{'fraction_non_alphanumeric': 0.03809494432277368, 'fraction_numerical': 0.010611061414602548, 'mean_word_length': 4.188860435339309, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Complex Langevin method: When can it be trusted?

We analyze to what extent the complex Langevin method, which is in principle capable of solving the so-called sign problems, can be considered as reliable. We give a formal derivation of the correctness and then point out various mathematical loopholes. The detailed study of some simple examples leads to practical suggestions about the application of the method.


I. INTRODUCTION
The complex Langevin method solves in principle the sign problems arising in simulations of various systems, in particular in QCD with finite chemical potential. After being proposed in the early 1980s by Klauder [1] and Parisi [2], it enjoyed a certain limited popularity [3,4], but very quickly certain problems were found. The first one was instability of the simulations with absence of convergence (runaways), the second one convergence to a wrong limit [5][6][7][8]. Nevertheless in recent years the method has been revived with sometimes impressive success [9][10][11][12][13][14][15][16]. In particular the use of adaptive stepsize has eliminated the problem of runaways [17].
But nagging problems remained due to the lack of clear criteria to decide when an apparently convergent simulation actually represented the truth. This was linked to the lack of a clear mathematical basis for the method, that would at the same time also provide criteria for its applicability. The purpose of the present paper is to clarify the situation at least to some extent. While we are still not able to close certain mathematical gaps and reach a complete analytic solution to the problems that have plagued the method, we give some strong numerical evidence that the method is correct in some cases and also suggest a plausible explanation for the failure in other cases; this leads to some pragmatic conclusions suggesting how to proceed in practice in a way that promises credible results.
The paper is organized as follows. In Sec. II we give a formal justification of the method, highlighting the assumptions underlying the derivation. In Sec. III three main questions raised by the formal arguments are listed. We then focus on one particular issue, boundary effects, in Sec. IV, and present detailed case studies in Sec. V. Tentative conclusions are given in Sec. VI.

II. FORMAL DERIVATIONS
For simplicity we concentrate here on models in which the fields take values in flat manifolds M = R n or M = T n , where T n is the n dimensional torus (S 1 ) n with coordinates (x 1 , . . . , x n ). The complications that arise when the fields live in nontrivial manifolds, as is of course the case in QCD, have been successfully dealt with in the literature (see for instance Ref. [18] for real, Refs. [10,12,19,20] for complex Langevin dynamics). But these complications are not really relevant for our discussion.
As is well known, the idea is to simulate a complex measure exp(−S)dx, with S a holomorphic function on a real manifold M, by setting up a stochastic process on the complexification M c of M, such that the expectation values of entire holomorphic observables O in this stochastic process converge to the ones with respect to the complex measure exp(−S)dx.

The complex Langevin equation (CLE) on M c is
where dw denotes the increment of the Wiener process and the equation is to be interpreted as a real stochastic process, namely dx =K x dt + dw, dy =K y dt, with K x = − Re∇ x S(x + iy), K y = − Im∇ x S(x + iy).
A slight generalization of Eq. (2) that has been considered and will play a role in this investigation is where dw R and dw I are independent Wiener processes, N I ≥ 0 and N R = N I + 1. This is usually referred to as complex noise. The introduction of a nonzero N I makes it possible to solve the Fokker-Planck equation (see below) numerically and also allows a random walk discretization of the complex Langevin process: where P ± are the transition probabilities and we have defined the steps ω such as to have the same δt in both sub-processes, to ensure correct evolution.
By Itô calculus, if f is a twice differentiable function on M c and is a solution of the complex Langevin equation (4), we have where L is the Langevin operator and f denotes the noise average of f corresponding to the stochastic process described by Eq. (4). In the standard way Eq. (4) leads to its dual Fokker-Planck equation (FPE) for the evolution of the probability density P (x, y; t), with L T is the formal adjoint (transpose) of L with respect to the bilinear (not hermitian!) pairing i.e., Lf, P = f, L T P . (14) Note that the FPE has the form of a continuity equation where (J x , J y ) is the probability current in the 2n dimensional space M c , given by We will also consider the evolution of a complex density ρ(x) on M under the following complex FPE where now the complex Fokker-Planck operator L T 0 is A slight generalization will be useful: For any y 0 ∈ M we consider a complex Fokker-Planck operator L T y0 given by L T y0 is the formal adjoint of The operators L T y0 act on suitable complex valued distributions (measures) on M, parameterized by the real variables (x 1 , · · · , x n ). But they do not allow a probabilistic interpretation, because they do not preserve positivity.
For any y 0 Eq. (17) with L T 0 replaced by L T y0 has the complex density as its (hopefully unique) stationary solution.
We next consider expectation values. Let O be an entire holomorphic observable with at most exponential growth; then we set and We would like to show that provided the initial conditions agree, which is true if we choose (for any y 0 ). In the limit t → ∞ the dependence on the initial condition should of course disappear by ergodicity.
The goal is to establish a connection between the 'expectation values' with respect to ρ and P for the class of observables chosen (entire holomorphic with at most exponential growth). The idea is to move the time evolution from the densities to the observables and make use of the Cauchy-Riemann (CR) equations. Formally (i.e. without worrying about boundary terms and existence questions) this works as follows: first we use the fact that we want to apply the complex FP operators L y0 only to functions that have analytic continuations to all of M c . On those analytic continuations we may act with the Langevin operator whose action on holomorphic functions agrees with that of L, since on such functions ∇ y = i∇ x and ∆ x = −∆ y so that the difference L −L vanishes.
The proliferation of Langevin/Fokker-Planck operators may be somewhat bewildering, but it is important to realize that L,L, L y0 are really all different operators: while L andL act on functions on M c (i.e. functions of (x 1 , y 1 , . . . , x n , y n )), agreeing on holomorphic functions, but disagreeing on general functions, L y0 acts on functions on M, i.e. functions of (x 1 , . . . , x n ).
We now useL to evolve the observables by the equation with the initial condition O(z; 0) = O(z), which is formally solved by In Eqs. (28, 29), because of the CR equations, the tilde may be dropped, and we will do so now. So we will have with its formal solution In fact Eq. (28) is also equivalent to the family of equations The first thing to notice is that O(z; t) will be holomorphic if O(z; 0) is. This can be seen as follows: Let∂ j (i = j, . . . , n) be the Cauchy-Riemann operators defined by∂ Applying this to both sides of Eq. (28) and using the fact that by the holomorphy of S∂ j commutes with L, we find this is just Eq. (28) again with O(z; t) replaced bȳ ∂ j O(z; t). Under the assumption that L generates a semigroup acting on∂ j O(z; t), Eq. (34) has a unique solution; since the initial condition is∂ j O(z; 0) = 0 we conclude that∂ j O(z, t) = 0 for all t ≥ 0, i.e. O(z; t) satisfies the CR equations for all t ≥ 0 and all j = 1, . . . , n. So O(z; t) is holomorphic in each component z j separately. By Hartogs' theorem [21,22] this implies joint holomorphy.
We now consider, for 0 ≤ τ ≤ t, and claim that it interpolates between the ρ and the P expectations: The first equality is obvious, while the second one can be seen as follows, using Eqs. (26,31), where we only had to assume that we can integrate by parts in x without worrying about boundary terms.
Our desired result follows if we can show that F (t, τ ) is independent of τ . To see this, we differentiate Integration by parts then shows that the two terms cancel, hence ∂ ∂τ F (t, τ ) = 0 and thus It is important to notice that this holds for all N I ; whereas the left hand side seems to depend on N I , the right hand side is manifestly independent of it.
If we knew in addition that with ρ(∞) given by Eq. (21) with y 0 = 0, we could now conclude that the expectation values of the Langevin process relax to the desired values; this convergence would follow if we knew that the spectrum of L T y0 lies in a half plane Re z ≤ 0 and 0 is a nondegenerate eigenvalue. But note that we do not really need convergence of P (x, y; t) for Eq. (40) to hold, since it will only be tested against analytic observables O.
Nevertheless the numerical evidence in many cases points to the existence of a unique stationary probability density P (x, y; ∞). The corresponding probability currents are divergenceless, but unlike the situation in the real Langevin process, they cannot vanish. A general feature of the stationary distribution that can be read off the FPE is the following: Assume that (x 0 , y 0 ) is a local stationary point of P , then , a local maximum of P can only occur where the divergence of the drift force is negative and a minimum where it is positive. For N I = 0 the conclusion is even stronger: where the divergence is negative (positive), there cannot be a local minimum (maximum) in x for fixed y. These properties provide some checks on numerical solutions.

III. QUESTIONS
There are three main questions raised by the formal arguments in the previous section: (1) Can the operators L,L, L y0 and their transposes be exponentiated; in more mathematical language: do these operators generate semigroups on some suitable space of functions?
(2) Are the various integrations by parts justified, which underlie the shifting of the time evolution from the measure to the observables and back, or are there boundary terms to worry about?
(3) Are the spectra of L, L y0 and their transposes entirely in the left half plane and is 0 a nondegenerate eigenvalue?
Concerning the first question, there are treatises (see for instance Refs. [23,24]) giving rather general sufficient conditions for the existence of a semigroup generated by differential operators of the general type considered here. Unfortunately it seems that the cases we have to deal with here are not covered by those general results; the main difficulties are (1) the strong growth of the drift given by the gradients of the action in some complex directions and (2) the fact that the drift is not always 'restoring'.
Question (1) for L, L T is intimately related to the question whether the stochastic process given by the complex Langevin equation exists for arbitrary long times. This is not obvious because typically in the classical (no noise) limit there are trajectories that go to infinity in finite time. While those trajectories occur only for a subset of measure zero of initial conditions, it is not obvious what happens after adding the noise. On the one hand, the noise will typically kick the process away from the unstable trajectories; on the other hand it may also kick it near the unstable trajectory, inducing very large excursions of the process. We will later illustrate this with some examples.
But let us say that the accumulated numerical evidence points not only to the existence of the process for arbitrarily large times, but also to the existence of a unique equilibrium measure for the process; unfortunately we could neither find results in the mathematical literature that would imply this, nor could we prove ourselves that this is the case.
Concerning the exponentiation ofL, L y0 , which should be easier, we still could not establish mathematically that it is possible on a space of functions containing the most obvious observables, such as exponentials.
There is a useful criterion for the existence of a bounded semigroup generated by an operator A on a Hilbert space [23]: A generates a bounded semigroup if it is dissipative, i.e. if A + A † ≤ 0. Unfortunately even in the simplest cases of a quadratic S the corresponding Langevin and FP operators are not dissipative. So they can at best generate exponentially bounded semigroups; if in addition the spectrum is in the left half plane, convergence to the equilibrium should still take place.
For the second question it would be necessary to have good control over the falloff of the solutions of the FPE in the imaginary directions: if we insert the observables O k (z) ≡ exp(ikz) into Eq. (39), we get for the Fourier transform (Fourier coefficients in the compact case) ρ(k; t) of the complex density ρ(x; t) ρ(k; t) = 1 (2π) n P (x, y; t)e ikx−ky dxdy.
This makes sense for all k only if P (x, y; t) decays more strongly than any exponential in imaginary direction. Our case studies described in the following sections indicate that this does not seem to be the case: in our first example the decay is probably exponential, but not stronger; in our second example the decay seems to be even weaker (O(|y| −r ) with r ≈ 2), so that exponentials cannot be used as observables. The remainder of the paper is mainly devoted to studying Question (2) in some toy models.
Finally let us remark that the third question is more difficult than the first one, and again the answer is not known rigorously. But again the numerical evidence strongly suggests a positive answer, depending on the model and the parameter values, in many interesting and relevant cases (see e.g. Refs. [12,14]).

IV. BOUNDARY EFFECTS
In this paper we are mainly concerned with Question (2). Even though we have very little analytic control, careful numerical studies reveal that, as remarked, the answer is generally no! In our case studies we find indications that the probability density P (x, y; t) indeed relaxes to an equilibrium density P (x, y; ∞), but that that limiting density decays at best like an exponential for |y| → ∞, in other cases only power-like. This limits the class of observables for which the integrations by parts can be performed without boundary terms.
But let us first take a closer look at the integrations by parts that occur in the formal arguments of Sec. II. The danger lies in a possibly insufficient falloff in the imaginary (y) directions, whereas in the real (x) direction we have either compactness or sufficient falloff due to the behaviour of the action S.
Let us remark that for N I > 0 the operators L and L T are uniformly strictly elliptic; this is important, because it implies regularity for any solution of the stationary FPE [25]. It is also to be expected that the semigroup exp(tL T ) has a smooth (even real analytic) kernel so that for t > 0 P (x, y; t) will be smooth (real analytic). This is supported by our numerical studies. So the problematic point is to show that the two terms on right-hand side of Eq. (38) actually cancel. This may fail because the observables O typically grow in imaginary direction, whereas the decay of P (x, y; t) (always assuming it exists) may be insufficient to compensate for it.
Let us see in a little more detail how the argument for the independence of N I may fail. For simplicity of presentation we consider the one-dimensional case (n = 1). We write the Langevin operator (10) as where L NI =0 is the Langevin operator for N I = 0, N R = 1. Then let us consider We use the formula and integration by parts to rewrite Eq. (44) as The term denoted by X collects possible boundary terms arising in the integration by parts. It vanishes only if the decay of P (x, y; t−τ ) is strong enough to offset any possible growth of O(x, y; τ ); otherwise it may either converge to a finite nonzero value or diverge. By the CR equations the first term vanishes, but the uncontrolled boundary term X remains.
Let us look at a simple example that shows how and when the formal argument fails. Consider (where the second term vanishes on account of the CR equations). By the formal argument this would be zero, being just a boundary term, but careful application of integration by parts, at first over the finite domain −Y − < y < Y + , gives Here it is clear that what matters is the combined asymptotic behaviour of P and O, and depending on the observable, X may be zero, finite, or divergent. Of course the form of the boundary terms is less simple when using integration by parts in Eqs. (38,44), but we expect that it is still the decay of the products like P O, O∇ y P , P ∇ y O that is relevant.
When trying to investigate the effects of the boundary numerically, one would in principle like to use the probability density obtained without a cutoff. In practice this is, however, not feasible, and we therefore introduce a cutoff in the imaginary direction, which is not sent to infinity. Such a device is necessary for the solution of the FPE, and even though the CLE does not require it, for the purpose of comparison we also introduce it there. This will, however, introduce additional problems with the formal arguments relying on the CR equations as well as integration by parts.
Concretely we proceed as follows: We restrict each y j to lie between −Y − and Y + and impose periodic boundary conditions on both the observables and the probability densities. This has a number of consequences. Firstly, observables O(x + iy) will in general not be continuous across the 'seam', where we identify y j = −Y − with y j = Y + . They can therefore not be interpreted as continuous functions and a priori the Itô formula (9) does not hold (it may still hold in the sense of distributions, which should be sufficient for our purposes). Furthermore, the jump across the seam will mean that the CR equations are no longer satisfied everywhere. For the evolved observables the CR equations cannot be expected to hold anywhere exactly, as the violation that occurred initially only at the boundary gets propagated everywhere by the Langevin evolution. Similarly, the drift in the FPE is expected to be discontinuous at the seam. Therefore we have to expect that P has a jump there as well; this in turn forces us to interpret the FPE in the sense of distributions.
One might wonder whether it would not be better to limit the fluctuations in the imaginary direction by introducing a smooth cutoff; but since such a smooth cutoff function will necessarily be nonholomorphic it will destroy the formal arguments even more; we therefore stick with the simplest choice of a periodic cutoff.
We conclude therefore that the introduction of a cutoff −Y − < y j < Y + and imposing periodic boundary conditions leads to a breakdown of the formal arguments given in Sec. II. Although it is difficult to quantify precisely the effect of this, it seems reasonable to expect that it is still the behavior of P O and similar products at large |y j | that determines what is happening. For the CLE it is also clear that a very large cutoff will practically not be felt, because the system very rarely will make contact with it. This is borne out by our numerics which clearly shows convergence to the limit of infinite cutoff. For the FPE, on the other hand the issue is less clear, because there are very large boundary terms arising from the gradients of the drift across the 'seam'. In any case, for the FPE we cannot directly compare with the cutoff-free results, because these do not exist.

A. The U(1) one-link model
To understand in more detail how boundary terms affect the behaviour of complex Langevin simulations, we studied in some detail the U(1) one-link model in the hopping (HDM) approximation that was already discussed in Ref. [12] for N I = 0.

The action is
The complex drift force is correspondingly and the two components of the drift read As discussed in Ref.
A special feature of this model is that for y = c the drift is purely in x direction. If in addition N I = 0, the Langevin process will never leave the line y = c if it starts there (we emphasize that the properties discussed in this paragraph do not hold for the full U(1) one-link model, which was studied in detail in Ref. [12]). It is therefore straightforward to find an explicit solution to the stationary FPE, It follows that this model is actually equivalent to one with a real action, once we shift y → y + c and replace β by a (a and c now embody the dependence on the parameters of the model). The numerics presented below show that the line y = c is an attractor for the Langevin process; this indicates that the solution (55) is unique (with proper normalization). These properties imply that for N I = 0 the dynamics is completely understood.
When N I > 0, the presence of the repulsive fixed point is responsible for the occurrence of large excursions, which are well known to be the scourge of complex Langevin simulations. For large y the drift terms dominate over the noise, and the Langevin process is essentially just a deterministic motion; the 'classical' trajectories are given by where the complex integration constant C is related to the starting point by It is easy to see that all trajectories, except those starting on the unstable trajectories (Im z(0) = π ) are attracted to the stable fixed point at z = ic.
As an illustration we show in Fig. 1 a scatter plot of a Langevin simulation clearly exhibiting the classical orbits. There were 500 update steps between consecutive points and we used N I = 1. To enhance the classical features, this simulation was done with the noise terms in the CLE (see Eqs. (59, 60) below) suppressed by a factor 10; this is of course equivalent to replacing ǫ by ǫ/100 while multiplying the force terms by a factor of 100. This scatter plot should capture some generic features that will also be present for other choices of the model parameters.
In order to numerically study the role of boundary terms and large |y| values at nonzero N I , we introduce a cutoff in imaginary direction, placed symmetrically around y = c, i.e. −Y + c < y < Y + c, and impose periodic boundary conditions. To see the effect of this cutoff, we compute numerically the expectation values of the observables exp(ikz) for k = ±1, ±2, for various values of N I and Y , both by simulating the Langevin process and by numerical solution of the FPE. Note that these observables grow exponentially at large y. The exact values are given by where I k (a) are the modified Bessel functions of the first kind.
The Langevin process is discretized in the usual way, y n+1 = y n + ǫK y (x n , y n ) + ǫN I η y,n where η x,n and η y,n are pseudorandom numbers with zero mean and variance 2. We use periodic boundary conditions in the imaginary direction, as stated above. We also use an adaptive step size [17], choosing ǫ such that the product of ǫ|K x + iK y | ≤ 0.001. To estimate the statistical error we run 100 trajectories with independent random starting points.
To solve the FPE numerically, we employ the periodicity in x and consider the Fourier decompositions with the inverse transformations given by  The FPE can be rewritten in terms of these modes as This equation is solved numerically with a discretized time step ǫ = 10 −5 and a spatial discretization in y of δ = √ ǫ. We vary the cutoff Y in the y direction from 10δ up to 500δ (1500δ in some cases). In the x direction, with −π < x ≤ π, we use 20 − 30 points. We found that convergence was reached after 10 6 time steps, corresponding to a Langevin time t ∼ 10. We found convergence for all the values of N I and Y studied. For N I = 0, we were not able to solve the FPE numerically, due to the singular behaviour, but the solution is known analytically, see Eq. (55).
We fix the parameters of the model to be In Fig. 2 we show some examples of the real probability distribution P (x, y, t) at large Langevin time t ∼ 20 for various N I and a given cutoff Y = 0.474. The y direction goes from left to right and the compact x direction from back to front. We observe that at small N I = 0.0001, the distribution resembles the analytic result at N I = 0. The distribution is very narrow in the y direction and boundary effects are not expected to play a role. Increasing N I results in a wider distribution, and boundary effects become clearly visible. The apparent non-smooth behaviour at the edge |y − c| = Y is to be expected from the discontinuity across the seam, discussed at the end of the previous section.
After obtaining the distribution, expectation values of observables follow from Eq. (22). In Tables  The following facts can be inferred from these results: (1) CLE and FPE give rather similar results, but sometimes they differ by several σ (statistical error of the CLE simulation).
(2) All the data show a clear dependence on N I , in contrast to the conclusion of the formal arguments.
For larger N I values both CLE and FPE give results different from the exact values.
(3) The best results are generally obtained for the smallest N I . In this case there is also the weakest Y dependence, in fact no Y dependence whatsoever for N I = 0.
Obviously the presence of the cutoff and periodic boundary conditions affects the CLE and FPE in a similar way. But the N I dependence shows the failure of the formal argument, even for Y = ∞. At least for k = ±1 one has a clear case of 'convergence to a wrong limit'. The data of the CLE with Y = ∞ and k = ±2 are actually not really converged, except for N I = 0. Observing the scatter plots for different Langevin times it appears unclear whether the observables exp(±2iz) really reach an equilibrium distribution -they tend to drift out to infinity, whereas their averages, while remaining small, suffer from huge fluctuations. On the other hand the observables exp(±iz) seem to reach a stable distribution for Y = ∞ and any N I .

B. The model of Guralnik and Pehlevan
Guralnik and Pehlevan [15] studied an instructive toy model on M = R and M c = C, called GP model henceforth. Its action is and was studied in connection with PT invariant but nonhermitian Hamiltionians (where PT indicates the combined action of parity and time reversal) [26].
The action (67) leads to the drift forces There is a stable fixed point at z = i and an unstable one at z = −i. The 'classical trajectories' obtained by leaving out the noise are given by Since this is a Möbius transformation from w = tanh(βt) to z(t) the trajectories are circles. They can be imagined to emerge from the unstable fixed point z = −i at t = −∞ and go to the stable fixed point z = i as t → ∞. Those classical trajectories again can be seen clearly in the large excursions, since there the noise becomes negligible. Fig. 5 shows the result from a Langevin simulation at N I = 1.0 and β = 1.0 (in this case there are 50000 update steps between two consecutive points).
Guralnik and Pehlevan give the exact results of the first three moments at β = 1, they are They solved the discretized Langevin equation numerically, using N I = 0.001, and obtained good agreement with those exact results. We did some more and probably longer simulations at N I = 0, 0.1, 0.5, 1.0. The results for Im z are shown in Fig. 6; again they show a clear dependence on N I , in conflict with the formal reasoning. While for small N I there is agreement with the exact result, for larger N I we again have convergence to the wrong limit.
The results for Re z 2 show a similar behaviour, except that at N I = 1.0 the result is −0.951(44), with a statistical error that is huge compared to the error 0.0007 found at N I = 0.5. The data for Im z 3 show an even more dramatic failure for larger N I : they diverge for N I = 1.0, whereas for N I = 0.5 fluctuations are large. Finally we measured Im e iz , which has an exact value of 2.2624. Here divergence becomes manifest already for N I = 0.5 (but it might occur for all N I > 0).
In order to have an idea of the equilibrium distribution for N I = 0 we show in Fig. 7 a scatter plot of 50000 configurations in the complex plane. Non-Gaussian behaviour is quite clear from this plot. Noticeable is the appearance of sharp edges of the distribution, possibly indicating jumps, but no δ functions like in the U (1) case in the hopping expansion. To obtain these pictures we sampled over 60000 points taken at equal intervals of 0.5 in Langevin time. Similar distributions have been observed in Ref.
[12] for the full U(1) one-link model. The non-Gaussian character is further demonstrated in Fig.  8: the histograms for both Re z and Im z deviate strongly from a Gaussian distribution.

A. Explanation of wrong results
We will now try to interpret these findings. They are quite analogous in the two model cases and we will try to reach some conclusions that can be generalized to more realistic models. The question of convergence vs. divergence apparently depends on the values of N I , but it is more plausible that there is no qualitative difference between the different positive values of N I , only the time needed to observe the asymptotic behaviour is different. On the other hand for N I = 0 there seems to be really a qualitative difference: the distributions develop discontinuities or even δ functions; but more important is the fact that they seem to drop very rapidly in the imaginary direction.
We tentatively conclude that for N I = 0 the systems relax to equilibrium measures that show at least exponential decay in imaginary direction (in our simple U (1) model this decay is of course much stronger -the measure is zero for |y| > c).
For N I > 0 the situation is less clear, but it seems that in the U (1) model we get a decay at least like exp(−|y|), but probably not stronger than any exponential. In the model of Guralnik and Pehlevan the data suggest a power-like decay (the power appears to be near 2).

To sum up our tentative conclusions:
There is a unique equilibrium distribution P (x, y; ∞) for the Langevin processes for any N I ≥ 0, but for N I > 0 it shows limited decay, especially in the imaginary direction.
The type of decay for N I > 0 depends on the model considered; but for some observables their growth in imaginary direction may conspire with the falloff of the equilibrium measure in such a way that we obtain convergence to a wrong limit. In any case one should not expect convergence of the mean value for all holomorphic observables.
For very small N I and limited simulation time, the behaviour is of course indistinguishable from the one at N I = 0, and one may reach a quasi-convergence to the right limit, even if an infinitely long simulation would diverge.
If we try to generalize boldly from our toy U (1) model to real lattice gauge theories, we expect that for N I > 0 and for most interesting observables, such as Wilson loops, Polyakov loops etc., we have to expect boundary terms contributing even as the boundary is sent to infinity. For multiply charged loops and N I > 0 the situation could be even worse: those boundary contributions may diverge as the boundary is moved to infinity. But these problems probably will not occur for N I = 0 and do not show up at very small values of N I either, at least if simulations are not run excessively long.
We realize that definite conclusions about the falloff of the probability density P (x, y) in the y direction have not yet been reached; we intend to return to a more detailed study of this question in a future paper.

B. Practical conclusions
The conclusions for practical applications are quite straightforward: (1) In complex Langevin simulations one should use N I = 0. If one wants to do a random walk simulation or an iterate of the Fokker-Planck operator, this is not possible, but one should make sure that N I ≪ 1.
(2) Always validate the simulations by comparison with other trusted computations for parameter values where other methods are available.
(3) Wilson or Polyakov loops of higher charge will generally have much larger fluctuations. In general they are not needed for physics applications, but they may be worth looking at because they can give information about the probability distribution.