SpinOrbit Force from Lattice QCD
Abstract
We present a first attempt to determine nucleonnucleon potentials in the parityodd sector, which appear in the , , , – channels, in lattice QCD simulations. These potentials are constructed from the NambuBetheSalpeter wave functions for and , which correspond to the and representation of the cubic group, respectively. We have found a large and attractive spinorbit potential in the isospintriplet 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 the , , and – channels. The strong attractive spinorbit force and a weak repulsive central force in spintriplet wave channels lead to an attraction in the channel, which is related to the wave neutron paring in neutron stars.
keywords:
Lattice QCD, nuclear force, spinorbit potential, scattering phase shift1 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 in describing properties of atomic nuclei and neutron stars (1); (2). In lattice QCD, the standard method to study hadronic interactions is the finite volume method (3) by calculating the scattering phase shift(4); (5); (6); (7); (8); (9). Recently, a new method to extract hadronic interactions from (lattice) QCD has been proposed, where nonlocal potential can be defined through the NambuBetheSalpeter (NBS) wave function. The method has been successfully applied to the nuclear forces (10); (11); (12). It has been extended to various other systems such as hyperonnucleon (YN), hyperonhyperon (YY), mesonbaryon, and the threenucleons(13). In Ref. (14), an explicit comparison between the finite volume method and the potential method is made in the case of scattering phase shifts, where a good agreement is obtained.
In the case of the nucleonnucleon (NN) system, for example, nonlocal potentials are defined through the Schrödinger equation for the NBS wave function. Below the pion production threshold, the nonlocal 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 spinindependent central potential , the spindependent central potential and the tensor potential , while the nexttoleading order (NLO) term is the spinorbit potential . Up to the NLO, there are altogether 8 independent local potentials, where denotes the total isospin(11). In the previous studies(10); (11); (12), , and have been determined from the NBS wave functions in and waves (the parityeven sector) at various lattice parameters. For the complete determination, however, we must employ the NBS wave functions in and waves (the parityodd sector) with nonzero relative momentum of the NN system, where relevant channels are , , and –. The corresponding potentials are given by
(1)  
(2)  
(3)  
with , and . The LS force has close relation to the spinorbit splittings of the nuclear spectra and the nuclear magic numbers (15). Large neutronneutron attraction due to the LS force in the – channel leads to the wave superfluidity in the stellar environment such as the neutron star interiors (16); (17); (18). It also affects the cooling properties of neutron stars (19).
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 nonzero angular momenta. In Section 4, we present numerical results of the potentials. The scattering phases and mixing parameters in the , , , and – channels calculated form these potentials are also presented.
2 NN potential in QCD
Below the inelastic threshold of the NN system (), we can define the nonlocal but energyindependent potential as(11)
(4) 
from the NambuBetheSalpeter (NBS) wave function in the center of mass (CM) frame defined by
(5) 
where with the nucleon mass , , denote local composite nucleon operators with spinor indices , in Dirac representation restricted to , is the state with baryonnumber and vanishing total momentum with the total energy . Due to the confinement of quarks and gluons, for large reduces to the relative wave function of the noninteracting nucleons, so that can be written as with being the asymptotic relative momentum between the nucleons. The identity leads to the derivative expansion up to the NLO order:
(6) 
where , and denote the tensor, the total spin and the (relative) orbital angular momentum operators, respectively, and 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.
The NBS wave function in eq. (5) and thus the nonlocal potential defined in eq. (4) depend on the choice of nucleon operators and . This is not surprising since the potential is not a physical observable and depends on how it is defined. On the other hand, it has been shown that the NBS wave function carries information of the scattering phase shift in its asymptotic behavior at large (20); (21); (22); (11). This property holds for an arbitrary choice of operators to define the NBS wave function, as long as basic properties in quantum field theories such as locality and unitarity are satisfied. Therefore, by construction, nonlocal potentials in different definitions give same and correct phase shifts through the Schrödinger equation in eq. (4). Despite that the above points have been already stated clearly and explicitly in our previous papers (see, e.g. (11)), similar remarks appear repeatedly in later literature (see e.g. (24)).
In our study, we take local interpolating operators for the nucleon as , and , to define the NBS wave function, where denotes the color indices. In our actual calculation, the derivative expansion of the nonlocal potential in eq. (6) is truncated at fixed order (the NLO in the present paper); then the resultant phase shifts are valid only in a certain energy interval which depends on the order of the truncation and the choice of the interpolating operator. The former dependence has been investigated for the NN case(23) and the case(14), while the latter dependence is still left for future studies.
3 NBS wave functions with nonzero angular momenta
The NBS wave function for the ground state can be extracted from nucleon fourpoint functions at large as
(7)  
where twonucleon source located at the timeslice is used to control the quantum numbers such as , and denotes the energy of the intermediate state . The summation over is performed to select states with the vanishing total spatial momentum.
To construct sources coupled to the parityodd states, we employ wall sources with nonzero momentum given by
(8) 
with , and , where denotes a plane wave with the spatial momentum parallel or antiparallel to a coordinate axis as
An element of cubic group acts on these six functions as permutation, where is the representation matrix of the cubic group whose explicit form can be generated by the basic matrices,
(9) 
with and being the rotations by +90 degrees around yaxis and zaxis, respectively, and for other are obtained by multiplying these two matrices in suitable orders. For instance, the representation matrix for , a rotation by +90 degrees around the xaxis, is obtained as . The spatial reflection corresponds to , , and , which leads to representation matrix . Analysis based on the group characters shows that is reduced to the direct sum of irreducible representations .
Together with the transformation property of the quark field where denotes the standard rotation matrix acting on upper two components in the Dirac representation, we have
(10)  
Therefore our momentum wall source
covers () for the spin singlet
sector and () for the spin triplet sector.
We introduce several projections to fix quantum numbers of the source operator . The projection for the total angular momentum is given by
(11)  
where and denote the dimension and the character of the irreducible representation of the cubic group, and the matrix is the projection matrix onto the irreducible representation . Similar considerations give the projection matrices for the parity
(12) 
as well as the total spin
(13)  
The projection matrices for are defined based on the subgroup of the cubic group which consists of multiples of as
(14) 
Since a product of these projection matrices
(15) 
has the property , the eigenvalues of are either 0 or 1. We diagonalize to obtain its eigenvectors with eigenvalue 1 as which is used to perform the projection of our twonucleon source as
(16) 
where is used to describe a possible degeneracy of a given set of quantum numbers . By construction, a state has conserved quantum numbers and in the twonucleon sector.
4 Numerical Results
4.1 Lattice Setup
In this study, we employ full QCD configurations generated by the CPPACS Collaboration on a lattice with the RG improved gauge action (Iwasaki action) at and with the improved Wilson quark (clover) action at and , which gives the lattice spacing fm, the spatial extension fm, the pion mass MeV and the nucleon mass MeV(25). The Dirichlet boundary condition along the temporal direction is employed to generate quark propagators to avoid contamination from backward propagation of nucleons with negative parity.
With the projection defined in the previous section, we obtain the NBS wave functions for in the spin singlet sector and for in the spin triplet sector. In order to improve statistics, we perform the measurement on source points by temporally shifting the location of the source, in addition to averages over the chargeconjugation/timereversal transformations. We further reduce noises by the average over the cubic group, as will be discussed later.
To construct the potentials, we use the timedependent method (12) 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 . To define and on the periodic lattice, we take the origin of to be the nearest periodic copy of the origin. On the spatial boundaries, i.e., or or , however, these operators are still illdefined. We therefore exclude data on these boundaries in our analysis. To extract potentials, we employ , which is determined from dependencies of potentials and phase shifts.
4.2 Modified timedependent method
In Ref. (12), the timedependent method has been proposed to extract the potential directly from the fourpoint functions. The method indeed gives more accurate and stable results, since it does not rely on the groundstate saturation at large . 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 with (26). In such a case, the formula in Ref. (12) receives a slight modification as
(17) 
where .
In our simulation, the dispersion relation for the nucleon can be fitted well with ( ) at MeV, showing no sign of higher order contributions in for () within statistical errors.
4.3 Extractions of potentials
The potential for the spinsinglet sector at NLO can be easily extracted from the equation
(18) 
for , dominated by , where , and we define an inner product with an average over the cubic group as , which reduces statistical noises of potentials. Note that here and in the following we use the fact that local potentials, , and , are invariant under the rotation in the cubic group. The result for is plotted in Fig.1 by green circles, which shows a strong repulsion at short distances.
For the spintriplet sector, three unknown functions up to NLO can be determined from the equation
(19) 
for three different sources, , , ( or ), dominated by , and –, respectively, where
In Fig. 1, we also plot (red), (black) and (blue), obtained from sources. (The result obtained form sources instead does not show a significant difference.) We observe that (i) the central potential is repulsive, (ii) the tensor potential is positive and weak compared to and , and (iii) the spinorbit potential is negative and strong. These features agree qualitatively well with those of the phenomenological potential in Ref. (27).
For both spinsinglet and spintriplet central potentials, there may be a very weak attractive pocket of less than a few MeV at medium distance fm). However, considering the statistical and systematic errors, its existence should be carefully examined in future studies.
4.4 Scattering phase shifts and effective potentials
For quantitative studies of the interactions, it is desirable to calculate not only the potential but also scattering phase shifts, since the potential is not a physical observable as mentioned above. In this section, we therefore investigate a nature of interactions, by calculating scattering phase shifts from the obtained potentials. In particular, we study whether the LS potential of Fig. 1 leads to attractive behaviors in the scattering phase shifts in the channel .
We calculate the scattering phase shifts by solving the Schrödinger equation with the above potentials, parameterized with multiGaussian forms, with for the central and spinorbit potentials, whereas for the tensor potential to mimic the short distance behavior, as shown in Fig 1. Here, a scaling parameter fm is introduced to simplify the notation. The uncorrelated fits are performed reasonably. The resultant fit parameters and /dof are given in Table 1.
channel  [MeV]  [MeV]  [MeV]  /dof  

