# Exploring Complex-Langevin Methods for Finite-Density QCD

###### Abstract

QCD at non-zero chemical potential () for quark number has a complex fermion determinant and thus standard simulation methods for lattice QCD cannot be applied. We therefore simulate this theory using the Complex-Langevin algorithm with Gauge Cooling in addition to adaptive methods, to prevent runaway behaviour. Simulations are performed at zero temperature on a lattice with 2 quarks which are light enough that is significantly larger than . Preliminary results are qualitatively as expected. The quark-number density is close to zero for , beyond which it increases, eventually reaching its saturation value of for sufficiently large. The chiral condensate decreases as is increased approaching zero at saturation, while the plaquette increases towards its quenched value. We have yet to observe the transition to nuclear matter at , presumably because the runs for between and saturation have yet to equilibrate.

Exploring Complex-Langevin Methods for Finite-Density QCD

J. B. Kogut

Department of Energy, Division of High Energy Physics, Washington, DC 20585, USA

and

Department of Physics – TQHN, University of Maryland, 82 Regents Drive, College Park, MD 20742, USA

E-mail: jbkogut@umd.edu

\abstract@cs

## 1 Introduction

QCD at a finite chemical potential for quark number has a complex action which prevents the direct application of simulation methods based on importance sampling. The Langevin equation is a stochastic differential equation for the evolution of the classical fields in a fictitious time, which does not rely on importance sampling. It is, in fact, a special case of the hybrid molecular-dynamics algorithm, where each trajectory consists of a single update.

The Langevin equation can be extended to complex actions by complexifying the fields [1, 2, 3, 4]. In the case of QCD this means promoting the gauge fields from to . Unfortunately, there is no proof that the long-time evolution of the fields under this Complex Langevin equation (CLE) provides a limiting value for observables. Even when this process does converge, the values it provides for observables are not guaranteed to be correct.

After successfully applying the CLE to spin models (see for example [5]), people were encouraged to apply it to lattice QCD at finite . Early attempts at applying the CLE to QCD were stymied by runaway behaviour, which was not corrected by adaptive methods. Recently it has been noted that at least some of this undesirable behaviour is due to the production of unbounded gauge transformations of compact gauge fields. Such behaviour can be controlled by gauge transforming to a gauge which minimizes the magnitudes of the gauge fields and hence their distance from the manifold[6]. This is called Gauge Cooling.

This has revived interest in the CLE for QCD at finite . These methods have been tested on simple models and for quark masses large enough that hopping-parameter methods can be applied [7, 8, 9, 10, 11, 12]. In addition, studies have been made of the conditions under which the CLE converges to the correct results [13, 14, 15, 16, 17, 18, 19]. QCD at finite and small masses has been simulated and the results compared with the heavy quark methods, for larger masses [20]. Very recently this has been extended to larger lattices at finite temperatures, where the transition from hadron/nuclear matter to a quark-gluon plasma is observed and results are compared with those from reweighting methods [21].

We are simulating QCD at zero temperature and finite for light quarks using the CLE, to test directly if it converges and produces believable results. We present preliminary results of our explorations.

## 2 Complex Langevin for finite density Lattice QCD

If is the gauge action after integrating out the quark fields, the Langevin equation for the evolution of the gauge fields in Langevin time is:

(2.0) |

where labels the links of the lattice, and . Here are the Gell-Mann matrices for . are Gaussian-distributed random numbers normalized so that:

(2.0) |

The complex-Langevin equation has the same form except that the s are now in . , now is

(2.0) |

where is the staggered Dirac operator. Note: backward links are represented by not . Note also that we have chosen to keep the noise-vector real. is gauge-covariant under , but not under . This means that gauge-cooling is non-trivial. Reference [18] indicates why this is not expected to change the physics. After taking , the cyclic properties of the trace are used to rearrange the fermion term so that it remains real for even after replacing the trace by a stochastic estimator.

To simulate the time evolution of the gauge fields we use the partial second-order formalism of Fukugita, Oyanagi and Ukawa. [22, 23, 24]

After each update, we gauge-fix iteratively to a gauge which minimizes the unitarity norm – gauge cooling [6]:

(2.0) |

where is the space-time volume of the lattice.

## 3 Zero-temperature simulations

### 3.1

For and infinite precision, Complex Langevin becomes Real Langevin. At 64-bit precision, roundoff allows the gauge fields to move (slowly) off the manifold. For , on an lattice we observe runaway solutions, even after Gauge Cooling!

