Long-Distance Nuclear Matrix Elements for Neutrinoless Double-Beta Decay from Lattice QCD

Neutrinoless double-beta ($0\nu\beta\beta$) decay is a heretofore unobserved process which, if observed, would imply that neutrinos are Majorana particles. Interpretations of the stringent experimental constraints on $0\nu\beta\beta$-decay half-lives require calculations of nuclear matrix elements. This work presents the first lattice quantum-chromodynamics (LQCD) calculation of the matrix element for $0\nu\beta\beta$ decay in a multi-nucleon system, specifically the $nn \rightarrow pp ee$ transition, mediated by a light left-handed Majorana neutrino propagating over nuclear-scale distances. This calculation is performed with quark masses corresponding to a pion mass of $m_\pi = 806$ MeV at a single lattice spacing and volume. The statistically cleaner $\Sigma^- \rightarrow \Sigma^+ ee$ transition is also computed in order to investigate various systematic uncertainties. The prospects for matching the results of LQCD calculations onto a nuclear effective field theory to determine a leading-order low-energy constant relevant for $0\nu\beta\beta$ decay with a light Majorana neutrino are investigated. This work, therefore, sets the stage for future calculations at physical values of the quark masses that, combined with effective field theory and nuclear many-body studies, will provide controlled theoretical inputs to experimental searches of $0\nu\beta\beta$ decay.


