Hamilton dynamics for the Lefschetz thimble integration akin to the complex Langevin method

The Lefschetz thimble method, i.e., the integration along the steepest descent cycles, is an idea to evade the sign problem by complexifying the theory. We discuss that such steepest descent cycles can be identified as ground-state wave-functions of a supersymmetric Hamilton dynamics, which is described with a framework akin to the complex Langevin method. We numerically construct the wave-functions on a grid using a toy model and confirm their well-localized behavior.


Introduction.
The first-principle approach to solve the theory is an ultimate goal of theoretical investigations. The quantum Monte Carlo simulation for the functional integral has been a powerful ab initio technique to reveal non-perturbative features of various physical systems. Successful applications include the lattice-QCD (Quantum Chromodynamics) simulation, the lattice Hubbard model, the path integral representation of spin systems, etc.
The Monte-Carlo algorithm is based on importance sampling, so it is demanded that the integrand should be positive semi-definite, i.e., e S ≥ 0 where S is the classical action. This positivity condition is often violated in systems of our interest and then we can no longer rely on Monte Carlo simulations [1,2]. In QCD a finite baryon or quark density introduces a mixture of Hermitian and anti-Hermitian terms in the action, and then e S acquires a complex phase. In repulsive Hubbard model away from half-filling or generally in fermionic systems with spin imbalance, the sign of the integrand may fluctuate. It is also a notorious problem of the complex phase that appears from the Berry curvature in the path integral representation of spin systems and this phase cannot be removed for frustrated situations such as the XY model on the Kagomé lattice.
Moreover, to approach real-time quantum phenomena, e S is an oscillating function by definition and the sign problem is unavoidable. Although the Monte Carlo simulation is useful to compute physical observables in equilibrium in the imaginary-time formalism, the analytical continuation is necessary to access the real-time information. In general, however, the analytical continuation is a quite pricey procedure, and some additional information on the system such as the pole and the branch-cut structures would be necessary.
Many ideas have been proposed so far to overcome the sign problem and, unfortunately, applicability of each method has been severely limited. Recently new techniques to complexify the theory are attracting more and more theoretical interest, which includes the path integral on Lefschetz thimbles [3][4][5][6][7][8][9][10][11][12][13][14] and the complex Langevin approach [15][16][17][18][19][20][21]. Except for several formal arguments, theoretical foundations for the complex Langevin method are not fully established, and not much is known about its reliability [16][17][18]. In the context of real-time quantum systems, the numerical simulation works for some initial density matrices [19,20]; however, it is recently reported in Ref. [21] that the real-time anharmonic oscillator at zero temperature converges to a wrong answer with unphysical width.
In contrast, the Lefschetz-thimble method has a solid mathematical foundation at least for finitely multiple integrals [3][4][5]; however, its practical applicability is still in the developing stage. This method decomposes the original integration cycle into several steepest descent ones, called Lefschetz thimbles, using complexified field variables. Picking up a single Lefschetz thimble, one can employ importance sampling [6][7][8] and one can evade the sign problem since the oscillatory factor in e S totally disappears on each Lefschetz thimble. On the other hand, zero-dimensional model studies have exemplified the importance of structures of multiple Lefschetz thimbles [10][11][12][13]. For further applications it is needed to deepen our understanding on more aspects of the Lefschetz thimbles.
The purpose of this Letter is to shed new light on the Lefschetz thimble method in a form of the Hamilton dynamics, which was first elucidated in Ref. [5]. In this reformulation, the Lefschetz thimbles can be identified as ground-state wave-functions of a supersymmetric topological quantum system. After reviewing this modified Lefschetz thimble method, for a quartic potential problem at zero dimension, we solve the Hamilton dynamics concretely to find the corresponding wave-functions. Based on our analytical and numerical observations, we discuss advantages of this reformulation for the numerical computation.

