Stable solvers for real-time Complex Langevin

This study explores the potential of modern implicit solvers for stochastic partial differential equations in the simulation of real-time complex Langevin dynamics. Not only do these methods offer asymptotic stability, rendering the issue of runaway solution moot, but they also allow us to simulate at comparatively largeLangevin time steps, leading to lower computational cost. We compare different ways of regularizing the underlying path integral and estimate the errors introduced due to the finite Langevin time. Based on that insight, we implement benchmark (non-)thermal simulations of the quantum anharmonic oscillator on the canonical Schwinger-Keldysh contour of short real-time extent.


Motivation
The sign problem (see ref. [1] for a mini-review) remains one of the central open challenges in modern theoretical physics and hinders progress in various different subfields. It underlies the challenges encountered in the study of transport properties of the quark-gluon-plasma [2][3][4][5][6][7][8], it is the central hurdle in the exploration of the QCD phase diagram at large Baryon density [9][10][11][12][13] and impacts the study of the thermodynamics of imbalanced Fermi gases [14][15][16], to name just three. The term sign problem refers to the fact that many strongly correlated quantum systems of phenomenological relevance can only be expressed through a path integral with complex valued Feynman weight. In turn, Monte-Carlo sampling methods, successful in case that the Feynman weight is purely real, become inapplicable and systemspecific strategies must be developed. One of the most technically challenging sign problems occurs in case of quantum systems formulated in Minkowski spacetime, where the Feynman weight amounts to a pure phase.
The sign problem has been shown to be NP-hard [17], which implies that a onesize-fits-all approach is unlikely to exist. Nevertheless, many examples are known in which the sign problem has been successfully overcome or at least tamed. An active research community (for a recent review see ref. [18]) is exploring multiple strategies. Among these are variants of reweighing, extrapolation from complex parameters and the reformulation of the system of interest in new degrees of freedom unaffected by a sign problem.

