Complex saddle points and the sign problem in complex Langevin simulation

We show that complex Langevin simulation converges to a wrong result, by relating it to the Lefschetz-thimble path integral, when the path-integral weight has different phases among dominant complex saddle points. Equilibrium solution of the complex Langevin equation forms local distributions around complex saddle points. Its ensemble average approximately becomes a direct sum of the average in each local distribution, where relative phases among them are dropped. We propose that by taking these phases into account through reweighting, we can solve the wrong convergence problem. However, this prescription may lead to a recurrence of the sign problem in the complex Langevin method for quantum many-body systems.


Introduction
Precise analysis of thermodynamic properties of a quantum many-body system, in particular, precise determination of its phase diagram is one of great challenges in theoretical physics. An ab initio simulation based on lattice field theory, in particular, so called Monte Carlo simulation is the most powerful tool for this. In many interesting cases, however, Monte Carlo simulation is hindered by the notorious sign problem. The importance sampling, using the Boltzmann weight e −S , breaks down when the action becomes complex. In hadron physics, lattice quantum chromodynamics (QCD) simulation suffers from the sign problem at finite quark densities [1,2], which is important to study quark matter inside neutron stars [3,4]. The sign problem occurs also in condensed matter systems [5,6,7]. Important examples are the fermionic Hubbard model away from half-filling, and geometric frustration in spin systems. A method to overcome the sign problem attracts a broad interest for application to the aforedescribed quantum many-body systems.
There have been a lot of attempts to tackle the sign problem. Among them, idea of complexification of the integration variables is one promising way to solve the sign problem. Theoretical attempts along this line are classified into two approaches, that is, the Lefschetz-thimble and the complex Langevin methods. The Picard-Lefschetz theory gives a generalization of the steepest descent method, and Lefschetz thimbles are steepest descent paths in the extended complex plane [8,9,10]. This method is formulated on rigorous mathematics, but it needs some approximation when applied to quantum many-body systems [11,12,13,14,15,16]. On the other hand, the complex Langevin method is an extension of the Langevin equation to a complex Boltzmann weight [17,18,19,20]. The numerical implementation of this is possible based on lattice field theory. The complex Langevin method has been widely applied from condensed matter to hadron physics [21,22,23,24,25,26]. There is a formal proof [27,28] on the correctness of the complex Langevin method, where it has been shown that the complex Langevin method correctly gives physical observables if the distribution obtained from the Langevin equation damps exponentially fast around infinities and singular points. This method is, however, known to give wrong results for some cases, where distribution does not show the exponentially fast decay, and thus the formal proof cannot be applied (For recent discussions, see also [29,30,31,32,33,34,35]). Therefore, it is important to unveil what properties of the classical action cause the wrong convergence of the complex Langevin method.
In this paper, we show within the semiclassical analysis that complex Langevin simulation converges to a wrong result, when path-integral weight at complex saddle points has different phases. This includes the case of the breakdown due to a singular drift term, e.g., the lattice QCD at finite density. We reveal that complex Langevin simulation breaks down more generic case where the Langevin drift term has no singular point. With the help of semiclassical analysis, we find that reweighting by the complex phase can partially solve the wrong convergence problem. However, the reweighting leads, in general, to a severe cancellation of the reweighting factor in many-body systems, which is nothing but a sign problem in terms of the complex Langevin method.

