41

沼津工業高等専門学校研究報告
第 57 号，2023 年 2 月

Numerical Simulation of Compact U(1) Lattice Gauge Theory
KOMA Yoshiaki ∗Abstract: Compact U(1) lattice gauge theory possesses a conﬁnement phase, which can be regarded
as a prototype of quantum chromodynamics QCD for describing strong interaction of quarks andgluons. We revisit this theory and perform Monte Carlo simulations to gain idea on the structureof the nonperturbative QCD vacuum.
Key Words: Compact U(1) lattice gauge theory, Lattice QCD simulations

1.

Introduction

Clarifying nonperturbative properties of quantum chromodynamics (QCD), a gauge ﬁeld theory ofSU(3) group symmetry for describing strong interaction of quarks and gluons, remains an important topic
in hadron physics. Monte Carlo simulations of QCD
deﬁned on a discrete spacetime, referred to as the lattice QCD simulations, oﬀer a powerful tool for thispurpose, which in principle make it possible for us tocompute directly the expectation values of physicalquantities involving nontrivial quantum eﬀects. The
problems are, however, that the numerical simulations
usually require a lot of computer resources for reducing numerical errors, and the numerical result itself
is not enough to understand the underlying physicalmechanism, even if one successfully obtains very precise numerical values. In this respect, complementarystudies are always useful to get an insight into the true
nature of nonperturbative QCD.
One of the mysterious phenomena that could be
explained by nonperturbative QCD is the fact that one
never observe an isolated quark in the accelerator experiment, which means that the quarks are somehow
conﬁned inside hadrons. This is called the quark conﬁnement problem, and is not yet satisfactorily solved.
For this problem the lattice QCD simulations provide
us with a numerical answer. The static potential extracted from the expectation value of the Wilson loop
shows that the functional form of the potential consists of a Coulombic term and a linearly rising term
as a function of the distance between the quark andthe antiquark. Although what kind of mechanism determines the functional form is not clariﬁed, the linear
term can explain why a single quark cannot be isolatedwith a ﬁnite amount of energy.
A plausible idea on the mechanism is that the
QCD vacuum acts as a dual superconductor for thecolor-electric ﬁeld associated with the quarks. In fact,
the static potential of the classical ﬂux-tube solution
in the dual Ginzburg-Landau (DGL) theory exhibits asimilar functional form [1], where the physical mechanism in this eﬀective theory originates from the dual
Meissner eﬀect. The static potential for the threequarks that we have successfully obtained by the lat-

∗ Division of Liberal Arts

tice QCD simulations [2] also supports the dual superconductor picture of the QCD vacuum. However,
it is as yet nontrivial why QCD has a property like adual superconductor. If this is the true nature of nonperturbative QCD, the theory should be connected
theoretically with the DGL theory.
In this report, we revisit numerical simulations of
the compact U(1) lattice gauge theory (U(1) LGT) [3,4]. This theory is not exactly the same as QCD, but
can be regarded as a small laboratory of QCD as it
possesses a conﬁnement phase. In addition, the reduction of gauge degrees of freedom drives down the costof simulations, allowing for computing various quantities. In the conﬁnement phase, many direct evidenceof supporting the dual superconductor picture havebeen found so far. The observed proﬁle of the electricﬁeld between a positive and a negative electric charges
has a one-dimensional ﬂux tube shape, accompaniedby circulating monopole supercurrent [5], which meansthat the dual Meissner eﬀect clearly occurs. Therefore, detailed comparisons with the U(1) LGT resultswith those of the DGL theory from various points ofview may give us some hints for the theoretical connection between QCD and the DGL theory.

2.

Numerical procedure

