Spin-Orbit Force from Lattice QCD

We present a first attempt to determine nucleon-nucleon potentials in the parity-odd sector, which appear in 1P1, 3P0, 3P1, 3P2-3F2 channels, in Nf=2 lattice QCD simulations. These potentials are constructed from the Nambu-Bethe-Salpeter wave functions for J^P=0^-, 1^- and 2^-, which correspond to A1^-, T1^- and T2^- + E^- representation of the cubic group, respectively. We have found a large and attractive spin-orbit potential VLS(r) in the isospin-triplet channel, which is qualitatively consistent with the phenomenological determination from the experimental scattering phase shifts. The potentials obtained from lattice QCD are used to calculate the scattering phase shifts in 1P1, 3P0, 3P1 and 3P2-3F2 channels. The strong attractive spin-orbit force and a weak repulsive central force in spin-triplet P-wave channels lead to an attraction in the 3P2 channel, which is related to the P-wave neutron paring in neutron stars.


Introduction
A study of the nuclear force in QCD is a first step toward the understanding of hadronic properties beyond single hadrons. In addition, the nuclear force plays a key role to describe properties of atomic nuclei and neutron stars [1,2]. Recently, a new method to extract hadronic interactions from (lattice) QCD, where non-local potential can be defined through the Nambu-Bethe-Salpeter (NBS) wave function, has been proposed and successfully applied to the nuclear force [3,4,5]. The method has been being extended to various other systems such as hyperon-nucleon (YN), hyperon-hyperon (YY), meson-baryon, and the three-nucleons [6].
In the case of the nucleon-nucleon (NN) system, for example, non-local potentials are defined through the Schrödinger equation for the NBS wave function. Below the pion production threshold, the non-local NN potential can be expanded by the number of derivatives with respect to its nonlocality of the relative coordinate. The leading order (LO) terms are the spin-independent central potential V 0 (r), the spin-dependent central potential V σ (r) and the tensor potential V T (r), while the next-to-leading order (NLO) term is the spin-orbit potential V LS (r). Up to the NLO, there are altogether 8 independent local potentials, V I=0,1 0,σ,T,LS where I denotes the total isospin [4]. In the previous studies [3,4,5], T have been determined from the NBS wave functions in S and D waves (the parity-even sector) at various lattice parameters. For the complete determination, however, we must employ the NBS wave functions in P and F waves (the parity-odd sector) with non-zero relative momentum of the NN system, where relevant channels are 1 P 1 , 3 P 0 , 3 P 1 and 3 P 2 -3 F 2 . The corresponding potentials are given by . The main purpose of this paper is to extract these four potentials, in particular, the spin-orbit (LS) force, using the NBS wave functions in the parity-odd sector. The LS force has close relation to the spin-orbit splittings of the nuclear spectra and the nuclear magic numbers [7]. Large neutronneutron attraction due to the LS force in the 3 P 2 -3 F 2 channel leads to the P -wave super fluidity in the stellar environment such as the neutron star interiors [8,9,10]. It also affects the cooling properties of neutron stars [11].
The present paper is organized as follows. In Section 2, we give a brief review of our method to obtain NN potentials on the lattice. Section 3 is devoted to the discussion on the NBS wave functions with non-zero angular momenta. In Section 4, we present numerical results of the potentials. The scattering phases and mixing parameters in 1 P 1 , 3 P 0 , 3 P 1 , and 3 P 2 -3 F 2 channels calculated form these potentials are also presented.

NN potential in QCD
Below the inelastic threshold of the NN system (W < 2m N + m π ), we can define the non-local but energy-independent potential as [4] from the Nambu-Bethe-Salpeter (NBS) wave function in the center of mass (CM) frame defined by where H 0 ≡ −∇ 2 /m N with the nucleon mass m N , p α ( x), n β ( y) denote local composite nucleon operators with spinor indices α, β in Dirac representation restricted to α, β = 0, 1, |B = 2, W is the state with baryon-number B = 2 and vanishing total momentum with the total energy W . Due to the confinement of quarks and gluons, φ( r; W ) for large | r| reduces to the relative wave function of the non-interacting nucleons, so that W can be written as W = 2 k 2 + m 2 N with k ≡ | k| being the asymptotic relative momentum between the nucleons. An identity U ( r, r ) = V ( r, ∇)δ( r − r ) leads to the derivative expansion up to the NLO order: where S 12 ≡ 3( σ 1 · r)( σ 2 · r)/r 2 − σ 1 · σ 2 , S ≡ ( σ 1 + σ 2 )/2 and L ≡ i r × ∇ denote the tensor, the total spin and the (relative) orbital angular momentum operators, respectively, and I = 0, 1 is the total isospin of the two nucleons.
These NN potentials can be extracted by solving Eq. (4) with Eq. (6) for NBS wave functions projected to appropriate quantum numbers.

