59

沼津工業高等専門学校研究報告
第 51 号，2017 年 1 月

Numerical evaluation of the gradient of the inter-quark potential
Yoshiaki KOMA∗
The inter-quark potential computed by lattice QCD Monte Carlo simulations contains valuable
information on the nonperturbative dynamics of quarks and gluons in quantum chromodynamics,
which is relevant to the understanding of the hadron structures. We investigate the functional formof the potential in detail not only by applying the χ2 ﬁtting but also by evaluating the gradient of
the potential, and clarify its behavior as a function of the inter-quark distance.
I.

INTRODUCTION

Recently, we have computed the potential energy for the various quark systems, the inter-quark
potential, by performing large scale numerical simulations of QCD using the discrete lattice formulation in
four-dimensional Euclidean space-time (lattice QCD
Monte Carlo simulations) [1–3]. QCD is an abbreviation of the quantum gauge ﬁeld theory associated with
the SU(3) group symmetry called “quantum chromodynamics” discovered in the 1970s [4], which is supposed to describe the dynamics of quarks and gluons
inside mesons and baryons, i.e., hadrons. We thenconfront with the problem of ﬁnding the appropriate functional form of the inter-quark potential, whichis relevant to the understanding of the hadron structures.
Given a set of experimental data, one often summarizes its behavior by ﬁtting it to a model function with adjustable parameters. The χ2 ﬁtting is
widely used for this purpose.

For the set of dis-crete data point (xi , yi ) (i = 1, 2, ..., np ) with thestandard deviation σi the model function yﬁt is assumed to be a function of x with a set of parameters�a = (a1 , a2 , ..., anpara ), which is determined so as tominimize the value of χ2 deﬁned by
2
np 

yi − yﬁt (xi ; �a)
.
χ2 =
σi
i=1

(1)

values associated with the statistical errors. However,a weakness of the χ2 ﬁtting may be that the model
function cannot be selected uniquely only by looking
at the value of χ2 . If there exists an analytical prediction from QCD on the behavior of the inter-quark
potential, the obtained functional form can be used
to examine the validity of QCD, but this is not thecase in QCD due to its nonperturbative nature, andthis is the reason why we have resorted to numericalsimulations.
In this report, we investigate the functional formof the inter-quark potential by applying not onlythe χ2 ﬁtting but also the additional method so as tocompensate the weakness of the χ2 ﬁtting. In particular, we look at the gradient of the inter-quark potential, that is, the derivative of the potential with respect
to the distance between quarks. We note, however,that this additional method is applicable only when
the numerical errors of the original data are highly
suppressed, otherwise the gradient will be noisy andrelevant information will be obscured. The numerical
data of our simulations satisfy this criterion by virtueof the simulation algorithm called the multi-level algorithm [5, 6].

II.

NUMERICAL PROCEDURES

Depending on the value of χ2 , one argues the goodness

We perform lattice QCD simulations within theof the ﬁtting and the validity of the model function.

quenched approximation in four dimensions with the

If the ﬁtting works well, the value of χ /ndf , wherelattice volume L3 ×T and the lattice spacing a by im-2

ndf = np − npara , will be of O(1).

posing periodic boundary conditions in all space-time

The χ ﬁtting can then be applied for the deter-

directions. The inter-quark potential is extracted frommination of the functional form of the inter-quark po-

the expectation value of the Polyakov loop correlation

tential, since the results are given by a set of numericalfunction (PLCF). In what follows we only explain how

2

to compute the inter-quark potential for the quarkantiquark system (mesons), but the similar procedure∗ Division of Liberal Arts / Physics groupis also applicable to the three-quark system (baryons).

60

沼津工業高等専門学校研究報告

We deﬁne a two-link correlator as

第 51 号

relatively large sub-lattice averages of its components

T(t, �r1 , �r2 )αβγδ ≡ U4 (t, �r1 )αβ U4 (t, �r2 )∗γδ ,(in our case it corresponds to the product of Ts within(2)

a sub-lattice), where the sub-lattice is deﬁned by di-which is a direct product of two temporal link vari-viding the lattice volume into several layers along theables at a time t separated by a distance r = |�r1−�r2 | andtime direction. During the computation of the sub-consists of a complex matrix with 3 = 81 components

lattice averages, the spatial links at the sub-lattice(Greek indices take the values from 1 to 3). The spa-boundaries are kept intact. The computation of thetial vectors �r1 and �r2 represent the spatial positions ofcorrelation function in this way is supported by thea quark and an antiquark, respectively. The two-link

transfer matrix formalism of quantum ﬁeld theory.

4

correlator acts on the color states in the 3 ⊗ 3̄ representation of the SU(3) group |n; �r1 , �r2 �αβ , which satis-

ﬁes T(t, �r1 , �r2 )αλγ� |n; �r1 , �r2 �αγ = e−En (r)aIII.

|n; �r1 , �r2 �λ� ,We carried out lattice QCD simulations using thewhere n is the principal quantum number (repeatedstandard Wilson gauge action at the coupling β = 6.00

indices are assumed to be summed over from 1 to 3).

