The QCD phase diagram in the limit of heavy quarks using complex Langevin dynamics

Complex Langevin simulations allow numerical studies of theories that exhibit a sign problem, such as QCD, and are thereby potentially suitable to determine the QCD phase diagram from first principles. Here we study QCD in the limit of heavy quarks for a wide range of temperatures and chemical potentials. Our results include an analysis of the adaptive gauge cooling technique, which prevents large excursions into the non-compact directions of the SL($3, \mathbb{C}$) manifold. We find that such excursions may appear spontaneously and change the statistical distribution of physical observables, which leads to disagreement with known results. Results whose excursions are sufficiently small are used to map the boundary line between confined and deconfined quark phases.


JHEP09(2016)087
correctly (in the presence of heavy quarks). We therefore have the opportunity to map the phase boundary, connecting the thermal transition at high temperature with the onset transition, where the quark density becomes non-zero, at large chemical potential. As such, it is a good test scenario to prepare for the realistic case of QCD with light quarks.
The CL method consists of stochastic explorations of a complexified configuration space, without the requirement of a positive weight [21][22][23][24]. It is precisely the method's complex nature that allows for a solution of the sign problem, even when it is severe [15,[25][26][27]. However, success is not guaranteed [28][29][30][31][32][33][34] and convergence to a wrong limit may occur. Based on the theoretical justification of the approach [35,36], these cases of incorrect convergence can be identified a posteriori. Here we employ the adaptive gauge cooling technique [17,37], which is necessary but not sufficient to avoid convergence to wrong limits. In addition, in the presence of a fermion determinant, the drift appearing in the CL equation is no longer holomorphic, which requires a reconsideration of the justification [35,36] and may lead again to incorrect convergence in practice [38][39][40]. However, all indications are that this is not an issue for the model considered in this paper [17,18]. We remark that applications of CL to full QCD can be found in refs. [18,41] and a comparison with multiparameter reweighting in ref. [42].
This paper proceeds as follows: in section 2 we review the complex Langevin method applied to lattice QCD and the gauge cooling technique. Section 3 presents the heavy dense (HD) approximation of QCD and lists the parameters used in our simulations. In section 4 we present results for observables related to the Polyakov loop and quark density as functions of the temperature and chemical potential as well as the resulting phase diagram. Issues related to instabilities and their relation to excursion into the non-compact directions of the Langevin equations are discussed in section 5. In section 6 we present a conclusion and an outlook for future work. Preliminary results have appeared in refs. [43][44][45][46].

Complex Langevin equation and gauge cooling
We consider QCD in the grand-canonical formulation, where the (quark) chemical potential µ couples to quark number. For an elementary introduction, see e.g. ref. [47]. After integrating out the bilinear quark fields, the partition function is written as where S YM is the Yang-Mills action, U are the gauge links, and M is the fermion matrix, depending on the chemical potential and the gauge links. Quantum expectation values can be evaluated using Langevin dynamics in a procedure known as stochastic quantisation [48]. In this scheme, expectation values are obtained as averages over a stochastic process by evolving dynamical variables over a fictitious time θ. Importantly, importance sampling does not enter in this formulation. On the lattice, for an SU(3) gauge theory with links U x,ν , a Langevin update, using a first-order discretisation in the Langevin time θ = n , reads [49] where λ a are the Gell-Mann matrices (with Tr λ a λ b = 2δ ab , the sum over a = 1, . . . , 8 is assumed) and η a x,ν are Gaussian white noise fields, which satisfy η a x,µ η b y,ν = 2 δ xy δ ab δ µν .
The dynamics is governed by the action S, which generates the drift where D a x,ν is the gauge group derivative The quark contribution leads to poles in the drift, namely where det M = 0 and M −1 does not exist. In some cases this affects the results negatively [38,40], but in HDQCD this is not the case, as far as is understood [18,39]. In order to avoid numerical instabilities and regulate large values of the drift, it is necessary to change the Langevin stepsize ε adaptively [16], based on the absolute value of the drift term K a x,ν . In theories that exhibit the sign problem the drift is complex, resulting in an exploration of a larger configuration space. This is how the sign problem is potentially evaded [15,[21][22][23][24]35]. In an SU(3) gauge theory, this procedure enlarges the gauge group to SL(3, C). The latter group, however, is not compact. Parametrising the gauge links as this implies that the gauge fields A a x,ν can now assume complex values. The extra degrees of freedom can lead to trajectories in which the imaginary parts of the gauge fields are not a small deformation. A measure of the distance from the unitary manifold can be given by unitarity norms etc., where Ω = N τ N 3 s is the four dimensional simulation volume. These norms are invariant under unitary gauge transformations, but not under general SL(3, C) transformations. They are exactly zero only if all links U x,ν are unitary.
Gauge cooling [17] is a procedure to reduce the distance to the unitary manifold via SL(3, C) gauge transformations. It consists of a sequence of gauge transformations which decrease the unitary norms d i in a steepest descent fashion Note that f a x is obtained via an infinitesimal gauge transformation of d 1 . In order to optimise the cooling procedure, the coefficient α is changed adaptively based on the absolute value of f a x [37]. Cooling is also stopped once the rate of change of the unitary norm is below a set target.