Lefschetz thimble and SUSY quantum mechanics.
Let us consider an N -dimensional real integral as a "quantum field theory" defined by a (complex) classical action S(x). In this theory our goal is to compute an expectation value of an "observable" O(x) defined by where x = (x (1) , x (2) , . . . , x (N ) ) ∈ R N and the normalization N is chosen such that 1 = 1. The starting point in our discussion is to reformulate this theory in an equivalent and more treatable way using a complexified representation: Here 2 ∈ R and dzdz represents the integration over the whole complex plane; i.e.
The choice of the generalized weight function P (z,z) may not be unique. Indeed, a trivial example is P (z,z) = N e S(z) i δ(z (i) −z (i) ). At the cost of complexifying the variables, nevertheless, it is often the case that P (z,z) could be endowed with more desirable properties for analytical and numerical computation than the original e S(x) .
A clear criterion to simplify the integral is to find P (z,z) such that the phase oscillation can be as much suppressed along integration paths as possible, while in the complex Langevin method P (z,z) is optimized to become a real probability. To suppress the phase oscillation, let us pick up a saddle point z σ satisfying S (z σ ) = 0. The steepest descent cycle or the 2/11 Lefschetz thimble J σ of the saddle point z σ is defined with a fictitious time t as This is a multi-dimensional generalization of the steepest descent path in complex analysis, which we will refer to as the downward path. The original integration path on the real axis in Eq. (1) can be deformed as a sum of contributions on J σ weighted with an integer m σ ; i.e., In mathematics it is established how to determine m σ from the intersection pattern between the steepest ascent (upward ) path from z σ and the original integration path [3][4][5]. It is important to note that Im S is a constant on each Lefschetz thimble for the application to the sign problem [6,8].
In the following, let us restrict ourselves to N = 1 for simplicity, because the generalization is straightforward. So far, the Lefschetz thimble is constructed as a line, and let us find a two-dimensional smooth distribution P (z,z) according to Ref. [5]. For that purpose, we define the "delta-functional one-form" δ(J σ ) supported on the Lefschetz thimble so that For instance, δ(R) = δ(y)dy. Such delta-functional forms δ(J σ ) (on a Kähler manifold) have a path-integral expression from the supersymmetric quantum mechanics [22][23][24] (see also Secs. 2.8 and 4 of Ref. [5] for more details in this context). Integration (4) can be represented as Here x, p are bosonic fields and π, ψ are fermonic ghost fields, and z(t) → z σ as t → −∞.
We should note that an integration in terms of z is promoted to the path integral on z(t) for t ≤ 0, while the observable and the weight O(z(0)) exp S(z(0)) are functions of z(0) only. Let us outline how these two expressions (4) and (5) are equivalent [5,[22][23][24]. We first integrate out p(t) to get the Dirac delta function, This delta function constrains the path integral on x(t) to a gradient-flow line defining Lefschetz thimbles. Since z(−∞) → z σ , this path integral for t < 0 gives a delta-functional support on J σ . However, the delta function produces an unwanted determinant factor. As is well-known, the path integral on ghost fields π(t), ψ(t) for t < 0 can eliminate that factor as Now, we obtain an integration over surface variables x(0), ψ(0), and denote them by x, ψ. Locally, the Lefschetz thimble J σ can be expressed as zeros of a certain function f , then we 3/11 can find that the path integral (5) eventually gives which is nothing but the local expression of the original integration (4). Going back to (5), this shows that the so-called residual sign problem comes from the fermionic surface term ψ 1 (0) + iψ 2 (0) because one can identify ψ i (0) = dx i as above. Importantly, with these added fields, p i , π i , ψ i , the action is BRST exact under a transformation;δx i = ψ i ,δψ i = 0,δπ i = −ip i ,δp i = 0. By definition the nilpotencyδ 2 = 0 is obvious. Thanks to the boundary fermionic operator in (5), the surface term is BRST-closed so long as the observables are holomorphic. This makes a sharp contrast to the complex Langevin method that could also acquire the BRST symmetry but it is violated by the surface term. Because of the BRST symmetry we can add any BRST exact terms without changing the original integral, and it is useful to insert εi 2 dt p 2 i . In summary, the effective Lagrangian that describes the fictitious time evolution is given by the following topological theory: which is nothing but a Legendre transform of an effective Hamiltonian: with [x i ,p j ] = iδ ij and {π i ,ψ j } = δ ij . The fermion number F =π 1ψ1 +π 2ψ2 is a conserved quantity of this Hamiltonian. After the time evolution from t = −∞ only the ground state with the lowest energy eigenvalue remains, so that the generalized weight is given by P (z,z)dzdz = Ψ(z,z) ∧ e S(z) dz, where Ψ(z,z) is the ground state wave-function and converges to δ(J σ ) in the limit ε i → +0. Note that the weight factor exp S(z) is necessary in this formula, since the wave function designates only the integration cycle J σ . We can further simplify this Hamilton problem by choosing ε = ε 1 = ε 2 . Performing the conjugate transformation Ψ = e −Re S/ε Ψ , the first derivative terms are eliminated as This describes supersymmetric quantum mechanics with the superpotential Re S [5]. Before applying this method to an interacting model, let us convince ourselves of perturbative correctness. For that purpose, we consider the simple Gaussian case: where x ∈ R is a one-component variable and ω ∈ C. One can regard this Gaussian integral as an elementary building-block of perturbative quantum field theory; in a non-interacting theory of the scalar field φ, the action is decomposed into for each Fourier mode. In this case we can immediately find that the bosonic wave-function should be the ground state of harmonic oscillator with the ground state energy |ω|. The 4/11 fermionic ground state energy −|ω| cancels the bosonic energy thanks to supersymmetry. As a result, unoccupied fermions point the direction of the Lefschetz thimble; i.e., the supersymmetric vacuum belongs to the F = 1 sector. Let us see this in the simplest example, ω ∈ R and ω > 0. The fermionic part of the Hamiltonian is ω 2 ([π 1 ,ψ 1 ] − [π 2 ,ψ 2 ]), and thus the 1 and 2 fermions are unoccupied and occupied, respectively. This leads to the fermionic ground state energy −ω, which cancels the bosonic ground state energy ω. Therefore, the unoccupied fermion is tangent to the Lefschetz thimble J = R.
The final result of P 0 (z,z) for ω ∈ C together with e S · e − 1 ε Re S is which reproduces the original integral (1) with the action (12). To see this for a polynomial O(x) it is sufficient to require z 2 = 1/ω, which is nothing but the free propagator and can be explicitly confirmed with Eq. (13). In order for the exponentially fast convergence of P 0 , the parameter ε needs to be 0 ≤ ε < 2. We here emphasize that the theory is equivalent to the original (1) for 0 ≤ ∀ε < 2 and the conventional Lefschetz thimble method is retrieved in the ε → 0 limit. Actually, in this limit of ε → 0, only a path of 2|ω|zz − ωz 2 −ωz 2 (≡ {2Im( √ ωz)} 2 ) = 0 contributes, which is nothing but a condition to guarantee Im S 0 (z) = 0 on the Lefschetz thimble. At finite ε, this restriction is smeared and P 0 may have a distribution around J with a width of order of ε where a complex phase arises in general. The non-positivity of P 0 is a big difference of this Lefschetz-thimble approach from the complex Langevin method, in which the distribution must be semi-positive definite.
In the model with quartic interaction, S = − ω 2 z 2 − λ 4 z 4 , the vacuum structure is drastically different from that of the Gaussian case (13). With λ = 0, there are three classical saddle points and the Morse index for each saddle point is 1. The Witten index is tr(−1) F = −3, and thus there are three supersymmetric vacua in this interacting model for any ε [25,26]. In the path integral expression (5), we can distinguish these three vacua by specifying the boundary condition at t = −∞, as we will discuss below in detail.
We define our zero-dimensional model with the quartic interaction by and hereafter, we will specifically choose the model parameters as This choice has been motivated by the application to real-time problems. As we already mentioned, Re ω corresponds to the width Γ or in the i prescription, and Im ω corresponds to −(k 2 − m 2 ). Therefore, the parameters (15) represent a situation at k 2 = m 2 + Γ in quantum field theory.
When ε is small enough, the Hamilton dynamics should become equivalent to the conventional formulation of the integral on the Lefschetz thimble, which means that P (z,z) should have a peak along the Lefschetz thimble only. This makes the analytic treatment much accessible since we do not have to solve the quantum mechanical problem in this zero-dimensional toy model. To identify the Lefschetz thimble, we should integrate the flow 5/11 (3) and find the upward and the downward paths. We show the numerical results in Fig. 1 for parameters slightly changed from Eq. (15). For our theory (14) three saddle points are located at z 0 = 0, z ± = ± −ω/λ. The latter is, for our choice of parameters (15), z ± = ±(0.897 + 0.372i), and S(z ± ) = −1/3. We can see that ω = 1 − i is a critical value at which the destinations of the downward flows from the saddle points change drastically, that is known as the Stokes phenomenon, as is clear from two panels for ω = 1 − 0.9i (left) and ω = 1 − 1.1i (right) in Fig. 1. In general we can show that the Stokes phenomenon occurs at Im (ω 2 /λ) = 0, and when λ is pure imaginary, this condition gives arg (ω) = −π/4. Importantly, not only the downward flows but also the upward flows change and the intersection number of the original and the upward paths changes accordingly [4,5]. In the case with ω = 1 − 0.9i only one upward path from z 0 crosses the real axis as seen in the left of Fig. 1, and so the Lefschetz thimble going through z 0 contributes to the integral. In the case with ω = 1 − 1.1i, on the other hand, three upward paths from z 0 and z ± all cross the real axis, and all three Lefschetz thimbles contribute to the integral.
Keeping the potential application to the real-time physics in mind, we need deeper understanding on the Stokes phenomenon. In fact, it occurs at k 2 = m 2 + Γ and so the thimble structure may fluctuate depending on the frequency k 0 , while in Euclidean theory it never happens because k 2 − m 2 is always negative (except for an unstable potential with m 2 < 0). It should be an interesting future problem to clarify the treatment of the Stokes phenomenon in the real-time systems [10].
For the same model with a different set of parameters, the Stokes phenomenon has been discussed in Refs. [27,28] and the probability distribution P (z,z) in the complex Langevin method has been numerically computed. The important insight obtained there is that P (z,z) in the complex Langevin method looks localized but has a power decay at large |z|, which causes a convergence problem. Hence, it would be an intriguing question how P (z,z) in the modified Lefschetz thimble method for ε = 0 should behave especially at large |z|.
As an application of the modified Lefschetz-thimble method with a regulator ε, let us compute the wave-function Ψ for the ε = ε 1 = ε 2 case, from which P (z,z) is to be constructed immediately. We can readily find the eigenstate in terms of ψ i for the last term in the Hamiltonian (11). By restricting ourselves to the F = 1 sector, we can define the effective potential in a form of 2 × 2 matrix-valued function that amounts to and then the Hamiltonian is Let us solve the ground state of the above H eff at finite ε, and we specifically adopt ε = 1 unless stated explicitly. When we solve the Hamilton dynamics, we should set initial conditions to select proper Lefschetz thimbles out. We choose the semiclassical ground state in the limit ε → +0 as our initial condition. Since Ψ = e ReS/ε Ψ and Ψ (zσ) → δ(J σ ) in ε → +0, Ψ (zσ) ∼ e ReS/ε δ(J σ ) for small ε. Let us consider the Lefschetz thimble around z 0 for instance, then its tangential direction at z 0 is given by x 2 = tan (π/8)x 1 as seen also from Fig. 1. The initial wave-function is thus proportional to The bosonic part is well localized at z 0 , which justifies the following choice of the initial wave function, Ψ (zσ) (t = −∞) = δ(z − z 0 )(− sin(π/8)dx 1 + cos(π/8)dx 2 ). Similarly, we can fix the initial wave function for z ± as Ψ (z±) = δ(z − z ± )(cos(π/8)dx 1 + sin(π/8)dx 2 ). For the numerical procedure we smeared the delta function in the initial wave function by a Gaussian as δ(z − z σ ) → exp{−20[(x 1 − x 1σ ) 2 + (x 2 − x 2σ ) 2 ]}. Then we discretized x 1 and x 2 from −2.5 to +2.5 with dx = 5 × 10 −2 . We then numerically integrate d dt Ψ (zσ) = −H eff Ψ (zσ) using the Euler method with dt = 10 −4 until the wave-function converges. The convergence is fast and stable and the wave-function hardly changes after t = 1 ∼ 2. We also mention that we utilized the Crank-Nicolson algorithm to improve the numerical stability when we compute the Laplacian.
With this prescription, we find three independent ground state wave-functions Ψ (zσ) = Ψ (zσ) 1 Fig. 2. The figures (a,b) in Fig. 2 show two components of Ψ (z0) and the figures (c,d) show those of Ψ (z+) . Since our system is symmetric under the reflection z → −z, we can easily read Ψ (z−) from the figures of Ψ (z+) . Remarkably, these ground-state wave functions are linearly independent from one another, which means that supersymmetry is unbroken. We have also verified that the Dyson-Schwinger equation, is satisfied within 1% accuracy. This clearly shows that the supersymmetric quantum mechanics provides a suitable framework to compute Lefschetz thimbles. Let us comment on the convergence of the wave-functions in the present modified Lefschetz thimble method. In the effective potential (16), the first term gives the dominant binding potential in our toy model for large |z| because the first term is a polynomial up to the sixth   Fig. 3 SUSY ground stateΨ (z0) , which is obtained starting with the saddle point z 0 but with the initial relative weight equal to both components. order, while the second term is up to the quadratic order. Therefore, the convergence in the present numerical approach is quite improved and there is no problem of the power-decay unlike the complex Langevin method. One might think that the remaining e iIm S is still oscillating, but this is not a problem practically. We have numerical verified that the effect of e iIm S is very small in the region where the profile of Ψ (zσ) is localized.
We also checked the robustness of the numerical results against small variations of the center position and the smearing width in initial conditions. If the smearing width of the initial wave function is changed, the overall normalization would be changed naturally, but once we normalize the wave-functions in the same manner (in Fig. 2 we normalized them 8/11 ) 2 ] = 1), then we eventually get the same result. Such insensitivity implies that each wave-function is well localized and the overlap at the saddle point is small. However, the overlap at the saddle point is not completely zero. To see this effect, let us skew the relative weight of the initial condition so that the initial wave-functions are not orthogonal. In Fig. 3, we showΨ (z0) 1 andΨ (z0) 2 starting with the same relative weights, namely,Ψ (z0) (−∞) = δ(z − z 0 )(dx 1 + dx 2 ) at the initial time, for the demonstration purpose. The resultΨ (z0) in Fig. 3 is given a natural interpretation as a superposition of three ground states shown in Fig. 2. Because of this overlap among wave functions, our prescription for the initial wave function needs further refinement in order to extract one Lefschetz thimble at ε = 1.
In order to see the importance of such refinements from another viewpoint, we checked behaviors of distributions P for small ε. Let us set ε = 0.2 with ω = 1 − i and λ = 1.5i, then we obtain Figs. 4 (a) and (b). Here, the normalization is given by d 2 x P = 1. Although it is localized around the saddle point z 0 = 0, its shape is still far from a one-dimensional line shown in Fig. 1. In order for comparison, we show the result also for ε = 0.2 with ω = 1 + i and λ = 1.5i in Figs. 4 (c) and (d). For this case, P is well localized to a one-dimensional line, which is nothing but the Lefschetz thimble J 0 at this parameter. We numerically observed that the weight function P is well localized around a Lefschetz thimble if Stokes phenomena do not happen and ε is sufficiently small. Expectation values of operators, such as z 2 , also give correct numbers within the numerical accuracy with such parameters. However, as ε becomes larger, wave functions spread as two-dimensional distributions, and expectation values of operators do not necessarily give correct values. Since the Dyson-Schwinger 9/11 equation (19) is satisfied, this problem must come from superpositions of the wave function with other ground states.
It is an important future study to clarify dependence on this initial condition more systematically, in order to take an appropriate linear combination of these wave-functions giving each Lefschetz thimble. This would be a key step to study how the Stokes phenomenon is realized at ε = 1. It also opens a new possibility to take into account multiple Lefschetz thimbles, since ground-state wave functions can be superposed by changing the initial condition.

