Lefschetz Thimbles and Wigner Functions

The major difficulty for the Feynman Path Integral Monte Carlo (PIMC) simulations of the quantum systems of particles is the so called “sign problem”, arising due to the fast oscillations of the path integral integrand depending on the complex-valued action. Our aim is to find universal tech-niques being able to solve this problem. The new method combines the basic ideas of the Metropolis and Hasting algorithms and is based on the Pi-card-Lefschetz theory and complex-valued version of Morse theory. The basic idea is to choose the Lefschetz thimbles as manifolds approaching the saddle point of the integrand. On this thimble the imaginary part of the com-plex-valued action remains constant. As a result the integrand on each thimble does not oscillate, so the “sign problem” disappears and the integral can be calculated much more effectively. The developed approach allows also finding saddle points in the complexified space of path integral integration. Some simple test calculations and comparisons with available analytical results have been carried out.


Introduction
One of the main difficulties for the Path Integral Monte Carlo (PIMC) simulation of the quantum systems of particles is the so called "sign problem". The "sign problem" arises in the simulations of the Wigner and Feynman path integrals, describing quantum systems and the finite density quantum chromodynamics due to the fast oscillations of the integrand defined by the complex-valued action. This integrand does not give a real and positive Boltzmann-like weight (for example, by the Wick rotation) to resort to the traditional Monte Carlo methods. The so-called reweighting algorithm is highly ineffective when the imaginary part of the action becomes very large, because one needs to take a sample from a configuration space, where the weights of nearby configurations have almost the same amplitudes but very different phases.
There have been many proposals to circumvent the "sign problem". Basic possible approach to this problem is to consider the variables, which are assumed to be real in the original formulation, to be complex and to extend the cycle of path-integration to a complex space in order to achieve better convergence.
So long under the typical physical conditions the integrand is holomorphic in the new complex variables and the final value of the path integral is unchanged by Cauchy's theorem. Making use of the Picard-Lefschetz theory and a complex version of Morse theory we can select the cycle approaching the saddle point at the path-integration, where the imaginary part of the complex action stays constant (Lefschetz thimbles) [1] [2] [3]. Since the imaginary part of the action is constant on each thimble, the "sign problem" disappears and the integral can be calculated much more effectively.
However, since different thimbles are strongly separated, one needs to develop method allowing incorporating contributions from all relevant thimbles. One of the natural ways to incorporate the relevant thimbles [4] [5] is to use a gradient flow of action starting from the original real space of integration and creating a new manifold at a finite flow time, which is equivalent to the integration over original real space. In particular, when the flow time approaches infinity, the new manifold is composed of Lefschetz thimbles. In practise, as long as the finite flow-time is large enough then the "sign problem" will be alleviated, as integrals turn into integrals of an oscillating function with decaying amplitude. However reducing the amount of flow time may also reduce the effectiveness against the "sign problem" and the multimodal problem simultaneously.
The alternative approach of sampling at calculation of the path-integrals on thimbles is based on the making use of the complexified Langevin equation [6] [7] [8]. Of course, application of the noise leads to departures from the thimble that accumulates and needs to be corrected. Numerically, this procedure can be made stable.
Besides, the Langevin algorithm, other algorithms have been proposed: an Hybrid Monte Carlo algorithm [9], that however is essentially as expensive as the Langevin algorithm; two Metropolis algorithms [4] [10] that are simpler and faster, but have the risk of poor acceptance for large systems and an alternative algorithm [11] that ensures a control of the thimble at the price of limited scalability.
In this work we present the new Metropolis-Hastings algorithm for searching critical points and subsequent generation of the Lefschetz thimbles. The algorithm allows to find the saddle points and to sample initial conditions for downward flows from vicinity of the saddle points. The developed approach combines the basic ideas of the Metropolis and Hasting algorithms. So to in-crease efficiency of numerical procedure we have separated the Markovian transition on the complex plane in two sub-steps: the proposal and the acceptance-rejection. The proposal probability allows proposing a new state for given one. Here we are free in choosing proposing probability as it affects only the efficiency of the sampling the main contribution to the integrals and does not change the final result of calculations. The acceptance distribution is the conditional probability to accept/reject the proposed state.
The integrals involved in our test calculations are one variable integral and can be performed analytically or independent numerically. However, it provides an interesting benchmark which can be seen as a limiting case of more realistic path integrals. It is non-trivial from the point of view of a Monte Carlo integration. For comparison with analytical results we present some calculations for the Airy function and some results on the Fourier transform of basic 1D factors comprising the Wiener path integral representation of the Wigner function in phase space. It also provides a case where different aspects of our approach can be clearly visualized.
This paper is organized as follows. In Section 2 we remind the basic ideas of

Wigner Approach to Quantum Mechanics
We consider a one-dimensional quantum-mechanical system consisting of one particle in potential field ( ) U q . The Hamilton function of this system is where p is the momentum of particle, q is its coordinate. Everywhere in this paper, we will assume that the potential field ( ) U q is analytical function. We assume that the system is in thermodynamic equilibrium with a thermostat. In other words, we consider a canonical ensemble with temperature T and volume V (V is one-dimensional).
The quantum canonical ensemble can be fully characterized by the Wigner function

