Simulating full QCD at nonzero density using the complex Langevin equation
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.
The determination of the phase diagram of finite density QCD is one of the great problems of theoretical physics today. One is interested in averages defined with the Euclidean path integral
where is the Yang-Mills action of the gauge fields and is the Dirac-matrix of the quark fields. Naive lattice simulations at 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 circumvent the problem, but these are of limited use pdf (), mostly being applicable for . An exception is the complex Langevin method parisi (), 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 aartsstam (); Aarts:2008wh (); Aarts:2010gr (); aj (); gaugecooling (), but in some cases also non-physical results are delivered Ambjorn:1985iw (); Ambjorn:1986fz (); aarts (); Pawlowski:2013gag (); Pawlowski:2013pje (). In this paper I demonstrate that the algorithm can be extended to full QCD with light quark masses on lattices with sufficiently small lattice spacings.
The complex Langevin method is based on setting up a complex Langevin equation (CLE) in an enlarged manifold, which is the complexification of the original field space parisi (). The original theory is recovered by taking expectation values of the analytically continued observables. For gauge theories this complexification is . 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’, using the Minkowskian formulation of the path integral bergesstam (); bergessexty (); opt (), or Yang-Mills theory with -term theta (). In this work I am concerned with finite 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 ajss (); guralnik (); haarmeas (); Duncan:2013wm (); Aarts:2013uza (), one can gain an insight whether the results are trustworthy using requirements such as the fast decay of the distributions.
Recently an important breakthrough in this field was the development of a ’gauge cooling’ algorithm for the CLE method gaugecooling (), where the gauge symmetry of the system is used to ensure a well localized distribution in the complexified field space, and thus convergence to the correct results.
In this work the CLE method is applied to the lattice discretization of full QCD, i.e. for the action
where is the Wilson plaquette action for the SU(3) link variables, and is the unimproved staggered fermion determinant for fermion flavors
where and indices represent spacetime coordinates, and are the staggered sign functions. Periodic (antiperiodic) boundary conditions are used in space (time) directions. The fermion matrix fulfills the symmetry condition:
with the “staggered matrix”, . This symmetry leads to . This means that the determinant becomes complex for , making a simulation based on importance sampling impossible. Without rooting (i.e. using 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 flavors is trivial, appears as a factor of a drift term. In this study I have chosen to use and .
The Langevin equation for the link variables is set up using the equation
Here are the generators of the gauge group, i.e. the Gell-Mann matrices. The drift force is determined by
with the left derivative
The drift term for the action (2) is written as
It has been suggested in logcut () that the drift corresponding to the action (2) might also include a term reflecting the branch cut of the complex logarithm on the negative real axis. This yet unclarified issue is the subject of ongoing research.
The drift term remains real only for . 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 batrouni (); fukugita1987 (), which is related to pseudofermionic variables, is introduced as follows. The drift of the link variables is calculated using
where the is a vector of Gaussian random numbers satisfying . For the calculation of the drift term one has to solve the linear system of equation . In terms of the solution the drift term is written as
This means that this algorithm uses the conjugate gradient (CG) algorithm once for every timestep for the solution of the linear system. In fukugita1987 () it was also shown that with a higher order algorithm one can get rid of part of the 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 complexified field space is desirable. A useful measure of the size of the distribution in imaginary directions of a link variable is the unitarity norm
where the equality is reached only for matrices. The enlarged gauge symmetry of the system can be used to decrease the unitary norm of the system, thus ensure convergence to the exact results opt (). Recently, we developed and tested in HQCD (see below) a procedure utilizing this freedom called gauge coolinggaugecooling () (reminiscent of stochastic gauge fixingstocgaugefix ()). The idea is the following: one uses gauge transformations
with to decrease the unitarity norm of the system. This can be accomplished by choosing the matrices in the direction 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 Aarts:2013uxa ().
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 simulation. Generally, we find also at , using sufficient cooling, that the level of unitarity norm stabilizes and allows one to obtain correct results (after extrapolation) for lattices with fine enough lattice spacings. As one observes, at low (below for ) the cooling is not effective enough to prevent the system from wandering off far from the manifold, and ’skirted’ distributions develop (as also observed in gaugecooling ()). Close to the continuum limit, however, the algorithm seems to be stable irrespective of the physical phase, as observed using cheaper HQCD simulations.
Observables are measured on ’slices’ of the phase diagram (meaning a scan using one variable while keeping the other fixed), 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 defined by
with the space-time volume, in units of the saturation density (which is reached when all available fermionic states on the lattice are filled). The density starts to increase right away, there is no sign of the Silver-Blaze phenomenon silverblaze () at this high temperature, as expected. Around the saturation is reached.
To measure the importance of the fermionic contribution to the weight of the system, we define the average sign of the determinant as
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 behaves similarly to the sign average in the non-quenched system as shown in Fig. 2.) In the saturation region the phase fluctuations 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 (defined 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 saturation, while the chiral condensate vanishes. The Polyakov loops have a peak at some nonzero (with the inverse Polyakov loop having a peak first), before they decay to zero, as the symmetry of the system is restored in the saturation region, where fermions no longer have an influence. Note that the critical of the system (the value for which the system is at the transition between confined and deconfined phases) is around 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 potentials oldhqcd (); feo (), see also owe (); deForcrand:2009dh (). In this approximation the spatial hoppings are dropped and the fermionic determinant simplifies considerably:
with the Polyakov loop
and the parameters and with the staggered mass , and the temporal extent of the lattice . Note that this is the ’symmetrized’ form of the determinant satisfying , otherwise the second factor could be dropped in the heavy-dense limit. The corresponding approximation for Wilson fermions was studied with the complex Langevin method in an earlier publication gaugecooling (). The HQCD approximation for one flavor of Wilson fermion amounts to substituting 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 validate the algorithm for full QCD, I compared simulations of HQCD (for details, see gaugecooling (); Aarts:2013uxa ()) to full QCD with staggered fermions, using the same mass parameter. In Fig. 4 the fermion density is compared, in Fig. 5, the Polyakov loops are compared.
One observes good agreement at the high mass of , as expected, since the HQCD expansion is based on the expansion of the fermion determinant using the small parameter . This of course does not prove that the results are fully reliable, but increases the confidence in the procedure, as the HQCD method was validated with reweighting at small gaugecooling (). 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 paper I have shown that finite 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 cooling is effective. This method avoids the sign and overlap problems, direct simulation results in the high density region are presented for the first time. The cost of the simulation depends on the volume 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 confidence 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 compare (), or the results gained using strong coupling expansions.
Acknowledgments. – I am indebted to Gert Aarts, Erhard Seiler and Ion-Olimpiu Stamatescu 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.
- (1) P. de Forcrand, PoS LAT 2009, 010 (2009) [arXiv:1005.0539]; G. Aarts, PoS LATTICE 2012 (2012) 017 [arXiv:1302.3028].
- (2) G. Parisi, Phys. Lett. 131 B (1983) 393.
- (3) G. Aarts and I. -O. Stamatescu, JHEP 0809 (2008) 018 [arXiv:0807.1597].
- (4) G. Aarts, Phys. Rev. Lett. 102 (2009) 131601 [arXiv:0810.2089 [hep-lat]].
- (5) G. Aarts and K. Splittorff, JHEP 1008, 017 (2010) [arXiv:1006.0332 [hep-lat]].
- (6) G. Aarts and F. A. James, JHEP 1201 (2012) 118 [arXiv:1112.4655].
- (7) E. Seiler, D. Sexty and I. -O. Stamatescu, Phys. Lett. B 723, 213 (2013) [arXiv:1211.3709 [hep-lat]].
- (8) J. Ambjorn and S. K. Yang, Phys. Lett. B 165, 140 (1985).
- (9) J. Ambjorn, M. Flensburg and C. Peterson, Nucl. Phys. B 275, 375 (1986).
- (10) G. Aarts and F. A. James, JHEP 1008 (2010) 020 [arXiv:1005.3468].
- (11) J. M. Pawlowski and C. Zielinski, Phys. Rev. D 87, 094509 (2013) [arXiv:1302.2249 [hep-lat]].
- (12) J. M. Pawlowski and C. Zielinski, Phys. Rev. D 87, 094503 (2013) [arXiv:1302.1622 [hep-lat]].
- (13) J. Berges and I. -O. Stamatescu, Phys. Rev. Lett. 95 (2005) 202003 [hep-lat/0508030].
- (14) J. Berges, S. .Borsanyi, D. Sexty and I. -O. Stamatescu, Phys. Rev. D 75 (2007) 045007 [hep-lat/0609058].
- (15) J. Berges and D. Sexty, Nucl. Phys. B 799 (2008) 306 [arXiv:0708.0779].
- (16) L. Bongiovanni, G. Aarts, E. Seiler, D. Sexty and I. -O. Stamatescu, arXiv:1311.1056 [hep-lat].
- (17) G. Aarts, E. Seiler and I. -O. Stamatescu, Phys. Rev. D 81 (2010) 054508 [arXiv:0912.3360]; G. Aarts, F. A. James, E. Seiler and I. -O. Stamatescu, Eur. Phys. J. C 71 (2011) 1756 [arXiv:1101.3270].
- (18) C. Pehlevan and G. Guralnik, Nucl. Phys. B 811 (2009) 519 [arXiv:0710.3756]; G. Guralnik and C. Pehlevan, Nucl. Phys. B 822 (2009) 349 [arXiv:0902.1503].
- (19) G. Aarts, F. A. James, J. M. Pawlowski, E. Seiler, D. Sexty and I. -O. Stamatescu, JHEP 1303, 073 (2013) [arXiv:1212.5231 [hep-lat]].
- (20) A. Duncan and M. Niedermaier, Annals Phys. 329, 93 (2013).
- (21) G. Aarts, P. Giudice and E. Seiler, arXiv:1306.3075 [hep-lat].
- (22) A. Mollgaard and K. Splittorff, arXiv:1309.4335 [hep-lat].
- (23) G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, B. Svetitsky and K. G. Wilson, Phys. Rev. D 32 (1985) 2736.
- (24) M. Fukugita, Y. Oyanagi and A. Ukawa, Phys. Rev. D 36 (1987) 824.
- (25) D. Zwanziger, Nucl. Phys. B 192 (1981) 259; P. Rossi, C. T. H. Davies and G. P. Lepage, Nucl. Phys. B 297 (1988) 287.
- (26) G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty and I. -O. Stamatescu, arXiv:1303.6425 [hep-lat].
- (27) T. D .Cohen, Phys. Rev. Lett. 91 (2003) 222001 [hep-ph/0307089].
- (28) I. Bender, T. Hashimoto, F. Karsch, V. Linke, A. Nakamura, M. Plewnia, I. O. Stamatescu and W. Wetzel, Nucl. Phys. Proc. Suppl. 26 (1992) 323; T. C. Blum, J. E. Hetrick and D. Toussaint, Phys. Rev. Lett. 76 (1996) 1019 [hep-lat/9509002]
- (29) R. De Pietri, A. Feo, E. Seiler and I. -O. Stamatescu, Phys. Rev. D 76 (2007) 114501 [arXiv:0705.3420].
- (30) M. Fromm, J. Langelage, S. Lottini and O. Philipsen, JHEP 1201 (2012) 042 [arXiv:1111.4953]; M. Fromm, J. Langelage, S. Lottini, M. Neuman and O. Philipsen, arXiv:1207.3005 [hep-lat].
- (31) P. de Forcrand and 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.