Hadron interactions in lattice QCD

Progress on the potential method, recently proposed to investigate hadron interactions in lattice QCD, is reviewed. The strategy to extract the potential in lattice QCD is explained in detail. The method is applied to extract $NN$ potentials, hyperon potentials and the meson-baryon potentials. A theoretical investigation is made to understand the origin of the repulsive core using the operator product expansion. Some recent extensions of the method are also discussed.


Introduction
Quarks and gluons, the fundamental degrees of freedom of Quantum Chromodynamics (QCD) never appear in Nature as asymptotic one-particle states. This phenomena is called confinement. Instead, only their bound states, i.e., mesons and baryons, are observed in experiment and these hadronic states are the asymptotic states whose interactions can be parametrized by the S-matrix. An ultimate goal for theoretical studies of the strong interaction, therefore, is to extract the properties of the hadronic S-matrix from QCD.
For a description of hadronic interactions, the nuclear force is one of the most fundamental quantities in nuclear physics. The origin of the nuclear force, however, is still one of the major unsolved problems in strong interaction physics even after the establishment of QCD. Although the nuclear force is not well understood theoretically, a large number of proton-proton and neutron-proton scattering data as well as deuteron properties have been experimentally accumulated and summarized as the nucleon-nucleon (NN) potential, which is designed to reproduce these experimental properties through the Schrödinger equation for the NN wave function. Once the potential is determined, it can be used to study systems with more than 2 nucleons by using various many-body techniques.
Phenomenological NN potentials which can fit the NN data precisely (e.g. more than 2000 data with χ 2 /dof ≃ 1 ) at T lab < 300 MeV are called the high-precision NN potentials. Those in coordinate space, some of which are shown in Fig. 1, reflect characteristic features of the NN interaction at different length scales [1,2,3,4,5]: (i) The long range part of the nuclear force (at relative distance r > 2 fm) is dominated by one-pion exchange introduced by Yukawa [9]. Because of the pion's Nambu-Goldstone character, it couples Figure 1: Three examples of the modern NN potential in 1 S 0 (spin-singlet and S-wave) channel: CD-Bonn [6], Reid93 [7] and Argonne v 18 [8]. Taken from Ref. [15].
to the spin-isospin density of the nucleon and hence leads to a strong spin-isospin dependent force, namely the tensor force.
(ii) The medium range part (1 fm < r < 2 fm) receives significant contributions from the exchange of two-pions (ππ) and/or 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.
(iii) The short range part (r < 1 fm) is best described by a strong repulsive core as originally introduced by Jastrow [10]. 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 [11,12,13].
It is then a challenge for the theoretical particle and nuclear physics communities to extract these properties of the NN interaction from first principle non-perturbative QCD calculations, in particular lattice QCD simulations. A theoretical framework suitable for such a purpose was first proposed by Lüscher [14]. For two hadrons in a finite box with the size L × L × L with 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 R is sufficiently smaller than the size of the box R < L/2, the behavior of the two-particle Nambu-Bethe-Salpeter (NBS) wave function ψ(r) in the interval R < |r| < L/2 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. Furthermore, it utilizes the finiteness of the lattice box effectively to extract the information of the on-shell scattering matrix and the phase shift.
A closely related but an alternative approach to the NN interactions from lattice QCD has been proposed recently [15,16,17]. The starting point is the same NBS wave function ψ(r) as discussed in Ref. [14]. Instead of looking at the wave function outside the range of the interaction, the authors consider the internal region |r| < R and define an energy-independent non-local potential U(r, r ′ ) from ψ(r) so that it obeys the Schrödinger type equation in a finite box. Since U(r, r ′ ) for strong interactions is localized in its spatial coordinates due to confinement of quarks and gluons, the potential receives only weak finite volume effect in a large box. Therefore, once U is determined and is appropriately extrapolated to L → ∞, one may simply use the Schrödinger equation in infinite space to calculate the scattering phase shifts and bound state spectra to compare with experimental data. A further advantage of utilizing the potential is that it is 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 in the quark mass corresponding to the formation of the NN bound state.
Since the recent progress for the study of the NN interaction by the first method has already been reviewed in Ref. [18], the recent progress for the second method is mainly considered in this review. In Sec. 2, the strategy of Ref. [15,16,17] to define the NN potential in QCD is explained in detail, and the lattice formulation is introduced in Sec. 3. Results of lattice QCD calculations for NN potentials are given in both quenched and full QCD in Sec. 4. Central potentials at the leading order of the velocity expansion is shown to reproduce qualitative features of the NN potential such as the repulsion at short distance and the attraction at medium to long distances. The tensor potential, which exists also at leading order, is extracted. Contrary to the case of the central potentials, it does not have a repulsive core. Higher order contributions in the velocity expansion are also investigated and shown to be small at low energy and low orbital angular momentum L. In Sec. 5, the method to extract the potential is applied to the hyperon-nucleon interactions such as NΞ and NΛ systems. Interactions between octet baryons in general are also investigated in the flavor SU(3) limit, where up, down and strange quark masses are all equal. In Sec. 6, we also consider a recent attempt to understand the origin of the repulsive core in the NN potential. Using the operator product expansion and renormalization group analysis in QCD, the potential derived from the NBS wave function in Sec. 4 is shown to have a repulsive core, whose functional form is also theoretically predicted. In Sec. 7, two extensions of the potential method are considered, together with explicit applications of these extensions to hadron interactions. One is the extension of the potential method to inelastic scattering, in order to investigate the ΛΛ system, while the other is the extraction of the potential from the time dependent NBS wave function in lattice QCD. With the latter method, the existence of the H-dibaryon is investigated in the flavor SU(3) limit. In Sec. 8, applications of the method to the three nucleon force, meson-baryon potentials and the potential in 2-color QCD are considered. Brief concluding remarks are given in Sec. 9.
2 Strategy to extract potential in QCD 2.1 Nambu-Bethe-Salpeter (NBS) wave function and its asymptotic behavior A key quantity to extract the potential from QCD is the equal time Nambu-Bethe-Salpeter wave function, defined by where |2N, W, s 1 s 2 is an eigen-state of QCD for two nucleons with total energy W = 2 k 2 + m 2 N and the total three-momentum P = 0, whose helicities are denoted by s 1 , s 2 . Here the local nucleon operator is given by where x = (x, t), the color indices are denoted by a, b, c, and α is the spinor index. The charge conjugation matrix in the spinor space is given by C = γ 2 γ 4 , and p, n denotes proton and neutron operators while u, d denote up and down quark operators. Note that ϕ W implicitly has two pairs of spinor-flavor indices which come from N α (r + x, t)N β (r, t) and two helicity indices s 1 and s 2 . The most remarkable property of the above NBS wave function is explained as follows. If the W is smaller than the threshold energy for one-pion production (i.e. W < 2m N + m π ), then its asymptotic behavior for large |x| can be evaluated [19,17]. The helicity component in the spin singlet channel (S = 0) is given by |s 1 s 2 = 1 √ 2 | + 1 2 , + 1 2 + | − 1 2 , − 1 2 , where the relative + sign is our convention. For this case, we have ϕ W (r) S=0 ≃ l,lz Z l,lz (S = 0)Y llz (Ω r ) sin(kr − lπ/2 + δ l0 (k)) kr e iδ l0 (k) where r = |r|, k = |k|, and δ lS (k) is the NN scattering phase shift in QCD with the total angular momentum l and the total spin S, which is determined by the unitarity of S-matrix in QCD below the inelastic threshold [17]. Here Y lm (Ω r ) is the spherical harmonic function with the solid angle Ω r of r.
The coefficient Z llz (0) has spinor components α, β and is given by where Z is the wave function renormalization for the nucleon operator N(x), D l mλ is the Wigner Dmatrix, and U αα (∇) and U ββ (−∇) are the 4×2 matrices acting on the 2×2 matrix χαβ(S, S z ). Explicitly we have For the spin triplet channel (S = 1), the asymptotic behavior of ϕ W is more involved but schematically written as ϕ W (r) S=1 ∝ Y llz (Ω r ) sin(kr − lπ/2 + δ l1 (k)) kr e iδ l1 (k) .
An explicit form of the asymptotic behavior for the NBS wave function in the triplet channel is given in appendix A. The asymptotic behaviors in eqs. (3) and (7) tell us that the NBS wave function at large separation r describes the scattering wave of the quantum mechanics whose phase shift agrees with the phase of the S-matrix in QCD. Therefore the NBS wave function satisfies the free Schrödinger equation at large r as at W < 2m N + m π , where µ = m N /2 is the reduced mass of the NN system. Note that these properties hold without using the non-relativistic approximation/expansion. In particular, only the upper components of the spinor indices for the NBS wave function (α = 1, 2 and β = 1, 2) are enough to reproduce all NN scattering phase shifts δ lS (k) with l = 0, 1, 2, 3, · · · and S = 0, 1. From these properties, (the upper spinor components of) the NBS wave function can be regarded as the "wave function" of the NN system at W < 2m N + m π . It is also noted that the equal-time constraint for the NBS wave function here is not a restriction to the extraction of physical observables such as the scattering phase shift, as evident from the fact that all informations on the scattering phase shifts are encoded in the asymptotic behavior of the equal-time NBS wave function. Moreover, the equal-time NBS wave function with non-zero total momentum (P = 0) is equivalent to the NBS wave function in the center of mass frame with non-zero time separation. In this case the 4-dimensional distance between two nucleon operators is always space-like.

Non-local potential from the NBS wave function
Since the NBS wave function satisfies the free Schrödinger equation at large r, one can define shortranged non-local potential by It is noted that the spinor indices α, β, γ, δ here run from 1 to 2, since all NN scattering phase shifts can be reproduced from the NBS wave function with α, β ∈ {1, 2} as discussed in the previous subsection. Therefore U αβ;γδ has 4×4 components, which can be determined from 4 components of ϕ W αβ for 4 different combinations of (s 1 , s 2 ). Note that since the NBS wave function ϕ W is multiplicatively renormalized, the potential U(x, y) is finite and does not depend on the particular renormalization scheme .
The non-local function U(x, y) is shown to be energy-independent as follows. Let V th be the space spanned by the wave function with W ≤ W th ≡ 2m N + m π : V th = {ϕ W |W ≤ W th }, and the projection operator to V th is defined as where ρ(W ) is the density of states at energy W , and N −1 (W 1 , W 2 ) is the inverse of the hermitian operator N(W 1 , W 2 ) defined by The non-local potential is then defined by It is easy to see that the above non-local potential satisfies eq.(9) at W ≤ W th as follows.
It should be noted that the non-local potential U(x, y) which satisfied eq.( 9) at W ≤ W th is not unique: For example, we can add the following term with arbitrary functions f W (x) to the non-local potential U(x, y) without affecting eq.(9) at W ≤ W th . The non-local potential U(x, y) in eq. (13) is energy independent by construction. Alternatively we can define a different non-local potential by which satisfies eq. (9) for all W . This potential, however, becomes long-ranged, due to the presence of inelastic contributions above W th . The extension of this method to non-elastic cases will be discussed in Sec. 7. In Ref. [18], it is claimed that the NBS wave function satisfies the Schrödinger equation with a non-local and energy-dependent potential. In general this is true but, as shown here, there is a scheme which makes the non-local potential energy-independent. This is sufficient for the strategy considered in this report.