Path Integral Representation
In general case, the density matrix in (2) cannot be calculated directly, since the kinetic and potential energy operators in Hamiltonian are non-commutative, so the statistical operator  . Therefore, we use the following procedure, leading to representation of the density matrix in form of path integrals (Winer path integral form).
Firstly, the statistical operator should be represented as product of 2K operators  , assuming that K is large number. Secondly, one should insert 2 1 K − unit operators 1 and replace each of them with completeness relation for states with a certain coordinate q : For large values of K, "high-temperature" statistical operators corresponding matrix elements can be easily calculated, using the completeness relations for states with a certain momentum p and the wave functions i q p pq =  . Finally, we obtain the representation of the Wigner function in form of 2K-dimensional integral over variables k q , and Journal of Applied Mathematics and Physics ξ : Here we assume that 2 is small and tends to zero for K → ∞ ; while the formula becomes exact.

Basics of the Complexification on Lefschetz Thimbles
According to the Morse theory [3] [4] [5] the regions of integration over each real k q in the path integrals (5) The ( ) which converge to , k q τ at l → −∞ , so that its intersection number with ( ) k q l σ ∈ ϒ is unity and vanishing otherwise, Strictly speaking this means that which holds in the homological sense. As consequence, for the critical points σ q satisfying

Numerical Algorithm
To do simulations we are going to combine the Monte Carlo method (MC) [16] [17] for searching the critical points q σ and the finite-difference methods for solving Equation (6) which can be rewritten as To increase efficiency of the numerical procedure we are going to separate the transition in two sub-steps: the proposal and the acceptance-rejection. The transition probability can be written as the product: ( ) ( ) ( ) The proposal distribution ( ) g ′ → q q is the conditional probability of proposing a state ′ q for given q . The acceptance distribution ( ) A ′ → q q is the conditional probability to accept the proposed state ′ q .
Inserting this relation in the previous equation, we have Then it is necessary to choose an acceptance that fulfills detailed balance. One This means that we always accept when the acceptance is bigger than 1 and we can accept or reject when the acceptance is smaller than 1.
It is important to notice that it is not clear, in a general problem, which distribution ( ) g ′ → q q one should use. It is a free parameter of the method which has to be adjusted to the particular problem "in hand". The probability

Results of Numerical Test Calculations
Before calculations of multidimensional integrals (5) we have to test the algorithm proposed above by calculations of some simpler integrals with known answer. It is reasonable to begin with consideration of low dimensional improper integrals of strongly oscillating functions.

The Airy Function
Let us consider a classic example: the Airy function is defined as the integral over the real axis C  and can be considered as a Fourier transform: The integrand is strongly oscillating function on C  , which makes a direct numerical evaluation infeasible. The left plot of Figure 1 shows lines of the constant imaginary part of the power in exponent in (14)  Results of the MC calculation at 2 i4 p = + are presented in Table 1

Short Time Wigner Path Integral
Now to test the developed approach let us consider elementary factor in the finite dimensional approximation of the path integral, which may be rewritten in the form like (see (5)): where, for example, 2 const = .
The left and right plots of Figure 3 are the contour plots of the image and real parts of the power of exponent in integrand in (15) respectively. Right plot of As before solution of the complex valued differential Equations (6) and (7) with MC initial conditions nearby both upper critical points allow to obtain averaged downward (red lines) and upward (blue lines) flows (see Figure 4). As

Short Time Feynman Path Integral
As the second example, we consider an elementary factor in discrete form of path integral (5) with imaginary parameter β : ( ) ( ) 2 const = . Two critical points and its Lefschetz thimble should be taken into account according to the explained above procedure (see red and blue lines on the right plot of Figure 5). In this case analytical estimations of the integral (16) are not available, while the independent numerical estimations due to the fast oscillations of the integrand in (16) is not reliable. The Monte Carlo calculations of the (16) demonstrate the fast convergence and will be considered later.

Conclusion
The main goal of this paper is to develop a new effective Monte Carlo method for numerical evaluation of a Feynman path integrals suffering to the "sign problem". This approach combines the basic ideas of the Metropolis and Hasting algorithms and is based on the Picard-Lefschetz theory and complex-valued version of Morse theory. Developed approach allows also simulating the path integral representation of the Wigner function. The basic ideas from mathematics come from Picard-Lefschetz theory and from Morse theory based on selection of the manifolds approaching the saddle points of the integrand, where the imaginary part of the complex-valued action stays constant (Lefschetz thimbles). Since the imaginary part of the action is constant on each thimble, the "sign problem" on the thimble disappears and the Feynman integrals can be calculated much more effectively. Some simple 1D test calculations and comparisons with available analytical results have been carried out. We hope that this method can also provide a new perspective in the path integration.