# Results on the heavy-dense QCD phase diagram using complex Langevin

###### Abstract

Complex Langevin simulations have been able to successfully reproduce results from Monte Carlo methods in the region where the sign problem is mild and make predictions when it is exponentially hard. We present here our study of the QCD phase diagram and the boundary between the confined and deconfined phases in the limit of heavy and dense quarks (HDQCD) for 3 different lattice volumes. We also briefly discuss instabilities encountered in our simulations.

Results on the heavy-dense QCD phase diagram using complex Langevin

Dénes Sexty

Department of Physics, Bergische Universität Wuppertal,

Gaussstraße 20, D-42119 Wuppertal, Germany

and Inst. for Theoretical Physics, Eötvös University,

Pázmány P. sétány 1/A, H-1117 Budapest, Hungary

E-mail: sexty@uni-wuppertal.de

\abstract@cs

## 1 Introduction

For the past three decades countless non-perturbative results of QCD have been obtained using the lattice formulation, such as decay constants, hadron masses and properties of the thermal transition to the quark-gluon plasma at zero density. However, a full picture of the QCD phase diagram at finite temperature and chemical potential, sketched in Figure 1, requires research at non-zero density as well.

A finite quark chemical potential in the Euclidean path integral formulation leads to the so-called sign problem, i.e., to a complex probability weight. In some situations, where the sign problem is considered “mild”, standard Monte-Carlo based methods, such as reweighting, and properties related to analyticity, like Taylor expansion, can still be applied [1]. However, these techniques cannot probe reliably regions where due to an exponentially hard overlap problem.

This difficulty is due to the quark determinant, which appears once the quark fields have been integrated out, since it is complex for real chemical potential ,

(1.0) |

where generically represents the gauge links.

We present here our results on the phase diagram of QCD in the heavy-dense approximation, to be explained below, obtained using complex Langevin simulations. We also discuss instabilities found in the simulations. This contribution is based on [2].

## 2 Complex Langevin equation

Stochastic quantisation [3] is a procedure that allows quantum expectation values to be evaluated as averages over a stochastic process by evolving the dynamical variables in a fictitious time using a Langevin equation,

(2.0) |

The gauge links are represented by , are the Gell-Mann matrices, is the stepsize and are white noise fields satisfying

(2.0) |

The drift is derived from the QCD action, , and the gauge group derivative is defined as

(2.0) |

We have used an adaptive Langevin stepsize, based on the absolute value of the drift term , in order to avoid numerical instabilities [4]. Quantum expectation values are computed as averages over the Langevin time after the system reaches equilibrium.

The sign problem is potentially evaded by allowing the system to explore a larger configuration space [5, 8, 9, 10]. For an SU() gauge theory the group is enlarged to SL(). In this context the gauge action is generalized by replacing for , which keeps the gauge action holomorphic.

The SL() manifold is not compact, which means the system might follow a trajectory where the imaginary parts of the gauge fields are no longer small deformations compared to the real ones. A measure of the distance from the unitary manifold is given by

(2.0) |

with being the lattice four-volume and the equality only holding if all are unitary. The gauge cooling [11] technique is employed to prevent the system from going too far from SU(). It uses gauge transformations to push the system as close as possible to the unitary manifold. In other words,

(2.0) |

where , similar to , is changed adaptively based on the absolute value of [12]. These transformations seek configurations that are closest to the SU() submanifold and gauge equivalent to those generated in the Langevin process. Gauge cooling has been shown not to affect the equilibrium probability distributions of the observables [13].

Direct simulations of full QCD [14, 15] and comparisons with the hopping expansion to all orders [16] and multi-parameter reweighting [17] are among recent works involving the complex Langevin method. Discussions regarding the role of the pole arising from the logarithm of the (fermion) determinant in the complex Langevin process can be found in [18, 19, 20].

## 3 Phase diagram of heavy-dense QCD

The heavy-dense approximation of QCD (HDQCD) [21, 9] is obtained by dropping the spatial hopping terms for the fermions, while the remainder can be treated exactly. The HDQCD fermion determinant reads

(3.0) |

where is the temporal extent of the lattice, is the number of degenerate quark flavours, and the Polyakov loop and its inverse are given by

(3.0) |