NBS wave functions with non-zero angular momenta
The NBS wave function for the ground state can be extracted from nucleon four-point functions at large t as where two-nucleon source J (J P ; t 0 ) located at the time-slice t = t 0 is used to control the quantum numbers such as J P , and W m denotes the energy of the intermediate state |m . The summation over y is performed to select states with the vanishing total spatial momentum. For p α (x) and n β (x), we employ, To construct sources coupled to the parity-odd states, we employ wall sources with non-zero momentum given by denotes a plane wave with the spatial momentum parallel or anti-parallel to a coordinate axis as An element g of cubic group O acts on these six functions as permutation, is the representation matrix of the cubic group whose explicit form can be generated by the basic matrices, with C 4y and C 4z being the rotations by +90 degree around y-axis and zaxis, respectively, and U (g) for other g ∈ O are obtained by multiplying these two matrices in suitable orders. For instance, the representation matrix for C 4x , a rotation by +90 degree around the x-axis, is obtained as . Analysis based on the group characters shows that U (g) is where S(g) denotes the standard rotation matrix acting on upper two components in the Dirac representation, we have for the spin triplet sector. 1 We introduce several projections to fix quantum numbers of the source operator J αβ (f i ). The projection for the total angular momentum J is given by where d (J) and χ (J) (g) denote the dimension and the character of the irreducible representation J of the cubic group, and the 24 × 24 matrix P (J) αβi;α β j is the projection matrix onto the irreducible representation J. Similar considerations gives the projection matrices for the parity 1 See Appendix A in Ref. [12] for a decomposition of a product of two irreducible representations of the cubic group and for the relation of irreducible representations between the cubic group and SO(3).
as well as the total spin S P (S=0) The projection matrices for J z are defined based on the C 4 subgroup of the cubic group which consists of multiples of C 4z as Since a product of these projection matrices 2 P (J,Jz,P,S) αβi;α β j ≡ P (J) P (Jz) P (P ) P (S) αβi;α β j .
has the property (P (J,Jz,P,S) ) 2 = P (J,Jz,P,S) , the eigenvalues of P (J,Jz,P,S) are either 0 or 1. We diagonalize P (J,Jz,P,S) to obtain its eigenvectors η where n is used to describe a possible degeneracy of a given set of quantum numbers J, J z , P, S. By construction, a state J (J P , J z , S)|0 has conserved quantum numbers J P , J z and S in the two-nucleon sector.

Lattice Setup
In this study, we employ N f = 2 full QCD configurations generated by CP-PACS Collaboration on a 16 3 × 32 lattice with the RG improved gauge action (Iwasaki action) at β = 1.95 and with the O(a)-improved Wilson quark (clover) action at κ = 0.1375 and C SW = 1.53, which gives the lattice spacing a = 0.1555(17) fm, the spatial extension L = 16a = 2.489(27) fm, the pion mass m π 1133 MeV and the nucleon mass m N 2158 MeV [13]. The Dirichlet boundary condition along the temporal direction is employed to generate quark propagators to avoid contamination from backward propagations of nucleons with negative parity.
With the projection defined in the previous section, we obtain the NBS wave functions for T − 1 in the spin singlet sector and for A − 2 ) in the spin triplet sector. In order to improve statistics, we perform the measurement on 32 source points by temporally shifting the location of the source, in addition to averages over the charge-conjugation/time-reversal transformations. We further reduce noises by the average over the cubic group, as will be discussed later.
To construct the potentials, we use the time-dependent method [5] with a slight modification to cope with the deviation of relativistic dispersion relation due to heavy quark mass, as explained in the following subsection. The nearest neighbor derivative is used to define the discretized laplacian, while the symmetric derivative is employed to define the operator L. To define S 12 and L on the periodic lattice, we take the origin of r to be the nearest periodic copy of the origin. On the spatial boundaries, i.e., x = ±L/2 or y = ±L/2 or z = ±L/2, however, these operators are still ill-defined. We therefore exclude data on these boundaries in our analysis. To extract potentials, we employ t − t 0 = 8, which is determined from t dependencies of potentials and phase shifts.

Modified time-dependent method
In Ref. [5], the time-dependent method has been proposed to extract the potential directly from the four-point functions. The method indeed gives more accurate and stable results, since it does not rely on the ground-state saturation at large t. We therefore employ this method also in this paper.
For a coarse lattice, however, the heavy quark may violate the relativistic dispersion relation of a single nucleon as E 2 m 2 N + α k 2 with α = 1 [14]. In such a case, the formula in Ref. [5] receives a slight modification as In our simulation, the dispersion relation for the nucleon can be fitted well with α = 0.88 (1)  of higher order contributions in k 2 for k 2 ≤ 1.25[GeV 2 ] (ka ≤ √ 5 × 2π/L) within statistical errors.

