Few-nucleon matrix elements in pionless effective field theory in a finite volume
MMIT-CTP/5275
Few-nucleon matrix elements in pionless effective field theory in a finite volume
W. Detmold and P. E. Shanahan
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA andThe NSF AI Institute for Artificial Intelligence and Fundamental Interactions (Dated: February 9, 2021)Pionless effective field theory in a finite volume (FVEFT π/ ) is investigated as a framework for theanalysis of multi-nucleon spectra and matrix elements calculated in lattice QCD (LQCD). By com-bining FVEFT π/ with the stochastic variational method, the spectra of nuclei with atomic number A ∈ { , } are matched to existing finite-volume LQCD calculations at heavier-than-physical quarkmasses corresponding to a pion mass m π = 806 MeV, thereby enabling infinite-volume bindingenergies to be determined using infinite-volume variational calculations. Based on the variationalwavefunctions that are constructed in this approach, the finite-volume matrix elements of variouslocal operators are computed in FVEFT π/ and matched to LQCD calculations of the correspond-ing QCD operators in the same volume, thereby determining the relevant one and two-body EFTcounterterms and enabling an extrapolation of the LQCD matrix elements to infinite volume. Asexamples, the scalar, tensor, and axial matrix elements are considered, as well as the magneticmoments and the isovector longitudinal momentum fraction. I. INTRODUCTION
Over the last 30 years, effective field theories have rev-olutionised nuclear physics, systematising the study ofnucleon-nucleon interactions and the properties of lightnuclei. Pionless effective field theory for few-nucleon sys-tems (EFT π/ ) in particular, which focuses on momentabelow the pion mass, has emerged as a powerful toolwith which to understand low-energy nuclear processesin many contexts [1–5] (see Ref. [6] for a recent review).Notably, in addition to its use in the analysis of experi-mental data, EFT π/ has found a key role in the analysisof lattice Quantum Chromodynamics (LQCD) calcula-tions of nuclear systems, providing a direct bridge be-tween QCD and nuclear physics.For example, in Ref. [7], LQCD calculations of A ∈{ , } nuclei at heavier-than-physical quark masses [8]were matched to auxiliary field diffusion Monte-Carlo cal-culations with EFT π/ interactions to constrain the twoand three-body counterterms of the EFT; these werethen used to make predictions for larger nuclei of atomicnumber A ≤
6. Further developments were presented inRef. [9], and this approach was extended to the next or-der in the EFT, and to still larger nuclei, in Refs. [10]and [11]. Ref. [12] presents studies of the quark-mass de-pendence of the magnetic moments and polarisabilities of A ∈ { , } nuclear systems, with both experimental re-sults and LQCD calculations at larger-than-physical val-ues of the pion mass used to constrain the countertermsof EFT π/ . Similarly, an early application of finite-volumeEFT π/ to electroweak matrix elements was presented inRef. [13] and extended in Ref. [14]. Furthermore, EFT π/ provided a powerful approach to analysing LQCD calcu-lations of second order weak processes [15–19].Recently, Eliyahu et al. [20] have taken the next stepsin this approach and used EFT π/ implemented via thestochastic variational method (SVM) [21, 22] in a finitecubic volume to analyse the binding energies of atomicnumber A ∈ { , } systems in that finite volume. Since the effects of a finite volume manifest in the infrared do-main, they can be captured in low-energy effective de-scriptions of QCD such as that provided by EFT π/ . SinceLQCD calculations are typically performed in multiplevolumes to enable an infinite-volume extrapolation, per-forming the EFT in the same volumes also maximisesthe constraining power of the LQCD results. As input,Ref. [20] utilised the NPLQCD collaboration’s LQCDcomputations of the ground-state energies and finite-volume energy shifts of these systems in 3 different lat-tice volumes with spatial extents L ∈ { . , . , . } fmat unphysical quark masses corresponding to the SU(3) f flavour-symmetric point where the up, down, and strangequark masses are degenerate and correspond to a pionmass m π = m K = 806 MeV [23]. Using EFT π/ in thesame volumes to determine counterterms, the bindingenergies were extrapolated to infinite volume. As thisexemplifies, the finite-volume pionless effective field the-ory (FVEFT π/ ) approach provides a powerful alternativeto L¨uscher’s method [24], in which finite-volume energiesfrom LQCD calculations are used to determine infinite-volume scattering phase shifts and bound-state energies.While L¨uscher’s method and its generalisations are modelindependent, the existing formalism is limited to two andthree-particle systems. The matching of LQCD results toEFT π/ , on the other hand, requires an underlying EFTbut can be applied to any system that the EFT can ad-dress.In this work, the application of FVEFT π/ is extended tonuclear matrix elements for the first time. After a briefsummary of the relevant EFT π/ Lagrangians in Sec. II,the SVM is introduced in Sec. III. Section IV presentsthe results of tuning the relevant two and three-bodycounterterms of FVEFT to reproduce the ground-stateenergies of A ∈ { , } nuclei computed in LQCD, paral-leling the analysis of Ref. [20]. Having determined thesecounterterms, Sec. V presents the tuning of the countert-erms describing the interactions of external currents toreproduce LQCD calculations of nuclear matrix elements, a r X i v : . [ nu c l - t h ] F e b thereby enabling an extrapolation of the finite-volumematrix elements to infinite volume. In particular, LQCDcalculations of scalar, tensor, and axial matrix elements,as well as magnetic moments and isovector longitudinalmomentum fractions, of A ∈ { , } states are investi-gated. To conclude, the outlook for these calculations isdiscussed in Sec. VI. II. PIONLESS EFFECTIVE FIELD THEORY
The pionless EFT Lagrangian describing the low-energy interactions of nucleons is given by L = L + L + L + . . . , (1)where L = N † (cid:18) iD + D M N (cid:19) N + . . . (2)contains the single-nucleon kinetic operator expanded inthe non-relativistic (NR) limit. Here, N represents thenucleon field, M N is the nucleon mass, and the ellipsis de-notes higher-order terms. The leading-order two-nucleoninteractions enter as L = − C S (cid:0) N T P i N (cid:1) † (cid:0) N T P i N (cid:1) − C T (cid:0) N T P a N (cid:1) † (cid:0) N T P a N (cid:1) + . . . (3)where P i ≡ √ σ σ i τ , P a ≡ √ σ τ τ a (4)are projectors onto spin-triplet and sin-singlet two-nucleon states respectively, and C S and C T are the rel-evant two-body low-energy constants (LECs). Here σ k and τ a are the Pauli matrices acting in spin and isospinspace, respectively. Eq. (3) can be re-expressed in a dif-ferent basis as [25] L = − (cid:104) C (cid:0) N † N (cid:1) + C (cid:0) N † (cid:126)σN (cid:1) (cid:105) + . . . , (5)where C T = C − C and C S = C + C . (6)Three-body interactions naively enter at higher order,but must be promoted to leading order as argued inRefs. [26, 27], and the relevant contribution to the La-grangian is L = − D N † N ) + . . . , (7)where D is the leading-order three-body LEC. A. Weak interactions
The weak decays and interactions of nuclear statesarise, after integrating out the weak gauge boson,through the effective Lagrangian (valid for energies E (cid:28) M W ) L W = − G F √ l µ + J − µ + h.c. + · · · , (8)where G F is the Fermi constant, l µ + involves a chargedlepton and neutrino, and the hadronic weak current canbe expressed in terms of the vector ( V µ ) and axial-vector( A µ ) currents as J i,µ = V i,µ − A i,µ , (9)with J ± µ = J ,µ ± i J ,µ . The isovector axial-vector cur-rent in EFT π/ is given by [28] A i,a = g A N † τ a σ i N + L ,A (cid:0) N T P i N (cid:1) † (cid:0) N T P a N (cid:1) + h.c. + . . . , (10)where the ellipsis denotes higher-order terms. In this ex-pression, g A is the nucleon axial charge, and the term pro-portional to the two-nucleon LEC L ,A provides the next-to-leading-order (NLO) corrections to the pp → de + ν fusion process, for example.The isoscalar axial current is similarly given by [29] A i, = − g A, N † σ i N − iL ,A (cid:15) ijk (cid:0) N T P j N (cid:1) † (cid:0) N T P k N (cid:1) + · · · , (11)where the isoscalar nucleon axial charge is g A, , and theLEC L ,A enters at the same order. B. Electromagnetic interactions
In the presence of an external electromagnetic (EM)field, the Lagrangians defined above are modified suchthat the derivatives are replaced by EM-covariant deriva-tives D µ = ∂ µ + iQA µ , where A µ is the vector potentialand Q is the electric charge operator, and also by theaddition of terms depending on the magnetic field B : L , EM + L , EM = J EMi B i , (12) Note that the normalisation of L ,A used here is the same as inRef. [28] but differs from Ref. [29], which uses a different projec-tor definition. Note that Ref. [29] uses a different projector definition than here,and correspondingly a different normalisation of L ,A . The notation of Ref. [30] is used. where the isoscalar and isovector currents coupling to themagnetic field are J EMi = e M N N † ( κ + τ κ ) σ i N − eL i(cid:15) ijk (cid:0) N T P k N (cid:1) † (cid:0) N T P j N (cid:1) + eL (cid:0) N T P i N (cid:1) † (cid:0) N T ¯ P N (cid:1) + h.c. , (13)where L and L are two-body LECs and κ = 12 ( κ p + κ n ) and κ = 12 ( κ p − κ n ) (14)are the isoscalar and isovector nucleon magnetic mo-ments. Note that electric field contributions and EMthree-body interactions enter at higher order. C. Scalar and tensor currents
The isovector and isoscalar scalar currents that arisefrom Higgs couplings and from potential dark matter in-teractions are given by S = g S, N † N − (cid:101) C S (cid:0) N T P i N (cid:1) † (cid:0) N T P i N (cid:1) , − (cid:101) C T (cid:0) N T P a N (cid:1) † (cid:0) N T P a N (cid:1) + . . . , (15) S a = g S, N † τ a N + i (cid:101) C V (cid:15) abc (cid:0) N T P b N (cid:1) † (cid:0) N T P c N (cid:1) + . . . . (16)Here g S, and g S, are the isoscalar and isovector one-body LECs that are related to the nucleon σ terms. Asdiscussed in Ref. [31], the two-body terms in the isoscalarscalar current are related to the corresponding terms inthe strong Lagrangian, Eq. (3). In particular, the LECs (cid:101) C S,T are the quark-mass–independent pieces of the La-grangian couplings C S,T .For the isoscalar and isovector antisymmetric tensorcurrents, the relevant EFT π/ expressions are T ij, = g T, (cid:15) ijk N † σ k N + iL ,T (cid:0) N T P i N (cid:1) † (cid:0) N T P j N (cid:1) + h.c. + . . . , (17) T ij,a = g T, (cid:15) ijk N † τ a σ k N + L ,T (cid:15) ijk (cid:0) N T P k N (cid:1) † (cid:0) N T P a N (cid:1) + h.c. + . . . , (18)where the one and two-body isoscalar (isovector) tensorLECs are g T, and L ,T ( g T, and L ,T ). D. Twist-two operators
The unpolarised twist-two operators that define mo-ments of parton distributions enter in EFT π/ as [32, 33] O µ ...µ n = (cid:104) x n (cid:105) v µ · · · v µ n N † N (cid:2) α n, N † N (cid:3) , (19) O µ ...µ n = (cid:104) x n (cid:105) v µ · · · v µ n N † τ N (cid:2) α n, N † N (cid:3) , (20)with subleading contributions from terms involvingderivatives that are suppressed in the power-counting,as well as from additional two-body terms that are notWigner SU(4) symmetric and are suppressed in the large N c limit. The subscripts 0 and 3 denotes the isoscalarand isovector combinations. Note that the isoscalar con-tributions arise from matching to both quark and gluonmatrix elements (which mix under QCD renormalisa-tion), and will give rise to finite-volume effects that arethe same in both cases. III. STOCHASTIC VARIATIONAL METHOD INA PERIODIC CUBIC VOLUME
In order to address finite-volume effects in few-nucleonsystems in EFT, few-body wavefunctions must be de-termined subject to the EFT interactions and the givenboundary conditions. There are multiple many-body ap-proaches that could be pursued for this task. Two ap-proaches that have been successfully applied are solvingthe 3-dimensional finite-volume Schr¨odinger equation viadiscretisation [8], and the stochastic variational method(SVM) for two and three-body systems [20]. The formerapproach works well for two-body interactions and waseffectively used in Ref. [8] to analyse hyperon-nucleon in-teractions where the effective range was not small com-pared to the spatial extent of the finite volume, L , and assuch the more direct L¨uscher method [24] could not beapplied. However, this coordinate-space–based approachscales poorly to larger systems.The SVM was introduced in nuclear physics [21] asa way to sample the possible spatial, spin, and isospin-wavefunctions for an A -nucleon system in a space thatis impractically large for an exhaustive approach, seeRefs. [22, 34] for reviews. This approach, detailed be-low (and applied to EFT π/ in Ref. [35], for example),involves the construction of a wavefunction by sequen-tial proposals of new stochastically-generated terms andthe optimisation of the linear coefficients of the terms bysolving the generalised eigenvalue problem of the varia-tional method. The SVM for a finite volume was firstintroduced in Ref. [36] where systems of bosons in pe-riodic cubic potentials were considered. Periodicity isimposed on the wavefunctions by considering all periodiccopies of the infinite-volume potential. The method wasfirst used for nuclei in Ref. [20], and a similar approachis used here. A. Finite volume Hamiltonian
The n -particle non-relativistic Hamiltonian that corre-sponds to the EFT π/ Lagrangian of Eq. (1) is H = − M N (cid:88) i ∇ i + (cid:88) i While in many applications of the SVM to nuclear sys-tems the angular momentum structure of the wavefunc-tion is tied to the spatial structure due to orbital mo-tion, in a cubic box orbital angular momentum is nota well-defined quantum number. As will be discussedbelow, one approach in this context is to build wavefunc-tions with particular transformation properties under thecubic group, again coupling spatial and spin degrees offreedom. However, since there are a finite number of ir-reducible representations of the cubic group, a simpler approach is to consider a factorisation of the spatial andspin-isospin wavefunctions.In this work, a trial wave function, Ψ ( N ) h , is built fromlinear combinations of symmetrised spatial wavefunctionsΨ sym L which satisfy periodic boundary conditions, cou-pled to the appropriate spin-flavour wavefunction | χ h (cid:105) ,i.e., Ψ ( N ) h ( x ) = N (cid:88) j =1 c j Ψ sym L ( A j , B j , d j ; x ) | χ h (cid:105) , (26)where the superscript ( N ) denotes the total number ofterms in the wavefunction, and the dependence of thespatial wavefunction on h is suppressed. The coordinate x = ( r , . . . , r n ) collects the spatial coordinates of the n particles with x j = r j . The spatial wavefunctions Ψ sym L are detailed in Sec. III B 1; the c j , j ∈ { , . . . , N } , arecoefficients, and the A j , B j and d j are the parameters ofthe j th spatial wavefunction included in the sum. Thespin-flavour wavefunction | χ h (cid:105) is a vector in spin-flavourspace for the given nucleus h ; the particular spin-flavourwavefunctions that are used in this work are given inSec. III B 3. 1. Shifted correlated Gaussian spatial wavefunctions To account for the anisotropy in the spatial wavefunc-tion due to the boundary conditions in a cubic volume,a shifted correlated-Gaussian basis for the trial wave-functions is used following the approach introduced inRef. [36]. States are constructed to be antisymmetricunder interchange of the spin-flavour degrees of freedomof pairs of nucleons and thus must have symmetric spatialwavefunctions under particle interchange.The basic Gaussian structure underlying these wavefunctions isΨ ( α ) ∞ ( A ( α ) , B ( α ) , d ( α ) ; x ( α ) ) = exp (cid:20) − x ( α ) T A ( α ) x ( α ) − 12 ( x ( α ) − d ( α ) ) T B ( α ) ( x ( α ) − d ( α ) ) (cid:21) , (27)where x ( α ) is an n -component vector collecting the α thCartesian component of the position of each particle. The n × n matrices A ( α ) and B ( α ) , and n -component vector d ( α ) , contain the parameters defining the wavefunction.The matrices A ( α ) are symmetric, containing n ( n − / B ( α ) are diagonal matrices with n real parameters. The finite-volume approach introducedby Yin and Blume [36] is implemented through sums of Note that χ is common to all terms in Eq. (26) in the currentimplementation. In other approaches for larger systems than willbe considered here, χ is also part of the stochastic sampling andwould be indexed by j [34]. periodic copies of the intrinsic wavefunction over shiftedvolumes to define a finite-volume wavefunctionΨ L ( A, B, d ; x ) = (cid:89) α ∈{ x,y,z } Ψ ( α ) L (cid:16) A ( α ) , B ( α ) , d ( α ) ; x ( α ) (cid:17) , (28)where A = diag { A ( x ) , A ( y ) , A ( z ) } is a block-diagonal3 n × n matrix that combines the A ( α ) matrices for eachdirection, and similarly B and d combine the B ( α ) and d ( α ) for each direction. The wavefunction for the α thdirection isΨ ( α ) L (cid:16) A ( α ) , B ( α ) , d ( α ) ; x ( α ) (cid:17) = (cid:88) b ( α ) Ψ ( α ) ∞ ( A ( α ) , B ( α ) , d ( α ) ; x ( α ) − b ( α ) L ) , (29)where the n -component vector b ( α ) has components b ( α ) j ∈ Z . The resulting wavefunction Ψ L ( A, B, d ; x ) sat-isfies the periodic constraintΨ L ( A, B, d ; x ) = Ψ L ( A, B, d ; x + n L ) (30)for all integer 3 n -tuples , n ∈ Z n . In order to sym-metrise the wavefunction under particle exchange, therows and columns of each A ( α ) and B ( α ) and rows of d ( α ) are interchanged under all n ! possible permutations, P , of particles. That isΨ sym L ( A, B, d ; x ) = (cid:88) P Ψ L ( A P , B P , d P ; x ) , (31)where A P is the permuted form of A and similarly for B P and d P . As discussed in Appendix B, the shiftedGaussian basis is able to describe scattering states atfinite volume as well as compact bound states. 2. Cubic harmonics A periodic spatial volume with identical extent in eachdirection has an underlying cubic symmetry and is in-variant under action of elements of the cubic group, H .By imposing cubic symmetry, wavefunctions that trans-form in particular representations of H can be con-structed, potentially allowing more efficient explorationof the space of correlated shifted Gaussians. Each termin the variational wavefunction can be constructed to re-spect the given transformation properties rather than re-lying on stochastic sampling of a sum of terms to discoverthe symmetry approximately. The H -covariant wave-function transforming in the representation R of H isgiven byΨ sym ,RL ( A, B, d ; x ) = (cid:88) p c ( R ) p Ψ sym L ( A p , B p , d p ; x p ) , (32) In practice, these sums are truncated as discussed in AppendixC. where p indexes the permutations of the Cartesian direc-tions, c ( R ) p are constants defining the representation, and A p is the appropriately block-permuted form of the ma-trix A and similarly for B p and d p . For the ground statesthat are considered here, the A (trivial) representationof H is assumed, for which c ( A ) p = 1. The utility of usingEq. (32) instead of Eq. (31), has been investigated. Over-all, it is found that N -term wavefunctions constructedfrom terms of the form Ψ sym ,RL are about a factor of fivebetter approximations than N -term wavefunctions con-structed from terms from Eq. (31), measured in termsof the number of wavefunction terms required to achieveconvergence within a given tolerance. However the costof evaluation of the matrix elements needed in the SVMis a factor of six slower using Eq. (32) than using Eq. (31).In the primary studies of this work, trial wavefunctionsare thus constructed using the simpler ansatz in Eq. (31). 3. Spin-flavour wavefunctions The simplest spin-flavour wavefunctions for the smallnuclei that are considered in this work are straightfor-ward to construct explicitly. In particular, the necessarystates are defined as | χ d,j z =+1 (cid:105) = 1 √ (cid:2)(cid:12)(cid:12) p ↑ n ↑ (cid:11) − (cid:12)(cid:12) n ↑ p ↑ (cid:11)(cid:3) , | χ d,j z =0 (cid:105) = 12 (cid:2)(cid:12)(cid:12) p ↑ n ↓ (cid:11) − (cid:12)(cid:12) n ↑ p ↓ (cid:11) + (cid:12)(cid:12) p ↓ n ↑ (cid:11) − (cid:12)(cid:12) n ↓ p ↑ (cid:11)(cid:3) , | χ pp (cid:105) = 1 √ (cid:2)(cid:12)(cid:12) p ↑ p ↓ (cid:11) − (cid:12)(cid:12) p ↓ p ↑ (cid:11)(cid:3) , | χ np,j =0 (cid:105) = 12 (cid:2)(cid:12)(cid:12) p ↑ n ↓ (cid:11) + (cid:12)(cid:12) n ↑ p ↓ (cid:11) − (cid:12)(cid:12) p ↓ n ↑ (cid:11) − (cid:12)(cid:12) n ↓ p ↑ (cid:11)(cid:3) , (cid:12)(cid:12) χ H ,j z =1 / (cid:11) = 1 √ (cid:2)(cid:12)(cid:12) n ↑ p ↑ n ↓ (cid:11) − (cid:12)(cid:12) n ↓ p ↑ n ↑ (cid:11) − (cid:12)(cid:12) p ↑ n ↑ n ↓ (cid:11) + (cid:12)(cid:12) p ↑ n ↓ n ↑ (cid:11) − (cid:12)(cid:12) n ↑ n ↓ p ↑ (cid:11) + (cid:12)(cid:12) n ↓ n ↑ p ↑ (cid:11)(cid:3) , (cid:12)(cid:12) χ He ,j z =1 / (cid:11) = 1 √ (cid:2)(cid:12)(cid:12) p ↑ n ↑ p ↓ (cid:11) − (cid:12)(cid:12) p ↓ n ↑ p ↑ (cid:11) − (cid:12)(cid:12) n ↑ p ↑ p ↓ (cid:11) + (cid:12)(cid:12) n ↑ p ↓ p ↑ (cid:11) − (cid:12)(cid:12) p ↑ p ↓ n ↑ (cid:11) + (cid:12)(cid:12) p ↓ p ↑ n ↑ (cid:11)(cid:3) , (33)where p ↑ ( ↓ ) and n ↑ ( ↓ ) denote proton and neutron statesof the given spin. C. Implementation of the stochastic variationalmethod The trial wavefunction Ψ ( N ) of Eq. (26) is constructedso as to minimize the bound which it provides on theground-state energy: E h ≤ (cid:82) Ψ ( N ) h ( x ) ∗ H Ψ ( N ) h ( x ) d x (cid:82) Ψ ( N ) h ( x ) ∗ Ψ ( N ) h ( x ) d x . (34)Because of the Gaussian structure of the trial wavefunc-tion, the various contributions to the Hamiltonian matrixelement, which are 3 n -dimensional integrals, can be eval-uated analytically [36], as shown in Appendix C.In the current application of the SVM, Ψ ( N ) is built upfrom 1 to N terms via an iterative procedure as follows:1. Given an M -term wavefunction (where M < N )defined by matrices A j , B j , d j for j ∈ { , . . . , M } , N proposal proposed candidates for the ( M + 1)thterm are constructed by randomly choosing matri-ces A M +1 , B M +1 and d M +1 . For simplicity of no-tation, the spatial wavefunction of the j th term isdenoted as Ψ j ( x ) ≡ Ψ sym L ( A j , B j , d j ; x ).2. For each candidate term Ψ M +1 ( x ), the normalisa-tion integrals[ N ( M +1) ] ij ≡ (cid:90) Ψ i ( x ) ∗ Ψ j ( x ) d x (35)and Hamiltonian matrix elements[ H ( M +1) ] ij ≡ (cid:90) Ψ i ( x ) ∗ (cid:104) χ h | H | χ h (cid:105) Ψ j ( x ) d x (36)are computed for { i, j } ∈ { , . . . , M + 1 } and usedto define matrices N ( M +1) and H ( M +1) , respec-tively. Note that only the additional ( M + 1)throw and column must be computed, given that the N ( M ) and H ( M ) matrices were computed in the pre-vious iteration.3. The generalised eigenvalue problem H ( M +1) c = λ N ( M +1) c (37)is solved for the eigenvalues λ ( M +1)0 ≤ λ ( M +1)1 ≤ . . . ≤ λ ( M +1) M +1 , and coefficient vectors c ( M +1) (cid:96) =( c , . . . , c M +1 ) for (cid:96) ∈ { , . . . , M + 1 } labelling theeigenvalue.4. The wavefunction which results in the small-est eigenvalue λ ( M +1)0 is selected from theset of N proposal candidates and added to theiteratively-constructed trial wavefunction to defineΨ ( M +1) ( x ).Naturally, the optimization at each step of this itera-tive procedure depends on the Hamiltonian H and henceon the LECs that define it; to enable optimization ofthe trial wavefunction across a broad range of LECs, inthe numerical study undertaken here the values of theLECs are varied for each step of optimization, cycling Note that each of the N proposal lowest eigenvalues λ ( M +1)0 issmaller than the lowest eigenvalue from the previous iterationˆ λ ( M )0 . h L = 3 . L = 4 . L = 6 . pp . . 3) 15 . . 8) 15 . . d . . 4) 22 . . 5) 19 . . H 65 . . 8) 63 . . 0) 53 . . through N couplings choices that span the relevant param-eter spaces. After initializing the procedure with the firsttrial wavefunction (the M = 1 term), for which the ma-trices N (1) and H (1) are single numbers and the gener-alised eigenvalue problem is trivial, additional terms areadded iteratively until the wavefunction has a fixed num-ber of terms, N . N must be taken large enough thatthe optimization procedure has converged by some defi-nition. Here, N is set by the criterion that repeated op-timizations based on different random seeds achieve thesame minimum energies within some tolerance, and thatadding some fixed number of additional terms to the trialwavefunctions does not alter the minimum found, withinthe same tolerance. Details of the optimisation proce-dure for the cases considered in this study are providedin Sec. IV. IV. GROUND STATES OF TWO ANDTHREE-NUCLEON SYSTEMS In this work, the finite-volume SVM is used tomatch the LECs of FVEFT π/ to the LQCD resultsfor two and three-nucleon systems which were ob-tained in Ref. [8], where nuclear states with SU(3) f -symmetric quark masses corresponding to m π = 806MeV were studied in three volumes of spatial extents L ∈ { . , . , . } fm. The binding energies extractedin that work, defined as ∆ E h = E h − AE p where A isthe atomic number of the state h , are tabulated in TableI. This matching procedure, used in this work as a firststep in the study of the matrix elements of various cur-rents in these states in this framework, closely mirrors,and reproduces, the analysis of Ref. [20]. A. Two-body states For each of the three finite volumes for which Ref. [8]provides LQCD data, variational wavefunctions were op-timised first for the pp and d two-body systems. Asdescribed in Sec. III C, the undetermined coefficient ofthe two-body potential in the Hamiltonian (which corre-sponds to the LEC C S in the case of the deuteron and C T in the case of pp ) is varied over N couplings = 10 differentchoices throughout the optimisation procedure, which arechosen to be evenly-spaced corresponding approximatelyto the plot range of Fig. 2; both two-body systems arethus optimised simultaneously. Taking N proposal = 30, □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○ - - - - - - □ △ ○ FIG. 1. Convergence of the eigenvalues λ ( N )0 to the ground-state energy of the diproton system as additional terms areadded to the variational wavefunction. The three colours cor-respond to wavefunctions optimized for the L = 3 . C T set to its optimised value after matching to the LQCDresults, as discussed in the main text. The results are shownfor one example of the regulator parameter corresponding to r = 0 . it is found for all two-body optimizations undertaken inthis work—at each finite volume, and for each of threechoices of the regulator Λ = √ /r corresponding to r ∈ { . , . , . } fm—that after 100 terms have beenadded to the wavefunction, the last 20 terms are within1% (typically, within 0.1%) of the final value for each ofthe values of the coupling that are used in optimization,and also that optimizations starting with different ran-dom seeds agree within that same tolerance. An exampleof this convergence is shown in Fig. 1.The ground-state d and pp binding energies resultingfrom the optimised variational wavefunction are shown asa function of the relevant undetermined LEC in Fig. 2,for all three finite volumes and for all choices of the reg-ulator parameter Λ which are studied. Comparing withthe LQCD results for the binding energies in each lat-tice volume, it is clear that for each value of Λ the re-sults in all volumes are consistent with a single consis-tent value of the relevant coupling, indicating that thereis no need to introduce higher order terms in EFT π/ . Thebest-fit values of the corresponding couplings, which de-pend on the regulator scale, are determined through acombined fit to the three volumes and are presented inTable II and Fig. 3. Note that the EFT π/ interactionproportional to C is suppressed by 1 /N c relative to theWigner-symmetric interaction with coefficient C in thelarge- N c limit; this hierarchy is born out in the fittedvalues of the couplings.Having determined the two-body couplings, wavefunc-tions optimised in infinite volume in exactly the same wayare used to determine the binding energies in the infinite-volume limit. Figure 4 shows the volume-dependence of r [fm] 0 . . . C − − − C − − − C S − − − C T − − − D π/ LECs determinedfrom matching the SVM calculations to the LQCD energyshifts for each value of the cutoff parameter r . For the two-body case, C S and C T are determined in terms of C , throughEq. (6). For D , three-nucleon optimisations were only per-formed for r = 0 . - - - - 100 0 - - - - - - - - - 100 0 - - - - - FIG. 2. Binding energies of the deuteron (upper panel) anddiproton (lower panel) systems as a function of the relevanttwo-body EFT π/ LECs. The three sets of curves correspond tothe three different choices of the regulator scale Λ = √ /r corresponding to r ∈ { . , . , . } fm (solid, dashed, anddotted) and the three colours correspond to the three differ-ent volumes: L = 3 . the binding energies of the deuteron and diproton sys-tems (in order to show a curve, binding energies com-puted with wavefunctions optimised at additional inter-mediate volumes are also shown). Extrapolations areshown for r = 0 . * * * * - - - - - - - - - - FIG. 3. The EFT π/ couplings C , determined by fitting tothe results of the LQCD calculations. The three sets of in-tersecting bands (blue for pp and red for d ) and correspond-ing ellipses show results obtained with wavefunctions opti-mised with the three different values of the regulator scaleΛ = √ /r studied here. The asterisks show the results ofan analogous analysis undertaken in Ref. [20], with differentvalues of the regulator r , as indicated on the figure. □ □ □□ □ □ □□□□◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ - - - - - □ □ FIG. 4. The volume-dependence of the deuteron and dipro-ton binding energies compared with the LQCD data whichwas used to determine the relevant LECs. The infinite-volumeextrapolations of the binding energies, computed as describedin the text, are shown in the rightmost sub-panel. values are indistinguishable. Although the values of theLECs depend on the value of the regulator r , the re-sulting finite and infinite-volume energies are regulator-independent. Table III lists the extrapolated binding en-ergies and compares them to Refs. [20, 23], with whichthey are in close agreement. B. Three-body states Having determined the two-nucleon couplings, an anal-ogous procedure is repeated for the triton to determinethe coefficient D of the leading-order three-nucleon cou-pling in Eq. (7). The stochastic optimisation of the three-body wavefunction is performed with the two-body LECs h L = ∞ Ref. [23] Ref. [20]pp -12.5(2.2) -15.9(3.8) -13.8(1.8) H -19.9(2.8) -19.5(4.8) -20.2(2.3) H .-60.2(6.5) -53.9(10.7) -58.2(4.7)TABLE III. The extrapolated infinite-volume binding ener-gies determined in the SVM approach for the two and three-body systems. Also shown are the binding energies deter-mined in the original LQCD study [23] and in Ref. [20] alsousing the SVM method. - - - - - - FIG. 5. The dependence of the triton binding energy onthe three-body coupling D in each of three finite volumesfor which there is LQCD data, L = 3 . r = 0 . fixed to their central values determined as discussed inthe previous subsection, for a single value of the Gaussianregulator parameter r = 0 . N couplings = 10 values of the three-body LEC D arecycled through in the wavefunction construction proce-dure, spanning the relevant coupling range (correspond-ing approximately to the range of the horizontal axis inFig. 5). The same convergence criteria as in the two-body case are satisfied after wavefunctions with N = 250terms have been constructed.Figure 5 shows the dependence of the triton binding en-ergy on the three-body coupling, D , for the optimal val-ues of the two-body couplings; the shaded bands aroundthe curves show the result of varying the two-body LECswithin their uncertainties for wavefunctions optimised ineach of the three volumes for which there is LQCD data.Figure 6 and Table III show the infinite-volume extrapo-lation of the triton binding energy. As for the two-bodysystems, the extrapolations and couplings reported hereare in close agreement with those of Ref. [20]. □ □ □◦ ◦ ◦ ◦ ◦ ◦ ◦ □□ - - - - - FIG. 6. The volume-dependence of the triton binding en-ergy, compared with the LQCD data. The infinite-volumeextrapolation of the binding energy, computed as describedin the text, is shown in the rightmost sub-panel. V. MATRIX ELEMENTS IN THE STOCHASTICVARIATIONAL METHOD Having determined finite-volume ground state (or inprinciple excited state) wavefunctions in the SVM, thosewavefunctions can be used to evaluate finite-volume ma-trix elements of operators in FVEFT π/ . The transitionmatrix element between an initial state Ψ i and final stateΨ f is given by (cid:104) Ψ f |O| Ψ i (cid:105) (cid:112) (cid:104) Ψ f | Ψ f (cid:105)(cid:104) Ψ i | Ψ i (cid:105) , (38)where O is a generic local EFT operator and bra-ketnotation is used for concision. Here, matrix elements of axial, electromagnetic, scalar,and tensor currents and the unpolarised twist-two oper-ators are studied. In each case, the relevant EFT π/ cur-rents of Sec. II are translated into operators acting on the n -body states. As with the nucleon-nucleon strong inter-actions, two-body contributions to the various currentsare regulated using the Gaussian approach as in Eq. (24)and rendered periodic using Eq. (25). Specifically, eachtwo-body current is implemented as (cid:2)(cid:0) N † ( r i )Σ N ( r i ) (cid:1) (cid:0) N ( r j ) † Σ (cid:48) N ( r j ) (cid:1) + h.c. (cid:3) g Λ ( r ij , L ) , (39)where Σ ( (cid:48) ) denotes a spin-isospin structure and g Λ ( r , L ),defined as in Eq. (24), implements a periodic regulatedform of the δ -function implied in local two-body EFTcurrents. Since the matrix elements that are consid-ered are for zero momentum transfer, the current is Matrix elements of non-local products of operators (such as thoserelevant for double- β decay) can also be approached in the finite-volume SVM as discussed in Sec. VI. integrated over the positions r i,j . For each current, X ∈ { A i,a , J EMi , S a , T ij,a , O ( n ) a } for a ∈ { , , , } , thezero–momentum-projected, regulated form is labelled as X .The evaluation of the relevant matrix elements fac-torises into a spin-isospin calculation that is specific toeach type of operator, and a calculation of the matrix el-ement of the spatial wavefunction. Since all currents thatare considered enter with the spatial dependence deter-mined by the Gaussian regulator function, these latterspatial matrix elements have a common form and aregiven for diagonal matrix elements in state h by h h (Λ , L ) = (cid:82) (cid:81) k d r k (cid:80) i 1= 2 (cid:32) (cid:101) L ,A h np ← d (Λ , L ) (cid:33) , (42) R H A, ≡ g A (cid:104) Ψ H |A , | Ψ H (cid:105)(cid:104) Ψ H | Ψ H (cid:105) = (cid:32) (cid:101) L ,A h H (Λ , L ) (cid:33) , (43)where the spin-flavour structure of the states used to ar-rive at these expressions are given in Sec. III.Figure 7(a) shows the constraints that the LQCD cal-culations [38] of the finite-volume matrix elements in thetwo channels place on the coupling combination (cid:101) L ,A .The consistency between the constraints in the two chan-nels suggests that higher-order terms in the axial cur-rent (two-body operators with derivative insertions [28],or three-body operators) are suppressed as their power-counting would suggest. Note that this consistency isregulator-dependent. Were this to persist for physicalquark masses, it would provide support for approachesto pp -fusion cross-section calculations that use tritium β -decay to constrain the relevant two-body LEC.The values of (cid:101) L ,A determined from each channelare scheme-dependent quantities, but can be combinedwith infinite-volume SVM wavefunctions to determinethe infinite-volume matrix elements. Figure 7(b) showsthe infinite-volume extrapolation for both channels, andthe extrapolated values are given in Table IV below.Analogous analysis of the isoscalar axial matrix ele-ments in the deuteron and He states allows for the de-termination of the two-body counterterm in Eq. (11).Ratios of the isoscalar axial current matrix element inthe deuteron and He states to that in the proton statecan be expressed as R d ; j z =1 A, ≡ − g A, (cid:104) Ψ d ; j z =1 |A , | Ψ d ; j z =1 (cid:105)(cid:104) Ψ d ; j z =1 | Ψ d ; j z =1 (cid:105) = 2 (cid:16) − (cid:101) L ,A h d (Λ , L ) (cid:17) , (44) R H A, ≡ − g A, (cid:104) Ψ H |A , | Ψ H (cid:105)(cid:104) Ψ H | Ψ H (cid:105) = (cid:18) − (cid:101) L ,A h He (Λ , L ) (cid:19) . (45)Figure 8 shows the constraints on (cid:101) L ,A obtained bymatching to the LQCD calculation of Ref. [39], andthe corresponding infinite-volume extrapolation of theLQCD matrix elements. The extrapolated values are alsoreported in Table IV. A mild tension is found between thevalues of (cid:101) L ,A extracted from each matrix element, indi-cating the potential need for higher-order terms in theEFT description. B. Magnetic moments The magnetic moments of light nuclei have been ex-tracted from the linear response of LQCD calculations to a constant background magnetic field oriented in the z -direction [40, 41]. In EFT, these quantities are de-termined by the couplings in Eq. (13). The differencesbetween the magnetic moments of the deuteron, H, and He states and the relevant naive shell-model predictionsin terms of proton and neutron magnetic moments canbe expressed as: δ ˆ µ d ≡ ˆ µ d − (ˆ µ p + ˆ µ n )= 2 M N e (cid:104) Ψ d ; j z =1 |J EM | Ψ d ; j z =1 (cid:105)(cid:104) Ψ d ; j z =1 | Ψ d ; j z =1 (cid:105) − κ = 2 M N L h d (Λ , L ) , (46) δ ˆ µ H ≡ ˆ µ H − ˆ µ n = 2 M N e (cid:104) Ψ H |J EM | Ψ H (cid:105)(cid:104) Ψ H | Ψ H (cid:105) − κ n = M N L + L ) h H (Λ , L ) , (47) δ ˆ µ He ≡ ˆ µ He − ˆ µ p = 2 M N e (cid:104) Ψ He |J EM | Ψ He (cid:105)(cid:104) Ψ He | Ψ He (cid:105) − κ p = − M N L − L ) h He (Λ , L ) , (48)where ˆ µ h is the magnetic moment of hadron h in nat-ural nuclear magnetons e/ M N defined using the nu-cleon mass at the quark masses of the lattice calculations, M N = 1 . δ ˆ µ h for h ∈{ d, H , He } in Ref. [40] are used to constrain the EFTcouplings L , . Since the magnetic moment differencesare dependent of various combinations of the couplings,the constraints take the form of bands in the L – L planeas shown in the figure. All three constraints are seen tobe consistent for a range of values of the couplings. Thisdetermination of the LECs allows for extrapolation of themagnetic moments to infinite volume as shown in Fig. 10,and produces the values shown in Table IV. C. Scalar matrix elements: nuclear σ terms Ratios of the matrix elements of the isoscalar andisovector scalar currents in hadron h to those in the pro-ton can be expressed as R hS, ≡ g S, (cid:104) Ψ h |S | Ψ h (cid:105)(cid:104) Ψ h | Ψ h (cid:105) = (cid:32) A h − f hS, g S, h h (Λ , L ) (cid:33) , (49) Note that in Ref. [40] the magnetic background field does notcouple to sea quarks, so only isovector quantities are calculatedcompletely; the error from this quenching of the magnetic field isignored here. In principle, the isovector np → dγ M1 transitioncan also be used to constrain L , but it is not used in this work. - - - - - - - - - - - - - - - - - - FIG. 9. The dependence of the various magnetic moment differences on the appropriate combinations of the two-bodycounterterms L , is shown in the upper row and the lower left panel. The lower right panel shows the combined constraintsimplied by agreement with the LQCD results of Ref. [40]. The details of the curves and points in the figure are as in Fig 7. ◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦□□ □□◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□ - FIG. 10. The infinite-volume extrapolations of the magneticmoment differences for d , H and He. and R hS, ≡ g S, (cid:104) Ψ h |S | Ψ h (cid:105)(cid:104) Ψ h | Ψ h (cid:105) = (cid:32) T h − f hS, g S, h h (Λ , L ) (cid:33) , (50) where A h denotes the atomic number of the nucleus, T h is its third component of isospin, and f hS, = (cid:101) C + (cid:101) C , h = d (cid:101) C − (cid:101) C , h = pp (cid:101) C − (cid:101) C , h = H , f hS, = (cid:40) (cid:101) C V , h = pp (cid:101) C V , h = H . (51)The quantity (cid:101)(cid:101) C V = (cid:101) C V /g S, is constrained by compar-ing Eq. (55) to LQCD calculations of ratios of isovectorscalar current matrix elements in different nuclei fromRef. [39]. The results of this comparison are shown inFig. 11. Determinations of (cid:101)(cid:101) C V from both the pp and H systems are consistent, with H providing a consider-ably more stringent constraint. The extracted values ofthe LEC are also used to extrapolate the LQCD matrixelements to infinite volume, as shown in the figure andpresented in Table IV.Similarly Fig. 12 compares LQCD calculations of theratio of the scalar isoscalar current matrix element in nu-clei to that in the proton to the expectations of Eq. (49)for the deuteron, diproton and H. The three states pro-vide sufficient information to constrain the LEC ratios (cid:101)(cid:101) C , = (cid:101) C , /g S, and the constrained values are used toextrapolate the LQCD matrix elements to infinite volumeas shown in Fig. 13 and presented in Table IV. As noted3in Sec. II C, the couplings (cid:101) C , that occur in the isoscalarscalar current are related to the Lagrangian counterterms C , in the limit of massless quarks. D. Tensor matrix elements For the tensor current, isoscalar matrix element ratiosto those in the nucleon are given by R hT, ≡ g T, (cid:104) Ψ h |T , | Ψ h (cid:105)(cid:104) Ψ h | Ψ h (cid:105) = (cid:0) S h − f hT, h h (Λ , L ) (cid:1) , (52)for h ∈ { d, H } , where S h is the third component of spinand f hT, = (cid:40)(cid:101) L ,T , h = d ; j z = 1 (cid:101) L ,T , h = H , (53)where (cid:101) L ,T = L ,T /g T, .For the isovector case, the corresponding ratio in H is R H T, ≡ g T, (cid:104) Ψ H |T , | Ψ H (cid:105)(cid:104) Ψ H | Ψ H (cid:105) = (cid:32) (cid:101) L ,T h H (Λ , L ) (cid:33) , (54)in terms of the LEC ratio (cid:101) L ,T = L ,T /g T, Figures 14 and 15 show the comparisons of Eqs. (52)and (54) to the respective LQCD calculations [39]. Inthe isoscalar case, consistency is seen in the values of (cid:101) L ,T that arise from comparison to either d or H matrixelement ratios from LQCD. As for the matrix elements ofother operators considered above, the constrained LECsenable infinite volume extrapolations of the matrix ele-ments which are presented in Table IV. E. Twist-two operators: the quark momentumfraction The first moment of the isovector unpolarised partondistribution has been computed in LQCD for nuclei with A ∈ { , } in Ref. [42], and correspond to the differencein the longitudinal momentum fractions carried by upand down quarks. In the finite-volume SVM, the ma-trix elements of the operators in Eq. (20) contain no spinstructure and so the finite-volume matching and extrapo-lation is analogous to that on the isovector scalar currentabove. In particular, ratios of the isovector matrix ele-ment in hadron h to the proton matrix elements are givenby R h O n , ≡ A h ( Z h − N h ) (cid:104) x n (cid:105) (cid:104) Ψ h |O ( n )3 | Ψ h (cid:105)(cid:104) Ψ h | Ψ h (cid:105) = (cid:18) α n, ( Z h − N h ) (cid:104) x n (cid:105) h h (Λ , L ) (cid:19) , (55) Quantity O State h ( → h (cid:48) ) O ( L = 4 . O ( L = ∞ ) δ ˆ µ ( h ) [40] d ( j z = ± 1) 0.011(80) 0.012(89) He -0.34(10) -0.35(10) H 0.45(16) 0.46(17) R hA, [38] H 0.979(10) 0.978(11) np ( → d ) 1.978(31) 1.975(36) R hA, [39] d . H 0 . R hT, [39] H 1 . np ( → d ) – -1.415(2) R hT, [39] d . H 0 . R hS, [39] pp . H 0 . R hS, [39] pp . d . H 2 . R h O n =1 , [42] pp . H 1 . m π = 806 MeV in a L = 4 . R d ← npT, , the LEC determinedfrom the H channel is used to make a prediction, as thereare no LQCD results available. for h ∈ { pp, H } .LQCD data in a calculation with L = 4 . α , , as shown in Fig. 16. The ex-trapolated infinite-volume matrix elements are reportedin Table IV. VI. DISCUSSION Finite-volume pionless effective field theory imple-mented through the stochastic variational method hasbeen used to extrapolate EFT wavefunctions and matrixelements for A ∈ { , } nuclei, matched to LQCD calcu-4 - - (a) ◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦□□ □□ ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□ (b) FIG. 11. (a) The dependence of the ratios of the pp (upper) and H (lower) isovector scalar matrix elements to that of theproton on the LEC ratio (cid:101)(cid:101) C V . (b) The infinite-volume extrapolation of these ratios after constraining the LEC ratio. Thedetails of the curves and points in the figure are as in Fig 7. - - - - - - - FIG. 12. The isoscalar scalar current matrix element ratios for d , pp and H as a function of the relevant combinations of LECratios (cid:101)(cid:101) C , . The lower right panel shows the resulting constraints on these LEC ratios. The details of the curves and points inthe figure are as in Fig 7. ◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦□□ □□ ◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦□□ □□ ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□ FIG. 13. The infinite-volume extrapolations of the isoscalarscalar current matrix element ratios for d , pp and H. lations in a finite lattice volume, to infinite volume. Thisnumerical approach can effectively describe bound-statesystems, can cleanly reproduce scattering states, and fur-thermore exhibits volume scaling that is consistent withthe predictions of the L¨uscher approach to high accuracy.To some degree, the method circumvents the complexi-ties of analytic approaches generalising that of L¨uscherfor two-body systems to larger number of particles. How-ever, as the atomic number of the system increases, thefinite-volume SVM scales relatively poorly, and it doesnot seem practical to extend much beyond A = 4 sys-tems or more than three-body interactions. Given thatthe use of the finite-volume aspect of the method is tai-lored to match LQCD calculations, for which increasing A is also costly, this is perhaps not a significant limi-tation: the finite-volume SVM can be used to determineEFT counterterms which can then be used in the infinite-volume SVM (or other many-body methods) to performcalculations for larger nuclei.For all of the matrix elements studied in this work,which include isoscalar and isovector scalar, axial, and tensor matrix elements, as well as magnetic moments andthe isovector longitudinal momentum fraction, it is foundthat for the large quark masses used in the LQCD calcu-lations, the lattice volume of L =4.5 fm as used for thecalculations is large enough that there are essentially nofinite-volume corrections. At lighter quark masses, how-ever, one might anticipate that larger lattice volumes willbe required to achieve the same behaviour. For almostall of the matrix elements investigated in this work, it isalso notable that although the constraints from LQCDcalculations of three-body systems are typically tighter,the EFT π/ LECs determined from the LQCD calculationsof two and three-body systems are consistent to withinone standard deviation (with the notable exception of theisoscalar axial LEC L ,A in Eq. (11), for which there is aslight tension). This indicates that higher-order terms inthe relevant currents are suppressed in the exponentialregulator scheme used in the SVM at these values of thequark masses. In some cases, this suppression has par-ticular consequence; for example, for the isovector axialmatrix elements, the LEC L ,A determined from mea-surements of tritium β -decay is used in calculations ofthe pp -fusion cross-section [43]. With the LECs deter-mined from LQCD calculations, predictions can be madefor other quantities for which there are no LQCD results;in Table IV, as an example, the np ← d tensor transitionmatrix element is predicted from the LEC determinedfrom the triton matrix element of the same current.Since the finite-volume SVM provides representationsof low-lying excited states as well as the ground statesthat have been the focus of this work, it can also beused to match matrix elements of second-order currentinsertions such as in double- β decay. In such processes, asum over excited nuclear states occurs for times betweenthose of the two currents, with the matrix elements ofinterest being (cid:88) n (cid:104) Ψ f |J | Ψ n (cid:105)(cid:104) Ψ n |J | Ψ i (cid:105) E − E n . (56)In the finite-volume EFT context this corresponds to in-clusion of the discrete states in principle up to energy-scale of the EFT cutoff. In order to accurately representthese contributions, care must be taken that these statesare equivalently well optimised. This will require signif-icant numerical effort, but is likely to be feasible. Ulti-mately, the finite-volume SVM appears to be a powerfultool to capitalise on LQCD calculations of systems withsmall A , which are approaching a novel era of systematiccontrol. ACKNOWLEDGMENTS We are grateful to J-W. Chen, Z. Davoudi, M. Illa,A. Parre˜no, M. J. Savage, and M. Wagman for discus-sions. The authors are supported in part by the U.S. De-partment of Energy, Office of Science, Office of Nuclear6 - (a) ◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦□□ □□ ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□ (b) FIG. 14. (a) The dependence of the d (upper) and H (lower) isoscalar tensor matrix elements on (cid:101) L ,T . (b) The dependenceof these matrix elements on the lattice extent, L . The details of the curves and points in the figure are as in Fig 7. - - - (a) ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□ (b) FIG. 15. (a) The dependence of the H isovector tensor transition matrix element on (cid:101) L ,T . (b) The dependence of this matrixelement on the lattice extent, L . The details of the curves and points in the figure are as in Fig 7. Physics under grant Contract Number DE-SC0011090and by the Carl G and Shirley Sontheimer ResearchFund. WD is also supported within the framework ofthe TMD Topical Collaboration of the U.S. Departmentof Energy, Office of Science, Office of Nuclear Physics,and by the SciDAC4 award DE-SC0018121. PES is ad-ditionally supported by the National Science Foundationunder EAGER grant 2035015, by the U.S. DOE EarlyCareer Award DE-SC0021006, and by a NEC research award. Appendix A: Alternate form for currents Currents are derived by first considering the relativis-tic form at the quark level, matching onto relativisticnucleon operators with the same C, P and T properties,and then performing a nonrelativistic reduction. In the7 - - (a) ◦ ◦ ◦ ◦◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦□□ □□ ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦□□ □□ (b) FIG. 16. (a) The dependence of the isovector twist-two matrix elements in pp and H on the two-body LEC, α , . (b) Thedependence of the matrix element on the lattice extent, L . The details of the curves and points in the figure are as in Fig 7. main text, the various currents were written using theprojectors in Eq. (4). They can also be written in termsof Pauli matrices in spin and isospin as follows, whereeach expression is given only to the order used in thiswork: A i,a = g A N † τ a σ i N − L ,A (cid:0) N † σ i N (cid:1) (cid:0) N † τ a N (cid:1) , (A1) A i, = − g A, N † σ i N + L ,A (cid:0) N † σ i N (cid:1) (cid:0) N † N (cid:1) , (A2) S = g S, N † N − (cid:101) C (cid:0) N † N (cid:1) (cid:0) N † N (cid:1) , − (cid:101) C (cid:0) N † σ i N (cid:1) (cid:0) N † σ i N (cid:1) , (A3) S a = g S, N † τ a N − (cid:101) C V (cid:0) N † τ a N (cid:1) (cid:0) N † N (cid:1) , (A4) T ij, = g T, (cid:15) ijk N † σ k N − L ,T (cid:15) ijk (cid:0) N † σ k N (cid:1) (cid:0) N † N (cid:1) , (A5) T ij ; a = g T, (cid:15) ijk N † τ a σ k N − L ,T (cid:15) ijk (cid:0) N † σ k N (cid:1) (cid:0) N † τ a N (cid:1) . (A6) Note that each of these terms is Hermitian so no Her-mitian conjugation is implied. For the scalar currents, (cid:101) C T = (cid:101) C − (cid:101) C and (cid:101) C S = (cid:101) C + (cid:101) C in Eq. (15), and { (cid:101) C , (cid:101) C } = { C S , C T } of Ref [31].Similarly, the two body part of Eq. (12) can be writtenas L , EM = − e L (cid:0) N † σ · B N (cid:1) (cid:0) N † τ N (cid:1) + e L (cid:0) N † σ · B N (cid:1) (cid:0) N † N (cid:1) . (A7) Appendix B: Scattering states for free particles andweak interactions The correlated shifted Gaussian basis is able to accu-rately describe low-energy finite-volume scattering statesas well as localised bound states. To demonstrate this,SVM approximations for non-interacting two particlestates are studied in this appendix.Figure 17(a) shows the energy eigenvalues obtained fora free two-nucleon system using wavefunctions approxi-mated using the shifted correlated Gaussian basis func-tions. The eigenvalues are shown in units of 2 M N L / π ,in which case the expectation is an integer-spaced spec-trum. States are symmetric under particle interchangeand should be ordered in terms of the sum of the squaredmomenta, N = | p | + | p | | p | , and should have degen-8 □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□△△△△△△△ △△△△△△△△△ △△△△△△△△△△△△△△△△△△△△△△△△○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○ □ △ ○ (a) □ □ □ □ □ □ □ □ (b) FIG. 17. (a) Energy eigenvalues, λ n for the free two-nucleon system obtained using the SVM in three different volumes, plottedin units of 2 M N L / π . The coloured regions indicate the expected multiplicity of eigenvalues. (b) Energy eigenvalues, λ fora weakly repulsive interaction obtained using the SVM in multiple volumes. The solid and dashed curves correspond to fitsusing Eq. (B1) with either { a, c } or { a } as fit parameters, respectively. eracies 1 , , , . . . corresponding to { p = (0 , , , p =(0 , , } for N = 0, { p = (1 , , , p = (0 , , } and permutations and sign changes for N = 1, { p =(1 , , , p = (0 , , } or { p = (1 , , , p = (0 , , } and permutations and sign changes for N = 2, and soon. As can be seen from the figure, the SVM is able tocleanly reproduce the low-energy part of the spectrum,including its degeneracies; with further numerical effortthis can be extended further. Similarly the free three-particle low-energy spectrum is well reproduced.The large L asymptotic behaviour of the ground stateis given by an expansion of the two-particle quantisationcondition derived by L¨uscher [24], namely λ = 4 πaM N L (cid:20) − c aL + c (cid:16) aL (cid:17) + . . . (cid:21) + O (cid:0) L − (cid:1) , (B1) where a is the scattering length and the geometric coef-ficients are c = − . c = +6 . C S = 31 MeV fm and C T = 0 in Eq. (3)) as a function of volume. Alsoshown are fits to this data using Eq. (B1) with the scat-tering length as a free parameter, as well as fits treatingboth the scattering length and the geometric constant c as free parameters. The latter fit returns a value of c fit1 = − . Appendix C: Matrix element formulae In this section, explicit formulae for the wavefunction integrals in Eqs. (35) and (36) are provided. In the expressionsbelow, Ψ sym L ( A i , B i , d i ; x ) is a symmetric n -body Gaussian wavefunction term defined in Eq. (31). The normalisationintegral of Eq. (35), corresponding to a cross-term between such n -body Gaussian wavefunction terms labelled by thesubscripts i and j respectively, can be expressed as[ N ] ij ≡ (cid:90) Ψ sym L ( A i , B i , d i ; x ) ∗ Ψ sym L ( A j , B j , d j ; x ) d x = (cid:88) P , P (cid:48) (cid:89) α ∈{ x,y,z } (cid:115) (2 π ) n Det[ C ( α ) i P ; j P (cid:48) ] | b ( α ) |≤ ˜ b (cid:88) b ( α ) exp (cid:34) − 12 Ω ( α ) i P ; j P (cid:48) (cid:35) , (C1)where Ξ ( α ) i P ; j P (cid:48) = LA ( α ) i P · b ( α ) + B ( α ) i P · ( L b ( α ) + d ( α ) i P ) + B j P (cid:48) · d ( α ) j P (cid:48) , (C2)Ω ( α ) i P ; j P (cid:48) = ( L b ( α ) ) · A ( α ) i P · ( L b ( α ) ) + ( L b ( α ) + d ( α ) i P ) · B ( α ) i P · ( L b ( α ) + d ( α ) i P ) + d ( α ) j P (cid:48) · B ( α ) j P (cid:48) · d ( α ) j P (cid:48) − Ξ ( α ) · [ C ( α ) i P ; j P (cid:48) ] − · Ξ ( α ) , (C3)9and C ( α ) i P ; j P (cid:48) = A ( α ) i P + A ( α ) j P (cid:48) + B ( α ) i P + B ( α ) j P (cid:48) , (C4)and where superscripts ( α ) and subscripts P on the wavefunction parameters A , B , and d extract the componentsof the parameters corresponding to the α direction, and permute the parameters for each of the n particles by thepermutation P , respectively. In the above equations, b ( α ) is a length- n vector; summing over all n -vectors correspondsto enforcing periodicity. In practice the infinite sum is truncated to vectors with a maximum norm ˜ b , and in thenumerical calculations in this work the cut is initially taken to be ˜ b = 3 for each term, and is iteratively increaseduntil the fractional change in the total summed expression from adding an additional term is less than 10 − .The spatial integrals involved in the Hamiltonian matrix elements are separated into the kinetic and two andthree-body potential terms as H = (cid:104) χ h | K + ( C + C σ · σ ) V + D V | χ h (cid:105) , (C5)where the spatial integrals for each of K , V , and V can be performed independently. The integral for the two-bodypotential term, again for n -body Gaussian wavefunction terms labelled by i and j , is[ V ] ij ≡ n (cid:88) a
12 Ω ( α ) i P ; j P (cid:48) (cid:35) × ˜ q (cid:88) q ( α ) = − ˜ q exp (cid:34) − ρ (cid:101) C ( α ) i P ; j P (cid:48) (cid:101) C ( α ) i P ; j P (cid:48) + 2 ρ (cid:16) [( C ( α ) i P ; j P (cid:48) ) − · Ξ ( α ) ] a − [( C ( α ) i P ; j P (cid:48) ) − · Ξ ( α ) ] b − Lq ( α ) (cid:17) (cid:35) , (C6)where the constant ρ = r = Λ is a re-scaling of the Gaussian regulator parameter, and (cid:101) C ( α ) i P ; j P (cid:48) = (cid:16) [ C ( α ) i P ; j P (cid:48) ] − aa + [ C ( α ) i P ; j P (cid:48) ] − bb − [ C ( α ) i P ; j P (cid:48) ] − ab − [ C ( α ) i P ; j P (cid:48) ] − ba (cid:17) − , (C7)where [ V ] a denotes the a th component of a vector V , and [ M ] − ab denotes the ( a, b ) component of the matrix M − .In this expression q ( α ) is an integer; combined with the sum over b ( α ) , summing over all values of q ( α ) corresponds toenforcing periodicity. In practice the infinite sum is truncated to integers with absolute value less than ˜ q = 40; thefractional change in the total summed expression from adding an additional term is less than 10 − .The relevant integral for the three-body potential term, for n -body Gaussian wavefunction terms labelled by i and j , can be expressed as[ V ] ij ≡ cyc (cid:88) a (cid:54) = b (cid:54) = c (cid:90) Ψ sym L ( A i , B i , d i ; x ) ∗ g Λ ( x a − x b , L ) g Λ ( x b − x c , L )Ψ sym L ( A j , B j , d j ; x ) d x = (cid:88) P , P (cid:48) cyc (cid:88) a (cid:54) = b (cid:54) = c (cid:89) α ∈{ x,y,z } (cid:115) (2 π ) n Det[ (cid:98) C ( α ) i P ; j P (cid:48) ] exp (cid:20) − (cid:16) d ( α ) i P · B ( α ) i P · d ( α ) i P + d ( α ) j P (cid:48) · B ( α ) j P (cid:48) · d ( α ) j P (cid:48) (cid:17)(cid:21) × | b ( α ) |≤ ˜ b (cid:88) b ( α ) exp (cid:20) − (cid:16) ( L b ( α ) ) · ( A ( α ) i P + B ( α ) i P ) · ( L b ( α ) ) + 2 d ( α ) i P · B ( α ) i P · ( L b ( α ) ) − Ξ ( α ) · [ (cid:98) C ( α ) i P ; j P (cid:48) ] − · Ξ ( α ) (cid:17)(cid:21) × ˜ q (cid:88) q ( α ) = − ˜ q exp (cid:34) − L r q ( α )2 + q ( α )2 L r P [ a,b ] v · [ (cid:98) C ( α ) i P ; j P (cid:48) ] − · P [ a,b ] v + q ( α ) Lr Ξ ( α ) · [ (cid:98) C ( α ) i P ; j P (cid:48) ] − · P [ a,b ] v (cid:35) × ˜ q (cid:88) t ( α ) = − ˜ q exp (cid:34) − L r t ( α )2 + t ( α )2 L r ∗ P [ b,c ] v · [ (cid:98) C ( α ) i P ; j P (cid:48) ] − · P [ b,c ] v + t ( α ) Lr Ξ ( α ) · [ (cid:98) C ( α ) i P ; j P (cid:48) ] − · P [ b,c ] v (cid:35) × exp (cid:20) t ( α ) q ( α ) L r P [ b,c ] v · [ (cid:98) C ( α ) i P ; j P (cid:48) ] − · P [ a,b ] v ) (cid:21) , (C8)0where (cid:80) cyc a (cid:54) = b (cid:54) = c indicates a sum over cyclic permutations of particles a , b and c , (cid:98) C ( α ) i P ; j P (cid:48) = C ( α ) i P ; j P (cid:48) + 1 r (cid:16) P [ a,b ] m + P [ b,c ] m (cid:17) , (C9)and vector and matrix projectors are defined component-wise as[ P [ a,b ] v ] c = δ ac − δ bc , (C10)[ P [ a,b ] m ] cd = δ cd ( δ ac + δ bc ) − δ ac δ bd − δ ad δ bc . (C11)As in Eqs. (C1) and (C6), the sums over b ( α ) , q ( α ) , and t ( α ) together enforce periodicity. In numerical evaluations ofEq. (C8), the same cut procedure for fixing ˜ b and ˜ q are used as for the evaluations of the previous expressions.Finally, the integral for the kinetic operator, for n -body Gaussian wavefunction terms labelled by i and j , is[ K ] ij ≡ − M N n (cid:88) a =1 (cid:90) Ψ sym L ( A i , B i , d i ; x ) ∗ ∇ a Ψ sym L ( A j , B j , d j ; x ) d x = 12 M N (cid:88) P , P (cid:48) (cid:88) α ∈{ x,y,z } (cid:115) (2 π ) n Det[ C ( α ) i P ; j P (cid:48) ] | b ( α ) |≤ ˜ b (cid:88) b ( α ) Θ ( α ) i P ; j P (cid:48) exp (cid:34) − 12 Ω ( α ) i P ; j P (cid:48) (cid:35) × β (cid:54) = α (cid:89) β ∈{ x,y,z } (cid:115) (2 π ) n Det[ C ( β ) i P ; j P (cid:48) ] | b ( β ) |≤ ˜ b (cid:88) b ( β ) exp (cid:34) − 12 Ω ( β ) i P ; j P (cid:48) (cid:35) , (C12)where (cid:126) = 1 is used andΘ ( α ) i P ; j P (cid:48) =Tr (cid:20) [( A ( α ) i P + B ( α ) i P ) · (cid:104) C ( α ) i P ; j P (cid:48) (cid:105) − · ( A ( α ) j P (cid:48) + B ( α ) j P (cid:48) ) (cid:21) − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( A ( α ) j P (cid:48) + B ( α ) j P (cid:48) ) · (cid:104) C ( α ) i P ; j P (cid:48) (cid:105) − · (cid:16) B ( α ) i P · ( L b ( α ) + d ( α ) i P ) + A ( α ) i P · ( L b ( α ) ) (cid:17) − ( A ( α ) i P + B ( α ) i P ) · (cid:104) C ( α ) i P ; j P (cid:48) (cid:105) − · B ( α ) j P (cid:48) · d ( α ) j P (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (C13) [1] D. B. Kaplan, M. J. Savage, and M. B. Wise, Phys. Lett. B424 , 390 (1998), arXiv:nucl-th/9801034 [nucl-th].[2] D. B. Kaplan, M. J. Savage, and M. B. Wise, Nucl. Phys. B534 , 329 (1998), arXiv:nucl-th/9802075 [nucl-th].[3] U. van Kolck, Nucl. Phys. A645 , 273 (1999), arXiv:nucl-th/9808007 [nucl-th].[4] J.-W. Chen, G. Rupak, and M. J. Savage, Nucl. Phys. A653 , 386 (1999), arXiv:nucl-th/9902056 [nucl-th].[5] P. F. Bedaque and U. van Kolck, Ann. Rev. Nucl. Part.Sci. , 339 (2002), arXiv:nucl-th/0203055.[6] H. W. Hammer, S. K¨onig, and U. van Kolck, Rev. Mod.Phys. , 025004 (2020), arXiv:1906.12122 [nucl-th].[7] N. Barnea, L. Contessi, D. Gazit, F. Pederiva, andU. van Kolck, Phys. Rev. Lett. , 052501 (2015),arXiv:1311.4966 [nucl-th].[8] S. Beane, E. Chang, S. Cohen, W. Detmold, H.-W.Lin, T. Luu, K. Orginos, A. Parreno, M. Savage, andA. Walker-Loud, Phys. Rev. Lett. , 172001 (2012),arXiv:1204.3606 [hep-lat].[9] J. Kirscher, N. Barnea, D. Gazit, F. Pederiva, and U. van Kolck, Phys. Rev. C , 054002 (2015), arXiv:1506.09048[nucl-th].[10] L. Contessi, A. Lovato, F. Pederiva, A. Roggero,J. Kirscher, and U. van Kolck, Phys. Lett. B , 839(2017), arXiv:1701.06516 [nucl-th].[11] A. Bansal, S. Binder, A. Ekstr¨om, G. Hagen, G. Jansen,and T. Papenbrock, Phys. Rev. C , 054301 (2018),arXiv:1712.10246 [nucl-th].[12] J. Kirscher, E. Pazy, J. Drachman, and N. Barnea, Phys.Rev. C , 024001 (2017), arXiv:1702.07268 [nucl-th].[13] W. Detmold and M. J. Savage, Nucl. Phys. A743 , 170(2004), arXiv:hep-lat/0403005 [hep-lat].[14] R. A. Brice˜no and Z. Davoudi, Phys. Rev. D88 , 94507(2013), arXiv:1204.1110 [hep-lat].[15] P. E. Shanahan, B. C. Tiburzi, M. L. Wagman, F. Win-ter, E. Chang, Z. Davoudi, W. Detmold, K. Orginos,and M. J. Savage, Phys. Rev. Lett. , 62003 (2017),arXiv:1701.03456 [hep-lat].[16] B. C. Tiburzi, M. L. Wagman, F. Winter, E. Chang,Z. Davoudi, W. Detmold, K. Orginos, M. J. Savage, and P. E. Shanahan, Phys. Rev. D96 , 54505 (2017),arXiv:1702.02929 [hep-lat].[17] Z. Davoudi and S. V. Kadam, (2020), arXiv:2012.02083[hep-lat].[18] Z. Davoudi and S. V. Kadam, Phys. Rev. D , 114521(2020), arXiv:2007.15542 [hep-lat].[19] R. A. Brice˜no, Z. Davoudi, M. T. Hansen, M. R.Schindler, and A. Baroni, Phys. Rev. D , 014509(2020), arXiv:1911.04036 [hep-lat].[20] M. Eliyahu, B. Bazak, and N. Barnea, (2019),arXiv:1912.07017 [nucl-th].[21] K. Varga and Y. Suzuki, Phys. Rev. C , 2885 (1995),arXiv:nucl-th/9508023.[22] J. Mitroy, S. Bubin, W. Horiuchi, Y. Suzuki, L. Adamow-icz, W. Cencek, K. Szalewicz, J. Komasa, D. Blume, andK. Varga, Rev. Mod. Phys. , 693 (2013).[23] S. R. Beane, E. Chang, S. D. Cohen, W. Detmold, H. W.Lin, T. C. Luu, K. Orginos, A. Parre˜no, M. J. Savage,and A. Walker-Loud (Nplqcd), Phys. Rev. D87 , 34506(2013), arXiv:1206.5219 [hep-lat].[24] M. Luscher, Commun. Math. Phys. , 153 (1986).[25] T. Mehen, I. W. Stewart, and M. B. Wise, Phys. Rev.Lett. , 931 (1999), arXiv:hep-ph/9902370.[26] P. F. Bedaque, H. Hammer, and U. van Kolck, Phys.Rev. Lett. , 463 (1999), arXiv:nucl-th/9809025.[27] P. F. Bedaque, H. Hammer, and U. van Kolck, Nucl.Phys. A , 444 (1999), arXiv:nucl-th/9811046.[28] M. Butler and J.-W. Chen, Phys. Lett. B520 , 87 (2001),arXiv:nucl-th/0101017 [nucl-th].[29] M. Butler and J.-W. Chen, Nucl. Phys. A675 , 575(2000), arXiv:nucl-th/9905059 [nucl-th].[30] G. Rupak, Nucl. Phys. A , 405 (2000), arXiv:nucl-th/9911018.[31] H. Krebs, E. Epelbaum, and U.-G. Meißner, Eur. Phys. J. A , 240 (2020), arXiv:2005.07433 [nucl-th].[32] J.-W. Chen and W. Detmold, Phys. Lett. B625 , 165(2005), arXiv:hep-ph/0412119 [hep-ph].[33] S. R. Beane and M. J. Savage, Nucl. Phys. A , 259(2005), arXiv:nucl-th/0412025.[34] Y. Suzuki and K. Varga, Stochastic variational approachto quantum-mechanical few body problems , Vol. 54 (1998)pp. 1–310.[35] V. Lensky, M. C. Birse, and N. R. Walet, Phys. Rev. C , 034003 (2016), arXiv:1605.03898 [nucl-th].[36] X. Y. Yin and D. Blume, Phys. Rev. A , 063609 (2013).[37] D. B. Kaplan, M. J. Savage, and M. B. Wise, Phys. Rev.C , 617 (1999), arXiv:nucl-th/9804032.[38] M. J. Savage, P. E. Shanahan, B. C. Tiburzi, M. L. Wag-man, F. Winter, S. R. Beane, E. Chang, Z. Davoudi,W. Detmold, and K. Orginos, Phys. Rev. Lett. ,062002 (2017), arXiv:1610.04545 [hep-lat].[39] E. Chang, Z. Davoudi, W. Detmold, A. S. Gambhir,K. Orginos, M. J. Savage, P. E. Shanahan, M. L. Wag-man, and F. Winter (NPLQCD), Phys. Rev. Lett. ,152002 (2018), arXiv:1712.03221 [hep-lat].[40] S. Beane, E. Chang, S. Cohen, W. Detmold, H. Lin,K. Orginos, A. Parreno, M. Savage, and B. Tiburzi,Phys. Rev. Lett. , 252001 (2014), arXiv:1409.3556[hep-lat].[41] S. R. Beane, E. Chang, W. Detmold, K. Orginos,A. Parre˜no, M. J. Savage, and B. C. Tiburzi (NPLQCD),Phys. Rev. Lett. , 132001 (2015), arXiv:1505.02422[hep-lat].[42] W. Detmold, M. Illa, D. J. Murphy, P. Oare, K. Orginos,P. E. Shanahan, M. L. Wagman, and F. Winter, (2020),arXiv:2009.05522 [hep-lat].[43] H. De-Leon, L. Platter, and D. Gazit, Phys. Rev. C100