Complex Langevin method and its failure
For simplicity, we discuss an oscillatory integral of one variable x, which can be extended to multiple integrals in a straightforward way, where Z is the normalization factor. The action S (x) is complex valued in general, which makes the Monte Carlo simulation of Eq. (1) difficult because of the sign problem. One proposal to calculate Eq. (1) for a complex valued action is the so-called complex Langevin method [18,19,20]. In this method, we solve the Langevin equations for complex values z = x + iy along the fictitious time direction θ, where η(θ) is real Gaussian noises satisfying η(θ) η = 0, and η(θ)η(θ ) η = 2δ(θ − θ ). Since the action S (x) is complex, the right-hand side of Eq. (2) is also complex. Thus, complexification of the variable z is unavoidable, which is the reason why this method is called the complex Langevin method. In the real Langevin method i.e., if S (x) is real, the ensemble average O(x η (θ)) η can be shown to converge to Eq. (1) as θ → ∞. It has been shown that the complex Langevin method also converges to Eq. (1) for the complex action S (x) when the tail of distribution obtained from Eq. (2) damps exponentially fast [27,28]. However, it has not been understood yet that what behavior is required to actions for the success of the complex Langevin simulation. We first show based on the semiclassical analysis that the complex Langevin method gives wrong results if there are several dominant saddle points with different complex phases. After that, we propose a new prescription to evade this breakdown. Ito calculus shows the following: If the expectation value of a holomorphic operator O(z η (θ)) converges as θ → ∞, the derivative of O(z η (θ)), O(z η (θ)) = O (z η (θ)), must satisfy the Dyson-Schwinger (DS) equation, Here the argument θ = ∞ of z is omitted. When = 0, the DS equation can be solved by complex saddle points S (z σ ) = 0 (σ ∈ Σ). Even at finite , the contour integral on the steepest descent path J σ around each z σ solves the DS equation [36,37,38], In general, any solutions of the DS equation are represented by a linear combination of the contour integrals on the steepest descent paths [38]. Therefore, ensemble average of a holomorphic operator O(z η ) at θ → ∞ can be represented as Here, d σ is a complex number. The Lefschetz-thimble method [9,10,39] is useful to connect those steepest descent integrals with the original one (1). If and only if d σ is an intersection number K σ , R between a steepest ascent path K σ and the original contour R, the original integral (1) is recovered. We analyze Eq. (5) in the semiclassical limit → +0. For this, we expand S (z) around each complex saddle point z σ as Let us first analyze the left hand side of Eq. (5). If Re ω σ > 0, the solution of the equation of motion (2) can converge into z σ as θ → ∞ in → 0. On the other hand, it cannot converge for Re ω σ ≤ 0. Then in the semiclassical approximation, we have where c σ ≥ 0, and c σ = 0 if Re ω σ ≤ 0. Next let us analyze the right hand side of Eq. (5). In the semiclassical approximation, the integral along the thimble becomes Also the denominator Z Z semi can be evaluated using the semiclassical analysis by setting O(z σ ) = 1 in the above discussions. Now, from the comparison of the both sides of Eq. (5) for arbitrary operators, we reach for dominantly contributing saddle points. Note that c σ ≥ 0. However, d σ needs to be an integer K σ , R to recover (1). These two statements contradict with one another in general.
As a result, we can conclude the following at least for semiclassical analysis: The complex Langevin method cannot reproduce the original integral (1) if there are several dominant saddle points with different complex phases. Equation (9) is completely sure for dominant saddle points. The above contradiction must be taken into account if the dominant saddle points have different complex phases. In other words, complex Langevin method may fail if there is some relative phase between the dominant saddle points. For subdominant saddle points, the ambiguity of Borel resummation of large order perturbations can give nontrivial cancellations [40,41,42,43,44,45,46,47]. Therefore, we cannot judge from our argument whether the complex Langevin method gives a correct result if there is only one dominant saddle point. This subtlety needs further studies. For a Gaussian action, it is easy to check that Eq. (9) is satisfied and the complex Langevin method works well. On the other hand, there exists a model with power-law tail, where the complex Langevin method does not work but there is only one dominant saddle point (see, e.g., [33]).
Let us give a few comments on previous studies. There is a formal proof [27,28] on the correctness of the complex Langevin method, but it relies on several nontrivial assumptions 1 . Combined with a recent study [34], they have shown that the formal proof breaks down if the complex Langevin distribution does not decay exponentially fast around infinities and singular points. Our analysis suggests without accessing details of the complex Langevin distribution that the breakdown happens if the dominant complex saddle points have different phases. Even in many-body systems, we can obtain saddle points by numerically solving Eq. (2) without random noises. This situation would naturally bring us to the conjecture that the complex Langevin distribution has a polynomial tail around infinities or singular points if several dominant saddle points contribute with different phases. It would be an important future study to check this conjecture in order to achieve a deeper understanding of the complex Langevin method.