Extractions of potentials
The potential for the spin-singlet sector at NLO can be easily extracted from the equation that ∂t 2 − ∂ ∂t , and we define an inner product with an average over the cubic group as F ( r), H( r) ≡ g∈O F * βα (g r)H αβ (g r), which reduces statistical noises of potentials. Note that here and in the following we use the fact that local potentials, V I C,S , V I T and V I LS , are invariant under the rotation g in the cubic group. The result for V I=0 C,S=0 (r) is plotted in Fig.1 by green circles, which shows a strong repulsion at short distances.
For the spin-triplet sector, three unknown functions up to NLO can be determined from the equation that In Fig. 1, we also plot V I=1 C;S=1 (r)(red), V I=1 T (r)(black) and V I=1 LS (r)(blue), sources instead does not show a significant difference.) We observe that (i) the central potential V I=1 C;S=1 (r) is repulsive, (ii) the tensor potential V I=1 T (r) is positive and weak compared to V I=1 C;S=1 (r) and V I=1 LS (r), and (iii) the spin-orbit potential V I=1 LS (r) is negative and strong. These features agree qualitatively well with those of the phenomenological potential in Ref. [15].
For both spin-singlet and spin-triplet central potentials, there may be a very weak attractive pocket of less than a few MeV at medium distance (r 1 fm). However, considering the statistical and systematic errors, its existence should be carefully examined in future studies.
We make a technical comment. We sometimes observe large condition numbers for eq. (19) (with three sources) near the spatial boundaries, which gives rise to points with large statistical errors at r 1 − 1.5 fm in Fig. 1.

Scattering phase shifts and effective potentials
We calculate the scattering phase shifts by solving the Schrödinger equation with the above potentials, parametrized with multi-gaussian forms, v(r) ≡ Ngauss i=1 a i exp(−ν i r 2 ) with N gauss = 3 for the central and spin-orbit potentials, whereas v(r) ≡ a 1 r exp(−ν 1 r 2 ) + a 2 r 3 exp(−ν 2 r 2 ) for the tensor potential to mimic the short distance behavior, as shown in Fig 1. The scattering observables are obtained from the long distance behaviors of linearly independent regular solutions, and are shown in Fig.2. Although the magnitude of the phase shifts obtained from our potentials are smaller than the experimental ones, general trends are well reproduced except for 3 P 0 case at low energies. The missing attraction in 3 P 0 channel is likely due to the weak tensor force V T caused by the large pion mass. Among others, the most interesting feature in Fig.2 is the attraction in the 3 P 2 channel, which is directly related to the paring correlation of the neutrons inside the neutron stars.
To obtain an intuitive understanding of the behavior of these phase shifts, we plot the potentials of the 1 P 1 , 3 P 0 , 3 P 1 and 3 P 2 channels, as defined in Eqs. (1)-(3) and below, in Fig.3. Indeed, one can see that V (r; 3 P 2 ) has a weak repulsive core surrounded by an attractive well; the attraction is driven by the strongly attractive LS force, V I=1 LS (r) in Fig.1. Note that we observe small violations of rotational invariance at short distance.  V(r; 1 P 1 ) V(r; 3 P 0 ) V(r, 3 P 1 ) V(r; 3 P 2 ) Figure 3: The potentials for 1 P 1 , 3 P 0 , 3 P 1 and 3 P 2 channels given in Eqs. (1)-(3) and below.

Conclusion
We have made a first attempt to determine NN potentials up to NLO in the parity-odd sector, which appears in 1 P 1 , 3 P 0 , 3 P 1 , 3 P 2 -3 F 2 channels. Using N f = 2 CP-PACS gauge configurations on a 16 3 × 32 lattice (a 0.16 fm and m π 1100 MeV), not only the central and the tensor potentials but also the spin-orbit potential have been derived for the first time. These potentials are constructed from NBS wave functions for J P = 0 − , 1 − , 2 − , which are generated by using the momentum wall sources with projections based on the representation theory of the cubic group.
We have observed that the qualitative behavior of the resultant potentials agree with those of phenomenological potentials: For the spin-singlet sector, the central potential V I=0 C;S=0 (r) is repulsive with a strong repulsive core at short distance. For the spin-triplet sector, (i) the central potential V I=1 C;S=1 (r) is also repulsive with a repulsive core at short distance, (ii) the tensor potential V I=1 T (r) is positive and quite weak, and (iii) the spin-orbit potential V I=1 LS (r) is negative and strong at short distance. We have then calculated scattering observables in 1 P 1 , 3 P 0 , 3 P 1 and 3 P 2 -3 F 2 channels, by solving Schrödinger equations with these potentials. It is interesting enough that we obtain, from the first principle lattice QCD approach, attractive phase shift driven by the strongly attractive LS force in the 3 P 2 channel, which has been known experimentally and has various implications in atomic nuclei and dense matter, though the magnitude of the phase shift is still small due to the large quark mass in our calculation.