Velocity expansion of the non-local potential
In principle, if one knows all NBS wave functions ϕ W , the non-local potential U can be constructed according to eq. (13). In practice, however, one can obtain only a few of them corresponding to the ground state as well as a few low lying excited states in lattice QCD simulations. Therefore, for practical applications, it is convenient to expand the non-local potential in terms of the velocity(derivative) with local functions as At the lowest few orders we have where r = |r|, σ i is the Pauli-matrix acting on the spin index of the i-th nucleon, S = ( σ 1 + σ 2 )/2 is the total spin, L = r × p is the angular momentum, and is the tensor operator. Each coefficient function is further decomposed into its flavor components as where τ i is the Pauli-matrix acting on the flavor index of the i-th nucleon. The form of the velocity expansion (18) agrees with the form determined by symmetries [20]. For the leading order of the velocity expansion, the local potential is given by which can be obtained from the NBS wave function at one value of W . Since S 12 = 0 for the spin single state, for example, we have

Remarks
There are several remarks on the extraction of the potential from QCD described in the previous sections.
First of all, it should be noted that the potential itself is not a physical observable, and it is therefore not unique. In particular the potential depends on the choice of the nucleon operator to define the NBS wave function. The local nucleon operator is one choice here but other definitions are equally possible, though the local operator is a convenient choice for the reduction formula of composite particles such as the nucleon [21,22,23]. A choice of the nucleon operator to define the NBS wave function is considered to be a "scheme" to define the potential. The potential is therefore a scheme dependent quantity, while physical observables such as the scattering phase shift and the binding energy of the deuteron are of course scheme independent.
Is such a scheme-dependent quantity useful ? The answer to this question is probably "yes", since the potential is useful to "understand" or "describe" the phenomena. For example, the repulsive core best summarizes the behavior of the NN scattering phase shift at larger energy in terms of the short distance behavior of the potential. A well-known example for a scheme-dependent but useful quantity is the running coupling in quantum field theory. Although the running coupling is of course scheme dependent, it is useful to understand the deep inelastic scattering data at high energy (i.e. asymptotic freedom).
Among different schemes, there exist of course good schemes. A good convergence of the perturbative expansion for a certain class of observables is a reasonable criteria for a good scheme in the case of the running coupling. In the case of the potential, a good convergence of the velocity expansion is important. In this respect, a completely local and energy-independent potential is the best one, and moreover the inverse scattering method tells us that it must be unique if it exists. We think that such a local and energy-independent potential exists only if no inelastic scattering appears at all values of energy (i.e. W th = ∞).
There exists a criticism that the repulsive core in the phenomenological potentials is meaningless since the low energy scattering data can be described equally well by other potentials without the repulsive core. In other words, the repulsive core can be removed by a unitary transformation of the wave function without changing the scattering phase shift at low energy. Even though this may be true, the resulting "potential" is highly non-local, and therefore is beyond the concept of a potential. This criticism thus corresponds to claiming that the asymptotic freedom is meaningless since it can be removed by some non-perturbative scheme of the "coupling", which is however beyond the concept of a coupling.
For the definition of the potential in QCD from the NBS wave function, the convergence of the velocity expansion can be checked by examining the energy (W ) dependence of the lower order potentials. For example, if we have ϕ Wn for n = 1, 2, · · · N, we can determine the N − 1 unknown local functions of the velocity expansion in N different ways. The variation among N different determinations gives an estimate of the size of the higher order terms neglected. Furthermore one of these higher order terms can be determined from ϕ Wn for n = 1, 2, · · · N. The convergence of the velocity expansion will be investigated explicitly in Sec. 4.
Let us consider what we have shown so far from a slightly different point of view. Instead of adopting the point of view that the "potential" can be defined in QCD, the analysis in this section shows that the use of quantum mechanics with potentials to describe the NN scattering can be justified in a more fundamental quantum field theory (QCD) through the NBS wave function, whose asymptotic behavior encodes phases of the S-matrix for the NN scattering. If the local approximation for the potential defined from the NBS wave function is good in the low energy region, one can use the local potential combined with quantum mechanics to investigate nuclear physics 1 .
It is now clear that there is no unique definition for the NN potential. Ref. [18,25,26], however, criticized that the NBS wave function is not "the correct wave function for two nucleons" and that its 1 It is more legitimated to show if three or more body potentials remain subdominant in the same framework in QCD before using the potential in nuclear physics. There exists another approach to investigate nuclei, by combining the chiral effective field theory with lattice calculations. (See Ref. [24] and references therein.) The chiral effective can easily incorporate may-body interactions order by order in the chiral expansion, though there are many free parameters to be determined from experimental data. On the other hand, the potential from lattice QCD does not contain such free parameters. Therefore it is interesting and important to determine some parameters of the chiral effective theory from the potential in lattice QCD. relation to the correct wave function is given by where N 0 (x, t) is "a free-field nucleon operator" and the ellipses denote "additional contributions from the tower of states of the same global quantum numbers". Thus 0|T {N 0 (x+r, 0)N 0 (x, 0)}|2N, W, s 1 , s 2 is considered to be "the correct wave function". In this claim it is not clear what is "a free-field nucleon operator" in the interacting quantum field theory. An asymptotic in or out field operator may be a candidate. If the asymptotic field is used for N 0 , however, the potential defined from the wave function identically vanishes for all r by construction. To be more fundamental, a concept of "the correct wave function" is doubtful. If some wave function were "correct", the potential would be uniquely defined from it. This clearly contradicts the fact discussed above that the potential is not an observable and therefore is not unique. This argument shows that the criticism of Ref. [18,25,26] is flawed.

Lattice formulation
In this section, we discuss the procedure to extract the NBS wave function from lattice QCD simulations. For this purpose, let us consider the correlation function on the lattice defined by where J (t 0 ) is a source operator which creates a two nucleon state of states and its explicit form will be considered later. By inserting a complete set and considering baryon number conservation, we have For a large time separation (t − t 0 ) → ∞, we have where W 0 is assumed to be the lowest energy of NN states. Since the source dependent term A 0 is just a multiplicative constant to the NBS wave function ϕ W 0 (r), the potential defined from ϕ W 0 (r) in the previous section is manifestly source-independent. Therefore the statement that the potential in this scheme is "source-dependent" in Ref. [27] is clearly wrong. In this extraction of the wave function, the ground state saturation for the correlation function F in eq. (26) is important. In principle, one can achieve this by taking a large t − t 0 . In practice, however, F becomes very noisy at large t − t 0 , so that the extraction of ϕ W 0 becomes difficult at large t − t 0 . Therefore it is crucial to find the region of t where the ground state saturation is approximately satisfied while the signal is still reasonably good. The choice of the source operator becomes important in order to have such a good t-region.

Choice of source operator
We can choose the source operatorJ to fix quantum numbers of the state |2N, W, s 1 , s 2 such as (J, J z ). Since a lattice QCD simulation is usually performed on the (finite) hyper-cubic lattice, we consider the cubic transformation group SO(3, Z) instead of the SO(3, R) as the symmetry of 3-dimensional space. Therefore the quantum number is classified in terms of the irreducible representation of SO(3, Z), Table 1: The number of each representation of SO(3, Z) which appears in the angular momentum L representation of SO(3, R). P = (−1) L represents an eigenvalue under parity transformation.
denoted by A 1 , A 2 , E, T 1 , T 2 whose dimensions are 1, 1, 2, 3, 3. A relation of irreducible representations between SO(3, Z) and SO(3, R) is given in table 1 for L ≤ 6, where L represents the angular momentum for the irreducible representation in SO(3, R). For example, the source operatorJ (t 0 ) in the A 1 representation with positive parity generates states with L = 0, 4, 6, · · · at t = t 0 , while the operator in the T 1 representation with negative parity produces states with L = 1, 3, 5, · · ·. For two nucleons, the total spin S becomes 1/2 ⊗ 1/2 = 1 ⊕ 0, which corresponds to T 1 (S = 1) and A 1 (S = 0) of the SO(3, Z). Therefore, the total representation J for a two nucleon system is determined by the product R 1 ⊗ R 2 , where R 1 = A 1 , A 2 , E, T 1 , T 2 for the orbital "angular momentum" while R 2 = A 1 , T 1 for the total spin. In table 2, the product R 1 ⊗ R 2 is decomposed into the direct sum of irreducible representations. Most of the results in this report use the wall source at t = t 0 defined by where α, β = 1, 2 are upper component spinor indices while f, g are flavor indices. Here N wall (t 0 ) is obtained by replacing the local quark field q(x) of N(x) by the wall quark field, with the Coulomb gauge fixing only at t = t 0 . Note that this gauge-dependence of the source operator disappears for the potential. All states created by the wall source have zero total momentum. Among them the state with zero relative momentum has the largest magnitude. The most important reason to employ the wall source is that the ground state saturation for the potential at long distance is better achieved for the wall source than other sources such as the smeared source. By construction, the source operatorJ wall (t 0 ) has zero orbital angular momentum at t = t 0 , which corresponds to the A 1 representation with positive parity. By the spin projection operator P (S) , e.g. P (S=0) = σ 2 and P (S=1,Sz=0) = σ 1 , we fix the J of the source as J (t 0 ; J P =+ , I) = P (S) βα J wall (t 0 ) αβ,f g (29) where P = ± is the parity and I = 1, 0 is the total isospin of the system. Since the nucleon is a fermion, an exchange of two nucleon operators in the source should give a minus sign. This fact fixes the total isospin once the total spin is given: (S, I) = (0, 1) or (1, 0). (Note that S, I = 0 are antisymmetric while S, I = 1 are symmetric under the exchange.) Since A + 1 ⊗ A 1 (S = 0) = A + 1 and A + 1 ⊗ T 1 (S = 1) = T + 1 , the state with either (J P , I) = (A + 1 , 1) for the spin-singlet or (J P , I) = (T + 1 , 0) for the spin-triplet is created at t = t 0 by the corresponding source operator. The NBS wave function extracted at t > t 0 has the same quantum numbers (J P , I) as they are conserved under QCD interactions. In addition the total spin S is conserved at t > t 0 for the two-nucleon system with equal up and down quark masses: Under the exchange of the two particles, the constraint (−1) S+1+I+1 P = −1 must be satisfied due to the fermionic nature of nucleon, while the parity P and the isospin I are conserved in this system. Therefore S is conserved. It is noted, however, that L is not conserved in general. While the state with (J P , I) = (A + 1 , 1) always has L = A + 1 even at t > t 0 , the one with (J P , I) = (T + 1 , 0) has both L = A + 1 and L = E + components 2 at t > t 0 , which corresponds to L = 0 and L = 2 in SO(3, R), respectively. Note that J or L in this report is used to represent the total or orbital quantum number of SO(3, Z) as well as SO(3, R), depending on the context. The orbital angular momentum L of the NBS wave function can be fixed to a particular value by the projection operator P (L) as ϕ W (r; J P , I, L, S) = P (L) P (S) ϕ W (r; J P , I) where ϕ W (r; J P , I) is extracted from for large t − t 0 . Here we also apply the total spin projection operator P (S) but this is redundant since the total spin S, already fixed by the source, is conserved as mentioned above. The projection operator P (L) to an arbitrary function φ(r) is defined in general by for L = A 1 , A 2 , E, T 1 , T 2 , where χ L denotes the character for the representation L of the cubic group SO(3, Z), g is one of 24 elements in SO(3, Z) and d L is the dimension of L.

