Applying Complex Langevin Simulations to Lattice QCD at Finite Density

We study the use of the complex-Langevin equation (CLE) to simulate lattice QCD at a finite chemical potential ($\mu$) for quark-number, which has a complex fermion determinant that prevents the use of standard simulation methods based on importance sampling. Recent enhancements to the CLE specific to lattice QCD inhibit runaway solutions which had foiled earlier attempts to use it for such simulations. However, it is not guaranteed to produce correct results. Our goal is to determine under what conditions the CLE yields correct values for the observables of interest. Zero temperature simulations indicate that for moderate couplings, good agreement with expected results are obtained for small $\mu$ and for $\mu$ large enough to reach saturation, and that this agreement improves as we go to weaker coupling. For intermediate $\mu$ values these simulations do not produce the correct physics. We compare our results with those of the phase-quenched approximation. Since there are indications that correct results might be obtained if the CLE trajectories remain close to the $SU(3)$ manifold, we study how the distance from this manifold depends on the quark mass and on the coupling. We find that this distance decreases with decreasing quark mass and as the coupling decreases, i.e. as the simulations approach the continuum limit.

the continuum limit.

I. INTRODUCTION
QCD at finite baryon/quark-number density describes nuclear matter, the constituent of neutron stars and the interiors of heavy nuclei. We are interested in calculating the phase diagram for nuclear matter. Knowing the properties of nuclear matter can yield an equationof-state and a better description of neutron stars. In addition it can yield information on the interaction of nuclear matter with particles passing through it. Nuclear matter at high temperatures undergoes a transition to a quark-gluon plasma. Such transitions can be observed in relativistic heavy-ion collisions. While for low densities this transition is expected to be a crossover, it is believed that for high densities it becomes a first-order transition. The critical end-point is the second-order transition (critical point) where this change occurs.
Finite density QCD has a sign problem which prevents the direct application of standard lattice QCD simulation methods, which rely on importance sampling. When finite density is implemented by use of a quark-number chemical potential, µ, the sign problem manifests itself by making the fermion determinant complex, with a real part of indefinite sign. Simulations using the complex Langevin equation (CLE) [1][2][3][4] can accommodate such complex actions. However, the CLE can only be shown to yield correct values for observables if the space over which the fields evolve is compact, the drift (force) term is holomorphic in the fields, and the solutions are ergodic [5][6][7][8][9][10][11][12][13].
For QCD, implementation of the CLE requires extending the gauge-field manifold from SU(3) to SL (3, C). Keeping the action holomorphic in the fields except on a space of measure zero requires making the action a function of the gauge fields and their inverses only. Then the drift term is meromorphic in the gauge fields, having poles at the zeros of the fermion determinant. Early attempts at simulations were frustrated by runaway behaviour which could not be controlled by adaptively decreasing the updating interval. Recently it was discovered that, for small enough couplings, this behaviour could be tamed by 'gaugecooling', gauge transforming the fields after each update to keep them as near as possible to the SU(3) manifold [14]. When this is done, the gauge fields appear to evolve over a compact manifold, at least for weak enough coupling. It remains to be determined under what conditions the CLE produces correct values for observables, despite the presence of poles in the drift term and/or in the operators defining these observables.
Preliminary work directed towards zero temperatures has been reported [25]. In heavy dense lattice QCD at finite chemical potential, the zeros of the fermion determinant (poles in the drift term) only affect the CLE results very close to the transition, provided one excludes the regions where the real part of the fermion determinant is negative, when necessary.
Alternatively, good results for heavy quarks can be obtained by restricting the length of CLE trajectories to keep them close to the SU(3) manifold. For lighter quarks, CLE simulations are found to produce good results at high temperatures (above the finite temperature phase transition), but the situation at lower temperatures is less clear.
We simulate zero temperature lattice QCD with 2 flavours (tastes) of staggered quarks at µ values from zero to saturation, using the CLE. Here we are interested in the phase transition from hadronic to nuclear matter which should occur at µ ∼ m N /3. Random matrix theories related to QCD at finite µ suggest that when CLE simulations fail, they produce the results of the phase-quenched model (the theory where the fermion determinant is replaced by its magnitude), which has a transition to a superfluid state with a pion-like condensate at µ ≈ m π /2 [26]. (Note that other random matrix CLE simulations seem more optimistic [27].) For this reason we also perform RHMC simulations of the phase-quenched theory at the same values of β = 6/g 2 and quark mass m over the same range of µ values and on the same lattice size as the full theory, for comparison. Hence it is important that we choose m such that m N /3 is significantly larger than m π /2. We perform our simulations at β = 5.6, m = 0.025 (in lattice units) on a 12 4 lattice and at β = 5.7, m = 0.025 on a 16 4 lattice. Preliminary results were reported at Lattice 2015-2018. See [28] and its references to our earlier talks.
What we find is that there is reasonable agreement with the expected values for observables at β = 5.6, for µ << m π /2 and µ >> m N /3, and significantly better agreement for β = 5.7. However, for µ in the transition region, these simulations fail to produce the expected physics. The transition to nuclear matter appears to start at µ even less than m π /2 instead of ∼ m N /3. There is some indication at β = 5.7 that the CLE might yield correct results for µ > ∼ m N /3. Since there are indications that the CLE can produce correct results when its trajectories remain close to the SU(3) manifold, we perform a systematic study of how the average distance to this manifold (measured using the 'unitarity norm') depends on the quark mass and the coupling. We find that this norm decreases with decreasing quark mass and with decreasing coupling, i.e. as we approach the continuum limit.
In section 2 we present the formulation of the CLE, we use. Section 3 describes our simulations at β = 5.6 and β = 5.7 and presents results. In section 4 we present our simulations to determine how the unitarity norm depends on quark mass m and lattice coupling g. Section 5 gives a summary, discussion and conclusions.

