The path from lattice QCD to the short-distance contribution to $0\nu\beta\beta$ decay with a light Majorana neutrino

Neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay of certain atomic isotopes, if observed, will have significant implications for physics of neutrinos and models of physics beyond the Standard Model. In the simplest scenario, if the mass of the light neutrino of the Standard Model has a Majorana component, it can mediate the decay. Systematic theoretical studies of the decay rate in this scenario, through effective field theories matched to \emph{ab initio} nuclear many-body calculations, are needed to draw conclusions about the hierarchy of neutrino masses, and to plan the design of future experiments. However, a recently identified short-distance contribution at leading order in the effective field theory amplitude of the subprocess $nn \to pp\,(ee)$ remains unknown, and only lattice quantum chromodynamics (QCD) can directly and reliably determine the associated low-energy constant. While the numerical computations of the correlation function for this process are underway with lattice QCD, the connection to the physical amplitude, and hence this short-distance contribution, is missing. A complete framework that enables this complex matching is developed in this paper. The complications arising from Euclidean and finite-volume nature of the corresponding correlation function are fully resolved, and the value of the formalism is demonstrated through a simple example. The result of this work, therefore, fills the gap between first-principle studies of the $nn \to pp\,(ee)$ amplitude from lattice QCD and those from effective field theory, and can be readily employed in the ongoing lattice-QCD studies of this process.