Leading order potential: spin-singlet case
We present the procedure to determine potentials at the reading order(LO): Since S 12 = 0 and σ 1 · σ 2 = −3 for the spin-singlet case, the LO central potential in the case of (J P , I) = (A + 1 , 1) state is extracted as where V I=1 X = V 0 X + V τ X in isospin space. The potential V C (r) (S,I)=(0,1) in the above is often referred to as the central potential for the 1 S 0 state, where the notation 2S+1 L J represents the orbital angular momentum L (see table 1), the total spin S and the total angular momentum J of J = L + S. It is noted, however, that in the leading order of the velocity expansion, the potential does not depend on the quantum number of the state J = L = A 1 . Moreover the A 1 state may contain L = 4, 6, · · · components other than L = 0, though the L = 0 component may dominate in the state. Therefore it is more precise to refer to V C (r) (S,I)=(0,1) as the spin-singlet (isospin-triplet) central potential determined from the state with J = L = A 1 . A possible difference of spin-singlet central potentials between this determination and others such as the one determined from J = L = E gives an estimate for contributions from higher order terms in the velocity expansion.

Leading order potential: spin-triplet case
Both the tensor potential V T and central potential V C appear at the LO for the spin-triplet case. Let us consider the determination from the (J P , I) = (T + 1 , 0) state. The Schrödinger equation for this state becomes where the spin-triplet central potential is given by The projections to A 1 and E components read The last quantity in eq. (38) is an approximation of the first line and a difference comes from T 1 and T 2 components, which are expected to be small. This approximate representation for Q is often employed in numerical simulations. Using these projections, V C and V T can be extracted as In numerical simulations, (α, β) = (2, 1) is mainly employed. If one neglects V T by putting V T = 0 in the above, one obtains the effective central potential for the spin-triplet (isospin-singlet) state as The difference between V C and V eff C is O(V 2 T ) in the second order perturbation for small V T .

A comparison with the finite volume method in lattice QCD
In this subsection we briefly compare the potential method with the direct extraction of the phase shift via the finite volume method in lattice QCD. First of all, by construction, the potential method gives the correct phase shift at where W is the total energy of the state from which the NBS wave function is defined, while phase shifts at other values of k are approximated ones obtained in the velocity expansion of the non-local potential. Secondly, the finite size correction to the potential is expected to be small. Indeed the finite volume method for the extraction of the phase shift in lattice QCD assumes that there is no finite size correction to the potential as long as the volume is large enough so that the interaction range of the potential is smaller than half of the lattice extension, L/2. Under this condition, there exists an asymptotic region in the periodic box where the scattering wave satisfies the free Schrödinger equation with a specific value of the energy, from which one can determine the phase shift at certain values of k in the infinite volume. This is Lüscher's finite volume formula for the phase shift [14].
Thirdly, we also expect that the quark mass dependence of the potential is much milder than that of physical observables such as the scattering length. While the scattering length is almost zero at the heavy quark mass region, it diverges when the bound state is formed at the lighter quark mass region. In this situation, the scattering length varies from zero to infinity as the quark mass changes [28]. Such a drastic change of the scattering length can easily be realized by a small change to the shape of the potential as a function of the quark mass.
Let us assume that higher order terms in the velocity expansion give negligible contributions at low energy so that the leading order local potential well reproduces the scattering phase shift. In this situation, some problems of the finite size method can be avoided by using the potential method. To extract the phase shift in the finite size method in lattice QCD, one has to assume that one particular angular momentum gives a dominant contribution among possible angular momenta in a given representation of the cubic group. For example, although a state in the A 1 representation contains not only an L = 0 contribution but also L = 4, 6, · · · contributions, one usually assumes that the L = 0 contribution dominates so that the energy shift in the finite volume is related to the scattering phase for the L = 0 state. In the case of the potential, on the other hand, such an assumption is unnecessary. One can determine the local potential in the velocity expansion from the A 1 state without specifying the dominant angular momentum. Once the potential is obtained, one can calculate the scattering phase shift for an arbitrary L by solving the Schrödinger equation in the infinite volume with the extracted potential. Furthermore, one can check the assumption made for the finite size method by comparing sizes of the scattering phases among different L's.
4 Results for nuclear potentials from lattice QCD

Quenched QCD results for (effective) central potentials
Let us show results in the quenched QCD, where creations and annihilations of virtual quark-antiquark pairs are all neglected. For the simulations, the standard plaquette gauge action is employed on a 32 4 lattice at the bare gauge coupling constant β = 6/g 2 = 5.7, which corresponds to the lattice spacing a ≃ 0.137 fm (1/a = 1.44(2) GeV), determined from the ρ meson mass in the chiral limit, and the physical size of the lattice L ≃ 4.4 fm [15]. As for the quark action, the standard Wilson fermion action is employed at three different values of the quark mass corresponding to the pion mass m π ≃ 731, 529, 380 MeV and the nucleon mass m N ≃ 1560, 1330, 1200 MeV, respectively.  The central potential in the spin-singlet channel and the effective central potential in the spintriplet channel reconstructed from the wave functions at m π ≃ 529 MeV are shown in Fig. 2(Right). These potentials reproduce the qualitative features of the phenomenological NN potentials, namely the repulsive core at short distance surrounded by the attractive well at medium and long distances. From this figure one observes that the interaction range of the potential is smaller than 1.5 fm. Therefore the box size L ≃ 4.4 fm is large enough to extract the phase shift by the finite size method, and furthermore the finite size corrections to the potentials themselves are expected to be small. Labels 1 S 0 and 3 S 1 of the potentials in the figure represent the fact that potentials are determined from A 1 wave functions, which are dominated by S wave components.
Instead of calculating the energy shift due to the finite size, one can extract the asymptotic momentum k, by fitting the NBS wave function ϕ(r) at large distance with the Green's function G(r; k 2 ) in a finite and periodic box for the Helmholtz equation (∇ 2 + k 2 )G(r; k 2 ) = −δ lat (r) with δ lat (r) being the periodic delta-function. Explicitly it is given by The asymptotic momentum k is related to the scattering phase shift δ 0 (k) or the scattering length a 0 for the S-states 3 as where Z 00 (1, q 2 ) with q = kL 2π is (the analytic continuation of) the generalized zeta-function Z 00 (s,   Table 3: Effective center of mass energies E obtained from the asymptotic momenta and the scattering length a 0 at different pion masses. Taken from Ref. [17]. -0.68(26) -0.97(37) 0.15(7) 0.23 (10) 16a ≃ 2.2 fm using the above form at m π = 529 MeV. This leads to the values of the effective energy E ≡ k 2 /m N , which can be translated to the scattering length a 0 by the Lüscher's formula (44). In Fig.4, we compare the NN central potentials in the spin-singlet channel for three different pion masses. As the pion mass decreases, the repulsion at short distance and the attraction at medium distance are enhanced simultaneously. In table 3, we give values of E and the S-wave scattering length a 0 , which show a net attraction of the NN interactions in both channels at these pion masses, though the absolute magnitudes of the scattering length a 0 are much smaller than the experimental values at the physical pion mass m π ≃ 140 MeV: a The above discrepancy is partly caused by the heavier pion masses and the absence of the dynamical quarks in quenched simulations. If we get closer to the physical pion mass in full QCD simulations, there should appear the "unitary region" where the NN scattering length shows the singularity associated with the formation of the di-nucleon bound state, so that a 0 changes sign [28]. Therefore the NN scattering length must become a non-linear function of the pion mass in this region. Unlike the scattering length, on the other hand, the NN potential does not necessarily have singular behavior in the unitary region, as demonstrated in the well-known quantum mechanical examples such as the low-energy scattering between ultracold atoms. To check this in QCD, it is of course important to study the NN potential in the full QCD simulations at lighter pion masses.
In addition to the above reasoning, there is a possibility that extracted values of k 2 have large systematic uncertainties caused by the contamination of the excited states at large distance for the wave functions. (and J z = S z = 0) states at m π ≃ 529 MeV. (Right) The same wave functions but the spherical harmonics components are removed from the non-A + 1 part. Taken from Ref. [17].
These NN scattering lengths extracted from the NBS wave function agree in sign but are much smaller in magnitude than the previous quenched results from the finite size method in smaller volume [29], while they disagree even in sign with the recent full QCD results form the finite size method (See Ref. [18] and references therein.).

Tensor potential
In Fig. 5(Left), we show the A 1 and non-A 1 components of the NBS wave function obtained from the J P = T + 1 (and J z = S z = 0) states at m π ≃ 529 MeV, according to eqs. (37) and (38). The A 1 wave function is multivalued as a function of r due to its angular dependence. For example, the (α, β) = (2, 1) component of the L = 2 part of the non-A 1 wave function is proportional to the spherical harmonics Y 20 (θ, φ) ∝ 3 cos 2 θ − 1.   Figure 6: (Left) The central potential V C (r) (1,0) and the tensor potential V T (r) obtained from the J P = T + 1 NBS wave function, together with the effective central potential V eff C (r) (1,0) , at m π ≃ 529 MeV. (Right) Pion mass dependence of the tensor potential. The lines are the four-parameter fit using one-pion-exchange + one-rho-exchange with Gaussian form factor. Taken from Ref. [17]. Fig. 6 (Left) are the central potential V C (r) (1,0) and tensor potential V T (r) together with the effective central potential V eff C (r) (1,0) , at the leading order of the velocity expansion as given in eqs. (39), (40) and (42), respectively.