II. COMPLEX LANGEVIN FOR FINITE DENSITY LATTICE QCD
If S(U) is the gauge action after integrating out the quark fields, the Langevin equation for the evolution of the gauge fields U in Langevin time t is: where l labels the links of the lattice, and η l = a η a l λ a . Here λ a are the Gell-Mann matrices for SU(3). η a l (t) are Gaussian-distributed random numbers normalized so that: We note in passing that the discretized Langevin Equation is the limiting case of the Hybrid Molecular Dynamics method where each trajectory has only a single update.
The complex-Langevin equation has the same form except that the Us are now in where M(U, µ) is the staggered Dirac operator. Note: backward links are represented by Note also that we have chosen to keep the noise term η real.
To simulate the time evolution of the gauge fields we use the partial second-order formalism of Fukugita, Oyanagi and Ukawa. [29][30][31] For an update of the fields by a 'time' increment dt, this gives: where γ = 1 2 + 1 4 dt and the Gaussian noise η is normalized such that: To proceed, we replace the spacetime trace with a stochastic estimator ξ where ξ is a vector over space-time and colour of Gaussian random numbers, normalized such that: which means, in particular, that the ξs in X 0 and X 1 are independent, unlike the ηs. After performing δ δU of ln(M) it is useful to use the cyclic property of the trace to rearrange the terms proportional to U and U −1 prior to introducing the stochastic estimators, so that this operator is antihermitian when µ = 0 and U is unitary. That way, in this special case, the complex Langevin equation becomes the real Langevin equation.
We apply adaptive updating, where if the force term becomes too large, dt is decreased to keep it under control. Because the Dirac operator is often ill-conditioned, we use 64-bit floating point precision throughout.
After each update, we adaptively gauge fix to the gauge which minimizes the unitarity norm:  theory. At µ = 0 the CLE value of this condensate at β = 5.6 is ≈ 6.6% too low, while at β = 5.7 it is ≈ 1.2% too low, a considerable improvement. As µ is increased from zero, the condensates for both βs start to fall for µ < m π /2 rather than for µ ≈ m N /3. Therefore, even though β = 5.7 shows an improvement over β = 5.6, it is still a worse approximation to the physics than is the phase-quenched theory. It is still unclear if going to a much weaker   theory look almost identical, except on their approach to saturation. Figure 4, shows an expanded version of the low µ region which indicates that the apparent agreement is only because the number densities are so close to zero for small µ. The results for the 2 theories are closer for β = 5.7 than for β = 5.6. This would be expected if the CLE produces phasequenched results in the weak-coupling limit, but also could indicate that it produces correct results in that limit. Note that it is difficult, if not impossible, to determine the positions of the transitions from these quark-number density graphs. Figure 5 shows the average unitarity norms for the CLE simulations at β = 5.6 and β = 5.7 as functions of µ. In both cases this norm decreases from its value at µ = 0 as µ is increased reaching a minimum around 0.5 or 0.6 before increasing to a maximum at saturation, which is the value for CLE simulations of the pure gauge theory at the same β value. The values of the unitarity norms for β = 5.7 lie below the values of those at β = 5.6 for the same µ. The minimum for β = 5.7 is ≈ 0.0054, far below the 0.1 which some have suggested is the approximate value below which the CLE can yield correct results. We thus suggest that, at β = 5.7, the CLE will produce correct results for µs around this rather broad minimum. Since at β = 5.7, the CLE appears to produce reliable results at saturation, it is reasonable to assume that observables measured by the CLE will have approximately correct values from near this minimum up to saturation. One might even hope that this range could be extended down to near µ = m N /3.
To test whether the assumption that the CLE is valid for β = 5.7, m = 0.025 for µ near

