Theoretical Foundation of the Nuclear Force in QCD
and Its Applications to Central and Tensor Forces
in Quenched Lattice QCD Simulations
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 well-understood theoretically, a large number of proton-proton and neutron-proton scattering data as well as deuteron properties have been accumulated and summarized e.g. in the Nijmegen database . To describe the elastic nucleon-nucleon () scattering at low-energies 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 : 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 Lippmann-Schwinger equation for the -matrix. Once the potential is determined, it can be used to study systems with more than 2 nucleons by using various many-body techniques.
Phenomenological potentials which can fit the data precisely (e.g. more than 2000 data points with ) at MeV are called the high-precision potentials: They include the potentials such as CD-Bonn , Argonne , and Nijm I, Nijm II and Reid93 . 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 :
The long range part of the nuclear force (the relative distance fm) is dominated by the one-pion exchange introduced by Yukawa . Because of the pion’s Nambu-Goldstone character, it couples to the spin-isospin density of the nucleon and hence leads to a strong spin-isospin dependent force, namely the tensor force.
The medium range part ( fm) receives significant contributions from the exchange of two-pions () and heavy mesons (, , and ). In particular, the spin-isospin 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 . 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 .
There is also a strong attractive spin-orbit 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 .
A repulsive core surrounded by an attractive well is in fact a common feature of the “effective” potential between composite particles. The Lenard-Jones potential between neutral atoms or molecules is a well-known 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 , exchange of non-linear pion field , and a combination of the Pauli principle with the one-gluon-exchange between quarks . 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 : 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 two-particle Bethe-Salpeter (BS) wave function in the interval under the periodic boundary conditions has sufficient information to relate the phase shift and the two-particle spectrum.
Lüscher’s method bypasses the difficulty to treat the real-time 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 Born-Oppenheimer picture, i.e., if one of the three quarks inside the baryon is infinitely heavy, one may define the potential between baryons a la Born-Oppenheimer . This is, however, not applicable to the nucleons with light quarks. The other employs the strong coupling limit, which has been proposed quite recently . Furthermore, it utilizes the finiteness of the lattice box effectively to extract the information of the on-shell scattering matrix and the phase shift. This approach has been applied to extract the scattering lengths in the quenched QCD simulations  and in the (2+1)-flavor QCD simulations with the mixed action .
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 energy-independent non-local 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 well-studied in connection with the BEC-BCS crossover in cold fermionic atoms , where the external magnetic field plays a role of the quark mass in QCD. For seminal suggestion on the rapid quark-mass dependence of the scattering length, see Ref. 23.
Since we consider the non-asymptotic region () of the wave function, the resultant potential and the -matrix are off-shell. 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 one-to-one 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 many-body problems but also to hyperon-nucleon, hyperon-hyperon and three-nucleon interactions which have much less experimental information than the systems. A first attempt to the hyperon-nucleon 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 two-body and many-body 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 field-theoretical 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 Okubo-Marshak decomposition) is reviewed. In Appendix C, matrix elements of the general potential are presented. In Appendix D, heat-kernel representation of the Green’s function is presented.
2 Non-local potential in quantum mechanics
2.1 Two-body force
To show the basic concept of the non-local potential in a finite box with the size , we start with a non-relativistic two-body problem described by the stationary Schrödinger equation:
where is the relative coordinate of the two spinless and non-relativistic particles, and is related to the discrete energy eigenvalues ( with being the reduced mass. The wave function obeys the periodic boundary condition. The non-local potential ‡)‡)‡)Here we use the standard term “non-local” in the sense that cannot be written as . is assumed to be energy-independent 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:
where is the non-relativistic kinetic energy operator satisfying . Since removes the non-interacting part of the wave function, is non-vanishing 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 energy-independent and non-local potential can be defined as
which leads to the Schrödinger equation Eq. (2.1) for . If we apply a unitary transformation to the wave function, , the non-local potential is modified as . Such unitary transformation does not affect the observables, while it changes the spatial structure of the wave function and the non-local 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 equal-time Bethe-Salpeter 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 low-energy scattering with sufficiently smaller than the intrinsic scale of the system or the scale of the non-locality of the potential. In such a case, the velocity expansion of in terms of its non-locality is useful : For example, a spin-independent potential with hermiticity, rotational invariance, parity symmetry, and time-reversal invariance can be expanded as
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
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  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 Many-body forces
For the interactions among composite particles, there are in principle many-body forces which take place in the system composed of more than two particles. The well-known example in nuclear physics is the Fujita-Miyazawa type three-body force acting among three nucleons [27, 28]. It is phenomenologically important for the extra binding of light nuclei  and for the extra repulsion in high density matter  and in elastic nucleus-nucleus scatterings .
The method to define the two-body potential from the relative wave function discussed above can be generalized to the many-body forces. Let us illustrate the procedure by considering the three-body system of spinless and distinguishable particles with equal mass . We consider the local potentials for both two-body and three-body forces just for simplicity. In the rest frame of the three-body system, we have
where and are the Jocobi coordinates. and are the kinetic energy operator with and being the reduced masses. is the total energy of the three-body system at rest. Because of the translational invariance, the two-body potential and the three-body potential are the functions of and .
If we know the wave function and the total energy on the left-hand side of Eq. (2.2), the three-body potential can be determined by the following procedure. We first consider the situation, , where , and are vanishingly small because of the assumed short-range 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 two-body system.
Once all the two-body potentials are determined, can be extracted from the wave function in the range, and , through the three-body equation Eq. (2.2). It is important to note that the three-body potential is always obtained together with the two-body 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 non-local potentials and to the -particle systems with different masses and internal degrees of freedom.
3 Non-local potential in field theory for spin 1/2 particles
3.1 Bethe-Salpeter wave function
In field theory, the best analogue of the two-particle wave function is the equal-time Bethe-Salpeter (BS) amplitude, so that we use the term “BS wave function” throughout this paper. Let us consider the following BS wave function for the 6-quark state with total energy and the total three-momentum in a finite box ,
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 “neutron-like” three-quarks located at point and “proton-like” three-quarks 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 non-relativistic case, we introduce the “effective center of mass energy”, . 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 short-ranged non-local potential as
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 two-nucleon system.
3.2 Interpolating operators
In Eq. (3.1), simplest interpolating operators for the neutron and the proton written in terms of the up-quark field and the down-quark fields would be
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 four-point 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 .
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 non-local 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 one-to-one 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
It is useful to classify the asymptotic two-particle 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 well-known relations given in Table I.
4.2 Okubo-Marshak decomposition
The general form of the potential in the two-component spinor space has been classified by Okubo and Marshak. 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 time-reversal invariance, fermi statistics and isospin invariance, the potential has a general decomposition
where is the projection operator to the iso-singlet () and iso-triplet :
Also, we define
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:
or in a more conventional notation,
The central and tensor potentials, and , in Eq. (4.11) are the leading-order (LO) terms of in the velocity expansion, while the spin-orbit potential, is the next-to-leading-order (NLO) term of .
4.3 Determination of the potentials
There are 8 unknown functions, , while we have 4 (2) diagonal and 1(0) off-diagonal 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 . Then, at most 16 independent (14 diagonal and 2 off-diagonal) 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 one-pion-exchange 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 pseudo-scalar coupling and the pseudo-vector coupling at low energy are related through . This is simply obtained by kinematics. On the other hand, chiral symmetry leads to the Goldberger-Treiman (GT) relation, , where is the nucleon axial-charge and MeV) is the pion decay constant.
With these relations, the OPEP reads
In quenched QCD without dynamical quarks, there arises a dipole ghost in the flavor-singlet 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
where with and being ghost parameters. The second term is the dipole ghost corresponding to the hairpin diagram with quark-line disconnected. Then the potential from the exchange reads 
where () is the pseudo-vector (pseudo-scalar) coupling of the flavor-singlet to the nucleon. Its magnitude does not necessarily be as large as the coupling . Note that the long range part of the potential has exponential fall-off instead of the Yukawa-type because of the dipole-term in Eq. (4.22).
Let us define a ratio between the central potential in the spin-singlet channel and that in the spin-triplet channel,
Since we have , and the similar relations for the isospin, the large behavior of has different sign and magnitude between the one-ghost-exchange and one-pion-exchange. 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 four-point correlator,
with the matrix element . The states created by the source have the conserved quantum numbers, (total angular momentum and its z-component) 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 :
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 S-state is then defined with the projection operator for the orbital angular momentum () and that for the spin ():