Complex Langevin Simulations of QCD at Finite Density -- Progress Report

We simulate lattice QCD at finite quark-number chemical potential to study nuclear matter, using the complex Langevin equation (CLE). The CLE is used because the fermion determinant is complex so that standard methods relying on importance sampling fail. Adaptive methods and gauge-cooling are used to prevent runaway solutions. Even then, the CLE is not guaranteed to give correct results. We are therefore performing extensive testing to determine under what, if any, conditions we can achieve reliable results. Our earlier simulations at $\beta=6/g^2=5.6$, $m=0.025$ on a $12^4$ lattice reproduced the expected phase structure but failed in the details. Our current simulations at $\beta=5.7$ on a $16^4$ lattice fail in similar ways while showing some improvement. We are therefore moving to even weaker couplings to see if the CLE might produce the correct results in the continuum (weak-coupling) limit, or, if it still fails, whether it might reproduce the results of the phase-quenched theory. We also discuss action (and other dynamics) modifications which might improve the performance of the CLE.


Introduction
We study the physics of nuclear matter, the constituent of neutron stars and the interiors of heavy nuclei. In addition, hot nuclear matter is produced in relativistic heavy-ion colliders, and was present in the early universe. Ignoring the important effects of the electromagnetic field, nuclear matter is QCD at finite quark-/baryon-number density.
Lattice QCD at a finite quark-number chemical potential µ has a complex fermion determinant. Hence standard simulation methods based on importance sampling fail. The Langevin approach does not rely on importance sampling and can be extended to complex actions [1][2][3][4]. For lattice QCD, this requires analytic continuation of the gauge fields from S U (3) to S L (3, C). However, complex Langevin (CLE) simulations cannot be guaranteed to work unless the trajectories are restricted to a finite domain and the drift (force) term is holomorphic in the fields. Detailed discussion of these requirements as they pertain to QCD at finite µ are described in [5][6][7][8][9][10][11]. In lattice QCD at weak coupling, potential runaway solutions are controlled by adaptively readjusting the updating 'time' increment and implementing gauge cooling -choosing a gauge which minimizes the average unitarity norm, which is a measure of the distance of the gauge fields from the S U(3) manifold [12].
For lattice QCD at finite µ, the gauge-cooled CLE trajectories do appear to be confined to compact domains. However, zeros of the fermion determinant produce poles in the drift term, so that it is meromorphic not holomorphic in the fields. Hence observables cannot be guaranteed to converge to the correct limits. Simulations of heavy-quark lattice QCD where the hopping parameter expansion is used, have been performed using the CLE [13][14][15][16][17]. More limited simulations have been performed for lighter quark masses [15,[18][19][20]. We have been testing the CLE for lattice QCD at finite µ at zero temperature, over a range of µ values from 0 to saturation.
Our earlier simulations at β = 6/g 2 = 5.6, m = 0.025 mainly on 12 4 lattices show some of the properties expected of finite density QCD, but fail in the details [21,22]. We are now extending these to weaker coupling, β = 5.7, on a 16 4 lattice, with m = 0.025. While this improves the agreement with known (RHMC) observables at µ = 0, the behaviour for µ > 0 through the transition region shows similar deficiencies to those at β = 5.6 albeit with some improvement. For µ values above the transition region observables from our β = 5.7 CLE simulations agree with expectations whereas β = 5.6 exhibits small discrepancies.
One important improvement over β = 5.6 is that the unitarity norm is significantly smaller for β = 5.7 than for β = 5.6. Since the experience of others has identified keeping this norm low as a a way of improving CLE results, this encourages one to try even weaker couplings requiring even larger lattices.
Recent results from random matrix theory indicate that, when the complex Langevin fails, it produces results consistent with the phase-quenched theory [23]. Preliminary indications from our simulations are that this is not true for lattice QCD at finite µ for these couplings. It is still an open question as to whether the CLE produces the correct results, the phase-quenched results or neither in the weak-coupling limit. Other random matrix CLE simulations seem more optimistic [24]. Here it is emphasized that gauge-cooling is essential for obtaining the correct results.

Complex Langevin Equation for finite density Lattice QCD
If S (U) is the gauge action obtained from 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 l λ a . Here λ a are the Gell-Mann matrices for S U(3). η a l (t) are Gaussian-distributed random numbers normalized so that: The complex-Langevin equation has the same form except that the Us are now in S L (3, C) where M(U, µ) is the staggered Dirac operator. Note: backward links are represented by U −1 not U † . Note also that we have chosen to keep the noise-vector η real. η is gauge-covariant under S U(3), but not under S L (3, C). This means that gauge-cooling is non-trivial. Reference [7] indicates why this is not expected to change the physics. After taking −iδS (U, µ)/δU l , the cyclic properties of the trace are used to rearrange the fermion term so that it remains real for µ = 0 even after replacing the trace by a stochastic estimator, at least for infinite precision. To simulate the time evolution of the gauge fields we use the partial second-order formalism of Fukugita, Oyanagi and Ukawa. [25][26][27] After each update, we gauge-fix iteratively to a gauge which minimizes the unitarity norm: where V is the space-time volume of the lattice. This is referred to as gauge cooling [12].
3 Zero Temperature Simulations on a 16 4 lattice at β = 5.7 Since our earlier CLE simulations of lattice QCD at β = 5.6, m = 0.025 on a 12 4 lattice correctly reproduced some aspects of the expected phase structure, but failed in the details, we are repeating such simulations at weaker coupling (while continuing our β = 5.6 runs). We are running simulations at β = 5.7, m = 0.025. In order to be sure that the µ = 0 theory is confined, we have increased our lattice size to 16 4 . Important µ values are m π /2 ≈ 0.194 and m N /3 ≈ 0.28. These were obtained from results published by the Columbia group [28,29]. We are running 2 or 3 million updates at each µ, from µ = 0 up to saturation. Our input dt = 0.01, and we perform 5 gauge-cooling steps after each update. After adaptive rescaling of dt this gives us from ≈ 80 to over 1000 equilibrated time-units per µ.
A recent paper by Bloch et al. [23] has demonstrated that when CLE simulations of a random matrix theory at finite µ, which is a model for QCD at finite µ, converges to the wrong limit, it produces the results of the phase-quenched theory. We therefore run RHMC simulations of phasequenched lattice QCD on the same size lattice at the same β over the same range of µ values. (For a description of our methods and the physics of phase-quenched QCD see for example [30,31].) These results are plotted on the same graphs as the CLE data, for comparison. Note that for these phase-quenched simulations we need to introduce a small symmetry breaking parameter λ. We will eventually need to run at more than one choice of this parameter and extrapolate to λ = 0.
The physics one expects is that all observables should remain at their µ = 0 values up to the transition. For the full theory, this transition should occur at µ just below m N /3, while for the phasequenched theory it should occur at µ ≈ m π /2. For the the phase-quenched theory, the transition should be second order, while for the full theory it is expected to be first order. Beyond this phase transition the chiral condensate ( ψ ψ ) should fall rapidly approaching zero at saturation, if not before. The quark number density should remain zero up to the transition. Above this it should rise, approaching 3 at saturation. The plaquette should remain constant up to the transition. Above this it should increase towards its saturation value. It is expected that, when µ is sufficiently large, the phase of the fermion determinant will cease to be important, and the observables of the full and phase-quenched theories should be identical. Figure 1 shows the chiral condensate, ψ ψ , from our β = 5.7 CLE simulations of the full theory and the phase-quenched results. We note that there is good agreement at µ = 0, which was not true at β = 5.6. However the condensate from our CLE simulations starts to fall (almost) immediately µ is increased from zero disagreeing with both the expected behaviour and the phase-quenched simulations. However, there is some indication that it falls more slowly than at β = 5.6. The full (CLE) and phasequenched results converge as µ is increased until by µ = 0.5 they are statistically indistinguishable. Similar convergence is observed for β = 5.6. Such large µ behaviour agrees with expectations.   Figure 2 shows the quark-number density from these simulations both of the full theory and of the phase-quenched approximation. While there is good agreement between the full (CLE) and phase quenched results for very small µ, where both are consistent with zero, the 2 diverge for 0.1 < µ < 0.5. The difference is small numerically, only because this quantity is small in this region. For µ ≥ 0.5 the results for the 2 theories are statistically indistinguishable. This contrasts with the β = 5.6 case where a small but significant difference persists up to saturation 1 .
In figure 3 we plot the plaquette measured in these simulations. The agreement between the 2 theories at µ = 0 is not as good as for the other observables, however, it is significantly better than for β = 5.6. The plaquette stays close to its µ = 0 value through the transition region. Again, the full (CLE) and phase-quenched results converge and remain together for µ ≥ 0.5 whereas for β = 5.6 a small difference persists up to saturation.
In summary, we see some improvement in the behaviour of the CLE for β = 5.7 over that for β = 5.6. However, at this stage, we cannot rule out the weak-coupling limit of the CLE producing phase-quenched results rather than the correct results.
Finally in figure 4 we show the average unitarity norm as a function of µ from our simulations on a 16 4 lattice at β = 5.7 and from our 12 4 runs at β = 5.6. We note that the unitarity norms for β = 5.7 are significantly smaller than those for β = 5.6. It has been suggested that keeping unitarity norms small increases the chance that the CLE will converge to the correct results, since there is less chance that the trajectories will encounter zeros of the fermion determinant (poles in the drift term). this suggests one should try even weaker couplings. Note also that both these graphs indicate that the unitarity norm has a deep minimum around µ = 0.5 or 0.6. We suggest that the CLE simulations  might converge to the correct distributions for µ above or even slightly below these minima, noting that this µ is approximately that for which the β = 5.7 full and phase-quenched results converge.

Summary and Discussion
We simulate lattice QCD at finite quark-number chemical potential using the complex Langevin equation (CLE) to determine whether this is a viable method for evading the sign problems of this theory. Our simulations are performed at β = 5.6, m = 0.025 on a 12 4 lattice and at weaker coupling β = 5.7, m = 0.025 on a 16 4 lattice. We compare our observables with those obtained in the phase-quenched approximation, where one uses only the magnitude of the fermion determinant, using exact (RHMC) simulations. For µ < m π /2, both theories should agree and the observables should remain fixed at their µ = 0 values. For the full theory the observables should remain constant up to µ ≈ m N /3. Appreciably beyond µ = m N /3, the phase of the determinant should be small enough that both theories should be in agreement.
At µ = 0 the chiral condensate for β = 5.6 shows sizable departures from the known value, while at β = 5.7 there is good agreement with the exact result. In both cases, the condensate starts to fall (almost) immediately µ > 0, which does not agree with either the full or phase-quenched expectations. However, the falloff for β = 5.7 does appear to be slightly less rapid than at β = 5.6. In the transition region the observables for both β values disagree both with what is expected and with those of the phase-quenched approximation. At small µs through the transition region, the plaquettes for both β = 5.6 and β = 5.7 differ from the known result for µ = 0 by a small but significant amount. However, the β = 5.7 plaquette is closer to the correct value than is the the β = 5.6 value. For large µs the chiral condensates for both βs converge to their phase-quenched values. In this regime the quark-number density and plaquette for the β = 5.7 simulations converge to the phase-quenched results while for β = 5.6 small but significant differences persist.
Hence, at β = 5.6 and β = 5.7 the CLE fails to produce correct results, but some improvement is seen as we go to weaker coupling. However, the results are worse than those of the phase-quenched approximation, so although it is possible that in the weak-coupling limit, the CLE produces correct results, it is also possible that it produces phase-quenched results. Therefore we are extending our simulations to even weaker couplings (β = 5.8, β = 5.9 and possibly β = 6.0) on a 32 4 lattice.
We also plan to research other (fermion) actions for which the CLE might exhibit better convergence. Because of flavour-symmetry (taste) breaking, the zeros of the 2-flavour determinant have order 1 2 instead of 2. This is likely to produce problems when calculating the chiral condensate or quark-number density, both of which have poles at these zeros. (Is this why the plaquettes appear to behave better over the transition region than the chiral condensate and quark-number densities?) This suggests using Wilson fermions whose exact vector flavour symmetry means that the 2-flavour determinant will have zeros of order 2 as in the continuum. Another possible action is one with an irrelevant chiral 4-fermion term, which moves the (pseudo-)Goldstone pions to the auxiliary field and raises the mass of the fermion fields, which produce the fermion determinant and have all the µ dependence [32,33]. This can be thought of as replacing the current quarks with pions and constituent quarks. Other possibilities include modifying the gauge-fixing [24], or modifying the fermion dynamics by adding irrelevant terms to the drift term [34,35]. Both are designed to reduce the size of the unitarity norms.
Since the CLE appears to perform better in the nuclear matter (non-zero quark-number density) phase, it might be possible to answer questions such as if and under what circumstances such exotic phases as colour-superconductors can be produced, even when the CLE fails in the region of the transition from hadronic to nuclear matter.
We also plan to apply the CLE to simulate QCD at finite µ and temperature, T , since it has been observed that the CLE works better in this domain. Here we will study the transition from hadronic and nuclear matter to a quark-gluon plasma, and search for the critical endpoint.