JHEP09(2016)087 3 Heavy dense QCD
We consider the heavy dense approximation of QCD (HDQCD) [9,15], in which the gluonic action is the standard Wilson Yang-Mills lattice action, while in the quark action spatial hopping terms are neglected but all chemical potential dependence, which resides in the temporal hopping terms, is retained. As mentioned earlier, the gluonic dynamics is contained without approximation.
The action S then consists of the gluonic term, x,ν is the standard plaquette and β the lattice gauge coupling, and minus the logarithm of the quark determinant in the HD approximation. The latter is obtained from the standard Wilson fermion action, by dropping the spatial hopping terms, such that where Γ ±ν = (1 ± γ ν )/2. Taking the determinant in Dirac space and in spacetime indices yields, for a single quark flavour (below we consider N f = 2 degenerate quarks), The power 2 originates from the gamma-matrix structure and the + sign from the antiperiodic boundary conditions. In this expression P (−1) x are the (inverse) Polyakov loops, 5) with N τ the number of time slices in the temporal direction. The temperature T is related to N τ via T = 1/(aN τ ), with a the lattice spacing. The parameter h = (2κ) Nτ , with κ the hopping parameter, arises from the hopping expansion and, finally, N f is the number of quark flavours. Important observables are the expectation value of the traced (inverse) Polyakov loops and the quark density, defined by Tr P x , (3.6)  Table 1. Parameters used in this study. The chemical potential µ is varied from 0 to 1.3µ 0 c , with µ 0 c = − ln(2κ). The lattice spacing is set using the gradient flow [50] and is approximate.
with V = N 3 s the spatial volume, and [15] (3.9) Here we used the notation We note here that P x and P −1 x are complex-valued for a given gauge configuration but that their expectation values are real, as they are related to the free energy of a single (anti) quark. Below we will also consider the symmetrised combination which is real for each SU(3) gauge link configuration. It is useful to consider the zero-temperature limit, N τ → ∞, at fixed lattice spacing. We take µ > 0 and first look at the density. The contribution from the anti-quarks, i.e. the second term in eq. (2.7), is exponentially suppressed. For the quark contribution, we write z as where we used that µ/T = µN τ , with µ expressed in lattice units after the equality sign. We see therefore that at zero temperature the density vanishes when µ < µ 0 c (Silver Blaze region [8,47]) and equals saturation density (n sat = 6N f ) when µ > µ 0 c , irrespective of the value of the Polyakov loop. Hence µ 0 c is the critical chemical potential for onset at T = 0, but the behaviour in the region µ > µ 0 c is a lattice artefact. For the Polyakov loop, we similarly note that at zero temperature and µ < µ 0 c , the quarks do not couple to the gauge fields and hence P = 0, as in the pure gauge theory, while when µ > µ 0 c , P has to be zero as well to ensure a finite determinant. Hence at T = 0, P = 0 for all µ, except possibly at µ = µ 0 c . The vanishing of P above onset is again due to the maximal number of quarks that can be placed on a finite lattice. For more discussion of these aspects, see e.g. ref. [19].
Simulation parameters are listed in table 1. In order to scan the phase diagram, a wide range of temperatures and chemical potentials is covered and a total of 880 ensembles with different combinations of N τ and µ were generated, for each of the three volumes. We use a fixed gauge coupling throughout this work, β = 5.8, and the estimate of the lattice spacing of a ∼ 0.15 fm has been obtained using the gradient flow [50]. Using a fixed lattice spacing yields an adequate coverage of the phase diagram at low temperature, with fixed lattice artefacts, but a poorer coverage at larger temperature.

