The complex Langevin method is extended to full QCD at non-zero chemical
potential. The use of gauge cooling stabilizes the simulations at small enough
lattice spacings. At large fermion mass the results are compared to the HQCD
approach, in which the spatial hoppings of fermionic variables are neglected,
and good agreement is found. The method allows simulations also at high
densities, all the way up to saturation.Comment: 5 pages, 5 figures, PLB version, minor change

Sexty, Denes

Physics Letters B

English

arXiv

AbstractThe complex Langevin method is extended to full QCD at non-zero chemical potential. The use of gauge cooling stabilizes the simulations at small enough lattice spacings. At large fermion mass the results are compared to the HQCD approach, in which the spatial hoppings of fermionic variables are neglected, and good agreement is found. The method allows simulations also at high densities, all the way up to saturation

Sexty, Dénes

Elsevier - Publisher Connector 

Simulating full QCD at nonzero density using the complex Langevin equation 

Dénes Sexty

Crossref

Simulating full QCD at nonzero density using The Complex Langevin Equation

Physics Letters B 729 (2014) 108–111Contents lists available at ScienceDirect
Physics Letters B
www.elsevier.com/locate/physletb
Simulating full QCD at nonzero density using the complex Langevin
equation
Dénes Sexty
Institut für Theoretische Physik, Universität Heidelberg, Germany
a r t i c l e i n f o a b s t r a c t
Article history:
Received 14 October 2013
Received in revised form 7 January 2014
Accepted 8 January 2014
Available online 10 January 2014
Editor: J.-P. Blaizot
The complex Langevin method is extended to full QCD at non-zero chemical potential. The use of gauge
cooling stabilizes the simulations at small enough lattice spacings. At large fermion mass the results are
compared to the HQCD approach, in which the spatial hoppings of fermionic variables are neglected,
and good agreement is found. The method allows simulations also at high densities, all the way up to
saturation.
© 2014 The Author. Published by Elsevier B.V. This is an open access article under the CC BY license
(http://creativecommons.org/licenses/by/3.0/). Funded by SCOAP3.The determination of the phase diagram of ﬁnite density QCD
is one of the great problems of theoretical physics today. One is
interested in averages deﬁned with the Euclidean path integral
〈
f [U ]〉= 1
Z
∫
DUe−Sg [U ] detM(μ,U ) f [U ], (1)
where Sg[U ] is the Yang–Mills action of the gauge ﬁelds and
M(μ,U ) is the Dirac-matrix of the quark ﬁelds. Naive lattice simu-
lations at μ = 0 using importance sampling are made unfeasible by
the fact that the determinant of the fermion matrix is a complex
number in general. Various methods have been invented to cir-
cumvent the problem, but these are of limited use [1], mostly be-
ing applicable for μ/T  1. An exception is the complex Langevin
method [2], which is not limited to small chemical potential. It
has been demonstrated that this method allows for the solution of
the sign problem in various systems [3–7], but in some cases also
non-physical results are delivered [8–12]. In this Letter I demon-
strate that the algorithm can be extended to full QCD with light
quark masses on lattices with suﬃciently small lattice spacings.
The complex Langevin method is based on setting up a com-
plex Langevin equation (CLE) in an enlarged manifold, which is the
complexiﬁcation of the original ﬁeld space [2]. The original the-
ory is recovered by taking expectation values of the analytically
continued observables. For SU(N) gauge theories this complexiﬁ-
cation is SL(N,C). This method can also be applied to other cases
where the action becomes complex, e.g. the case of real time
evolution, where the complexity of the action is much ‘larger’, us-
ing the Minkowskian formulation of the path integral [13–15], or
Yang–Mills theory with Θ-term [16]. In this work I am concerned
E-mail address: d.sexty@thphys.uni-heidelberg.de.http://dx.doi.org/10.1016/j.physletb.2014.01.019
0370-2693/© 2014 The Author. Published by Elsevier B.V. This is an open access article
SCOAP3.with ﬁnite density physics, where the complexity of the action is
present at non-zero chemical potential. The analytic understanding
of the breakdowns and successes of the complex Langevin method
has improved in the last few years [17–21], one can gain an in-
sight whether the results are trustworthy using requirements such
as the fast decay of the distributions.
Recently an important breakthrough in this ﬁeld was the de-
velopment of a ‘gauge cooling’ algorithm for the CLE method [7],
where the gauge symmetry of the system is used to ensure a well
localized distribution in the complexiﬁed ﬁeld space, and thus con-
vergence to the correct results.
In this work the CLE method is applied to the lattice discretiza-
tion of full QCD, i.e. for the action
Seff [U ] = Sg[U ] − NF4 lndetM(μ,U ) (2)
where Sg[U ] is the Wilson plaquette action for the SU(3) link
variables, and M(μ,U ) is the unimproved staggered fermion de-
terminant for NF fermion ﬂavors
M(μ,U )xy =mδxy +
∑
ν
ην(x)
2a
[
eδν4μUν(x)δx+aν ,y
− e−δν4μU−1ν (x− aν)δx−aν ,y
]
, (3)
where x and y indices represent spacetime coordinates, and ημ(x)
are the staggered sign functions. Periodic (antiperiodic) boundary
conditions are used in space (time) directions. The fermion matrix
fulﬁlls the symmetry condition:
xM(μ,U )xyy = M†
(−μ∗,U)yx (4)
with the “staggered γ5 matrix”, x = (−1)x1+x2+x3+x4 . This symme-
try leads to detM(−μ∗,U ) = (detM(μ,U ))∗ . This means that the
determinant becomes complex for Reμ = 0, making a simulationunder the CC BY license (http://creativecommons.org/licenses/by/3.0/). Funded by
D. Sexty / Physics Letters B 729 (2014) 108–111 109based on importance sampling impossible. Without rooting (i.e. us-
ing NF < 4 by taking a root of the fermion determinant in the path
integral), the staggered determinant describes 4 tastes of fermions.
In the Langevin dynamics (see below) the implementation of any
(not necessarily integer) number of ﬂavors is trivial, NF appears as
a factor of a drift term. In this study I have chosen to use NF = 4
and NF = 2.
The Langevin equation for the link variables is set up using the
equation
Ux,ν (τ + ) = Rx,ν(τ )Ux,ν (τ ), (5)
with
Rx,ν(τ ) = exp
[
i
∑
a
λa(Kaxν +
√
ηaxν)
]
. (6)
Here λa are the generators of the gauge group, i.e. the Gell-Mann
matrices. The drift force is determined by
Kaxν = −Daxν Seff [U ] (7)
with the left derivative
Daxν f (U ) = ∂α f
(
eiαλaUx,ν
)∣∣
α=0. (8)
The drift term for the action (2) is written as
Kaxν = −Daxν Sg[U ] + NF
4
Tr
[
M−1(μ,U )DaxνM(μ,U )
]
. (9)
It has been suggested in [22] that the drift corresponding to the
action (2) might also include a term reﬂecting the branch cut of
the complex logarithm on the negative real axis. This yet unclari-
ﬁed issue is the subject of ongoing research.
The drift term remains real only for μ = 0. Since the explicit
calculation of the inverse of the fermion matrix is quite costly, this
naive algorithm is feasible only for small lattice sizes. As a cost
effective alternative, the bilinear noise scheme [23,24], which is
related to pseudofermionic variables, is introduced as follows. The
drift of the link variables is calculated using
Kaxν = −Daxν Sg[U ] + NF
4
η+M−1DaxνMη, (10)
where the η is a vector of Gaussian random numbers satisfying
〈η∗xηy〉 = δxy . For the calculation of the drift term one has to solve
the linear system of equation M+ψ = η. In terms of the solution ψ
the drift term is written as
Kaxν = −Daxν Sg[U ] + NF
4
ψ+DaxνMη. (11)
This means that this algorithm uses the conjugate gradient (CG)
algorithm once for every timestep for the solution of the linear
system. In [24] it was also shown that with a higher order algo-
rithm one can get rid of part of the O (2) corrections from the
Fokker–Planck equation.
From previous studies of the complex Langevin equation one
learns the heuristic approach that a well localized distribution of
the variables in the complexiﬁed ﬁeld space is desirable. A useful
measure of the size of the distribution in imaginary directions of a
link variable is the unitarity norm
Tr
((
UU+ − 1)2) 0, (12)
where the equality is reached only for SU(N) matrices. The en-
larged gauge symmetry of the system can be used to decrease the
unitary norm of the system, thus ensure convergence to the ex-
act results [15]. Recently, we developed and tested in HQCD (see
below) a procedure utilizing this freedom called gauge cooling [7]Fig. 1. Unitarity norm as a function of Langevin time with and without cooling for
several values of the Langevin timestep  .
(reminiscent of stochastic gauge ﬁxing [25]). The idea is the fol-
lowing: one uses gauge transformations
Ux,ν → Ω(x)Ux,νΩ−1(x+ aν) (13)
with Ω(x) ∈ SL(N,C) to decrease the unitarity norm of the system.
This can be accomplished by choosing the Ω(x) matrices in the di-
rection of the steepest descent of the unitarity norm. Advanced
versions of the algorithm ensuring faster decay of the unitarity
norm use adaptive stepsize and Fourier acceleration [26].
The consequence of using the bilinear noise scheme is that an
imaginary part of the drift term is generated already at zero μ,
as the drift is real only on the average. Using a smaller Langevin
step allows the system to better approximate the drift term within
a given Langevin time-window, therefore the resulting equilibrium
unitarity norm of the simulation should vanish in the zero  limit,
see Fig. 1. Without gauge cooling the non-unitarities generated by
the noise term (or rounding errors in the case of the exact inverse
algorithm) would grow exponentially, breaking down the simula-
tion. Generally, we ﬁnd also at μ = 0, using suﬃcient cooling,
that the level of unitarity norm stabilizes and allows one to ob-
tain correct results (after  → 0 extrapolation) for lattices with
ﬁne enough lattice spacings. As one observes, at low β (below
β ≈ 5.0–5.3 for NF = 4) the cooling is not effective enough to pre-
vent the system from wandering off far from the SU(3) manifold,
and ‘skirted’ distributions develop (as also observed in [7]). Close
to the continuum limit, however, the algorithm seems to be sta-
ble irrespective of the physical phase, as observed using cheaper
HQCD simulations.
Observables are measured on ‘slices’ of the T , μ phase diagram
(meaning a scan using one variable while keeping the other ﬁxed),
to gain insight in the behavior of the system. On Fig. 2 a horizontal
slice at high temperature is shown. The density of the fermions in
the system is measured, as deﬁned by
〈n〉 = 1
Ω
∂ ln Z
∂μ
, (14)
with Ω the space–time volume, in units of the saturation density
(which is reached when all available fermionic states on the lattice
are ﬁlled). The density starts to increase right away, there is no
sign of the Silver–Blaze phenomenon [27] at this high temperature,
as expected. Around μ/T = 8 the saturation is reached.
To measure the importance of the fermionic contribution to the
weight of the system, we deﬁne the average sign of the determi-
nant as
〈
e2iϕ
〉=
〈
detM(μ)
〉
. (15)detM(−μ)
110 D. Sexty / Physics Letters B 729 (2014) 108–111Fig. 2. Average phase factor and density as a function of the chemical potential.
Fig. 3. The fermion density, the chiral condensate (deﬁned by 〈∂ ln Z/∂m〉/Ω) and
the trace of the Polyakov loop and its inverse as a function of the chemical potential.
Since the calculation of the determinant is very costly, it is only
measured on small lattices, see Fig. 2. (For the Langevin dynamics
the calculation of the determinant is not needed.) One sees that
even on this small lattice the average sign is close to zero in a
big range of the physically interesting region, making reweighting
unfeasible. (The feasibility of the reweighting is controlled by the
sign average in the phasequenched system (where the determinant
in the measure is substituted with its absolute value), which be-
haves similarly to the sign average in the non-quenched system as
shown in Fig. 2.) In the saturation region the phase ﬂuctuations of
the determinant vanish again, as the necessary energy to create a
hole in the sea of fermions requires more energy than is available
in the thermal bath, thus the fermions decouple from the system.
In Fig. 3 we show the fermionic observables as well as the
trace of the Polyakov loops (deﬁned in (17)) and its inverse again
on a horizontal slice of the phase diagram. The expected physical
scenario is realized: the density of the fermions grows until sat-
uration, while the chiral condensate vanishes. The Polyakov loops
have a peak at some nonzero μ (with the inverse Polyakov loop
having a peak ﬁrst), before they decay to zero, as the Z3 sym-
metry of the system is restored in the saturation region, where
fermions no longer have an inﬂuence. Note that the critical β of
the system (the value for which the system is at the transition be-
tween conﬁned and deconﬁned phases) is around βc ≈ 5.5 for the
parameters used in Fig. 3, so the slice is slightly above the critical
temperature. To reach smaller temperatures, lattices using larger
temporal extent are needed.
A well known approximation to full QCD is heavy quark QCD
(HQCD), which is valid for heavy quarks and large chemical po-
tentials [28,29], see also [30,31]. In this approximation the spatialFig. 4. Comparison of the average densities measured in HQCD and in full QCD with
staggered fermions.
Fig. 5. Comparison of the average Polyakov loops measured in HQCD and in full QCD
with staggered fermions.
hoppings are dropped and the fermionic determinant simpliﬁes
considerably:
det
(
M(μ,U )
)=∏
x
det(1+ C Px)det
(
1+ C ′P−1x
)
(16)
with the Polyakov loop
Px =
NT −1∏
τ=0
U (τ ,x),4, (17)
and the parameters C = eμNT /(2m)NT and C ′ = e−μNT /(2m)NT
with the staggered mass m, and the temporal extent of the lat-
tice NT . Note that this is the ‘symmetrized’ form of the deter-
minant satisfying detM(−μ) = (detM(μ))∗ , otherwise the second
factor could be dropped in the heavy-dense limit. The correspond-
ing approximation for Wilson fermions was studied with the com-
plex Langevin method in an earlier publication [7]. The HQCD
approximation for one ﬂavor of Wilson fermion amounts to sub-
stituting m = 1/(4κ) in Eq. (16), as well as taking the square of
the right hand side of (16).
Increasing the quark mass, the HQCD approach will become a
better and better approximation of full QCD. To test at which mass
scale will the HQCD become quantitatively accurate, and to vali-
date the algorithm for full QCD, I compared simulations of HQCD
(for details, see [7,26]) to full QCD with staggered fermions, using
the same mass parameter. In Fig. 4 the fermion density is com-
pared, in Fig. 5, the Polyakov loops are compared.
One observes good agreement at the high mass of am = 4, as
expected, since the HQCD expansion is based on the expansion of
D. Sexty / Physics Letters B 729 (2014) 108–111 111the fermion determinant using the small parameter 1/am. This of
course does not prove that the results are fully reliable, but in-
creases the conﬁdence in the procedure, as the HQCD method was
validated with reweighting at small μ [7]. At smaller masses the
results are quantitatively different, but the qualitative behavior is
very similar, where the biggest effect on the density and on the
Polyakov loop seems to be a rescaling of the chemical potential.
In this Letter I have shown that ﬁnite density simulations of full
QCD using the CLE with gauge cooling all the way up to saturation
are feasible using small enough lattice spacings, where the cool-
ing is effective. This method avoids the sign and overlap problems,
direct simulation results in the high density region are presented
for the ﬁrst time. The cost of the simulation depends on the vol-
ume similarly to a hybrid Monte Carlo simulation, as the inversion
of the fermion matrix is the main numerical cost. In particular the
cost increases polynomially with the volume, in contrast with the
exponentially costly reweighting approach.
The results correctly reproduce the saturation physics and are
found to agree with HQCD for large quark masses. To increase
the conﬁdence in the reliability of the results, further checks are
needed in the regions where different approaches are available,
such as results at small chemical potentials [32], or the results
gained using strong coupling expansions.
Acknowledgements
I am indebted to Gert Aarts, Erhard Seiler and Ion-Olimpiu Sta-
matescu for many discussions and collaboration on related work.
A large part of the numerical calculations for this project was done
on the bwGRiD (http://www.bw-grid.de), member of the German
D-Grid initiative, funded by BMBF and MWFK Baden-Württemberg.
References
[1] P. de Forcrand, PoS LAT 2009 (2009) 010, arXiv:1005.0539;
G. Aarts, PoS LATTICE 2012 (2012) 017, arXiv:1302.3028.
[2] G. Parisi, Phys. Lett. B 131 (1983) 393.
[3] G. Aarts, I.-O. Stamatescu, J. High Energy Phys. 0809 (2008) 018, arXiv:
0807.1597.
[4] G. Aarts, Phys. Rev. Lett. 102 (2009) 131601, arXiv:0810.2089 [hep-lat].
[5] G. Aarts, K. Splittorff, J. High Energy Phys. 1008 (2010) 017, arXiv:1006.0332
[hep-lat].[6] G. Aarts, F.A. James, J. High Energy Phys. 1201 (2012) 118, arXiv:1112.4655.
[7] E. Seiler, D. Sexty, I.-O. Stamatescu, Phys. Lett. B 723 (2013) 213, arXiv:
1211.3709 [hep-lat].
[8] J. Ambjorn, S.K. Yang, Phys. Lett. B 165 (1985) 140.
[9] J. Ambjorn, M. Flensburg, C. Peterson, Nucl. Phys. B 275 (1986) 375.
[10] G. Aarts, F.A. James, J. High Energy Phys. 1008 (2010) 020, arXiv:1005.3468.
[11] J.M. Pawlowski, C. Zielinski, Phys. Rev. D 87 (2013) 094509, arXiv:1302.2249
[hep-lat].
[12] J.M. Pawlowski, C. Zielinski, Phys. Rev. D 87 (2013) 094503, arXiv:1302.1622
[hep-lat].
[13] J. Berges, I.-O. Stamatescu, Phys. Rev. Lett. 95 (2005) 202003, arXiv:hep-lat/
0508030.
[14] J. Berges, S. Borsanyi, D. Sexty, I.-O. Stamatescu, Phys. Rev. D 75 (2007) 045007,
arXiv:hep-lat/0609058.
[15] J. Berges, D. Sexty, Nucl. Phys. B 799 (2008) 306, arXiv:0708.0779.
[16] L. Bongiovanni, G. Aarts, E. Seiler, D. Sexty, I.-O. Stamatescu, arXiv:1311.1056
[hep-lat].
[17] G. Aarts, E. Seiler, I.-O. Stamatescu, Phys. Rev. D 81 (2010) 054508, arXiv:
0912.3360;
G. Aarts, F.A. James, E. Seiler, I.-O. Stamatescu, Eur. Phys. J. C 71 (2011) 1756,
arXiv:1101.3270.
[18] C. Pehlevan, G. Guralnik, Nucl. Phys. B 811 (2009) 519, arXiv:0710.3756;
G. Guralnik, C. Pehlevan, Nucl. Phys. B 822 (2009) 349, arXiv:0902.1503.
[19] G. Aarts, F.A. James, J.M. Pawlowski, E. Seiler, D. Sexty, I.-O. Stamatescu, J. High
Energy Phys. 1303 (2013) 073, arXiv:1212.5231 [hep-lat].
[20] A. Duncan, M. Niedermaier, Ann. Phys. 329 (2013) 93.
[21] G. Aarts, P. Giudice, E. Seiler, arXiv:1306.3075 [hep-lat].
[22] A. Mollgaard, K. Splittorff, arXiv:1309.4335 [hep-lat].
[23] G.G. Batrouni, G.R. Katz, A.S. Kronfeld, G.P. Lepage, B. Svetitsky, K.G. Wilson,
Phys. Rev. D 32 (1985) 2736.
[24] M. Fukugita, Y. Oyanagi, A. Ukawa, Phys. Rev. D 36 (1987) 824.
[25] D. Zwanziger, Nucl. Phys. B 192 (1981) 259;
P. Rossi, C.T.H. Davies, G.P. Lepage, Nucl. Phys. B 297 (1988) 287.
[26] G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty, I.-O. Stamatescu, arXiv:1303.6425
[hep-lat].
[27] T.D. Cohen, Phys. Rev. Lett. 91 (2003) 222001, arXiv:hep-ph/0307089.
[28] I. Bender, T. Hashimoto, F. Karsch, V. Linke, A. Nakamura, M. Plewnia, I.O. Sta-
matescu, W. Wetzel, Nucl. Phys. Proc. Suppl. 26 (1992) 323;
T.C. Blum, J.E. Hetrick, D. Toussaint, Phys. Rev. Lett. 76 (1996) 1019, arXiv:
hep-lat/9509002.
[29] R. De Pietri, A. Feo, E. Seiler, I.-O. Stamatescu, Phys. Rev. D 76 (2007) 114501,
arXiv:0705.3420.
[30] M. Fromm, J. Langelage, S. Lottini, O. Philipsen, J. High Energy Phys. 1201
(2012) 042, arXiv:1111.4953;
M. Fromm, J. Langelage, S. Lottini, M. Neuman, O. Philipsen, arXiv:1207.3005
[hep-lat].
[31] P. de Forcrand, M. Fromm, Phys. Rev. Lett. 104 (2010) 112005, arXiv:0907.1915
[hep-lat].
[32] S. Borsányi, Z. Fodor, S.D. Katz, D. Sexty, in preparation.


file:///data/remote/core/dit/data/elsevier/pdf/24a/aHR0cDovL2FwaS5lbHNldmllci5jb20vY29udGVudC9hcnRpY2xlL3BpaS9zMDM3MDI2OTMxNDAwMDMxOA==.pdf

Simulating full QCD at nonzero density using the complex Langevin equation

Abstract

Similar works

Full text

Available Versions

Elsevier - Publisher Connector

Crossref

Elsevier - Publisher Connector