Dynamical stabilisation of complex Langevin simulations of QCD

The ability to describe strongly interacting matter at finite temperature and baryon density provides the means to determine, for instance, the equation of state of QCD at non-zero baryon chemical potential. From a theoretical point of view, direct lattice simulations are hindered by the numerical sign problem, which prevents the use of traditional methods based on importance sampling. Despite recent successes, simulations using the complex Langevin method have been shown to exhibit instabilities, which cause convergence to wrong results. We introduce and discuss the method of Dynamic Stabilisation (DS), a modification of the complex Langevin process aimed at solving these instabilities. We present results of DS being applied to the heavy-dense approximation of QCD, as well as QCD with staggered fermions at zero chemical potential and finite chemical potential at high temperature. Our findings show that DS can successfully deal with the aforementioned instabilities, opening the way for further progress.


Introduction
Strongly interacting matter at finite baryon number density and temperature has been, and remains, an active research subject to understand QCD under extreme conditions. Features of QCD are typically studied in thermodynamic equilibrium, where the theory has two external parameters: the temperature T and and baryon chemical potential μ B . Varying those allows the exploration of the QCD phase diagram in the T -μ B plane. Known phases include ordinary nuclear matter and the quark-gluon plasma (QGP), with a colour superconducting phase expected at large μ B . Of great appeal are also the boundaries that mark the transition between these phases. This phase diagram has a fascinating structure, which is of a e-mail: pyfelipe@uw.edu b e-mail: jaeger@cp3.sdu.dk significance for the study of hot and/or dense systems, such as the early universe and heavy-ion collisions.
Heavy-ion collisions have been successfully used to investigate the high temperature behaviour of QCD at the relativistic heavy ion collider (RHIC) and the large hadron collider (LHC). These facilities, together with future ones, namely the Facility for Antiproton and Ion Research (FAIR) and the nuclotron-based ion collider facility (NICA), will further explore the phase diagram of QCD. They will allow the study of hadronic interactions under extreme conditions, such as higher baryonic density or very high temperatures.
From a theoretical perspective, some insight, at high temperature or density, can be gained from perturbation theory. A full picture of the phase diagram, however, requires non-perturbative methods. Recent lattice results at non-zero temperature include [1,2]. Typically, lattice QCD simulations at finite baryon/quark density are carried out using the grand canonical ensemble, with the chemical potential introduced as conjugate variable to the appropriate number density (quark, baryon, etc). At finite quark chemical potential, the simulations have to overcome the infamous sign problem-a complex weight in the Euclidean path integral. This imposes severe limitations on the applicability of standard numerical methods [3,4]. Many approaches to deal with the sign problem have been proposed, including the complex Langevin method [5][6][7][8], strong coupling expansions [9][10][11], Lefschetz thimbles [12][13][14][15][16], holomorphic gradient flow [17], density of states [18][19][20][21] and sign-optimized manifolds [22].
The complex Langevin (CL) method is an extension of the stochastic quantisation technique [23] to a complexified configuration space, without requiring a positive weight [5,7,8]. The complex nature of the method allows the circumvention of the sign problem, even when it is severe [24][25][26]. However, convergence to wrong limits has been observed both at Euclidean time [27][28][29][30][31], and real time [32,33]. These cases of incorrect convergence can be identified a posteriori, based on the theoretical justification of the method [34][35][36][37]. Further discussions on the criteria for correct convergence of com-plex Langevin can be found in Ref. [38,39]. Moreover, gauge cooling (GC) [40] has improved the convergence of complex Langevin simulations for gauge theories. The effects of gauge cooling on the complex Langevin method have been studied analytically in [41]. Investigations of gauge cooling in random matrix theories has been performed in [42,43].
Complex Langevin simulations, combined with gauge cooling, have successfully been used in QCD with a hopping expansion to all orders [44], with fully dynamical staggered fermions [45] and to map the phase diagram of QCD in the heavy-dense limit (HDQCD) [46]. In that work, we noticed that, despite the use of gauge cooling, instabilities might appear during the simulations. Here, we introduce and elaborate on our method of dynamic stabilisation (DS), which has been constructed to deal with these instabilities. This paper is organised as follows: in Sect. 2 we review the complex Langevin method. Section 3 motivates and introduces the method of DS. Tests of this procedure, applied to QCD in the limit of heavy-dense quarks (HDQCD) [47,48] are discussed in Sects. 4 and 5. Section 6 shows the outcome of applying dynamic stabilisation to simulations with staggered quarks at zero chemical potential and at finite chemical potential and high temperatures. We summarise our findings in Sect. 7. Appendices A and B review the HDQCD approximation and the staggered formulation of lattice quarks, which have been used in our investigations.