I. INTRODUCTION
Neutrinos are the most poorly understood particles within the Standard Model.In the original conception of the Standard Model, they were presumed to be massless until the discovery of neutrino oscillations [1,2], which showed that the masses of at least two of the neutrino mass eigenstates are nonzero.The physical mechanism that generates neutrino masses, however, is still uncertain.If neutrinos are their own antiparticles, their masses could arise through a Majorana mass term Here, (ν iL ) C = C νT iL with C being the chargeconjugation matrix, and ν iL is a left-handed neutrino field for each of the mass eigenstates labelled by i ∈ {1, 2, 3}.These mass eigenstates are related to the flavor eigenstates ℓ ∈ {ν e , ν µ , ν τ } via ν ℓL = ℓ U ℓ,i ν iL , where U ℓ,i are the elements of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [3,4].Alternatively, if the Standard Model is extended to include yet-to-be-observed right-handed neutrinos ν iR , Dirac mass terms −m i νiL ν iR + h.c.arise naturally, for example, through a Yukawa coupling to the Higgs field analogously to that for the charged leptons.
Resolving whether neutrinos are their own antiparti-cles, that is, whether terms such as those in Eq. ( 1) are present, is one of the major open problems of modern particle physics.Since Eq. ( 1) permits lepton-number violation by two units, experimental probes of the Majorana nature of the neutrino search for processes that create and destroy neutrinos in pairs.Neutrinoful double-beta (2νββ) decay consists of two simultaneous electroweak nuclear transitions in the combined reaction nn → ppeeν eL νeL , where two neutrons (n) decay into two protons (p), two electrons (e), and two antineutrinos (ν eL ).This process is the rarest experimentally observed Standard Model process [5,6], and only occurs at measurable rates in nuclei that are stable against single-beta decay but favor a double-beta decay.If neutrinos are Majorana, then the two outgoing antineutrinos could mutually annihilate, resulting in a neutrinoless double-beta (0νββ)-decay nn → ppee , which could, in principle, occur in the same nuclei that can undergo 2νββ decay.Numerous experiments have searched for 0νββ decay [7][8][9][10] but, to date, none has conclusively shown that it occurs.At present, the most stringent bound on a 0νββ decay half-life is T 0νββ 2.3 × 10 26 yr at 90% C.L. for 136 Xe from the KamLAND-Zen experiment [9].
In any theory with a Majorana mass term as in Eq. ( 1), 0νββ decay can be induced via a light lefthanded neutrino propagating between two Standard-Model electroweak vertices, as depicted at the quark level in Fig. 1.Since the left-handed neutrino is nearly massless, the electroweak interactions can be widely separated (up to the diameter of the nucleus undergoing decay), so the resultant interactions are termed long distance.Beside this minimal extension of the Standard Model, many beyond-the-Standard-Model theories that allow for lepton-number violation, generate short-distance sixfermion (4-quark-2-electron) effective operators that can also induce 0νββ decay [11,12].The contributions of these operators to the π − → π + ee transition have been studied in Refs.[13,14] with the lattice-quantumchromodynamics (LQCD) framework.This work will not consider such short-distance scenarios and focuses on the long-distance mechanism in Fig. 1.
In the light Majorana-neutrino exchange mechanism, the necessity of a helicity flip between the electroweak-current insertions implies that the amplitude for 0νββ decay is proportional to the effective 0νββ neutrino mass defined as m ββ = i U 2 ei m i and to a hadronic or nuclear matrix element.A variety of nuclear models have been used to estimate the matrix elements in experimentally relevant nuclei, and significant differences exist between the values predicted by those models [10,15].The resultant model uncertainty can be roughly estimated (but not bounded) by the spread among model predictions and amounts to a factor of three or more.This results in large uncertainties when extracting a bound on m ββ from experimental constraints on half-lives.Reducing these uncertainties is crucial for interpreting experimental searches for 0νββ decay [16,17].
LQCD is a well-established non-perturbative technique for numerically evaluating hadronic and nuclear quantities rooted in quantum chromodynamics (QCD), the theory of the strong force [18,19].It, therefore, offers a firstprinciples method for determining hadronic and nuclear matrix elements relevant to 0νββ decay and has been previously used to study 2νββ decay [20,21].Nonetheless, the complexity of LQCD computations grows rapidly with baryon number, so initial calculations relevant to 0νββ decay have focused on the mesonic π − → π + ee transition [22,23] as a subprocess in a nuclear decay [11,[24][25][26].This work extends the approach developed in mesonic calculations to baryonic systems, the Σ − → Σ + ee and nn → ppee transitions, with the latter relevant to experimental studies in nuclei.The nn → pp transition cannot occur in free space due the the unbound initial state and the dominance of the single-beta decay mode.However, the transition amplitude is well defined and calculable with LQCD even in the absence of a nuclear medium, providing a promising avenue to isolate few-nucleon contributions to the full amplitude in large nuclei.The Σ − → Σ + transition also does not correspond to an experimentally observable decay mode (being much slower than the first-order weak decay of the Σ − to nπ − ); it is studied here to understand systematic uncertainties in the LQCD calculations more thoroughly than can be done with the nn → pp transition alone.
By themselves, the LQCD calculations presented here are not sufficient to determine nuclear matrix elements of phenomenological relevance but require connection to nuclear effective field theories (EFTs).EFTs provide a low-energy description of nuclear processes, including both neutrinoless and neutrinoful double-beta decay, in terms of a set of low-energy constants (LECs) that are a priori unknown parameters [27][28][29][30][31][32][33].Matching a (finitevolume) LQCD calculation of the nn → ppee transition amplitude to that expressed within a nuclear EFT allows the relevant LECs to be extracted [24,25,[34][35][36][37]. Once the systematic uncertainties associated with the present LQCD calculation are fully controlled in future studies, the constrained EFT can be used with many-body methods to calculate nuclear matrix elements in larger nuclei, hence reducing the model uncertainty that currently limits interpretation of experimental results.The present work explores prospects for the matching procedure to extract a leading-order LEC appearing in the pionless-EFT description of the nn → ppee process.

II. THEORETICAL AND COMPUTATIONAL APPROACH
This section presents the details of the theoretical and computational approach of this work.After introducing the physical 0νββ decay amplitude with a light Majorana neutrino in Sec.II A, Sec.II B demonstrates how such an amplitude can in principle be extracted from appropriate two-and four-point correlation functions in LQCD.A more thorough discussion of the exact mapping between the two quantities will be left to Sec.IV.
A. 0νββ decay amplitude in the long-range scenario At energies well below the electroweak scale, the Hamiltonian for single-β decay is given by where G F is the Fermi constant and V ud is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element encoding the down quark (d) to up quark (u) transition [38,39].At second order in perturbation theory, this interaction gives rise to a bi-local matrix element of the form [40] ⟨f |S (2) |i⟩ ≡ where S (2) is the second-order contribution to the weak interaction S-matrix, and ū1L and ū2L are the spinors of the outgoing left-handed electrons with momenta p 1 = (E 1 , p 1 ) and p 2 = (E 2 , p 2 ) and state normalization factors N e1 and N e2 , respectively.The quark-level lefthanded weak current is and T denotes the time-ordering operation.The neutrino propagator is given by where ) is the left-handed projector.In the last line, the neutrino propagator has factored into a product of Dirac matrices and a massless bosonic propagator, neglecting the neutrino mass compared with momenta characteristic of the hadronic scale.Finally, initial and final hadronic states are denoted by |N i ⟩ and |N f ⟩, and are assigned four-momenta p i,f = (E i,f , p i,f ), respectively.Importantly, the spatial momenta of the electrons are set to zero throughout, i.e., p 1 = p 2 = 0.The S-matrix element in Eq. ( 5) can be simply written as ⟨f |S (2) where the hadronic and leptonic tensors are defined via with Here, ū1,2 are the spinors corresponding to the outgoing electrons at rest.Equation ( 9) can be further processed by inserting a complete set of intermediate hadronic states |n⟩ (with energy Ẽn and momentum pn ) between the currents according to 1 = , then using x to transform the Heisenberg-picture currents back to the spacetime origin (with P 0 and P being the energy and momentum operators, respectively).One can then insert the form of the neutrino propagator in Eq. ( 8) and perform integrations over the spacetime coordinates and over the neutrino propagator to arrive at ⟨f |S (2) with Note that the sum over states |n⟩ involves an implicit integration over the total three-momentum of the intermediate state.
Considering that the expression in parentheses in the numerator of Eq. ( 14) is symmetric under the exchange of µ and ν indices, only the symmetric part of Γ µν contributes to the matrix element.Therefore, one may replace γ µ γ ν with γ {µ γ ν} = g µν in Eq. ( 14), giving Γ {µν} = g µν Γ with Γ = ū1 (1 + γ 5 )C ūT 2 .Taking u 1 and u 2 to have opposite spins (as is required by Pauli exclusion when outgoing momenta vanish), one can show that Γ = 1 up to normalization factors accounted for by N e1 and N e2 .
Finally, defining the amplitude one obtains This quantity encapsulates all of the strong-interaction dynamics of the 0νββ decay and is the target of the LQCD calculations discussed below.

B. 0νββ decay from LQCD correlation functions
LQCD calculations are performed in Euclidean spacetime to enable Monte Carlo methods.As a result, correlation functions and matrix elements extracted from them are defined in Euclidean spacetime.There are subtleties in the connection between Euclidean and Minkowski matrix elements of time-separated currents when on-shell intermediate states are produced [36,[41][42][43].Nonetheless, as will be discussed later, such states can be avoided in the present calculations; hence Euclidean and Minkowski matrix elements may be related simply by a phase from Wick rotation.As a result, this work will not distinguish Euclidean from Minkowski quantities hereafter but will state the relation between them when necessary.Furthermore, the LQCD study of this work is performed in the isospin limit, corresponding to degenerate up and down quark masses, and does not incorporate electromagnetic interactions.Additionally the electron mass is neglected, m e = 0, and consequently E i = E f ≡ E 0 in the 0νββ processes that are studied here.The formalism below is adapted to such a limit.Finally, all quantities are assumed to be defined in an infinite continuous spacetime volume throughout this section.The extension to a discretized finite volume is presented in Secs.III and IV.
To proceed, one can define two-point and four-point (Euclidean) correlation functions, which are calculable in LQCD (once spacetime is compactified and discretized).O i and O f are source and sink interpolating operators with the necessary quantum numbers to create the initial and final hadronic states for a given transition.A similar two-point function to Eq. ( 17) can be formed using the final-state interpolating operators but is equivalent to C 2 (t ′ ) in the isospin limit.Concrete choices for the interpolating operators will be discussed in Sec.III.
The integrals over the spatial coordinates project the final state and the two currents to zero momentum, so without loss of generality, the source interpolating operator O † i is placed at the spatial origin.After integrating over spatial coordinates as noted, the correlation functions only depend on the relative (Euclidean) time separations defined as The spectral decomposition of the bi-local matrix element in Eq. ( 18) is given by Here, ∆E n0 ≡ E n − E 0 denotes energy splitting between the ground state of the source interpolating operator and the nth excited state with the same quantum numbers, while ∆ E n0 = Ẽn − E 0 denotes energy splitting between the source ground state and the nth state with the quantum numbers of the intermediate hadronic system.Contributions from backwards-propagating states have been neglected (i.e., an infinite temporal extent is assumed).The factors A, B, and C are constants with respect to Euclidean time, expressible in terms of various excited-state matrix elements.The subleading terms represented by the ellipsis decay at least as quickly as e −∆E20tsrc or e −∆E20t snk .Similarly, the spectral decomposition of the two-point function takes the form where D is constant with respect to Euclidean time.
The connection to the amplitude in Eq. ( 16) is clearest for the ratio of four-point and two-point functions, which can be expressed as As indicated, this ratio depends on the three relative operator-time separations.It then follows that the (Euclidean) amplitude is given by

III. LQCD CALCULATION
The LQCD calculation in this work is performed on an ensemble of 12,136 QCD gauge-field configurations separated by 10 trajectories.The ensemble has a lattice spacing of a = 0.145 fm and a volume of (L/a) 3 × (T /a) = 32 3 ×48.Furthermore, sea quarks are implemented at the SU (3) flavor-symmetric point with degenerate up, down, and strange quark masses corresponding to a pion mass of m π = 806 MeV.The details of the gauge and fermion actions and the hybrid Monte Carlo scheme used to generate the ensemble are described in Ref. [44], with the same action used in other studies of few-baryon systems [20,[44][45][46][47][48][49][50][51][52][53][54][55].Of particular importance for this calculation, the proton, neutron, Σ 0 , Σ ± , and Λ are all degenerate, with a common mass of 1.64 GeV [44].

A. Interpolating operators
The single-baryon interpolating operators used in this work are where the superscript σ is a free spinor index, C = iγ 2 γ 4 is the Euclidean charge conjugation matrix, and P + = (1+γ 4 )/2 is the positive-parity projector. 1The color and spin contractions implicit in the preceding expressions are defined explicitly for an arbitrary set of three quarks (q i ∈ {u, d, s}) and products of Dirac matrices (Γ 1 , Γ 2 ) via where (α, β, σ, δ) and (a, b, c) are spin and color indices, respectively, and the square brackets visually isolate the diquark interpolating operator.The projection of all quarks to positive parity is appropriate for the large 1 The relation between the Minkowski and Euclidean γ matrices according to the convention of this work are 0 as given in Ref. [19].
FIG. 2. A schematic depiction of the four-point correlation function for the nn → pp transition used in this work.Quark propagators (solid lines) were constructed from a zeromomentum wall source and from point sinks.Extended propagators, defined in Eq. ( 31), are denoted by orange lines while the regular, spectator, propagators are shown in black.The neutrino propagator (dashed line) between ty and tx is given in Eq. (30).
quark masses used in the present calculation.The dinucleon interpolating operators are defined as where N ∈ {n, p} and the additional Cγ 5 couples the nucleon spins into the required spin-singlet combination.
For the nn → pp transition, the source and sink operators in the four-point function are

B. Propagator computation
The two-point correlation functions in Eq. ( 17) were computed with a wall source and a point sink.For the four-point correlation functions in Eq. ( 18), propagators were computed originating from both the source and the sink and contracted at the two operator positions x and y, as shown in Fig. 2. While Eq. ( 18) requires summing over all sink interpolating-operator positions, computing propagators from every point at the sink would be prohibitively expensive.Therefore, only a sparse grid of 4 3 sink points (corresponding to a sparsening factor of (L/a)/4 = 8 in each direction) was used.As studied in Ref. [56], this sparse grid corresponds to a partial three-momentum projection and does not modify the low energy spectrum. 2n each configuration, spatial grids of point sinks were constructed on every eighth timeslice.The computationally cheaper zero-momentum wall sources were computed on every timeslice in order to study the effects of varying source-sink separation.A total of 432 propagators were computed on each configuration.
Since all quarks in the interpolating operators in Eqs.(24) to (27) are projected to positive parity, only six (out of twelve) spin-color components of each propagator needed to be computed.The wall sources with zero three-momentum were constructed in Coulomb gauge with gauge fixing performed in GLU [57].Propagators were computed using the QPhiX inverters [58]. 3he bosonic propagator associated with the neutrino is defined in a finite periodic Euclidean spacetime in the LQCD calculation.Furthermore, the contribution from the spatial zero momentum is subtracted from the propagator: where the sum runs over non-zero finite-volume momenta and is truncated at |q| ≤ π/a to regulate the ultraviolet divergence at x = y.This form of the propagator is chosen to make matching to the nuclear EFT seamless [36]. 4he removal of the zero mode ensures that all intermediate states will be at a higher energy than the initial and final states for the volume used in this work, since the minimum neutrino energy is |q| = 2π/L.This approach avoids the difficulties of four-point correlation functions growing exponentially in operator separation times that affected π − → π + ee calculations with a massless intermediate state [22,23].

C. Contractions
The four-point correlation function is computationally expensive due to the number of Wick contractions involved and the sums over the sink and both current positions.First, extended propagators S αβ,µ ab (x) were built at the current insertion points x and y via where (J µ ) δζ is the Dirac structure of the weak current, the propagator S ζβ eb (x|x i ) originates at the source, S αδ ae (x f |x) is constructed from the propagator from the sink S δα ea (x|x f ) by γ 5 -hermiticity, and dependence on x i and x f is left implicit on the left-hand side.Then, at fixed operator times t x , t y , two extended propagators were combined with the bosonic propagator D(x − y) (without any spinor or color indices) to obtain a fourquark tensor with the discrete 3D Fourier transform F[f ](p; t) = x e ip•x f (x, t) computed efficiently using the fast Fourier transform implemented via the FFTW library [59] as in Ref. [23].
The tensor in Eq. ( 32) was then contracted with the spectator quark propagators connecting the source and sink interpolating operators as prescribed by Wick's theorem to form the four-point nn → pp and Σ − → Σ + correlation functions. 5The Σ − → Σ + correlation function is explicitly given as and the nn → pp correlation function includes N u !N d != (4!) 2 = 576 terms in the square brackets, each with three additional spectator quark propagators.Due to the link smearing and improvement in the gauge action and the clover term in the fermion action, time separations of at least three lattice units are required between the current-insertion points and either source or sink locations to avoid contamination from contact terms.Subject to this constraint, the four-point correlation function for the nn → pp transition was computed at all operator insertions for source-sink separations ranging from 6a to 16a, beyond which the statistical noise became prohibitively large.For the Σ − → Σ + transition, where the statistical noise was milder, contractions were computed for all separations less than T /2 = 24a.

D. Renormalization
The left-handed electroweak-current insertion J µ = 1 2 ūγ µ (1 − γ 5 )d is the difference of vector and axial-vector insertions.The local lattice currents for these two contributions renormalize separately, so the renormalized current insertion has the form Due to the interference between the two insertions of these terms in the four-point function, the renormalization factors (or at least the relative renormalization Z V /Z A ) are included at the time the correlation functions are computed.The renormalization factors for the action parameters used in this work have been computed in Ref. [60]: E. Extraction of matrix elements

Analysis of two-point functions
The ground-state energies m n = m Σ and E nn = E pp are extracted from the respective two-point functions given in Eq. (17). Figure 3 shows the effective-mass functions for the Σ and nn correlation functions, where aE eff (t ′ ) = ln (C 2 (t ′ )/C 2 (t ′ + a)).Results for fitting the effective mass to a constant using correlated χ 2 minimization are given on the right of Fig. 3 as a function of the minimum time used in the fit.For Σ, t ′ min ∈ {10, . . ., 19}.For nn, t ′ min ∈ {9, . . ., 13} and a cut of t ′ max = 16 is imposed to restrict to points where the statistical noise for the two-point function remains below 30% of the central value.Fits with smaller values of t ′ min were conducted but resulted in poor fit quality (χ 2 /dof > 2, where dof denotes the number of degrees of freedom) and are therefore not shown.Stability at the level of one standard deviation is observed for the masses extracted from different fits.The horizontal bands show the result of combined averages and uncertainties using weights based on the Akaike Information Criterion (AIC) [61].The final results for the masses in lattice units are Note that the interpolating operators used in this work are different from those used in previous studies but yield masses consistent with these earlier studies [44,54,55,62].At the level of precision achieved in this study, the dineutron is consistent with either a bound state or a scattering state.t /a

Analysis of four-point functions
The extraction of nuclear matrix elements from a LQCD calculation of the ratio R(t snk , t, t src ) defined in Eq. ( 21) requires controlling excited-state contributions from the source and sink in Eq. ( 22), followed by extrapolation and integration over the current separation as in Eq. (23).A two-step analysis procedure is used.First, for fixed current separations, the Euclidean time dependence is modeled with respect to the source and the sink locations to remove excited-state contributions.The output of the first step is therefore   21) as a function of temporal separation t between the two currents.The upper data are for nn → pp, while the lower data are for Σ − → Σ + .In both cases, the colors distinguish different values of current-sink separations t snk as indicated in the legends.For clarity, the points at fixed (t, t snk ) have been slightly offset in the horizontal direction as tsrc varies.
Second, the integral in Eq. ( 23) must be evaluated to determine the amplitude A i→f .Equation (22) shows that R i→f (t) decays as q,n e −(|q|+∆En0)|t| .As shown concretely below, at the present statistical precision and at finite lattice spacing, the sum can be well approximated by a single exponential where E (R) and A (R) are an effective energy gap and amplitude associated with the asymptotic ratio R i→f (t).Departures from this behavior, arising from the full spectrum of states in the sum q,n e −(|q|+∆En0)|t| are expected at short times.However as discussed above, the short-time data (t/a ≤ 2) are sensitive to details of the lattice discretization and are excluded from this analysis; subsequent calculations at finer lattice spacings will likely reveal additional contributions to the amplitude from these higher-energy states.Since these cannot be resolved in the current study, however, the required integral in Eq. ( 23) can be approximated as LQCD results for the ratios R Σ − →Σ + (t snk , t, t src ) and R nn→pp (t snk , t, t src ) are shown in Fig. 4, displayed as a function of the temporal separation between the currents.An alternative view of the data, focusing on the source and sink separations, is given in Fig. 5 for both Σ − → Σ + and nn → pp.As expected from the spectral decomposition, excited-state contamination is generically present from both the source and the sink.The one exception is for the source-time dependence of R nn→pp , which at the present level of precision is statistically consistent with a constant. 6 First-stage fits: R i→f (t snk , t, t src ) → R i→f (t).For fixed current separation t, the data are fit to Eq. ( 22).For Σ − → Σ + , only the leading contributions proportional to e −∆E10t snk and e −∆E10tsrc are retained (with unknowns A, B, and ∆E 10 ).For nn → pp, only the contribution proportional to e −∆E10t snk is included (with unknowns B and ∆E 10 ), as no dependence on t src is observed within uncertainties.Examples of the resulting fits are shown by the solid black curves in Fig. 5.The limiting value of R i→f (t snk , t, t src ) → R i→f (t) emerging from the fit is shown by the common horizontal line.The fit displayed in the upper row of Fig. 5 for Σ − → Σ + has χ 2 /dof of 1.04 for dof = 240; the fit in the lower row for nn → pp has χ 2 /dof of 1.15 for dof = 72.Fits of similar quality are obtained for each fixed temporal separation 6 While the fits appear to control excited state contamination well, there is always the possibility of low-lying excited states distorting the results of LQCD calculations, and this concern is of particular importance in the nn → pp transition due to the dense low-lying spectrum in nuclear systems [54].Further study with a variety of interpolating operators would be beneficial to confirm the plateau values observed in this work.t snk /a  The ratio of four-point and two-point correlation functions, R i→f (t snk , t, tsrc), defined in Eq. ( 21) for a fixed current separation t/a = 3.The left (right) column shows the dependence on the sink-current (source-current) separation.To show the simultaneous dependence on both t snk and tsrc, the same data appear in both columns, and matching points appear in the same color on the left and right.The solid black curves show the result of a correlated fit to the all the data displayed for a given process (Σ − → Σ + or nn → pp).In each row, the limiting value of R i→f (t) determined from the fit is shown by the common horizontal line.
of the currents, yielding R i→f (t) as a function of t.The results of this process are shown in Fig. 6.
To verify stability of the fitting procedure, the values of (t min snk , t min src ) included in the fit are varied for each fixed t with t min snk and t min src varied independently in {3, 4, 5, 6}, which modulates the size of excited-state effects.To account for any variation in the output values for R i→f (t), the results at fixed t are combined using model averaging with AIC weights [61] to yield the black points in Fig. 6.
Second-stage fits: R i→f (t) → A i→f .As shown in Fig. 6, R i→f (t) is saturated by a single decaying exponential for t/a ≥ 3 in both panels.This statement is illustrated in Fig. 7, which shows the effective energy and effective amplitude For Σ − → Σ + , both quantities exhibit clear plateaus before statistical noise begins to dominate at large times.
For nn → pp, the data are noisier but consistent with a constant.The data for R i→f (t) are fit to Eq. ( 39), varying t min ∈ {3, 4, . . ., 7} and t max ∈ {t min + 3, t min + 4, . . ., t max max } to check for stability, where the variations in t max extend to t max max = 10 for the nn → pp transition.For the Σ − → Σ + transition, the data were clean enough to allow t max max to be extended to 17, and a single Pull FIG. 6.The asymptotic ratio R i→f (t) shown on a logarithmic scale for Σ − → Σ + (upper panel) and nn → pp (lower panel).Each cluster of colored points represents fit results at fixed t (with varying tsrc and t snk ) such as those shown in Fig. 5.The results at each fixed t are combined using model averaging with weights based on the AIC to yield the black points.The gray line and error band show the result of second-stage fits to model the dependence on the currentcurrent separation for t/a ≥ 3. The bottom of each panel displays the pull, i.e., the difference between the fit and data in units of the uncertainty.Points excluded from the fit appear in light gray.
exponential still sufficed for the second-stage fit. 7Results are combined using weights based on the AIC, with the final posterior values for E (R) and A (R) indicated by the horizontal bands in Fig. 7. Due to correlations, the 7 The statistically cleaner Σ − → Σ + channel provides a useful check on the systematic uncertainties of the second stage of the nn → pp analysis; fits to R Σ − →Σ + (t) for t max max = 10 are consistent within uncertainties with those with t max max = 17, and consequently, fits with t max max = 10 were also deemed sufficient for the nn → pp case.
eff and amplitude eff from fits to the ratio R i→f (t) for Σ − → Σ + (upper panel) and nn → pp (lower panel).The horizontal lines and error bands show the final posterior results from fits to the exponential decay in Eq. (39).The amplitudes have been re-scaled by arbitrary factors for ease of visualization.
uncertainty in E (R) is somewhat smaller than suggested visually by E (R) eff in Fig. 7.The gray bands in Fig. 6 show the fit results against the data for R i→f (t).The posterior values for E (R) and A (R) can then be used to evaluate the integral in Eq. (40).
The final values for the renormalized amplitudes are where the final uncertainties include both statistical uncertainties and systematic uncertainties from the model averaging as well as the uncertainty arising from Z V in Eq. (35).The renormalized amplitude for Σ − → Σ + is determined with a fractional uncertainty of roughly 10%, of which the dominant uncertainties are the ratio ) and Z V (≈ 5%).The relative breakdown in similar for nn → pp.The small (few-percent) uncertainty in the ratio Z A /Z V is neglected in this work, since it would require recomputing all of the contractions while propagating this uncertainty.

IV. PROSPECTS FOR NUCLEAR EFT MATCHING
Direct LQCD calculations of 0νββ amplitudes in experimentally relevant nuclear isotopes are beyond the reach of the current computational paradigm.The reasons include a substantial increase in complexity of quarklevel nuclear correlation functions with increasing atomic number, a severe signal-to-noise degradation of correlation functions as a function of Euclidean time and atomic number, and nuclear excitation gaps that are small compared to the QCD scale which thus demand unrealistically precise spectral resolution.As a result, nuclearstructure calculations based on nucleonic degrees of freedom, and nuclear-level Hamiltonians and currents, will be the primary method to access phenomenologically relevant nuclear matrix elements for the forseeable future.These Hamiltonians and currents can be systematically constructed from few-nucleon EFTs, assuming the existence of reliable power-counting schemes.Nonetheless, such a program is limited by the lack of knowledge of input interactions at the few-nucleon level, particularly for the 0νββ process, which has not yet been observed, and importantly, does not occur naturally in few-nucleon systems.As a result, fully controlled LQCD input at or near the physical values of the quark masses will be crucial in order to constrain unknown low-energy constants (LECs) of the EFTs.Pionless EFT is a commonly used theoretical framework for studying few-nucleon processes at low energies [27][28][29][30].Pionless EFT was applied to the 0νββ decay in Refs.[24][25][26] to determine the amplitude for nn → ppee process at the lowest EFT orders.Nonetheless, it was found that the EFT amplitude is undetermined for the long-range scenario even at leading order due to the presence of an unknown short-distance LEC, called g N N ν (µ), which characterizes the strength of the four-nucleon-two-electron contact interaction at a given renormalization scale, µ.Later studies provided various estimates of this coupling based on a dispersive analysis [63,64] and large-N c considerations [65].However, there remain significant model dependence and uncertainty in these determinations, which have been shown to lead to an amplified uncertainty in the nuclear matrix elements in larger nuclear isotopes [66].Ultimately, LQCD will be able to provide a first-principles determination of this LEC.Such calculations, nonetheless, provide the values of matrix elements in a Euclidean finite spacetime volume, which need to be connected to the physical amplitudes in the corresponding EFT.
Such a formalism for the case of leading-order pionless EFT was developed in Ref. [36].Explicitly, the amplitude, defined in Eq. ( 15), can be related to the leadingorder LEC of the EFT by the following matching relation: Here, pi and pf are the nonrelativistic binding momenta defined as pi,f = m n E i,f for energy shifts E i,f and the dependence of A nn→pp on these momenta has been made explicit.M(p) denotes the elastic two-nucleon scattering amplitude in the spin-singlet channel, which can be approximated by an effective-range expansion: with scattering length a and effective range r. g N N ν (µ) in Eq. ( 45) is a dimensionless constant related to the LEC g and J ∞ (µ) is a known function given by with γ E being Euler's constant [24][25][26].Furthermore, R and δJ V are two finite-volume functions, whose forms are given in Refs.[36,37].Compared with the matching relation in Eq. ( 28) of Ref. [36] which connects the absolute values of the left and right-hand sides of Eq. ( 45), this work resolves the sign ambiguity in this equation so as to allow for a unique constraint to be placed on the LEC g N N ν (µ).In the isospin limit where pi = pf ≡ p, the relation can be simplified as Despite the relation described in this section, and the LQCD results obtained for A nn→pp (p) in this work, several caveats preclude a rigorous determination of g N N ν (µ) via Eq.( 49) at the present time.First and foremost, the LQCD matrix element here is obtained at unphysically large quark masses.Clearly, it is the value of g N N ν (µ) with the physical quark masses that is of phenomenological interest and, a priori, the quark-mass dependence of such an LEC is unknown.Therefore, an attempt to constrain g N N ν (µ) or the renormalization-scale-independent quantity at the quark masses of this work will likely have little bearing on the physical value of the coupling.Nonetheless, one may still obtain an estimate of the value of this LEC at the quark-mass value of this work, in which case the corresponding values of two-nucleon scattering parameters need to be used in the matching relation.To date, there are two classes of LQCD computations of low-energy two-nucleon spectra and scattering parameters at m π ≈ 800 MeV via the use of Lüscher's finite-volume formalism.The earlier computations involve asymmetric two-nucleon correlation functions, and point to the existence of rather deep bound states in both the spin-singlet and spin-triplet two-nucleon channels [44,50,52,55,67,68].These were subsequently used to constrain the relevant LECs in electromagnetic and weak reactions of two-nucleon systems at various pion masses and allowed preliminary extrapolations to the physical point [20,21,45,69,70].However, at the finite-volume ground-state two-nucleon energy, which sets the kinematics of the amplitude in this work, the pionless EFT converges poorly when using the values for the effective range and scattering length in those studies.Therefore, obtaining the desired 0νββ-decay amplitude using those results requires extensions of the current leading-order matching formalism, or the use of alternate power-counting schemes.The other set of calculations at m π ≈ 800 MeV build symmetric correlation functions to enable accessing the low-lying spectra via a variational method.These lead to upper bounds on ground-state energies that are also consistent with less bound or unbound two-nucleon systems within uncertainties [54,62,71].No bound states are seen in complementary studies using the Bethe-Salpeter potential method [72,73].While the associated scattering length and effective range for these bounds allow the use of the leading-order matching formalism here, it is non-trivial to turn variational bounds on the energies to bounds on the desired LEC of the EFT, given the nonlinearity of the matching relation.
Despite these caveats, the matching to the EFT amplitude using the above calculation of A nn→pp , leads to gNN ν (µ = m π = 806 MeV) values that differ by a factor of four depending on whether the non-variational determinations of two-nucleon energy and scattering parameters or those from the variational studies are used (assuming the variational bounds are saturated).In both cases, the extracted values are within an order of magnitude of the phenomenological estimate of Ref. [64].Consequently, increasingly controlled determinations of the two-nucleon quantities that are input to the matching relation are needed for a robust determination of this LEC.For calculations with physical quark masses, such two-nucleon quantities are well determined phenomenologically, which would ease the matching procedure.
Improving on this situation thus requires calculations of A nn→pp and the finite-volume two-nucleon spectrum at or near the physical quark masses.A point worth emphasizing is that the pionless EFT converges at the finite-volume ground-state energy of the spin-singlet twonucleon system, provided that the lattice volume is sufficiently large, hence putting another requirement on future calculations.For an exploration of the impact of volume on the determination of g N N ν (µ) at the physical values of quark masses, see Ref. [37].

V. SUMMARY AND CONCLUSION
Within the coming few decades, the sensitivity of experimental neutrinoless double-beta decay searches is projected to increase by several orders of magnitude, corresponding to an order of magnitude decrease in the effective 0νββ masses that can be probed [16].Given current best estimates of nuclear matrix elements, these experiments will likely-but not definitively-be sensitive to the entirety of the parameter space for the inverted hierarchy of neutrino masses.These searches thus have a large discovery potential but also present the possibility of definitively ruling out the Majorana nature of the neutrino if they find no such decays and if neutrino oscillation experiments confirm the inverted mass hierarchy.Thus, either positive or negative results in next-generation experiments will shed crucial light on this problem provided that the dominant mode of decay is via the exchange of a light Majorana neutrino and that the corresponding nuclear matrix elements can be computed accurately to extract m ββ from measured (bounds on) half-lives.
Starting with the low-energy constants from nuclear effective field theories, nuclear many-body theories can provide ab initio calculations of binding energies and 0νββ matrix elements in light to moderate (A ≲ 48) nuclei [74,75].For heavier nuclei (16 ≲ A ≲ 132), EFT-based approximations to nuclear physics can predict 0νββ half-lives with more control than the nuclear models currently used [76][77][78].As such, determining these low-energy constants in the timescales relevant for these next-generation experiments is of substantial importance to the nuclear-and particle-physics communities [16,17].
This work presents the first LQCD calculation of the long-distance 0νββ-decay amplitude of a nuclear system, yielding the result a 2 A nn→pp = 0.078 (16) (51) on a single LQCD ensemble with a lattice spacing of a = 0.145 fm, a lattice volume of (L/a) 3 ×T /a = 32 3 ×48, and quark masses corresponding to a pion mass of m π = 806 MeV.The baryonic transition A Σ − →Σ + was also determined for the first time.While this calculation was performed at quark masses that are too large to match to experiment directly, it shows that the relevant matrix elements are calculable in LQCD in multi-baryon systems.This work further discusses prospects for the determination of the leading-order pionless-EFT LEC g N N ν from the LQCD matrix element.Repeating this calculation at lighter quark masses will be non-trivial due to the exponentially worsening signal-to-noise problem as the lightquark masses decrease, a problem especially challenging in multi-baryon systems.However, such calculations are important, as they are the only way to determine experimentally relevant values for the LECs of the nuclear EFTs in a model independent way.

1 / 2 >FIG. 1 .
FIG.1.The quark-level diagram responsible for the longdistance contribution to neutrinoless double-beta decay, corresponding to light left-handed Majorana-neutrino exchange between two W bosons.

FIG. 4 .
FIG. 4.The ratio of four-point and two-point correlation functions defined in Eq. (21) as a function of temporal separation t between the two currents.The upper data are for nn → pp, while the lower data are for Σ − → Σ + .In both cases, the colors distinguish different values of current-sink separations t snk as indicated in the legends.For clarity, the points at fixed (t, t snk ) have been slightly offset in the horizontal direction as tsrc varies.
FIG. 5.The ratio of four-point and two-point correlation functions, R i→f (t snk , t, tsrc), defined in Eq. (21) for a fixed current separation t/a = 3.The left (right) column shows the dependence on the sink-current (source-current) separation.To show the simultaneous dependence on both t snk and tsrc, the same data appear in both columns, and matching points appear in the same color on the left and right.The solid black curves show the result of a correlated fit to the all the data displayed for a given process (Σ − → Σ + or nn → pp).In each row, the limiting value of R i→f (t) determined from the fit is shown by the common horizontal line.