Over the last 30 years, effective field theories have revolutionized nuclear physics, systematizing the study of nucleon-nucleon interactions and the properties of light nuclei. Pionless effective field theory for few-nucleon systems (EFT π ) in particular, which focuses on momenta below the pion mass, has emerged as a powerful tool with which to understand low-energy nuclear processes in many contexts [1][2][3][4][5] (see Ref. [6] for a recent review). Notably, in addition to its use in the analysis of experimental data, EFT π has found a key role in the analysis of lattice Quantum Chromodynamics (LQCD) calculations of nuclear systems, providing a direct bridge between QCD and nuclear physics.
For example, in Ref. [7], LQCD calculations of A ∈ f2; 3g nuclei at heavier-than-physical quark masses [8] are matched to auxiliary field diffusion Monte Carlo calculations with EFT π interactions to constrain the two-and three-body counterterms of the EFT; these were then used to make predictions for larger nuclei of atomic number A ≤ 6. Further developments were presented in Ref. [9], and this approach was extended to the next order in the EFT, and to still larger nuclei, in Refs. [10,11]. Reference [12] presents studies of the quark-mass dependence of the magnetic moments and polarizabilities of A ∈ f2; 3g nuclear systems, with both experimental results and LQCD calculations at larger-than-physical values of the pion mass used to constrain the counterterms of EFT π . Similarly, an early application of finite-volume EFT π to electroweak matrix elements was presented in Ref. [13] and extended in Ref. [14], and Ref. [15] implemented numerical studies of two-and three-body systems in discretized EFT π . Furthermore, EFT π provided a powerful approach to analyzing LQCD calculations of second order weak processes [16][17][18][19][20].
Recently, Eliyahu et al. [21] have taken the next steps in this approach and used EFT π implemented via the stochastic variational method (SVM) [22,23] in a finite cubic volume to analyze the binding energies of atomic number A ∈ f2; 3g systems in that finite volume. Since the effects of a finite volume manifest in the infrared domain, they can be captured in low-energy effective descriptions of QCD such as that provided by EFT π . 1 Since LQCD calculations are typically performed in multiple volumes to enable an infinite-volume extrapolation, performing the EFT in the same volumes also maximizes the constraining power of the LQCD results. As input, Ref. [21] utilizes the NPLQCD Collaboration's LQCD computations of the ground-state energies and finite-volume energy shifts of these systems in three different lattice volumes with spatial extents L ∈ f3.4; 4.5; 6.7g fm at unphysical quark masses corresponding to the SUð3Þ f flavor-symmetric point where the up, down, and strange quark masses are degenerate and correspond to a pion mass m π ¼ m K ¼ 806 MeV [25]. Using EFT π in the same volumes to determine counterterms, the binding energies were extrapolated to infinite volume. As this exemplifies, the finite-volume pionless effective field theory (FVEFT π ) approach provides a powerful alternative to Lüscher's method [26], in which finite-volume energies from LQCD calculations are used to determine infinite-volume scattering phase shifts and bound-state energies. While Lüscher's method and its generalizations are model independent, the existing formalism is limited to two-and three-particle [27] energies and matrix elements [28][29][30]. 2 The matching of LQCD results to EFT π , on the other hand, requires an underlying EFT but can be applied to any system that the EFT can address.
In this work, the application of FVEFT π is extended to nuclear matrix elements for the first time. After a brief summary of the relevant EFT π Lagrangians in Sec. II, the SVM is introduced in Sec. III. Section IV presents the results of tuning the relevant two-and three-body counterterms of FVEFT to reproduce the ground-state energies of A ∈ f2; 3g nuclei computed in LQCD, paralleling the analysis of Ref. [21]. Having determined these counterterms, Sec. V presents the tuning of the counterterms describing the interactions of external currents to reproduce LQCD calculations of nuclear matrix elements, thereby enabling an extrapolation of the finite-volume matrix elements to infinite volume. In particular, LQCD calculations of scalar, tensor, and axial matrix elements, as well as magnetic moments and isovector longitudinal momentum fractions, of A ∈ f2; 3g states are investigated. To conclude, the outlook for these calculations is discussed in Sec. VI.
The pionless EFT Lagrangian describing the low-energy interactions of nucleons is given by
where
contains the single-nucleon kinetic operator expanded in the nonrelativistic limit. Here, N represents the nucleon field, M N is the nucleon mass, and the ellipsis denotes higher-order terms. The leading-order two-nucleon interactions enter as
where
are projectors onto spin-triplet and spin-singlet twonucleon states, respectively, and C S and C T are the relevant two-body low-energy constants (LECs). Here, σ k and τ a are the Pauli matrices acting in spin and isospin space, respectively. Equation ( 3) can be reexpressed in a different basis as [40]
where
Three-body interactions naively enter at higher order but must be promoted to leading order as argued in Refs. [41,42], and the relevant contribution to the Lagrangian is
where D 0 is the leading-order three-body LEC.
The weak decays and interactions of nuclear states arise, after integrating out the weak gauge boson, through the effective Lagrangian (valid for energies E ≪ M W )
where G F is the Fermi constant, l μ þ involves a charged lepton and neutrino, and the hadronic weak current can be expressed in terms of the vector (V μ ) and axial-vector (A μ ) currents as
with J AE μ ¼ J 1;μ AE iJ 2;μ . The isovector axial-vector current in EFT π is given by [43] foot_1
where the ellipsis denotes higher-order terms. In this expression, g A is the nucleon axial charge, and the term proportional to the two-nucleon LEC L 1;A provides the next-to-leading-order corrections to the pp → de þ ν fusion process, for example. The isoscalar axial current is similarly given by [44] foot_2
where the isoscalar nucleon axial charge is g A;0 and the LEC L 2;A enters at the same order.
In the presence of an external electromagnetic (EM) field, the Lagrangians defined above are modified such that the derivatives are replaced by EM-covariant derivatives
where A μ is the vector potential and Q is the electric charge operator and also by the addition of terms depending on the magnetic field B,foot_3
where the isoscalar and isovector currents coupling to the magnetic field are
where L 1 and L 2 are two-body LECs and
are the isoscalar and isovector nucleon magnetic moments. Note that electric field contributions and EM three-body interactions enter at higher order.
The isovector and isoscalar scalar currents that arise from Higgs couplings and from potential dark matter interactions are given by S 0 ¼ g S;0 N † N -CS ðN T P i NÞ † ðN T P i NÞ;
Here, g S;0 and g S;3 are the isoscalar and isovector one-body LECs that are related to the nucleon σ terms. As discussed in Ref. [46], the two-body terms in the isoscalar scalar current are related to the corresponding terms in the strong Lagrangian, Eq. ( 3). In particular, the LECs CS;T are the quark-mass-independent pieces of the Lagrangian couplings C S;T .
For the isoscalar and isovector antisymmetric tensor currents, the relevant EFT π expressions are
where the one-and two-body isoscalar (isovector) tensor LECs are g T;0 and L 2;T (g T;3 and L 1;T ).
The unpolarized twist-2 operators that define moments of parton distributions enter in EFT π as [47,48]
where the explicit Lorentz indices of the operators are suppressed on the left-hand side of the definitions for compactness of notation. Subleading contributions from terms involving derivatives are suppressed in the power counting, and contributions from additional two-body terms that are not Wigner SU(4) symmetric are suppressed in the large N c limit. The subscripts 0 and 3 denote the isoscalar and isovector combinations. Note that the isoscalar contributions arise from matching to both quark and gluon matrix elements (which mix under QCD renormalization) and will give rise to finite-volume effects that are the same in both cases.
In order to address finite-volume effects in few-nucleon systems in EFT, few-body wave functions must be determined subject to the EFT interactions and the given boundary conditions. There are multiple many-body approaches that could be pursued for this task. Two approaches that have been successfully applied are solving the three-dimensional finite-volume Schrödinger equation via discretization [8] and the SVM for two-and three-body systems [21]. The former approach works well for twobody interactions and was effectively used in Ref. [8] to analyze hyperon-nucleon interactions where the effective range was not small compared to the spatial extent of the finite volume, L, and as such the more direct Lüscher method [26] could not be applied. However, this coordinate space-based approach scales poorly to larger systems.
The SVM was introduced in nuclear physics [22] as a way to sample the possible spatial, spin, and isospin wave functions for an A-nucleon system in a space that is impractically large for an exhaustive approach; see Refs. [23,49] for reviews. This approach, detailed below (and applied to EFT π in Ref. [50], for example) involves the construction of a wave function by sequential proposals of new stochastically generated terms and the optimization of the linear coefficients of the terms by solving the generalized eigenvalue problem of the variational method. The SVM for a finite volume was first introduced in Ref. [51] in which systems of bosons in periodic cubic potentials were considered. Periodicity is imposed on the wave functions by considering all periodic copies of the infinite-volume potential. The method was first used for nuclei in Ref. [21], and a similar approach is used here.
The n-particle nonrelativistic Hamiltonian that corresponds to the EFT π Lagrangian of Eq. ( 1) is
where i; j; k ∈ f1; …; ng label the particle, r ij ¼ r ir j is the displacement between particles i and j, and ∇ 2 i denotes the Laplacian for particle i. The two-and three-particle potentials are
and
where the interactions have been regulated using Gaussian smearing. This smearing function is given by
where r ¼ ðr ðxÞ ; r ðyÞ ; r ðzÞ Þ and is dependant on the regulator parameter Λ (also commonly expressed in terms of a length scale r 0 , related as Λ ¼ ffiffi ffi 2 p =r 0 ). In a finite volume, periodicity can be imposed by replacing g Λ ðrÞ by a regulator which is periodic in each of the spatial directions,
in which the sums run over all periodic copies of the finite volume.
While in many applications of the SVM to nuclear systems the angular momentum structure of the wave function is tied to the spatial structure due to orbital motion, in a cubic box, orbital angular momentum is not a well-defined quantum number. As will be discussed below, one approach in this context is to build wave functions with particular transformation properties under the cubic group, again coupling spatial and spin degrees of freedom. However, since there is a finite number of irreducible representations of the cubic group, a simpler approach is to consider a factorization of the spatial and spin-isospin wave functions.
In this work, a trial wave function, Ψ
where the superscript ðNÞ denotes the total number of terms in the wave function and the dependence of the spatial wave function on h is suppressed. The coordinate x ¼ ðr 1 ; …; r n Þ collects the spatial coordinates of the n particles with x j ¼ r j . The spatial wave functions Ψ sym L are detailed in Sec. III B 1; the c j , j ∈ f1; …; Ng, are coefficients, and the A j , B j , and d j are the parameters of the jth spatial wave function included in the sum. The spin-flavor wave function jχ h i is a vectorfoot_4 in spin-flavor space for the given nucleus h; the particular spin-flavor wave functions that are used in this work are given in Sec. III B 3.
To account for the anisotropy in the spatial wave function 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 in Ref. [51]. States are constructed to be antisymmetric under interchange of the spin-flavor degrees of freedom of pairs of nucleons and thus must have symmetric spatial wave functions under particle interchange.
The basic Gaussian structure underlying these wave functions is
for all integer 3n-tuples,foot_5 n ∈ Z 3n . In order to symmetrize the wave function under particle exchange, the rows and columns of each A ðαÞ and B ðαÞ and rows of d ðαÞ are interchanged under all n! possible permutations, P, of particles. That is,
where A P is the permuted form of A and similarly for B P and d P . As discussed in Appendix B, the shifted Gaussian basis is able to describe scattering states at finite volume as well as compact bound states.
A periodic spatial volume with identical extent in each direction has an underlying cubic symmetry and is invariant under action of elements of the cubic group, H 3 . By imposing cubic symmetry, wave functions that transform in particular representations of H 3 can be constructed, potentially allowing more efficient exploration of the space of correlated shifted Gaussians. Each term in the variational wave function can be constructed to respect the given transformation properties rather than relying on stochastic sampling of a sum of terms to discover the symmetry approximately. The H 3 -covariant wave function transforming in the representation R of H 3 is given by
where p indexes the permutations of the Cartesian directions, c ðRÞ p are constants defining the representation, and A p is the appropriately block-permuted form of the matrix A and similarly for B p and d p . For the ground states that are considered here, the A 1 (trivial) representation of H 3 is assumed, for which c
The utility of using Eq. ( 32), instead of Eq. ( 31), has been investigated. Overall, it is found that N-term wave functions constructed from terms of the form Ψ sym;R L are about a factor of 5 better approximations than N-term wave functions constructed from terms from Eq. (31), measured in terms of the number of wave function terms required to achieve convergence within a given tolerance. However, the cost of evaluation of the matrix elements needed in the SVM is a factor of 6 slower using Eq. ( 32) than using Eq. (31). In the primary studies of this work, trial wave functions are thus constructed using the simpler ansatz in Eq. (31).
The simplest spin-flavor wave functions for the small nuclei that are considered in this work are straightforward to construct explicitly. In particular, the necessary states are defined as
where p ↑ð↓Þ and n ↑ð↓Þ denote proton and neutron states of the given spin.
The trial wave function Ψ ðNÞ of Eq. ( 26) is constructed so as to minimize the bound which it provides on the groundstate energy:
Because of the Gaussian structure of the trial wave function, the various contributions to the Hamiltonian matrix element, which are 3n-dimensional integrals, can be evaluated analytically [51], as shown in Appendix C.
In the current application of the SVM, Ψ ðNÞ is built up from 1 to N terms via an iterative procedure as follows:
(1) Given an M-term wave function (where M < N) defined by matrices A j , B j , d j for j ∈ f1; …; Mg, N proposal proposed candidates for the (M þ 1)th term are constructed by randomly choosing matrices A Mþ1 , B Mþ1 and d Mþ1 . For simplicity of notation, the spatial wave function of the jth term is denoted as Ψ j ðxÞ ≡ Ψ sym L ðA j ; B j ; d j ; xÞ. (2) For each candidate term Ψ Mþ1 ðxÞ, the normalization integrals
and Hamiltonian matrix elements
are computed for fi; jg ∈ f1; …; M þ 1g and used to define matrices N ðMþ1Þ and H ðMþ1Þ , respectively. Note that only the additional (M þ 1)th row and column must be computed, given that the N ðMÞ and H ðMÞ matrices were computed in the previous iteration. (3) The generalized eigenvalue problem
is solved for the eigenvalues λ
Mþ1 , and coefficient vectors c , is selected from the set of N proposal candidates 8 and added to the iteratively constructed trial wave function to define Ψ ðMþ1Þ ðxÞ. Naturally, the optimization at each step of this iterative procedure depends on the Hamiltonian H and hence on the LECs that define it; to enable optimization of the trial wave function across a broad range of LECs, in the numerical study undertaken here, the values of the LECs are varied for each step of optimization, cycling through N couplings choices that span the relevant parameter spaces. After initializing the procedure with the first trial wave function (the M ¼ 1 term), for which the matrices N ð1Þ and H ð1Þ are single numbers and the generalized eigenvalue problem is trivial, additional terms are added iteratively until the wave function has a fixed number of terms, N. N must be taken large enough that the optimization procedure has converged by some definition. Here, N is set by the criterion that repeated optimizations based on different random seeds achieve the same minimum energies within some tolerance and that adding some fixed number of additional terms to the trial wave functions does not alter the minimum found, within the same tolerance. Details of the optimization procedure for the cases considered in this study are provided in Sec. IV.
In this work, the finite-volume SVM is used to match the LECs of FVEFT = π to the LQCD results for two-and threenucleon systems which were obtained in Ref. [8], where nuclear states with SUð3Þ f -symmetric quark masses corresponding to m π ¼ 806 MeV were studied in three volumes of spatial extents L ∈ f3.4; 4.5; 6.7g fm. The binding energies extracted in that work, defined as ΔE h ¼ E h -AE p where A is the atomic number of the state h, are 8 Note that each of the N proposal lowest eigenvalues λ ðMþ1Þ 0 is smaller than the lowest eigenvalue from the previous iteration λðMÞ 0 .
tabulated in Table I. This matching procedure, used in this work as a first step in the study of the matrix elements of various currents in these states in this framework, closely mirrors, and reproduces, the analysis of Ref. [21].
For each of the three finite volumes for which Ref. [8] provides LQCD data, variational wave functions were optimized first for the pp and d two-body systems. As described in Sec. III C, the undetermined coefficient of the two-body potential in the Hamiltonian (which corresponds 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 different choices throughout the optimization procedure, which are chosen to be evenly spaced corresponding approximately to the plot range of Fig. 2; both two-body systems are thus optimized simultaneously. Taking N proposal ¼ 30, it is found for all two-body optimizations undertaken in this work-at each finite volume and for each of three choices of the regulator Λ ¼ ffiffi ffi 2 p =r 0 corresponding to r 0 ∈ f0.2; 0.3; 0.4g fm-that after 100 terms have been added to the wave function the last 20 terms are within 1% (typically, within 0.1%) of the final value for each of the values of the coupling that are used in optimization and also that optimizations starting with different random seeds agree within that same tolerance. An example of this convergence is shown in Fig. 1.
The ground-state d and pp binding energies resulting from the optimized variational wave function are shown as a function of the relevant undetermined LEC in Fig. 2, for all three finite volumes and for all choices of the regulator parameter Λ which are studied. Comparing with the LQCD results for the binding energies in each lattice volume, it is clear that for each value of Λ the results in all volumes are consistent with a single consistent value of the relevant coupling, indicating that there is no need to introduce higher-order terms in EFT π . The best-fit values of the corresponding couplings, which depend on the regulator scale, are determined through a combined fit to the three volumes and are presented in Table II and Fig. 3. Note that the EFT π interaction proportional to C 1 is suppressed by 1=N 2 c relative to the Wigner-symmetric interaction with coefficient C 0 in the large-N c limit; this hierarchy is born out in the fitted values of the couplings. to the ground-state energy of the diproton system as additional terms are added to the variational wave function. The three colors correspond to wave functions optimized for the L ¼ 3.4 (blue), 4.5 (orange), and 6.7 fm (green) volumes, with the LEC C T set to its optimized value after matching to the LQCD results, as discussed in the main text. The results are shown for one example of the regulator parameter corresponding to r 0 ¼ 0.2 fm. For each volume, the results of optimization procedures starting with two random seeds are shown. =r 0 corresponding to r 0 ∈ f0.2; 0.3; 0.4g fm (solid, dashed, and dotted), and the three colors correspond to the three different volumes: L ¼ 3.4 (blue), 4.5 (orange), and 6.7 fm (green). The horizontal bands show the values of the binding energies for each volume from the LQCD calculations of Ref. [8].
Having determined the two-body couplings, wave functions optimized in infinite volume in exactly the same way are used to determine the binding energies in the infinitevolume limit. Figure 4 shows the volume dependence of the binding energies of the deuteron and diproton systems (in order to show a curve, binding energies computed with wave functions optimized at additional intermediate volumes are also shown). Extrapolations are shown for r 0 ¼ 0.2 fm, but the extrapolations for other values are indistinguishable. Although the values of the LECs depend on the value of the regulator r 0 , the resulting finite-and infinite-volume energies are regulator independent. Table III lists the extrapolated binding energies and compares them to Refs. [21,25], with which they are in close agreement.
Having determined the two-nucleon couplings, an analogous procedure is repeated for the triton to determine the coefficient D 0 of the leading-order three-nucleon coupling in Eq. (7). The stochastic optimization of the three-body wave function is performed with the two-body LECs fixed TABLE II. The two-and three-body EFT = π LECs determined from matching the SVM calculations to the LQCD energy shifts for each value of the cutoff parameter r 0 . For the two-body case, C S and C T are determined in terms of C 0;1 through Eq. ( 6). For D 0 , three-nucleon optimizations were only performed for r 0 ¼ 0.2 fm. =r 0 studied here. The asterisks show the results of an analogous analysis undertaken in Ref. [21], with different values of the regulator r 0 , as indicated on the figure . FIG. 4. The volume dependence of the deuteron and diproton binding energies compared with the LQCD data which was used to determine the relevant LECs. The infinite-volume extrapolations of the binding energies, computed as described in the text, are shown in the rightmost subpanel. TABLE III. The extrapolated infinite-volume binding energies determined in the SVM approach for the two-and three-body systems. Also shown are the binding energies determined in the original LQCD study [25] and in Ref. [21] to their central values determined as discussed in the previous subsection, for a single value of the Gaussian regulator parameter r 0 ¼ 0.2 fm. As for the two-body case, N couplings ¼ 10 values of the three-body LEC D 0 are cycled through in the wave function construction procedure, spanning the relevant coupling range (corresponding approximately to the range of the horizontal axis in Fig. 5). The same convergence criteria as in the two-body case are satisfied after wave functions with N ¼ 250 terms have been constructed. Figure 5 shows the dependence of the triton binding energy on the three-body coupling, D 0 , for the optimal values of the two-body couplings; the shaded bands around the curves show the result of varying the twobody LECs within their uncertainties for wave functions optimized in each of the three volumes for which there are LQCD data. Figure 6 and Table III show the infinitevolume extrapolation of the triton binding energy. As for the two-body systems, the extrapolations and couplings reported here are in close agreement with those of Ref. [21].
Having determined finite-volume ground-state (or in principle excited state) wave functions in the SVM, those wave functions can be used to evaluate finite-volume matrix elements of operators in FVEFT π . The transition matrix element between an initial state Ψ i and final state Ψ f is given by where O is a generic local EFT operator and bra-ket notation is used for concision. 9Here, matrix elements of axial, electromagnetic, scalar, and tensor currents and the unpolarized twist-2 operators are studied. In each case, the relevant EFT π currents of Sec. II are translated into operators acting on the n-body states. As with the nucleon-nucleon strong interactions, two-body contributions to the various currents are regulated using the Gaussian approach as in Eq. ( 24) and rendered periodic using Eq. (25). Specifically, each two-body current is implemented as ½ðN † ðr i ÞΣNðr i ÞÞðNðr j Þ † Σ 0 Nðr j ÞÞ þ H:c:g Λ ðr ij ; LÞ; ð39Þ where Σ ð0Þ denotes a spin-isospin structure and g Λ ðr; LÞ, defined as in Eq. ( 24), implements a periodic regulated form of the δ-function implied in local two-body EFT currents. Since the matrix elements that are considered are for zero momentum transfer, the current is integrated over the positions r i;j . For each current, X ∈ fA i;a ; J EM i ; S a ; T ij;a ; O ðnÞ a g for a ∈ f0; 1; 2; 3g, the zero-momentumprojected, regulated form is labeled as X .
The evaluation of the relevant matrix elements factorizes into a spin-isospin calculation that is specific to each type of operator and a calculation of the matrix element of the spatial wave function. Since all currents that are considered enter with the spatial dependence determined by the Gaussian regulator function, these latter spatial matrix elements have a common form and are given for diagonal matrix elements in state h by
where R n ¼ fr 1 ; …; r n g indicates dependence on the coordinates of each of the n particles. For transition matrix elements between states a and b, the corresponding expression is
The relevant one-and two-body LECs of currents in FVEFT = π can be tuned such that the EFT π matrix elements determined in this way reproduce the matrix elements (or their ratios to the proton matrix element) determined in LQCD in a particular lattice volume or set of volumes. In what follows, the LQCD matrix elements in atomic number A ¼ f2; 3g systems calculated at m π ¼ 806 MeV on the L ¼ 4.5 fm ensemble, discussed above, are used. The EFT π counterterms determined in this way are specified in the Gaussian-regulated scheme and should not be compared with the corresponding counterterms determined in dimensional regularization. Indeed, for the current purposes, the extraction of counterterms can be viewed simply as an intermediate step in extracting matrix elements at infinite volume. The infinite-volume-extrapolated matrix elements can be matched to EFT π regulated in more common schemes, such as dimensional regularization or the power-divergent subtraction scheme [52], to determine the two-body LECs for comparison to other extractions. A. Axial matrix elements: Proton-proton fusion, tritium β-decay, and isoscalar charges
In order to evaluate the isovector axial current matrix elements, the EFT current in Eq. (A1) is used. With the proton axial matrix element determining g A (up to exponentially small volume effects), ratios of matrix elements of the relevant current in both two-and three-body states to that of the proton can be used to determine the two-body coupling ratio L1;
where the spin-flavor structure of the states used to arrive at these expressions is given in Sec. III. Figure 7(a) shows the constraints that the LQCD calculations [53] of the finite-volume matrix elements in the two channels place on the coupling combination L1;A . The consistency between the constraints in the two channels suggests that higher-order terms in the axial current (two-body operators with derivative insertions [43] or three-body operators) are suppressed as their power counting would suggest, for this choice of regulator scheme and scale. Were this to persist for physical quark masses, it would provide support for approaches to pp-fusion cross section calculations that use tritium β-decay to constrain the relevant twobody LEC.
The values of L1;A determined from each channel are scheme-dependent quantities but can be combined with infinite-volume SVM wave functions to determine the infinite-volume matrix elements. Figure 7(b) shows the infinite-volume extrapolation for both channels, and the extrapolated values are given in Table IV below.
Analogous analysis of the isoscalar axial matrix elements in the deuteron and 3 He states allows for the determination of the two-body counterterm in Eq. (11). Ratios of the isoscalar axial current matrix element in the deuteron and 3 He states to that in the proton state can be expressed as
Figure 8 shows the constraints on L2;A obtained by matching to the LQCD calculation of Ref. [54] and the corresponding infinite-volume extrapolation of the LQCD matrix elements. The extrapolated values are also reported in Table IV. A mild tension is found between the values of L2;A extracted from each matrix element, indicating the potential need for higher-order terms in the EFT description.
The magnetic moments of light nuclei have been extracted from the linear response of LQCD calculations to a constant background magnetic field oriented in the z-direction [55,57]. In EFT, these quantities are determined by the couplings in Eq. ( 13). The differences between the magnetic moments of the deuteron, 3 H, and 3 He states and the relevant naive shell-model predictions in terms of proton and neutron magnetic moments can be expressed as
where μh is the magnetic moment of hadron h in natural nuclear magnetons e=2M N defined using the nucleon mass at the quark masses of the lattice calculations, M N ¼ 1.634ð18Þ GeV [25].
FIG. 9. The dependence of the various magnetic moment differences on the appropriate combinations of the two-body counterterms L 1;2 is shown in the upper row and the lower left panel. The lower right panel shows the combined constraints implied by agreement with the LQCD results of Ref. [55]. The details of the curves and points in the figure are as in Fig. 7.
FIG. 10. The infinite-volume extrapolations of the magnetic moment differences for d, 3 H, and 3 He. In Fig. 9, the LQCD calculations of δμ h for h ∈ fd; 3 H; 3 Heg in Ref. [55] are used to constrain the EFT couplings L 1;2 . 10 Since the magnetic moment differences are dependent of various combinations of the couplings, the constraints take the form of bands in the L 1 -L 2 plane as shown in the figure. All three constraints are seen to be consistent for a range of values of the couplings. This determination of the LECs allows for extrapolation of the magnetic 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 and isovector scalar currents in hadron h to those in the proton can be expressed as
and
where A h denotes the atomic number of the nucleus, T h 3 is its third component of isospin, and
The quantity CV ¼ CV =g S;3 is constrained by comparing Eq. ( 55) to LQCD calculations of ratios of isovector scalar current matrix elements in different nuclei from Ref. [54]. The results of this comparison are shown in Fig. 11. Determinations of CV from both the pp and 3 H systems are consistent, with 3 H providing a considerably more stringent constraint. The extracted values of the LEC are also used to extrapolate the LQCD matrix elements to infinite volume, as shown in the figure and presented in Table IV.
Similarly, Fig. 12 compares LQCD calculations of the ratio of the scalar isoscalar current matrix element in nuclei to that in the proton to the expectations of Eq. ( 49) for the deuteron, diproton, and 3 H. The three states provide sufficient information to constrain the LEC ratios C0;1 ¼ C0;1 =g S;0 , and the constrained values are used to extrapolate the LQCD matrix elements to infinite volume as shown in Fig. 13 and presented in Table IV. As noted in Sec. II C, the couplings C0;1 that occur in the isoscalar scalar current are related to the Lagrangian counterterms C 0;1 in the limit of massless quarks.
For the tensor current, isoscalar matrix element ratios to those in the nucleon are given by
for h ∈ fd; 3 Hg, where S h 3 is the third component of spin and FIG. 13. The infinite-volume extrapolations of the isoscalar scalar current matrix element ratios for d, pp, and 3 H. 10 Note that in Ref. [55] the magnetic background field does not couple to sea quarks, so only isovector quantities are calculated completely; the error from this quenching of the magnetic field is ignored here. In principle, the isovector np → dγ M1 transition can also be used to constrain L 1 , but it is not used in this work.
where L2;T ¼ L 2;T =g T;0 . For the isovector case, the corresponding ratio in 3 H is 52) and ( 54) to the respective LQCD calculations [54]. In the isoscalar case, consistency is seen in the values of L2;T that arise from comparison to either d or 3 H matrix element ratios from LQCD. As for the matrix elements of other operators considered above, the constrained LECs enable infinite-volume extrapolations of the matrix elements, which are presented in Table IV.
The first moment of the isovector unpolarized parton distribution has been computed in LQCD for nuclei with A ∈ f2; 3g in Ref. [56] and corresponds to the difference in the longitudinal momentum fractions carried by up and down quarks. In the finite-volume SVM, the matrix elements of the operators in Eq. ( 20) contain no spin structure, and so the finite-volume matching and extrapolation are analogous to those on the isovector scalar current above. In particular, ratios of the isovector matrix element in hadron h to the proton matrix elements are given by
for h ∈ fpp; 3 Hg. LQCD data in a calculation with L ¼ 4.5 fm constrain the single-nucleon isovector momentum fraction and twonucleon counterterm α 1;3 , as shown in Fig. 16. The extrapolated infinite-volume matrix elements are reported in Table IV.
Finite-volume pionless effective field theory implemented through the stochastic variational method has been used to extrapolate EFT wave functions and matrix elements for A ∈ f2; 3g nuclei, matched to LQCD calculations in a finite lattice volume, to infinite volume. This numerical approach can effectively describe bound-state systems, can cleanly reproduce scattering states, and furthermore exhibits volume scaling that is consistent with the predictions of the Lüscher approach to high accuracy. To some degree, the method circumvents the complexities of analytic approaches generalizing that of Lüscher for twobody systems to larger number of particles. However, as the atomic number of the system increases, the finite-volume SVM scales relatively poorly, and it does not seem practical to extend much beyond A ¼ 4 systems or more than threebody interactions. Given that the use of the finite-volume aspect of the method is tailored to match LQCD calculations, for which increasing A is also costly, this is perhaps not a significant limitation; the finite-volume SVM can be used to determine EFT counterterms, which can then be used in the infinite-volume SVM (or other many-body methods) to perform calculations 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 and the isovector longitudinal momentum fraction, it is found that for the large quark masses used in the LQCD calculations the lattice volume of L ¼ 4.5 fm as used for the calculations is large enough that there are essentially no finitevolume corrections. At lighter quark masses, however, one might anticipate that larger lattice volumes will be required to achieve the same behavior. For almost all of the matrix elements investigated in this work, it is also notable that, although the constraints from LQCD calculations of threebody systems are typically tighter, the EFT π LECs determined from the LQCD calculations of two-and three-body systems are consistent to within one standard deviation [with the notable exception of the isoscalar axial LEC L 2;A in Eq. (11), for which there is a slight tension]. This indicates that higher-order terms in the relevant currents are suppressed in the exponential regulator scheme used in the SVM at these values of the quark masses. In some cases, this suppression has particular consequence; for example, for the isovector axial matrix elements, the LEC L 1;A determined from measurements of tritium β-decay is used in calculations of the pp-fusion cross section [58]. With the LECs determined from LQCD calculations, predictions can be made for other quantities for which there are no LQCD results; in Table IV, as an example, the np ← d tensor transition matrix element is predicted from the LEC determined from the triton matrix element of the same current. Since the finite-volume SVM provides representations of low-lying excited states as well as the ground states that have been the focus of this work, it can also be used to match matrix elements of second-order current insertions such as in double-β decay. In such processes, a sum over excited nuclear states occurs for times between those of the two currents, with the matrix elements of interest being
In the finite-volume EFT context, this corresponds to inclusion of the discrete states in principle up to energyscale of the EFT cutoff. In order to accurately represent these contributions, care must be taken that these states are equivalently well optimized. This will require significant numerical effort but is likely to be feasible. Ultimately, the finite-volume SVM appears to be a powerful tool to capitalize on LQCD calculations of systems with small A, which are approaching a novel era of systematic control.
and sign changes for N ¼ 1, fp 1 ¼ ð1;1;0Þ;p 2 ¼ ð0;0;0Þg or fp 1 ¼ ð1; 0; 0Þ; p 2 ¼ ð0; 1; 0Þg and permutations and sign changes for N ¼ 2, and so on. As can be seen from the figure, the SVM is able to cleanly reproduce the low-energy part of the spectrum, including its degeneracies; with further numerical effort, this can be extended further. Similarly, the free three-particle low-energy spectrum is well reproduced. The large L asymptotic behavior of the ground state is given by an expansion of the two-particle quantization condition derived by Lüscher [26], namely
where a is the scattering length and the geometric coefficients are c 1 ¼ -2.837297 and c 2 ¼ þ6.375183. Figure 17(b) shows the two-body ground-state energy extracted from the SVM with a small repulsive coupling [i.e., the energy in the isoscalar channel with C S ¼ 31 MeV fm 3 and C T ¼ 0 in Eq. ( 3)] as a function of volume. Also shown are fits to this data using Eq. (B1) with the scattering length as a free parameter, as well as fits treating both the scattering length and the geometric constant c 1 as free parameters. The latter fit returns a value of c fit 1 ¼ -2.63ð19Þ that is consistent with the actual value, showing that the SVM is correctly able to reproduce the expected finite-volume asymptotic behavior. Similar agreement with the expect asymptotic behavior is found for small attractive interactions.
In this section, explicit formulas for the wave function integrals in Eqs. (35) and (36) are provided. In the expressions below, Ψ sym L ðA i ; B i ; d i ; xÞ is a symmetric n-body Gaussian wave function term defined in Eq. (31). The normalization integral of Eq. ( 35), corresponding to a cross-term between such n-body Gaussian wave function terms labeled by the subscripts i and j, respectively, can be expressed as
and where superscripts ðαÞ and subscripts P on the wave function parameters A, B, and d extract the components of the parameters corresponding to the α direction and permute the parameters for each of the n particles by the permutation P, respectively. In the above equations, b ðαÞ is a length-n vector; summing over all n-vectors corresponds to enforcing periodicity. In practice, the infinite sum is truncated to vectors with a maximum norm b, and in the numerical calculations in this work, the cut is initially taken to be b ¼ 3 for each term and is iteratively increased until the fractional change in the total summed expression from adding an additional term is less than 10 -10 .
The spatial integrals involved in the Hamiltonian matrix elements are separated into the kinetic and two-and threebody potential terms as
where the spatial integrals for each of K, V 2 , and V 3 can be performed independently. The integral for the two-body potential term, again for n-body Gaussian wave function terms labeled by i and j, is
where ½V a denotes the ath component of a vector V and ½M -1 ab denotes the ða; bÞ component of the matrix M -1 . In this expression, q ðαÞ is an integer; combined with the sum over b ðαÞ , summing over all values of q ðαÞ corresponds to enforcing periodicity. In practice, the infinite sum is truncated to integers with absolute value less than q ¼ 40; the fractional change in the total summed expression from adding an additional term is less than 10 -10 .
The relevant integral for the three-body potential term, for n-body Gaussian wave function terms labeled by i and j, can be expressed as
As in Eqs. (C1) and (C6), the sums over b ðαÞ , q ðαÞ , and t ðαÞ together enforce periodicity. In numerical evaluations of Eq. (C8), the same cut procedure for fixing b and q is used as for the evaluations of the previous expressions. Finally, the integral for the kinetic operator, for n-body Gaussian wave function terms labeled by i and j, is
The Lüscher approach has been extended to n-particle systems under certain assumptions in Refs.[31][32][33][34][35][36][37][38][39].
Note that the normalization of L 1;A used here is the same as in Ref.[43] but differs from Ref.[44], which uses a different projector definition.
Note that Ref.[44] uses a different projector definition than here and correspondingly a different normalization of L 2;A .
The notation of Ref.[45] is used.
Note that χ is common to all terms in Eq.(26) in the current implementation. In other approaches for larger systems than will be considered here, χ is also part of the stochastic sampling and would be indexed by j[49].
In practice, these sums are truncated as discussed in Appendix C.
Matrix elements of nonlocal products of operators (such as those relevant for double-β decay) can also be approached in the finite-volume SVM as discussed in Sec. VI.