on the L3 × T = 243 × 24 lattice with the multi-level

The energies En (r) are positive, which are common

algorithm. The lattice spacing a = 0.093 [fm], which

to all color components of |n; �r1 , �r2 �αβ . The multi-

is determined by the Sommer scale r0 = 0.50 [fm] [7].

plication rule of the two-link correlators for the two

One Monte Carlo update consisted of one heat-bathneighboring time at t and t + a is{T(t, �r1 , �r2 )T(t + a, �r1 , �r2 )}αβγδ

= T(t, �r1 , �r2 )αλγ� T(t + a, �r1 , �r2 )λβ�δ .

NUMERICAL RESULTS

and ﬁve over-relaxation steps. The number of sublattice Nsub = 6, which corresponds to Ntsl = 4, sothat aNtsl = 0.372 [fm]. The number of internal up-(3)

dates Niupd = 10000 and the number of gauge con-

By using the two-link correlators, the PLCF is

ﬁgurations Ncnf = 200. The numerical errors of thethen constructed asinter-quark potential and its derivatives are estimatedby the jackknife method.

TrP (�r1 )TrP (�r2 )∗In Fig. 1, we plot the quark-antiquark potential= {T(0, �r1,�r2 )T(a, �r1,�r2 )· · ·T(T −a, �r1,�r2 )}ααγγ . (4)as a function of the inter-quark distance r/a. We

The expectation value is evaluated by inserting theﬁnd that the most of the data fall into one curve.

complete set of eigenstates at all times t = 0, a, ..., T −The labels (i, j, k) represent the base vectors along

a as

∞


the quark-antiquark direction. For instance, the data
(5)