arXiv:2105.02735v1 [hep-lat] 6 May 2021
The present study sets out to contribute to ongoing efforts to beat the sign problem by considering the complexification of system degrees of freedom. This research field has a long history, giving birth to two major promising currents: the Lefshets thimble approach [19] and complex Langevin [20]. In the former, one identifies manifolds in the complex plane, the so-called thimbles, on which the imaginary part of the classical action remains constant and thus ordinary Monte-Carlo sampling may commence. A residual sign problem persists as one has to average over different thimbles. Together with the computationally demanding task of locating the thimbles, these challenges constitute two areas of active research interest.
Complex Langevin on the other hand is based on the concept of stochastic quantization [21,22]. Quantum and statistical fluctuations of a system are represented by noise in an additional (d + 1) + 1 temporal dimension, which reproduces the correlation functions in (d + 1) dimensions. In practice one is required to evolve the field degrees of freedom by a stochastic partial differential equation (SDE) in an additional, so-called Langevin time, generating representations of the quantum system along the way. Expectation values of observables are estimated by taking the mean over these field configurations. Langevin stochastic quantization has proven successful in systems in which naive Monte-Carlo methods are also applicable [23]. It has furthermore been shown to correctly simulate several systems with complex weights and thus complexified field degrees of freedom [20]. Even though straight forward in principle, it has been realized early on by the community that in its standard formulation, the complex Langevin approach suffers from three major shortcomings, which have to be addressed, before the method can serve as a reliable tool in the precision study of strongly correlated quantum systems.
The three main challenges affecting complex Langevin identified in the literature are its stability, the ergodicity in the presence of non-holomorphic actions and most crucially the convergence to incorrect results (for recent insight see e.g. [24,25]. We believe that it is paramount to disentangle each of these issues, in order to be able to solve them one-by-one. Hence we focus in this study solely on the question of stability, returning to the remaining two in future work. (I.e. in order to remain in the parameter range where the complex Langevin method itself is known to converge to the correct results, we limit ourselves to a short real-time extent in this study.) The question of stability in complex Langevin is intimately connected to the wellknown phenomenon of runaway solutions. In general, such divergent behavior can arise from two sources. Either the complex Langevin method itself does not converge to a finite result, or the numerical methods used to implement the discrete Langevin time evolution introduce artifacts, which in turn give rise to unphysical divergencies. In order to make progress on understanding the former, we must disentangle numerical artifacts from methods artifacts.
Observed early on [26], runaways are now commonly treated by deploying adaptive step-size prescriptions [27] in the solution of the stochastic Langevin dynamics. One motivation for our work is the fact that even though adaptive step-size has proven to alleviate the problem of runaways in many systems in practice, it does not prevent their occurrence in principle. I.e. runaway solutions may appear even if adaptive step-size is deployed (see e.g. [28]).
In this paper, we approach the stability of complex Langevin from the point of view of the stiffness of the underlying SDEs. While no precise definition of stiffness exists, we take the pragmatic view that it refers to systems in which naive explicit time-stepping prescriptions fail to recover the correct solution. Surveying the landscape of complex Langevin implementations, we find that a majority of studies rely on the simple forward Euler discretization. Early on, improvements in the spirit of deterministic Runge-Kutta methods have been proposed [29], but to our knowledge only a single study [30] has embraced a higher-order method that takes into account the stochastic character of the Langevin evolution equation.
Our study aims at bringing to the table some of the progress made in the solution of SDEs in other fields. In particular, we propose to deploy implicit solvers, which are designed with stiff problems in mind. Besides the simple Euler-Maruyama scheme, which we use extensively in this paper, we will discuss what ingredients are needed in order to set up higher-order schemes for SDEs, compared to the case of purely deterministic equations.
Once stable solvers are available, we can proceed to investigate the stability and accuracy of the complex Langevin method itself. Our goal lies in simulating real-time physics, which in the continuum requires a form of regularization. After exploring different ways how a regularization may be incorporated in the Langevin evolution, we implement high accuracy simulations of the (0+1) dimensional anharmonic oscillator in thermal equilibrium and as a genuine initial value problem from a Gaussian density matrix.
The paper is organized in the following way: we start in section 2 with an introduction of the equations underlying the complex Langevin approach and discuss some explicit and implicit SDE solvers for their solution. In the following, we prepare the grounds for numerical simulations by introducing and discretizing the anharmonic oscillator model in section 3.1. We discuss different ways to regularize its path integral section 3.2 and will learn how to describe the errors made by a finite step size in the Langevin evolution. Armed with this insight, we carry out benchmark complex Langevin simulations of the quantum anharmonic oscillator on the canonical Schwinger-Keldysh contour at short real-times both in thermal equilibrium and in a non-equilibrium setting in section 4. We close with a summary and outlook in section 5.

Complex Langevin and SDE solvers
The task at hand is to compute quantum statistical expectation values of an observable O. Conventionally such expectation values are formulated in terms of a Feynman path integral where Z = Dφ exp[iS [φ]] denotes the partition function.
In Stochastic Quantization, we obtain the expectation values from the evolution of the system in an artificial Langevin time τ L (for an in-depth review of the approach, see Ref. [31]). The Langevin-like evolution equation for the field φ(x, τ L ) in its simplest form consists of a drift term, derived from its classical action S[φ], as well as a Gaussian noise term η(x, τ L ) The quantity x may e.g. refer to a four-vector x = (x 0 , x) with Minkowski (realtime) x 0 and the 3 spatial dimensions x. It is important to keep in mind that the physical time (x 0 ) and the fictitious Langevin time τ L are not related to each other. Note that the delta function in the correlator of the noise encompasses all dimensions of x. This prescription places an independent stochastic process at each space-time point.
Due to the complex drift term, the field degrees complexify and we can rewrite the evolution equations instead in terms of the real and imaginary part of the field as such that eq. (2) turns into two coupled but real-valued equations for the real-and imaginary part of the field degrees of freedom We have here used the standard construction in which the noise term η(c, τ L ) is real. These are the stochastic partial differential evolution equations we will solve in the subsequent sections. The first central task is to find appropriate numerical solvers to accommodate these equations.

Numerical schemes
Stochastic partial differential equations (SDE) are a central modern tool in the modeling of various phenomena in science, technology and in particular finance.
Most of the equations arising in these research fields do not lend themselves to an analytic treatment and thus require numerical solvers. To this end, the past two decades have seen vigorous research activity in the development of accurate and efficient algorithms. One major impulse towards these developments can be found in the by now classic book by Klöden and Platen [32], which not only contains a comprehensive survey of both explicit and implicit SDE solvers but also provides a pedagogic introduction into the underlying Ito and Stratonovic calculus. One central message of the book states that the series expansions, commonly used to set up deterministic discretization schemes, need to be amended by additional terms in the stochastic case due to the different scaling properties of stochastic variables. In particular, it is shown that deterministic algorithms may yield much lower convergence rates in an SDE setting than for the PDEs they were originally designed for.
The goal of this section is to introduce some of these numerical schemes and to explicitly match the complex Langevin equations to the mathematical notation used in the literature, preparing us for a straightforward implementation through standard libraries, such as the SDE solver package found in the Julia language.
Let us formulate a general stochastic differential equation for N stochastic variables φ i in Langevin time τ L , enumerated by the superscript j.
Its evolution is governed by a diffusion term conventionally denoted by a j (φ, τ L ), which may depend on all other stochastic variables, as well as the Langevin time explicitly. We incorporate N independent Wiener processes dW j which obey the standard relations They affect the dynamical degrees of freedom φ j via the mixing matrix b(φ, τ L ). This general non-constant noise coefficient matrix may have a non-trivial dependence on Langevin time and the stochastic variables.
In the concrete case of stochastic quantization, we replace the discrete parameter j with the combination of an discrete index for field degrees per spacetime point and the continuous parameter x. Hence the sum over j turns into a combined sum and integral j → j dx. Kronecker deltas remain for discrete indices, while we have to introduce Delta functions for spacetime δ jk → δ jk δ(x − x ). This leads us to the expression which can be matched to our complex Langevin eq. (2) using Note that in its standard form, the CL noise coefficients are constant in φ and τ L . As a first step, we need to take care of the drift term, which originates from a discretized action. To this end, one discretizes physical spacetime on which the field degrees live, turning continuous x into a discrete set of coordinates x m . The integral in the classical action may be approximated using a Newton-Cotes formula with the weights ω m = ω(x m ), which yields the following drift and noise terms for the continuous Langevin time SDE This expression for the diffusion term a and the noise coefficient matrix b can be used straightforwardly in numerical schemes for stochastic differential equations. In the remainder of this section, we discuss some of these schemes in more detail.
The most common implementation of the CL dynamics deploys the simple Euler-Maruyama (EM) scheme. Discretizing Langevin time in equidistant steps of size ∆τ L we introduce the discrete Wiener process increment ∆W j λ = W j λ+1 − W j λ . The subscript λ denotes at which Langevin time step the random variables are evaluated. The update step is then given by the following expression, where summation over repeated indices is implied Here we actually refer to a whole class of EM schemes, which differ by the choice of a single real-valued parameter θ. It controls the level of implicitness. For θ = 0 one recovers the fully explicit forward EM scheme, while for θ = 1 the implicit variant ensues. The choice of θ = 1/2 is special, as it refers to a semi-implicit Crank-Nicholson-like implementation of the EM scheme.
In contrast to the deterministic Euler schemes, the different variants of the EM scheme for a non-trivial noise term share a numerical accuracy of strong order O( √ ∆τ L ). I.e. in general they perform worse than in the deterministic case, which is a common ailment afflicting the direct application of deterministic schemes to SDEs. In the case of simple CL dynamics we are fortunate however, in that the noise term remains trivial and thus the simple EM scheme can be shown to be of strong order O(∆τ L ).
While the same numerical accuracy is shared among the different members of the EM family of schemes, their numerical stability varies significantly. It is well known that the forward EM scheme is at best conditionally stable, while the fully implicit scheme θ = 1 is robust against instabilities, being what is called in the literature L-stable. Similarly, it can be shown (c.f. Crank-Nicolson) that the semi-implicit scheme θ = 1/2 is also unconditionally asymptotically stable [33].
We stress that stability and accuracy are two separate qualities of a scheme, where stability only refers to the ability of the numerical solver to follow the true solution within the limitations placed by the accuracy of the scheme. Unconditional stability however also guarantees that as long as the true solution remains bounded, the scheme will not produce divergent runaway solutions. This property is what leads us to propose the deployment of (semi-)implicit solvers for complex Langevin, as it allows us to disentangle possible breakdown of the stochastic quantization prescription from a breakdown of the numerical solver.
In general, an implicit scheme is more costly than its explicit cousin at each individual update step. We need to solve a non-linear system of equations arising from the drift term a j,n (φ λ+1 ) in eq. (10), which is commonly implemented by a variant of Newton's method. As we will see, the favorable stability properties may however allow one to choose larger step sizes in Langevin time, leading to an overall reduction in computation cost.
Similarly as for deterministic differential equations, one may improve on the simple Euler schemes by developing Runge-Kutta solvers for SDEs, usually referred to as SRK schemes. They offer a straightforward way to increase the accuracy of the solver by combining approximations of the stochastic variables at intermediate steps within one update interval. The simplest of these is the Runge-Kutta Milstein scheme [32,34] of strong order O(∆τ L ) for diagonal noise Again we have indicated a whole family of schemes, whose implicitness is governed by the θ parameter. Note that these schemes differ from a naive application of the deterministic second-order Runge-Kutta (RK2) prescription through the presence of a term quadratic in the Wiener process in the second line of eq. (11). In case of simple CL eq. (2) with constant real noise coefficients, the Milstein scheme reduces to the EM scheme, reaffirming that for trivial noise the simplest algorithm already offers order 1.0 accuracy. Let us also touch on higher-order schemes. The next order one can reach is 1.5 [32], at which the SRK prescription for constant additive noise reads Here we use vector notation and we have explicitly chosen θ = 1/2 for simplicity of the presentation. In this equation, we find that a genuinely new contribution arises even for trivial noise. It consists of a combination of the drift term together with Wiener processes dW and dZ. The latter one refers to additional independent processes with the same mean and variance as the dW as well as dW j dZ j = 0. By incorporating the proper contributions arising from Ito's lemma in the series expansions underlying these Runge-Kutta schemes one may thus construct consecutive improvements to the naive EM scheme. In our study, we will draw upon the implementation of the above-mentioned schemes through the SDE module in the DifferentialEquations.jl [35,36] library provided in the Julia language. The concrete implementations of eqs. (10) to (12) differ slightly due to performance improvements outlined in the literature (see the documentation of [35]), which however has no effect on their stability and accuracy properties.
All the methods listed above can be implemented with an adaptive step-size prescription. This offers two concrete benefits. On the one hand, the stability properties of a simulation can be improved, as the step size is adapted to fulfill the Courant-Friedrichs-Lewy stability condition at each update. For stiff problems and explicit solvers, this approach is limited in practice by the step size becoming so small that the number of steps along Langevin time grows beyond available computational power. On the other hand adaptive step size also allows us to increase the step size at intermediate times to reduce the computational burden while staying within a predefined accuracy tolerance for the update step. One drawback of adaptive step size is the fact that an analytic investigation of the properties of the solver becomes more involved. We will thus deploy adaptive step size in all simulations except those where we study the finite time discretization artifacts and look at the corrections from the discretized Fokker-Planck equation.
Many different adaptive step prescriptions are deployed in the literature. One of the more sophisticated approaches implemented e.g. in the Julia library compares updates of solvers of different order and takes the difference as an error estimate [37]. The step size then is chosen to keep this error estimate below a pre-defined threshold. More simply we may monitor the size of the drift term in the Langevin equation and adjust the time step such that the change induced by the drift term remains below a certain threshold. We have found that some adaptive step-size algorithms implemented in the literature do not contain a limit on the maximum step size, which may spoil the accuracy of the outcome and required us to implement such an upper limit by hand.

On the issue of large excursions
Having reviewed different explicit and implicit prescriptions for the solution of the complex Langevin SDE, we may now explore how these methods fare in addressing the issue of stability. A common challenge that plagues complex Langevin simulations is the occurrence of large excursions. While the overwhelming majority of trajectories contributing to the final expectation value are located in a well-contained area around the origin, some paths are found to venture significantly further out into the complex plane. A simple example of this behavior can be found in the system of a single degree of freedom, evolving in the potential V (φ) = iφ 4 . On average it leads to paths that stay within around 2 dimensionless units from the origin, however excursions up to |φ| ∼ 10 sporadically occur.
These excursions can be understood by inspecting the flow field −4iφ 3 , as shown in fig. 1 along the line where φ R = φ I . As φ I φ R the flow lines tilt upward, while for φ I φ R they tilt downwards. As one moves exactly on top of the line, the field lines will keep going straight out towards infinity. This is a property of the continuum theory and not an artifact of the numerical solution.
In principle, this is not a problem, since the noise term makes sure that one never stays on this line indefinitely. However, as the size of φ increases, the size of the flow 4iφ 3 also increases significantly compared to the noise term. In turn, after the noise kicks the system away from the diverging path, it now follows a path dominated by the drift term with only a small contribution of the noise term. Such a path tends to go out to even larger values of |φ| until it eventually returns to the dominating region. A representative example is shown in fig. 1 as the green solid line.
It has been understood that such excursions constitute one of the reasons for the stability issues of numerical implementations of complex Langevin. Let us have a look at how well the true path is recovered by the simple EM schemes for different settings of implicitness for a fixed step size. In general, the explicit method is prone to overshooting the correct trajectory (red solid line), while the fully implicit method also fails to stay close to it but does so by undershooting the correct result (violet solid line). The overshooting of the explicit method easily leads to divergent behavior as the errors accumulate. For the implicit method on the other hand the simulation remains stable. Note however that its accuracy still suffers due to the deviation from the true trajectory. The semi-implicit EM scheme on the other hand combines the best of both worlds, as it offers the unconditional stability of the implicit scheme and limits the undershooting to a minimum, as can be seen in the dashed black line in fig. 1. Note that due to the large flow the actual Langevin time spend in one of these excursions is very small compared to the total length of the trajectories needed to accumulate reasonable statistical uncertainties.
We can understand the behavior seen in fig. 1 already from an inspection of the discretization prescription in the free theory. There the drift term is linear (a(φ) = iM φ) and we can rewrite the EM scheme of eq. (10) as Taking the expectation value of the field at τ L = (λ + 1)∆τ L we get For θ = 0 the fraction simplifies to |1 + i∆τ L M | > 1, which induces an increase in the magnitude of the field value for each step. On the other hand for θ = 1 one finds |1 − i∆τ L M | −1 < 1, which represents shrinkage of the magnitude. For the special case of (θ = 1 2 ) we obtain a semi-implicit scheme with actly preserves the magnitude of the expectation value of the field. We chose the above example to illustrate a key qualitative difference between solvers of varying degrees of implicitness. It should be mentioned that for the specific scenario shown here, the differences are only sizeable since we do not use an adaptive step size. A ∆τ L = 10 −4 allows for stable dynamics close to the origin where the drift is small. However at |φ| ∼ 10 we have a relatively large drift term of around ∆τ L 4φ 3 = 0.4. In practice, using adaptive step size one would reduce ∆τ L along the excursion, keeping the deviation from the exact solution small. The conclusion that explicit schemes accumulate errors according to overshooting of the true trajectory and implicit schemes according to undershooting it remains unchanged.
3 Towards stable real-time simulations of the quantum anharmonic oscillator

Formulating and discretizing the model
Our main goal in this study is to implement stable simulations of the early realtime dynamics of a strongly coupled quantum anharmonic oscillator. This system amounts to a (0+1)d field theory prototype and has been studied in the literature in detail [38,39], establishing itself as a benchmark for the success of different real-time approaches. Real-time expectation values for an observable O arising in a system that evolves from a mixed initial state ρ can be described via the Schwinger-Keldysh closed time-path formalism. In its canonical implementation, it deals with field degrees of freedom placed on a time contour, with both a forward-facing branch along the time axis (housing φ + ) and a backward branch (housing φ − ), both of which are attached at the initial time to the density matrix ρ(φ 1 , φ 2 ) Let us consider the case of thermal equilibrium with inverse temperature β = 1/T . The density matrix takes on the standard Boltzmann form ρ ∝ exp[−βH] and we may conveniently absorb the sampling over initial conditions into a path integral along the imaginary time axis, compactified to a length of β The corresponding Schwinger Keldysh contour now contains three parts: the forward and backward real-time branch, as well as the Euclidean contour all of which contribute with their own action and which we will summarize as iS The real-time action for the anharmonic oscillator explicitly reads where λ refers to the coupling constant. We give all of our results in units of m, which is the same as setting m = 1. The first step to take is to discretize the time coordinate x 0 ∈ C along the Schwinger-Keldysh contour. To this end, we introduce a contour parameter ξ ∈ R which on the forward and backward branch refers to a real-valued time, while on the Euclidean branch refers to a negative imaginary time.
The time integral in the action becomes a line integral over the contour parameter, discretized into N C steps a j . The fields, evaluated at the discrete real-time steps are denoted by φ( j a j ) = φ j . We deploy the trapezoidal rule for the action integral, which amounts to averaging over the left and right Riemann sums. Consistently we approximate the derivative by finite differences choosing forward difference at the point j and backward differences at the point j + 1, leading us to the standard expression To implement the complex Langevin equations of motion, we need to calculate the drift term i δS δφj . Using the discretization scheme above we obtain The factor 1 1 2 (|aj |+|aj−1|) , which we have included here explicitly in the drift term must then also be consistently included in the noise term b jk = 2 1 2 (|a k |+|a k−1 |) δ jk . Note however that via a rescaling of the drift term, one may drop the factors of a if one consistently drops them also from all the Kronecker deltas associate with functional derivatives. At this stage we have only discretized the physical coordinates, leaving us with a continuous Langevin time prescription to stochastically quantize the anharmonic oscillator.
The discretization of the action already introduces numerical artifacts in the solution of the continuous-time complex Langevin equation. In order to make sure that we do not misinterpret such errors as arising from the finite Langevin time discretization in an actual simulation, we take a closer look at them here.
Analogous to constructing the transfer matrix operator we can go backward from the discretized path integral to the corresponding operator expressions while leaving the real-time step size a j finite. In that case, we have to deal with the fact that the Campbell-Baker-Hausdorff formula gives non-trivial contributions when decomposing the action integral into individual exponentials. Let us define the exponentiated Hamiltonian of the system via the following matrix elements According to our eq. (18) the RHS contains the potential evaluated at neighboring values of the field φ j+1 and φ j , which requires that the potential operator acts on both the left and right state. Since only one complete set of momentum eigenstates is involved in transforming the kinetic term back to its operator form we end up Had we considered just one single potential term in eq. (18) the corresponding operator expressions would have turned out to be Using as parameters λ = 24 and m = 1, similar to what we will deploy in the actual complex Langevin simulations in the following sections, we can now study the effects of the finite real-time spacing explicitly. To this end we compute the forward correlator φ(x 0 )φ(0) using matrix mechanics in the truncated Hilbert space spanned by the 32 lowest-lying energy eigenstates of the harmonic oscillator, according to the different effective Hamilton operatorsĤ defined above.
As can be seen in fig. 2, we find a characteristic artifact introduced by the finite real-time steps. Its main manifestation is the appearance of a non-zero value of the imaginary part of the correlator at the origin. Neighboring values of the imaginary part are correspondingly also shifted away from their true values. When using the O(a) discretization of eq. (22) this effect is of the order of 12% for a = 0.08 (with respect to the maximum value the imaginary part of the correlator takes on.) The strength of the effect scales as expected linearly with a so that we obtain 6% deviation for a = 0.04 and 4% deviation for a = 0.08/3. Switching from eq. (22) to (23) changes the sign in the shift of the imaginary part, indicating that the discretization amounts to a complex phase factor. Combining the opposing phase factors in the symmetric formulation of eq. (21) cancels out the effect to a significant extent, reducing the deviation from the continuum results to 0.2% already for a = 0.08. We will make sure to keep these discretization artifacts below the percent level in the Complex Langevin simulations in the following.

Regularizing the model
To make the continuum theory well-defined on the canonical Schwinger-Keldysh real-time contour, we need to introduce an infinitesimal damping term into the otherwise purely oscillatory behavior of the Feynman path integral weight. Otherwise, the continuous Langevin-time evolution will not be able to converge to a finite result. In an analytic setting, this is conventionally achieved by introducing by hand an additional term R > 0 in the action such that the new effective action reads S = S + iR(φ, ). A simple example of such a regulator is R = 1 2 φ 2 is e.g. discussed in [40].
It is long known that such a regulator also controls the rate of convergence in a complex Langevin simulation [38,41]. Similar to an analytic computation where the correct result is obtained by setting → 0 only a posteriori, a simulation operates at finite , which, depending on its chosen value may significantly distort the computed expectation values. Only after an extrapolation over several different simulations do we recover the true solution. It goes without saying that the closer we can simulate to the correct solution, the less troubled the extrapolation procedure will be. Therefore we will attempt to deploy as small a regulator as possible.
In this section, we will briefly discuss the classic approach to introduce a tilt in the Schwinger-Keldysh contour [38] and report on our observation that the implicit scheme itself offers a regularization by construction.

Tilted Schwinger-Keldysh contour
One way to regularize our model is to introduce an imaginary tilt in the Schwinger-Keldysh contour, as deployed e.g. in ref. [38] and shown in the leftmost panel of fig. 3. Since in the thermal setting we are interested in correlators on the forward contour (the values of the mixed correlators are related via the KMS relation) one tries to keep the tilt on the forward branch small. The backward branch on the other hand may tilt downward more steeply, as long as it reaches the negative imaginary axis before −iβ. To be more concrete, the tilted Schwinger-Keldysh contour ( fig. 3 a) has two distinct parts. Part one (C 1 ) is tilted under an angle α from 0 to t max −sin(α)βi. Part two (C 2 ) is tilted such that it arrives at the imaginary axis at the point −iβ, which due to periodic boundary conditions coincides with the starting point of the first part of the contour.
The information about the shape of the real-time contour is fully contained in the choice of the, in general complex, time step a j . It denotes the distance between point j and j + 1 on the contour and appears explicitly in the complex Langevin Figure 3 Three different realizations of the thermal Schwinger-Keldysh contour for a system with temperature T = 1 β . The leftmost setting (a) corresponds to the contour adopted e.g. in ref. [38]. Our first goal is to be able to remove the tilt on the forward contour (b) and ultimately in preparation for the non-equilibrium setting to move both the forward and backward branch very close to the real-time axis (c). We find that the inherent regularization of the implicit solver allows us to realize scenarios (b) and (c) in practice.
drift term in eq. (19). A tilt of the real-time contour manifests itself as a non-zero imaginary part in a j such that We see that if a i has a negative imaginary part, the overall prefactor becomes (−ia I j ) = +i|Im(a i )|, turning the corresponding R > 0 into a positive quantity. Coming to the conclusion that such a positive term R successfully acts as a regulator however is not as straightforward as it appears at first sight. In the case of a complex Feynman weight, the field themselves becomes complexified and R exhibits both a real-and imaginary part. In the free theory ref. [40] has shown that a positive R = 1 2 φ 2 allows us to take a well defined late Langevin-time limit of the complex Langevin dynamics with a regularization of the two-point function that amounts to φ(k)φ(−k) = i/(k 2 − m 2 + i ) supporting the downward tilt of the Schwinger-Keldysh contour as an appropriate regulator.
Reducing the tilt of C 1 reduces the strength of the regularization. It is well known and we have reconfirmed in our numerical experiments that concurrently the stochastic dynamics become more and more stiff. I.e. when deploying an explicit scheme the probability to encounter runaway solutions increases significantly as the tilt is reduced. In the classic work of [38], an explicit solver was combined with a 0.01β tilt of the contour. While this tilt is small at the early times considered in previous and also this study, it already introduces a deviation from the true solution in the unequal-time correlation functions, which goes beyond the statistical errorbars of the simulation. With the goal of extending CL simulations to later real-times in the future, we will be urged to reduce the tilt even further.
As an example let us carry out a simulation using a similar setup as in [38], with a tilt of 0.01β in the anharmonic oscillator action (eq. (18)). In order to remain in the region where complex Langevin converges to the correct result, we select as maximum real-time extent x max 0 = 0.5. As a solver the general Euler-Maruyama scheme (eq. (10)) with adaptive step-size [1] is chosen. The θ value in 10 is set to θ = 1 2 corresponding to a semi-implicit scheme, which, as we have seen in section 2.2, preserves the magnitude of the expectation value throughout the simulation well. We average over a total of 500 trajectories, each of which reaches a total Langevin time of τ L m = 100. Observables are read out every δτ L m = 0.1 in Langevin time.
We plot the real-and imaginary part of the unequal time correlation function G ++ (x 0 ) = φ(0)φ(x 0 ) − φ(0) φ(x 0 ) on the forward branch vs. the contour parameter ξ in the top panel of fig. 4. While a small effect, we can already distinguish between the analytic solution on the real-time axis (black solid) and the solution on the tilted contour ( green solid ) within the precision of our simulation. The analytic solution here is obtained again using matrix mechanics in the truncated Hilbert space spanned by the 32 energy eigenstates of the harmonic oscillator.
As shown in the magnified insets, close to x 0 = 0.5 the tilt leads to a visible deviation from the true solution. I.e. such a tilt does affect the solution at early times and will become sizeable once the simulation can be extended to a phenomenologically relevant real-time extent.
In the lower panel of fig. 4 the field expectation value φ and the equal time correlation function φ 2 are plotted vs. the contour parameter along both branches of the contour. We find that they agree with the constant value predicted by the true solution. I.e. as is known, these quantities are less susceptible to the tilt as the unequal-time correlation function.

Regularization via an implicit scheme
Previous studies and the preceding subsection have shown that introducing a large enough tilt in the Schwinger-Keldysh contour, as in section 3.2.1, allows us to regularize the oscillatory behavior of the path integral. Depending on the size of the tilt it does so effectively enough for even an explicit solver to capture the ensuing complex Langevin dynamics. The price to pay is a systematic deviation of the correlation function from the result on the real-time axis which grows with the maximum extent of the forward contour. [1] All simulations based on the implicit scheme can be carried out without adaptive step-size. The numerical cost in that case will simply be higher, as an overall smaller step-size is needed to reach the same accuracy. Nevertheless, compared to simulating with an explicit scheme, we can deploy a much larger step size. The implicit scheme already works well with ∆τL = 10 −3 , while the explicit scheme requires us to go to ∆τL = 10 −5 . In this study, our goal is to explore the potential of implicit solvers for complex Langevin. We have already seen how their simplest formulation, in form of the EM scheme, avoids the occurrence of runaway solutions in section 2.2. We will return to the evolution equations and show that the formulation of the general EM scheme harbors additional terms, which play a role in the regularization of the path integral itself.
Let us focus on the simple but relevant case of the free theory here with V (φ) = 1 2 m 2 φ 2 . The update step of the general EM scheme reads where j = ∆τL |ωj | . To simplify the derivation below, we assume without loss of generality that the step size i = is constant along the contour. In the free theory we may in addition write the action in a simple matrix form as ∂S λ ∂φj = M φ λ j . Substituting this into eq. (25) yields The explicit entries in M ij are obtained via eq. (19) as In order to proceed, we bring the implicit part of the update over to the RHS and assume that is sufficiently small to expand the inverse matrix. The relevant quantitative criterion here is that the magnitude of the eigenvalues of θM are smaller than 1, i.e., the max [|λ 1 |, |λ 2 |, ...] < 1. In turn, we obtain Let us truncate the expansion at second order in and focus on the contributions to the drift term The correction to the drift term of second order in may be absorbed into an effective action for the general EM scheme The expression for S θ tells us that the difference between the EM scheme with finite θ and the fully explicit one lies in the presence of one additional term. It is proportional to the complex unit i and both depend on the time step and the implicitness parameter. Similar to the regulator term from a tilted contour in eq. (24) it is positive and thus leads to a damping of the oscillations of the path integral. eq. (31) thus constitutes a new means of regularization unavailable to explicit solvers. The above argument is further supported by numerical tests, which show that the regularization becomes weaker as the Langevin time step is reduced. Intuitively it also agrees with the behavior of the numerical solvers we discussed in the context of large excursions in section 2.2. The term proportional to 2 in eq. (30) features a minus sign, which leads to the stable undershooting of the true solution shown in fig. 1. In turn, it is this correction that will prevent the Langevin dynamics from diverging in the late Langevin time limit, realizing the role of a regularizer in the underlying path integral.

Finite Langevin time step errors
We have argued in the preceding section that implicit solvers provide a novel intrinsic regularization of the underlying complex path integral. This regularization depends on the implicitness parameter but more importantly depends on the finite Langevin time step ∆τ L . It tells us that for finite step size our system remains well defined but that when moving towards continuous Langevin time, the dynamics will become more difficult to tame, as we concurrently remove our regularization. This conundrum can be avoided if it is possible to analytically correct for the finite Langevin step size corrections in our observables. Then we may choose a small but not too small value of ∆τ L (depending on the parameters of the system) and carry out the simulation in a well-defined manner, accounting for the difference to the continuous Langevin solution a posteriori. In this section, we set out to derive such correction terms.
Our strategy is as follows: As a first step, we follow [29] and show that the effects of an implicit solver scheme at finite Langevin step size ∆τ L can be cast in the language of an effective action for the Fokker-Planck equation. In order to exploit the wellestablished methods underlying the derivation of the Fokker-Planck equation from Langevin dynamics, we restrict ourselves to a scenario with purely real Feynman weights, i.e. the imaginary time one. We continue in a second step to guess how the effective action obtained in the real case generalizes to the complex case. This heuristic step is supported by numerical evidence, which confirms that it allows us to correct numerical artifacts introduced by finite ∆τ L in practice.
Let us again focus on the simple but relevant free theory with V (φ) = 1 2 m 2 φ 2 . The update step of the general EM scheme, now for the Euclidean action reads where j = ∆τL |ωj | and the negative sign in the drift term arises from the Wick rotation into imaginary time. For θ = 0 these dynamics have been investigated in ref. [29].
To simplify the derivation below, we assume without loss of generality that the step size i = is constant along the contour. Remember that in the free theory can write the action in a simple matrix form as ∂S λ ∂φi = M φ λ . Substituting all into eq. (32) yields Let us bring the implicit part of the update over to the RHS and take is small enough to expand the first term in parentheses The above expression, up to order 5/2 , can be written in index notation, using We will now derive the corresponding Fokker-Plank equation and the effective action based on the above update prescription. The standard approach (see e.g. [42]) is to rewrite the probability distribution for φ, denoted as P[φ] at discrete Langevin time step λ + 1 in terms of its values at step λ using a delta-distribution. The argument of the delta distribution contains the Langevin update step from λ to λ + 1 and is averaged over the ensemble After expanding the delta function in powers of f j and integrating over φ one arrives at the Kramers-Moyal expansion for the discretized stochastic process, A Fokker-Planck equation may be obtained by considering terms up to the order 2 , which are encoded in the correlation functions of the update term f . To make these explicit we use the following properties of the noise η j = 0, η j η k = 2δ jk , η j η k η l = 0 and η j η k η l η m = 4 (δ jk δ lm + δ jl δ km + δ jm δ kl ), which leads to the following four expressions: , To the lowest order in we obtain the following Fokker-Planck equation Since it is the equilibrium distribution, which is of main interest to us, let us set the LHS to zero. To be more concise we will rewrite M jk = −∇ j S k and use eq. (40) to make the replacement ∇ j P = −S j P + O( ) within the curly brackets of eq. (41), consistent with the order of the approximation. The new terms at this order then give Reexpressed as a modified action we arrive at the intermediate result In the above expression, we see that the θ parameter governs the size and sign of a real-valued addition to the action. For θ > 1 2 the contribution is positive and for θ < 1 2 it is negative, distinguishing clearly between the implicit regime and the explicit regime. Note that similar to our discussion of the large excursions, the semi-implicit case of θ = 1 2 is special, as it cancels all corrections associated with the S 2 k term. The derivation outlined above cannot be translated one-to-one into the complex case. We would need to instead express the complex Langevin evolution in terms of the real-valued joint probability distribution of the real-and imaginary part of the complexified fields. Doing so, we were unable to derive a similarly closed-form as eq. (48). The structure of the correction terms obtained in the real case however invites a heuristic generalization to the complex domain using the replacement −S → iS, which leads to Let us find out whether this expression describes the dynamics of the complex Langevin simulation in practice. Similar to the discussion for the real-valued case in [29], we can attempt to counteract the effects introduced by a finite Langevin step size in the action by a redefinition of the fields. In our case the leading order change in fields according to eq. (50) amounts tõ where S j is nothing but the drift term. Note that S jj in the free theory is just a constant, so that acting with one more derivative on it makes that term vanish. Thus the redefinition of the fields changes the action to first order in such that it cancels the S 2 j contribution in eq. (50) I.e. if eq. (50) is the correct generalization then observables evaluated in terms of φ i instead of φ should show reduced deviations from the continuous Langevin time result. For the equal time two-point function we e.g. obtain the following corrected expression For later reference, we denote the correction term linear in as Σ. In order to assess the validity of the above arguments, let us simulate the harmonic oscillator on the real-time contour c) of fig. 3, i.e. on a contour without tilt in the real-time branch, up to a maximum extent of x max 0 = 0.5. Deploying an equidistant real-time spacing on the forward and backward branch of the contour and a Langevin time step of ∆τ L = 10 −2 , we compute the difference between the numerical result and the analytic solution ∆φ 2 = φ 2 CL − φ 2 QM as the colored boxes in the top panel of fig. 5 vs. the contour parameter ξ. We find characteristic features in both the realand imaginary part of this quantity. The artifacts introduced by the implicit solver at finite Langevin time show opposite sign in the imaginary part and same sign in the real-part comparing the forward and backward branch. On the Euclidean time interval, only the real-part receives significant modifications.
Interestingly, the correction term (φ −φ) 2 taken from eq. (53), when plotted as the colored triangles in fig. 3 already follows the behavior of the deviations in a qualitative fashion. At the same time, we observe that it appears to consistently over-predict those artifacts. Limiting ourselves to the corrections linear in Langevin step size , shown in gray, we find that they capture the artifacts even more accurately. This difference between the linear and quadratic terms in to us hints at the need to include higher-order corrections in the expansion of eq. (50) to arrive at a reliable correction term beyond leading order. To conclude, we find that the linear correction terms derived from a heuristic generalization of the robust result in eq. (49) capture the discrete dynamics of our complex Langevin simulation in a qualitative fashion, lending numerical support to eq. (50).
The correction terms obtained in eq. (50) contain the first and second derivative of the action S k and S kk . We can gain additional insight into what role they play Figure 5 Comparison of the finite Langevin step-size artifacts in the equal-time correlator ∆φ 2 = φ 2 CL − φ 2 QM to the estimates of that error based on eq. (53). The filled triangles denote the full estimate including the terms proportional to 2 , which provide the correct qualitative behavior but systematically overestimate ∆φ 2 . On the other hand the leading order expression Σ, proportional to , shown as filled circles captures the error even quantitatively within the statistical uncertainty. The position along the contour is parametrized by ξ, which for ξ < 1 points to real-time values and for 1 < ξ < 2 refers to imaginary times. (top) Estimation of the errors in the free theory (harmonic oscillator) using ∆τ L = 10 −2 , as well as (bottom) for the anharmonic oscillator at λ = 24 with ∆τ L = 10 −3 . In both cases, the data is based on 1000 separate trajectories each of total length τ L = 200. from the following considerations based on the translation invariance of the integrals over the fields. Let us start by stating the fact that where we have shifted the integral in the second line, such that the constant φ 0 has been moved into the exponent. While the distribution obtained from CL is different, the expectation values of the shifted system will still be the same. We can now exploit the presence of φ 0 in the integrand and the weight to derive relations between different n-point correlation functions. Using the first and second derivative with respect to φ 0 we have Using the first derivative with n = 1 and the second derivative with n = 0 we obtain the following two expressions respectively For the free case, S jj is a constant. The first term is particularly interesting as it tells us that for continuous Langevin time, the term φ j (−iS j ) should be constant and correspond to the normalization of the system. In the presence of discrete Langevin time steps, we found that it is a term proportional to φ j (−iS j ) , which describes the corrections and which are not constant as shown in fig. 5. We thus interpret the corrections in eq. (53) as counteracting in part the deviations from the correct normalization of the continuum theory. When we derived the modifications to the Fokker-Planck equation in the free theory, we were able to express them in terms of the quantity S j . We may ask whether this expression also holds in the interacting theory. To this end, we carry out simulations of the anharmonic oscillator at short real-times on the untilted Schwinger-Keldysh contour up to x max 0 = 0.5, where complex Langevin is known to converge to the correct solution (for more details see section 4.1). The same solver as for the harmonic oscillator is deployed and we choose a Langevin step-size of ∆τ L = 10 −3 . In the lower panel of fig. 5 we plot the resulting deviations from the analytic solution (filled boxes), compared to the naive application of eq. (53) to the interacting theory (filled triangles). Again we observe that the expression up to second order in slightly overestimates the artifacts but that restricting us to the linear term in allows us to capture the discretization errors within uncertainties.
With an expression at hand that allows us to correct the finite Langevin step size corrections for small values of , we are able to exploit the regularization properties of the implicit solvers in practice. What remains for each explicit system is to choose a step-size and θ parameter, keeping in mind the trade-off between regularization artifacts and numerical cost. Having too small of a step size in the implicit scheme will reduce the effect of the regulator, which in turn will lead to the appearance of large excursions. Even though these excursions do not represent a problem in principle (the approach is inherently stable) they may lead to high computational cost if a fixed accuracy goal is prescribed. Using an intermediate step sizes ∼ 10 −3 appears to give the best trade-off for the interacting systems considered in this study. The finite step size provides an effective regulator to the path integral and the finite step-size artifacts can be remedied by the correct procedure discussed above.
We emphasize that the implicit EM scheme provides enough of an intrinsic regularization that we may forego a tilting of the Schwinger-Keldysh contour all together. I.e. we gain access to the fields very close to the actual forward and backward real-time branch of the canonical Schwinger-Keldysh contour, which is particularly useful in the study of non-equilibrium field theory, in which the forward and backward correlators are not related via the KMS relation.
Armed with the insight laid out in the previous sections we are now ready to carry out stable simulations of the quantum anharmonic oscillator at short real times.