Discussions and Conclusions.
In this Letter, the supersymmetric reformulation of the Lefschetz thimble integration was studied and its practical computation was discussed. The Lefschetz thimbles are now regarded as the ground sates of a supersymmetric Hamiltonian. Our computational scheme for the Hamiltonian system is essentially the same with that for the Fokker-Planck equation in the complex Langevin method. Those supersymmetric wavefunctions are numerically computed for a zero-dimensional toy model with the classical action S = S 0 + S int where S 0 = − ω 2 x 2 and S int = − λ 4 x 4 . Since the Morse indices of all the saddle points are the same, the number of saddle points must be the same as that of linearly independent ground states.
For the Gaussian model with S 0 only, the above-mentioned situation is clearly true, which is explicitly checked by computing the wave-function analytically and constructing one supersymmetric ground sate. From this almost trivial example, one can learn an important lesson about the two-dimensional smooth distribution P 0 (z,z). In the "semi-classical" limit of ε → 0, we recover a delta functional support along the Lefschetz thimble in the original formulation. For non-zero ε, the phase oscillation arises away from the Lefschetz thimble, and the formulation nevertheless reproduces the correct expectation value.
For the interacting model with S int , we performed numerical computations for the supersymmetric quantum mechanics, and confirmed the existence of three linearly independent ground states by restricting ourselves to the F = 1 sector. Since the wave-function in the F = 1 sector consists of two components, we must set them in the proper initial conditions. We started from a localized wave-function in the vicinity of each saddle point, and our numerical computation shows a remarkable stability under small modifications on the initial conditions. This reflects the fact that all the saddle points are attractive unlike the complex Langevin method. At the same time we also found substantial dependence on the initial relative weight of these two components. This clearly indicates that we need a careful refinement of the initial condition to use the modified Lefschetz thimble method for numerical simulations, but it also opens a new possibility for some convenient scheme to take into account multiple Lefschetz thimbles.
In the case of the Fokker-Planck system, the ground state is often uniquely determined and the initial condition dependence does not appear as in the case in ordinary quantum mechanics without any special symmetry. The Fokker-Planck operator of the complex Langevin method can be written as a functional integral in a similar manner, and it also shows the same type of the BRST symmetry. There are, however, several important differences in these two formalisms: First of all, the diffusion term in the complex Langevin equation does not give the gradient flow because the sign is different. Because of this difference, the BRST invariance cannot be promoted to the supersymmetric system satisfying 2H = {Q,Q} with 10/11 supercharges Q,Q. Moreover, the fermionic sector should be restricted to F = 2 instead of F = 1. Therefore, it is not straightforward to relate these two formalisms yet. If one could establish a firm relation between two methods, it would be a great progress and the present modified Lefschetz thimble could provide us with an opportunity for a fully complementary approach to the sign problem together with the complex Langevin method.