Prescription
Let us propose a prescription to circumvent this inconsistency. We denote the equilibrium distribution of the complex Langevin method by P. The expectation value is given as In the semiclassical limit, P will be represented, by using a sum of distributions P σ localized at complex saddle points z σ , as P = σ P σ , which gives the expectation value (7) i.e, By defining (nonholomorphic) functions χ σ satisfying χ σ P τ δ στ P τ , we define a phase function Θ by If P σ does not overlap with others, χ σ can be chosen as the characteristic function of supp(P σ ). This is not true in general, and we must find χ σ satisfying the condition with a good approximation. The expectation value of a holomorphic operator O(z) is given, by reweighting with Θ, as Even when the random noise or equivalently correction is included, so long as P is well localized around each saddle point, this prescription seems to work nicely. Note that this replacement does not break DS equations (3) so far as the semiclassical analysis is valid. Now our question is "What d σ , or c σ , is adopted in the complex Langevin method?" If Re ω σ ≤ 0 means K σ , R = 0, the following d σ is consistent with c σ ≥ 0: If this is true, the complex Langevin method gives an extension of the so-called phase quenched approximation to include complex saddle points: We adopt it as a working hypothesis in the following sections, although this is not the unique solution for consistency. Using this hypothesis, the phase function Θ is given by and the reweighting formula (13) is available for practical use 2 .

Numerical simulation
We test our proposal by applying it to two models with and without a singular drift term. We numerically solved (2) with the fictitious time step ε = 5.0×10 −6 and 5.0×10 −7 for the models with and without the singular drift term, respectively. We adopted a higher order algorithm [49]. Errors were estimated by using the jackknife method, and each quantity is computed by using 5.0 × 10 5 configurations. Below we set = 1.

One-site fermion model
First, as a nontrivial example with a singular drift term, we analyze a one-site fermion model. This is the simplest model to suffer from the sign problem same as that in lattice QCD simulations [50,51,52,53]. To demonstrate that the modified complex Langevin method can simulate the Silver Blaze like feature [54] in the one-site model is a good landmark to show its applicability to the sign problem in many-body systems.
After introducing a Hubbard-Stratonovich field ϕ, we explicitly integrate out the original fermionic fields. The partition function reads [50] Z = dϕ bg e −S (ϕ bg ) , with the action, where ϕ bg = β 0 dτϕ(τ)/β is the zero Matsubara mode of ϕ. U(> 0), µ, and β = 1/T are the on-site repulsive interaction, chemical potential, and inverse temperature, respectively. We dropped nonzero Matsubara modes of ϕ, since they do not couple to µ [50]. The auxiliary field ϕ bg is related to the fermion number density n by where we used the equation of motion to obtain the last expression. The integral (17) is analytically calculable, but instead, we shall apply the complex Langevin method. Due to the logarithmic term, the action has infinitely many saddle points, which appear in the period of 2πT . Since the Lefschetz-thimble method is valid even with these logarithmic singularities [55, 56], all the discussions in previous sections are available in order to conclude the failure of the complex Langevin method.
In the large βU limit, the saddle points ϕ m are given as [50] with m ∈ Z. The classical action at ϕ m reads [50] S 0 − βU In Eqs. (22) and (23), we have calculated only the m-dependent leading terms in the large βU expansion. In this model, β −1 plays a role of but the classical action (18) depends on it in a nontrivial way. Therefore, Re[S m ] becomes a function of β and m, which make difficult to judge the dominance of saddle points. According to Eq. (22), saddle points with m 2 βU would give dominant contributions in the large βU limit. Thus, the zero temperature limit, which corresponds to the classical limit in Sec. 2, is not described by the unique saddle point and the condition d σ = K σ , R is not trivially recovered. According to Eq. (23), these different saddle points have different complex phases, and thus the complex Langevin simulation may fail except for special cases.
We show the fermion number n as a function of the chemical potential µ in Fig. 1. The satndard complex Langevin method predicts the wrong linear µ-dependence. This wrong behavior is also obtained from the mean field or the one-thimble approximation [50]. To find a reweighting factor, we use approximate expressions on the saddle points in the leading order of the large βU expansion in Eqs. (21)-(23) [50]. The saddle points are in between singular points of the logarithm ϕ s m = i µ + U 2 + T (2πm + 1/2), namely, Re ϕ s m < Re ϕ m < Re ϕ s m+1 . The distribution generated by solving the complex Langevin equation is localized arournd ϕ m and decays by a power law as it is getting close to ϕ s m along the real part direction. For the imaginary part direction the distribution exponentially decays. Then we put χ m (ϕ bg , ϕ bg ) = θ(Re ϕ s m < Re ϕ bg < Re ϕ s m+1 ). The residual sign coming from ω σ turns out to be negligible for βU = 30. Now Θ is given explicitly as Θ = m∈Z e −2πi(µ/U+1/2)m θ(2m − 1 < β Re ϕ bg /π < 2m + 1), (24) where θ(x) is the step function. We also show the fermion number after reweighting in Fig. 1. The result becomes much better, but we may still need an improvement of Θ for exact agreement. The number density seems to linearly decrease in each plateaux as chemical potential increases. This behavior is incorrect from the view point of the thermodynamics stability since the compressibility must be non-negative. There might exist the physics not included in our weighting prescription.
We show the average phase factor Θ η as a function of the chemical potential µ in Fig. 2. It becomes small near jumping points of n at µ/U = 0 and 1, and is getting close to one near the half filling µ/U = 1/2. If we apply the conventional reweighting by the Monte Carlo method to the original integral (17), however, the severe sign problem appears for every µ/U > −1/2 [50]. The cancellation of the phase function in the modified complex Langevin method is milder than that in the reweighting by the Monte Carlo method. The same cancellation may happen near phase transition points in many-body systems. If Θ η becomes exponentially small as the system size increases, it is also true that the sign problem is still obstinate in the complex Langevin method.