This model exhibits the sign and Silver Blaze [22] problems, but is computationally cheaper to simulate than full QCD. Other studies include combinations with strong-coupling expansions [23, 24], reweighting [25], histogram [26], density of states methods [27] and mean-field methods [28]. Other complex Langevin studies of this model have led to important algorithmic improvements [9, 4, 11, 16].

Our study consisted of an extensive scan of the plane with a total of ensembles using different combinations of and for each of the three different lattice volumes, described in Table 1. The critical chemical potential indicates the transition to higher densities at zero temperature (). We have used fixed gauge coupling and hopping parameter , resulting in an approximate lattice spacing of fm, which has been determined using the gradient flow [29].

From the analysis of the distance from the SU() manifold and distributions of observables we have learned that the complex Langevin method combined with gauge cooling produces correct results as long as [2]. Therefore we present here simulation results for which the distance is less than .

fm | ||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

28 | 24 | 20 | 16 | 14 | 12 | 10 | 8 | 7 | 6 | 5 | 4 | 3 | 2 | |

[MeV] | 48 | 56 | 67 | 84 | 96 | 112 | 134 | 168 | 192 | 224 | 268 | 336 | 447 | 671 |

The main observable we have used in this work is the average of the traced (inverse) Polyakov loop, defined by

(3.0) | ||||

(3.0) |

It is worth noting that and are complex-valued but their expectation values are real, as they are proportional to with being the free energy of a single (anti) quark. We consider the symmetrised combination, which is real for SU() gauge configurations,

(3.0) |

In [2] results for the density and its susceptibility can be found as well.

The left panel of Figure 2 shows a D plot of the average symmetrised Polyakov loop as function of and for a volume of and the parameters given in Table 1. Each black point represents the average from an individual simulation and the cubic spline surface is plotted to guide the eye. The Polyakov loop exhibits both the transition to high densities, which is induced by the quark dynamics, and the thermal deconfinement transition, driven by gluonic dynamics.

For the determination of the boundary between the different phases we have studied the Binder cumulant [30] of ,

(3.0) |

The result is shown on the right hand side of Figure 2, where a clear separation between the confined, with , and deconfined, , phases is visible.

Figure 3 shows the deconfinement transition in the plane, along with a polynomial fit, for the lattice volumes of , and . The points and error bars represent the region where the Binder cumulant increases from to .

Finite size effects are clearly visible in Figure 3, but are less pronounced for the two larger volumes. The polynomials shown are of the form , which enforces that . Comparisons with higher order polynomials have revealed that this is sufficient for our data set. A study of different functions, including an expected non-analiticity at , along with the fitted values for the polynomials can be found in [2].

In principle, the Binder cumulant can be used to determine the order of the phase transition, since at the transition point its value depends only on the universality class [30]. However, due to our choice of fixed lattice spacing our values for the temperature were constrained by being an integer number. A deeper analysis of the volume dependence would require a more precise determination of the critical temperature as function of , with smaller error bars.

## 4 Instabilities

During our analysis we have found cases where after a certain Langevin time the distribution of the observables would change, typically becoming wider or shifting their average values. Based on the formal justification [10, 31] and comparisons with reweighting [25] we concluded that these instabilities do not reflect the original theory.

We have found evidence that a large distance from SU() is correlated with the appearance of the instabilities. Therefore, we have imposed a conservative cutoff on the distance, i.e., exclude data points where . A more detailed discussion is found in [2].

## 5 Conclusion

We have presented here our results on the phase diagram of QCD in the limit of heavy quarks, using complex Langevin simulations. The combination with gauge cooling allowed us to perform ab-initio simulations on the whole plane. We have mapped the phase boundary among different phases using the Binder cumulant for the symmetrised Polyakov loop and described the boundary line by a polynomial.

Despite some instabilities during the Langevin process, which happen when the system is too far from the unitary submanifold even after gauge cooling is applied, we have been able to collect enough data for a reliable analysis.

The prospect for the future would be to improve the precision in the determination of the phase boundary, and also the order of the transition. For that it will be necessary to change the temporal extent and the lattice spacing simultaneously. Furthermore, a better control of the Langevin process is important and work in this direction is under development.

## Acknowledgments