Shown in
Note that V eff C (r) contains the effect of V T (r) implicitly as higher order effects through processes such as 3 S 1 → 3 D 1 → 3 S 1 . At the physical pion mass, V eff C (r) is expected to obtain sufficient attraction from the tensor potential, which causes the appearance of a bound deuteron in the spin-triplet (and flavor-singlet) channel while an absence of the bound dineutron in the spin-singlet (and flavor-triplet) channel. The difference between V C (r) (1,0) and V eff C (r) in Fig. 6 (Left) is still small in this quenched simulation due to relatively large pion mass. This is also consistent with the small scattering length in the previous subsection.
The tensor potential in Fig. 6 (Left) is negative for the whole range of r within statistical errors and has a minimum around 0.4 fm. If the tensor potential receives a significant contribution from one-pion exchange as expected from the meson theory, V T (r) is rather sensitive to the change of the pion mass. As shown in Fig. 6 (Right), it is indeed the case: Attraction of V T (r) is substantially enhanced as the pion mass decreases.
The central and tensor potentials obtained from lattice QCD are given at discrete data points. For practical applications to nuclear physics, however, it is more convenient to parameterize the lattice results by known functions. Such a fit for V T (r) is given by the form of one-pion-exchange + one-rhoexchange with Gaussian form factors as where b 1,2,3,4 are the fitting parameters while m π (m ρ ) is taken to be the pion mass (the rho meson mass) calculated at each pion mass. The fit line for each pion mass is drawn in Fig. 6 (Right). It may be worth mentioning that the pion-nucleon coupling constant extracted from the parameter b 3 in the case of the lightest pion mass (m π = 380 MeV) gives g 2 πN /(4π) = 12.1(2.7), which is encouragingly close to the empirical value.

Convergence of the velocity expansion
The potentials are derived so far at the leading order of the velocity expansion. It is therefore important to investigate the convergence of the velocity expansion: How good is the leading order approximation ? How small are higher order contributions ? If the non-locality of the NN potentials were absent, the leading order approximation for the potentials would give exact results at all energies. The non-locality of the potentials therefore becomes manifest in the energy dependence of the potentials. So far the LO potentials are extracted with periodic boundary conditions in the spatial directions for quark fields. This leads to the lowest values of the effective center of mass energy E almost zero. To study the energy dependence, the leading order local potentials at E ≃ 45 MeV, realized by antiperiodic boundary conditions in the spatial directions, are calculated in quenched QCD at m π ≃ 529 MeV and L ≃ 4.4 fm [30,31,32]. In this case, 4 types of "momentum-wall" sources, defined by (46) are employed, where f (x) = cos((±x ± y + z)π/L). Note that f (x) = 1 corresponds to the wall source used in the periodic boundary condition. These momentum-wall sources induce L = T + 2 as well as L = A + 1 states. In Fig. 7(Left), the spin-singlet potential V C (r) (S,I)=(0,1) obtained from the L = A + 1 state at E ≃ 45 MeV (red circles) is compared with that at E ≃ 0 MeV (blue circles), while a comparison is made for the spin-triplet potentials in Fig. 8, V C (r) (S,I)=(1,0) (left) and V T (r) (right). Good agreements between results at two energies seen in these figures indicate that higher order contributions are rather small in this energy interval. In other words, these local potentials obtained at E ≃ 0 MeV can be safely used to describe the NN scattering phase shift in both spin-singlet and -triplet channels between E = 0 MeV and E = 45 MeV at this pion mass in quenched QCD.
Non-locality of the potential may become manifest also in its angular momentum dependence, since the orbital angular momentum L = r × p contains one derivative. In Fig. 7 (Right), the spin-singlet potential V C (r) (S,I)=(0,1) obtained from the L = T + 2 state, whose main component has L = 2, is compared to the one from the L = A + 1 state, whose main component has L = 0. In this case local potentials are determined at the same energy, E ≃ 45 MeV, but different orbital angular momentum. Although the statistical errors are rather large in the case of L = T 2 , a good agreement between the two is again Quenched QCD result observed, suggesting that the L dependence of the potential is small for the spin-singlet state. By these comparisons, it is observed that both energy and orbital angular momentum dependences for local potentials are very weak within statistical errors. We therefore conclude that contributions from higher order terms in the velocity expansion are small and that the LO local potentials in the expansion obtained at E ≃ 0 MeV and L = 0 are good approximations for the non-local potentials at least up to the energy E ≃ 45 MeV and orbital angular momentum L = 2.
Hereafter "potential" in this report means the local potential at the leading order, unless otherwise stated.

Full QCD results
Needless to say, it is important to repeat calculations of NN potentials in full QCD on larger volumes at lighter pion masses. The PACS-CS collaboration is performing 2 + 1 flavor QD simulations, which cover the physical pion mass [33,34]. Gauge configurations are generated with the Iwasaki gauge action   MeV, which is compared with the previous quenched results at comparable pion mass m π ≃ 731 MeV but at a ≃ 0.137 fm, given in Fig. 9(Right). Both the repulsive core at short distance and the tensor potential become significantly enhanced in full QCD. The attraction at medium distances tends to be shifted to the outer region, while its magnitude remains almost unchanged. These differences may be caused by dynamical quark effects. For a definite conclusion on this point, however, a more controlled comparison at the similar lattice spacing is needed.
In Fig. 10(Left), the spin-singlet central potential V C (r) (0,1) determined from the orbital A 1 channel is plotted at three pion masses, while the spin-triplet central potential V C (r) (1,0) and the tensor potential V T (r) from the orbital A + 1 − T + 2 couple channel are given in Fig. 11. As in the quenched QCD, the repulsive cores at short distance, the attractive pocket at medium distance and the strength of the tensor potential are all enhanced as pion mass decreases.
The phase shifts of the NN scattering for 1 S 0 obtained from the above V C (r) (0,1) are given in Fig. 10(Right). At low energy, the phase shift increases due to the attraction at medium distance, while at high energy it decreases as a consequence of the repulsive core at short distance. The shape of the scattering phase shift as a function of energy is qualitatively similar to but is much smaller than in magnitude the experimental one, plotted by the black solid line in Fig. 10(Right). As discussed before, the pion mass dependence is so large for the scattering phase shift that full QCD simulations at physical pion mass are needed to reproduce the experimental behavior.
It is noted here that ground state saturation has to be achieved to an accuracy of around 1 MeV, which is about 0.05 % of the total mass of the two-nucleon system, in order to determine the shift of the potential (E = k 2 /m N ) at the same accuracy. The value of E has a strong influence on the value of the scattering length. Such a high precision is not yet attained in full QCD calculations, since significantly larger t and, accordingly, larger statistics are required. An alternative method to overcome this difficulty will be discussed in Sec. 7 .

Hyperon Interactions
Hyperon(Y ) potentials (hyperon-nucleon and hyperon-hyperon) serve as the starting point in studying the hyper-nuclei physics. Properties of these potentials can also determine structures in the core of neutron stars. In spite of their importance, only a limited knowledge for hyperon potentials is available so far, since experimental data such as scattering phase shifts are difficult to obtain, due to the short lifetime of hyperons. Therefore it is important to calculate hyperon potentials in lattice QCD by using the potential method.

Quenched result for N Ξ potentials
Since all octet-baryons and decuplet Ω are stable in the strong interaction, there are many hyperon potentials in 2 + 1 flavor QCD. The method for the NN potentials can be straightforwardly applied to the I = 1 NΞ channel, since pΞ is simply obtained from pn by replacing a d-quark in the neutron with the s-quark and it does not have strong decay into other channels. Unstable channels such as the I = 0 NΞ, which can decay into ΛΛ via the strong interaction, will be discussed later. In addition, experimentally, not much information has been available on the NΞ interaction except for a few studies: a recent report gives the upper limit of elastic and inelastic cross sections [36] while earlier publications suggest weak attractions of Ξ− nuclear interactions [37,38,39]. The Ξ−nucleus interactions will be soon studied as one of the day-one experiments at J-PARC[40] via (K − , K + ) reaction with a nuclear target.
Ref. [41] gives the first quenched result for I = 1 NΞ potentials. Lattice parameters are the same as for the quenched NN potential. In addition to two values of the light quark mass, the quenched strange quark is introduced and is fixed to one value. The potential is calculated at (m π , m N , m Ξ ) = (511(1)MeV, 1300(4)MeV, 1419(4)MeV) and (368(1)MeV, 1167(7)MeV, 1383(6)MeV), using the NBS wave function with the interpolation operators defined by Since both p and Ξ 0 have (I, I z ) = 1/2, 1/2, the pΞ 0 system has I = 1 with the strangeness S = −2.
The left (right) of Fig. 12 gives the (effective) central potential of the pΞ 0 system obtained from the L = A 1 representation for the spin-singlet (triplet) at m π = 511 MeV and 368 MeV. Potentials in the I = 1 NΞ system for both channels show a repulsive core at r ≤ 0.5 fm surrounded by an attractive well, similar to the NN systems. In contrast to the NN case, however, the repulsive core of the pΞ 0 potential in the spin-singlet channel is substantially stronger than in the triplet channel. The attraction in the medium to long distance region( 0.6 fm ≤ r ≤ 1.2 fm ) is similar in both channels. The height of the repulsive core increases as the light quark mass decreases, while a significant difference is not seen for the attraction in the medium to long distance within statistical errors. Potentials in Fig. 12 are weakly attractive on the whole in both spin channels at both pion masses, in spite of the repulsive core at short distance, though the attraction in the triplet is a little stronger than that in the singlet.
The solid lines in Fig. 12 are the one-pion exchange potential (OPEP), given by with (m π , m N ) = (368MeV, 1167MeV), where the pseudo-vector πΞΞ coupling f πΞΞ is related to the πNN coupling as f πΞΞ = −f πN N (1 − 2α) with the parameter α = F/(F + D), and g πN N = f πN N mπ 2m N . The empirical vales, α ≃ 0.36 and g πN N /(4π) ≃ 14.0, are used for the plot. Unlike the NN potential, the OPEP in the present case has opposite sign between the spin-singlet channel and spin-triplet channel. In addition, the absolute magnitude is smaller due to the factor 1 − 2α. No clear signature of the OPEP at long distance (r ≥ 1.2 fm) is observed in Fig. 12 within statistical errors. Furthermore, there is clear departure from the OPEP at medium distance (0.6 fm ≤ r ≤ 1.2 fm) in both channels. These observations may suggest an existence of state-independent attraction.