2173(268)  762(62)  236(65)  11(2)  2.1(0.3)  0.6(0.1)  1.7(0.8)  
421(122)  233(74)  397(16)  11.5(0.4)  1.3(0.2)  3.9(0.1)  1.1(1.0)  
711(11)  16(5)  —  2.6(0.2)  0.5(0.1)  —  0.8(0.5)  
45(17)  181(5)  315(12)  0.4(0.1)  1.4(0.2)  5.3(0.3)  3.6(0.7) 
The scattering observables are obtained from the long distance behaviors of linearly independent regular solutions, and are shown in Fig.2. The inner error is statistical, while the outer one is statistical and systematic combined in quadrature. Here, to estimate the systematic error, we take into account the uncertainty arising from the truncation of the derivative expansion and from the choice of fitting functions for the potentials. To estimate systematic errors associated with the truncation of the derivative expansion, we calculate phase shifts also at , and take differences of central values between and 7 as systematic errors. A dependence of phase shifts on a choice of fitting functions for the potentials is estimated by changing the fitting function to a Yukawatype. It turns out that the former dominates the systematic error except that the latter dominates in the channel. Although the magnitude of the phase shifts obtained from our potentials are smaller than the experimental ones, general trends are well reproduced except for the case at low energies. The missing attraction in the channel is likely due to the weak tensor force caused by the large pion mass. Among others, the most interesting feature in Fig.2 is the attraction in the 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 , , and channels in Fig.3, as defined in Eqs. (1)(3) and below. Indeed, one can see that has a weak repulsive core surrounded by an attractive well; the attraction is driven by the strongly attractive LS force, in Fig.1.
We give a comment on the reliable energy region of these phase shifts. Through the timedependent method, these phase shift are obtained based on the NBS wave functions at the energy points MeV with . Except for these energy points and the point for where the value of the phase shift vanishes by definition, the phase shifts are obtained by assuming that the derivative expansion is converged so that the truncated potentials in Eq. (6) do not depend on the energy. Because the contribution from MeV gradually dominates the intermediate states in as increases, the most reliable energy region of the phase shift is around MeV. Reliability for MeV can be explicitly examined by enlarging the spatial volume. Reliability for MeV can be examined by changing relative weight of excited states, which can be done either by studying the dependence of the truncated potentials or by changing the twonucleon source operator. Note that these arguments apply to any choices of interpolating fields.
5 Conclusion
We have made a first attempt to determine NN potentials up to NLO in the parityodd sector, which appears in the , , , – channels. Using CPPACS gauge configurations on a lattice ( fm and MeV), not only the central and the tensor potentials but also the spinorbit potential have been derived for the first time. These potentials are constructed from NBS wave functions for , 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 spinsinglet sector, the central potential is repulsive with a strong repulsive core at short distance. For the spintriplet sector, (i) the central potential is also repulsive with a repulsive core at short distance, (ii) the tensor potential is positive and quite weak, and (iii) the spinorbit potential is negative and strong at short distance.
We have then calculated scattering observables in the , , and – 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 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. This, together with the missing attraction in the channel, indicates an importance to carry out simulations at and around the physical quark mass.
Numerical calculations in this report are performed on the University of Tsukuba Supercomputersystem (T2K). This work is supported in part by the JSPS GrantinAid for Scientific Research (23540321, 24740144, 24740146), the JSPS GrantinAid for Scientific Research on Innovative Areas(Nos.2004: 20105001, 20105003) and SPIRE (Strategic Program for Innovative Research). We are also grateful for authors and maintainer of CPS++(28), a modified version of which is used for this work. We thank the CPPACS Collaboration and ILDG/JLDG for providing the 2flavor QCD gauge configurations(25); (29).
Footnotes
 journal: Physics Letters B
 See Appendix A in Ref.(23) 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 .
 The result does not depend on the order of multiplications, since these projection matrices are mutually commutative with each other.
