Theoretical Foundation of the Nuclear Force in QCD
and Its Applications to Central and Tensor Forces
in Quenched Lattice QCD Simulations
1 Introduction
The origin of the nuclear force is one of the major unsolved problems in particle and nuclear physics even after the establishment of the quantum chromodynamics (QCD). Although the nuclear force is still not wellunderstood theoretically, a large number of protonproton and neutronproton scattering data as well as deuteron properties have been accumulated and summarized e.g. in the Nijmegen database [1]. To describe the elastic nucleonnucleon () scattering at lowenergies below the pion production threshold together with the deuteron properties, the notion of the potential (either in the coordinate space or in the momentum space) turns out to be very useful [2]: it can be determined phenomenologically to reproduce the scattering phase shifts and bound state properties either through the Schrödinger equation for the wave function or through the LippmannSchwinger equation for the matrix. Once the potential is determined, it can be used to study systems with more than 2 nucleons by using various manybody techniques.
Phenomenological potentials which can fit the data precisely (e.g. more than 2000 data points with ) at MeV are called the highprecision potentials: They include the potentials such as CDBonn [3], Argonne [4], and Nijm I, Nijm II and Reid93 [5]. Also systematic low energy construction of the nuclear force on the basis of the chiral perturbation theory is being developed [6, 7].
The phenomenological potentials in the coordinate space are known to reflect some characteristic features of the interaction at different length scales [2]:

The long range part of the nuclear force (the relative distance fm) is dominated by the onepion exchange introduced by Yukawa [8]. Because of the pion’s NambuGoldstone character, it couples to the spinisospin density of the nucleon and hence leads to a strong spinisospin dependent force, namely the tensor force.

The medium range part ( fm) receives significant contributions from the exchange of twopions () and heavy mesons (, , and ). In particular, the spinisospin independent attraction of about 50 – 100 MeV in this region plays an essential role for the binding of atomic nuclei.

The short range part ( fm) is best described by a strong repulsive core as originally introduced by Jastrow [9]. Such a short range repulsion is important for the stability of atomic nuclei against collapse, for determining the maximum mass of neutron stars, and for igniting the Type II supernova explosions [10].