Full and quenched QCD results for N Λ potentials
Spectroscopic studies of the Λ and Σ hypernuclei, carried out both experimentally and theoretically, suggest that the Λ-nucleus interaction is attractive while the Σ-nucleus interaction is repulsive. If this is the case, the Λ particle would be the first strange baryon instead of Σ − to appear in the core of the neutron stars [42]. It is therefore interesting and important to investigate the nature of the NΛ interaction in lattice QCD, by calculating the NΛ potential using the method of this report. Since the Λ is the lightest hyperon, the NΛ potential can be calculated as in the case of NN potentials.  Figure 13: (Left) The spin-singlet central potential for NΛ obtained from the orbital A + 1 channel in 2+1 flavor QCD at m π ≃ 414 MeV (red) and 699 MeV (green). (Right) The spin-triplet central potential and the tensor potential for NΛ obtained from the orbital A + 1 − T + 2 coupled channel in 2+1 flavor QCD at m π ≃ 414 MeV (red and blue) and 699 MeV (green and magenta). Taken from Ref. [43].  Figure 14: (Left) The spin-singlet central potential for NΛ obtained from the orbital A + 1 channel in quenched QCD at m π ≃ 407 MeV (red) and 512 MeV (green). (Right) The spin-triplet central potential and the tensor potential for NΛ obtained from the orbital A + 1 − T + 2 coupled channel in quenched QCD at m π ≃ 407 MeV (red and blue) and 512 MeV (green and magenta). Taken from Ref. [43].
In Ref. [43], the NΛ potentials are calculated in both full and quenched QCD. The 2+1 flavor full QCD gauge configurations generated by the PACS-CS collaboration are employed for the calculations of the potentials on a 32 3 × 64 lattice at a = 0.091(1) fm, while in the quenched calculation, the potentials are obtained on a 32 3 × 48 lattice at a = 0142(1) fm. Numerical values for some hadron masses for these calculations are given in Table 4, together with some lattice parameters. Fig. 13 shows the NΛ potentials obtained from 2+1 flavor QCD calculations as a function of r at m π ≃ 699 MeV and 414 MeV. The spin-singlet central potential obtained from the J = A 1 channel is plotted on the left, while the spin-triplet central potential and the tensor potential obtained from the J = T + 1 channel are given on the right. The central potential multiplied by a volume factor (r 2 V C (r)) is also shown in the left panel in addition to the V C (r) itself in the right panel, in order to compare the strength of the repulsion between two quark masses.

Flavor SU(3) limit
In order to unravel the nature of the various channels in the hyperon interactions, it is more convenient to consider an idealized flavor SU(3) symmetric world, where u, d and s quarks are all degenerate with a common finite mass. In the flavor SU(3) limit, one may capture essential features of the interaction, in particular, the short range force without contamination from the quark mass difference. Calculations in the SU(3) limit allow us to extract potentials for irreducible flavor multiplets: Potentials between asymptotic baryon states are obtained by the recombination of the multiplets with suitable Clebsh-Gordan coefficients.
In the flavor SU(3) limit, the ground state baryon belongs to the flavor-octet with spin 1/2, and two baryon states with a given angular momentum are labeled by the irreducible representation of SU(3) as where "symmetric" and "anti-symmetric" stand for the symmetry under the exchange of the flavor for two baryons. For the system with orbital S-wave, the Pauli principle imposes 27, 8 and 1 to be spin-singlet ( 1 S 0 ), while 10, 10 and 8 to be spin-triplet ( 3 S 1 ). The NBS wave function is defined by from which the corresponding (effective) central potential is given as   Table 5. For example, symmetric NN belongs to 27, which therefore can be considered as the NN spin-singlet potential in the flavor SU(3) symmetric limit. Similarly V (10) , V (10) and V (8a) can be considered as some BB potentials of the particle basis in the SU(3) symmetric limit, while V (1) and V (8s) are always mixtures of different BB potentials in the particle basis. Fig. 15 shows V (27) (r) and V (10) (r), which correspond to spin-singlet and spin-triplet NN potentials, respectively. Both have a repulsive core at short distance with an attractive pocket around 0.6 fm. These qualtative features are consistent with the previous results found for the NN potential in both quenched and full QCD. The upper-right panel of Fig. 16 shows that V (10) (r) has a stronger repulsive core and  Fig. 16 has a very strong repulsive core among all 6 channels, while V (8a) (r) in the lower-right panel has a very weak repulsive core. In contrast to all other cases, V (1) (r) shows attraction instead of repulsion at all distances, as shown in the lower-left panel.
Above features are consistent with what has been observed in a phenomenological quark model [46]. In particular, the potential in the 8 s channel in this quark model becomes strongly repulsive at short distance since the six quarks cannot occupy the same orbital state due to the Pauli exclusion for quarks. On the other hand, the potential in the 1 channel does not suffer from the quark Pauli exclusion and can become attractive due to the short-range gluon exchange. Such agreements between the lattice data and the phenomenological model suggest that the quark Pauli exclusion plays an essential role for the repulsive core in BB systems.
The BB potentials in the baryon basis can be obtained from those in the SU(3) basis by a unitary rotation as where U is a unitary matrix which rotates the flavor basis |X to the baryon basis |i , i.e. |i = U iX |X . The explicit forms of the unitary matrix U in terms of the CG coefficients are found in Appendix B.
In Fig. 17, as characteristic examples, let us show the spin-singlet potentials for S=−2, I=0 channel determined from the orbital A + 1 representation at m π = 835 MeV. To obtain V ij (r), the potentials in Then the right hand side of Eq. (52) is used to obtain the potentials in the baryon basis. The left panel of Fig. 17 shows the diagonal part of the potentials. The strong repulsion in the 8 s channel is reflected most in the ΣΣ(I=0) potential due to its largest CG coefficient among three channels. The strong attraction in the 1 channel is reflected most in the NΞ(I=0) potential due to its largest CG coefficient. Nevertheless, all three diagonal potentials have a repulsive core originating from the 8 s component. The right panel of Fig. 17 shows the off-diagonal parts of the potentials which are comparable in magnitude to the diagonal ones. Since the off-diagonal parts are not negligible in the baryon basis, fully coupled channel analysis is necessary to study observables. A similar situation holds even in (2+1)-flavors where the strange quark is heavier than up and down quarks: The SU(3) basis with approximately diagonal potentials is useful for obtaining essential features of the BB interactions, while the baryon basis with substantial magnitude of the off-diagonal potentials is necessary for practical applications. Other potentials in the baryon basis are given in Ref. [45]. Since the 8 s state does not couple to the spin-triplet channel, the repulsive cores in the spin-triplet channel are relatively small. The off-diagonal potentials are not generally small: For example, the NΛ-NΣ potential in the spin-triplet channel is comparable in magnitude at short distances with the diagonal NΛ-NΛ and NΣ-NΣ potentials. Although all quark masses of 3 flavors are degenerate and rather heavy in these simulations, the coupled channel potentials in the baryon basis may give useful hints for the behavior of hyperons (Λ, Σ and Ξ) in hyper-nuclei and in neutron stars [47,48].
The flavor singlet channel has attraction for all distances, which might produce the bound state, the H-dibaryon, in this channel. The present data, however, are not sufficient to make a definite conclusion on the H-dibaryon, since a single lattice with small extension L ≃ 2 fm is employed. In order to investigate whether the H-dibaryon exists or not in the flavor SU(3) limit, data on several different volumes are needed. Such a study on the H-dibaryon will be discussed in Sec. 7.
In order to extend the study in the flavor SU(3) limit to the real world where the strange quark is much heavier than light quarks, the potential method used so far has to be extended to more general cases, which will also be considered in Sec. 7.

Origin of repulsive core
As seen in the previous sections, lattice QCD calculations show that the NN potential defined through the NBS wave function has not only the attraction at medium to long distance that has long been well understood in terms of pion and other heavier meson exchanges, but also a characteristic repulsive core at short distance, whose origin is still theoretically unclear. Furthermore, the BB potentials in the flavor SU(3) limit show several different behaviors at short distance: some have a stronger/weaker repulsive core than NN while the singlet has an attractive core. In this section, recent attempts [49,50,51] to theoretically understand the short distance behavior of the potential in terms of the operator product expansion (OPE) is explained.

OPE and repulsive core
Let us first explain the basic idea. We consider the equal time NBS wave function defined by where |E is some eigenstate of a certain system with total energy E, and O A , O B are some operators of this system. (We suppress other quantum numbers of the state |E for simplicity.) The OPE reads Suppose that the coefficient function of the OPE behaves in the small r(= |r|) limit as where θ, φ are the angles of r, the NBS wave function becomes where The potential at short distances can be calculated from this expression. For example, in the case of the Ising field theory in two dimensions, the OPE for the spin field σ is given by where O 1 (x) (=:ψψ(x) : in terms of free fermion fields) is an operator of dimension one. This leads to where |E is a two-particle state with energy E = 2 √ k 2 + m 2 . From this expression the potential becomes in the r → 0 limit. The OPE predicts not only the r −2 behavior of the potential at short distance but also its coefficient −3/16. Furthermore the potential at short distance does not depend on the energy of the state in this example [30,52].
In QCD the dominant terms at short distance have α C = 0. Among these terms, we assume that C has the largest contribution such that β C > β C ′ for ∀C ′ = C. Since such dominant operators with α C = 0 mainly couple to the zero angular momentum (L = 0) state, let us consider the NBS wave function with L = 0. Applying ∇ 2 to this wave function, we obtain the following classification of the short distance behavior of the potential.
(1) β C = 0: The potential at short distance is energy independent and becomes which is attractive for β C > 0 and repulsive for β C < 0.
(2) β C = 0: In this case the potential becomes where β C ′ < 0 is the second largest exponent. The sign of the potential at short distance depends on the sign of D C ′ (E)/D C (E).
On the lattice, we do not expect divergence at r = 0 due to lattice artifacts at short distance. The above classification holds at a ≪ r ≪ 1/Λ QCD , while the potential becomes finite even at r = 0 on the lattice.
Since QCD is an asymptotic free theory, the 1-loop calculation for anomalous dimensions becomes exact at short distance. The OPE in QCD is written as where g (m) is the renormalized coupling constant (quark mass) at scale µ. In the limit that r = |y| = e −t R → 0 ( t → ∞ with fixed R), the renormalization group analysis leads to where β (1) = 1 16π 2 11 − 2N f 3 is the QCD beta-function at 1-loop, and γ C,(1) Here γ (1) X is the 1-loop anomalous dimension of the operator O X . An appearance of D C AB (R, 0, 0, µ) on the right-hand side tells us that it is enough to know the OPE only at tree level. From the above expression, β C is given by β C = γ C,(1) AB 2β (1) .

Two flavor case
We first consider the OPE for N f = 2 QCD [49]. The 1-loop calculations show that the largest value of β C is always zero for both spin-singlet and spin-triplet channels for N f = 2 QCD and that the second largest value of β C ′ is given by where S = 0, 1 denotes the total spin. This corresponds to the case (2) in the previous subsection. Therefore the OPE and renormalization group analysis in QCD predicts the universal functional form of the NN central potential at short distance as X 27 8 s 1 10 10 8 a γ (X) 0 6 12 0 0 4 Non-relativistic op. yes no yes yes yes yes Table 6: The largest value of β C in unite of 1/(33 − 2N f ) of 3-flavor QCD for each representation. The last line indicates that the operator corresponding to the largest value of β C exists or not in the non-relativistic limit.
which is a little weaker than a 1/r 2 singularity, while for the tensor potential we have at the 1-loop order. The OPE, however, can not tell whether the potential at short distance is repulsive or attractive, which is determined by the sign of the coefficient. If D X (E) and D Y (E) are evaluated by the nonrelativistic quark model at the leading order, we obtain For both cases, the ratio has positive sign, which gives repulsion at short distance, the repulsive core.

