Stabilizing complex Langevin for real-time gauge theories with an anisotropic kernel

The complex Langevin (CL) method is a promising approach to overcome the sign problem that occurs in real-time formulations of quantum field theories. Using the Schwinger-Keldysh formalism, we study SU($N_c$) gauge theories with CL. We observe that current stabilization techniques are insufficient to obtain correct results. Therefore, we revise the discretization of the CL equations on complex time contours, find a time reflection symmetric formulation and introduce a novel anisotropic kernel that enables CL simulations on discretized complex time paths. Applying it to SU(2) Yang-Mills theory in 3+1 dimensions, we obtain unprecedentedly stable results that we validate using additional observables and that can be systematically improved. For the first time, we are able to simulate non-Abelian gauge theory on time contours whose real-time extent exceeds its inverse temperature. Thus, our approach may pave the way towards an ab-initio real-time framework of QCD in and out of equilibrium with a potentially large impact on the phenomenology of heavy-ion collisions.


Contents
In recent years a number of systems have been studied using the CL method such as spin models [34,35], Bose gases at finite chemical potential [36,37] and scalar fields in realtime [38][39][40][41][42]. One of the major applications of CL has been the investigation of QCD at finite chemical potential [28,29,[43][44][45][46][47]. In contrast, apart from a few pioneering studies [39,48,49], non-Abelian gauge theories in real time have received less attention, since the sign problem associated with real-time paths has been widely believed to be particularly severe. The main problems with CL are numerical instabilities and issues with wrong convergence, i.e. the stochastic process approaches a wrong stationary solution. These problems are notoriously difficult to solve and a number of methods to reduce them have been proposed in recent years. Often numerical instabilities are related to the discretization of the Langevin equation and can be reduced by adaptive step sizes [50] or implicit solvers [41]. Problems with wrong convergence are mitigated using kernels [35,[51][52][53][54][55] and, especially in the case of gauge theories, through gauge cooling [27,56] and dynamical stabilization [57], which have led to substantial progress in the study of the QCD phase diagram.
As an important step towards a genuine real-time formulation of QCD, here we investigate lattice simulations of Yang-Mills theory on real-time contours using the CL approach. Since earlier studies have suffered from severe problems regarding wrong convergence [39], we address these issues using modern stabilization methods that are usually applied to QCD at finite density. Although applying these methods leads to remarkable improvements regarding stability on regularized complex time paths [49], they turn out to be insufficient in the approach towards the actual Schwinger-Keldysh contour. Therefore, we revisit the basic formulation of the CL equations for complex time contours and develop a new anisotropic kernel, which in numerical simulations of SU(2) Yang-Mills theory in 3+1 dimensions leads to unprecedentedly stable results with correct convergence. In particular, we find that the stabilizing effect of the kernel leads to metastable regions in Langevin time, which can be extended systematically at the expense of increased lattice sizes. The metastable regions are long-lived enough that, for the first time, we are able to obtain correct results on complex time paths whose real-time extent exceeds the inverse temperature. Our new CL equations thus mark a substantial conceptual progress that may allow us to calculate real-time observables in QCD from first principles with important applications in heavy-ion physics. This paper is organized as follows. We first revisit the formulation of CL on complex time paths for a simple quantum mechanical toy model in Sec. 2. To resolve ambiguities in the CL equation, our formulation is based on parameterizing the complex time contour. We extend our approach to lattice gauge theory in Sec. 3. Motivated by the continuum limit of the Schwinger-Keldysh contour, we exploit the kernel freedom of CL to arrive at new CL update equations using an anisotropic kernel. In Sec. 4 we demonstrate and validate that this approach leads to remarkably stable results when applied to SU(2) Yang-Mills theory in 3+1 dimensions. We conclude this work with a summary of our findings and an outlook for our future studies in Sec. 5. Additional details on CL formulations with contour parametrization, their discretized equations and useful stabilization techniques can be found in the Appendices.

Real-time complex Langevin for a simple system
In this section we discuss the CL method, its application to a quantum mechanical model on a complex time contour and resulting ambiguities from a naive formulation. We resolve these issues in Sec. 2.4 where we also employ the kernel freedom to reproduce the Minkowksi and the Euclidean cases. An unambiguous discretized CL equation is introduced in App. A.

Introducing the Complex Langevin method
The dynamics of quantum mechanical systems can be described using Feynman's path integral approach. In its real-time path integral formulation, one typically encounters oscillatory integrals for expectation values of the form for an observable O(x) where S(x) is the classical action of the theory. While such integrals can be solved analytically for simple systems, the numerical computation of expectation values for more complicated systems is typically not feasible with standard methods such as Monte Carlo integration or numerical integration methods based on discretization. The former is not applicable since exp [iS(x)] can not be interpreted as a probability density, and the latter fails due to the highly oscillatory nature of the integrand, also known as the sign problem. The complex Langevin method is an approach to overcome this issue. It introduces an artificial auxiliary time coordinate θ known as the Langevin time and uses complexified degrees of freedom, i.e., x → z(θ) ∈ C. The evolution along θ is described by a complex stochastic process known as the CL equation and is given by where dS/dz is called the drift term. The action S for complex arguments z is understood to be the analytic continuation of S(x). The real-valued noise term η(θ) is Gaussian distributed and satisfies ⟨η(θ)⟩ = 0, ⟨η(θ)η(θ ′ )⟩ = 2δ(θ − θ ′ ).

(2.3)
If the stochastic process described by the CL equation converges (see [32,33] for discussions on the assumptions), the distribution of the stochastic variable z approaches the stationary solution of the corresponding complex Fokker-Planck equation. In this case, the calculation of expectation values in Eq. (2.1) can be equivalently replaced by sampling the stochastic process at large Langevin times θ In this way, the CL method can circumvent the sign problem. However, additional ambiguities are encountered when we consider the real-time path integral formulation. In the following we discuss and resolve them for a simple quantum mechanical system before continuing with more complicated gauge field theories in the following sections.

Simple model on a complex time contour
As a concrete example, we consider a quantum mechanical model in thermal equilibrium with dependence on real time within the Schwinger-Keldysh formalism. The action is given by where x(t) is the trajectory of a particle with mass m = 1 in a potential V . The complex contour path C that is integrated over in Eq. (2.5) denotes the Schwinger-Keldysh contour [22,23] and is assumed to be continuous, starting at t = 0 and ending at t = −iβ without crossings. Here, β denotes the inverse temperature of the system. A sketch of the contour is visualized as the blue curve shown in Fig. 1a. The Schwinger-Keldysh formalism allows us to calculate expectation values via the path integral Here, x + and x − denote the trajectories on the forward and backward real-time paths C + and C − , respectively, while x E is defined on the Euclidean (purely imaginary) part of the contour C E . All times on C − are considered later than on C + and both real-time branches are separated infinitesimally along the imaginary time direction. At the start and endpoints of the contour, x(t) satisfies periodic boundary conditions In our model (2.5), the real-time paths in Eq. (2.6) lead to highly oscillatory integrals similar to Eq. (2.1).

Ambiguities for complex time contours
Our goal is here to formulate the CL equation for the model in Sec. 2.2 by complexifying the trajectory of the particle and introducing the auxiliary Langevin time θ Compared to the stochastic process described by Eq. (2.2), we have introduced an additional dependence on the physical time t ∈ C. A straightforward but naive generalization of Eq. (2.2) to the Schwinger-Keldysh time contour reads where | θ is a shorthand notation for the replacement z(t) → z(θ, t) after differentiation, i.e., denoting that the drift term is evaluated for fields at Langevin time θ. The issue with this formulation is that the noise correlator at t − t ′ for some points t, t ′ on the complex contour is ambiguous because the Dirac delta distribution is usually defined for real-valued arguments. Thus, the noise correlator must be adapted in order to obtain an unambiguous CL equation for complex time contours. However, there are at least two special cases, where this ambiguity does not arise: the Minkowski and the Euclidean time contours. The CL equation for Minkowski time t ∈ R generalizes Eq. (2.2) in a natural way with the Minkowski action On the other hand, for a Euclidean time contour t = −iτ with τ ∈ [0, β] we may write where S E denotes the Euclidean action Since there appear no complex terms in the Euclidean case, the stochastic process stays real and the CL equation reduces to the real Langevin equation. Evidently, both systems can be formulated without ambiguities, because the arguments in the noise correlator are always real-valued. Thus, our goal is to find a consistent formulation of CL which not only correctly reproduces both Minkowski and Euclidean contours as limiting cases, but is well-defined for any parametrizable complex time contour such as the Schwinger-Keldysh contour.

Parametrizing the time contour & kernel freedom
The ambiguities from the Dirac delta distribution in the noise correlator can be resolved by choosing an explicit bijective parametrization of the time contour C with λ → t(λ) on where t 0 and t 1 are the start and end points. An unambiguous CL equation, written in terms of the curve parametrization, is given by where the noise satisfies Here, the noise correlator is well-defined since the curve parameters λ, λ ′ are real-valued.
In the following we demonstrate that the parameterized CL equation correctly reproduces the Minkowksi and the Euclidean case and is invariant under reparameterization. The action in terms of an integral along the contour parameter λ reads where we have kept derivatives of the trajectories x in terms of t instead of the curve parameter λ for brevity. Computing the variation yields the drift term where t = t(λ) and t ′ = t(λ ′ ). The above relation between the correlators suggests the following transformation between η(θ, λ) and η(θ, t): The above CL equation includes additional dt/dλ factors in front of the drift and noise terms when compared to the original stochastic process defined in Eq. (2.10). However, different CL equations can be grouped into equivalence classes, where two equations are equivalent if the stationary solution in the limit θ → ∞ is described by the same probability distribution. The correspondence between Langevin equations and their stationary Fokker-Planck equation is only unique up to the so-called kernel freedom (see, e.g., chapter 4 of [58]). More precisely, all Langevin equations of the following form approach the same stationary solution of the corresponding Fokker-Planck equations: where Γ(t, t ′ ) denotes a field-independent kernel which is required to be factorizable Comparing Eq. (2.24) with the original formulation in Eq. (2.10), we notice that the additional factors correspond to a field-independent kernel where we have used Eq. (2.20) and S = −iS E due to dt = −idτ . Similarly, the noise correlator reads The parametrized Euclidean CL equation can thus be written as The above equation is equivalent to the original Euclidean process in Eq. (2.12) up to a kernel Having shown that the parameterized CL equation in Eq. (2.15) correctly reproduces both the Minkowski and the Euclidean case and, by virtue of the parameterization, has an unambiguous noise correlator for general complex time contours t(λ) ∈ C, we posit that it is the right approach to formulate the CL method for the Schwinger-Keldysh contour.
Moreover, we emphasize that one can similarly show that CL equations with different parametrizations are related by a kernel and thereby yield the same stationary solution in the case of convergence. Starting from Eq. (2.15) and introducing a reparametrization λ → ξ(λ), t(λ) → t(ξ) = t(ξ(λ)) with dξ/dλ > 0, we can perform the same steps as before to find the reparameterized CL equation In order to numerically simulate the stochastic process defined by the CL equation, we have to discretize the trajectory x(t) along the time contour. We choose a time-reversal symmetric discretization of the action Eq. (2.17) where t k = t(λ k ) ∈ C, x k = x(t k ) for k = 0, . . . , N t , and one identifies x 0 ≡ x Nt due to the periodic boundary conditions. With the time steps a t,k = t k+1 − t k , and parameterizing the time path by its arc length with the step a λ,k = λ k+1 − λ k = |a t,k |, we arrive in App. A at the discrete CL update step in the Euler-Maruyama scheme where ϵ denotes the Langevin time step and the discrete noise satisfies This is consistent with the update step used in the literature [41], where CL was applied to an anharmonic oscillator on a Schwinger-Keldysh contour.
As mentioned in Ref. [41], the Eq. (2.36) is related to a discrete CL equation without a contour parametrization by interpreting the additional factors as a rescaling of the Langevin time step This can be understood as a discretized version of the kernel freedom discussed above. Consequently, both equations will converge to the same stationary solution in the limit of a small Langevin time step ϵ → 0. However, we want to emphasize that the instabilities CL suffers from might be mitigated or even aggravated by some kerneled CL equations, as the equivalence only holds for ϵ → 0 at large Langevin times θ → ∞. In Sec. 3.5 we will exploit the kernel freedom for gauge theories in a similar way to effectively stabilize the CL equation for real-time Yang-Mills theories.

Revisiting the CL method for real-time Yang-Mills theory
In this section we formulate the CL method for non-Abelian gauge fields on a complex Schwinger-Keldysh time contour C , depicted as the blue curve in Fig. 1a. The Yang-Mills action is given by with the non-Abelian field strength tensor and totally antisymmetric structure constants f abc . The coupling constant is given by g.
Throughout this work we assume implicit summation over Lorentz indices µ, ν = 0, 1, . . . , d and color indices a, b = 1, . . . , N 2 c − 1. Employing the Schwinger-Keldysh formalism, our goal is to compute expectation values using the CL method, for gauge fields satisfying periodic boundary conditions

Complex Langevin for Yang-Mills theories on complex time contours
Applying the CL method to SU(N c ) gauge theories requires the complexification of the corresponding Lie algebra su(N c ) → sl(N c , C). The CL equation for Yang-Mills theory on complex time contours can be naively formulated as where x, x ′ ∈ R (d) denote spatial coordinates and the subscript | θ implies that the drift term is evaluated for fields at Langevin time θ. The direction and color dependent noise term η a µ (θ, x) is governed by a Gaussian distribution in each degree of freedom. We face the same issues with complex times t, t ′ ∈ C as in Sec. 2.3 because of the Dirac distribution in the noise correlator. Hence, the equation above is only consistent for Minkowski time contours where we have t ∈ R and replace S YM with the action in Minkowski time S M . The Minkowski CL equations read with the Minkowski action Similarly, there is a natural way to write down the evolution equation in Euclidean time where S E denotes the Yang-Mills action in Euclidean time Here summation over the Euclidean indices m, n = 1, . . . , d Having resolved the ambiguities of the Dirac distribution for the toy model in Sec. 2, we follow the same strategy to arrive at a consistent formulation for non-Abelian gauge theories on a general complex time contour. As before, we demand that the stationary solutions for Minkowksi and Euclidean time contours are retained. In the following section this is achieved in analogy to Sec. 2.4 by introducing a contour parametrization. However, we will additionally need to take into account that the coordinate transformation λ → t(λ) also induces a coordinate change for the gauge fields A µ .

Parameterizing the time contour
In the spirit of Sec. 2.4, we parameterize the time contour λ → t(λ) by a real-valued curve parameter λ. The CL equations for gauge fields, written in terms of the contour parameter, read where µ, ν ∈ {λ, x, y, z} and we have set the number of spatial dimensions to d = 3 for definiteness. Different from the naive formulation in Eq. (3.5), the contour parameter formulation in Eq. (3.10) involves an unambiguous Dirac delta distribution δ(λ − λ ′ ). Note that the evolution equation for µ = λ is written in terms of the λ-component A a λ (λ, x) instead of the temporal component A a t (t, x). As explained in App. B, we can write Eq. (3.10) in terms of the functional derivative with respect to fields in physical time To relate Eq. (3.11) to the aforementioned Minkowski and Euclidean cases, one needs kernel freedom. The notion of kernels discussed in Sec. 2.4 can be generalized to gauge theories. In particular, a kerneled CL equation for gauge theories is of the form

Discretizing the CL equation for Yang-Mills theory
Following the same approach as in Sec. 2, the next step is to discretize the contour parameterized CL equation for gauge fields. We proceed by approximating the underlying space-time as a discrete lattice. For the spatial part we choose a regular cubic lattice with spacing a s for each time slice. Each spatial lattice site can thus be represented by a vector x = a s 3 i=1 n iêi with unit vectorsê i parallel to the lattice axes and integer-valued coordinates n i ∈ {0, 1, . . . , N s − 1}, where N s denotes the number of spatial sites in one direction. Since we consider general complex time contours t(λ), we choose an inhomogeneous discretization for the time direction. The discrete time contour is defined by t k = t(λ k ) for k ∈ {0, 1, . . . N t } with a discrete contour parameter λ k . Here, N t is the number of temporal lattice sites along the contour. The time and parameter spacings are given by a t,k = t k+1 − t k and a λ,k = λ k+1 − λ k for k < N t . We use periodic boundary conditions for the degrees of freedom.
In contrast to our toy model, more care has to be taken regarding the degrees of freedom, namely the gauge fields A a µ (x). Lattice gauge theory provides a gauge symmetric discretization scheme, which uses group-valued gauge links as degrees of freedom. The set of gauge links {U x,µ } is given by Wilson lines along the lattice edges, i.e. connecting neighboring lattice sites x and x + a µêµ (no summation over µ), and where links with a negative index U x,−µ point into the opposite direction, i.e. U x,−µ = U −1 x−µ,µ . Thus, our goal is to recast the continuous CL evolution in Eq. (3.11) in terms of gauge links. Although our approach closely follows [39], there are a few subtleties regarding the discretized time contour, which motivate a careful re-derivation of the CL equations for gauge links. Here, we only summarize the results and refer to App. C for a detailed derivation.
Let us first introduce plaquettes U x,µν defined as 1 × 1 Wilson loops which, in the limit of small lattice spacings, approximate the field strength tensor at the mid-point of the face spanned by directionsμ andν This approximation allows us to discretize the Yang-Mills action in Eq. (3.1) in terms of plaquettes, which yields the Wilson action with the coupling constants 19) and the averaged time-step in the spatial plaquette term The latter guarantees time reversal symmetry and quadratic accuracy of the Wilson action even for general time contours. The drift term entering the CL equation is the functional derivative of the Wilson action and reads (see App. C) where t a are the traceless Hermitian generators of the SU(N c ) group andÃ a x,µ = ga µ A a x,µ . Taking all subtleties into account, we transform the CL equation in (3.11) to an update equation for gauge links on complex time contours in App. C, arriving at with the Langevin time step ϵ and the dimensionless noise correlator

Comparison to earlier approaches
Our re-derivation of the CL equations for gauge links shows that for a general time contour, additional prefactors in front of the drift and noise terms appear. It is then instructive to compare our update step to the one used in [39], which in our notation reads where the noise term satisfies the same correlator as in Eq. (3.24). Comparing this update step to Eqs. (3.22) and (3.23), we see that the main differences are contour and lattice spacing dependent prefactors. Similarly to our discussion at the end of Sec. 3.2, these can be absorbed into a CL kernel. Therefore, both schemes are equivalent due to the kernel freedom provided that a convergent solution is reached for θ → ∞. Additionally, our update scheme reduces to Eq. (3.25) if we choose a homogeneously discretized contour with a λ,k =ā λ,k = a s .
In Ref. [39] it was reported that the CL update procedure in Eq. (3.25) is prone to an instability towards wrong results, which we will refer to as wrong convergence in the following. The severity of such an instability is tightly correlated to the value of the tilt angle α of the complex time contour in Fig. 1. The smaller the angle, the earlier one finds wrong convergence, as we will discuss in Sec. 4. This makes it impossible to obtain correct results for sufficiently small tilt angles with the update step of Eq. (3.25).
To mitigate this unstable behaviour, we investigate known stabilization methods that were designed to improve convergence and stability in CL simulations of QCD at finite chemical potential. In particular, we conduct thorough tests of adaptive stepsize (AS) [50], gauge cooling (GC) [56] and dynamical stabilization (DS) [57]. Their impact on the CL dynamics for real-time Yang-Mills theory will be shown in Sec. 4.2 while details of these stabilization techniques and of how we use them are revisited in App. D. Consistent with [49], our conclusion is that while they improve the simulations, they are insufficient to obtain stable and correct results for small enough tilt angles. Therefore, another approach is required to obtain correct expectation values. In the next subsection we motivate the introduction of a novel anisotropic kernel which shows great potential by systematically improving on the stability of our simulations as shown in Sec. 4.

Regularization of the path integral and a new anisotropic kernel
We will now introduce another approach to stabilize real-time CL simulations which is well motivated by the regularization of the path integral for lattice gauge theories. It was pointed out that the path integral for the Wilson action on a Minkowski time contour has an ill-defined continuum limit and thereby does not yield a unitary theory [59][60][61]. It was proposed in Ref. [60] to resolve this issue through path integral contour deformation. The ill-defined continuum limit for the Wilson action was traced back to the character expansion of compact variables in terms of modified Bessel functions in Ref. [61].
A possible way to regularize the path integral on the Minkowski contour is to multiply the lattice-spacing dependent coupling constant by a phase factor [61] β 0 → e +iα β 0 , where α > 0. In the context of real-time simulations, this regularization can be understood as an infinitesimal Wick rotation of the positive branch C + of the Schwinger-Keldysh time contour. For the negative branch C − one has to invert the phase factor e +iα → e −iα . The resulting tilted contour is visualized as an orange line in Fig. 1a. This prescription translates to the following replacements for the temporal lattice spacing along the time contour: The regularization parameter α then corresponds to the tilt angle of the contour. In [61] it was further emphasized, that the order of limits of the temporal continuum limit a t,k → 0 and the subsequent limit to small tilt angles α → 0 has to be respected to ensure the convergence of the path integral in Eq. (3.3).
The introduced CL formalism for complex time paths via contour parametrization can be utilized to account for this regularization by simply using a tilted Schwinger-Keldysh contour with tilt angle α. Hence, we implicitly introduce the phase factors in the update steps in Eqs. (3.22) and (3.23) which are not explicitly dependent on α. An interesting aspect of our update scheme is that issues may arise when we undo the discretization of the time contour and take the continuous time limit. More specifically, taking the limit a t,k ∼ a λ,k → 0 for fixed spatial lattice spacing a s and fixed Langevin step ϵ is potentially problematic due to the blow-up of the factors a s /ā λ,k in the spatial update step.
Fortunately, we can exploit the kernel freedom of CL to circumvent this problem. The main idea is to rescale the Langevin time step for the spatial update in Eq. (3.23) via which for a λ,k < a s effectively slows the Langevin evolution of the spatial links. Similarly, we may perform an analogous rescaling for the temporal update in Eq.
This contour-dependent rescaling of ϵ can be understood as a field-independent kernel transformation applied to our contour-parameterized CL update scheme. 1 The new kerneled update steps are then given by which is well-behaved in the limit of a λ,k → 0. Comparing this kerneled update scheme to the more traditional scheme of Eq. (3.25), we see that our modification amounts to an anisotropic kernel, which treats temporal and spatial links differently. More specifically, the Langevin update for temporal links is slowed down by a factor of (a λ,k /a s ) 2 compared to the update of spatial links. Similar to the toy model in Sec. 2, we choose a λ,k = |a t,k | for practical simulations, although any other parameterization would also be admissible. In Sec. 4.4 we show that increasing the number of temporal lattice sites N t , and thereby decreasing the temporal lattice spacing a t,k , leads to improved stability of our simulations. We want to highlight that this systematic behaviour is precisely what is needed to obtain correct results for the continuum limit for the following reasons. We have already remarked that the correct order of limits is followed by calculating a t,k → 0 of the path integral first for finite α and subsequently taking the limit to vanishing tilts α → 0. In practice we can utilize the continuum limit via the anisotropic kernel to reach better stability of the simulation, which at the same time allows us to simulate for decreasing tilt angles. Our approach is therefore well-suited to sample the path integrals of SU(N c ) gauge theories on the Schwinger-Keldysh contour.

Stabilizing real-time Yang-Mills simulations
We test the impact of our kernel on the stability of the CL method by simulating SU(2) Yang-Mills theory on a 3+1 dimensional lattice. To reach improved stability we utilize modern stabilization techniques summarized in App. D in conjunction with our novel anisotropic kernel in the CL update step (3.32) and (3.33). Compared to the traditional update step of Eq. (3.25), our kernel introduces an explicit lattice spacing dependence of our CL step. This allows us to systematically improve the stability of the simulations by increasing the number of temporal lattice sites N t .
The use of stabilization methods combined with our anisotropic kernel alleviates problems of wrong convergence to a large degree, but does not fully resolve them. In particular, we are not able to achieve correct convergence for arbitrarily large Langevin times θ → ∞. Similar to what has been demonstrated in [39], the general behaviour of the CL evolution is as follows: the stochastic process first undergoes a transition away from the initial state towards a region of correct convergence. Depending on the physical and numerical parameters of the system, the evolution resides in this region for a certain time, which we term the metastable region. However, as the evolution continues, we observe that the process moves towards another stable solution that exhibits wrong convergence, i.e. the expectation values eventually converge to wrong results. This particular instability, which drives the process away from the region of correct convergence, can be strongly reduced by our new update equations. We find that the anisotropic kernel allows us to systematically extend the metastable region by increasing the resolution along the complex time contour. Hence, we are able to extract correct expectation values by taking samples only from the metastable region, which can in principle be made arbitrarily large at the cost of an increased lattice size. Additional statistics can be gained by sampling over independent simulations. However, we will not follow this strategy in this paper since here our goal is to assess the CL method with our new kernel.

Numerical setup
To compare our results with Ref. [39], we consider an isosceles contour as depicted in Fig. 1b. The forward and backward parts C + , C − of the isosceles contour connect t = 0 with t = max[Re(t)] − iβ/2 and t = max[Re(t)] − iβ/2 with t = −iβ respectively, where β is the inverse temperature 2 . The maximal real-time extent max[Re(t)] of the isosceles contour can be related to the tilt angle α via . (4.1) The Euclidean contour is realized in the limit α → π/2 where tan (α) → ∞. Large but finite values of tan (α) lead to small real-time extents, whereas for tan α ≤ 0.5, the realtime extent exceeds the inverse temperature max[Re(t)] ≥ β. Although we only present results for the isosceles contour, we note that with our new update scheme, we observed similar improvements of stability for discretized Schwinger-Keldysh contours depicted as the orange line in Fig. 1a.
The calculation of expectation values using the CL method is achieved by sampling the observable from the stochastic process of the link configuration. We initialize our simulations with homogeneous configurations of identity matrices 3 and evolve them sufficiently until the system equilibrates. Once a metastable equilibrium is reached, observables fluctuate around constant values. In practice, we use the time-translation invariance of thermal equilibrium to estimate the boundaries of the metastable region. As we will show below, we compare a time-independent observable O on the complex time path to the results from a (stable) Euclidean simulation where the value can be calculated to high accuracy. For instance, we set the beginning and the end of the metastable region at the points when the observable roughly coincides with the Euclidean simulation and when a plateau has formed around which the observable fluctuates. We have checked that for a sufficiently long metastable region, our exact choice for its starting and ending Langevin times becomes less important.
In general, observables at different Langevin times are correlated. To obtain unbiased expectation values, we sample data points that are sufficiently far apart with respect to the autocorrelation (Langevin) time. We approximate the autocorrelation function R O associated with an observable O with and find the autocorrelation time T O by a parameter fit with an exponential decay. The standard deviation of the observable at Langevin time θ is denoted by σ θ . In most figures of this section where we show the evolution of an observable over θ, we rescale the Langevin time by the autocorrelation time. This allows a fair comparison of different CL simulations. Additionally, we use simple moving averages to reduce the overall noise of our data, which allows us to better see the systematic behaviour in our plots. We consider systems with inverse temperature β = 4 in lattice units, spatial lattice spacing a s = 1 and coupling constant g = 1 on a 3+1 dimensional lattice with N s = 4 and N t = 16 if not stated otherwise. In our current setup the temporal lattice spacing a t is chosen such as to match a given N t , tilt angle and β in lattice units. For the stated standard values this implies a bare spacing of a t = a s /4 for a Euclidean contour. The stabilization methods discussed in App. D have additional numerical parameters that must be chosen appropriately. We set these parameters by hand to optimize the stability of the simulations while keeping possible biases small. We use a bound parameter B = (2π) 3 for the adaptive step size (AS), which removes runaway instabilities without distorting the results. For gauge cooling (GC), the number of gauge cooling steps is chosen between one and five with a force parameter within α GC ∈ [10 −4 , 5 · 10 −3 ]. For dynamical stabilization (DS) we observed that the parameter α DS ∈ [1, 10 3 ] leads to a stabilizing effect for sufficiently large tilt angles α, but introduces a small bias in the obtained results. We observed that we lose the stabilizing effect for smaller α DS and introduce a significant bias for larger α DS . Note that DS always introduces a bias due to its modifications to the CL equation, as will become evident from Fig. 2a and Table 1. We will therefore refrain from combining this stabilization method with our kernel and only use it for comparisons.

Effectiveness of the new CL update scheme
We assess the effect of the introduction of the anisotropic kernel in Eqs.
We compare our simulation results with those obtained from the traditional update scheme Eq. (3.25) originally used in [39]. Focusing on time-independent observables has the remarkable benefit that they can be computed from real Langevin simulations on Euclidean contours where no sign problem occurs and the gauge links remain unitary. As a result, the Euclidean simulations do not exhibit any instabilities towards wrong convergence and reliably produce reference data that allows us to check the validity of the CL results. The numerical results of our simulations are summarized in Figs. 2, 3 and Table 1. In the table, the expectation values are computed in the metastable region after a short transition time. We take samples in half-steps of the autocorrelation time T O to calculate the mean values and the statistical errors. In the case of unstable simulations, where the autocorrelation time becomes very small due to large fluctuations, we simply average over the whole ensemble of samples. The CL simulations are carried out for different tilt angles α of the isosceles contour. In Ref. [39] it was shown that the severity of the instabilities towards wrong convergence regions depends on the tilt angle, as smaller tilts lead to an earlier transition towards the wrong distribution. The dotted lines in Fig. 2a reproduce the values presented in Ref. [39]. We observe that the Langevin process initially approaches the correct value for the observable tan(α)  for some angles but then transitions and converges to another wrong expectation value.
The transition towards wrong convergence can be mitigated in part by using stabilization methods such as gauge cooling (GC) and dynamical stabilization (DS). For a large tilt angle tan(α) = 2.0 we see a promising stabilization effect when the configuration is iteratively gauge transformed via the GC procedure. However, for smaller tilts we found that GC is not sufficient and can even lead to numerical instabilities associated with the blow up of the gauge gradient. Confirming previous studies [49], a combination of GC and DS can stabilize tilt angles down to tan(α) = 1.0 but introduces a bias due to the DS penalty term, which leads to a slight offset of the solid blue line in Fig. 2a and of the corresponding expectation value given in Table 1 when compared to the Euclidean value. We further observe that DS breaks down for tilt angles as small as tan(α) = 0.5 due to the rapid growth of the unitarity norm (as we will see in the next subsection), which leads to a dominating penalty term. In this case, we found no admissible choice of the DS parameter α DS such that DS stabilizes the simulation without significantly distorting the expectation value.
The results of our simulations with the anisotropic kernel are shown in Fig. 2b and at the bottom of Table 1. In order to stabilize the simulations for small tilt angles, we increase the number of temporal lattice sites N t . In addition, we perform gauge cooling after each update step, but do not use dynamical stabilization. As can be seen from the expectation values in Table 1, our anisotropic kernel yields results in very good agreement with the Euclidean case with the same temporal resolution for all tested tilt angles and, in contrast to the GC+DS simulations, does not induce a bias.
It is interesting to point out that merely increasing the number of lattice sites N t without using our kernel does not lead to improvements regarding stability. This is shown in Fig. 3 for two tilt angles. In both cases, we observe no improvement of stability by choosing finer discretizations of the time contour. In contrast, an extended metastable region is only achieved in the figure using the anisotropic kernel of Eqs. (3.32) and (3.33) in combination with a higher temporal resolution.

Dyson-Schwinger equations and unitarity norm
In the previous section we have shown that we can successfully reproduce Euclidean data by stabilizing CL using the anisotropic kernel. In order to gain further trust in the simulation results, we check if the Dyson-Schwinger equations hold for the plaquette variables. For spatial plaquettes, this relation is given by (see Ref. [39]) In Fig. 4 we show the numerical results of the left-hand side (LHS) and right-hand side (RHS) of Eq. (4.4) as a function of Langevin time θ. 5 We show simulations without stabilization in the left panels and results with our kernel in the right panels (as in Fig. 2b). It has been remarked in Ref. [39] that the Dyson-Schwinger equations approximately hold even for the wrong convergence region. However, we find a clear difference between the regions of correct and wrong convergence regarding the RHS. We observe that for regions of wrong convergence, as shown in the left panels of Fig. 4, the RHS of the equation is governed by fluctuations around the same values as obtained for the LHS but with a much larger variance. In contrast, stabilized simulations using GC and our kernel (right panels of Fig. 4) indicate that both LHS and RHS coincide with little noise as long as we remain in the metastable region with correct convergence. Remarkably, as can be seen from unstablized simulations with tan (α) = 1.0 in the upper left panel, the stochastic process initially transitions through a short lived metastable region with small fluctuations. Since all of our simulations show the same pattern, we conclude that the fluctuations of the RHS of the Dyson-Schwinger equations can be taken as an indicator of a metastable region. It has been demonstrated in simple models [62] that there is a connection between wrong convergence and the appearance of boundary terms. 6 Due to our results we therefore expect 5 In order to obtain a simple measure to what degree the Dyson-Schwinger equations are satisfied, we average over all lattice sites and perform the expectation value for small ranges of Langevin time (∆θ ≪ T O ). 6 See also the discussion in Ref. [63] about the relation between stationary solutions of CL equations and complexified solutions to the Dyson-Schwinger equations.
no boundary terms to emerge in the Dyson-Schwinger equations in the metastable region and plan to investigate this more carefully in the future.
As an additional observable for validating our simulations, we compute the unitarity norm Since this quantity is non-negative and vanishes identically for purely unitary gauge link configurations, it can be used as a measure of "distance" from the unitary subgroup. In other systems, small values of F [U ] have been empirically shown to be associated with correct convergence [27,56]. Our numerical results for F [U ] as a function of Langevin time are shown in Fig. 5: starting from an initially unitary configuration in SU(N c ), the CL evolution drives the system away from unitarity into the non-compact directions of the complexified gauge group SL(N c , C). This leads to a growing F [U ] with increasing Langevin time. Without stabilization (left panels of Fig. 5), the unitarity norm grows quickly and plateaus at a large value after a short time. The particular Langevin time at which the plateau is reached approximately coincides with the time where the stochastic process has transitioned towards wrong convergence (compare with Fig. 2a). On the other hand, using gauge cooling and the anisotropic kernel, we observe a drastically slowed down growth of the unitarity norm (right panels of Fig. 5). Although the growth is already strongly reduced by gauge cooling alone, we find that our anisotropic kernel can reduce the increase of F [U ] over time even further. Hence, our simulation results indicate that also the unitarity norm can be used as a validating observable for the metastable region, at least on a qualitative level. For instance, the steady growth of the norm for tan α = 0.5 suggests that the metastable region will end eventually. Interestingly, for tan α = 1 the norm remains on the same level even at late Langevin times in our simulations. This could be an indication that in this case our simulations can be extended to later Langevin times without leaving the metastable region.

Systematics of the anisotropic kernel approach
The anisotropic kernel effectively increases the Langevin time spent in the metastable region with increasing number of lattice points along the time contour at fixed maximal real-time extent and inverse temperature β. To better understand this systematic behaviour, we perform multiple simulations at a fixed tilt angle, while varying the discretization of the time contour. The results are shown in the upper panels of Fig. 6 for the average spatial plaquette on an isosceles contour with tan(α) = 0.5 and tan(α) = 0.4. In the latter case, the real-time extent exceeds the inverse temperature β (see Eq. (4.1)). We observe that the width of the plateau, where the observable fluctuates around the correct value, increases relative to the autocorrelation time T O for increasing N t . Simulations using the anisotropic kernel are, however, not indefinitely stable. We can see that the stochastic process deviates from the correct values at large Langevin times. Nevertheless, it is apparent that the extent of the metastable region can be systematically enlarged in order to calculate expectation values to sufficient accuracy as shown in Table 1. We emphasize that this systematic    Figure 7: Normalized histograms (distributions) of the real and imaginary parts of the drift term D a x,µ S W ≡ δS W /δÃ a x,µ rescaled by iϵ with ϵ = 10 −4 for each degree of freedom without averaging over the lattice. The drift term was recorded in the metastable region of the simulations on an isosceles contour with tilt angle tan(α) = 0.625 using the indicated stabilization techniques and number of temporal lattice sites N t . behaviour of our kerneled CL update step is also seen in the evolution of the unitarity norm depicted in the lower panels of Fig. 6. With increasing N t its values remain small for a longer Langevin time region, which indicates prolonged stability of the CL simulations as discussed in Sec. 4.3. Consequently, our approach allows us -for the first time -to simulate non-Abelian gauge theories with larger real-time extent than the inverse temperature, albeit in a metastable manner. 7 We emphasize that here we go to large N t and thus anisotropic lattices in order to demonstrate the systematic growth of the metastable region. However, for practical purposes it will be sufficient to use discretizations with smaller N t such as N t = 256 or even smaller for tan α = 0.5 in Fig. 6. In this case sufficient statistics to compute observables can be achieved by performing multiple independent simulations in addition to averaging over the metastable region. This also reduces artifacts that stem from the anisotropic lattice spacings and simplifies a proper renormalization procedure, which we will report elsewhere.
Another important aspect of the convergence of CL towards the correct probability density function is indicated by a sufficiently fast decay of the drift D a x,µ S W ≡ δS W /δÃ a x,µ that enters the CL equation [32,33]. Intuitively, the histograms show how fast the stochastic process strays into the complex configuration space, which leads to instabilities or wrong convergence. Figure 7 shows the normalized histograms of the real and imaginary parts of the drift term for the simulation of an isosceles contour with tilt angle tan(α) = 0.625. The values are recorded during the metastable region without averaging over the lattice. We observe that our anisotropic kernel gradually contracts the histogram with growing N t . We emphasize that this behavior cannot be observed by increasing N t without the use of our kernel and the GC method alone only slightly narrows the histograms for the simulated system. Our kernel leads to localized histograms for N t = 256 without any skirts. This suggests that the criterion of correctness of CL, that is discussed in great detail in Ref. [64], may be satisfied by the introduction of our kernel for the observed system.

Conclusion
In this work, we have revisited the CL method for non-Abelian gauge theories on complex time paths, as required for simulations in physical time. We have stated a time-reversal symmetric and unambiguous formulation of CL equations both in the continuum and numerically on the lattice. We found that the introduction of an anisotropic kernel allows for systematic improvement of the stability of our simulations. In particular, our first objective was to obtain a consistent formulation of the CL equation. We have resolved the ambiguities of the traditional CL equation for a simple quantum mechanical toy model which we then generalized to Yang-Mills theory on complex time contours. The main issue is that commonly used CL equations are ambiguous on 7 Note that this approach typically requires more computational resources for decreasing tilt angles as we need finer discretizations to reach large enough plateaus of the metastable region. As the kernel is controlled by adapting the temporal lattice spacing, the memory requirements only grow linearly with Nt. Computational time grows more quickly, since the autocorrelation time also grows with Nt. Therefore, more update steps are required on a larger lattice to sample enough uncorrelated data points. complex time paths due to Dirac delta distributions of the time coordinates, which appear in the noise correlators. Therefore, we introduced a contour parameter formulation of the CL equation. Exploiting the kernel freedom, we have shown that it is parametrization independent and consistent with the known CL equations for time paths along the real and imaginary axis, respectively. For numerical simulations, the obtained CL equations were discretized in a time-symmetric manner such that they correctly approach the corresponding continuum limit for small lattice spacings. We additionally exploited the kernel freedom for Yang-Mills theory to obtain a novel CL update scheme given by Eqs. (3.32) and (3.33) that differs from the traditional equations in terms of an anisotropic kernel.
We have then demonstrated for an SU(2) gauge theory in 3+1 dimensions that our new CL equations lead to remarkable improvements regarding stability and convergence. Previously, the main problem has been that CL simulations with earlier update equations suffered from severe instabilities and converged to wrong results. In particular, these issues become significant for simulations of Yang-Mills theory on tilted time contours, which, in the limit of vanishing tilt angle, approximate the Schwinger-Keldysh time contour required for Yang-Mills theory in physical time. We have shown here that recently developed stabilization techniques such as adaptive step size, gauge cooling and dynamical stabilization mitigate these problems but are insufficient when the tilt angle of the discretized time path becomes too small. Dynamical stabilization even introduces a bias that becomes highly problematic at small tilt angles. Our new CL update scheme paves the way towards resolving these issues. Even in the case of small tilt angles, the CL evolution exhibits a metastable region of correct convergence, whose lifetime can be prolonged systematically by increasing the resolution along the complex time contour.
Using the expectation value of the average spatial plaquette as an example, we have demonstrated that we obtain correct results for any of the previously studied tilt angles by taking samples from this region. We validated our results in multiple ways: Exploiting the time invariance of the observable, we compared our results with stable Langevin simulations along the Euclidean time contour. Moreover, we have studied the unitarity norm, Dyson-Schwinger equations, and the distribution of the drift term entering the CL update step. We have found that in the metastable region, the unitarity norm remains small and that both sides of the Dyson-Schwinger equations coincide with high accuracy as opposed to regions of wrong convergence which exhibit much larger fluctuations. Our anisotropic kernel additionally narrows the distribution of the drift term, which stabilizes the CL evolution. This procedure can be extended systematically to larger metastable regions or to smaller angles. For the first time with Yang-Mills simulations, this enabled us to obtain correct results on a tilted time contour whose real-time extent exceeds its inverse temperature.
Our novel approach introduces a powerful and promising tool towards computing realtime observables and non-equilibrium dynamics of gauge theories in CL on a Schwinger-Keldysh time contour. In perspective, this may allow us to calculate transport coefficients, viscosities and spectral functions in QCD directly from first principles, which have important applications in heavy-ion collisions and beyond. Conceptually, our method is indeed able to yield results for unprecedentedly small tilt angles, potentially allowing us to perform the limit to a Schwinger-Keldysh time contour. However, this requires large lattices and high resolution in the temporal direction. We will investigate the possibilities of our framework by extending our analysis to additional observables, the assessment of boundary terms, the extraction of unequal-time correlation functions and the approach towards a continuous Schwinger-Keldysh time contour in forthcoming studies. Furthermore, the extraction of physical observables from lattice simulations will also necessitate the determination of the physical lattice spacing and the renormalized lattice anisotropy. Beyond real-time applications, the revision of the discretized CL equation in this work and the introduction of a novel anisotropic kernel may also benefit other applications of the CL method, in particular QCD at finite chemical potential. We intend to investigate this exciting direction in the future.
its variation reads where the λ-spacings are a λ,k = λ k+1 − λ k and the time steps are a t,k = t k+1 − t k . In order to define a quantity that approaches the functional derivative in the limit of infinitesimal time steps, we write the discrete variation as which is the discrete analogue of Eq. (2.18). This allows us to read off the discrete functional derivative which in the limit of a λ,k → 0 yields the correct result of Eq. (2.19). Note that the discrete functional derivative is related to the usual derivative of S latt via A natural choice for the parameter of the time path is the arc length. This amounts to choosing for all k and yields We proceed by discretizing the noise term of the CL equation along the contour. For a small spacing of the contour parameter a λ,k → 0, we may approximate the two-point correlation function in the following way: where δ kk ′ is the Kronecker delta and the discrete (with respect to λ) noise correlator reads The factor (|a t,k | + |a t,k−1 |)/2 corresponds to the averaged time step centered around t k .
We can now formulate the CL equation with a discretized contour parameter for the considered toy model as Finally, we discretize the Langevin time θ in steps of ϵ which yields the Euler-Maruyama scheme The correlator of the discrete noise is given by

B Parameterizing the time contour for gauge theories
In this Appendix we provide the details of our considerations in Sec. 3.2 for gauge theories in a contour parameter formulation. Let us recall Eq. (3.10): with µ, ν ∈ {λ, x, y, z}. Note that it involves an unambiguous Dirac delta distribution δ(λ − λ ′ ) and is written in terms of the λ-component A a λ (λ, x). The Yang-Mills action can be written using an integral over the contour parameter We view the parameterization t(λ) as a change of the time coordinate t(λ) → λ. As such, the contraction of Lorentz indices must be carried out with the appropriate metric and the components of the gauge field transform according to The functional derivative of the action can be rewritten in terms of t-derivatives. A comparison of the integrands of the variation of S YM .

(B.6)
We can therefore write Eq. (B.1) in terms of the functional derivative with respect to fields in physical time (stated as Eq. (3.11) in the main text) In analogy to Sec. 2.4, we now consider a Minkowski time contour t(λ) ∈ R to show that we can reproduce Eq. (3.6) with our contour parameter CL equation (B.1). For this we can relate the noise correlator ⟨η a µ (θ, λ, x)η b ν (θ ′ , λ ′ , x ′ )⟩ to an expression using a Dirac distribution for real-valued time arguments in an unambiguous manner. This naturally leads to the transformation behaviour of the noise field which yields the original correlator Note that even though the spatial component of the gauge field does not acquire a factor when transforming from λ to t (see Eq. (B.4)), the spatial noise field η i (θ, t, x) must be multiplied with dt/dλ to reproduce the correlator in the Minkowski case. Taking the Minkowski action S M into account, we can now put the pieces together and obtain (B.10) We observe that for parametrized Minkowski time contours, the CL equation in Eq. (B.10) represents a kerneled version of the original Minkowski formulation in Eq. (3.6). The kernel which relates these equations reads  (3.23), which correctly approach the continuum limit for small lattice spacings.

C.1 Approximating the Yang-Mills action
A general Wilson line along an arbitrary path 8 C is given by The path C is parameterized by x µ (s) with s ∈ [0, 1] with startpoint x µ (0) and endpoint x µ (1), and t a are the traceless Hermitian generators of SU(N c ). The symbol P denotes path ordering defined by Under a general gauge transformation Ω(x) ∈ SL(N c , C) the Wilson line transforms according to In the limit of small lattice spacings, gauge links can be approximated by matrix exponentials of the gauge fields evaluated at the mid-points of lattice edges where subscript x is a shorthand for the lattice site at time t k and position x. More compactly, we write where no sum over µ is implied. Links with a negative index U x,−µ point into the opposite direction, i.e. U x,−µ = U −1 x−µ,µ . Plaquettes U x,µν are defined as 1 × 1 Wilson loops which, in the limit of small lattice spacings, approximate the field strength tensor at the mid-point of the face spanned by directionsμ andν Using the above approximation, we discretize the Yang-Mills action in Eq. (3.1) in terms of plaquettes, which yields the Wilson action Here we have used the averaged time-step in the spatial plaquette term. This guarantees time reversal symmetry and quadratic accuracy of the Wilson action even for general time contours. Adopting a similar notation to [39], the Wilson action can be written more compactly as where we introduce the coupling constants (C.14) As before, the index k refers to the time slice index associated with the lattice site x.

C.2 Relating the drift terms
The drift terms appearing in the continuum CL evolution in Eq. (3.11)  where the variation of the plaquette is given by δU x,µν = δU x,µ U x+µ,ν U −1 x+ν,µ U −1 x,ν + U x,µ δU x+µ,ν U −1 x+ν,µ U −1 x,ν + U x,µ U x+µ,ν δU −1 x+ν,µ U −1 x,ν + U x,µ U x+µ,ν U −1 x+ν,µ δU −1 x,ν . (C.20) We then reorder terms in the trace and rename indices to write the variation as The sums over |i| and |j| run over positive and negative orientations of the direction, i.e. given some arbitrary expression X i In the above expression, a single link U x,µ of the set of gauge links is perturbed by a group element e iεt a close to the unit element.
which is known as wrong convergence (see e.g. [39] where this occurs for real-time lattice gauge theories). There may be two root causes for wrong convergence. Firstly, the spectrum of the corresponding Fokker-Planck (FP) operator is required to be negative semi-definite. Otherwise there is no guarantee that the stationary solution describes the desired path integral. Secondly, a wrong stationary solution may result from non-vanishing boundary terms that spoil the criterion of correctness in the derivation of the CL method [33,62,65]. Hence, the issues of wrong convergence of the stochastic process may in some cases not be primarily of numerical nature, i.e., due to the discretization of the Langevin equation, 9 but rather of a more fundamental origin [32]. For lattice gauge theories, some of these issues can be mitigated using modern stabilization techniques, which we summarize in this part of the appendix.
Adaptive step size (AS) The run away instability can be removed by the introduction of an adaptive Langevin time step ϵ [50]. We modify ϵ by introducing an upper bound B of the drift term K a x,µ = δS W /δÃ a x,µ by reducing the step size if the latter becomes too large. This is done with the substitution We set the bound parameter B small enough that it regularizes potentially large drift terms, which are responsible for runaway instabilities. Moreover, B is tuned in such a way that it is sufficiently larger than the average of the maximal drift term in the metastable region in order to avoid large biases of the simulation results.
Gauge cooling (GC) It has been empirically shown that the minimization of a functional known as the unitarity norm F [U ] can mitigate wrong convergence instabilities. The unitarity norm measures the non-unitarity of the link configuration, or more intuitively, it measures how "far away" the complexified links are from the unitary subgroup. Gauge cooling exploits the gauge freedom of the complexified system under SL(N c , C) transformations, which allows us to minimize the unitarity norm using gauge transformations. For the minimization, we adopt the gauge cooling procedure introduced in [56] and developed further in [27]. As a non-unitarity measure we use a version of the unitarity norm which differs from the original formulation by the inclusion of the square. We found that this modification turns out to be more efficient in our simulations. The minimization process is done by gauge transforming the links according to where the gauge transformation V x = exp(i α GC γ x ) (with γ x = t a γ a x ∈ su(N, C)) is determined by a gradient descent scheme. The gauge gradient is calculated via (see also Eq. (C.26)) D a x,µ f [U ] = lim where f [U ′ (δ)] is given by the functional evaluated for the modified link configuration U ′ , which is altered at x and µ such that U x,µ is replaced by U ′ x,µ = exp (iϵt a ) U x,µ . We obtain a suitable gauge transformation by computing the gauge gradient of F [U V ] and choosing V such that the gradient points towards the steepest decent. We find where we use Gauge transforming the link configuration during the CL simulation amounts to adding another term to the FP operator of the FP equation. In [66] it was shown that this additional term vanishes when applied on gauge invariant observables and does not bias the results. Gauge cooling is iteratively performed after each CL update step. The number of gauge cooling steps and the magnitude of α GC depend on the lattice size, the initial configuration and the Langevin time step of the simulation. In [27] it was further observed that gauge cooling converges slower for larger lattices, and adaptive GC steps were introduced. However, we find that these adaptive techniques do not have a significant advantage for the purposes of our real-time Yang-Mills simulations.
Dynamical stabilization (DS) As a third stabilization technique we use dynamical stabilization [57]. It introduces an additional term in the CL update equation that counteracts large imaginary drifts. The drift term is substituted by K a x,µ →K a x,µ = K a x,µ + iα DS M a x , (D.8) where the penalty term M a x acts on the imaginary part of the drift term and is proportional to the local non-unitary of the configuration. This leads to a reduction of the unitarity norm, but can not be understood as an admissible gauge transformation. Nevertheless, this method has resulted in practical advancements in the calculation of the QCD equation of state at finite density [47]. In contrast to gauge cooling, DS has (as of yet) no rigorous justification.