For , on a lattice without gauge cooling, we observe runaway solutions. With gauge cooling, the trajectory moves slowly off the manifold. We perform 100,000 updates with input . Adaptively rescaling to keep the drift(force) term under control is reduced to , so the total run takes Langevin time units, at the end of which the unitarity norm . This we can probably tolerate, especially since we expect it to improve with weaker couplings and larger lattices. Figure 2 shows the time evolution of the unitarity norm with and without gauge-cooling.

For this run plaquette , whereas the RHMC algorithm gives , and , RHMC . This is reasonable agreement for a short run with an inexact algorithm.

### 3.2

We simulate on a lattice at , with . Potentially important values include and (masses from HEMCGC collaboration [25, 26, 27]). The first is the position of the expected transition for the phase-quenched approximation. The second is the approximate position of the expected transition to nuclear matter. We start with a limited number of values to probe the various regions of the zero-temperature phase diagram. The values we choose are , , , , , and . In each case we start the simulation from an ordered start and use 5-step gauge-cooling.

The first thing we look for, is evidence that the trajectories for a given set of parameters are restricted to a compact region of the manifold. Without this it is (almost) impossible for these simulations to produce meaningful results. If the simulations do converge to a limiting distribution, one must then address the question as to whether this is the correct limit.

At we have performed sufficient updates for the unitarity norm to level off. The unitarity norm appears to have leveled off, indicating that the system is evolving over a compact domain of The evolution of this norm over the trajectory is shown in figure 2. The quark number density . Hence the system has reached saturation where , as expected for large . This is where each site is occupied with 3 quarks in a colour singlet state (nucleon). The chiral condensate – small and consistent with zero as expected. The plaquette , consistent with the idea that, at saturation, the quarks are frozen out and the system approximates quenched QCD. The quenched plaquette at , . The total trajectory length time-units. In fact, after equilibration, .

To date, of the other values we are simulating, and appear to have equilibrated. appears to be close to equilibrating. , and have yet to equilibrate. We believe that this is because these runs need more updates, since their unitarity norms have yet to reach values achieved for the simulations.

We present ‘data’ for the quark-number densities (figure 4), and chiral condensates (figure 4) as functions of with the understanding that the points at , and are expected to change to become closer to the values at as the system equilibrates. This is because, since we start each run with all gauge links on the manifold, then as the system equilibrates, the gauge links move away from this manifold. Because not all s are equilibrated, we have not included error-bars in these figures.

These preliminary results (each point represents 90,000 – 500,000 sweeps/updates of the lattice) agree qualitatively with our expectations. The quark-number density remains close to zero for . For larger s it becomes non-zero, increasing towards its saturation value of as is increased. The chiral condensate decreases monotonically from its value as is increased, approaching zero at saturation. We will need to wait until each point has equilibrated to where it is clear that the gauge fields are varying over a compact region in the manifold, before we can get truly quantitative results. This takes longer for since it takes more iterations to invert the Dirac operator as increases, until close to saturation, and is smaller. Since each run is starting from the manifold, we expect metastability for due to the presence of a supposed first-order transition at . This will also slow down the approach to equilibrium just above the transition.

Because the runs for just above the transition have yet to equilibrate, we have been unable to observe this transition to nuclear matter.

## 4 Summary, discussion and outlook

We apply Complex-Langevin simulations with gauge cooling to lattice QCD at finite quark-number chemical-potential () at zero temperature. Our current simulations are on a lattice with , , . Preliminary results look promising, but more simulations are needed. Adaptive updating with gauge cooling does appear to stabilize the algorithm. However, this only appears to work provided the gauge coupling is not too strong. Since gauge cooling does not completely fix the gauge – the unitarity norm is invariant under gauge transformations – it is possible that further gauge fixing might improve the situation. Fixing to Landau gauge in the subgroup suggests itself.

We need answers to the following. Do these simulations converge and converge to the correct limit? Do we observe a phase transition to nuclear matter at ? Is there a spurious transition at ? Do these simulations produce the expected 2-flavour colour-superconductor at large ()?.

Smaller masses are needed – (?). We also need larger lattices, weaker coupling… Our current code is inefficient serial code which needs improvement. When we are convinced that the algorithm works we will parallelize our code.