Given an action of the theory S[φ] of ﬁeld variables φ, the expectation value of a quantum operator
O(φ) is schematically computed as�O� = Dφ O(φ) e−βS[φ] ,(1)
where the value β controls the strength of interaction, corresponding to the inverse temperature in statistical mechanics. In practice, we deﬁne the theory
on a discrete lattice of spacetime and generate a sequence of independent ﬁeld conﬁgurations {φi } with
the Boltzmann weight of exp(−βS[φ]) by the MonteCarlo method. We then evaluate the average of theoperator for the number of conﬁgurations Nconf , which
ﬁnally reduces to the expectation value �O� by taking
the limit Nconf → ∞. This is the basic methodologyof the lattice simulations. Of course, it is impossible
to prepare inﬁnite number of the ﬁeld conﬁgurations
with inﬁnite lattice volume with inﬁnitesimal latticespacing, and hence statistical and systematic errors of
the numerical results need to be controlled properly.

42

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

In the U(1) LGT, the action consists of link variables (variables between two nearest sites) deﬁned byUμ (m) = exp(iθμ (m)) ∈ U(1), where μ and m denotea direction and a site labels. The link variables arecompact as they are unchanged under the transformation θμ → θμ + 2πnμ with an integer nμ . The most
simple action of the U(1) LGT is called the Wilson
gauge action of the formSμν (m) ={1 − Re Uμν (m)} , (2)S[U ] =m, μ<ν

1.0
E

0.8

E

0.6
0.4
0.2

m, μ<ν

0.0
0.0

where Uμν (m) = Uμ (m)Uν (m + μ̂)Uμ† (m + ν̂)Uν† (m)are called the plaquette variables. We consider four
dimensional spacetime with the lattice volume V =L3 T , where L and T are the space and time extentsof the lattice, and impose periodic boundary conditions in all directions. The number of link variables is
Nlink = 4L3 T . A sequence of independent link variables is generated by using a mixture of heatbath (HB)
and over-relaxation (OR) methods. In the weak coupling limit as β → 0, the action reduces to the continuum action for the ﬁeld strength part of quantum
electrodynamics (QED).
3.

第 57 号

1.0

1.5

2.0

β

Fig. 1: The average action per plaquette E as a functionof β values on a 164 lattice.

1.0

ACR

0.8

0.6

Numerical results

β = 0.80
β = 1.00
β = 1.20

0.4

In order to test the validity of our OpenMP codewith double precision arithmetic, we ﬁrst compute theaverage action per plaquette �E� = �S�/Nplaq (expectation value of the action density) on a 164 lattice asin the works [3, 4], where an update of link variables
is achieved by 1 HB followed by 3 OR. The numberof plaquette variables is Nplaq = 6L3 T . In Fig. 1, weplot the result as a function of the β values. Each data
point is given by the average of Nconf = 50 conﬁgurations, which are separated by 100 updates. We obtainthe same result as in [3, 4], supporting the validityof our setup. The present result is obtained by using the Mac mini (2018). The phase transition occurs
around β � 1.01. In order to achieve an eﬃcient link
update, we have checked the acceptance ratio (ACR)
in advance as a function of the trial numbers Nhit as
in Fig. 2, and then decided to take Nhit = 4, whichguarantees that about 95% of whole link variables are
replaced by the new ones in 1 HB process.
We next examine the probability density function (PDF) of the action current density and its relation to the monopole current [6]. We construct theAbelian plaquette variables from the phase of the linkvariables as

0

2

4

6

8

10

Nhit

Fig. 2: The acceptance ratio (ACR) of link variables (average of 80 updates) as a function of the trial numbers Nhitat β = 0.80, 1.00, and 1.20 on a 164 lattice.