belong to (1, 0, 0) are for the case that the quark andwhere w0 = 1 is guaranteed automatically. We ﬁnallyaxis), while the data belong to (1, 1, 0), (1, 1, 1), andobtain the ground state potential V (r) ≡ E0 (r) by(2, 1, 0) are for the cases that the quark and antiquark�TrP (�r1 )TrP (�r2 )∗ � =V (r) = −wn e−En (r)T ,antiquark are separated along the x, y, and z axes (on-n=0

are separated along the oﬀ-axis. The unit Euclidean√
distances for these oﬀ-axis settings are r/a = 2,√√
3, and 5, respectively. Note that the tree-level1ln�TrP (�r1 )TrP (�r2 )∗ �+O(e−(E1 −E0 )T ), (6)Twhere the terms of O(e−(E1 −E0 )T ) are always negli-

improvement is applied to the inter-quark distance be-gible at zero temperature. Although it is impossible

fore plotting the data by using the three-dimensional

to compute the expectation values of the PLCF accu-

Green’s function in inﬁnite volume, which helps to re-

rately within the ordinary simulations as they are ex-duce the lattice discretization eﬀect [7]. The dotted

tremely small at long distances, it is possible to over-line is the model function in the form

come this problem by using the multi-level algorithm.
The idea of the multi-level algorithm is to com-c
Vﬁt (r) = − + σr + μ ,rpute the correlation function, which may have an ex-(7)which is ﬁtted to the on-axis (1, 0, 0) data. We obtain

tremely small expectation value, from the product of2

61(1,0,0)(1,1,0)(1,1,1)(2,1,0)2
c/r + σa F(r/a)aV(r/a)Numerical evaluation of the gradient of the inter-quark potential

2(1,0,0)(1,1,0)(1,1,1)(2,1,0)fit (1,0,0)

r/a

r/a

FIG. 1: The quark-antiquark potential as a function of theinter-quark distance r/a.

FIG. 2: The ﬁrst derivative of the quark-antiquark potential with respect to the inter-quark distance r/a.

c = 0.2826(7), σa2 = 0.04558(7), and μa = 0.7540(5)Therefore we expect

for the ﬁtting range r/a = 1.855 ∼ 9.974 (correspond-

V (r + h) + V (r − h) − 2V (r)
h2
2
h
h4
= V �� + 2 f (4) + 2 f (6) + · · ·
4!
6!
2c
h2
h4
= − 3 (1 + 2 + 4 + · · · )
r
r
r
2c
h2 −1
= − 3 (1 − 2 ) ,rr

ing to r/a = 2 ∼ 10 before the tree-level improvement)with np = 9 and χ2 = 3.1 (χ2 /ndf = 0.52). We call c

the Coulombic coeﬃcient, σ the string tension, and μ
just a constant. We then ﬁnd that the oﬀ-axis data
behave nicely up to the distances longer than the halfof the spatial lattice volume r/a = 12, at where theon-axis data already suﬀer from ﬁnite volume eﬀects.

and obtain

In order to conﬁrm the validity of this functionalc=−form, we further investigate the gradient of the potential data. In Fig. 2, we plot the ﬁrst derivative

Eq. (7) isinter-quark distance r/a,

c=−

(8)

r3 ��
V (r) .
2 ﬁt

(13)Note that the O(h2 ) terms present in the ordinary

where h/a = 1 for the on-axis (1, 0, 0) data, and h/a =√ √ √2, 3, 5 for the oﬀ-axis (1, 1, 0), (1, 1, 1), (2, 1, 0)

second diﬀerence is properly removed by the factor
2(1 − hr2 ) in Eq. (12).

data, respectively. If the model function in Eq. (7) is

In Fig. 3, we plot the Coulombic coeﬃcient cvalid, the ﬁrst derivative F should behave ascdVﬁt= 2 +σ .
drrr3
h2 V (r + h)+V (r − h)−2V (r)) . (12)(1− 2 )(2rh2

The corresponding continuum expression derived fromof the quark-antiquark potential with respect to theV (r + h) − V (r)h,
F (r + ) =
2
h

(11)

computed by Eq. (12) as a function of r/a.

We

ﬁnd that the value of c is approximately constant(9)

around c � 0.3. The two solid lines correspond to thec = 0.282 and c = 0.315, which are the minimum and

The data seem to fall into this curve (dotted line) andapproach the value σa = 0.04558 at long distances

the maximum values obtained by the χ2 ﬁtting of the(solid line).

potential for the on-axis and the oﬀ-axis settings. The

2values of c in the range r/a = 3 ∼ 8 lie within these

It is possible to extract the Coulombic coeﬃcient c directly from the second derivative of the po-

two lines. The deviation at long distances r/a > 9 is

tential. For the model function in Eq. (7) its deriva-

due to ﬁnite volume eﬀects, which is expected to dis-

tives with respect to r behave asappear for larger lattice volumes. If we consider the

2c
4!c
6!c
(4)
(6)
��
= − 3 , Vﬁt = − 5 , Vﬁt = − 7 , ... .
Vﬁt
r

r

r

deviation at short distance r/a < 3 seriously, we may(10)

conclude that the term −c/r in Eq. (7) is not enough3

62

沼津工業高等専門学校研究報告

第 51 号

Eq. (7) is(1,0,0)(1,1,0)(1,1,1)(2,1,0)

1 d2(rVﬁt (r)) .
2 dr2

(15)

c

σa2 =

We ﬁnd that the values of σa2 behave as constant,which lie on the solid line corresponding to the value
σa2 = 0.04558 obtained by the χ2 ﬁtting.

Wethen conclude that the linear-rising term is certainlypresent in the quark-antiquark potential. The deviar/ation at long distances is due to ﬁnite volume eﬀects,
while the deviation at short distances indicates that

FIG. 3: The Coulombic coeﬃcient c extracted from thesecond derivative of the potential.

there exist the terms which do not behave as constanteven after multiplying r to the potential.

(1,0,0)(1,1,0)(1,1,1)(2,1,0)SUMMARY

σa

2

IV.

We have investigated the functional form of theinter-quark potential obtained by the lattice QCD
simulations not only by applying the χ2 ﬁtting but
also by evaluating the gradient of the potential. We
have conﬁrmed that the quark-antiquark potential can

r/a

be described mostly by the sum of the Coulombic term

FIG. 4: The string tension σa2 extracted from the second
derivative of the potential after multiplying r.

and the linear term, which is actually the well-known
result, but we have further clariﬁed the r dependence

to describe the short distance behavior of the poten-of the Coulombic coeﬃcient and the string tension.

tial. A possible explanation may be that the c is not

The results indicate that some modiﬁcations to theconstant but dependent on r. It is interesting to in-

Coulombic term are needed for the perfect description

vestigate the short distance behavior in detail, whichof the potential including the data at short distances.

however is beyond the scope of the present report.

We plan to apply similar analyses presented in thisreport to the determination of the functional form of

In Fig. 4, we plot the string tension σa2 directly

the three-quark potential [1–3].

extracted from the second derivative of the potential

The author is grateful to Miho Koma (Nihon

after multiplying r,

University) for fruitful collaboration. The computa-

1 (r + h)V (r + h)+(r − h)V (r − h) − 2rV (r)
.
σa =
2
h2
(14)

tion was performed on the supercomputer NEC SX-

2

ACE at Cybermedia Center, Osaka University. This
work is partially supported by JSPS KAKENHI Grant

The corresponding continuum expression derived fromNumber 24740176.

[1] Y. Koma and M. Koma, PoS LAT2013 (2013) 469.
[2] Y. Koma and M. Koma, PoS LAT2014 (2014) 352.
[3] M. Koma and Y. Koma, PoS LAT2015 (2015) 324.
[4] W. Greiner and A. Schäfer, “Quantum chromodynamics,” Berlin, Germany: Springer (1994).

[5] M. Lüscher and P. Weisz, JHEP 09 (2001) 010.
[6] M. Lüscher and P. Weisz, JHEP 07 (2002) 049.
[7] S. Necco and R. Sommer, Nucl. Phys. B622 (2002) 328.

4