Hadron Interactions from lattice QCD

We review our strategy to study hadron interactions from lattice QCD using newly proposed potential method. We first explain our strategy in the case of nuclear potentials and its application to nuclear physics. We then discuss the origin of the repulsive core, by adding strange quarks to the system. We also explore a possibility for H-dibaryon to exist in flavor SU(3) limit of lattice QCD. We conclude the paper with an application of our strategy to investigate the maximum mass of neutron stars.


Introduction
As numerical simulations in lattice QCD mature, we can attack more difficult or complicated problems in strong interactions.Calculations of hadronic matrix elements of electroweak operators are predictions from lattice QCD, which can be used to determine parameters of the standard model.Another interesting application is to investigate interactions among hadrons using lattice QCD.Since hadrons are bound states of quarks and gluons, their interactions are residual effects of strong interactions among them inside hadrons, and therefore they are more complicated quantities in lattice QCD than properties of an isolated hadron such as masses or decay constants.In addition, the euclidean nature of lattice QCD makes problems of hadron interactions more complicated and difficult.Thanks to the finite volume technique [1], however, scattering lengths and phase shifts of hadrons can in principle be extracted from lattice QCD.Some recent progress on these quantities is reviewed in this paper.
I also consider a more complicated but important quantity of hadron interactions, the force between nucleons (the nuclear force).In 1935 Yukawa introduced a virtual particle, the pion, to account for the nuclear force [2], by which protons and neutrons are confined in nuclei.Since then enormous efforts have been made to understand the nucleon-nucleon (NN) interaction at low energies both theoretically and experimentally.In Fig. 1, I present modern NN potentials, which are characterized by the following features [3,4].At long distances (r ≥ 2 fm ) there exists weak attraction, which is well understood and is dominated by the one pion exchange (OPE), as first pointed out by Yukawa.At medium distances (1 fm ≤ r ≤ 2 fm), contributions from the exchange of multi-pions and heavy mesons such as ρ, ω and σ lead to slightly stronger attraction.At short distances (r ≤ 1 fm), attraction turns into repulsion, and it becomes stronger and stronger as r gets smaller, forming the strong repulsive core [5], which is essential not only for describing the Figure 1: Three examples of the modern NN potential in 1 S 0 (spin singlet and S-wave) channel: Bonn [8], Reid93 [9] and AV18 [10].NN scattering data, but also for the stability and saturation of atomic nuclei, for determining the maximum mass of neutron stars, and for igniting Type II supernova explosions [6].Although the origin of the repulsive core must be related to the quark-gluon structure of nucleons, it remains one of the most fundamental problems in nuclear physics for a long time [7].It is a great challenge for us to derive the nuclear potential including the repulsive core from lattice QCD.Recent progress on this issue is explained in this paper.

Conventional method in lattice QCD
In this section, some recent results for the scattering length and phase shift for pions are discussed.