Phase diagram
We have performed an extensive scan of T -µ plane, to determine the phase structure by direct simulation. 1 In order to track the reliability [35,36] of the results we measured the unitarity norms (2.7), studied the distributions of observables, and compared with results obtained with reweighting [12], where applicable. From this analysis, we inferred that complex Langevin dynamics in combination with gauge cooling produces correct results, provided that the unitarity norm does not become too large, d 2 O(0.1). In light of these observations we present here only simulation data for which the unitarity norm is smaller than 0.03. In this regime we can extract physical information on the phase boundary of HDQCD. We come back to larger unitarity norms in section 5. Figure 1 shows the quark density n and the symmetrised Polyakov loop P s = 1 2 P +P −1 as functions of the temperature and chemical potential on the spatial volume of 10 3 . The plotted surfaces are cubic splines to guide the eye and each black point represents the average from an individual simulation. Other parameters are given in table 1. We have used the lattice spacing of a ∼ 0.15 fm to convert the temperature to physical units and expressed the chemical potential in terms of µ 0 c . The Polyakov loop shows both the thermal deconfinement transition, driven by gluonic dynamics, and the transition to high densities, driven by quark dynamics. The region where µ > µ 0 c is a lattice artefact and the Polyakov loop drops again to zero at low temperature, as explained above. At higher temperature, the Polyakov loop is nonzero for all chemical potentials. At low temperature the quark density rises sharply at µ = µ 0 c to saturation density (n sat = 12). This behaviour is smoothened out at higher temperature. The density only rises slowly as µ increases from zero; for heavy quarks, the quark number susceptibility at µ = 0 is exponentially suppressed. Figure 2. Susceptibility of the quark density n (left) and symmetrised Polyakov loop 1 2 P +P −1 (right). In both cases peak heights have been cut, resulting in white plateaus. Figure 2 shows the susceptibilities for the aforementioned observables, which outline the corresponding transitions. Note that in both cases dominant peaks are not shown, to improve visibility. In principle the phase boundary can be determined from these susceptibilities. A better signal, however, is obtained by employing the Binder cumulant B [51], which for an observable O is defined as