Stable CL simulations at short real-times
In this section, we present numerical results of simulating real-time complex Langevin on the canonical Schwinger-Keldysh contour with short time extent of x max 0 = 0.5 using the implicit EM scheme. We will start out with a system in thermal equilibrium which is formulated on contour c) of fig. 3. As a second example, we take a look at a system with Gaussian initial conditions, where only the forward and backward real-time branch of the contour remains and a Euclidean branch is absent.

Dynamics in thermal equilibrium
Our simulation uses the same parameters as adopted in the classic work of ref. [38], i.e. λ = 24 and m = 1. To discretize the real-time contour c) in fig. 3 for a temperature T = 1/β = 1 and real-time extent of x max 0 = 0.5, we use 16 points for the forward and backward branch each and an additional 32 points along the negative imaginary time axis. This choice of an equidistant |a| = 0.031 guarantees that the finite time spacing artifacts to the correlation functions remain at the permille level.
To regularize the path integral we deploy the general EM scheme with its implicitness parameter set to θ = 0.6. We use the adaptive step size prescription of the Julia stochastic processes library with a maximum step size of ∆τ L = 0.005. This choice provides an efficient enough regularization to avoid costly excursions, while at the same time the deviation from the continuous Langevin result remains smaller than our statistical uncertainty. For the computation of the correlation functions of interest, field configurations are collected based on 500 different trajectories. We read out observables on each of them in intervals of δτ L = 0.1 up to a total Langevin time of τ L = 100.
In the top panel of fig. 6 we show the unequal time correlation function . The corresponding continuum Langevin time solution is plotted as solid black curve, which is obtained from a matrix mechanics computation based on the truncated Hilbert space spanned by the lowest 32 energy eigenstates of the harmonic oscillator. The magnified insets confirm that our solution accurately reproduces the continuum solution on the forward and backward branch of the canonical Schwinger-Keldysch contour up to these early real-times.
In the lower panel of fig. 6 we plot the Euclidean correlator G E (x 0 = −i(ξ − 1)) = φ(0)φ(ξ) . It features a vanishing imaginary part and a real part which correctly exhibits a symmetry around x 0 = −iβ/2, corresponding to the contour parameter ξ = 1.5 here. Again the continuum solution from matrix mechanics is given as solid black line and we find excellent agreement. Let us take a look at another set of observables, which have been discussed in the literature. The top panel of fig. 7 we show the field expectation value φ (filled box and triangle) and the equal time correlation function φ 2 (filled star and cross) along the whole extent of the simulation contour parametrized by ξ. Within the statistical uncertainties of our simulation, we find full agreement with the continuoustime Langevin solution. Should one be interested in higher precision results, one will eventually find minute differences from the continuum result, similar to those shown in the lower panel of fig. 5.
We foresee that the new insight obtained in eq. (53) will help us in future studies to distinguish artifacts arising from finite Langevin-time discretization from those connected to a convergence to the wrong result. One concrete example is the equaltime correlation function, whose deviation from a constant value has previously been taken as an indication for the arrival at an unphysical solution. If the errors from a finite ∆τ L are accounted for, can the remaining deviation be unambiguously associated with wrong convergence.