where md = m + (1̂ + 2̂ + 3̂ + 4̂)/2 denotes a dual site.
We deﬁne the action current density on the dual links
so as to be comparable with the location of monopole
current such as1|�μαβγ |(Sαβ (m+μ̂)+Sαβ (m+μ̂+γ̂)) . (5)12Note that this satisﬁes � md ,μ Sμ (md )�/Nlink = �E�.
In Fig. 3, we plot the normal PDF and the PDF juston the monopole current on a 164 lattice, where three
β values at 0.80, 1.00 (conﬁnement phase), and 1.20
(deconﬁnement phase) are selected. We ﬁnd an interesting behavior in the conﬁnement phase such that thenormal PDF is distorted by the presence of monopole
current. A remarkable point may be that the peak
position of the PDF on the monopole current is insensitive to the change of β values, although the heightdecreases as the β increases. There is a small peakeven at β = 1.20 in the deconﬁnement phase, although
its contribution to the PDF is negligible.
We then compute the static potential between
a positive and a negative electric charges. A usualSμ (md ) =θμν (m) = θμ (m)+θν (m+ μ̂)−θμ (m+ ν̂)−θν (m) , (3)and decompose this into a regular part θ̄μν ∈ [−π, π)
and a singular part nμν (m) = 0, ±1, ±2 as θμν (m) =θ̄μν + 2πnμν (m). The boundary of nμν (m) is thenidentiﬁed as the monopole current
1
kμ (md ) = − �μνρσ ∂ν nρσ (m + μ̂) ,
2

0.5

(4)
2

43

Numerical Simulation of Compact U(1) Lattice Gauge Theory42.0
β = 0.80

normalon monopole

1.5
aV (r/a)

PDF

3

2

Vﬁt (r)1.0(1, 0, 0)

1

0.5

0
0.00

0.0

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

0.25

0.50

0.75

1.00

1.25

1.50

0

2

4

Sμ (md )
4
normalon monopole

310

Fﬁt (r)(1, 0, 0)0.4(1, 1, 0)
a2 F (r/a)

PDF

8

0.5
β = 1.00

2

1

0.3

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

0.2
0.1

0
0.00

0.25

0.50

0.75

1.00

1.25

0.0

1.50

0

Sμ (md )

2

4

6

8

10

r/a

Fig. 4: The static potential (upper) and the force (lower)
at β = 1.00 on a 244 lattice for Ntsl = 3 and Niupd = 2000
with Nconf = 20. The labels (1, 0, 0), (1, 1, 0), (1, 1, 1), and(2, 1, 0) specify the unit vectors pointing from a positive
to a negative charges.

4
β = 1.20

normalon monopole

3
PDF

6
r/a

2

In Ref. [5], we computed the PLCF by using themultilevel algorithm [7] and obtained precise results
on a 164 lattice. The multilevel algorithm is applicable to any lattice gauge theories as long as a local
gauge action, such as the Wilson action, is simulated.
In this algorithm, the lattice volume is divided into
sublattices along the time direction and the parts ofoperator are simulated with internal updates of linkvariables, where the spatial links at the boundariesof the sublattices remain intact. The eﬃcient use of
the algorithm requires a tune of the temporal size ofthe sublattice Ntsl and the number of internal updates
Niupd . In Ref. [5], we performed simulation with thetemporal size of the sublattice Ntsl = 2 within single
precision arithmetic.
In the present study, we have examined the temporal size dependence of the numerical errors on thepotential and the force at β = 1.00 on a 244 latticewith Ntsl = 2, 3, and 4, where the static charges are
separated by the distances along not only on axis but
also oﬀ axis. Note that the force is computed from thepotential by

1

0
0.00

0.25

0.50

0.75

1.00

1.25

1.50

Sμ (md )

Fig. 3: The probability density function (PDF) of the action current density at β = 0.80 (top), 1.00 (middle), and1.20 (bottom) on a 164 lattice. The circle symbols are forthe normal action current density, while the square just on
the monopole current.