JHEP09(2016)087
(4.1) Let O be zero in one phase and nonzero in another, and assume that the higher moments are governed by Gaussian fluctuations. It is then easy to see that where in the latter case it is assumed that The Binder cumulant for the symmetrised Polyakov loop expectation value P s is shown in figure 3. The separation between the confined phase, with P s = 0, and deconfined phase, with P s = 0, is clearly visible. At low temperature, the transition can easily be identified, due to the adequate coverage of the parameter space and the relatively sharp transition. At higher temperature, the setup with fixed lattice spacing does not have sufficient resolution to determine the thermal transition with precision. Nevertheless, a clear phase boundary is seen to emerge. To identify the transition between both phases, we determine the parameters for which the Binder cumulant reaches 1/3, and the results are shown in figure 4. The uncertainties are estimated by taking half the distance between neighbouring points in both T and µ directions. As mentioned above, the resolution in the temperature direction is limited due to having only integer N τ values, which leads to large discretisation effects for the thermal transition. The transition to higher densities can be mapped out with much more precision. To parameterise the transition temperature as a function of the chemical potential, we have fitted the estimates for T c (µ) to a number of fitting functions. Using the notation we considered an expansion around x ∼ 0, i.e., fit A: where we used that T c (µ) is an even function of µ [52]. Given that due to the lattice setup the transition is better determined around x 1 than around 0, and that T c (µ 0 c ) = 0, we  Table 2. Fit parameters and reduced χ 2 for fits A and B, see eqs. (4.4), (4.5), used to describe the chemical potential dependence of the transition temperature, T c (µ), for three spatial volumes.
have considered a power series around x = 1 as well, namely fit B: The expansion parameters {a k } and {b k } are trivially related, provided that k a k = 0 emerges from the fit. Finally, to take into account nonanalytic behaviour around x = 1, as required by the Clausius-Clapeyron relation (∂T c (µ)/∂µ → ∞ at µ = µ 0 c ), we included one additional term and used fit C: with 0 < α < 1.
The left panel of figure 4 shows fits A, B and C with n = 2 for our largest volume of V = 10 3 . The non-analytic behaviour at µ ≈ µ 0 c is not evident from our data; our lowest temperature is still away from 0. Hence treating α in fit C as a fit parameter does not yield additional information on the transition line and we do not consider C any further. Fits A and B are seen to be compatible with each other, indicating that T c (µ 0 c ) = 0 emerges without imposing it. The fit coefficients and the corresponding reduced χ 2 can be found in table 2 for fit A and B, for n = 2. Note that a rough estimate for T c (µ = 0) in MeV is given by a 0 ∼ b 1 + b 2 , which sets the scale of the coefficients.
On the right-hand side of figure 4 we compare three different polynomials for fit B with n = 2, 3 and 4, again for V = 10 3 . Higher-orders polynomials result in an almost identical curve as the fourth-order polynomial fit (n = 2). Hence adding more parameters does not result in an improved fit. Fits B with n = 2 for all three volumes studied here are shown in figure 5. We observe clear finite-size effects, especially for the smallest simulation box (6 3 ). A much smaller trend can be seen in the two larger volumes. The main limitation, however, comes from the discretisation at high temperature, as discussed above.
The Binder cumulant is in principle suitable to determine the order of the phase transition, as its value at the transition point only depends on the universality class [51]. Further analyses of the volume dependence would, however, require a more precise determination of T c as a function of µ throughout the phase diagram, with smaller uncertainties.

Instabilities
In our simulations we encounter instabilities, complicating the analysis. These result in a widening of the distribution of observables during the Langevin process and affect susceptibilities and other quantities significantly. Based on the formal justification [35,36] and a comparison with reweighting [12], one can conclude that the wider distributions do not reflect the original theory. In this section we describe some of these features.
In figure 6 we show examples of the Langevin time evolution of the real part of the Polyakov loop P and the unitarity norm d 2 , at low (left) and high (right) temperature. We observe two distinct segments, characterised by a small unitarity norm and controlled fluctuations in the initial part, followed by larger fluctuations and unitarity norm afterwards. At the higher temperature, this also leads to a tunnelling transition for the Polyakov loop, from around 0.2 to 0.