Complex Langevin
We study QCD by employing the method of stochastic quantisation [23], with which quantum expectation values can be computed using Langevin dynamics. These expectation values are evaluated as averages over a stochastic process, in which dynamical variables are evolved over a fictitious time θ . Notably, importance sampling does not enter in this formulation. The partition function for lattice QCD in the grand canonical ensemble, where the (quark) chemical potential μ couples to the quark number, is where, for the second equality, the bilinear quark fields have been integrated out. U represents the gauge links, S YM is the Yang-Mills action, S F = d 4 xψ Mψ is the fermion action with M being the fermion matrix, which depends on the gauge links and the chemical potential, and S = S YM − ln det M. We consider a SU(3) gauge theory with links U x,ν , defined on a lattice of spatial volume N 3 s and temporal extent N τ . A Langevin update, using a first-order discretisation scheme in the Langevin time θ = nε, is given by [52]: (2) where λ a are the Gell-Mann matrices, with Tr[λ a λ b ] = 2δ ab , and η a x,ν are Gaussian white noise fields satisfying The Langevin drift, K a x,ν , is obtained from the action S, where D a x,ν is the gauge group derivative The trace in the fermionic contribution to the drift is usually evaluated stochastically, see Sect. 6 for more details and alternative methods. Poles may appear in the drift in the presence of quarks, when det M = 0 and M −1 does not exist. In some situations this has a negative impact on the results [30,53,54], but, as far as understood, this is not the case in HDQCD [44,55]. For further reference, we refer the reader to the extended discussion on the issues arising from the branch cuts of the logarithm of the determinant [53][54][55][56]. In [57] it was clarified that it is the drift's behaviour around the poles, rather than the branch cuts, that affects the reliability of the complex Langevin method. It is also necessary to employ adaptive algorithms to change the Langevin step size ε, in order to avoid runaways due to large values of the drift [58].
When the sign problem is present, the Langevin drift is complex. This results in the exploration of a larger configuration space. The sign problem is circumvented by allowing the gauge links to take values in enlarged manifolds [5][6][7][8]24,34,48]. In the case of QCD, the gauge group extends from SU(3) to SL(3, C). The "distance" from the unitary manifold can be used to identify these trajectories. A possible measurement of this distance is given by the unitarity norm: It has been shown that simulations in which the unitarity norm is kept under control lead to reliable results, matching exact ones or results from different methods, when available [40,46]. One procedure to reduce the distance to the unitary manifold is known as gauge cooling [40]. It consists of a sequence of SL(3, C) gauge transformations, designed to decrease the unitarity norm in a steepest descent style: The transformation parameters, f a x , are obtained by requiring that the first variation of d with respect to a gauge transformation is negative semi-definite. The coefficient α can be changed adaptively to optimise the cooling procedure [59]. A variable number of gauge cooling steps, depending on the rate of change of the unitarity norm, can be applied between successive Langevin steps [60].
In our studies involving the heavy-dense limit of QCD (HDQCD) [47,48], we have considered the expectation value of the traced (inverse) Polyakov loops, where V is the spatial volume. The average Polyakov loop is an order parameter for Yang-Mills theories, as it is related to the free energy of a single quark by P ∼ e −F q /T . In the presence of dynamical quarks, it is no longer an order parameter. However, it still provides information on whether quarks are free or confined within hadrons. Another useful observable is the average phase of the quark determinant, measured in a phase quenched ensemble, When the sign problem is mild, the phase is not expected to vary much, leading to an average close to unity. On the other hand, in situations with severe sign problems, e 2iφ can average out to zero. When dealing with fully dynamical quarks, we have studied the chiral condensate: This is an order parameter only for massless quarks, but like the Polyakov loop, still provides information on quark confinement in general. The data was first presented in Ref. [46] 3 Dynamic stabilisation We found that even with a large number of gauge cooling steps, instabilities still may appear [46] in simulations of HDQCD, which we briefly review in Appendix A. These change the distribution of the observables during the Langevin process and lead to wrong results. Figure 1 shows the Langevin time evolution of the Polyakov loop and of the unitarity norm. This situation has a very mild sign problem, with average phase e 2iφ = 0.9978(2) − 0.0003(57)i, and thus results from reweighting are reliable. We observe two distinct regions: one is characterised by a sufficiently small unitarity norm and agreement between gauge cooling and reweighting results. At a larger Langevin time, i.e. θ 50, the agreement disappears as the unitarity norm becomes too large. It has been concluded in Ref. [46] that a large unitarity norm is an indicator of these instabilities, with 0.03 being a conservative threshold, after which results become unreliable.
To keep the unitarity norm under control we have developed a new technique-DS-which consists of adding a SU(3) gauge invariant force to the Langevin drift. This force is designed to grow rapidly with the unitarity norm d and to be directed towards the SU(3) manifold. One possible implementation is given by the substitution: with the new term 123 We remark that our choice for the additional force, which acts equally in all four directions, is not unique. The parameter α DS allows us to control the strength of the this force. A similar strategy has been used successfully for nonrelativistic fermions in one dimension [61,62]. We point out that M a x is not invariant under general SL(3, C) gauge transformations, but it is with respect to SU(3) transformations. Moreover, it is not holomorphic, since it is constructed to be a function of only the non-unitary part of the gauge links, i.e., of the combination UU † . This is necessary to make M a x scale with the unitarity norm, such that explorations of the non-unitary directions can be controlled. Therefore, it cannot be obtained from a derivative of the action. This invalidates the standard justification for the validity of complex Langevin [34,35], which require a holomorphic Langevin drift. Nevertheless, numerical evidence of the convergence to the correct limit of CL simulations with dynamic stabilisation will be shown in Sects. 5 and 6.
A naïve expansion of M a x in powers of the lattice spacing is possible if one writes the SL(3, C) gauge links as U The continuum behaviour will be discussed in Sect. 5, where we show results for different gauge couplings. By construction, the DS drift is purely imaginary and thus acts only on the imaginary parts of the Langevin drift. Checks of how the Langevin and DS drifts behave in a situation with a severe sign problem are shown in Sect. 4. An initial result of complex Langevin simulations using one step of gauge cooling 1 and DS is shown in Fig. 2, where a comparison of results from DS, gauge cooling 2 and reweighting is shown. We have used the same parameters of Fig. 1 and found agreement with reweighting for the entire length of the simulation. Figure 2 also demonstrates that it is possible to stabilise complex Langevin simulations in a way that gauge cooling alone is not able to, allowing for longer simulation times and thus smaller statistical errors.