method is to extract it from the expectation valueof the Wilson loop, but here, we extract it from thecorrelation function of the Polyakov loop operator
T� t), which we call the PLCF, as
P (m)� = t=1 U4 (m,V (r) = −1
ln�P (m� 1 )P (m� 2 )� + O(e−(ΔE)T ) , (6)T� 2 |. It is important to note that thewhere r = |m
� 1 −msystematic error is much smaller than that from the
Wilson loop once the expectation value of the PLCF
is computed accurately.

F (r +3h
V (r + h) − V (r)
)=
,
2
h

(7)

44

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

第 57 号0.4
σa2 = 0.17
(1, 0, 0)

0.3

Fﬁt (r) =gives the values c = 0.2699(52) and σ = 0.1713(20)with χ2 /Ndf = 3.2, where the data in the range
r/a = [1.5, 5.5] are taken into account. The slopeof the linear term σ, which we call the string tension,from Eq. (8) agrees with that from Eq. (9). The small
discrepancy in the Coulombic coeﬃcient may reﬂect asystematic discretization eﬀect.
If the functional form of Eq. (8) is appropriate
to describe the behavior of the potential, the string
tension can also be extracted by evaluating(1, 1, 1)
σa2(9)(1, 1, 0)0.2(2, 1, 0)

0.1

0.0

c
+σ ,
r2

0

2

4

6

8

10

r/a

Fig. 5: The direct extraction of the string tension.

σ=
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)
data, respectively.
Our limited observation suggests that Ntsl = 3
is an optimal choice at β = 1.00 among the three.
The results with Ntsl = 3 are then plotted in Fig. 4,
where the number of internal updates is set to be
Niupd = 2000 with the conﬁgurations Nconf = 20. The
standard jack-knife method is used to evaluate the statistical errors. We used the Mac Pro (Late2013) in thiscomputation. In fact, the behavior of the potentials
with Ntsl = 2 and with Ntsl = 3 are not so diﬀerent
with each other, but the force exhibits a better behavior with Ntsl = 3. If one is interested in the long
distance behavior of the observables, of course, thenumber of Niupd should be increased. We ﬁnd that
there is no large discrepancy between the data fromthe on-axis and from the oﬀ-axis data, which meansthat the rotational invariance is reasonably restoredat β = 1.00.
The ﬁtting of the potential for the on-axis data to
the functional form of a Coulomb and a linear terms,

1 d2(rV (r)) ,
2 dr2

(10)

where we take the diﬀerence instead of the derivativeon the lattice. The result is plotted in Fig. 5. Although the repetition of taking diﬀerence of the data
increases the numerical error, we observe a plateau
around σa2 ∼ 0.17, which agrees with those from thepotential in Eq. (8) and the force in Eq. (9).

4.

Summary

gives the values c = 0.2363(71), σ = 0.1736(27), andμ = 0.5211(96) with χ2 /Ndf = 1.6, where the data in
the range r/a = [1, 6] are taken into account. Similarly, the ﬁtting of the force for the on-axis data to
the functional form (derivative of Eq. (8) with respectto r),We have performed numerical simulation of theU(1) LGT, and have investigated some of basic physical quantities involving the action and the static
potential in order to test our numerical code with
OpenMP. The peculiar structure of the PDF of the action current density that we have found in the conﬁnement phase indicates that the monopoles surely play
an important role for the conﬁnement mechanism. In
order to make a quantitative comparison between theU(1) LGT and the DGL theory, we further need toaccumulate numerical data. It is also important todiscriminate the lattice discretization eﬀect at shortdistances. The use of the lattice Green function inthree dimensions will be helpful for this purpose. As
demonstrated in [8] there is a possibility that the dynamics of monopoles can be transformed into a scalarﬁeld theory like the Abelian Higgs model. The lattice
results will be useful for the quantitative evaluation ofsuch an idea.
The author is grateful to Miho Koma for fruitfulcollaboration.

[1] Y. Koma, E. M. Ilgenfritz, T. Suzuki and H. Toki,
Phys. Rev. D 64, 014015 (2001)[2] Y. Koma and M. Koma, Phys. Rev. D 95, 094513(2017)[3] M. Creutz, L. Jacobs and C. Rebbi, Phys. Rev. D 20,1915 (1979)[4] K. J. M. Moriarty, Phys. Rev. D 25, 2185 (1982)[5] Y. Koma, M. Koma and P. Majumdar, Nucl. Phys. B
692, 209 (2004)[6] T. A. DeGrand and D. Toussaint, Phys. Rev. D 22,2478 (1980)[7] M. Lüscher, P. Weisz, JHEP 09, 010 (2001)[8] K. Bardakci and S. Samuel, Phys. Rev. D 18 (1978),2849

cVﬁt (r) = − + σr + μ ,r(8)4