JHEP09(2016)087
N τ = 20, µ = 0.5 100 < θ < 250 330 < θ < 500 Reweighting N τ = 4, µ = 0.7 20 < θ < 60 100 < θ < 500 Reweighting  The data in figure 6 is analysed further in table 3. We have determined expectation values for the Polyakov loop, and its susceptibility χ P and Binder cumulant in each of the two intervals where the Polyakov loop fluctuates consistently around a certain value. Results obtained with reweighting are shown as well. We note that the observables are, within the statistical error, in agreement with the latter in the first interval, but not in the second one. The apparent agreement B ∼ 0 at low temperature for the entire interval mostly reflects that P ∼ 0 throughout, and hence the susceptibility is a more sensitive measure of accuracy. In figure 7 we compare histograms for both scenarios. A Gaussian fit is added to guide the eye. For the region with larger unitarity norms, the distribution is broader, with a larger tail. At high temperature, there is in addition a shift of the mean. We conclude that the region with smaller unitarity norm leads to acceptable results, while those with a larger value do not.
The behaviour described above has been seen for different chemical potentials and temperatures, but in all cases widening of the distributions coincided with a severe change  in unitarity norm. We have checked that using smaller stepsizes does not prevent these transition from occurring. The inability to control the unitarity norm on coarser lattices was already noted in ref. [17].
To check the behaviour closer to the continuum limit, we have performed additional simulations with larger gauge coupling, β = 6.0 and 6.2. Figure 8 shows the real part of the Polyakov loop for an identical setup as in figure 6 and table 3. Simulations at low temperature (N τ = 20) are shown on the left and at high temperature (N τ = 4) on the right. On the finer lattices and at low temperature, the unitarity norms remain practically 0 for the entire simulation time. At the higher temperature, the unitarity norm still rises, but with a smaller exponent. Once the unitarity norm becomes too large, fluctuations become significantly larger and skirts emerge, as in the case discussed above. This behaviour can be seen in figure 9, which shows the histograms for the high-temperature runs for the two larger β values on a 10 3 lattices. Hence we conclude that the instabilities are still present on finer lattices, but that they set in later (at high temperature) or only appear beyond the length of the Langevin trajectory (at low temperature).
In order to maintain the volume of the lattice in physical units, simulations with a larger gauge coupling require larger simulation volumes to compensate the smaller lattice spacing. In the preceding section, we have found that employing a gauge coupling of β = 5. 8 Figure 9. Histograms of the real part of the Polyakov loop before and after the rise of the unitarity norm, for the larger gauge couplings of β = 6.0 (left) and 6.2 (right), at high temperature (N τ = 4, µ = 0.7). a compromise between simulation costs and the ability to extract reliable information on the phase boundary of HDQCD. The Langevin time when the unitarity norm starts rising varies considerably for different setups. In most cases that happens sufficiently after the thermalisation stage, which leaves enough data points to allow us to extract observables and perform a subsequent analysis. However, since the amount of available data suitable for analysis differs greatly between ensembles, we find different uncertainties in each setup, including the integrated auto-correlation time [53]. In order to implement these findings, we have made sure that for the results presented in the previous section, only simulation data for which the unitarity norm d 2 is smaller than 0.03 were included.
To check the sensitivity with respect to changes in the cutoff imposed on the unitarity norm d 2 and the robustness of physical observables, we show in figure 10 the dependence of the quark density and the Polyakov loop on the maximally allowed unitarity norm, for N τ = 4 and µ = 0.7 on a 10 3 lattice. We observe that the obervables are stable and JHEP09(2016)087 independent of the cutoff over a wide range, up to d 2 ∼ 0.5. A clear transition is visible for larger values of the cutoff, which coincides with the widening of the distributions, discussed above. Note that the change in the statistical uncertainties at larger d 2 cutoff follows from the decrease of data points available in the analysis. Similar behaviour is seen at other parameter values. We conclude that observables are robust under changes in the cutoff imposed on the unitarity norm. In the previous section we have conservatively chosen a small value for the cutoff, i.e. d 2 < 0.03, to stay sufficiently away from the wrong behaviour observed for unitarity norm of O(1).

Conclusion
We have studied the phase diagram of QCD in the presence of heavy quarks, using complex Langevin simulations. Combining gauge cooling with a careful monitoring of the Langevin process, we have shown that it is possible to perform ab-initio simulations in the entire T -µ plane. The phase boundary between the confined and the deconfined phases was determined via the Binder cumulant of the symmetrised Polyakov loop and the resulting line can be fitted in terms of simple polynomials. In our setup, in which the lattice spacing is fixed and temperature is varied by changing the temporal extent, the main uncertainty occurs at high temperature, where discretisation effects are severe. The transition at low temperature, however, can be determined with more precision.
During the Langevin process, we observed instabilities, which take configurations far away from the SU(3) submanifold, even in the presence of gauge cooling. These events lead to incorrect convergence. By monitoring the unitarity norm, we found that it is nevertheless possible to collect sufficient simulation data which can be used in a reliable manner. There are strong indications that this situation will improve on finer lattices.
As an outlook, we note that in order to determine the phase boundary, and the order of the transition, throughout the T -µ plane with more precision, it will be necessary to vary both the lattice spacing and the temporal extent simultaneously, both in the model considered here as in full QCD. Besides this, an important additional step is a better control on the Langevin process and work in this direction is currently under development.