Double-well potential model
Next, we consider a model without a singular drift term, whose action is given by with α > 0. This action has three saddle points on the complex plane. Only two of them have positive Re ω σ , and contribute to the semiclassical analysis. These two saddle points (z 1 and z 2 ) are, respectively, located on the first and second quadrant planes (Re z 1 > 0 and Re z 2 < 0). They have different complex phases, except when α = 0. The complex Langevin simulation may fail at finite α from our semiclassical analysis given in Sec. 2.
The distribution of this model seems to have the power law behavior. We show the partially integrated distribution P y (y) = dx P(x, y)/ dxdy P(x, y), and its fifth moment y 5 P y in Fig. 3. The distribution may behave as P y ∼ y −5 at y → ∞. The power law implies that the expectation value of a higher power of z, e.g., z n (n ≥4) diverges 3 . Thus the complex Langevin simulation apparently breaks down, as expected. Remark here that since z n (n ≥1) does not satisfy the DS equation (3) if the power law exponent is true, our argument based on the DS equation in Sec. 2 is no longer available. Nevertheless the complex Langevin method actually breaks down, and our prescription works well for lower dimensional operators as is seen in the following. We show the expectation values of iz and (iz) 2 as a function of α in Figs. 4 and 5. The complex Langevin simulation converges to a wrong result (red squares.) Based on our prescription, we put χ 1 (z, z) = θ(Re z) and χ 2 (z, z) = θ(−Re z). The result of the reweighting is also shown in Figs. 4 and 5 with blue circles. The reweighting works perfectly, and we resolve the wrong convergence problem. This is also true for (iz) 3 . For a diverging higher power of z, (iz) n (n ≥4), our prescription does not work, and the expectation values suffer from the large fluctuations before and after reweighting. · · · · · · · · · · · · · · · · · · · · Á Á Á Á Á Á Á · · · · · · · · · · · · · · · · · · · · Á Á Á Á Á Á Á

Concluding remarks
We have analytically shown within the semiclassical approximation that complex Langevin method gives wrong results, when there are several dominant saddle points with different complex phases. Since the interference of these complex phases is an essential ingredient to understand the Silver Blaze phenomenon [50], the usual complex Langevin method might not be reliable in order to tackle the cold and dense nuclear matters. Moreover, this interference is also of great importance in order to study the dynamical phenomena, such as a particle production, using the real-time path integral [58,59,60]. For more general situation where the semiclassical analysis breaks down, we need further study to show the failure of the complex Langevin method.
The next problem is to modify the distribution so as to reproduce the expectation values in the original theory. We proposed a reweighting prescription by introducing a working hypothesis, which is consistent with the semiclassical analysis. This must be justified or revised in future study. Also the correct treatment of subdominant saddle points must be clarified.
Our prescription is numerically confirmed for two models with and without singular drift terms. In particular, we succeeded to simulate the nonanalytic behavior of the one-site fermion model at low temperatures.
If our prescription were proven or revised, the modified complex Langevin method could provide a way to perform numerical simulations on multiple Lefschetz thimbles. However, it requires us to get complete knowledge on complex saddle points to assign correct phase function. Furthermore, our prescription causes a large cancellation of relative phases among saddle points, although it is somewhat milder than that of the conventional reweighting by the Monte Carlo method. This implies the sign problem possibly occurs in the modified complex Langevin method. To find more efficient prescription must be an important future study.