Two-neutrino double-beta decay in pionless effective field theory from a Euclidean finite-volume correlation function

Two-neutrino double-beta decay of certain nuclear isotopes is one of the rarest Standard Model processes observed in nature. Its neutrinoless counterpart is an exotic lepton-number nonconserving process that is widely searched for to determine if the neutrinos are Majorana fermions. In order to connect the rate of these processes to the Standard Model and beyond the Standard Model interactions, it is essential that the corresponding nuclear matrix elements are constrained reliably from theory. Lattice quantum chromodynamics (LQCD) and low-energy effective field theories (EFTs) are expected to play an essential role in constraining the matrix element of the two-nucleon subprocess, which could in turn provide the input into ab initio nuclear-structure calculations in larger isotopes. Focusing on the two-neutrino process $nn \to pp \, (ee \bar{\nu}_e\bar{\nu}_e)$, the amplitude is constructed in this work in pionless EFT at next-to-leading order, demonstrating the emergence of a renormalization-scale independent amplitude and the absence of any new low-energy constant at this order beyond those present in the single-weak process. Most importantly, it is shown how a LQCD four-point correlation function in Euclidean and finite-volume spacetime can be used to constrain the Minkowski infinite-volume amplitude in the EFT. The same formalism is provided for the related single-weak process, which is an input to the double-$\beta$ decay formalism. The LQCD-EFT matching procedure outlined for the double-weak amplitude paves the road toward constraining the two-nucleon matrix element entering the neutrinoless double-beta decay amplitude with a light Majorana neutrino.