Introduction.-The lepton-number violating process (A, Z) → (A, Z + 2) + ee, with A and Z being, respectively, the atomic and proton numbers of a parent nucleus, if observed, will mark a major discovery. Beyond its confirmation of the presence of a Majorana component to the neutrino mass [1], which would be unlike any other fermions in the Standard Model (SM), our knowledge of beyond-SM (BSM) mechanisms that may be responsible for this decay can only be enhanced by combining theoretical calculations of the rate, and other decay observables, with experimental findings [2][3][4]. Furthermore, planned experimental endeavors will crucially benefit from theoretical predictions of the expected rates in various isotopes given the BSM scenarios considered [2,3,[5][6][7]. A widely considered scenario is a minimal extension of the SM in which the light neutrinos of the SM are promoted to Majorana neutrinos, which by virtue of being their own anti-particles, can be emitted and reabsorbed by the nucleus undergoing the decay. Given the estimated masses of the neutrinos in the sub-eV to a few-eV range [8], the corresponding nuclear matrix element is long range in nature and receives contributions from intermediate nuclear states. Nonetheless, the currents involved in the transition are the weak current of the SM.
Despite the long-range nature of the process, recent nuclear effective field theory (EFT) analyses of the elementary subprocess nn → pp (ee) have revealed a shortdistance contribution to the amplitude at leading order (LO), with a low-energy constant (LEC) of the corresponding isotensor contact operator that absorbs the ultraviolet (UV) scale dependence of the amplitude through Renormalization Group (RG) [9,10]. As such a subprocess cannot be observed in free space, and given the program that has been formed around the use of nuclear EFTs to systematically improve the ab initio nuclear structure calculations of the nuclear matrix elements [11][12][13][14][15] toward experimentally-relevant isotopes, the missing knowledge of the value of such a shortdistance contribution appears to impede progress, and has promoted indirect but less certain estimations [10].
Lattice QCD, which numerically solves QCD on a finite grid in an Euclidean spacetime, has the promise of reliably constraining the EFTs of 0νββ in the few-nucleon sector [16][17][18], and has already demonstrated its reach and capability in constraining pionic matrix elements for lepton-number violating processes π − → π + (ee) and π − π − → ee within the light-neutrino scenario [19][20][21], the π − → π + (ee) process within a heavy-scale scenario [22], as well as the two-neutrino double-β decay (2νββ) of a two-nucleon state [23,24] (the latter yet at unphysically large quark masses due to the computational cost). Lattice-QCD matrix elements for these processes, however, lack certain complexities compared with the desired nn → pp (ee) process with a light Majorana neutrino, whose determination is the key to matching to the EFT descriptions developed in recent years [9,10,25]. While the numerical evaluations of the matrix elements are underway, the interpretation of these matrix elements in terms of the physical amplitude, and their matching to EFTs have so far been missing from the course of developments. In this paper, such a framework will be developed and presented for the first time, and the steps involved in the matching will be described. This frame-work, along with a realistic example that will be outlined, therefore, demonstrate how the results of this work can be used in the ongoing and upcoming studies to obtain the short-distance LEC of the EFT from the lattice QCD output.
The matching framework of this work builds upon major developments in recent years in accessing local and non-local transition amplitudes in hadronic physics from the corresponding finite-volume matrix elements in Euclidean spacetime obtained with lattice QCD [26][27][28][29][30][31][32][33][34][35][36], and in particular, a recent work on developing a similar formalism for the two-neutrino process nn → pp (eeν eνe ) [37]. Nonetheless, the neutrinoless process involves additional complexities due to a propagating neutrino in the intermediate state, requiring new components to be included in the matching condition. The matching involves constructing the relation between (Minkowski and infinite-volume) physical amplitudes in the EFT and the Minkowski finite-volume matrix element in lattice QCD as the first ingredient, as well as the relation between the latter and its Euclidean counterpart, which is complicated by the presence of on-shell multi-particle intermediate states, as the second ingredient. The Euclidean finite-volume matrix element is what is provided by lattice QCD, and hence the connection to the physical amplitude is built with these two ingredients.
EFT amplitude at leading order.-In a SM EFT of 0νββ decay [25,[38][39][40][41][42][43][44], the lepton-number (L) violating operator with the lowest mass dimension is a Majorana mass term, Here, C = iγ 2 γ 0 denotes the charge conjugation matrix, ν L is the left-handed (electron) neutrino field, m ββ = i U 2 ei m i is the effective neutrino mass, with U ei being the elements of the Pontecorvo-Maki-Nakagawa-Sato (PMNS) matrix [45,46]. m i is the mass of the neutrino flavor i. The contributions from higher-dimensional operators to the 0νββ decay can be also significant (considering the suppression from the small m ββ factor form the lowest-dimensional Majorana mass operator). Nonetheless, the scenario in which the the decay is induced by the operator in Eq. (1) is the most natural possibility and will be the focus of this work. While the decay can only proceed in certain nuclear media, the subprocess to be studied is nn → pp (ee), in which the two d quarks in the initial two-neutron state interact with the left-handed W boson of the SM, the emitted antineutrino would turn into a neutrino through the operator in Eq. (1) as it propagates, and gets absorbed by the two u quarks mediated by another W boson (and two electrons will be emitted in the final states). Since quarks are bound to nucleons and nucleons interact via the non-perturbative strong force, to relate the rate of the decay in nn → pp (ee) to the underlying SM EFT, one needs to map this problem to a nuclear EFT, and constrain the EFT using a direct calculation of the matrix element with lattice QCD.
The nuclear EFT considered here is the pionless EFT [47][48][49][50][51], where the Lagrangian of free and strongly interacting nucleons can be organized as (2) Here, ∂ t is the time derivative and ∇ is the spatial gradient operator. N = ( p n ) is an isospin doublet comprised of the proton, p, and the neutron, n, fields, each with mass M . Isospin symmetry will be assumed throughout as the isospin-breaking effects are small and contribute at higher orders than considered. P i ≡ 1 √ 8 σ 2 τ 2 τ i is a projector for the isotriplet channel, and the ellipsis denotes higher-order terms in a momentum expansion. A similar interacting term can be written for the isosinglet channel, however, it will not be of use in the following LO study of the decay amplitude. The effective Lagrangian for the charged-current (CC) weak interaction is given by [52,53], v and S are the nucleon velocity and spin, respectively (v = (1, 0) and S = (0, σ 2 ) in the nucleon's rest frame), τ + = (τ 1 + iτ 2 )/2 where τ i are isospin Pauli metrices, and g A is the nucleon's axial charge. The leptonic current contains the left-handed electron, e L , and neutrino, ν L , fields. Last but not least, one can construct a contact ∆L = 2 four-nucleon-two-electron operator in the EFT: where P + = (P 1 + iP 2 )/2 [9,10,54]. 1 While naive dimensional analysis suggests that this operator must contribute at a high order, RG considerations require promoting this operator to leading order, as derived in Refs. [9,10], and will be discussed shortly. The full transition amplitude for the nn → pp (ee) process is not separable to the hadronic and leptonic amplitudes as in the two-neutrino process given the presence of a neutrino that propagates between the two weak currents. Nonetheless, the contribution from final-state electrons (as well as constants proportional to G F and V ud ) 1 Our sign convention for L (∆L=2) N is opposite to that in Ref. [10], but it reproduces the amplitude in the same reference. Our convention is in agreement with Refs. [9,54].
can still be separated from a hadronic amplitude that includes the hadronic matrix element that is convoluted by the neutrino propagator. This latter contribution is what one would evaluate in lattice QCD and match to nuclear EFTs. We assume a simple kinematic in which the total three-momenta of the system is zero, and the electrons are at rest, each having energy E 1 = E 2 = m e , where m e is the electron's mass. Furthermore, at the LO in the EFT, two further simplifications arise: i) only s-wave interactions of the nucleons contrinue, ii) the amplitude only receives contributions from a static neutrino potential, and contributions from the small non-zero neutrino mass in the denominator of the neutrino propagator, as well as radiative neutrinos, can be ignored. The mixed hadronic-leptonic amplitude can then be written as nn→pp denotes contributions in which the neutrino propagates between two external nucleons. These will not matter for matching to the lattice-QCD matrix element. On the other hand, M (Int.) nn→pp denotes contributions in which the neutrino propagates between two nucleons dressed by strong interactions on both sides. This amplitude depends upon the short-distance LEC g N N ν through: Here, E i and E f denote the energy of the incoming twoneutron state and the outgoing two-proton state, respectively. M is the LO strong-interaction scattering amplitude of the isotriplet channel, and J ∞ is a function representing the s-channel two-loop diagram with an exchanged Majorana neutrino: The integral is divergent in the UV and in the dimensional regularization scheme is regularized to where L(E i , E f ; µ) ≡ ln The UV divergence of the loop function necessities the introduction of a counterterm at the same order, i.e., g N N ν .
that can be constrained with lattice QCD, as will be shown shortly. g N N ν (µ) can then be determined using the value of C 0 and J ∞ at a given renormalization scale µ.
Matching between finite and infinite volume.-Keeping the Minkowski signature of spacetime intact, we now consider a finite spatial volume with cubic geometry and with extent L along each Cartesian coordinates with periodic boundary conditions. The time direction is assumed to be infinite. The correlation function at LO in the EFT is depicted in Fig. 1, and can be written as Here, B † nn and B pp are the matrix elements of initial-and final-state interpolating operators between vacuum and on-shell "in" and "out" two-nucleon states, respectively. The ellipsis denotes terms that will not matter for the matching relation below, see Ref. [37] for further detail in a similar process. F is a finite-volume function defined as with The discretized energy eigenvalues of the two-nucleon system in a finite volume, E m , are obtained from the "quantization condition" F −1 (E m ) = 0 [55,56]. Finally, a new finite-volume function δJ V , corresponding to the two-loop diagram with the exchanged neutrino propagator needs to be evaluated: where in the summations, k 1 , k 2 ∈ 2π L Z 3 . This sumintegral difference can be evaluated numerically for given values of E i and E f , the detail of which is presented in the Appendix. The requirement k 1 = k 2 removes the zero spatial-momentum mode of the neutrino in the loop to render the finite-volume sum finite. Correspondingly, the finite-volume correlation function in lattice QCD will need to implement a zero-mode regulated neutrino propagator to match to this expression. Such a treatment of the infrared singularities in a finite volume is customary in the lattice QCD+QED studies of hadronic masses [57][58][59], decay amplitudes [60][61][62], and two-hadron scattering [63,64]. Other prescriptions have been proposed and implemented for lattice-QCD calculation of the π − → π + (ee) amplitude [20,21]. These can be generalized to be applicable to the nn → pp (ee) process following similar steps as those in this work.
To proceed with finding the matching relation, one notes that the finite-volume correlation function in Eq. (10) has the same general structure as that for the two-neutrino process obtained in Ref. [37]. As a result, all steps introduced in Ref. [37] can be closely followed to obtain the matching relation between finite and infinitevolume matrix elements. In particular, upon Fourier transforming Eq. (10) with E i and E f to form the correlation function in the mixed time-momentum representation, and comparing it against the same correlation function that is obtained from a direct four-point function upon inserting complete sets of intermediate finitevolume states between the currents, one arrives at where R(E n ) = lim E→En (E − E n ) F(E), and E n i(f ) denotes a finite-volume energy of the initial (final) twonucleon state. T Here, J =qτ + γ µ (1 − γ 5 )q with q = ( u d ), which can be implemented in lattice-QCD calculations. At the hadronic level, it matches to N † τ + (v µ − 2g A S µ )N in Eq. (3). Nonetheless, being a quark-level current means that T (M) L also incorporates the contact ∆L = 2 interaction in Eq. (4). S ν (z 0 , z) denotes the Minkowski finitevolume propagator of a Majorana neutrino, with its zero spatial-momentum mode removed: Minkowski to Euclidean matching.-The quantity T (M) L in Eq. (15), whose connection to the physical amplitude was established in Eq. (14), is defined with a Minkowski signature. On the other hand, with lattice QCD only Euclidean correlation functions can be evaluated. Unfortunately in the case of non-local matrix elements, generally one cannot obtain the former from the latter upon an analytical continuation [35]. To appreciate the subtlety involved, and to introduce a procedure that, nonetheless, allows constructing the Minkowski matrix element from its counterpart in Euclidean spacetime, one may consider the following correlation function: which can be computed directly with lattice QCD. 2 τ ≡ iz 0 is the Euclidean time, and the subscript (E) on quantities is introduced to distinguish Euclidean forms from Minkowski counterparts. In particular, S (E) ν is the Euclidean neutrino propagator in a finite volume with its zero spatial-momentum mode removed, It is now clear that simply integrating over the Euclidean time with weight e E1τ can be problematic if on-shell intermediate states are allowed. Here, E 1 is the energy of the first or the second electron depending on the time ordering. This can be seen by expressing the Heisenberg-picture operator in Euclidean spacetime as J (E) (τ, z) = eP 0τ −iP ·z J (E) (0) e −P0τ +iP ·z , whereP 0 andP are energy (Hamiltonian) and momentum operators, respectively, and upon inserting a complete set of single-and multi-particle states in the volume between the two currents. Without loss of generality, we assume that J (E) (0) is the same as its Minkowski counterpart, i.e., it has no phase relative to the Minkowski current at the origin. The Euclidean superscript of the Schrödingerpicture currents will therefore be dropped. It then becomes clear that for those values of intermediate-states energies and momenta such that the integration over Euclidean time with e E1τ will be divergent.
Here, E * m are the finite-volume energy eigenvalues of the intermediate spin-triplet two-nucleon state with total momentum P * m , and we assume that threeparticle intermediate states with on-shell kinematics are not possible given the initial-state energy. The problematic contributions satisfying conditions i and ii can be subtracted from Eq. (17), leaving the rest to read The spectral decomposition of T (E) ≥ L has, therefore, exactly the same form as the Minkowski counterpart upon an overall i factor. Here, (21) where it is assumed that there are N (N ) states satisfying condition i (ii) above, and The remaining contributions arising from on-shell intermediate states, called T (E) < L , can be formed separately with the knowledge of the single-current matrix elements in a finite volume between the initial (final) and intermediate states: Eqs. (23) and (20) can now be combined to construct the desired Minkowski quantity T whose relation to the physical nn → pp (ee) amplitude at LO in the EFT was already established in Eq. (14). This completes the matching framework that relates M (Int.) nn→pp in Eq. (6), and hence the new shortdistance LEC g N N ν , to the lattice-QCD correlation function G (E) L (τ ) in Eq. (17). It should be noted that the single-current matrix elements required for this matching relation, i.e., those appearing in Eq. (22), can themselves be evaluated with lattice QCD, and can be matched to the physical amplitude for the single-β transition amplitude [29,37].
Discussion and outlook.-Given significant progress in lattice QCD studies of nuclear matrix elements in recent years [17,23,[65][66][67][68][69], albeit yet with unphysical quark masses, it is expected that lattice QCD will be able to evaluate the four-point correlation function in Eq. (17) at the physical values of the quark masses, along with the required two-and three-point functions that allow the construction of the finite-volume Minkowski amplitude in Eq. (24). This can then be used in Eq. (14) to constrain the physical EFT amplitude, hence the unknown short-distance contribution. The practicality of the method, however, relies on the presence of only a finite (and few) number of on-shell intermediate states that are composed of no more than two hadrons, such that the Minkowski amplitude can be constructed. One can estimate the expected nature and the number of intermediate states by examining a plausible example. Let us take L = 8 fm to ensure the validity of the finite-volume formalism used with physical quark masses, up to exponentially suppressed contributions [70,71]. The finitevolume spectrum of the two-nucleon isotriplet channel at rest arising from singularities of the function in Eq. (11) can be determined using the experimentally known phase shifts [72], giving the ground-state energy E ni ≈ −2.6 MeV (which polynomially approaches zero as L → ∞). A simple kinematic can be considered for the transition amplitude such that E i (= E ni ) = E f (= E n f ), and where the currents carry zero energy and momentum so that the final-state two-nucleon system remains at rest. Given the available total energy, and the quantum numbers of the currents, the only allowed intermediate state is the twonucleon isotriplet channel at rest, whose low-lying spectrum in this volume is E * m ≈ {−5.6, 10.5, · · · } MeV. While it may appear that the ground state of this system constitutes an on-shell intermediate state, requiring construction of the Minkowski amplitude through an evaluation of the isosinglet to isotriplet matrix element, one must note that since the zero spatial momentum is not allowed for the neutrino propagation in the finite volume, none of the on-shell conditions in Eq. (19) can be satisfied with the kinematics considered (noting that the minimum allowed energy of an on-shell neutrino in this volume is |P * | = 2π/L ≈ 155.0 MeV). As a result, T L (τ ), and Eq. (14) can be readily used to obtain the physical amplitude from the lattice-QCD four-point function G (E) L (τ ). The example described demonstrates that obtaining the physical amplitude of the nn → pp (ee) process from lattice QCD is even more straightforward than its twoneutrino counterpart, as in the latter there is a larger kinematic phase space allowed for on-shell intermediate states-a feature that is shared with the nuclear structure calculations of the corresponding nuclear matrix elements [73]. Without experimental input, such nuclear many-body calculations, however, are unable to isolate and quantify the short-distance contribution to the matrix element, motivating the need for first-principles determination of the matrix elements in the few-nucleon sector, and a thorough matching to nuclear EFTs. The current framework, therefore, takes an essential step in enabling such a matching program in the upcoming years. Besides its application in the 0νββ process, the formalism outlined in this paper will find its use in a range of hadronic processes that consist of single-or two-hadron initial, intermediate, and final states, and where a light lepton (or photon) propagator is present, such as in the semi-leptonic rare decays of the kaon [74]. The steps involved in performing the two-loop sum-integral difference defined in Eq. (13) for an expedited convergence will be outlined in this section. To simplify the notation, the conventions n ≡ |n| and n 2 ≡ |n| 2 are used for any three-vector n.
The two-loop integral involving the neutrino propagator is given by This integral is divergent in the UV region of integrating variables and must be regulated. While for the discussion of the physical amplitude in the main text, dimensional regularization is a natural choice as presented in Eq. (9) [10], for the evaluation of the sum-integral difference, a cutoff regulator Λ proves most useful. As the physical amplitude, as well as the matching condition are UV convergent, both choices can be used in the matching framework. In particular, with the cutoff regularization, J ∞ evaluates to In a finite volume with cubic geometry and spatial extent L along each Cartesian coordinate and with periodic boundary conditions, the analog of Eq. (25) is given by replacing integrals with sums over quantized three-momenta k = 2πn/L with n ∈ Z 3 : Here, the i terms are dropped from the denominators since discrete sums are defined over non-singular values of k 1 and k 2 . Equation (27) differs from Eq. (25) by power-law correction in 1/L which can be isolated from the difference To evaluate δJ V , let us first convert the summation variable in Eq. (27) from k 1(2) to n 1(2) and rescale p 1(2) as p 1(2) = p 1(2) L/2π. Next, one can observe that the UV divergence in Eq. (27) is the same as that occurred in the sum when p 2 1 = p 2 2 = 0. Using a cutoff regulator Λ, this sum reads The upper bound on the sum over n 1 indicates that only integer triplets that satisfy n 1 ≤Λ (= ΛL/2π) must be included. The sum over n 2 is left unbounded. Now adding and subtracting J V (0, 0) and upon using Eq. (26), Eq. (28) becomes where To evaluate the lattice sums in Eqs. (31)- (34), one can use the method of tail-singularity separation (TSS) described in Ref. [75]. In this method, the sum is split into two pieces: one containing the singular contributions and the other containing a power-law tail which is sufficiently smooth such that it can be approximated by its integral counterpart. As an example, let us sketch out the details for evaluating the lattice sum in Eq. (33). The TSS scheme can be achieved by introducing exponential factors containing a small positive number, α, and rewriting X 3 as The smooth function containing a power-law tail is obtain by gathering 1 − e −αn 2 and 1 − e −α(n 2 −p 2 ) factors, which is then approximated, up to O(e −π 2 /α ) corrections, by n → d 3 n with values at the poles removed. This gives 2 1 − e −α(n 2 −p 2 ) n 2 (n 2 −p 2 ) − 2α p 2 sinh (αp 2 ) + O(e −π 2 /α ).
TSS can now be used on X 6 (n 1 ,p 2 ) and X Here, P denotes the Cauchy principal value of the radial integration in n 2 for the pole n 2 2 =p 2 . Using Eqs. (40) and (41) in Eq. (38), the outer sum over n 1 is then split into two parts: terms containing exponentially suppressed terms in α and terms independent of α. The former can be summed directly while the latter, which stems from evaluating the α-independent integrals in Eq. (40) and (41), can be calculated using TSS just like Eq. (36). As an example, with this procedure, Eq. (38) forp 2 = 1 evaluates to X 6 (1, 1) = 264, which agrees with Eq. (A10) of Ref. [63] up to three significant figures. Arbitrary accuracy can be achieved by decreasing the value of α and increasing the number of integer-triplets used in the sums with increasing magnitude.