References
 E. Epelbaum, H. W. Hammer, U. G. Meissner, Rev. Mod. Phys. 81 (2009) 17731825.
 R. Machleidt and D. R. Entem, Phys. Rept. 503 (2011) 1 [arXiv:1105.2919 [nuclth]].
 M. Lüscher, Nucl. Phys. B 354 (1991) 531.
 M. Fukugita, Y. Kuramashi, H. Mino, M. Okawa and A. Ukawa, Phys. Rev. Lett. 73 (1994) 2176 [heplat/9407012].
 S. R. Beane, K. Orginos and M. J. Savage, Int. J. Mod. Phys. E 17 (2008) 1157 [arXiv:0805.4629 [heplat]]
 S. R. Beane, W. Detmold, K. Orginos and M. J. Savage, Prog. Part. Nucl. Phys. 66, (2011) 1 [arXiv:1004.2935 [heplat]].
 T. Yamazaki et al. [PACSCS Collaboration], Phys. Rev. D 81, 111504 (2010) [arXiv:0912.1383 [heplat]].
 T. Yamazaki, K. i. Ishikawa, Y. Kuramashi and A. Ukawa, Phys. Rev. D 86, 074514 (2012) [arXiv:1207.4277 [heplat]].
 S. R. Beane, E. Chang, S. D. Cohen, W. Detmold, P. Junnarkar, H. W. Lin, T. C. Luu and K. Orginos et al., Phys. Rev. C 88, 024003 (2013) [arXiv:1301.5790 [heplat]].
 N. Ishii, S. Aoki and T. Hatsuda, Phys. Rev. Lett. 99 (2007) 022001.
 S. Aoki, T. Hatsuda and N. Ishii, Prog. Theor. Phys. 123 (2010) 89 [arXiv:0909.5585 [heplat]].
 N. Ishii et al. [HAL QCD Collaboration], Phys. Lett. B 712 (2012) 437 [arXiv:1203.3642 [heplat]].
 Reviewed in, S. Aoki et al. [HAL QCD Collaboration], PTEP 2012 (2012) 01A105 [arXiv:1206.5088 [heplat]].
 T. Kurth, N. Ishii, T. Doi, S. Aoki and T. Hatsuda, JHEP 1312, 015 (2013) [arXiv:1305.4462 [heplat]].
 A. Bohr and B. R. Mottelson, Nuclear Structure. (World Scientific Publishing Co., 1998). Sect. 25.
 R. Tamagaki, Prog. Theor. Phys. 44, 905 (1970). T. Takatsuka and R. Tamagaki, 46 (1971) 114. T. Takatsuka, 47, 1062 (1972), ibid. Prog. Theor. Phys. 48, 1517 (1972).
 M. Hoffberg, A.E. Glassgold, R.W. Richardson, M. Rudermann, Phys. Rev. Lett. 24 (1970) 775.
 M. Baldo, O. Elgaroey, L. Engvik, M. HjorthJensen and H. J. Schulze, Phys. Rev. C 58 (1998) 1921 [nuclth/9806097].
 D. Page, J. M. Lattimer, M. Prakash and A. W. Seiner, arXiv:1302.6626 [astroph.HE].
 C. J. D. Lin, G. Martinelli, C. T. Sachrajda and M. Testa, Nucl. Phys. B 619 (2001) 467 [heplat/0104006].
 S. Aoki et al. [CPPACS Collaboration], Phys. Rev. D 71 (2005) 094504 [heplat/0503025].
 N. Ishizuka, PoS LATTICE2009 (2009) 119.
 K. Murano, N. Ishii, S. Aoki, T. Hatsuda and , Prog. Theor. Phys. 125 (2011) 1225 [arXiv:1103.0619 [heplat]].
 M. C. Birse, arXiv:1208.4807 [nuclth].
 A. Ali Khan et al. [CPPACS Collaboration], Phys. Rev. D 65 (2002) 054505 [Erratumibid. D 67 (2003) 05990 ].
 A. X. ElKhadra, A. S. Kronfeld and P. B. Mackenzie, Phys. Rev. D 55 (1997) 3933 [heplat/9604004].
 R. B. Wiringa, V. G. J. Stoks and R. Schiavilla, Phys. Rev. C 51 (1995) 38.
 Columbia Physics System (CPS), http://qcdoc.phys.columbia.edu/cps.html

http://www.lqcd.org/ildg
http://www.jldg.org/  J. J. de Swart, C. P. F. Terheggen and V. G. J. Stoks, nuclth/9509032.