Gauge cooling in complex Langevin for QCD with heavy quarks

We employ a new method,"gauge cooling", to stabilize complex Langevin simulations of QCD with heavy quarks. The results are checked against results obtained with reweigthing; we find agreement within the estimated errors. The method allows us to go to previously unaccessible high densities.

(Dated: May 5, 2014) We employ a new method, "gauge cooling", to stabilize complex Langevin simulations of QCD with heavy quarks. The results are checked against results obtained with reweigthing; we find agreement within the estimated errors. The method allows us to go to previously unaccessible high densities.
PACS numbers: 11.15.Ha,11.15Bt Introduction. -An important problem of particle physics is to derive the phase diagram of hot and dense QCD from first principles. The difficulty lies in the fact that the action becomes complex at finite density. Various methods have been employed to get around this problem: Taylor expansions at µ = 0 and analytic continuation from imaginary µ, reweighting etc., but these have a limited range of applicability [1]. A general method that recently has received a lot of attention is the Complex Langevin Equation (CLE) [2][3][4][5][6][7][8] . It has had some successes as well as some failures; the failures seem to be related to the problem that the system wants to spend "too much time too far out" in the complexified configuration space; this problem and possible ways to deal with it were discussed in [8,9].
In QCD at finite density the enlarged gauge freedom in the complexified field space can actually open a new way to limit the dangerous large excursions by employing nonunitary gauge transformations. This freedom was used in a somehwat different way for a SU (2) toy model [9]; the procedure used there, however, does not generalize to a lattice model in an obvious way.
Here we introduce the method of gauge cooling (g.c.), designed to keep the link variables as close as possible to the unitaries, thereby preventing the dangerous large excursions. We test the idea first for a simple Polyakov loop model, where we have exact results to compare our data with; then we apply it to a heavy quark approximation of QCD with chemical potential (HQCD) (cf. [7]).
Complex Langevin for gauge models -The complex Langevin method [10,11] for a complex action S is based on setting up a stochastic process on the complexification of the configuration space. The longtime average of holomorphic observables is then supposed to give the correct average corresponding to the complex weight exp(−S). In lattice models of QCD the configuration space is SU (3) #links (see [12] for the Real Langevin approach); after complexification this becomes SL(3, C) #links [13] λ a , a = 1 . . . , 8 are the Gell-Mann matrices; dw x,µ,a are independent Wiener increments, normalized as D a,x,µ is a differential operator (derivation) acting on where {U (δ)} means the variable U x,µ has been replaced by exp(iδλ a )U x,µ with all other variables unchanged. The first term on the rhs of Eq.(1) (the drift term) is gauge covariant and transverse to the gauge orbits since S is gauge invariant, but the noise term contains components along the gauge orbits. Since the gauge group is noncompact, the process (1) may go far from the unitary submanifold.
An Euler discretization of Eq.(1) (see also [6]) gives where K a,x,µ = D a,x,µ S is the drift force and η are independent Gaussian noises satisfying Gauge cooling. -It is a well known problem of the CLE method that the system may drift too far out into the imaginary directions, causing numerical problems. To quantify the distance from unitarity we use the 'unitarity norm' Gauge invariance is the freedom to transform variables for all links x, x +μ, where initially V x , V x+μ are unitary, but after complexification V x , V x+μ ∈ SL(3, C). We use this freedom to minimize the unitarity norm, without changing the observables. Concretely, we intersperse standard 'dynamical' Langevin sweeps with several 'g.c. sweeps'. These are deterministic moves essentially in the direction opposite to the gradient of F .
A maximally nonunitary gauge transformation at site y in direction a is given by (α ∈ R) for all µ, with all other link variables remaining unchanged. The 'gauge gradient' of F in the a direction at the lattice site y is then The g.c. updates of the configuration are given by whereα = ǫα; α determines the strength of the g.c. force, whereas ǫ is a discretization parameter as in Eq.(4). Note that even ifα is not small, Eq.(10) is still a gauge transformation; it just might not be optimal for reducing F .
Polyakov loop model. -The Polyakov loop model is given by a 1D lattice consisting of N links with periodic boundary conditions. Analytically it reduces to a onelink integral, but it is a useful laboratory to check the effect of g.c.. The action is given by where we allow β 1,2 to be complex. Here we choose β 1 = β + κe µ , β 2 = β * + κe −µ ; S will in general be complex. The observables of interest are tr P k = tr(U 1 . . . U N ) k , k = ±1, ±2, ±3. The effect of the g.c. on tr P is shown in Fig.1, where we vary α of Eq.(10). In Fig. 2 we show the effect of varying α on the distribution of the values of tr P. We see that with too small α the distribution has a 'skirt' of slow decay. From [8] we know that slow decay typically leads to incorrect results; a small skirt, however, does not lead to appreciable deviations from the exact values. The main point is that enough cooling leads to correct results, see Fig.1.
We also measured tr U ±2 and tr U ±3 with similar results, for N up to 1024; it turns out that increasing N requires an increase in cooling.  Heavy quark QCD (HQCD). -The model was first explored using the CLE without g.c. in [7] on small lattices and for moderate µ. It is defined by dropping all spatial hopping terms of the quarks; this amounts to simplifying the Wilson fermion determinant of lattice QCD with nonzero chemical potential µ to HQCD is described by the action where S G = P Re tr(U P ) is the Wilson plaquette action and antiperiodic b.c. in time are chosen for the fermions. We have  and similarly for the second factor in Eq.(12), where A related model has been studied numerically by a reweighting (RW) technique to deal with the complex density in [14] (and more recently in [15]), where for larger µ and larger lattices the sign problem causes bad signal to noise ratios. See [14] for the general setting. We simulate the model by the CLE on various lattices, with κ = 0.12, various values of β and µ, using Eq.(4) with ǫ = 10 −5 − 2 × 10 −4 . In some cases we use adaptive control for the dynamical stepsize and the number of cooling steps. The system is thermalized after a cold start up to Langevin time t = 10 before taking averages.
In Fig.3 we show the Langevin evolution of the unitarity norm: it is seen how g.c. stabilizes the process. We see that g.c. is needed even for µ = 0 where the process is real, but has an instability pushing it away from the unitary submanifold in the absence of cooling.
In Fig.4 we show the baryon density n/n sat as a function of the chemical potential µ on an 8 3 × 6 lattice. We find saturation for µ 2, suggesting that our results are correct also for fairly large µ. The figure also shows the average phase of exp(−S), defined by (see [7]). With growing density the average phase drops to almost zero and rises again in the saturation region. Fig.5 shows P and P ′ vs. µ. Both reach a maximum where the density takes off (µ ≈ 1.4), then drop again to nearly zero (see also [16][17][18]); the curve for P ′ is shifted slightly to the left with respect to P . The analytic curves for β = 0 in the figure exhibit a similar behavior; they are calculated essentially as in [14]. This behavior is due to the fact that for small as well as large µ the determinants in Eqs. (12) are dominated by the constant terms, leading to small expectation values of P and P ′ . In physics terms, to get a nonzero answer one has to form a localized colorless bound state between the dynamical quarks and the external charge provided by P or P ′ . For small µ > 0, when typically there is only 1 quark present, it can form a 'meson' with P ′ but not with P ; for larger µ, when typically 2 quarks are present, P , but not P ′ can form a 'baryon' with them. The situation is slightly complicated by the fact that up to 6 quarks can exist at every site, so P ′ can also form a colorless state with 4 and P with 5 quarks. With no quarks or in the completely filled state it is not possible to form such bound states. The simple picture given here ignores fluctuations, including a small probability, damped by µ, of finding antiquarks, so at µ = 0 P and P ′ are small but nonzero, whereas at very large µ they will go to zero. Thus both quantities will have a peak, but P will peak later than P ′ .
Comparison of CLE with RW. -To check the CLE results we also simulated HQCD with a RW technique, in the region where this is feasible (cf. [14]). We compare the results for the Polyakov loops and for plaquette   Fig.6. The results agree within the errors (large errors only occur for RW). Note that around µ = 1.2 the RW method breaks down, whereas CLE with g.c. works fine for arbitrarily large µ.

expectations in
In Fig.7 we show the expectation values of Polyakov loops at µ = 0.85 vs. β, both obtained from RW and CLE. The figure shows that CLE works for large β and also allows us to cross over into the confining region. Deeper in the confining region, however, we see small differences between the RW and CLE data, apparently related to an emergence of a 'skirt' of the distribution.
Concluding remarks. -The instabilities that arise without or with too weak g.c. are apparently caused by (1) the existence of repulsive fixed points, pushing the configurations exponentially fast away from the unitary submanifold and (2) roundoff and other numerical errors. Each gauge orbit contains configurations arbitrarily far from the unitaries, and there roundoff errors pile up to affect also the observables, violating gauge invariance.
The modified process remedies these problems; it thereby also avoids the slow decay of the equilibrium distribution. The CLE method with g.c. is tested successfully for the Polyakov loop model. It allows, apparently for the first time, to simulate a real QCD like gauge model at finite density all the way into the saturation region. The method does not suffer from any sign or overlap problem. There is no fundamental obstacle to extend the method to full QCD.