Scattering length of pions
In Fig. 2, recent theoretical estimates of the S-wave ππ scattering lengths, a 0 0 and a 2 0 in units of 1/M π , are summarized [11], where L and I of the scattering length a I L represent the total angular momentum and the total isospin, respectively.The band bounded by two black solid lines represents the theoretically allowed region (universal band) [11].The three black dots are predictions from 2-flavor chiral perturbation theory (ChPT) at tree, 1-loop and 2-loop levels from left to right.Other results are obtained by using the formulae [12] M π (2a where M 4 π α 3,4 are corrections to ChPT from dispersion integrals.They are approximated at the physical point as M 4  π α 3 = 0.135 + 0.77(a 0 0 − 0.220) − 1.50(a 2 0 + 0.0444), ( M 4 π α 4 = 0.061 + 0.48(a 0 0 − 0.220) − 0.26(a 2 0 + 0.0444).
(2.4) l3,4 are the low energy constants in N f = 2 ChPT and defined by Here m is the quark mass and F is the pion decay constant in the chiral limit.Once l3 and l4 are known, the above expressions give a 0 0 and a 2 0 .The red ellipse in the figure corresponds to l3 = 2.9(2.4) from the mass spectrum of the pseudoscalar octet and l4 = 4.4(2) from the scalar form factor of the pion.The narrow azure strip indicates the region allowed if l3 is treated as free parameter while l4 is fixed to the above value.There exist several lattice determinations of l3 and l4 .The MILC collaboration [13], using 2+1 flavors of staggered dynamical quarks, obtain the low energy constants L 4,5,6,8 of N f = 3 ChPT, which are translated to l3 = 0.8(2.3)and l4 = 4.0(6) by standard one loop formulae.The results for a 0 0 and a 2 0 are given by a green ellipse in the figure, which is consistent with the previous estimate (the red ellipse).The result of N f = 2 dynamical Wilson quarks, l3 = 3.0(5)(1) [14], with free l4 leads to the narrow purple strip in the figure, which intersects the red ellipse.
These lattice results for a 0 0 and a 2 0 are indirectly obtained by using l3 and l4 , which have been determined from the quark mass dependences of M π and F π .There exist several direct estimates on a 2 0 using the finite volume method [1].The horizontal orange band in Fig. 2 indicates the direct result for a 2 0 by the NPLQCD collaboration [16] with domain-wall valence quarks on a 2+1 staggered sea.Recently they have reported a more accurate result, M π a 2 0 = −0.0443(4)[17], which is consistent with the red ellipse.So far ChPT predictions and all indirect /direct lattice results agree with each other.However a direct evaluation of a 0 0 is missing.Thus it is a great challenge for lattice QCD to extract a 0 0 directly by the finite volume method.

The ρ meson decay width
A less difficult than a 0 0 but still challenging problem is the determination of the P-wave scattering phase shift for the I = 1 two-pion system, from which the ρ meson decay width is extracted.Applying the finite volume method in the laboratory system to this problem, the CP-PACS collaboration has calculated the decay width [18] in lattice QCD with a renormalization group improved gauge action and N f = 2 dynamical clover quarks on a 12 3 × 24 lattice at m π /m ρ ≃ 0.41 and the lattice spacing 1/a ≃ 0.92 GeV.In order to realize the kinematics such that the energy of the two pions is close to m ρ , one pion has a non-zero momentum p = (2π/L)e 3 and the other is at rest, while the ρ meson has the same non-zero momentum.Energies of these states for non-interacting hadrons are W 0 1 = m 2 π + p 2 + m π for the two pions and W 0 2 = m 2 ρ + p 2 for the ρ.At m π /m ρ ≃ 0.41 on π + p 2 ≃ 1.47m ρ .The hadron interaction shifts the energies from W 0 n to W n , which are related to the two-pion scattering phase shift δ in the infinite volume through the Rummukainen-Gottlieb formula [19].
To extract two energy levels W n close to each other, a 2 × 2 matrix of the time correlation function has been constructed, where ρ 3 (t) is an interpolating operator for the neutral ρ meson with the momentum p and the polarization vector parallel to p, and (ππ)(t) is an interpolating operator for two pions, (ππ for large t.To keep symmetry between source and sink in G(t) U( 1) noises are introduced.The total number of quark propagators per configuration becomes 520, including 10 noises times 2 source points [18].
From the energy levels W n extracted, the invariant mass and the momentum k are given by π , assuming the continuum dispersion relation, while they become cosh( In both cases it is observed that δ > 0 (attractive) at √ s < m ρ = 0.86(1) while δ < 0 (repulsive) at √ s > m ρ .This property confirms the existence of a resonance at a mass around m ρ [20].
In Fig. 3, sin 2 δ , which is proportional to the scattering cross section of the two-pion system, is plotted against the invariant mass √ s.The two lines in the figure represent fits of the data by tan δ = g ρππ 6π where M R is the resonance mass, plotted also in the upper panel of Fig. 3, and g ρππ is the effective coupling of ρππ.If g ρππ were obtained at several quark masses, it could be extrapolated to the physical quark mass.Since this calculation is performed only at one quark mass, however, it is assumed that g ρππ is independent of the quark mass, so that the ρ meson decay width at the physical quark mass is estimated by where The fit result of g ρππ gives Γ ph = 162(35) MeV (continuum) and Γ ph = 140( 27)MeV (lattice), which are consistent with the experiment value, Γ = 150 MeV.
Since the first attempt is successful, the next step is the chiral extrapolation.The ρ meson decay width is a good bench mark quantity which shows the dynamical nature of QCD vividly.Therefore collaboration groups which have full QCD configurations at small quark masses are encouraged to calculate this quantity.

Potentials for heavy hadrons
In the following two sections, extractions of the potential between hadrons from lattice QCD are discussed.A straight-forward way to define a potential is to calculate an energy of a twohadron system as a function of the distance between the two hadrons.A difficulty exists in the definition of the distance between two hadrons, since the two hadrons are moving around changing their relative distance.To overcome this difficulty, one (infinitely) heavy quark may be introduced in each hadron, so that the distance between two hadrons is defined by the distance between the two heavy static quarks, which do not move in space for all the time.This definition, similar to the static quark potential from a Wilson loop or Polyakov lines, is indeed employed to calculate potentials between two heavy hadrons in lattice QCD.

Baryon-Baryon
In Ref. [21], the energy of two heavy baryons has been measured as a function of the distance r between two heavy quarks Q, using quenched QCD with the plaquette gauge action and the Wilson quark action at a ≃ 0.19 fm on a 20 3 × 24 lattice.The binding energy of two baryons, V BB − 2V B , is plotted as a function of r in Fig. 4 at m π ≃ 500 MeV for the light quark, where V BB is an energy for two baryons while V B is an energy for one baryon.The differences of the three figures are explained in the caption.Surprisingly, the binding energy is very small and shows almost no dependence on r for all three cases.Thus, a repulsive core is not seen for the heavy baryon potential in this calculation.

Meson-Meson and others
In Ref. [22], the binding energy of two heavy-light mesons has been computed, each of which is made of one static quark and one light quark, using quenched QCD with the DBW2 gauge action at a ≃ 0.1 fm and the Wilson quark action at m π ≃ 400 MeV.The central potential for total spin S and isospin I of two light quarks in a meson-meson system is decomposed as where σ i (τ i ) acts on the spin (isospin) of the light quark q i , and r is the distance between the two static quarks.Results for V X (r) (X = 1, σ , τ, σ τ) are plotted in Fig. 5.All potentials show attraction at short distance (r ≤ 0.2 fm ).In Ref. [23], a meson-meson potential and a meson-baryon potential has been calculated with the same setup of the baryon-baryon case in the previous subsection.As shown in Fig. 6, however, the dependence on r is very small in both cases.In particular, no short-distance attraction is observed for the two heavy-light mesons, contrary to the results in Fig. 5.There are several candidates which might cause the difference: the light quark mass, the lattice spacing, the lattice volume, the statistics or the way to extract the binding energy.Further investigations are necessary for definite conclusions about the static quark approach to potentials between hadrons.

Potentials from wave functions
In Ref. [24], the nucleon-nucleon (NN) potential has been extracted in lattice QCD by a totally new and different method, which is explained in this section.

Wave functions and potentials
The starting point of the new method for the NN potential is a wave function defined by where N(x,t) is an interpolating field of the nucleon, |2N; E is a 2N state with energy E < E inelastic with the inelastic threshold energy E inelastic of the 2N system.This wave function was first analyzed for a spin model in Ref. [25], and has been employed in lattice QCD to calculate the ππ scattering length [26].This wave function has the following properties.For large enough r = |r|, is a free non-relativistic Hamiltonian with the reduced mass µ = m N /2.This means that the interaction between two nucleons vanishes for a large separation.Moreover, for example in the case of the 1 S 0 channel, it can be shown that ϕ S (r, E) ∼ e iδ 0 (k) sin(kr where δ 0 (k) is the S wave scattering phase shift and k is determined by E = k 2 /(2µ).In Ref. [26,27], it is shown for a two-pion system that ϕ S (r, E) = e ikr + ik 4πr H(k; k) j 0 (kr) + P where π + p 2 and M(p, k) is the (off-shell) scattering amplitude.A similar but a little more complicated expression can be shown also for the NN system [28].For large r, if E < E inelastic , the contribution from I others vanishes exponentially and eq.( 4.3) approaches eq.( 4.2).
The potential is defined from the wave function through the Schrödinger equation which symbolically gives

Lattice calculation
A corresponding lattice wave function in a finite box at E ≃ 0 is given explicitly by where N i α = ε abc ( t q a Cγ 5 τ 2 q b )q i,c α , P 1 = 1 and P 0 = τ 2 are isospin projections, P 1 = 1 and P 0 = σ 2 are spin projections, and the summation over x gives zero total momentum.The summation over the discrete rotations R of the cubic group O implies that the state belongs to the A + 1 representation of the cubic group, which is expected to couple to an L = 0 ground state as well as L ≥ 4 excited states of the rotation group in the continuum theory.
The wave function without projections has been extracted from the 4-point nucleon correlator as where A n = 2N; E n |J NN (t 0 )|0 and J NN (t 0 ) = P I i j P S αβ N i α N j β .In N , a wall source is employed by replacing q(x,t 0 ) with Q(t 0 ) = ∑ x q(x,t 0 ) after Coulomb gauge fixing.For large enough t the wave function for E = E 0 is obtained.
In a finite volume, energy levels are shifted from those in the infinite volume as ) due to interactions.From these shifts, the scattering phase shift can be extracted [1].

Numerical simulations
In the actual calculation, the NN scattering for L = 0 is considered.There are two channels, the spin singlet (S = 0) channel 1 S 0 and the spint triplet (S = 1) channel 3 S 1 , where the standard notation 2S+1 L J is used.Since only a central potential appears for 1 S 0 , the definition (4.5) directly gives V NN (r) = V C (r).On the other hand, in the case of 3 S 1 , the potential becomes V NN (r) = V C (r) + V T (r)S 12 (r), where S 12 (r) = 3(σ 1 • r)(σ 2 • r) − (σ 1 • σ 2 ) is the tensor operator.The tensor potential induces mixing between 3 S 1 and 3 D 1 states, so that the definition (4.5) gives the so-called effective central potential, which includes the mixing effect from the 3 D 1 state as a second order perturbation.

Results
In the left panel of Fig. 7, the NN wave function normalized to 1 at r ≃ 2.2 fm is plotted as a function of r for 1 S 0 and 3 S 1 at m π ≃ 527 MeV.In both cases the wave function shows an increase at 0.5 fm < r < 1.5 fm, suggesting attraction, while it decreases at r < 0.5 fm, indicating the existence of the repulsive core.In the right panel of Fig. 7 the central (effective central) NN potential extracted from the wave function by eq.( 4.5) at m π ≃ 527 MeV is plotted as a function of r for 1 S 0 ( 3 S 1 ).Interestingly the potential obtained qualitatively agrees with the NN potential determined from scattering experiments: weak attraction at long distance, a little stronger attraction at intermediate distance and strong repulsion at short distance (the repulsive core).The solid line is the Yukawa potential (One Pion Exchange Potential) given by which agrees well with the data at long distance, with g 2 πN /(4π) ≃ 14.0 from experiments, m π ≃ 0.527 GeV and m N ≃ 1.34 GeV from lattice data.Thus, in some sense, the Yukawa theory for the nuclear force at r > 1 fm is confirmed by lattice QCD.
The quark mass dependence of the NN potential for 1 S 0 is shown in Fig. 8.As the quark mass decreases, the repulsive core at short distance gets stronger, and at the same time, attraction at intermediate distances also becomes a little stronger [29].In the future it will be interesting to see if this quark mass dependence remains in full QCD.

Energy dependence of the potential
Since the wave function depends on the energy E = p 2 /(2µ), the potential defined by eq.(4.5) is also energy dependent in general: where J, the total angular momentum, is fixed, and r = |r|.To make our argument simpler, the spin degrees of freedom are not considered here.Since V J (r, E) carries more information than the scattering phase shift δ (E) does, V J (r, E) is redundant and therefore not physical.For the potential defined from the wave function to be physically meaningful, its energy dependence must be weak for some range of small E. In Fig. 9, the ππ wave function and the corresponding potential in the S channel are shown at r = (x, y, 0) in quenched QCD at m π /m ρ ≃ 0.51 for p = 0 (left) and p = (2π/L, 0, 0) (right) on an L 3 box [30] .Although the wave functions at the two different energies are very different, the potentials look similar.For the ππ system, the energy dependence of the potential is not so strong at low energy.Note, however, that the wave function (and therefore the potential) at p = (2π/L, 0, 0) is calculated at equal time in the laboratory system and is transformed back to the center of mass system, so that the relative time between the two pions is x dependent: t = vγx.

A unique local potential
In this subsection a method to extract a unique potential from wave functions is proposed [28].As mentioned in the previous subsection, the potential seems E dependent in general.In addition, it may also depend on the choice of the interpolating operator N(x,t) in (4.1): one can use a different Ñ(x,t) unless it changes the asymptotic behavior (the phase shift).To begin with, it is argued that, by introducing a non-local potential, the energy dependence of the potential can be removed.The non-local potential is defined by the following equation, where all quantities are taken to be real.For notational simplicity, ∑ is used even for the case that a variable is continuous and the superscript J for the angular momentum is omitted here.Since {ϕ(r, E)} E form a complete set, there exists an inverse such that Using this we obtain (5.4) In the actual simulations, it is impossible to obtain ϕ(r, E) for all E. If ϕ(r, E)'s are obtained for (5.5) The coefficient U k can be obtained as where φ −1 j,k is the inverse of φ j,k = d k φ (r, E j ).Note that d becomes a difference on a lattice.It is now argued that the non-local potential U obtained above can be transformed to a local E independent potential V .The Schrödinger equations for them become where H 0 + U is not hermitian unless U (r, r ′ ) = U (r ′ , r), while H 0 + V is hermitian.Note here that the eigenvalues E in both equations must be equal to ensure that both wave functions give the same scattering phase shift (and the bound state spectrum if any).Existence and uniqueness of V are suggested by the inverse scattering theory, which tells us that a local potential for a fixed J is uniquely reconstructed from the scattering phase shift for all E and informations of possible bound states.Since the eigenfunctions |ψ, E of the hermitian operator can be constructed as The compatibility of non-local and local potentials leads to which, together with the completeness of |ϕ, E , implies (5.12) If the number of degrees of freedom for E and r is denoted as #E = #r = N on a finite lattice, eq.(5.12) gives N × N constraints, while the unkown Λ and V have N × N and N components, respectively.Therefore the determination of both Λ and V from eq.(5.12) seems an ill-posed problem.However, N(= #E) components of Λ can be taken freely, since the eigenvalue equation (5.8) does not depend on C E , the norm of the eigenfunction, for all E's.Using this freedom, eq.(5.12) is enough to determine Λ and V .Eq.(5.12) is the master equation, which determines the unique local energy independent potential V , once a non-local potential U is obtained from a wave function defined through a particular choice of operators.This conceptually solves the uniqueness (or E dependence) problem of "the potential from a wave function".The results in the previous section correspond to the 0th order solution of this equation: V = U 0 and Λ = 1.

Future perspectives
To understand interaction properties of hyperons, baryons which include at least one strange quark, is an important subjects in the nuclear physics.Hyperon-nucleon (Y N) and hyperon-hyperon (YY ) interactions are relevant to structures of the neutron-star core and the existence/absence of Hdibaryon states.Properties of hypernuclei, nuclei which contain hyperons, will be also studied, as a project at J-PARC (Japan Proton Accelerator Research Complex).However, Y N and YY interactions are poorly known both theoretically and experimentally so far.Lattice QCD calculations of the NN potential in the previous section can be extended to Y N and YY potentials.
In Fig. 10, the wave function and the corresponding potential for pΞ 0 are plotted as a function of r in quenched QCD [31].The lattice parameters are the same as in the case of the NN potential in the previous section.The light quark mass corresponds to m π ≃ 370 MeV while the strange quark mass is tuned to reproduce the physical K meson mass, m K ≃ 550 MeV.Qualitative features of the Y N potential are similar to those of the NN potential.Weak attraction appears at long and intermediate distances while a strong repulsive core shows up at short distance.However the spin dependence of the potential is stronger than in the NN case.In particular, the repulsive core in the 1 S 0 channel is much stronger than that in the 3 S 1 channel.Although these results are still preliminary, they are interesting and encouraging as a first step.
Lattice QCD calculations of potentials between hadrons from wave functions have just begun and the first result of the NN potential surprisingly reproduces all the known features of the NN potential such as weak attraction at long distance, a little stronger attraction at intermediate distance and a strong repulsive core at short distance.For a quantitative comparison between lattice results and experimental ones, however, chiral and continuum extrapolations are necessary to remove systematic errors.The inclusion of dynamical quark effects in the NN potential will be the most exciting improvement in future calculations.Now a door is open for us to vast fields in nuclear physics with lattice QCD.

Figure 3 :
Figure 3: Scattering phase shift sin 2 δ (lower panel), positions of m ρ and resonance mass M R (upper panel) in lattice units."Cont" refers to results obtained with the continuum dispersion relation while "Lat" to those obtained with the lattice dispersion relation.a 12 3 × 24 lattice, the invariant mass of two free pions takes the value √ s ≃ 0.97 × m ρ , which is much closer to m ρ than that in the center of mass system, E = 2 m 2π + p 2 ≃ 1.47m ρ .The hadron interaction shifts the energies from W 0 n to W n , which are related to the two-pion scattering phase shift δ in the infinite volume through the Rummukainen-Gottlieb formula[19].To extract two energy levels W n close to each other, a 2 × 2 matrix of the time correlation function using a lattice dispersion relation.Then W 1 leads to tan δ = 0.0791(9) at √ s = 0.776(8) (continuum dispersion) or tan δ = 0.074(7) at √ s = 0.800(8) (lattice dispersion), and W 2 gives tan δ = −0.38(11)at √ s = 0.95(3) (continuum) or tan δ = −0.52(14)at √ s = 0.98(3) (lattice).

Figure 4 :
Figure 4: The binding energy, V BB − 2V B (GeV), as a function of r (fm).Left: All the light-quark flavors are different.Middle: One pair of quark flavors is identical.Right: Two pairs are identical.

Figure 5 :
Figure 5: The central meson-meson potential (MeV) as a function of r in lattice units.

Figure 6 :
Figure 6: The meson-meson potential (left) and the meson-baryon potential as a function of r at m π ≃ 700 MeV.Triangles (circles) denote results without (with) light quark exchange diagrams.

Figure 7 :
Figure 7: Figures from [24].Left: The normalized NN wave function for the singlet(circles) and the triplet(triangles) as a function of r (fm) at m π ≃ 527 MeV.Right: The central (effective central) NN potential (MeV) as a function of r (fm) for the singlet (triplet) at m π ≃ 527 MeV.Numerical simulations have been made in quenched QCD on a 32 4 lattice with the plaquette gauge action and the Wilson quark action, at a ≃ 0.137 fm from the ρ meson mass.The number of configurations is 2000 for m π ≃ 370 MeV and 527 MeV, and 1000 for m π ≃ 732 MeV.A Dirichlet boundary condition (DBC) in time and a periodic B.C. in space are employed, and the wall source with Coulomb gauge fixing is placed at t 0 = 5, to avoid an influence of the DBC.Calculations have been performed on a Blue Gene/L at KEK, which has 57.3 TFlops peak performance.It took about 4000 hours of 512 Nodes (a half-rack, 2.87TFlops peak) with 34-48% sustained speed to complete them.