I. INTRODUCTION
In certain nuclear isotopes with an even number of neutrons and protons, the rate of single-β decay is suppressed compared with a double-β decay. Consequently, two of the neutrons in a parent nucleus can convert to two protons while emitting two electrons and two antineutrinos: (A, Z) → (A, Z + 2) + eeν eνe . Here, A and Z are the atomic and proton numbers of the parent nucleus, respectively. Since the theoretical prediction of this process and the first estimation of its rate in 1930s [1], such a decay has been observed conclusively in a dozen isotopes ranging from 48 Ca to 238 U [2,3]. This decay is one of the rarest Standard Model processes in nature, with experimental half-lives ranging from ∼ 10 19 to ∼ 10 24 years [2,3]. The exotic beyond the Standard Model counterpart of this process, namely the neutrinoless mode (A, Z) → (A, Z + 2) + ee, remains at the center of vigorous theoretical studies and experimental searches [4,5], as its observation will unveil the nature of neutrinos and provides a mechanism for lepton-number violation, hence baryon-number violation, in the universe.
The two-neutrino mode (2νββ), nonetheless, continues to gain much attention, both theoretically and experimentally, for several reasons. First, this decay mode is a dominant background for the much less probable neutrinoless process which could occur in the same isotopes. Therefore, accurate constraints on its decay rate, and on the spectral shape of electron energies emitted in the decay, are crucial for better understanding of the background in neutrinoless double-β (0νββ) decay searches [3,5,6]. Second, the measured decay rates can be converted to constraints on the corresponding nuclear matrix elements (MEs), under the assumption that only the Standard a davoudi@umd.edu b ksaurabh@umd.edu arXiv:2007.15542v1 [hep-lat] 30 Jul 2020 Model weak interactions are in play. Subsequently, these MEs can be compared against theoretical determinations to test the validity of the nuclear-structure models used [7][8][9][10]. This, along with theoretical calculations of the single-β decay MEs, can refine these models, and inform the calculations of the exotic 0νββ decay process, for which the role of theory is crucial in constraining the new-physics mechanisms underlying the decay [4,5]. Last but not least, the increased sensitivity of 0νββ decay searches in recent years has subsequently led to a number of high-statistics measurements of the two-neutrino mode [9,[11][12][13][14][15], opening a window for probing potential beyond the Standard Model scenarios through this decay mode as well (see e.g., Ref. [16] for constraints on right-handed currents from the energy distribution and angular correlation of outgoing electrons in the 2νββ decay). Such investigations require accurate nuclear MEs to be computed within the Standard Model, augmented with reliable uncertainties.
Computing nuclear MEs for the 2νββ process (and its neutrinoless counterpart) is a rather challenging task given the quantum many-body nature of the nuclear isotopes used, as well as uncertainties in nuclear interactions and currents entering the calculations [17,18]. This manuscript takes a bottom-up approach to this problem, laying out a possible framework for providing the Standard Model input for the ME involved at the microscopic level, namely that in the nn → pp (eeν eνe ) process, using a combined LQCD and EFT approach. While such a process will not occur unless embedded in certain nuclear backgrounds, an apparent hierarchy of nuclear interactions and nuclear responses to external probes at low energies indicates that the most significant contribution to the nuclear decay may arise from single and correlated two-nucleon couplings, and corrections to this picture can be evaluated within an EFT power counting systematically. The contributions to the ME at each order in the EFT depend upon the low-energy constants (LECs) for interactions and currents, which can ideally be constrained by matching to the Standard Model determination of the same ME using LQCD methodology. Setting up this step is at the core of the present work. Upon successful completion of LQCD determinations of the double-weak two-nucleon ME [19,20], the formalism presented can be employed to constrain EFTs that will be used in systematic ab initio calculations of the decay rate in relevant isotopes. The compatibility of the EFT description in the few-nucleon sector with ab initio nuclear many-body methods, and the validity of the adopted EFT power counting when applied to a large nuclear isotope, will require further study and are not the focus of the current work. 1 Pionless EFT [22][23][24][25], a nuclear EFT relevant for processes with intrinsic momenta well below the pion mass, has shown notable success in studies of light nuclei and their electromagnetic and weak response, both at physical and unphysical quark masses (the latter accessible from LQCD), see Ref. [26] for a recent review. It is known, in particular, that within the Kaplan-Savage-Wise power counting of the pionless and pionfull EFT [22,23], a better convergence is observed as the EFT order is increased for observables in the two-nucleon (NN) spin-singlet channel [27,28]. In a recent study, pionless EFT has further been applied to the problem of 0νββ decay in the two-nucleon system up to next-to-leading order (NLO), and is found to require a new short-distance coupling, or contact LEC, already at leading order (LO) to fully renormalize the two-nucleon ∆I = 2 transition amplitude, where I denotes the total isospin [29,30]. This feature is not expected to occur for the two-neutrino mode, since the 1/q 2 -type neutrino potential induced by the exchange of a light neutrino among the nucleons is absent in the two-neutrino process. Here, q denotes the momentum of the neutrino (i.e., antineutrino) exchanged. Nonetheless, there has been no explicit calculation of the nn → pp (eeν eνe ) amplitude in the pionless EFT with nucleonic degrees of freedom (DOF) to date.
The first attempt at constructing and constraining the nuclear ME for this process was presented in Refs. [19,20], in which a dibaryon formulation of the pionless EFT [31,32] was considered. In this case, a contact LEC, called h 2,S in this reference, was naturally introduced to account for the conversion of a spin-singlet dibaryon with I 3 = −1 to that with I 3 = 1 at tree level, similar to the role of the so-called l 1,A coupling of the dibaryon formalism that turns a spin-singlet dibaryon with I 3 = 0 to a spin-triplet dibaryon with S 3 = 0 at tree level. Here, I 3 (S 3 ) is the third component of the total isospin (spin) operator. Importantly, the first LQCD computation of the relevant ME at a pion mass of ≈ 800 MeV was performed at the same time, and led to a constraint on this new short-distance coupling, 2 demonstrating its roughly equal contribution to the ME when compared with the l 1,A effects in the deuteron-pole contribution. It is valuable to obtain an EFT amplitude within the nucleonic formulation, which is a more suitable framework in connecting to nuclearstructure calculations of the ME in larger isotopes. Such an amplitude is calculated in the current work in Sec. IV A for the 2νββ decay, and as a recap, in Sec. III A for the single-β decay process (for which the result had been known previously in the context of general single-weak processes in the two-nucleon systems, e.g., in Refs. [33][34][35]). It is shown that in both cases, the axial coupling of the nucleon and the correlated coupling of the two-nucleon system to an axial-vector current are sufficient to obtain a renormalization-scale independent amplitude at NLO. Such a conclusion for the power counting can be verified through simultaneous LQCD computations of both MEs.
Perhaps an even more significant component of this work is to provide the framework for obtaining the hadronic amplitude for the two-nucleon double-β decay from four-point correlation functions of two nucleons obtained with LQCD in a finite Euclidean spacetime. It is well known that multi-hadron amplitudes cannot be directly accessed from LQCD Euclidean correlation functions [36], and a mapping is required to constrain the amplitudes evaluated at eigenenergies of the system in a finite volume. These finite-volume energies, obtained from hadronic two-point functions, are sufficient for obtaining elastic and inelastic amplitudes with strong interactions in the two-and three-hadron sectors up to the next inelastic thresholds, as introduced in Refs. [37,38] and extended in Refs. . To constrain transition amplitudes involving external currents, finite-volume MEs, obtained from LQCD three-point functions, are additionally required. Such formalisms are developed in Refs. [45,[70][71][72][73][74][75][76][77]. For nonlocal MEs derived from four-point functions, such as those relevant for rare kaon-decays and Compton amplitudes [78][79][80][81][82][83][84][85][86][87][88], and hadronic double-β decays [19,20,[89][90][91][92], in addition to identifying the volume dependence of MEs, another crucial step is involved. Explicitly, due to the dependence of the time-ordered product of currents present in a four-point function on the spacetime signature, Euclidean and Minkowski MEs are fundamentally different when intermediate states with on-shell kinematics are present. These must be subtracted from the Euclidean correlation function, and be separately constructed from the knowledge of appropriate two-and three-point functions. These contributions must then be added back to the non-problematic contributions to obtain the complete Minkowski ME in a finite volume, which is used to construct the physical amplitude, see e.g., Ref. [88].
Various components of the formalisms mentioned above, such as identifying volume corrections to external and intermediate two-hadron states, as well as reconstructing the Minkowski amplitude from the Euclidean four-point function, are relevant for the mapping of the nucleonic ME in the 2νββ transition as well. Nonetheless, none of the work developed so far can be applied directly to the problem at hand, mainly due to new single-body contributions to the two-body nonlocal ME. 3 Furthermore, while the FV matching for the single-weak processes, such as pp fusion, has been presented in the past within the pionless-EFT framework [45], such a matching for the EFT amplitude of the 2νββ transition with nucleonic DOF is developed for the first time in the current work, as outlined in Secs. IV B, IV C, and IV D below. This parallels a related mapping with dibaryon formalism presented in Ref. [20] using a different approach. Finally, a new derivation is presented for the previously known result for the single-weak amplitude in the two-nucleon sector, combining the pionless-EFT result for the correlation function in a finite volume with the derivation of Ref. [77] for a general mapping for 2 → 2 single-weak processes, as outlined in Secs. III B and III C. The conclusion and outlook of this work are presented in Sec. V.

II. SETTING UP THE FORMALISM
In pionless EFT, the hadronic Lagrangian 4 is arranged according to the number of nucleons, where L (n) is comprised of n−nucleon operators, and ellipsis denotes higher-nucleon operators. The single-nucleon Lagrangian is the nonrelativistic kinetic-energy operator of the nucleon at LO, and ellipsis denotes relativistic corrections. Here, ∂ t is the time derivative and ∇ is the spatial gradient operator. N = (p, n) T 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. The Lagrangian with two-nucleon operators is given by where the overhead arrow indicates which nucleon fields are being acted by the derivative operator, and ellipsis denotes higher-derivative operators. Here and everywhere, repeated indices are summed over where i = 1, 2, 3. P i and P i are the spin-isospin projection operators for spin-singlet ( 1 S 0 ) and spin-triplet ( 3 S 1 ) channels, respectively, with definition and normalization where σ i (τ i ) are the Pauli matrices acting in spin (isospin) space. The term proportional to C 0 ( C 0 ) corresponds to the LO contact interaction in the 1 S 0 ( 3 S 1 ) channel while the terms proportional to C 2 ( C 2 ) describe the NLO momentum-dependent interactions. The couplings in the EFT are generally renormalization-scale dependent, and their Renormalization Group flow is determined by requiring the physical observables such as amplitudes to be scale independent order by order in the EFT. The renormalization-scale dependence of the couplings in Eq. (3) is given in Eq. (29) below. The effective Lagrangian for the charged-current (CC) weak interaction is given by where G F is the Fermi's constant. The leptonic current contains electron, e, and neutrino, ν, fields, and the hadronic current can be written in terms of vector and axial contributions as The superscript (subscript) in the hadronic (leptonic) current denotes isovector components while the subscript (superscript) denotes the spacetime vector components. The vector current mediates Fermi transitions, while the axial current governs Gamow-Teller transitions, which correspond to different isospin (∆I) and spin (∆S) selection rules. The focus of this work is on isovector transitions that lead to (double) β decays. Further, only low-energy processes are considered such that pionless EFT is a proper description. As a result, the only relevant currents are the spatial component of the axial-vector isovector currents. 5 In fact, it is the nuclear ME of such currents, i.e., Gamow-Teller MEs, that are not constrained precisely phenomenologically and their determination will be the subject of a LQCD-EFT matching program. At LO in pionless EFT, the nonrelativistic one-body operator with such quantum numbers is where τ + = (τ 1 + iτ 2 )/ √ 2 and g A is the nucleon axial charge. k = 1, 2, 3 is the spatial Lorentz index. The momentum-independent two-body axial-vector/isovector current is where P + = (P 1 + iP 2 )/ √ 2. This current contributes to the β-decay amplitude at NLO, as shown below, leading to a Renormalization-Group equation for L 1,A , obtained in terms of stronginteraction couplings in Eq. (3) as well as the nucleon's axial charge, see Eqs. (56) and (57). All fields introduced in this subsection have spatial and temporal dependence.

A. Infinite-volume amplitudes
The nuclear weak-decay processes involve nonperturbative strong interactions between nucleons, governed by the terms in Eq. (3), and perturbative weak interactions between nucleons and leptons from L CC in Eq. (5). The Feynman amplitude of a weak process involving n-lepton currents is obtained by considering the nth order term in the perturbative expansion of the S-matrix with respect to the CC interaction. This amplitude can be decomposed into leptonic and hadronic parts as where X denotes the transition considered, and all the amplitudes introduced have appropriate dependence on on-shell four-momenta of initial and final states that are suppressed for brevity. The leptonic amplitude, M lep. X , involves the kinematics of the electron(s) and neutrino(s) and the hadronic amplitude, M X , is the amplitude involving initial and final nucleonic states. The aim of this paper is to relate the MEs obtained in a finite Euclidean spacetime to this physical hadronic amplitude for the case of single-and double-β decays. As a result, we turn our focus on the amplitude M X and will not carry out the straightforward convolution 6 of the hadronic contribution with the leptonic amplitude. This convolution is more involved for the case of 0νββ decay process, which includes a Majorana neutrino propagator, and will be studied elsewhere [93]. Given its nonperturbative nature, the hadronic amplitude M X is obtained by solving the Schrödinger's equation with the strong interaction potentials that act between nucleons in initial, final, and when relevant, intermediate nuclear states. 7 The β-decay processes considered in this work involve two-nucleon initial and final states and Gamow-Teller transitions, and pionless EFT will be used to construct the amplitudes up to and including NLO. Thus, as mentioned above, the only hadronic currents that contribute are those given in Eqs. (8) and (9). The Hamiltonian for these processes is given byĤ whereĤ S is the strong Hamiltonian withĤ 0 andV S being the free two-nucleon Hamiltonian and the strong interaction potential, respectively. These are derived straightforwardly from the strong-interaction Lagrangian density in pionless EFT presented in Eq. (1). The weak-Hamiltonian operatorĤ CC can be understood as the product of a leptonic part and a hadronic part, and the latter is the only relevant part when acting on the state of two nucleons. In particular, assuming a final-state kinematic such that each hadronic current transfers zero spatial momenta (corresponding to pairs of electron and neutrino in a n-β decay carrying away no spatial momenta), and at the EFT order considered in this work, the hadronic part,Â CC , can be written aŝ whereÂ are constructed from the one-body and two-body current operators, respectively. The one-body operator is extended to two-nucleon Hilbert space by taking a tensor product (⊗) 8 with the identity operator, 1, in the nonparticipating nucleon space. For a given N N → N N process, Lorentz index k determines the change in spin along the spin quantization axis. In order to treat strong interactions nonperturbatively, one needs to find the Feynman amplitude between two-nucleon eigenstates ofĤ S , which are constructed from the eigenstates ofĤ 0 [94]. In the center-of-mass (CM) frame of the two-nucleon system, theĤ 0 eigenstate in channel N N ( 1 S 0 or 3 S 1 ) with relative momentum p and energy E is given bŷ wherep is the relative-momentum operator in the CM frame, and p 2 = p · p. This state is normalized according to consistent with the normalization convention for nonrelativistic two-body states with zero total momentum. Using the free retarded Green's function with > 0, and the T-matrix for strong interaction theĤ S eigenstate in two-nucleon channel N N with energy E and relative momentum p is given by The hadronic Feynman amplitude, M X , between the initial and final two-nucleon states with CM energies E i and E f , respectively, is then given by where E − f is defined through an advanced Green's function, i.e., i → −i in Eq. (16), and J X is the hadronic part of the T-matrix for the weak interaction. For the nth order weak process, J X is constructed using n insertions ofÂ CC and strong retarded Green's function, G S (E + ), defined as The hadronic part of the T-matrix corresponding to the first-and second-order weak processes of this paper will be introduced in Eqs. (51) and (70), respectively.
Equations (17), (18), and (19) imply that the hadronic Feynman amplitude can be expressed in terms of MEs ofV S andÂ CC between free eigenstates, defined in Eq. (14). For a given two-nucleon channel, these states can be constructed using the spin-isospin projection operators in Eq. (4). From the Lagrangian in Eq. (3), the MEs ofV S between such momentum eigenstates are given by for LO and NLO interactions, respectively. Similar relations can be written for the spin-triplet channel with the corresponding couplings. One can similarly obtain the MEs of the one-body and two-body weak hadronic Hamiltonian using Eq. (8), (9), and (13), which are only nonvanishing when taken between the spin-singlet and spin-triplet states: Using Eqs. (21) and (22) and inserting a complete set of free single-and multi-particle states in the iterative sum in Eq. (17), one can obtain the T S MEs up to NLO: Similar expressions can be formed for the spin-triplet channel upon replacements 1 S 0 → 3 S 1 , T S → T S , C 0 → C 0 , and C 2 → C 2 . I 0 is an ultraviolet (UV) divergent integral that is regularized with the power-divergence subtraction (PDS) scheme introduced in Ref. [22], where p = √ M E, and µ is the renormalization scale. In a Feynman-diagram expansion, this divergence corresponds to the two-nucleon s-channel loop, see Fig. 2.
The on-shell N N → N N elastic scattering amplitude is defined as with a similar expression for the amplitude in the spin-triplet channel upon replacements 1 S 0 → 3 S 1 , M → M, and T S → T S . From the requirement of renormalization-scale invariance of this amplitude, the scale dependence of the strong interaction couplings can be deduced: with similar relations for the spin-triplet channel upon replacements C 0 → C 0 , and C 2 → C 2 .
In the remainder of this paper, two simplifying conventions will be applied to notation. First, the N N channel index will be dropped from the states, and their assignment must be deduced from the ME under consideration. Instead, when helpful given the context, the kinematic dependence of the states will be displayed. Second, the energies are meant to be arising from a retarded propagator unless specified, hence the superscript '+' on energies will be dropped for brevity.

B. Finite-volume methodology
In this section, the well-known formalism for relating two-body elastic scattering amplitudes will be reviewed, setting the stage for the derivation of the matching relation between finitevolume three-and four-point functions in the two-nucleon systems and the corresponding hadronic transition amplitudes. The derivation presented, drawn primarily from the approach of Refs. [40,76,77,88], provides a suitable ground for the finite-volume analysis of later sections. In all finitevolume quantities considered, the temporal extent of spacetime is taken to be infinite, while the spatial extents are taken to be finite. The spatial geometry is taken to be cubic with extents L, and periodic boundary conditions are assumed on the fields. As a result, while the energy remains a continuous variable, each Cartesian component of the spatial momentum of the non-interacting nucleons is discretized in units of 2π/L.
In comparing assumptions of this work with the conditions on LQCD correlation functions, several differences can be highlighted. First, the spacetime is discretized in LQCD computations. In the following, therefore, it is assumed that the continuum limit of LQCD quantities are taken before matching to physical amplitudes through the matching relations provided in this work. Second, unless specified, the matching relations do not make any reference to the Euclidean spacetime and correlation functions/MEs are represented in Minkowski spacetime in both finite and infinite volume. Nonetheless, quantities computed with LQCD correspond to Euclidean spacetime. The distinction between Minkowski and Euclidean correlation functions becomes important only for the four-point functions, where the currents are separated in time, and the analytic structure of the MEs may differ significantly between different time signatures [88]. As a result, despite the twoand three-point functions, a straightforward analytic continuation does not transform quantities from one time signature to the other. Consequently, extra steps must be taken to use the matching relations of this work to extract physical amplitudes from LQCD four-point functions. These steps are outlined in Sec. IV D, following and extending the protocol of Ref. [88] for nonlocal MEs of single-hadron states.
Consider the two-nucleon systems in the spin-singlet channel. Similar relations can be obtained identically for the spin-triplet, or any two-hadron channels for that matter, below the lowest three-hadron inelastic threshold. Nonrelativistic kinematic is assumed throughout although generalization to relativistic kinematic is straightforward and known [23]. The relations presented are general for all partial-wave amplitudes (that besides possible mixing of the physical amplitudes in the spin-singlet channel, get further mixed in a finite cubic volume). However, the s-wave limit will be made explicit when necessary, which is relevant for two-nucleon processes of this work that occur at low energies. The key idea of relating the finite-volume MEs to infinite-volume amplitudes is to equate two equivalent definitions of the n-point correlation function in a finite volume [76,77], as will become evident shortly.
Consider the momentum-space finite-volume two-point correlation function of two nucleons projected to zero spatial momentum, FIG. 3. Diagrammatic representation of the finite-volume two-point correlation function corresponding to the expansion in Eq. (35). The black dots are the interpolating operators for the two-nucleon state in the spin-singlet channel. All other components are introduced in Fig. 1. The s-channel loops are evaluated as a sum over discretized momenta, as discussed in the texts.
which can be expressed, alternatively, in the mixed momentum-time representation: In these equations, B † and B are two-nucleon source and sink interpolating operators for the spinsinglet channel, respectively. The operators in the second line have spacetime dependences and x ≡ (x 0 , x) and y ≡ (y 0 , y). P denotes the total four-momentum of the two-nucleon system and P = (E, P = 0). The Minkowski metric g µν = diag(1, −1, −1, −1) is assumed throughout. T denotes the Minkowski time-ordering operator, and without loss of generality, it is assumed that that y 0 > x 0 . The subscript L on the spatial integral denotes the integral is performed over a finite cubic volume, and the [· · · ] L is used to denote the finite-volume nature of the ME enclosed. In going from the second to the third line, the relation B(x) = e iP 0 x 0 −iP ·x B(0)e −iP 0 x 0 +iP ·x is used, whereP 0 andP are, respectively, energy and momentum operators acting on the adjacent states. Furthermore, a complete set of discretized finite-volume two-nucleon states are inserted, with n denoting the state index. These finite-volume states in the CM frame are characterized by the total energy of two nucleons, E n , and an L argument to remind the finite-volume nature of these states. They are chosen to satisfy the normalization condition to be compared with the normalization of the infinite-volume states in Eq. (15). Note that for on-shell states, E n = p 2 n /M , where p n ≡ |p n | is the corresponding interacting momentum of each nucleon in the CM frame.
The goal now is to express C L (P ) in terms of a diagrammatic representation illustrated in Fig. 3, with the building blocks provided in Fig. 1. One can then perform the inverse Fourier integral in Eq. (31) to arrive at an equivalent form that can be compared against Eq. (33). The correlation function C L (P ) contains an infinite coupled sum with s-channel two-nucleon loops, iI V 0 , as depicted in Fig. 1-g, connected via the two-nucleon Bethe-Salpeter kernels, K, as shown in Fig. 1-a for pionless EFT. This expansion can be written in the compact notation of Refs. [69,88] as Here, all kinematic dependence of the functions are suppressed, but these will be made explicit shortly.B † andB are the zero-spatial-momentum projected source and sink interpolating operators for the two-nucleon system in the The symbol ⊗ should be interpreted as an sum/integral convolution of the adjacent functions. Explicitly, for two arbitrary functions χ and ξ that are smooth in q = (q 0 , q): Since χ and ξ are smooth functions, which is the case for interpolating fields and the kernels in Eq. (35) below inelastic thresholds, the only difference between the sum and the corresponding infinite-volume integral (beyond exponentially suppressed corrections) arise from singularities of 1/(E − q 2 M ) function in q. Note that i is now redundant since the sum runs over discrete values of momenta. Subsequently, Eq. (37) can be written as where the first term, χ ⊗ I 0 ⊗ ξ, is defined in the same way as the right-hand side of Eq. (36) except for the replacement 1 On the other hand, the second term, χ F 0 ξ, can be interpreted as the ordinary product of three functions, that are in general matrices in the angular momentum basis. In this basis, MEs of functions should be understood as projection of the original functions to a given partial-wave component with respect to the angular variable defined by q. Furthermore, the F 0 function, defined as the difference between the sum and integral, projects the functions adjacent to it to on-shell values of momentum, corresponding to |q| → √ M E ≡ p. For the s-wave projection that is a good approximation at low-energies, the finite-volume function F 0 can be written in a simple form where the sum-integral difference notation is defined as Furthermore, c lm is a function closely related to Lüscher's Z function: where n is a Cartesian vector with integer components, and n denotes direction of n. In the remainder of this paper, the product of functions can be interpreted either as an ordinary product of scalar quantities in the s-wave approximation, or a matrix product of quantities expressed in the angular-momentum basis, in which case the generalized form of F 0 should be used, see Ref. [37,38,40]. For the spin-triplet channel, the spin quantum number must be taken into account, giving rise to the finite-volume function derived in Ref. [47]. The s-wave approximation of this function, nonetheless, is the same at that presented in Eq. (39). Expressions for F 0 generalized to moving frames, hadrons with unequal masses, systems with relativistic kinematic, volumes with elongated sides, and general twisted BCs exist, see Refs. [55] for a review. 9 Now using Eq. (38) to expand and rearrange Eq. (35), it is straightforward to see that C L (P ) can be written as The first term is obtained by collecting all terms with I 0 in Eq. (35) to produces the infinite-volume two-point correlation function. All the functions in the second term are evaluated at the on-shell values of energy, i.e., they are solely a function of only one kinematic variable, E. B and B † are the endcap functions which are built out of infinite-volume quantities. Explicitly, iM is the full on-shell N N → N N scattering amplitude in the spin-singlet channel, as shown in Fig. 2, and is given by Given that the product of the functions in the second term of Eq. (42) is now an ordinary (matrix) product, these terms can be summed to all orders in F 0 : where a new finite-volume function F is defined for convenience in later discussions: Equation (45) has only singularities at interacting energies of two-nucleon system in the finite volume [76], E n . These arise from Lüscher's 'quantization condition': For a generic case with no s-wave approximation, the condition reads det F −1 0 + M = 0, where determinant is taken in the angular-momentum space. In any case, a cutoff is required on the partial waves included, as otherwise the relation is not of practical use. In summary, Lüscher's quantization condition provides a constraint on the physical elastic amplitude of two nucleons at the finite-volume eigenenergies of the two-nucleon system, which are accessible from LQCD computations of Euclidean two-point correlation functions in a finite volume. The condition is valid up to exponential corrections that go as e −L/R , where R ∼ M −1 π is the range of strong interactions and M π denotes the mass of the pion.
Having arrived at the form in Eq. (45) for C L (P ), one can now perform integration over E in Eq. (31) to obtain where the generalized Lellouch-Lüscher residue matrix is given by: Note that in order to arrive at this result in Minkowski space, the i prescription in the correlation functions must be recovered, i.e., E → E + i . Now comparing Eq. (48) with Eq. (33) for each value of E n , one obtains the following matching condition: This relation, and its spin-singlet counterpart, will be used in later sections to simplify the matching condition for three-and four-point functions obtained from the same source and sink two-nucleon interpolating operators.

III. SINGLE-β DECAY PROCESS
Low-energy single-weak processes in the two-nucleon sector, including pp fusion: pp → de + ν e , neutrino(antineutrino)-induced disintegration of the deuteron: ν(ν)d → npν(ν) andνd → e + nn, and muon capture on the deuteron: µd → nnν µ have been studied in the past within the framework of pionless EFT [33][34][35][95][96][97]. The Gamow-Teller ME in these processes is dominated by the single-nucleon contribution. At LO, the ME is characterized by the nucleon's axial charge, g A . The two-nucleon contribution, characterized by the L 1,A LEC in Eq. (9), is only a few percents of the full ME. It, nonetheless, constitutes the dominant source of theoretical uncertainty in the determinations of relevant cross sections [98]. A precise LQCD determination of the L 1,A LEC at the physical values of the quark masses will be a critical goal of the next-generation nuclear LQCD studies in the coming years. The two-nucleon contribution to Gamow-Teller ME can also be constrained from the half-life of the tritium [99,100]. With the finite-volume LQCD technology developed in this work, simultaneous fits to single-and double-β decay processes in the twonucleon sector will be possible, which could provide better constraints on this unknown LEC, see also Refs. [19,20,101]. Alternatively, the constraint on L 1,A from single-weak processes could be used to evaluate the significance of the higher-order two-weak two-nucleon contribution in the double-β decay, hence testing the validity of the EFT power counting of this work. Although the single-weak process nn → np( 3 S 1 ) does not occur in free space, its Gamow-Teller ME in the isospinsymmetric limit is the same as that in all the single-weak processes mentioned above. Since this process constitutes a subprocess of the double-β decay of the two-neutron system, the subscript nn → np is adopted for the amplitudes and correlation functions below, but the formalism is general for any 3 S 1 → 1 S 0 or time-reversed transitions in two-nucleon systems. Since these will be ingredients to the double-β decay formalism, single-weak amplitudes and their matching to finite-volume correlation functions will be studied in detail in this section. The derivation of the matching relation is new and follows that of Ref. [77] for general 2(J) → 2 processes. The result obtained is in agreement with the earlier results on this problem presented in Refs. [45,71].
A. Physical single-β decay amplitude Following the formalism presented in Sec. II A, the single-β decay amplitude will be constructed in this section at NLO in pionless EFT. The hadronic transition considered in Eq. (19) is X ≡ nn → np. For simplicity, we consider the initial two-nucleon state to be at rest with energy E i and further set the four-momentum carried away by the weak current to (E, 0), so the two-nucleon system in the final state remains unboosted with energy E f , where E f = E i − E. The initial neutronneutron state is a spin-singlet state while the final neutron-proton state is a spin-triplet state with m = 0, ±1, where m is the eigenvalue of σ 3 . This state arise from the (k =)mth component of the FIG. 4. Diagrams contributing to the single-β decay hadronic amplitude in Eq. (52) at NLO in pionless EFT. All building blocks of this amplitude are introduced in Fig. 1. The ( ) symbol under a given diagram indicates that a counterpart of the diagram must be also considered in which the diagram is reversed and the replacement 1 S 0 ↔ 3 S 1 is applied to all components except the initial and final two-nucleon states. For an on-shell amplitude, the external legs must be evaluated at on-shell kinematics. currents in Eqs. (8) and (9). Given the rotational symmetry, transition amplitudes with different m values are equal. Thus, we fix m = 0 for the final state and omit the m index from the notation.
The hadronic part of the T-matrix for the nn → np transition, a first-order process in the weak Hamiltonian, is given by whereÂ CC is given in Eq. (12). Substituting this in Eq. (19) and using Eq. (18), one can decompose the amplitude into different terms depending on the number of T S insertions. Diagrammatically, the contributions to the hadronic amplitude are shown in Fig. 4. This amplitude can be computed up to NLO in strong interactions by inserting a complete set of hadronic states between the operators and using Eqs. (16), (23), (24), (25), and (26). A straightforward calculation gives where p i(f ) = M E i(f ) p i(f ) , with the wide hat used to denote the directionality of the threevector. The amplitude is conveniently split into three terms, the free-neutron decay contribution, a divergence-free (DF) term and a divergent term, with divergence refers to the behavior of the terms in the E i → E f limit. Only the second term, the divergence-free amplitude, is relevant for the finite-volume calculation of this process, as shown in Sec. III B. It is given by Here, I 1 is a loop introduced in Fig. 1-i with three propagators, and is given by Furthermore, in arriving at Eq. (53), the identity is used. The L 1,A is a renormalization-scale independent combination of strong and weak couplings, 10 This is because, first of all, from Eqs. (54) and (27), I 1 (E i , E f ) can be seen to be µ independent. Secondly, divergent terms in Eq. (52) are also renormalization-scale independent, as the g A coupling and the N N → N N elastic scattering amplitudes are scale independent. Thus, for the physical amplitude in Eq. (52) to be independent of the renormalization scale, L 1,A should satisfy This, along with Eq. (29) and its spin-triplet counterpart, determines the µ dependence of L 1,A . This result was first obtained in Refs. [34,35].

B. Finite-volume correlation function
The first step in establishing a matching relation for the infinite-volume single-β decay amplitude calculated in Sec. III A is to construct the momentum-space three-point function in a finite volume, C (J ) L (P i , P f ). P i(f ) denotes the four-momentum of the two-nucleon state in the spin-singlet (triplet) channel at rest with P i(f ) = (E i(f ) , P i(f ) = 0). This correlation function can be expanded diagrammatically as shown in Fig. 5. The result can be expressed as Here,B † ( B † ) andB( B ) are two-nucleon source and sink interpolating operators for the spin-singlet (triplet) channel, respectively, all projected to zero spatial momentum. I V 0 , K , K, as well as ⊗ product notation, are defined in Sec. II B. I V 1 is the finite-volume counterpart of the I 1 loop defined in Sec. III A, i.e., the replacement d 3 q (2π) 3 → 1 L 3 q must be applied. It is important to expand the correlation function in Eq. (58) in fixed orders in the EFT to ensure the renormalization-scale independence of physical observables extracted from the correlation function, including the energy eigenvalues and the finite-volume MEs. For convenience, this simple form will be carried out for now but a pionless EFT expansion at NLO will be performed shortly.
The goal is to separate the infinite-volume correlation function and to re-arrange purely finitevolume contributions in a manner similar to that outlined in Sec. II B. Firstly as was shown in Eq. (38), the sum convolution of functions adjacent to I V 0 can be separated into two contributions, the infinite-volume piece and a remnant finite-volume piece in which ⊗ products turn into ordinary (matrix) products of functions, and where the F 0 function puts any adjacent functions on shell. The sum convolution of functions adjacent to I V 1 for two generic left and right functions χ and ξ, defined as FIG. 5. Diagrammatic representation of the finite-volume correlation function with a single insertion of the weak current corresponding to the expansion in Eq. (58). The black dots are the interpolating operators for the two-nucleon channels in spin-singlet or triplet states. All other components are introduced in Fig. 1. The loops are evaluated as a sum over discretized momenta, as discussed in the texts.
can proceed similarly, except now the left and right sides of I V 1 have, in general, different on-shell kinematic. The general case of the convolution of functions with arbitrary momentum dependence adjacent to a (relativistic) I V 1 is worked out in Ref. [77]. Here, the situation is simpler, as at NLO in the pionless-EFT expansion of the three-point function, there are only two relevant scenarios to consider. First, one encounters two contact (momentum-independent) kernels or current on both sides of the I V 1 loop, in which case the convolution becomes trivially the ordinary product of (on-shell) functions. The function I V 1 in such ordinary products is that defined in Eq. (54) upon replacements where I 0 and F 0 are defined in Eqs. (27) and (39), respectively. It is evident that I V 1 is UV convergent and is hence scale independent. Note also that this function is regular in the limit E i → E f . The other possibility for the convolution of I V 1 with two adjacent functions at NLO is when there is one momentumdependent kernel C 2 q 2 (or C 2 q 2 ) on either side of I V 1 (with q being the summed momentum), in which case the application of the finite-volume counterpart of the identity in Eq. (55) leads to ordinary product of on-shell functions.
With these considerations, the expansion in Eq. (58) can be organized as follows. First, collecting all terms with only infinite-volume loop contributions isolates the infinite-volume three-point function. Second, one can collect all terms that contain at least one F 0 in the ∞ n=0 ⊗ K ⊗ I V Here, ellipsis denotes all other terms in Eq. (58) that neither contribute to C (J ) ∞ nor to the finitevolume terms in which factors of F and F are both present. The reason for such a rearrangement of the terms in the finite-volume correlation function will become evident shortly. M DF,V nn→np is the finite-volume counterpart of the single-β decay divergence-free amplitude M DF nn→np defined in Eq. (53), in which I 1 is replaced with Here, F 1 = I V 1 − I 1 , or in terms of the sum-integral notation defined in Eq. (40), 11 Note that we have overloaded the I V 1 symbol, used both as a sum with an arbitrary summand (made of left and right functions) times the corresponding propagators, and as the sum evaluated with just the propagators. When I V 1 is not convoluted with adjacent functions with the ⊗ product sign, it is meant to be a function on its own as defined in this paragraph.
where F 0 is defined in Eq. (39). Finally, it should be noted that all terms in Eq. (60) are evaluated at the on-shell kinematic, given their proximity to F 0 functions, and given the convolution rules with the I V 1 loop explained above. Eq. (60) will be used in the next subsection to derive the matching relation between the finite-volume ME and the physical hadronic amplitude for the single-β decay process.

C. The matching relation
To obtain the infinite-volume hadronic transition amplitude for the nn → np process from the corresponding ME in a finite volume, once again one could inspect the finite-volume correlation function. In momentum space, the correlation function is given by: while the momentum-time representation of this correlation functions can be written in different ways: In the last equality y 0 > z 0 > x 0 is assumed, complete sets of finite-volume states are inserted between the operators, and the operators in the Heisenberg picture are expressed in the Schrödinger picture. Now Eq. (60) can be used in Eq. (64) to perform the energy integrations. This gives Here, the only contributions to the energy Fourier integrals arise from the poles of the F and F functions, E n f and E n i , which are the finite-volume energy eigenvalues of the spin-triplet and spinsinglet two-nucleon systems at rest, respectively, see Eq. (47). R ( R) is the residue of F ( F) at the corresponding finite-volume energies, as defined in Eq. (49). This step makes it clear why only the contributions with both factors of F and F were collected in the correlation function explicitly, as the inverse Fourier integrals of the remaining contributions vanish. For a comprehensive proof of the absence of any poles beside the finite-volume two-nucleon energies in the finite-volume correlation function, see Ref. [77].
Having arrived at Eq. (67), it can be compared with its equivalent form in Eq. (66). For each E n i and E n f , one obtains: Multiplying this equation by its complex conjugate and using Eq. (50), one finally arrives at the relation that connects the hadronic ME in a finite volume with the divergence-free part of the physical amplitude: This equation is the main result of this section. Note that the divergence-free part of the physical hadronic amplitude is embedded in M DF,V nn→np , see Eq. (61). To construct the full hadronic transition amplitude, the divergent parts must be added to M DF nn→np . The divergent amplitudes, nonetheless, are comprised of contributions that can be determined from the nucleon single-β decay amplitude (g A coupling), as well as the strong-interaction two-nucleon scattering amplitudes. This relation, therefore, provides a means to constrain the L 1,A LEC from a LQCD calculation of the finite-volume ME of the axial-vector current in the two-nucleon system, see Ref. [101] for a first calculation along this direction.

IV. DOUBLE-β DECAY PROCESS
The analysis of the single-β decay in infinite and finite volume can now be extended to the case of the double-β decay. While the overall strategy is the same as that outlined in Sec. III, the doubleβ decay problem presents new features originating from the presence of two spacetime-displaced axial-vector currents in the MEs. This, first of all, makes distinct contributions arising from different time ordering of the currents for general kinematics. Second, and more significantly, it necessitates further steps to be taken to connect the Minkowski finite-volume ME to that calculated in Euclidean spacetime, which is the case with LQCD computations. The first two subsections follow the procedure of Sec. III, while the last subsection addresses this latter point. It is found that no new contact two-nucleon two-current LEC is needed to renormalize the physical hadronic amplitude at NLO, therefore the only two-nucleon axial-current LEC that is to be constrained from the matching relation is L 1,A . LQCD computations of this ME will be able to test the validity of this power counting.

A. Physical double-β decay amplitude
In this section, the hadronic amplitude for the double-β decay will be calculated. This is the amplitude given in Eq. (19) for X = nn → pp. The two currents are the same for this transition, and in the spin-isospin symmetric limit, all three spin components of the intermediate spin-triplet state contribute equally. In the following, only the amplitude for transitions with the k = 3 component of the currents in Eq. (13) is considered, and in the full amplitude this contribution must be multiplied by a factor of 3. Once again, the currents are assumed to carry out no spatial momenta and their kinematics is given by Q 1 = (E 1 , 0) and Q 2 = (E 2 , 0). The initial and final two-nucleon states are at rest, which constrains the intermediate states between the current to have zero spatial momentum. The hadronic contribution to the T-matrix is given by: Here,Â CC is the hadronic contribution to the weak Hamiltonian defined in Eq. (12), and G S is the strong retarded Green's function defined in Eq. (20). There will be two possibilities for the intermediate states, one with total energy E * 1 = E i − E 1 and one with total energy E * 2 = E i − E 2 .
FIG. 6. Diagrammatic representation of the hadronic transition amplitude contributing to the nn → pp process. All components of the diagrams are introduced in Fig. 1, where the external legs are evaluated at the on-shell kinematics. The ( ) symbol under a given diagram indicates that the reversed diagrams must be included as well (without changing the initial and final states). An exact set of diagrams must be computed but with the first and the second currents reversed in order, which is equivalent to setting E 1 ↔ E 2 in these diagrams, where E 1 and E 2 are the energies carried out by each of the two currents, see the discussions in the text.
These should be considered separately in the time-independent Lippmann-Schwinger formalism that is adopted for computing the amplitudes. The two contributions correspond to different time orderings of the currents in a time-dependent quantum field theory approach that is used in the definition of the finite-volume correlation functions. J nn→pp can now be used in Eq. (19) to obtain the hadronic amplitude up to NLO in pionless EFT. The diagrammatic expansion of this amplitude is shown in Fig. 6. Following the method described in Sec. II A, one obtains indicates that the counterpart of every term where E 1 is exchanged with E 2 must be included as well. Similar to Eq. (52), terms in Eq. (71) have been grouped into a divergence-free nn → pp amplitude, as well as divergent amplitudes which are singular in the limit E 1 , E 2 → 0. The divergent part contains the divergence-free single-weak transition amplitude, which is given by Eq. (53), 12 as well as N N → N N elastic scattering amplitudes for both spin-singlet and spin-triplet FIG. 7. Diagrammatic representation of the finite-volume correlation function with two insertions of the weak current corresponding to the expansion in Eq. (75). The black dots are the interpolating operators for the spin-singlet two-nucleon states. All other components are introduced in Fig. 1. The loops are evaluated as a sum over discretized momenta, as discussed in the texts. An exact set of diagrams must be computed but with the first and the second currents reversed in order, which is equivalent to setting E 1 ↔ E 2 in these diagrams, where E 1 and E 2 are the energies carried out by each of the two currents, see the discussions in the text.
channels. The divergence-free part of the nn → pp amplitude, M DF nn→pp , is given by where I 1 is give by Eq. (54) and the new type of loop arising from four propagators and two weak currents, as shown in Fig. 1-j, is given by and a useful identity is used to simplify the expressions: As the only new ingredient of the double-β decay amplitude compared with the single-β decay is M DF nn→pp , it is this quantity that is aimed to be constrained from a finite-volume matching relation in the next subsections. At the order considered in the EFT, the physical hadronic amplitude in Eq. (71) is evidently renormalization-scale independent (note that I 2 loop is UV convergent), and as a result, no new LEC beyond those present in the single-β decay amplitude is needed to renormalize the amplitude, as stressed before. A contact two-weak two-nucleon operator, in particular, should appear at next-to-NLO using a naive-dimensional analysis, and no new UV divergence at the previous order is observed to require its promotion to a lower order.

B. Finite-volume correlation function
The relevant finite-volume correlation function for the double-β decay process in momentum space can be represented by diagrams shown in Fig. 7, where an expansion at NLO in the pionless EFT can be performed later. These diagrams correspond to the following expansion Here, the kinematic dependence of the functions is suppressed, except for the dependence of the interpolating functions on total initial and final momenta. Ellipsis indicates that a similar expansion with E 1 ↔ E 2 must be included. All ingredients in Eq. (75) are introduced in the previous sections, except for ⊗I V 2 ⊗, where I V 2 is the finite-volume counterpart of the loop function with two insertions of the single-body weak current. For generic left and right functions χ and ξ, ⊗-product sign is defined as and similarly for E * 2 . For the analysis of this work at NLO in the EFT, the left and right convoluting functions are either momentum independent kernels/current, in which case the ⊗ sign trivially becomes an ordinary product, or one of the convoluting functions is proportional to q 2 and one is constant, in which case the finite-volume counterpart of the identity in Eq. (74) (in which q , I 1 → I V 1 , and I 2 → I V 2 ) can be used to return to the ordinary product of functions. When appearing in an ordinary product, I V 2 is simply given by Eq. (73) upon replacements just mentioned. Note that this function is regular in the limit E i → E f .
The expansion in Eq. (75) can now be turned into a useful form, in a similar fashion to Sec. III B. The idea is to first isolate the infinite-volume contribution, which arises from separating the loop functions to an infinite-volume integral and a sum-integral difference, and collecting all terms with exclusively infinite-volume loops. Next, given the identity in Eq. (38), factors of F(E f ) and F(E i ) can be formed out of ∞ n=0 ⊗ K ⊗ I V 0 n adjacent to B and B † , up to additional contributions that do not concern us. The F functions enforce on-shell kinematic on the adjacent function, which combined with identities on the ⊗I V 1 ⊗ and ⊗I V 2 ⊗ as described before, ensure that all functions in such a term are evaluated on-shell. These features, along with straightforward algebra, lead to where now ellipsis denotes any other terms left out while expanding Eq. (75) in terms of infiniteand finite-volume contributions. The reason behind this peculiar rearrangement of the terms, as discussed in Sec. III C, is to only keep contributions with poles in E i and E f (which are the finite-volume energy eigenvalues of two nucleons in the spin-singlet channel). This fact will be used in the matching relation that will be derived shortly. M DF,V nn→np in Eq. (77) is related to the physical divergence-free transition amplitude, M DF nn→np , which are defined in Eqs. (61) and (53) Here, F 2 = I V 2 − I 2 , or in terms of the sum-integral notation defined in Eq. (40), where F 1 is defined in Eq. (62). Equation (78) provides a direct relation between the finite-volume quantities and the physical divergence-free double-β decay amplitude.

C. The matching relation
The momentum-space four-point correlation function in the finite volume that was constructed in the previous subsection can be written as The time-momentum representation of the same correlation function is written in different forms following similar steps as in the previous sections. In the last equality, it is it assumed that y 0 > z 0 > 0 > x 0 or y 0 > 0 > z 0 > x 0 . The result in Eq. (77) for C can now be input in Eq. (81). Integration over E i and E f leads to nonvanishing contributions from the only poles of the function, which are the finite-volume eigenenergies in the spin-singlet channel, E n i and E n f , respectively. Note that these are solutions to quantization condition in Eq. (47). This gives rise to Here, R is the residue of the finite-volume function F evaluated at the corresponding finite-volume energy, as defined in Eq. (49). Equating this for each E n i and E n f with Eq. (83) and using Eq. (50), one finally obtains the desired matching relation for the double-β decay amplitude: with the finite-volume divergence-free amplitudes for single-and double-β decays defined in Eqs. (61) and (78), respectively. This relation is a main result of this work as it connects the (Minkowski) finite-volume ME of two time-ordered currents to the physical (Minkowski) amplitudes in infinite volume. The only subtlety is that LQCD calculations do not directly have access to the former, i.e., the left-hand side of this equation. Instead, they provide the Euclidean counterpart of this ME. It is, therefore, necessary to describe how to construct the Minkowski finite-volume ME from its Euclidean counterpart. Once this construction is carried out, Eq. (85) enables constraining the physical double-β decay amplitude.

D. Constructing the Minkowski correlation function from its Euclidean counterpart
The quantity derived in Eq. (85) in connection to the physical hadronic amplitude for double-β decay is defined in Minkowski spacetime. In particular, the ME on the left-hand side of Eq. (85) has the following spectral representation: where the kinematic dependence of T (M) L on the left-hand side is left implicit. m is an index corresponding to the tower of finite-volume energy eigenstates with quantum numbers of J nn→np |E i , L (or identically E i , L|J nn→np ). In contrary, such a spectral decomposition is ill-defined in Euclidean spacetime as integration over the Euclidean time τ ≡ z Extending the procedure of Ref. [88] to nonlocal two-nucleon MEs of this work, one can start by defining the time-momentum representation of the Euclidean correlation function: which can be directly accessed via a LQCD computation. The initial and final states are at zero spatial momentum and hence the p i and p f dependence of the function is dropped. T (E) means time ordering must be implemented with respect to the Euclidean-time variable. Moreover, the Heisenberg-picture operator in Euclidean spacetime satisfies J pp→np (0) e −P 0 τ +iP ·z , whereP 0 andP are energy (Hamiltonian) and momentum operators, respectively. In the following, without loss of generality we assume that J (E) pp→np (0) is the same as its Minkowski counterpart, i.e., it has no phase relative to the Minkowski current at origin. We thus drop the Euclidean superscript on the Schrödinger-picture currents. Now assuming that there are N lowest-lying intermediate states (indexed as m = 0, · · · , N − 1) with energy E * m ≤ E f +E 1 , and N lowest-lying intermediate states (indexed as m = 0, · · · , N −1) with energy E * m ≤ E i − E 1 , the contribution from all states except the lowest N or N states in the theory can be directly analytically continued to Minkowski space. With this observation, one can break the Minkowski ME T (M) L to two parts as following: where the first term, defined as cannot be directly accessed via the Euclidean four-point correlation function, and can only be reconstructed from the knowledge of N or N (whichever larger) finite-volume two-nucleon eigenenergies and the MEs of the weak current between initial/final two-nucleon states and the N or N allowed intermediate states. These must be evaluated using LQCD computations of two-and three-point correlation functions in a finite volume. The second term in Eq. (88) needs the values of the fourpoint function evaluated with LQCD, as well as the two-and three-point functions that need to be evaluated in obtaining the first term. Explicitly, where G L is the Euclidean bi-local MEs defined in Eq. (87) and G < L is: with the c m function defined in terms of finite-volume single-β decay MEs satisfying either E * m ≤ E f + E 1 or E * m ≤ E i − E 1 , corresponding to the first and second sum in Eq. (91), respectively: Given all these ingredients, the physical (infinite Minkowski spacetime) double-β decay amplitude can be extracted from a (finite) set of relevant two-, three-, and four-point correlation functions in a Euclidean finite spacetime, quantities that are all accessible from LQCD. In fact, as pointed out in Ref. [88], it will be essential that the matching relation, including Minkowski spacetime construction and the finite-volume corrections to the ME, are all fit to LQCD data on the required n-point functions simultaneously to insure the exact cancelations among various contributions that lead to the correct analytic structure of the infinite-volume amplitude.
To close, it should be noted that a closely-related procedure was implemented in Refs. [19,20] to isolate the deuteron-pole contribution from the Euclidean four-point correlation function, resulting in the desired double-β decay ME. Nonetheless, the initial and final states in that study were bound nuclear states given larger-than-physical values of the quark masses used in the study, with only one bound on-shell intermediate state (the deuteron). This, therefore, required a simpler formalism to extract the physical ME with no need for a FV matching. The complete framework of this work, therefore, will be necessary in future LQCD calculations of the same process at lighter values of the quark masses, where two-nucleon states will no longer be bound. The limitations of the applicability of the formalism presented are its reliance on a finite-order EFT expansion, which can be straightforwardly improved by including higher-order contributions in the analysis, and its validity only below three-hadron thresholds. This latter limitation will be more challenging to relax, but progress is plausible given the finite-volume technology developed in recent years for the three-hadron problem from LQCD, see Ref. [69] for a review.

V. CONCLUSION AND OUTLOOK
Lattice quantum chromodynamics coupled with applicable effective field theories is becoming a critical component of modern and systematic approaches to understanding and predicting certain nuclear phenomena from first principles by supplementing quantum many-body methods in nuclear physics [102][103][104][105][106][107][108][109]. This paper provides the elementary building blocks for implementing such a perspective in the case of single-and double-weak processes in nuclei, with a focus on the lessdeveloped process of two-neutrino double-β decay. This work consists of two components: Within the framework of pionless EFT with nucleonic degrees of freedom, which is a good description of the nn → pp (eeν eνe ) process at low energies, the physical amplitude is constructed for the first time in this work. While this process is not expected to occur in nature given the short lifetime of the neutron compared with the inverse rate of the double-weak process, such a transition is still the primary process occurring in isotopes undergoing a double-β decay of ordinary type. As a result, any constraint on the relevant nuclear matrix element at the two-nucleon level serves as the microscopic input to quantum many-body calculations of the rate. Investigating various contributions to the amplitude, including those with electron-neutrino emission on the external nucleons, are shown to confirm the working principles of the pionless EFT as a renormalized approach to nuclear observables when applied to a double-weak amplitude.
No new low-energy constants beyond the nucleon (g A ) and the two-nucleon (L 1,A ) axialvector couplings are shown to be needed to guarantee a renormalization-scale independent amplitude at next-to-leading order. This is in contrary to the case of pionless EFT with dibaryon DOF in Refs. [19,20], where a new short-distance two-nucleon two-weak current coupling, h 2,S , was introduced, and its effect on the amplitude at unphysically large values of the quark masses (determined with LQCD) was shown to be comparable to the contributions from the l 1,A coupling (of the dibaryon formalism). By examining the consistency of L 1,A values when constrained from single-and double-weak processes independently, future analysis of LQCD results will determine the validity of the pionless EFT power counting with nucleonic DOF as derived in this paper, and may have implications for 'quenching' of g A in double-β decays.
Despite the lack of experimental constraints on the process, LQCD can still give access to the nonlocal correlation functions of two-nucleon systems with two insertions of the weak current, hence determining the relevant MEs contributing to the physical amplitude. There is, however, a nontrivial procedure involved to perform this matching, as developed for the first time in this paper for this problem. The procedure builds upon the previously developed finite-volume formalisms for the local MEs of two nucleons and nonlocal MEs of a single hadron. It requires the determination of the finite-volume two-nucleon spectrum in the spin-singlet and spin-triplet channels, and the three-point and four-point correlation functions corresponding to the single-β and double-β decay processes from LQCD for a range of kinematics. The reason is that the Euclidean and Minkowski correlation functions can not be mapped to each other straightforwardly when on-shell multi-hadron intermediate states are present, whose effects must be isolated and treated separately [88]. Furthermore, by assuming that only two-nucleon intermediate states can go on shell, power-law volume corrections to the amplitude are derived within the pionless-EFT approach at NLO, completing the matching relation.
This formalism brings us one step closer to the ultimate goal of constraining the physical amplitude for the neutrinoless process nn → pp (ee) with the light Majorana exchange, or in turn the new short-distance LEC in pionless EFT that is present at LO [29]. Obtaining this constraint from a Euclidean four-point correlation function involves the convolution of the nuclear ME with the neutrino propagator, further complicating the relation to the Minkowski correlation function. Nonetheless, upon slight modifications, identification of the finite-volume corrections to the reconstructed Minkowski amplitude follows the same procedure as outlined in this work. These steps will be presented in future work [93].
Finally, an interesting observation is regarding the utility of Eq. (85) even when the amplitude is computed using a quantum simulator or a quantum computer. In a quantum simulation, the most natural approach is to construct the states and evolve and measure them using analog or digital simulation protocols. As a result, no Monte Carlo importance sampling of quantum field configurations is required to demand a Wick rotation to Euclidean spacetime, as is the case in LQCD computations. Nonetheless, a quantum processor does not provide infinite computing resources, requiring truncations on the physical Hilbert space of the fields. Among the unavoidable truncations is that of the spatial extent of the system, bringing in the need to determine the infinite-volume limit of observables from finite-volume quantities. When it comes to multi-hadron observables, approaching such a limit can be involved. One approach is to take the ordered double limit L → ∞ and → 0 in the finite-volume Minkowski amplitude, which as shown in Ref. [110] in given examples, may require large spatial volumes to allow convergence to the infinite-volume value. Another approach employs the finite-volume formalisms such as that developed in this work, e.g., Eq. (85), to provide the mapping between the Minkowski but finite-volume ME and the physical (Minkowski and infinite-volume) amplitude at a given volume. The latter circumvents the need for the involved procedure described in Sec. IV D, but is limited to low-energy kinematics below three-hadron inelastic thresholds, which, however, can be generalized. This is an example of the applicability of a range of theoretical developments in LQCD for hadronic and nuclear physics to the quantum simulation of these strongly-interacting systems in nature.