There is also a strong attractive spinorbit force in the isospin 1 channel at medium and short distances. This leads to the neutron pairing in neutron matter and hence the neutron superfluidity inside neutron stars [10].
A repulsive core surrounded by an attractive well is in fact a common feature of the “effective” potential between composite particles. The LenardJones potential between neutral atoms or molecules is a wellknown example in atomic physics. The potential between He nuclei is a typical example in nuclear physics. The origin of the repulsive cores in these examples is known to be the Pauli exclusion among electrons or among nucleons. The same idea, however, is not directly applicable to the potential, because the quark has not only spin and flavor but also color which allows six quarks to occupy the same state without violating the Pauli principle. To account for the repulsive core of the force, therefore, various ideas have been proposed as summarized in Ref. 11: exchange of the neutral meson [12], exchange of nonlinear pion field [13], and a combination of the Pauli principle with the onegluonexchange between quarks [14]. Despite all these efforts, convincing account of the nuclear force has not yet been obtained.
In this situation, it is highly desirable to study the interactions from the first principle lattice QCD simulations. A theoretical framework suitable for such purpose was first proposed by Lüscher [15]: For two hadrons in a finite box with the size in periodic boundary conditions, an exact relation between the energy spectra in the box and the elastic scattering phase shift at these energies was derived: If the range of the hadron interaction is sufficiently smaller than the size of the box , the behavior of the twoparticle BetheSalpeter (BS) wave function in the interval under the periodic boundary conditions has sufficient information to relate the phase shift and the twoparticle spectrum.
Lüscher’s method bypasses the difficulty to treat the realtime scattering process on the Euclidean lattice.^{*)}^{*)}*) There are several studies of the interactions on the lattice, which do not rely on Lüscher’s method. One uses the BornOppenheimer picture, i.e., if one of the three quarks inside the baryon is infinitely heavy, one may define the potential between baryons a la BornOppenheimer [16]. This is, however, not applicable to the nucleons with light quarks. The other employs the strong coupling limit, which has been proposed quite recently [17]. Furthermore, it utilizes the finiteness of the lattice box effectively to extract the information of the onshell scattering matrix and the phase shift. This approach has been applied to extract the scattering lengths in the quenched QCD simulations [18] and in the (2+1)flavor QCD simulations with the mixed action [19].
Recently, the present authors proposed a closely related but an alternative approach to the interactions from lattice QCD [20, 21]. The starting point is the same BS wave function as discussed in Ref. 15. Instead of looking at the wave function outside the range of the interaction, we consider the internal region and define the energyindependent nonlocal potential from so that it obeys the Schrödinger type equation in a finite box. Since for strong interaction is localized in its spatial coordinates due to confinement of quarks and gluons, the potential receives finite volume effect only weakly in a large box. Therefore, once is determined and is appropriately extrapolated to , one may simply use the Schrödinger equation in the infinite space to calculate the scattering phase shifts and bound state spectra to compare with experimental data. Further advantage of utilizing the potential is that it would be a smooth function of the quark masses so that it is relatively easy to handle on the lattice. This is in sharp contrast to the scattering length which shows a singular behavior around the quark mass corresponding to the formation of the bound state.^{†)}^{†)}†) Similar situation is wellstudied in connection with the BECBCS crossover in cold fermionic atoms [22], where the external magnetic field plays a role of the quark mass in QCD. For seminal suggestion on the rapid quarkmass dependence of the scattering length, see Ref. 23.
Since we consider the nonasymptotic region () of the wave function, the resultant potential and the matrix are offshell. Therefore, they depend on the nucleon interpolating operator adopted to define the BS wave function. This is in a sense an advantage, since one can establish a onetoone correspondence between the nucleon interpolating operator and the potential in QCD, which is not attainable in phenomenological potentials. It also implies that the potential on the lattice and the phenomenological potentials are equivalent only in the sense that they give the same observables, so that the comparison of their spatial structures should be made only qualitatively.
The purpose of this paper is twofold: First, we will present a theoretical foundation of our method to extract the potentials from lattice QCD. Then, we will give a full account of the application of the method to the quenched lattice QCD simulations. Once our method in lattice QCD is proved to work in the system, it will have various applications not only to nuclear manybody problems but also to hyperonnucleon, hyperonhyperon and threenucleon interactions which have much less experimental information than the systems. A first attempt to the hyperonnucleon potential has been already reported in Ref. 24, and more on hyperons will appear in the future publications.
This paper is organized as follows. In §2, we illustrate the derivation of the twobody and manybody potentials from the wave function in quantum mechanics. In §3, the idea in the previous section is generalized to the interaction of composite particles in field theory. In §4, we classify the general structure of the potential in the velocity expansion and show the procedure to determine each term. In §5, the method to determine the potential from the lattice QCD data is discussed in detail for the effective central potential at low energy. We also discuss the method to extract the tensor potential in our approach. In §6, potentials obtained from the quenched lattice QCD simulations are presented. Section 7 is devoted to summary and concluding remarks. In Appendix A, a fieldtheoretical derivation of the asymptotic BS wave function at large distance is presented. In Appendix B, the way to make general decomposition of the potential (the OkuboMarshak decomposition[25]) is reviewed. In Appendix C, matrix elements of the general potential are presented. In Appendix D, heatkernel representation of the Green’s function is presented.
2 Nonlocal potential in quantum mechanics
2.1 Twobody force
To show the basic concept of the nonlocal potential in a finite box with the size , we start with a nonrelativistic twobody problem described by the stationary Schrödinger equation:
(2.1) 
where is the relative coordinate of the two spinless and nonrelativistic particles, and is related to the discrete energy eigenvalues ( with being the reduced mass. The wave function obeys the periodic boundary condition. The nonlocal potential ^{‡)}^{‡)}‡)Here we use the standard term “nonlocal” in the sense that cannot be written as . is assumed to be energyindependent and Hermitian, , so that the discrete energy eigenvalues are real and corresponding eigenfunctions can be made orthonormal. For the scattering states (bound states) in the infinite volume, we have (). On the other hand, negative in the finite volume does not necessarily imply the existence of the bound state at .
We consider the potential whose spatial extension is sufficiently small in the sense that is exponentially suppressed for with being smaller than . We define the “inner region” by . Then, the wave function in the “outer region” satisfies the Helmholtz equation, , with the periodic boundary condition.
Let us consider the following inverse problem: Suppose we have no information about except that it is smooth and short ranged, while we know linearly independent wave functions and associated energy in a finite box for .^{§)}^{§)}§)This is more luxurious situation than the usual inverse scattering problem where only the scattering phase shifts in the outer region are available. Now, we introduce the following function:
(2.2) 
where is the nonrelativistic kinetic energy operator satisfying . Since removes the noninteracting part of the wave function, is nonvanishing only in the inner region irrespective of the sign of .
By taking into account the fact that may not be orthonormal, we introduce the norm kernel , so that the projection operator to the space spanned by the wave functions with reads . Then, an energyindependent and nonlocal potential can be defined as
(2.3) 
which leads to the Schrödinger equation Eq. (2.1) for . If we apply a unitary transformation to the wave function, , the nonlocal potential is modified as . Such unitary transformation does not affect the observables, while it changes the spatial structure of the wave function and the nonlocal potential.
If are all real and , the potential becomes a hermitian operator in the subspace . Otherwise, the hermiticity is not obvious and should be checked case by case. In field theory discussed later, corresponds to the equaltime BetheSalpeter amplitude in a finite box and corresponds to the threshold energy of inelastic channels.
In practice, the potential defined in Eq. (2.1) has limited use, because the number of states satisfying the condition is not generally large for lattice QCD in a finite box. This problem can be evaded when we focus on the lowenergy scattering with sufficiently smaller than the intrinsic scale of the system or the scale of the nonlocality of the potential. In such a case, the velocity expansion of in terms of its nonlocality is useful [26]: For example, a spinindependent potential with hermiticity, rotational invariance, parity symmetry, and timereversal invariance can be expanded as
(2.4)  
(2.5) 
where and with . Each coefficient of the expansion is the local potential and can be determined successively by the wave functions at low energies: For example, if we have five wave functions corresponding to , we obtain
(2.6) 
Pretending that and are independent of each other, Eq. (2.1) for can be solved algebraically to obtain and . Hermiticity of the potential can be checked by the consistency among the local potentials thus determined. Stability of the potentials against the number of wave functions introduced can be also checked.
An advantage of defining the potential from the wave functions in the “inner region” is that the effect of the periodic boundary condition is exponentially suppressed for finite range interactions: Then one can first make appropriate extrapolation of or to , and then solve the Schrödinger equation using the extrapolated potential to calculate the observables such as the phase shifts and binding energies in the infinite volume.^{¶)}^{¶)}¶)Strictly speaking, the local potentials with higher derivatives must be treated as perturbation to keep the Schrödinger equation as a second order differential equation. This is in contrast to the approach by Lüscher [15] in which the wave functions in the “outer region” suffering from the boundary conditions is ingeniously utilized to probe the scattering observables. Apparently, the two approaches are the opposite sides of a same coin.
2.2 Manybody forces
For the interactions among composite particles, there are in principle manybody forces which take place in the system composed of more than two particles. The wellknown example in nuclear physics is the FujitaMiyazawa type threebody force acting among three nucleons [27, 28]. It is phenomenologically important for the extra binding of light nuclei [29] and for the extra repulsion in high density matter [30] and in elastic nucleusnucleus scatterings [31].
The method to define the twobody potential from the relative wave function discussed above can be generalized to the manybody forces. Let us illustrate the procedure by considering the threebody system of spinless and distinguishable particles with equal mass . We consider the local potentials for both twobody and threebody forces just for simplicity. In the rest frame of the threebody system, we have
(2.7) 
where and are the Jocobi coordinates. and are the kinetic energy operator with and being the reduced masses. is the total energy of the threebody system at rest. Because of the translational invariance, the twobody potential and the threebody potential are the functions of and .
If we know the wave function and the total energy on the lefthand side of Eq. (2.2), the threebody potential can be determined by the following procedure. We first consider the situation, , where , and are vanishingly small because of the assumed shortrange nature of the potentials. Then, can be determined by changing within the range . One can carry out similar procedure to determine and . Alternatively, one may determine from the genuine twobody system.
Once all the twobody potentials are determined, can be extracted from the wave function in the range, and , through the threebody equation Eq. (2.2). It is important to note that the threebody potential is always obtained together with the twobody potential: they are closely tied through the wave function. If one makes the unitary transformation of the wave function, both and are changed simultaneously.
The above procedure can be formally generalized to the nonlocal potentials and to the particle systems with different masses and internal degrees of freedom.
3 Nonlocal potential in field theory for spin 1/2 particles
3.1 BetheSalpeter wave function
In field theory, the best analogue of the twoparticle wave function is the equaltime BetheSalpeter (BS) amplitude, so that we use the term “BS wave function” throughout this paper. Let us consider the following BS wave function for the 6quark state with total energy and the total threemomentum in a finite box ,
(3.1) 
where the relative coordinate is denoted as . The local composite operators for the proton and the neutron are denoted by and with spinor indices and . The QCD vacuum is denoted by , while the state is a QCD eigenstate with baryon number 2 and with the same quantum numbers as the pn system. One should keep in mind that is not a simple superposition of a product state , since there are complicated exchanges of quarks and gluons between the two composite particles. The stationary BS wave function may be regarded as a probability amplitude in to have “neutronlike” threequarks located at point and “protonlike” threequarks located at point .
The spatial extent of the interaction in QCD is short ranged and is exponentially suppressed beyond the distance fm. Therefore, the spatial part of the BS wave function in the “outer region” () satisfies the Helmholtz equation, , up to an exponentially small correction. Here the “asymptotic momentum” is related to the total energy through the relation, . To make a formal resemblance with the nonrelativistic case, we introduce the “effective center of mass energy”, [15]. As shown in Appendix A, using the unitarity of the matrix, we can show that the asymptotic behaviour of the BS wave function at large is identical to that of the scattering wave in the quantum mechanics, with the identification that the phase of the matrix is the scattering phase shift of the BS wave function.
Now, we apply the same logic as the quantum mechanical case in §2.1. The threshold of the pion production is chosen to be . Namely, is a function which has a support only in the inner region as long as stays below the threshold. Thus we can define the shortranged nonlocal potential as
(3.2)  
(3.3)  
(3.4) 
where and . By construction, the solution of Eq. (3.1) with extrapolated to reproduces the correct BS wave function in the asymptotic region, and hence the phase shifts and binding energies of the twonucleon system.
3.2 Interpolating operators
In Eq. (3.1), simplest interpolating operators for the neutron and the proton written in terms of the upquark field and the downquark fields would be
(3.5)  
(3.6) 
where , and the color indices are denoted by , and . The charge conjugation matrix in the spinor space is denoted by .
As shown in Appendix A, local operators such as given in Eqs. (3.2) and (3.2) are most convenient for relating the BS wave function to the fourpoint Green’s function and the scattering observables at . Closely related observation was obtained long time ago by Nishijima, Zimmermann and Hagg who derived the generalized reduction formula for local composite fields [33].
In principle, one may choose any composite operators with the same quantum numbers as the nucleon to define the BS wave function.^{∥)}^{∥)}∥)In practice, however, we had better restrict ourselves to consider only local composite operators for the nucleon, since it is very difficult, although not entirely impossible, to derive the reduction formula for nonlocal composite operators without violating the causality of relativistic theories. Different operators give different BS wave functions and different potentials, although they lead to the same observables such as the phase shifts and binding energies. This is quite analogous to the situation in quantum mechanics that the unitary transformation of the wave function changes the structure of the potential while the observables are not modified. A theoretical advantage of our approach based on lattice QCD is that we can unambiguously trace the onetoone correspondence between the potential and the interpolating operator in QCD. This is in contrast to the phenomenological potentials where connections to QCD operators are not attainable.
4 General form of the potential
In the previous section, we illustrated the procedure to define the potential between the neutron and the proton, which has spinor indices running from 1 to 4. In order to derive the general structure of the potential at low energies, we restrict ourselves to consider only the upper components of these spinor indices in the following sections.
4.1 Symmetry of the two nucleon system
0  1  
0  1  0  1  
odd  even  even  odd  
,  
,  
,  
,  
⋮  ⋮  ⋮  ⋮  ⋮  ⋮  ⋮ 
It is useful to classify the asymptotic twoparticle states by the orbital angular momentum (), the total spin () and the total angular momentum () together with the total isospin . Using the standard notation, , and taking into account constraints due to Pauli principle, we have the wellknown relations given in Table I.
4.2 OkuboMarshak decomposition
The general form of the potential in the twocomponent spinor space has been classified by Okubo and Marshak[25]. We leave the derivation in Appendix B and recapitulate only the results here. By using the helmiticity, translational invariance in space and time, Galilei invariance, rotational invariance, parity and timereversal invariance, fermi statistics and isospin invariance, the potential has a general decomposition
(4.1)  
where is the projection operator to the isosinglet () and isotriplet :
(4.3) 
Also, we define
(4.4)  
(4.5)  
(4.6)  
(4.7)  
(4.8) 
with . The anticommutators in Eq. (4.2) are necessary to make the potential hermitian, since do not commute with the scalar potentials ().
If we keep the terms only up to the first order in , we obtain the conventional form of the potential at low energies commonly used in nuclear physics:
(4.9) 
or in a more conventional notation,
(4.11)  
The central and tensor potentials, and , in Eq. (4.11) are the leadingorder (LO) terms of in the velocity expansion, while the spinorbit potential, is the nexttoleadingorder (NLO) term of .
4.3 Determination of the potentials
For given , and , the matrix elements of the LO and NLO potentials up to in Eq. (4.11) have the following structure (see Appendix C and also see Ref. 26):
(4.12)  
(4.13)  
(4.14) 
with
(4.15)  
(4.16)  
(4.17) 
There are 8 unknown functions, , while we have 4 (2) diagonal and 1(0) offdiagonal matrix elements at each for as seen from Table I. On the lattice, it is relatively unambiguous to extract information for using the irreducible representations of the cubic group [15]. Then, at most 16 independent (14 diagonal and 2 offdiagonal) information as seen in Table I are obtained for 8 unknowns , so that each can be determined in two different ways.
4.4 Long range part of the potential
In QCD with dynamical quarks, the lightest hadron is the pion. Therefore, the longest range interaction between the nucleons is dictated by the onepionexchange potential (OPEP). For later purpose, let us here summarize several features of OPEP with special care about its chiral behavior.
First of all, the equivalence theorem implies that the pseudoscalar coupling and the pseudovector coupling at low energy are related through . This is simply obtained by kinematics. On the other hand, chiral symmetry leads to the GoldbergerTreiman (GT) relation, , where is the nucleon axialcharge and MeV) is the pion decay constant.
With these relations, the OPEP reads
(4.18)  
(4.19)  
(4.21)  
Here we have used the equivalence theorem to obtain Eq. (4.19) from Eq. (4.18) and use the GT relation to obtain Eq. (4.21) from Eq. (4.19). and in Eq. (4.21) are the values in the chiral limit.
In quenched QCD without dynamical quarks, there arises a dipole ghost in the flavorsinglet channel (the channel in the case of two flavors) which couples to the nucleons [34, 35]. The propagator in the quenched approximation is written as
(4.22) 
where with and being ghost parameters. The second term is the dipole ghost corresponding to the hairpin diagram with quarkline disconnected. Then the potential from the exchange reads [35]
where () is the pseudovector (pseudoscalar) coupling of the flavorsinglet to the nucleon. Its magnitude does not necessarily be as large as the coupling [36]. Note that the long range part of the potential has exponential falloff instead of the Yukawatype because of the dipoleterm in Eq. (4.22).
Let us define a ratio between the central potential in the spinsinglet channel and that in the spintriplet channel,
(4.25) 
Since we have , and the similar relations for the isospin, the large behavior of has different sign and magnitude between the oneghostexchange and onepionexchange. Therefore can be used as a tool to identify the ghost contribution at large distance as will be discussed in §6.5.
5 Central and tensor forces in lattice QCD
5.1 BS wave function on the lattice
To define the BS wave function on the lattice with the lattice spacing and the spatial lattice volume , we start from the fourpoint correlator,
(5.1)  
(5.2)  
(5.3)  
(D.11)  
(D.12)  
(D.13)  
(D.14)  
(D.15)  
(D.16)  
(D.17)  
(D.18)  
(D.19)  
(D.20)  
(D.21)  
(D.22)  
(D.23)  
(D.24)  
(D.25)  
(D.26)  
(D.27)  
(D.28)  
(D.29)  
(D.30)  
(D.31)  
(D.32) 
with the matrix element . The states created by the source have the conserved quantum numbers, (total angular momentum and its zcomponent) and (parity). For studying the nuclear force in the () channel and the ( and ) channel, we adopt a wall source located at with the Coulomb gauge fixing only at :
(5.4) 
where and are obtained by replacing the local quark fields and in Eqs. (3.2) and (3.2) by the wall quark fields,
(5.5) 
By construction, the source operator Eq. (5.4) has zero orbital angular momentum at , so that states with fixed are obtained by the spin projection with , e.g. and . Note that the and are not separately conserved: Therefore, the state created by the source becomes a mixture of the and at later time .
The BS wave function in the orbital Sstate is then defined with the projection operator for the orbital angular momentum () and that for the spin ():