Extension to three flavors
The above calculation has been extended to N f = 3 QCD [51]. In the 3-flavor case, some channel may become attractive at short distance since the Pauli exclusion principle is less significant than in the 2flavor case. Indeed the lattice QCD calculations in the flavor SU(3) limit shows the attractive potential for the singlet channel, as seen in the previous section. The largest value of β C is given for each X representation of the flavor SU(3) in unit of 1/(33 − 2N f ) in the table 6, where we define From the table, we observe that the largest value of β C is zero in the 27, 10 and 10 channels. This is consistent with the nucleon case in the previous subsection, which belong to 27 (spin-singlet) and 10(spin-triplet). These three channels correspond to the case (2), so that the potentials are given by eq. (63). On the other hand, the largest value of β C becomes positive in the 8 s , 8 a and 1 channels, which correspond to the case (1). Therefore the (effective) central potential becomes attractive at short distance as where m B is the octet baryon mass. The attractive core of the potential in the flavor singlet channel agrees with the behavior of the potential found for the numerical simulation of lattice QCD in the previous section, while for other two channels, 8 s and 8 a , the prediction by the OPE disagrees with the lattice QCD results: The potential in the 8 s channel is most repulsive among 6 channels and the potential in the 8 a channel still has a repulsive core, which is however weaker than others. The disagreement between the OPE and the lattice QCD result for the 8 s channel may be understood by the fact that no local 6 quark operator exists for this channel in the non-relativistic limit, as shown in the table 6: the 8 s operator with the largest positive β C has a very small coefficient at low energy, so that the other operators with zero or negative β C may still dominate at distance scales comparable to the lattice spacing a = 0.1 − 0.2 fm. For the 8 a case, the weakest repulsive core in the lattice QCD simulation suggests that the attraction from the leading operator with the positive β C may be cancelled by other contributions from the sub-leading operators with zero or negative β C at the distance scale comparable to the lattice spacing of the simulations. It is therefore important to confirm the prediction from the OPE, by investigating the behavior of the repulsive core for each channel in the flavor SU(3) limit at finer lattice spacings and hopefully in the continuum limit.

Extensions
In this section, recent extensions of the potential method are considered.

Inelastic scattering
The potential method discussed so far is shown to be quite successful in order to describe elastic hadron interactions. Hadron interactions in general, however, lead to inelastic scatterings as the total energy of the system increases. In order to extract hadron interactions which describe such inelastic scatterings from lattice QCD, an extension of the potential method is considered in this subsection.
Let us first discuss the case of A + B → C + D scattering where A, B, C, D represent some 1-particle states. This is a simplified version of the octet baryon scattering in the strangeness S = −2 and isospin I = 0 channel, where ΛΛ, NΞ and ΣΣ appear as asymptotic states of the strong interaction if the total energy is larger than 2m Σ .
We here assume k is the total energy of the system, and E X k = m 2 X + k 2 . In this situation, the QCD eigen-state with the quantum numbers of the AB state and center of mass energy W is expressed in general as where We define the following NBS wave functions, Using the partial wave decomposition such that 4 the NBS wave function of the 2-channel system behaves for large r as where δ i l (W ) is the scattering phase shift , whereas θ(W ) is the mixing angle. This expression shows that the NBS wave functions for large r agree with scattering waves described by two scattering phases δ i l (W ) (i = 1, 2) and one mixing angle θ(W ). Because of this property, these wave functions satisfy (∇ 2 + k 2 )ϕ AB (r, k) = 0, (∇ 2 + q 2 )ϕ CD (r, q) = 0 (80) for r → ∞.
Let us now consider QCD in the finite volume V . In the finite volume, |AB, W and |CD, W are no longer eigen-states of the hamiltonian. True eigenvalues are shifted from W to W i = W + O(V −1 ) (i = 1, 2). By diagonalization method in lattice QCD simulations, it is relatively easy to determine W 1 and W 2 . With these values Lüscher's finite volume formula gives two conditions, which, however, are insufficient to determine three observables, δ 1 l , δ 2 l and θ. (See [53,54,55] for recent proposals to overcome this difficulty.) An alternative approach to extract three observables, δ 1 l , δ 2 l and θ, has been proposed in lattice QCD through the above NBS wave functions [56,57]. We consider the NBS wave functions at two different values of energy, W 1 and W 2 , in the finite volume: We then define the coupled channel non-local potentials from the coupled channel Schrödinger equation as for i = 1, 2. As before the velocity expansion is introduced as and at the leading order of the expansion, we have These equations for i = 1, 2 can be solved as Once we obtain the coupled channel local potentials V XY,V Z (x), we solve the coupled channel Schrödinger equation in infinite volume with some appropriate boundary condition such that the incoming wave has a definite ℓ and consists of the AB state only , in order to extract three observables for each ℓ (δ 1 l (W ), δ 2 l (W ) and θ(W )) at all values of W . Of course, since V XY,V Z is the leading order approximation in the velocity expansion of U XY,V Z (x; y), results for three observables δ 1 l (W ), δ 2 l (W ) and θ(W ) at W = W 1 , W 2 are also approximate ones and might be different from the exact values. By performing an additional extraction of V XY,V Z (x) at (W 3 , W 4 ) = (W 1 , W 2 ), we can test how good the leading order approximation is.
The method considered above can be generalized to inelastic scattering where a number of particles is not conserved. For illustration, let us consider the scattering A + B → A + B and A + B → A + B + C where the total energy W satisfies m A + m B + m C < W < m A + m B + 2m C .
The following NBS wave functions at the center of mass system are used: where and 1/µ AB = 1/m B + 1/m C . Here y = r B − r C is a relative coordinate between B and C with the reduced mass µ BC , while x = r A − R BC is the one between A and the center of mass of B and C with R BC = (m B r B + m C r C )/(m B + m C ). We define the non-local potential from the coupled channel equations as where with another reduced mass defined by 1/µ A,BC = 1/m A + 1/(m B + m C ). We consider the following velocity expansions  where the hermiticity of the non-local potentials gives V AB,ABC (x, y) = V ABC,AB (x, y).
At the leading order of the velocity expansions, the coupled channel equations become By considering two values of energy such that W = W 1 , W 2 , we can determine V ABC,AB and V ABC,ABC from the second equation as Using the hermiticity relation V AB,ABC (x, y) = V ABC,AB (x, y), we can extract V AB,AB from the first equation as for W = W 1 , W 2 . A difference of V AB,AB (x) between two estimates at W 1 and W 2 gives an estimate for higher order contributions in the velocity expansions.
Once we obtain V AB,AB , V AB,ABC = V ABC,AB and V ABC,ABC , we can solve the coupled channel Schrödinger equations in the infinite volume, in order to extract physical observables. As W increases and becomes larger than m A +m B +nm C , the inelastic scattering A+B → A+B +nC becomes possible. As in the case of A+B → A+B +C in the above, we can define the coupled channel potentials including this channel, though calculations of the NBS wave functions for multi-hadron operators become more and more difficult in practice.

Coupled channels with S = −2 and I = 0
As an application of the method in the previous subsection, let us consider BB potentials for the S = −2 and I = 0 channel, which consist of the ΛΛ, NΞ and ΣΣ components in terms of low-lying octet baryons. Mass differences of these components are quite small such that 2m Λ = 2232 Mev, m N + m Σ = 2257 MeV and 2m Σ = 2386 MeV. Using diagonalized source operators, the NBS wave functions at three different values of energy, for i = 0, 1, 2, are extracted in lattice QCD simulations, where AB = ΛΛ, NΞ and ΣΣ, and k i AB satisfies Using the notation where 1/µ AB = 1/m A + 1/m B , the coupled channel 3 × 3 potential matrix is given by Here the last factor is the inverse of the 3 × 3 matrix ϕ W i CD (r, k i CD ) with indices i and CD. Gauge configurations generated on a 16 3 × 32 lattice at a ≃ 0.12 fm ( therefore L ≃ 1.9 fm) in 2+1flavor full QCD simulation are employed to calculate the coupled channel potentials at three different values of the light quark mass with the fixed bare strange quark mass [58]. Quark propagators are calculated with the spatial wall source at t 0 with the Dirichlet boundary condition in time at t = t 0 +16. The wall source is placed at 16 different time slices on each gauge configuration, in order to enhance signals, together with the average over forward and backward propagations in time. Corresponding hadron masses and number of gauge configurations are given in table 7.
The coupled channel potential matrix V AB,CD from the NBS wave function for Set 1 is shown in Figure 18. The flavor dependence of the height of the repulsive core at short distance region is observed. In particular, the ΣΣ potential has the strongest repulsive core of these three channels. It is interesting to see that off-diagonal parts of the potential matrix roughly satisfy the hermiticity relation V AB,CD = V CD,AB within statistical errors. In addition the off-diagonal parts are similar in magnitude for V ΛΛ,ΣΣ and V N Ξ,ΣΣ with the diagonal parts, but V ΛΛ,N Ξ is much smaller.  Figure 19: Transition potentials in the flavor SU(3) IR basis. Red, blue and green symbols correspond to results of Set1, Set2 and Set3, respectively. The result of the flavor SU(3) symmetric limit at the same strange quark mass is also plotted with brown symbols [45]. Taken from Ref. [58].
In order to compare the results of the potential matrix calculated in three configuration sets, the potentials from the particle basis are transformed to those in the flavor SU(3) irreducible representation (IR) basis as where U is a unitary transformation matrix whose explicit form is given in Appendix B. The potential matrix in the IR basis is convenient and a good measure of SU(3) breaking effects by comparing three configuration sets since it is diagonal in the SU(3) symmetric limit.
In Figure 19, the results of the potential matrix in the IR basis are compared between different configuration sets, together with the one in the flavor SU(3) symmetric limit. As the pion mass decreases, the repulsive core in the V 27,27 potential increases. The V 1,27 and V 8,27 transition potentials are consistent with zero within statistical errors. On the other hand, it is noteworthy that the flavor SU(3) symmetry breaking effect becomes manifest in the V 1,8 transition potential.