Dependence of observables on α DS
The complexity of gauge theories makes it difficult to predict the effect of the control parameter α DS on the Langevin dynamics. However, two limiting cases can be expected: for small α DS the DS drift becomes very small, essentially not affecting the dynamics. For large values of α DS , the DS force heavily suppresses excursions into the non-unitary directions of SL(3, C), which can be interpreted as a gradual reunitarisation of the gauge links. We illustrate the effect of different α DS on complex Langevin simulations of HDQCD, see appendix A, in two cases. The first scenario corresponds to an average phase of the quark determinant close to unity, i.e., when the sign problem is mild and comparisons with reweighting are possible. In the second case, the average phase is very small, indicating a severe sign problem. Both scenarios have been simulated with inverse coupling β = 5.8 and hopping parameter κ = 0.04. Additionally, one gauge cooling step has been applied between consecutive Langevin updates.
Results for the first scenario are shown in Fig. 3. These simulations use a volume of = 10 3 ×4 and chemical potential μ = 0.7. Complex Langevin results have been extrapolated to zero step size. For sufficiently large α DS our results are compatible with both reweighting and phase quenched results. With our resolution, it is not possible to discriminate between these two cases. Figure 3 seems to indicate that one could choose an arbitrarily large α DS . However, it is necessary to keep in mind that, for these parameters, these simulations have the average phase of the fermion determinant close to unity, and a very mild sign problem. When the sign problem is severe, indicated by a highly oscillating phase of the fermion determinant, reweighting cannot be applied reliably. In principle, complex Langevin combined with gauge cooling is a viable option to simulate these regions of the phase diagram. However, as seen in Fig. 1, we have observed disagreement with reweighting, when the unitarity norm becomes too large. Figure 4 shows the Polyakov loop as function of the Langevin time for a scenario with severe sign problem, without and with DS. These simulations were carried out a volume of = 8 3 × 20 and μ = 2.45. In this case the average phase is e 2iφ = − 0.0042(35) − 0.0047 (35)i, indicating a very small overlap between the full and phase quenched models. The Polyakov loop in the simulation without DS changes to a different value for θ 20 when the unitarity norm exceeds O(0.1). We use the region before the unitarity norm rises as a reference point to test dynamic stabilisation. Due to this small sampling region, the statistical uncertainties of the gauge cooling simulations are comparatively large. The results are compatible for a wide region of α DS , as shown in Fig. 5 3 , where the results from gauge cooling are indicated by green bands. We find disagreement when α DS is outside a certain window. This can be understood as follows: for α DS very small the DS drift is too small to be effective; on the other hand, large values of the control parameter cause a heavy suppression of the exploration of the non-unitary directions. Figure 5 shows the existence of a region in α DS , which agrees with GC. Our data suggests that   4 . In other words, the region of least sensitivity to α DS seems to provide the best estimate. However, a cross-check with another method is necessary to verify that these values provide the correct result. Figure 6 shows the average unitarity norm as a function of α DS for both previously studied scenarions, with average phases close to unity and close to zero (parameters shown in the figure). It is visible that DS is able to restrict the exploration of SL(3, C) to submanifolds whose distance to SU(3) decrease with α DS . When the sign problem is severe, a plateau in the unitarity norm seems to emerge for large α DS ; while  For the remainder of this section we study histograms of the drift of Eq. (13), as they are relevant in the context of the criteria for correctness [34][35][36][37]: a heavy-tailed distribution leads to incorrect results. In Fig. 7 we show the histograms of the DS drift for four choices for the control parameter α DS . These simulations correspond to the scenario with a severe sign problem, i.e. = 8 3 × 20 and μ = 2.45. Larger drifts are less frequent than smaller ones. For large enough values of α DS the histograms become more localised distributions.
For α DS = 10 0 we observe larger values of the product α DS εM a x , due to the larger unitarity norm. We remind the reader that M a x is a function of the combination UU † (see Eq. 14), similar to the unitarity norm. As α DS increases, the unitarity norm decreases and then plateaus (seen in Fig. 6), and so does M a x . Intuitively, for very large α DS the DS drift The suppression of the imaginary part of the Langevin drift can be seen in the histogram of Fig. 8. The real part of the drift is plotted in Fig. 9 and essentially remains unchanged by DS, once α DS is sufficiently large, as evident in the inset. For too small values of α DS , the system can explore a large region of SL(3, C), which causes a different behaviour and lead to convergence to a wrong limit.  8 show histograms with compact distributions, with no skirts or "tails". This indicates an absence of boundary terms, which would spoil the proof of convergence [63]. A more in-depth analysis, however, is necessary to verify whether the lack of exponential fall-off in this situation breaks the proof of convergence.
In fig. 10 we study the influence of DS in the Langevin drift, using the same simulation parameters as in the previous plots. We show the ratio of the added DS drift over the total drift, i.e. |M|/(|K |+|M|). As the drifts are complex without a definite sign, we use the absolute values. The figure shows that even for very large values of α DS the ratio never exceeds 7%. For typical values this ratio is smaller than 1%. This is a good indication that DS accounts for a small change of the overall drift.