Non-equilibrium dynamics
Having confirmed the efficacy of the implicit solver for the simulation of the early real-time dynamics of the anharmonic oscillator in thermal equilibrium the next step is to move to an out-of-equilibrium setting.
We follow ref. [38] and choose a Gaussian initial density matrix. Its form allows us to incorporate the information about initial conditions into a modification of the action of the system on the first and last point on the Schwinger-Keldysh contour. Note that here the contour consists only of a forward and backward real-time branch, which are not connected via periodic boundary conditions. The most general form of the Gaussian density matrix [43] leads to the following expression for the system action The five independent parameters, which specify the Gaussian initial state, represent the initial values of the field expectation value, the two-point correlation function and their derivatives The subscript c refers to the connected correlator, in which the expectation value of the field and its derivatives are subtracted. The drift term of the discretized complex Langevin dynamics (eq. (19)) is affected by the Gaussian initial density matrix only at the boundaries of the contour. Consistent with the trapezoidal rule underlying the discretization of the action integral, we choose forward derivatives at the starting point and backward derivatives when considering the endpoint of the contour. No changes are needed at intermediate contour steps. Similar to [38] we set η = 0 andφ 0 = 0, which leads to the following two explicit terms to implement at the boundary The Langevin equation remains in its standard form for the field degrees on the forward and backward branch respectively with no changes to the noise term. As in previous studies in the literature, we deploy here m = 1 and a relatively small coupling of λ = 1. This choice leaves us safely in the regime where CLE converges to the right solution. The field starts out at a finite expectation value φ 0 = φ(t = 0) = 1 at restφ 0 = 0. The spread in the values of the initial field, encoded in the correlation function is set symmetrically σ = 1 to a value of ζ = 1. Mixing terms between field and derivatives vanish via η = 0. We distribute 32 points along each of the two real-time branches to cover the maximum time extent of x max 0 = 0.5. Similar to the thermal case we deploy the EM solver with θ = 0.6 implicitness parameter using the Julia adaptive step size prescription with a maximum Langevin step size of ∆τ L = 0.005. Statistics are collected on 500 different trajectories of length τ L = 100, reading out observables on intervals δτ L = 0.1. Comparisons to matrix mechanics are also available in this scenario, however, the energy eigenfunctions of the harmonic oscillator are not well suited for truncating this particular Hilbert space. Instead, we discretize the Hamiltonian in the coordinate basis using 1024 points in the distance range x ∈ [−10, 10], the result of which will be shown as solid lines in the subsequent plots.
Our out-of-equilibrium simulation results are collected in fig. 8. Now with time translational invariance gone, we can follow the non-trivial behavior of the field expectation value φ , plotted as solid squares and triangles in the upper panel along the contour, parameterized by ξ. The initial conditions of φ 0 = 1 as well as unit variance manifest themselves in the value φ 2 (ξ = 0) = 2. Our results agree within statistical uncertainties with the analytic solution on the real-time axis.
The most interesting unequal-time correlation functions in the out-of-equilibrium scenario are the forward and backward quantities G +− = G > and G −+ = G < , which provide access to the spectral function of the system. Plotted in the lower panel of fig. 8, we find that here the statistical error remains larger than in the thermal case at the same collected statistics, indicating the presence of larger excursions in the Langevin dynamics as in the more strongly coupled thermal case. Compared to the continuous Langevin time solution from matrix mechanics, the numerical solution again shows excellent agreement.