Time dependent method
One of the practical difficulties to extract the NBS wave function and the potential from the correlation function Eq.(24) is to achieve the ground state saturation in numerical simulations at large but finite t − t 0 with reasonably small statistical errors. While the stability of the potential against t − t 0 has been confirmed within statistical errors in numerical simulations reviewed in the report, the determination of W for the ground state suffers from systematic errors due to contaminations of possible excited states, which can be seen as follows. There exist 3 different methods to determine W . The most well-known method is to determine W from the t − t 0 dependence of the correlation function eq.(24) summed over r to pick up the zero momentum state. On the other hand, one can determine k 2 of W by fitting the r dependence of the NBS wave function with its expected asymptotic behavior at large r or by reading off the constant shift of the Laplacian part of the potential from zero at large r. Although the latter two methods usually give consistent results within statistical errors, the first method (the t dependence method) sometimes leads to a result different from those determined by the latter two (called the r dependent method together) at the value of t−t 0 usually employed in numerical simulations. Although, in principle, the increase of t − t 0 is needed in order to see an agreement between t and r dependence methods, it is difficult in practice due to larger statistical errors at larger t − t 0 for the two-baryon system.
In order to overcome this practical difficulty, the method to extract the potential from the NBS wave function has been modified as follows. Let us consider the correlation function Eq.(24) again: If t is large enough so that contributions from O(e −W th t ) terms can be neglected 5 , we have where where the velocity expansion is introduced in the last line and higher other than the leading order terms are then omitted. The leading order potential is therefore given by or where R(r, t) = F (r, t)/e −2m N t . Here it is assumed that O(e −W th t ) contributions can be neglected at large t. The non-relativistic formula for V LO (r) above can be easily generalized to the case that masses of two particles are different, by the replacement that R(r, t) = F (r, t)/e −(m A +m B )t . Note also that the potential extracted in this method automatically satisfies that V LO (r) → 0 as r → 0 without the constant shift. This property may be used to check whether this extraction works correctly or not. On the lattice, the t derivative should be approximated by the t difference. In practice, one may adopt a particular method for the t difference, in order to reduce statistical as well as systematic errors for V LO (r).
The non-relativistic approximation can be removed by using the second order derivative in t as as long as O(e −W th t ) contributions are negligibly small. For this method to apply, two particles should have the same mass. Statistical errors of the second order difference on the lattice must be kept small in numerical simulations.
One may introduce a more general correlation function as Using this new quantity, we have from which the non-local potential is extracted as HereF −1 (x, y, t) is the approximated inverse of the hermitian operator F (x, y, t), defined bỹ where λ n (t) and v(x, t) are an eigenvalue and a corresponding eigenvector of F (x, y, t), respectively, and zero eigenvalues are removed in the summation. Since the modified potential also satisfies the same Schrödinger equation ∀{c n }, the non-local potential is NOT unique, and U(x, y) is scheme dependent, as discussed before.

Bound H dibaryon in flavor SU(3) limit
As an application of the method in the previous subsection, let us consider the singlet potential in the flavor SU(3) limit in order to investigate whether the bound H dibaryon exists or not in this case. At the leading order of the velocity expansion, the central potential is defined in this method by To check the qualitative consistency with previous results, the central potential in the 27-plet channel is plotted in Fig. 20 obtained in three different lattice volumes with L = 1.94, 2.90, 3.87 fm at m ps = 1015 MeV and (t − t 0 )/a = 10. This is the case corresponding to the spin-singlet NN potential. Compared with statistical errors, the L dependence is found to be negligible. The t dependence is also small as long as (t − t 0 )/a ≥ 9. As expected, the potential approaches zero automatically for large r. The figure shows a repulsive core at short distance surrounded by an attractive well at medium and long distances, which is qualitatively consistent with previous results in quenched and full QCD simulations. Shown in Fig. 21(Left) and Fig. 21(Right) are the volume and the quark mass dependences of the central potential in the flavor-singlet channel V (1) C (r), respectively, at (t−t 0 )/a = 10 where the potentials do not have appreciable change with respect to the choice of t. The flavor-singlet potential is shown to have an "attractive core" and to be well localized in space. Because of the latter property, no significant volume dependence of the potential is observed within the statistical errors, as seen in Fig. 21(Left). As the quark mass decreases in Fig. 21(Right) , the long range part of the attraction tends to increase.
The resultant potential is fitted by the following analytic function composed of an attractive Gaussian core plus a long range (Yukawa) 2 attraction: With the five parameters, b 1,2,3,4,5 , the lattice results can be fitted reasonably well with χ 2 /dof ≃ 1. The fitted result for L = 3.87 fm is shown by the dashed line in Fig. 21(Left). Solving the Schrödinger equation with the fitted potential in infinite volume, the energies and the wave functions are obtained at the present quark masses in the flavor SU(3) limit. It turns out that, at each quark mass, there is only one bound state with binding energy 30-40 MeV. In Fig. 22(Left), the energy and the root-mean-squared (rms) distance of the bound state are plotted in the case of (t − t 0 )/a = 9, 10, 11 at m ps = 673 MeV and L = 3.87 fm, where errors are estimated by the jackknife method. Although the statistical error increases as t increases, we observe small changes of central values, which are considered as the systematic errors. Fig. 22(Right) shows the energy and the rms distance of the bound state at each quark mass obtained from the potential at L = 3.87 fm and (t − t 0 )/a = 10. Despite the fact that the potential has quark mass dependence, the resultant binding energies of the H-dibaryon are insensitive in the present range of the quark masses. This is due to the fact that the increase of the attraction toward the lighter quark mass is partially compensated by the increase of the kinetic energy for the lighter baryon mass. It is noted that there appears no bound state for the potential of the 27-plet channel in the present range of the quark masses.
The final results of the binding energy B H and the rms distance r 2 are summarized below, where the 1st and 2nd parentheses correspond to statistical errors and systematic errors from the t-dependence, respectively. Recently the existence of H-dibaryon is also suggested by a direct calculation of its binding energy in 2+1 full QCD simulations [60], where B H = 16.6(2.1)(4.6) MeV is reported in the L → ∞ extrapolation at m π ≃ 389 MeV .

Other applications
In the last section, some other applications of the potential method are reviewed.

Three nucleon force
Recent precise calculations of few-nucleon systems clearly indicate that the 2 nucleon force alone is insufficient to understand the nuclei, which calls for three (and/or more) nucleon forces. Actually, the three nucleon force (TNF) is supposed to play an important and nontrivial role in various phenomena in nuclear and astrophysics. For the binding energies of light nuclei, the attractive TNF is required to reproduce the experimental data. On the other hand, the repulsive TNF is necessary to reproduce the empirical saturation density of symmetric nuclear matter. For the equation of state(EoS) of asymmetric nuclear matter, a repulsive TNF is required to explain the observed maximum neutron star mass.
Pioneered by Fujita-Miyazawa [61], the TNF has been mainly studied from the two-pion exchange picture with the ∆-excitation. In addition, the repulsive TNF is often introduced phenomenologically [62]. Recently, the TNF based on chiral effective field theory is developing [63], but the unknown low-energy constants can be obtained only by the fitting to the experimental data. Since the TNF originates from the fact that the nucleon is not a fundamental particle, it is essential to study the TNF from fundamental degrees of freedom(DoF) , i.e., quarks and gluons.
In Ref. [64,65] such first-principle calculations of the TNF in lattice QCD have been reported. If the potential method is applied to the three nucleon (3N) system, a straightforward calculation is impossible due to the significantly enlarged DoF. In Ref. [64,65], two different approaches have been considered.
The NBS wave function for 3N is defined by where r 12 ≡ x 1 − x 2 , r 123 ≡ x 3 − (x 1 + x 2 )/2 are the Jacobi coordinates, and |W is the 3N state with energy W . At the leading order of the velocity expansion, the NBS wave function satisfies denotes the potential between (i, j)-pair, V TNF (r 12 , r 123 ) the TNF, µ 12 = m N /2, µ 123 = 2m N /3 the reduced masses. If ϕ W (r 12 , r 123 ) is calculated for all r 12 , r 123 and all V 2N,ij (x i − x j ) are available by lattice calculations, V TNF (r 12 , r 123 ) can be extracted. Unfortunately, this is not the case: Since both r 12 and r 123 have L 3 DoF, the calculation cost is more expensive by a factor of L 3 compared to the 2N system. Furthermore, the number of diagrams to be calculated in the Wick contraction tends to diverge with a factor of N u ! × N d ! (N u,d are numbers of u,d quarks in the system). It is also noted that not all 2N potentials are available in lattice QCD at this moment: Only parity-even 2N potentials have been obtained so far.
The first method in Ref. [64] to avoid these problems is to consider the effective 2N potential in the 3N system by taking the summation over the location of the spectator nucleon N(x 3 ), ϕ W (r 12 ) = x 3 ϕ W (r 12 , r 123 ) = r 123 ϕ W (r 12 , r 123 ).
The effective potential between N(x 1 ) and N(x 2 ) is then defined by T potential. Taken from Ref. [64].
In this calculation, the DoF of r 123 is integrated out beforehand, and thus the calculation cost is reduced by a factor of ∼ 1/L 3 , compared to the straightforward calculation. The difference V eff ( r) − V 2N ( r) can be considered to be the "finite density effect" in the 3N system. Part of this effect is attributed to the genuine 2N potential with the nontrivial 3N correlation, while another originates from the genuine TNF.
The triton channel( I = 1/2, J P = 1/2 + ) is studied as the 3N system. Since the spectator nucleon is projected to the S-wave, the possible quantum numbers between the (effective) 2N are only 2S+1 L J = 1 S 0 , 3 S 1 , 3 D 1 . Gauge configurations in 2-flavor QCD on a 16 3 × 32 lattice at a ≃ 0.16 fm [66] are employed for the calculation at m π ≃ 1.13 GeV and m N ≃ 2.15 GeV. Fig. 23(Left) show results for V eff (r) in the triton channel at t − t 0 = 8, where the constant shift by energy is not included for the central potentials. It is noteworthy that V eff (r) is obtained with good precision even though the signal to noise ratio is expected to be worse for more quarks in the system. Fig. 23(Right) gives V eff (r) − V 2N (r) for the tensor potential, which is free from the constant shift. The difference is consistent with zero within a few MeV statistical errors. Similar results are reported for central potentials as well in Ref. [64]. There is no indication of a TNF effect. A possible explanation is that the TNF effect is suppressed at heavy quark mass. Basically similar results are obtained, however, for lighter pion masses(m π ≃ 0.7 GeV and 0.57 GeV) [64]. Another possibility is that the TNF effect is suppressed by the summation over the location of the spectator nucleon. While the TNF effect is expected to be enhanced when all three nucleons are close to each other, such 3-dimensional spacial configurations have small contributions in the spectator summations.
In order to assess this possibility, a second method has been investigated in Ref. [64,65], where the linear setup with r 123 = 0 is used for the 3N wave function. In this case, the third nucleon is attached to (1, 2)-nucleon pair with only S-wave. Considering the total 3N quantum numbers of I = 1/2, J P = 1/2 + , the wave function can be completely spanned by only three bases, which can be labeled by the quantum numbers of (1, 2)-pair as 1 S 0 , 3 S 1 , 3 D 1 . Therefore, the Schrödinger equation is simplified to the 3 × 3 coupled channel equations with the bases of ϕ1 S 0 , ϕ3 S 1 , ϕ3 D 1 . Even in this case the subtraction of V 2N remains nontrivial: the parity-odd potentials, which must be subtracted, are not available in lattice QCD at this moment. The subtraction problem of parity-odd potentials can be avoided for triton by using the symmetric wave function, Combined with the Pauli principle, it is automatically guaranteed that any 2N-pairs couple with even parity only, since this wave function is anti-symmetric in spin/isospin spaces for any 2N-pairs. Therefore the TNF can be extracted unambiguously in this channel, without the information of parity-odd 2N potentials.
The same gauge configurations used for the effective 2N potential study are employed in the numerical simulations. Fig. 24(Left) gives each wave function of ϕ S = 1 , ψ3 D 1 as a function of r = |r 12 /2| in the triton channel at t − t 0 = 8. Among these three ϕ S dominates the wave function, since ϕ S contains the component for which all three nucleons are in S-wave.
By subtracting the V 2N from the total potentials in the 3N system, the TNF is detemined. Fig. 24 (Right) shows results for the scalar/isoscalar TNF, where the r-independent shift of V 2N is not included, and thus about O(10) MeV systematic error is understood. There are various physical implications in Fig. 24 (Right). At the long distance region of r, the TNF is small as is expected. At the short distance region, the indication of a repulsive TNF is observed. Recalling that the repulsive short-range TNF is phenomenologically required to explain the saturation density of nuclear matter, etc., this is a very encouraging result. Of course, further study is necessary to confirm this result, e.g., the study of the ground state saturation, the evaluation of the constant shift by energies, the examination of the discretization error.