Continuum behaviour of dynamic stabilisation
To check the continuum limit of the DS, we have performed three simulations at different gauge couplings, specifically β = 5.4, 5.8 and 6.2, in a lattice of volume 8 3 ×20, κ = 0.04, μ = 2.45 and α DS = 10 3 . As in the previous section, we have applied one gauge cooling step between subsequent Langevin updates. The resulting histograms for the DS drift are shown in Fig. 11. As in the previous section, the histogram for the DS drift does not fall off exponentially, but this figure also indicates an absence of boundary terms. At finer lattices, larger values of the dynamic stabilisation drift M a x occur less frequently. Furthermore, the change in gauge coupling has a small effect on the real part of the Langevin drift.
The deconfinement transition of the heavy dense approximation of QCD was studied in Refs. [40,59]. The gauge coupling was varied in the interval 5.4 ≤ β ≤ 6.2. These simulations have a lattice of volume 6 3 × 6, chemical potential of μ = 0.85 and hopping parameter of κ = 0.12. It was found that the average plaquette from complex Langevin with Fig. 11 Histogram of the drift added by DS for different values of the gauge coupling β Fig. 12 The average spatial plaquette as function of β. Our simulations use a volume 6 3 × 6, μ = 0.85 and κ = 0.12. The data points have been slightly shifted for better readability just gauge cooling disagrees with reweighting for β 5.5. We have investigated whether DS can remedy this discrepancy. As in our previous studies, we added one step of gauge cooling between consecutive Langevin updates. Figure 12 shows the spatial plaquette as function of the inverse coupling. We present results generated using reweighting (RW), DS with α DS = 10 3 , and gauge cooling (GC). For comparison we add simulations with just gauge cooling with and without a cut to small unitarity norm (d < 0.03). We find good agreement between complex Langevin simulations using DS and the reweighting results from refs. [40,59,64] in both confined and deconfined phases. In Table 1 we further list the average spatial plaquette between simulations that used reweighting, dynamic stabilisation and just gauge cooling. The discrepancy between reweighting and gauge cooling results is clearly visible. On 123 Table 1 The average value for the spatial plaquette, from HDQCD simulations at 6 3 × 6, κ = 0.12 and μ = 0.85, using reweighting, dynamic stabilisation and gauge cooling. Reweighting data have been taken from [40,59,64]  In order to evaluate the fermionic contribution to the Langevin drift of Eq. (4) we employ a bilinear noise scheme and the conjugate gradient method to calculate the trace and inverse, respectively. The fermion matrix we use is that of staggered quarks, reviewed in Appendix B. One characteristic of the bilinear noise scheme is that, at μ = 0, the drift is real only on average [45]. Therefore, a non-zero unitarity norm is expected even for vanishing chemical potential. This can cause simulations to diverge, and has been addressed in Ref. [31] where the necessary elements of the inverse were calculated exactly. It has been pointed out that solutions in the confined phase are wrong, showing that stability does not imply correctness. We have investigated whether DS is able to successfully keep the unitarity norm under control, by comparing complex Langevin and hybrid Monte-Carlo (HMC) simulations 5 . We have used four different lattice volumes, 6 4 , 8 4 , 10 4 and 12 4 .
First, we have identified a suitable value for the control parameter α DS following the procedure in Sect. 4, i.e. using Eq. (18). After finding the optimal values for α DS for each lattice size, we extrapolated the results to zero Langevin step size. We have performed studies with four degenerate quark flavours of mass m = 0.025 and inverse coupling β = 5.6. We have analysed the average values of the plaquette and (unrenormalised) chiral condensate. For these parameters, the chiral symmetry is broken, as indicated by ψψ = 0. The results for the chiral condensate at zero chemical potential for = 6 4 are shown in Fig. 13, where the green band indicates the result from the HMC run. The Langevin simulations had an average step size of ∼ 4 × 10 −5 , which leads to approximately 3000-5000 independent configurations, including an auto-correlation analysis as proposed in [65]. We find very good agreement between hybrid Monte-Carlo and complex Langevin simulations. The two leftmost points have a unitarity norm larger than 0.03 and are thus not taken into account. We note here that in none of our simulations of staggered quarks the unitarity norm was zero. A behaviour similar to the points with μ = 2.45 of Fig. 6 is, in fact, observed. We used our empirical criterion of d < 0.03 [46] as threshold for reliability of the results. Figure 14 displays a comparison of the plaquette between complex Langevin and HMC for the lattice volume of 12 4 . A straight line has been fitted to the points generated by the Langevin simulations to extrapolate to zero step size, as the integration scheme is of first order [52]. We find clear agreement within the quoted uncertainties. Results for the average plaquette and chiral condensate after extrapolation to zero step size can be found in Table 2. The table shows excellent agreement between HMC and CL simulations for both observables in all four volumes considered. Recent works on complex Langevin and gauge cooling applied to staggered fermions include [66,67]. There, a discrepancy between the CLE and exact results is reported for V = 12 4 at the same inverse coupling and quark mass used here, but with two flavours of staggered fermions. This tension could be removed by using dynamic stabilization and careful extrapolation to zero step size.
6.2 Staggered quarks at μ = 0 We have carried out a qualitative simulation of staggered quarks at high temperatures spanning a wide range of μ. The chemical potentials vary from μ = 0 until saturation, where the entire lattice is filled with quarks. We use lattices with a spatial volume of V = 12 3 for two different temperatures, N τ = 2 and 4, with two degenerate quark flavours of mass m = 0.025. The inverse coupling is fixed to β = 5.6. With these input parameters, the pion and nucleon masses are m π ≈ 0.42 and m N ≈ 0.93 in lattice units [68]. At high temperatures, the inversion of the fermion matrix is numerically cheap and converges quickly even for large μ. At lower temperatures, however, we have seen that the inversion becomes more expensive, as the number of iterations easily exceed 10 4 . Further work using more state-of-the-art inverters and algorithms are under way and will enable simulations at lower temperatures. As before, we employ one step of gauge cooling and add a DS force with a control parameter of α DS = 10 3 . The chiral condensate, shown in Fig. 15, is not extrapolated to zero step size, but serves a proof of principle. Also shown in the plot are vertical lines indicating the regions of pion (μ = m π /2) and baryon condensation (μ = m N /3).