We thank Erhard Seiler and Ion-Olimpiu Stamatescu for collaboration on related topics. We are grateful for the computing resources made available by HPC Wales. We acknowledge the STFC grants ST/L000369/1, ST/K000411/1, ST/H008845/1, ST/K005804/1 and ST/K005790/1, the STFC DiRAC HPC Facility (www.dirac.ac.uk), STFC grant ST/L000369/1, the Royal Society and the Wolfson Foundation. FA is grateful for the support through the Brazilian government programme “Science without Borders” under scholarship number BEX 9463/13-5.

## References

- [1] P. de Forcrand, PoS LAT2009 (2009) 010, [1005.0539].
- [2] G. Aarts, F. Attanasio, B. Jäger and D. Sexty, JHEP 09 (2016) 087, [1606.05561].
- [3] G. Parisi and Y.-s. Wu, Sci. Sin. 24 (1981) 483.
- [4] G. Aarts, F. A. James, E. Seiler and I.-O. Stamatescu, Phys. Lett. B687 (2010) 154–159, [0912.0617].
- [5] J. R. Klauder, Acta Phys. Austriaca Suppl. 25 (1983) 251–281.
- [6] J. R. Klauder, J. Phys. A16 (1983) L317.
- [7] J. R. Klauder, Phys. Rev. A29 (1984) 2036–2047.
- [8] G. Parisi, Phys. Lett. B131 (1983) 393–395.
- [9] G. Aarts and I.-O. Stamatescu, JHEP 09 (2008) 018, [0807.1597].
- [10] G. Aarts, E. Seiler and I.-O. Stamatescu, Phys. Rev. D81 (2010) 054508, [0912.3360].
- [11] E. Seiler, D. Sexty and I.-O. Stamatescu, Phys. Lett. B723 (2013) 213–216, [1211.3709].
- [12] G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty and I.-O. Stamatescu, Eur. Phys. J. A49 (2013) 89, [1303.6425].
- [13] K. Nagata, J. Nishimura and S. Shimasaki, PTEP 2016 (2016) 013B01, [1508.02377].
- [14] D. Sexty, Phys. Lett. B729 (2014) 108–111, [1307.7748].
- [15] D. K. Sinclair and J. B. Kogut, PoS LATTICE2015 (2016) 153, [1510.06367].
- [16] G. Aarts, E. Seiler, D. Sexty and I.-O. Stamatescu, Phys. Rev. D90 (2014) 114505, [1408.3770].
- [17] Z. Fodor, S. D. Katz, D. Sexty and C. Török, Phys. Rev. D92 (2015) 094516, [1508.05260].
- [18] K. Splittorff, Phys. Rev. D91 (2015) 034507, [1412.0502].
- [19] A. Mollgaard and K. Splittorff, Phys. Rev. D91 (2015) 036007, [1412.2729].
- [20] J. Nishimura and S. Shimasaki, Phys. Rev. D92 (2015) 011501, [1504.08359].
- [21] I. Bender, T. Hashimoto, F. Karsch, V. Linke, A. Nakamura, M. Plewnia et al., Nucl. Phys. Proc. Suppl. 26 (1992) 323–325.
- [22] T. D. Cohen, Phys. Rev. Lett. 91 (2003) 222001, [hep-ph/0307089].
- [23] P. de Forcrand, J. Langelage, O. Philipsen and W. Unger, Phys. Rev. Lett. 113 (2014) 152002, [1406.4397].
- [24] J. Glesaaen, M. Neuman and O. Philipsen, JHEP 03 (2016) 100, [1512.05195].
- [25] R. De Pietri, A. Feo, E. Seiler and I.-O. Stamatescu, Phys. Rev. D76 (2007) 114501, [0705.3420].
- [26] S. Ejiri, Eur. Phys. J. A49 (2013) 86, [1306.0295].
- [27] N. Garron and K. Langfeld, [1605.02709].
- [28] T. Rindlisbacher and P. de Forcrand, JHEP 02 (2016) 051, [1509.00087].
- [29] S. Borsanyi et al., JHEP 09 (2012) 010, [1203.4469].
- [30] K. Binder, Z. Phys. B43 (1981) 119–140.
- [31] G. Aarts, F. A. James, E. Seiler and I.-O. Stamatescu, Eur. Phys. J. C71 (2011) 1756, [1101.3270].