IV. DEPENDENCE OF THE UNITARITY NORM ON QUARK MASS AND COUPLING
It has been observed that the CLE is better behaved if its trajectories remain near to the SU(3) manifold, i.e. if the unitarity norm remains small. At least for large quark masses, small means less than ∼ 0.1. It is therefore of interest to determine how the unitarity norm depends on the parameters of the theory.
First we consider how the unitarity norm depends on the quark mass m. In the previous section, we have observed that, for fixed m and β, the local maximum for small µ occurs at µ = 0. The global maximum occurs at saturation. Since this maximum is the value for the pure gauge theory at this β, it does not depend on m. Therefore we shall determine the m dependence of the unitarity norm at µ = 0, which is thus relevant for small µ.
Since the unitarity norm at µ = 0 appears to decrease with increasing β and we shall be interested in β ≥ 5.6, we study the m dependence with β fixed at β = 5.6. We perform updates.
In figure 8 we plot the unitarity norm versus 1/m at β = 5.6. We note that it decreases by almost an order of magnitude as m is decreased from infinity to 0.01.
For a given β, the unitarity norm has its maximum for µ at saturation, for which it is mass independent. Hence we choose to perform a CLE at saturation for each β, that is we simulate the pure gauge theory for that β, knowing that this will give the upper bound to the unitarity norm for that β. Moreover, since there are no quarks, we do not have to worry that we should really change m when we change β so as to keep on a line of constant physics. Since there are no quarks, these simulations are fast, and for given β one can use a much smaller lattice than for the same β with light quarks. Without quarks, the drift term is holomorphic in the gauge fields, so the CLE should produce results correct to order dt 2 , provided the fields evolve on a compact domain.
We perform CLE simulations for a selection of βs in the range 5.  figure 9. These decrease as a β increases, that is as the coupling g decreases. In fact over the range of βs considered, this norm decreases by more than an order of magnitude.
In figure 10 we plot the average plaquette from our CLE simulations and compare it to the 'exact' value from a Monte-Carlo simulation. Except for β = 5.6 these values are close and get closer as β is increased. This difference for the lower βs is too large to be due to statistical errors or the inexact nature of Langevin simulations. (This is surprising since the drift term is holomorphic in the fields, and the region spanned by these CLE For µ > m N /3 we have not seen any sign of new exotic phases, such as a coloursuperconducting phase. Nor is there any indication of a transition to quark-matter below saturation. There is also no indication of a difference between the full theory and its phasequenched approximation. Such a difference might be expected because of the existence of a pion-like superfluid phase in the latter. Because there is some indication that the CLE is more likely to produce correct results if the trajectories stay close to the SU(3) manifold, we have performed CLE simulations over a range of quark masses at a fixed coupling, and over a range of couplings at infinite quark mass (pure SU(3) gauge theory). What we find is that the average distance of the trajectories from SU(3) as determined by the unitarity norm decreases as the coupling and quark mass are decreased, that is as we approach the continuum limit. This gives us some hope that the CLE might produce the correct physics in this limit. However, random matrix theories suggest that it might produce phase-quenched results. Since the simulations described in section 3 do not rule out either possibility, further simulations at weaker couplings, which require larger lattices, are needed.
Other methods have been suggested to try and obtain correct results from CLE simulations. One is to add terms to the action which cause the CLE to avoid the poles, with coefficients, which when taken to zero, yield the original action [36]. The question then is can these coefficients be taken small enough to allow them to be continued to zero, without them losing their effectiveness in avoiding the poles. Preliminary results on very small lattice look promising, but it remains to be seen if this method will work on larger lattices.
A second method is based on the observation that one needs to keep the unitarity norm small for the CLE to work. This is achieved by changing the dynamics of the CLE by adding a force in the direction of decreasing unitarity norm to the drift term, with a coefficient which can be made arbitrarily small [37]. This additional force should be irrelevant (in the re-normalization group sense) so that it will vanish as the lattice spacing goes to zero.
It should also vanish when the gauge fields lie on the SU(3) manifold. Such a force will not be holomorphic in the gauge fields, so that adding it to the drift term could completely destroy any relationship between the CLE and the physics contained in the original functional integral, so careful testing is needed.
The use of other actions should be investigated. For example, use of Wilson fermions should be studied, since they will preserve the continuum order of the zeros of the fermion determinant, which is equal to the number of flavours. For staggered fermions this order is decreased by a factor of 4 by 'taste' breaking, and is only recovered in the continuum limit.
Application of the CLE to finite temperature QCD at finite µ, in particular with regard to the effect of finite µ on the transition of hadronic/nuclear matter to a quark-gluon plasma, should be studied. It has been observed that the CLE should work better in this case, since finite temperatures move the zeros of the fermion determinant away from the CLE trajectories. In the case of 2-flavour staggered quarks, our studies indicate that this is only likely to work (at least for m = 0.025) if β = 5.6 lies in the low temperature phase. This requires N t > ∼ 12.