Summary and outlook
Dynamic stabilisation was introduced to deal with instabilities found in complex Langevin simulations, especially when the inverse coupling is small (typically β 5.5) or the unitarity norm rises steadily. The method is based on adding a non-holomorphic drift to the complex Langevin dynamics to keep simulations in the vicinity of the SU(3) manifold. We have studied the dependence of the observables on the control parameter α DS and have presented a criterion to tune it appropriately. We also found numerical evidence that the DS drift decreases when the lattice spacing is reduced and has a localised distribution. DS improved results on the deconfinement transition for HDQCD, previously shown in [40], where a discrepancy between reweighting and complex Langevin was observed. We find good agreement with reweighting for all gauge couplings in both confined and deconfined phases.
We presented a study of complex Langevin simulations of QCD with naïve staggered fermions at vanishing chemical potential. After extrapolating the Langevin results to zero step size, we found excellent agreement between complex Langevin and hybrid Monte-Carlo simulations for the plaquette and chiral condensate for four different lattice volumes, despite DS adding a non-holomorphic drift. Our findings rectify the discrepancy found in earlier studies in [66,67]. For μ > 0, we were able to observe changes in the chiral condensate as the chemical potential increases at high temperatures. In those cases, DS kept the unitarity norm under control and allowed for long simulations. However, the extent of these studies were limited, as simulations at lower temperatures showed a numerical difficulty arising from the inversion of the fermion matrix.
More analytical work on the justification of DS is desirable, as DS formally violates the proof of convergence of CL. Nevertheless, numerical evidence clearly shows no difference between HMC and CL, even at a sub-permille level. This needs to be confirmed at non-zero μ, by comparing with other approaches, such as those mentioned in Sect. 1.