We will also investigate whether we can make a fully second-order version. In addition we will investigate whether it makes sense to make a complex extension of hybrid molecular-dynamics. Would complex hybrid molecular-dynamics be expected to converge to the correct limit, if the complex langevin does? If so, such an algorithm would expect to be faster and have smaller errors than the complex-langevin methods.

A precise measurement of the value of at the transition to nuclear matter () would yield the binding energy/nucleon () in the absence of electromagnetic interactions, since . However, since , this will be a formidable task. More accessible nuclear physics will be to study the propagation of hadrons in the nuclear-matter medium. If our 2-flavour simulations are successful, we will also simulate and try to observe the 3-flavour colour-superconductor with its colour-flavour locking.

We also plan to simulate at high temperatures, near to the finite-temperature phase transition, and look for the critical endpoint. We will try to determine how good is the resonance-gas model. Also planned are investigations of the phase structure of -flavour QCD with independent chemical potentials for the and quarks.

## Acknowledgements

These simulations are being performed on PCs belonging to Argonne’s HEP Division, Edison and Carver at NERSC, and Blues at LCRC, Argonne.

## References

- [1] G. Parisi, Phys. Lett. B 131, 393 (1983).
- [2] J. R. Klauder, Acta Phys. Austriaca Suppl. 25, 251 (1983).
- [3] J. R. Klauder, J. Phys. A 16, L317 (1983).
- [4] J. R. Klauder, Phys. Rev. A 29, 2036 (1984).
- [5] F. Karsch and H. W. Wyld, Phys. Rev. Lett. 55, 2242 (1985).
- [6] E. Seiler, D. Sexty and I. O. Stamatescu, Phys. Lett. B 723, 213 (2013) [arXiv:1211.3709 [hep-lat]].
- [7] G. Aarts and I. O. Stamatescu, JHEP 0809, 018 (2008) [arXiv:0807.1597 [hep-lat]].
- [8] G. Aarts, Phys. Rev. Lett. 102, 131601 (2009) [arXiv:0810.2089 [hep-lat]].
- [9] G. Aarts and F. A. James, JHEP 1201, 118 (2012) [arXiv:1112.4655 [hep-lat]].
- [10] G. Aarts and K. Splittorff, JHEP 1008, 017 (2010) [arXiv:1006.0332 [hep-lat]].
- [11] G. Aarts, F. Attanasio, B. Jäger, E. Seiler, D. Sexty and I. O. Stamatescu, PoS LATTICE 2014, 200 (2014) [arXiv:1411.2632 [hep-lat]].
- [12] G. Aarts, E. Seiler, D. Sexty and I. O. Stamatescu, Phys. Rev. D 90, no. 11, 114505 (2014) [arXiv:1408.3770 [hep-lat]].
- [13] G. Aarts, F. A. James, E. Seiler and I. O. Stamatescu, Eur. Phys. J. C 71, 1756 (2011) [arXiv:1101.3270 [hep-lat]].
- [14] G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty and I. O. Stamatescu, Eur. Phys. J. A 49, 89 (2013) [arXiv:1303.6425 [hep-lat]].
- [15] 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]].
- [16] G. Aarts, P. Giudice and E. Seiler, Annals Phys. 337, 238 (2013) [arXiv:1306.3075 [hep-lat]].
- [17] J. Nishimura and S. Shimasaki, Phys. Rev. D 92, no. 1, 011501 (2015) [arXiv:1504.08359 [hep-lat]].
- [18] K. Nagata, J. Nishimura and S. Shimasaki, arXiv:1508.02377 [hep-lat].
- [19] H. Makino, H. Suzuki and D. Takeda, arXiv:1503.00417 [hep-lat].
- [20] D. Sexty, Phys. Lett. B 729, 108 (2014) [arXiv:1307.7748 [hep-lat]].
- [21] Z. Fodor, S. D. Katz, D. Sexty and C. Török, arXiv:1508.05260 [hep-lat].
- [22] A. Ukawa and M. Fukugita, Phys. Rev. Lett. 55 (1985) 1854.
- [23] M. Fukugita, Y. Oyanagi and A. Ukawa, Phys. Rev. D 36 (1987) 824.
- [24] M. Fukugita and A. Ukawa, Phys. Rev. D 38, 1971 (1988).
- [25] K. M. Bitar et al., Phys. Rev. D 42, 3794 (1990).
- [26] K. M. Bitar et al., Phys. Rev. D 49, 6026 (1994) [hep-lat/9311027].
- [27] K. M. Bitar et al., Phys. Rev. Lett. 65, 2106 (1990).