Summary and Outlook
In this study, we have explored and showcased the potential of implicit solvers in real-time complex Langevin simulations. With the intention to disentangle the issue of numerical artifacts, such as runaway trajectories, from foundational issues, such as the convergence to wrong results, our focus in this paper remained restricted solely to early real-times.
Two central benefits of the implicit solvers were laid out in detail. On the one hand, the implicit solvers can be shown to be unconditionally asymptotically stable, preventing the occurrence of runaway trajectories, as long as the underlying complex Langevin dynamics remain finite. Using the Langevin dynamics of the free theory as a simple but relevant example, we showed in section 2 that the difference between implicit and explicit methods lies in the accumulation of errors that either undershoot or overshoot the true trajectory. While the undershoot in the implicit case also leads to a reduction in accuracy of the solution, it manages to prevent the occurrence of runaways.
In section 3 we carried out a comparison of the update prescription for the explicit and implicit EM scheme, which revealed that the effect of the latter can be captured in one additional term in an effective action. That term takes the form of a regulator +iR and depends on the implicitness parameter θ, as well as Langevin step size ∆τ L . Since R > 0, it indeed dampens the oscillations in the underlying path integral. We conclude that this additional term provides an intrinsic regularization of the path integral unavailable to the explicit solvers.
Subsequently, we analyzed the finite Langevin time discretization artifacts in terms of an effective action in the Fokker-Planck equation for the case of a purely real path integral. We then heuristically generalized the result to the complex case and provided numerical support that our educated guess indeed captures the numerical artifacts introduced due to finite Langevin time steps in the free theory and even the strongly coupled interacting case. This correction formula allows us to exploit the inherent regularization properties of the implicit solvers in practice, as we may now simulate the system at a small but finite Langevin step size ∆τ L in a well-defined manner and correct for the effect of the regulator a posteriori.
The first three sections have provided us with insight into the regularization properties of different numerical schemes, insight into the effects of finite real-time discretization and we have derived the form of finite Langevin time steps artifacts that allow us to compensate for the effect of the regulator. We thus proceeded in section 4 to carry out benchmark numerical simulations of the anharmonic oscillator in (0 + 1)d on the canonical Schwinger-Keldysh contour without tilt and maximum real-time extent of x max 0 = 0.5. Both in the thermal case and in a scenario with Gaussian non-thermal initial conditions, we find excellent agreement between the complex Langevin simulation and the analytic solution from matrix mechanics. The fact that the implicit solver gives access to the backward path on the real-time axis allows us for the first time to compute the actual forward and backward correlators G +− = G > and G −+ = G < , whose difference encodes the phenomenologically relevant spectral function of the system.
We believe that the availability of implicit solvers and an improved understanding of discretization artifacts will help to improve the reliability of the complex Langevin approach and provides new momentum to attack the pressing open challenges associated with it. The stability and regularization properties of the implicit schemes offer benefits in other applications of complex Langevin beyond real-time simulations, such as the treatment of strongly interacting systems at finite chemical potential (for a recent review on CL and the QCD phase diagram see e.g. [13]).
Many different paths forward exist. One aspect we are following up on is the role of regularization in the path integral for the success of complex Langevin convergence. When we introduce a tilt in the Schwinger-Keldysh contour it led us in eq. (24) to a regulator term that incorporates all terms of the action. We may instead ask how the system reacts to introducing a regulator on individual terms in the action by modifying in either the kinetic, mass or self-interaction term the lattice spacing from a j → a j − iκ with κ > 0. In the following we will thus work with the contour b) of fig. 3, where the forward branch is located on the real-time axis and only the backward branch tilts downwards to intersect with the imaginary time axis at β. We have seen that for x max 0 = 0.5 the complex Langevin approach in the strongly coupled thermal scenario with λ = 24 converges to the correct solution given by matrix mechanics. Extending the contour to later real-times, we encounter significant deviations already at x max 0 = 0.8. A prominent characteristic of the incorrect solution is an artificial downward shift in the real-part of the unequal-time correlation function, as shown by the red data points in fig. 9. In addition, the curvature of the imaginary part of the correlator, given as open squares also deviates from the true solution beyond statistical uncertainty.  fig. 3). The direct simulation based on the explicit EM scheme converges to an incorrect result given by the red data points. No improvement is observed for regularizing the φ 4 term (blue). The correct solution is recovered when regularizing the momentum term (brown).
In the regime 0.75 < x max 0 1 we observe that this incorrect convergence can be overcome by a choice of regularization on the forward branch. Interestingly, when introducing an imaginary part in the a j 's associated with the interaction term (blue filled square and open circle) the simulation outcome remains unchanged. On the other hand, modifying the kinetic term with a small imaginary part a j − i × 10 −3 leads to a significant improvement as indicated by the brown-filled circle and open triangles, which agree with the analytic solution from matrix mechanics. On the other hand, such a regularization based strategy fails to achieve its purpose, once the real-time extent of the Schwinger-Keldysh contour goes beyond unity in units of the mass. Our goal in future work is to gain a systematic understanding of how the regularization achieves to recover the correct results, possibly by studying the associated Fokker-Planck equation in low-dimensional models.
Furthermore, the availability of implicit and in particular higher-order solvers benefits the systematic exploration of kernels for the Langevin dynamics (for a modern perspective on CL kernels see e.g. [44]). In the real-valued case, kernels can be used to improve the convergence properties of the stochastic quantization procedure. In complex Langevin, they have been studied with mixed success as means to remedy the convergence to wrong solutions. Robust numerical SDE solvers (c.f. the Runge-Kutta Milstein scheme of eq. (11)), which can accommodate nontrivial kernels with Langevin-time and field dependencies will allow us to explore a much broader class of kernels than before in future studies.