Meson-baryon interactions
The potential method can be naturally extended to the meson-baryon systems and the meson-meson systems. In this subsection, two applications of the potential method to the meson-baryon system are discussed.
The first application is the study of the KN interaction in the I(J P ) = 0(1/2 − ) and 1(1/2 − ) channels by the potential method. These channels may be relevant for the possible exotic state Θ + , whose existence is still controversial.
The KN potentials in isospin I = 0 and I = 1 channels have been calculated in 2 + 1 full QCD simulations, employing 700 gauge configurations on a 16 3 × 32 lattice at a = 0.121(1) fm and (m π , m K , m N ) = (871(1), 912(2), 1796 (7)) in unit of MeV [67]. Fig. 25   channels. The large r behavior of the NBS wave functions in both channels do not show a sign of a bound state, though more detailed analysis is needed with larger volumes for a definite conclusion. On the other hand, the small r behavior of the NBS wave functions suggests a repulsive interaction at short distance (r < 0.3 fm). The repulsion in the I = 1 channel seems to be stronger than that in the I = 0 channel. The LO potential V (r) for the KN state without the constant energy shift E is shown in Fig. 26 in the I = 0 (left) and I = 1 (right). As expected from the NBS wave functions in Fig. 25, the repulsive interactions are observed at short distance in both channels, while the attractive well appears at the medium distance (0.4 < r < 0.8 fm) in the I = 0 channel. These results indicate that there are no bound states in I(J π ) = 0(1/2 − ) and 1(1/2 − ) states at m π ≃ 870 MeV.
The second application is to study the charmonium-nucleon interaction using the potential method. Since charmonium does not share the same quark flavor with the nucleon, the charmonium-nucleon interaction is mainly induced by the genuine QCD effect of multi-gluon exchanges. Theoretical studies based on QCD suggest that the cc-N interaction is weakly attractive. It is argued that the cc-nucleus (A) bound system may be realized for the mass number A ≥ 3 if the attraction between the charmonium and the nucleon is sufficiently strong [68,69]. Precise information on the cc-N potential V ccN (r) is  Figure 27: (Left) The effective central potential in the s-wave η c -N system at m π = 640 MeV. The solid line is the fit by the Yukawa potential while the dashed line is the one by the phenomenological potential adopted in Ref. [68]. (Right) The volume dependence of the η c -N potential. Taken from Ref. [71].
therefore indispensable for exploring nuclear-bound charmonium states such as the η c -3 He or J/ψ-3 He bound state in few body calculations [70]. In Ref. [72], the charmonium-nucleon potentials are calculated in quenched QCD on 16 3 × 48 and 32 3 × 48 lattices at a ≃ 0.94 fm at three different values of the light quark mass corresponding to (m π , m N ) ≃ (640, 1430), (720, 1520), (870, 1700) in unit of MeV and one fixed value of the charm quark mass corresponding to m ηc ≃ 2920 MeV and m J/Ψ ≃ 3000 MeV.
The effective central η c -N potential, evaluated from the NBS wave function with measured values of E and µ, is shown in Fig. 27(Left). The η c -N potential clearly exhibits an entirely attractive interaction between charmonium and the nucleon without any repulsion at all distances. Absence of the short range repulsion (the repulsive core) may be related to absence of the Pauli exclusion between the heavy quarkonium and the light hadron. The interaction is exponentially screened in the long distance region r ≥ 1 fm. The fit of the potential with the Yukawa form −γe −αr /r gives γ ∼ 0.1 and α ∼ 0.6 GeV (solid line in Fig. 27), which are compared with the phenomenological values γ = 0.6 and α = 0.6 GeV in Ref. [68] (dashed line). The strength of the Yukawa potential γ is six times smaller than the phenomenological value, while the Yukawa screening parameter α obtained from lattice QCD is comparable with the phenomenological one. The cc-N potential is found to be rather weak in lattice QCD.
As shown in Fig. 27(Right), there is no significant difference between potentials at two different spatial sizes (La ≈ 3.0 and 1.5 fm). The size dependence of the η c -N potential seems small since the η c -N potential is quickly screened to zero and turns out to be short ranged.
No large quark-mass dependence is observed in Fig. 28(Left). This may be understood by the argument that the cc-N interaction is mainly governed by multi-gluon exchanges, which do not depend explicitly on the quark mass. Taking a closer look at the inset of Fig. 28(Left), however, one finds that the attractive interaction in the η c -N system tends to get slightly weaker as the light quark mass decreases.
The spin-independent part of the central J/ψ-N potential is shown at m π = 640 MeV in Fig. 28(Right), together with the η c -N potential for comparison. While no qualitative difference between the η c -N and J/ψ-N potentials is observed, the attractive interaction in the J/ψ-N potential is a little stronger than the η c -N potential. The attraction, however, is still not strong enough to form a bound state in the J/ψ-N channel.

Two-color QCD and potentials
In Ref. [72], potentials between baryons in two-color (SU (2)) QCD are investigated on the lattice. In two-color QCD, two quarks (diquark) form a "baryon", whose local interpolating operator is given by where ε ab is the 2 × 2 anti-symmetric tensor, i, j are flavor indices and Γ = C, Cγ 5 , Cγ µ , Cγ µ γ 5 with the charge conjugation matrix C. The NBS wave function is then defined by ϕ W,Γ ij,kl (r)e −W t = 0|T {D ij,Γ (x + r, t)D kl,Γ (x, t)}|W where |W is an eigenstate of two color QCD with total energy W and "baryon" number 2. The potential derived from this wave function is denoted as V Γ ij,kl (r). The potentials have been calculated in SU(2) quenched QCD on a 24 3 × 48 lattice at a ≃ 0.1 fm determined from the string tension with the assumption that √ σ = 440 MeV. The lightest "baryon" corresponds to the scalar diquark state (Γ = C). Therefore the NBS wave function is evaluated for Γ = C at four values of the quark mass, which give m S = 1.044(2), 0.836(2), 0.618(2), 0.377(3) in lattice unit with m S being the scalar diquark mass. The potential is then extracted from the A 1 state with Γ = C, and is therefore denoted by V ij,kl (r). Fig. 29 shows the LO central potentials plotted as functions of r for (i, j, k, l) = (1, 2, 3, 4)(left) and (1, 2, 1, 2)(right), where quark-exchange diagrams are absent for the former while they exist for the latter .
The potential V 12,34 (r) has an attractive interaction at all length scales, which becomes stronger as m S decreases. The r-dependence is monotonic at all values of m S . At large r (r ≥ 4 in lattice units), the potential has small m S dependence, while it has stronger m S dependence at small r (r ≤ 4). The magnitude of the potential decreases as m S increases and finally the potentials at m S = 0.836 and 1.044 coincide with each other.
For V 12,12 , on the other hand, strong repulsions appear at short distance. The m S dependence is not monotonic: The potential is a smooth function of r at small m S , while it has a pocket at intermediate distance for m S = 1.044. As m S decreases, the repulsive core rapidly increases and the attractive pocket disappears.
The qualitative difference between the two cases suggests that contributions from quark exchange diagrams, which more or less represent the Pauli exclusion effect, are responsible for the repulsive core. On the other hand, the attraction may be explained by gluon exchanges, which are the main part of the potential V 12,34 (r). More details analysis can be found in Ref. [72].

Conclusion
In this report, we have reviewed the progress on the study of hadron interactions via the potential method in lattice QCD. The key quantity of the method is the NBS wave function of the two particle state. Not only is the asymptotic behavior of the NBS wave related to the scattering phase shift (the phase of the S-matrix) in QCD, but also the non-local but energy-independent potential (or the interaction kernel) can be extracted from the NBS wave function in the non-asymptotic region. By construction the potential correctly reproduces the scattering phase shift at all energies below the inelastic threshold. In practice the non-local potential is approximated by the velocity expansion in terms of local functions, so that physical observables such as the scattering phase shift are approximately calculated. It is also possible to check the accuracy of the approximation.
Using the local potential approximation at the leading order in the velocity expansion, nucleonnucleon, hyperon-nucleon and hyperon-hyperon interactions have already been investigated successfully. The BB potential in the flavor SU(3) limit shows the attraction in the flavor singlet channel, which is strong enough to form a bound state, H-dibaryon in this limit. The repulsive core is analyzed in terms of the operator product expansion and the renormalization group.
The potential method is extended to the case that the total energy is above the inelastic threshold. This method seems to work for the BB interaction with S = −2 and I = 0, where the coupled channel potentials among ΛΛ, NΞ and ΣΣ are obtained. The potential method is also applied to the three nucleon force, the meson-baryon interactions and two color QCD.
Of course, more careful studies of systematic errors such as the finite volume effect, dynamical quark effect, quark mass dependence and the lattice spacing effect are needed, before applying potentials obtained in lattice QCD in order to investigations of nuclei and hyper-nuclei.
Finally I stress here that the potential derived from the NBS wave function adds a new tool to investigate hadron interactions in lattice QCD, in addition to the standard finite volume method proposed by Lüscher [14]. In particular, the hadron scattering above the inelastic threshold can be treated in lattice QCD by the potential method. Other extensions of this method will also be looked for.