Lefschetz Thimbles and Wigner Functions

Vladimir S. Filinov, Alexander S. Larkin

Joint Institute for High Temperatures, Russian Academy of Sciences, Moscow, Russia.

DOI: 10.4236/jamp.2020.87098                

Abstract

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 techniques being able to solve this problem. The new method 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. 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 complex-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.

Keywords

Wigner Function, Path Integral, Quantum Monte Carlo, Picard-Lefschetz Theory, Morse Theory

1. 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 increase 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 the Wigner formulation of quantum mechanics. Section 3 is devoted to the brief deriving the path integral representation of the Wigner function. In Section 4 we explain the basic ideas of the complexification of the variables in path integrals and introduce in complex-valued space of the path integration the manyfold approaching the saddle point, where the imaginary part of the path integral action is constant (Lefschetz thimble). Path integration on the thimble allows reducing the “sign problem”. In Section 5 we derive the Metropolis-Hasting algorithm allowing finding the saddle points in complex-valued space of integration. Vicinities of the saddle pints are used as initial conditions to generate the Lefschetz thimbles. Section 6 presents results of the test calculations for 1D case with known analytical answer. Conclusion of the work is given in Section 7.

2. Wigner Approach to Quantum Mechanics

We consider a one-dimensional quantum-mechanical system consisting of one particle in potential field $U\left(q\right)$. The Hamilton function of this system is

$H\left(p\mathrm{,}q\right)=\frac{p^2}{2m}+U\left(q\right)\mathrm{,}$ (1)

where p is the momentum of particle, q is its coordinate. Everywhere in this paper, we will assume that the potential field $U\left(q\right)$ 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 $W\left(p\mathrm{,}q\right)$, which is essentially a density matrix in the mixed coordinate-momentum representation and defined as a Fourier transform of the density matrix:

$W\left(p\mathrm{,}q\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}\text{d}\xi {\text{e}}^{\frac{\text{i}}{\hslash }p\xi }\langle q-\xi /2|{\text{e}}^{-\beta \hat{H}}|q+\xi /2\rangle \mathrm{,}$ (2)

where $\beta =1/\left(kT\right)$, $\hat{H}$ is Hamiltonian, obtained from (1) by replacing $p,q$ by momentum and coordinate operators $\hat{p},\hat{q}$ respectively. The Wigner function $W\left(p\mathrm{,}q\right)$ can be formally considered as quantum generalization of classical distribution function in the phase space. This interpretation is noncompletely correct, because $W\left(p\mathrm{,}q\right)$ may be non-positive. However momentum and coordinate distribution functions, as well as average values of physical quantities depending on p and q, can be obtained from the Wigner function in classical-like way:

$\begin{array}{l}F\left(p\right)={\displaystyle \underset{V}{\int }}\text{d}qW\left(p\mathrm{,}q\right)\mathrm{,}\\ W\left(q\right)={\displaystyle \int }\frac{\text{d}p}{2\pi \hslash }W\left(p\mathrm{,}q\right)\mathrm{,}\\ \langle \hat{A}\rangle ={\displaystyle \int }\frac{\text{d}p\text{d}q}{2\pi \hslash }W\left(p\mathrm{,}q\right)A\left(p\mathrm{,}q\right)\mathrm{,}\end{array}$ (3)

where $A\left(p\mathrm{,}q\right)$ is so-called Weyl symbol of operator $\hat{A}$.

3. 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 $exp\left(-\beta \hat{H}\right)$ does not split into $exp\left[-\beta {\hat{p}}^2/\left(2m\right)\right]$ and $exp\left[-\beta U\left(\hat{q}\right)\right]$. 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 $exp\left(-\beta \hat{H}\right)$ should be represented as product of 2K operators $exp\left[-\beta /\left(2K\right)\hat{H}\right]$, assuming that K is large number. Secondly, one should insert $2K-1$ unit operators $\hat{1}$ and replace each of them with completeness relation for states with a certain coordinate $|q\rangle$ :

$\begin{array}{l}\langle q-\xi /2|{\text{e}}^{-\beta \hat{H}}|q+\xi /2\rangle \\ =\langle q-\xi /2|{\text{e}}^{-\frac{\beta }{2K}\hat{H}}\hat{1}\cdots \hat{1}{\text{e}}^{-\frac{\beta }{2K}\hat{H}}|q+\xi /2\rangle \\ ={\displaystyle \int }\cdots {\displaystyle \int }\text{d}{q}_{-K+1}\cdots \text{d}{q}_{K-1}\left[{\displaystyle \underset{k=-K}{\overset{K-1}{\prod }}}\langle q_{k+1}|{\text{e}}^{-\frac{\beta }{2K}\hat{H}}|q_k\rangle \right]\mathrm{.}\end{array}$ (4)

For large values of K, “high-temperature” statistical operators $exp\left(-\beta /\left(2K\right)\hat{H}\right)$ can be decomposed into the product of the operators $exp\left[-\beta /\left(2K\right){\hat{p}}^2/\left(2m\right)\right]$ and $exp\left[-1/2\beta /\left(2K\right)U\left(\hat{q}\right)\right]$ with accuracy $O\left[{\left(\beta /\left(2K\right)\right)}^2\right]$. After this step,

corresponding matrix elements can be easily calculated, using the completeness relations for states with a certain momentum $|p\rangle$ and the wave functions $\langle q|p\rangle =\text{i}pq/\hslash$. Finally, we obtain the representation of the Wigner function in form of 2K-dimensional integral over variables $q_k$, $k=\mathrm{0,}\pm \mathrm{1,}\cdots \mathrm{,}\pm \left(K-1\right)$ and $\xi$ :

$W\left(p\mathrm{,}q\right)={\left(\frac{2\pi \hslash \Delta \tau }{m}\right)}^{-K}{\displaystyle \int }\text{d}\xi {\displaystyle \int }\cdots {\displaystyle \int }\text{d}{q}_{-K+1}\cdots \text{d}{q}_{K-1}{\text{e}}^{-{\Phi }_{2K+1}\left(q_{-K}\mathrm{,}\cdots \mathrm{,}{q}_K\right)}+O\left(\Delta \tau \right)\mathrm{,}$ (5)

${\Phi }_{2K+1}\left(q_{-K}\mathrm{,}\cdots \mathrm{,}{q}_K\right)=-\frac{\text{i}}{\hslash }p\xi +\frac{\Delta \tau }{\hslash }{\displaystyle \underset{k=-K}{\overset{K-1}{\sum }}}\left[\frac{m}{2}{\left(\frac{q_{k+1}-q_k}{\Delta \tau }\right)}^2+\frac{U\left(q_{k+1}\right)+U\left(q_k\right)}{2}\right]\mathrm{.}$

Here we assume that $q_{\pm K}=q\mp \xi /2$. Parameter $\Delta \tau =\beta \hslash /\left(2K\right)$ is small and tends to zero for $K\to \infty$ ; while the formula becomes exact.

4. Basics of the Complexification on Lefschetz Thimbles

According to the Morse theory [3] [4] [5] the regions of integration over each real $q_k$ in the path integrals (5) is equivalent to the integration over complex valued set of Lefschetz thimbles which is homologically equivalent due to the Cauchy’s integral theorem to the integration on real cycle $C_{ℝ}$. Assuming that $q_k$ takes the complex values $q_k\in ℂ$ and the action $\Phi \left[q_k\left(\tau \right)\right]$ is extended to a holomorphic function of q let us consider the set $\Sigma$ of critical points (saddle

points) $q_{k\mathrm{,}\sigma }$, which satisfy condition ${\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}|}_{\bar{q}={\bar{q}}_{k,\sigma }}=0$. The real Morse function in our case can be defined as $h\equiv -\Re \left\{\Phi \left[q_k\right]\right\}$ and the associate gradient (downward) flow equations are given by [2] [12]:

$\frac{\text{d}q}{\text{d}l}=\overline{\frac{\partial \Phi \left(q\right)}{\partial q}}\mathrm{,}l\in ℝ\mathrm{.}$ (6)

The Morse function h is always strictly decreasing along a flow. Associated with a critical point $q_{k\mathrm{,}\sigma }$, a Lefschetz thimble ${\Upsilon }_{\sigma }$ [12] is defined by the union of all downward flows, which trace back to $q_{k\mathrm{,}\sigma }$ at $l\to -\infty$. Let us note that if $q_k\left(l\right)$ equals a critical point at some l, then the flow equation implies that $q_k\left(l\right)$ is constant for all l. So a non-constant flow can only reach a critical point at $l\to -\infty$.

One can also introduce another submanifold ${\Pi }_{\tau }$ of by the union of all upward flows, satisfying equation with opposite sign to Equation (6) [12] [13] [14]

$\frac{\text{d}}{\text{d}l}{\tilde{q}}_k\left(l\right)=-\frac{\partial \bar{\Phi }\left[{\overline{\stackrel{˜}{q}}}_k\right]}{\partial {\overline{\stackrel{˜}{q}}}_k}\mathrm{,}l\in {ℝ}^1$ (7)

which converge to $q_{k\mathrm{,}\tau }$ at $l\to -\infty$, so that its intersection number with $q_k\left(l\right)\in {\Upsilon }_{\sigma }$ is unity and vanishing otherwise, $\langle {\Upsilon }_{\sigma }\mathrm{,}{\Pi }_{\tau }\rangle ={\delta }_{\sigma \mathrm{,}\tau }$.

Strictly speaking this means that $q_k\left(l\right)\in {\Upsilon }_{\sigma }$ if $q_k\left(l\right)$ is the solution of Equation (6) and for any positive (maybe very small) $ϵ>0$ there exists the positive (maybe very large) $L>0$ that for all negative l such that $-l>L$ we have $\Vert q_k\left(l\right)-q_{k\mathrm{,}\sigma }\Vert <ϵ$. This allows in numerical simulation to use the following approximation. Due to the restrictions on computational time at solving Equation (6), (7) it is necessary to exclude the small $ϵ$ vicinity of critical point, where the parameter $l\to -\infty$. So in numerical simulations the staring points for downward flow have to be chosen outside of a small vicinity of $q_{k\mathrm{,}\sigma }$ and the averaging results of calculations by the Monte Carlo method can be done over ensemble of the downward flows related to decreasing small $ϵ$. Algorithm of this MC approach, sampling the main contribution in integrals like (5) will be discussed below.

Then, according to Morse theory [2] [12] [15], it follows that

$C_{ℝ}={\displaystyle \underset{\sigma \in \Sigma }{\sum }}{n}_{\sigma }{\Upsilon }_{\sigma }\mathrm{,}{n}_{\sigma }=\langle C_{ℝ}|{\Pi }_{\tau }\rangle$ (8)

which holds in the homological sense. As consequence, for the critical points $q_{\sigma }$ satisfying $\Re \left[-\Phi \left(q_{\sigma }\right)\right]\ge max\Re \left(-\Phi \left(q\right)\right)$, $q\in C_{ℝ}$, it holds that $\langle C_{ℝ}|{\Pi }_{\sigma }\rangle =0$ and the associated thimbles do not contribute to the path integration. On the other hand, it holds that if $\langle C_{ℝ},{\Pi }_{\sigma }\rangle =1$ the associated thimbles contribute with the relative weights proportional to ${\text{e}}^{-\Re \left(\Phi \left[q\left(\tau \right)\right]\right)}$. Now, for example, the momentum distribution function Equation (3) can be given by the formula

$F\left(p\right)={\displaystyle \underset{\sigma \in \Sigma }{\sum }}{n}_{\sigma }{\text{e}}^{-\text{i}\Im \left(\Phi \left[q\left(\tau \right)\right]\right)}{\displaystyle {\int }_{{\Upsilon }_{\sigma }}}Dq\left(\tau \right){\text{e}}^{-\Re \left(\Phi \left[q\left(\tau \right)\right]\right)}\mathrm{.}$ (9)

5. Numerical Algorithm

To do simulations we are going to combine the Monte Carlo method (MC) [16] [17] for searching the critical points $q_{\sigma }$ and the finite-difference methods for solving Equation (6) with initial conditions ${\bar{q}}_{\sigma }$ obtained by MC method in the small $ϵ$ vicinity of $q_{\sigma }$.

Used here MC method is based on the Metropolis-Hastings algorithm [16] [17], which resides in designing a Markov process (by constructing transition probabilities $P\left(\bar{q}\to {\bar{q}}^{\prime }\right)$ ), such that its stationary distribution to be equal to $w\left(\bar{q}\right)$. The derivation of the algorithm starts with the condition of detailed balance:

$w\left(\bar{q}\right)P\left(\bar{q}\to {\bar{q}}^{\prime }\right)=w\left({\bar{q}}^{\prime }\right)P\left({\bar{q}}^{\prime }\to \bar{q}\right)$ (10)

which can be rewritten as

$\frac{P\left(\bar{q}\to {\bar{q}}^{\prime }\right)}{P\left({\bar{q}}^{\prime }\to \bar{q}\right)}=\frac{w\left({\bar{q}}^{\prime }\right)}{w\left(\bar{q}\right)}\mathrm{.}$ (11)

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: $P\left(\bar{q}\to {\bar{q}}^{\prime }\right)=g\left(\bar{q}\to {\bar{q}}^{\prime }\right)A\left(\bar{q}\to {\bar{q}}^{\prime }\right)$. The proposal distribution $g\left({\bar{q}}^{\prime }\to \bar{q}\right)$ is the conditional probability of proposing a state ${\bar{q}}^{\prime }$ for given $\bar{q}$. The acceptance distribution $A\left(\bar{q}\to {\bar{q}}^{\prime }\right)$ is the conditional probability to accept the proposed state ${\bar{q}}^{\prime }$.

Inserting this relation in the previous equation, we have

$\frac{A\left(\bar{q}\to {\bar{q}}^{\prime }\right)}{A\left({\bar{q}}^{\prime }\to \bar{q}\right)}=\frac{w\left({\bar{q}}^{\prime }\right)}{w\left(\bar{q}\right)}\frac{g\left({\bar{q}}^{\prime }\to \bar{q}\right)}{g\left(\bar{q}\to {\bar{q}}^{\prime }\right)}\mathrm{.}$ (12)

Then it is necessary to choose an acceptance that fulfills detailed balance. One common choice is the Metropolis’s suggestion:

$A\left(\bar{q}\to {\bar{q}}^{\prime }\right)=min\left(\mathrm{1,}\frac{w\left({\bar{q}}^{\prime }\right)}{w\left(\bar{q}\right)}\frac{g\left({\bar{q}}^{\prime }\to \bar{q}\right)}{g\left(\bar{q}\to {\bar{q}}^{\prime }\right)}\right)\mathrm{.}$ (13)

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\left(\bar{q}\to {\bar{q}}^{\prime }\right)$ one should use. It is a free parameter of the method which has to be adjusted to the particular problem “in hand”. The probability $g\left(\bar{q}\to {\bar{q}}^{\prime }\right)$ affects only the efficiency of sampling the main contribution to the integrals and does not change the final result of calculations. For optimization of the MC finding the main contribution to the integral (9) the choice of the $g\left(\bar{q}\to {\bar{q}}^{\prime }\right)$ may be the following $g\left(\bar{q}\to {\bar{q}}^{\prime }\right)={\text{e}}^{-\beta \Re \left(\Phi \left[{\bar{q}}^{\prime }\right]\right)}/{\text{e}}^{-\beta \Re \left(\Phi \left[\bar{q}\right]\right)}$ with free appropriate fit parameter $\beta$. To optimize the MC search of the critical points ${\bar{q}}_{k\mathrm{,}\sigma }$ ( ${\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}|}_{\bar{q}={\bar{q}}_{k,\sigma }}=0$ ) we define the probability $w\left(\bar{q}\right)$ as

$w\left(\bar{q}\right)=exp\left(-b{\left|\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}\right|}^2\right)$ with parameter $b\ge 1$. The ideal acceptance rate,

which is the fraction of proposed samples that is accepted during the last N samples, have to be in interval of 23% - 50%.

The Metropolis-Hastings algorithm consists in the following steps:

1) Initialization: pick an initial state point $\bar{q}$ at random.

2) Randomly pick a state ${\bar{q}}^{\prime }$, according to probability $g\left(\bar{q}\to {\bar{q}}^{\prime }\right)$.

3) Accept the state according to the probability $A\left(\bar{q}\to {\bar{q}}^{\prime }\right)$. If not accepted, that means that ${\bar{q}}^{\prime }=\bar{q}$, and so there is no need to update anything. Else, the system transits to ${\bar{q}}^{\prime }$.

4) Go to 2 until many M states were generated to “forget” initial $\bar{q}$ and to obtain the average position of the critical point $\langle \bar{q}\rangle$.

5) If $w\left(\langle \bar{q}\rangle \right)\ge 0.95$ carry out several iterations by complex-valued extension of the Newton’s method to produce better approximations to the roots (or zeroes) of the complex-valued function ${\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}|}_{\bar{q}={\bar{q}}_{\sigma }}=0$. Else go to 2.

6) Save the state ${\bar{q}}_{\sigma }=\langle \bar{q}\rangle$.

7) At random pick an initial point $\delta \bar{q}$ at the small vicinity of the zero point of the complex-valued space (maybe in any given quarter).

8) Analogously randomly pick a new state $\delta \overline{q^{\prime }}$ according to the probability $g\left({\bar{q}}_{\sigma }+\delta q\to {\bar{q}}_{\sigma }+\delta q^{\prime }\right)$.

9) Go to 8 many times $M^{\prime }$ to “forget” initial $\delta \bar{q}$ state and to obtain the average value of the $\langle \delta \overline{q^{\prime }}\rangle$.

10) Solve equation for upward flow ${\Pi }_{\tau }$ with initial conditions of ${\bar{q}}_{\sigma }+\langle \delta q^{\prime }\rangle$ to check if $n_{\sigma }=\langle C_{ℝ}\mathrm{,}{\Pi }_{\sigma }\rangle =1$. If $n_{\sigma }=1$ go to 11, otherwise go to 2.

11) Solve Equation (6) with initial conditions of ${\bar{q}}_{\sigma }+\langle \delta q^{\prime }\rangle$ until modulus of the integrand in (9) will be smaller of several order of magnitude of its initial value and calculate the integral sum related to (9).

12) Go to 2 many times and calculate the integrals (9).

Function $F\left(p\right)$ (9) has be averaged over the set of downward flows-solutions of the Equation (6) with initial conditions $q$ provided by MC procedure from neighborhood of $q_{\sigma }$.

6. 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.

6.1. 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:

$\text{Ai}\left(p\right)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}\text{d}xexp\left[\text{i}\left(px+\frac{x^3}{3}\right)\right]\mathrm{,}$ (14)

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) on complex plane, while increasing and lowering values are denoted by the red and blue regions. The critical points are in the left upper and the right bottom quarters of the complex plane. We can deform the integration path in the complex plane of variable $z=x+\text{i}y$, as long as the new path belongs to the original relative homology class, which connects regions of strongly decaying modulus of the integrand at infinity (blue regions on the centre plot of Figure 1). Right panel of this figure shows the contour plot of the MC probability $w\left(\bar{q}\right)=exp\left(-b{\left|\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}\right|}^2\right)$ with two red circles at its maximum values at the critical points.

Testing the proposed approach starts from finding critical points by suggested MC method. It turns out that the Markovian chain traveling on the whole complex plane always stabilizes in the vicinity of the critical point $z_{\sigma }=-1.11+\text{i}1.79$ (point in the left upper quarter of the complex plane in Figure 1) ignoring the second critical point ( $z_{\sigma }=+1.11-\text{i}1.79$ ) in the favoure of the first one. So to

Figure 1. (Color online) The contour plot of the imaginary part (left panel) and the real part (central panel) of the power in exponent in (14) on complex plane for $p=2+4\text{i}$. (Right plot) Contour plot of the probability $w\left(\bar{q}\right)=exp\left(-b{\left|\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}\right|}^2\right)$.

force the Markovian chain to stabilize in the vicinity of the second critical point we have to use the special restrictions. Reason of this behavior of the Markovian chain is the asymmetry of the contour plots of the real and imaginary parts of the power in exponent in (14) (see Figure 1).

Solution of the complex valued differential Equations (6 and 7) with MC initial conditions nearby both critical points allow to obtain averaged downward (red lines) and upward (blue lines) flows (see both plots of Figure 2). Right plot of this figure shows in detail the downward and upward flows from different quarters of the small enough vicinities of the critical points. Let us note that initial conditions for red downward and blue upward flows were the same to test their fast converge to the related ${\Upsilon }_{\sigma }$ and ${\Pi }_{\sigma }$ respectively. Let us note that the power low grows on the complex plane of the right part of differential Equations (6) and (7) results in limitations on the “time” l of obtained solutions at needed given accuracy.

As only the blue line for red point crosses the real axis ( $\langle C_{ℝ}|{\Pi }_{\sigma }\rangle =1$ ) we calculate integral defined the Airy function along the red Lefschetz thimble. As we mentioned above the Markovian chain traveling on the whole complex plane prefer namely this critical point ignoring the second one. The reason of this interesting fact has been further investigated.

Results of the MC calculation at $p=2+\text{i}4$ are presented in Table 1 as well as results of some additional calculations for $p=0+\text{i}4$ and $p=4+\text{i}0$. Comparison of obtained MC results with well known values of Airy function demonstrates a good enough agreement. Discrepancy between exact values of the Airy function and related values obtained by proposed procedure can be explained by approximations used in transitions from initial integral to its Lefschetz thimbles representation accounting for only the main contribution to the contour integrals.

Figure 2. (Color online) (Left plot) The averaged downward flows are lines-1 ( $\in {\Upsilon }_{\sigma }$ ) and the upward flows are lines-2 ( $\in {\Pi }_{\sigma }$ ). The critical points are points 3 and 4 ( $\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}=0$ ) of the integrand in (14). Red critical point and the associated thimbles contribute to the contour integral of the Airy function (14) as $\langle C_{ℝ},{\Pi }_{\sigma }\rangle =1$, while the blue point not, as $\langle C_{ℝ},{\Pi }_{\sigma }\rangle =0$. (Right plot) Details of initial conditions obtained by MC method in different quarters of small vicinity of the critical point.

Table 1. The MC results versus the exact Airy function.

6.2. 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)):

$I_k\left(p_k\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}\text{d}{q}_{k}exp\left\{\text{i}\left(p_k{q}_k+\text{i}\left({\left(q_k-const\right)}^2+q_k^4/4\right)\right)\right\}\mathrm{,}$ (15)

where, for example, $const=2$.

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 Figure 3 presents contour plot of the probability, $w\left(\bar{q}\right)=exp\left(-b{\left|\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}\right|}^2\right)$ which with conditional probability $g\left(\bar{q}\to {\bar{q}}^{\prime }\right)$ allows finding by MC method the critical points and provides the initial conditions for solution of the Equations (6) and (7). According to the Monte Carlo simulation the critical points are $z_1\approx 0.5712+\text{i}0.8786$, $z_2\approx -0.5712+\text{i}0.8786$ and $z_3\approx -1.77+\text{i}0.0$, which agree with the regions of the saddle-like behaviour of the contour lines on left and right plots of Figure 3. Here the Markovian chain traveling on the whole complex plane always stabilizes in the vicinity of the critical point $z_1$ and $z_2$ ignoring the point $z_3$.

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 both blue lines for red point cross the real axis ( $\langle C_{ℝ}|{\Pi }_{\sigma }\rangle =1$ ) we calculate sum of the contributions to the integral (15) along the both red Lefschetz thimbles with opposite sign due to the different thimble orientation. According to the Morse theory contribution of the bottom critical point has to be ignored as the related blue line does not cross real axis (here not shown) and the maximum of the integrand on the real axis is smaller than the value of the integrand at this critical point.

As before details the MC initial conditions and related downward and upward flows are presented by right plot of Figure 4. Let us note that red and blue lines start at the same MC initial points. Reason of asymmetry in behavior of the red and blue lines is in asymmetry of the real and imaginary parts of the power in exponent (see Figure 4). Let us stress the fast convergence of the upward and downward flows to the related limit lines ${\Upsilon }_{\sigma }$ and ${\Pi }_{\sigma }$. Comparison of the Monte Carlo result for $I_k=0.01306+\text{i}0.00006$ with independent exact calculations

Figure 3. (Color online) The contour plot of imaginary part(left panel) and real part (central panel) of the of the power in exponent in (15) on complex plane for $const=2$. (Right plot) Contour plot of the probability $w\left(\bar{q}\right)=exp\left(-b{\left|\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}\right|}^2\right)$.

Figure 4. (Color online) The averaged downward flows are lines-1 ( $\in {\Upsilon }_{\sigma }$ ) and the upward flows are lines-2 ( $\in {\Pi }_{\sigma }$ ) at $const=2$ and $const=2$ and $p=2+\text{i}4$. The critical points of the integrand in (15) are points 3 ( $\partial \bar{\Phi }\left[\bar{q}\right]/\partial \bar{q}=0$ ). Red critical point and the associated thimbles contribute to the contour integral (15) as $\langle C_{ℝ},{\Pi }_{\sigma }\rangle =1$. (Right plot) Details of MC initial conditions in different quarters of small vicinity of the critical points.

$I_k=0.01402+\text{i}0$ demonstrate a good enough accuracy of the developed approach.

6.3. Short Time Feynman Path Integral

As the second example, we consider an elementary factor in discrete form of path integral (5) with imaginary parameter $\beta$ :

$I\left(p\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}\text{d}xexp\left\{\text{i}\left(px+{\left(x-const\right)}^2+\frac{x^3}{3}\right)\right\}\mathrm{.}$ (16)

Contour plots on Figure 5 show imaginary and real parts of the complex-valued action $\Phi \left(z\right)=\text{i}\left(pz+{\left(z-const\right)}^2+\frac{z^3}{3}\right)$, when $p=2+4\text{i}$ and

$const=2$. 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.

Figure 5. The contour plot of phase (left panel) and modulus (central panel) of the integrand from (16) on complex-valued plane for $p=2+4\text{i}$ and $const=2$. Averaged downward and upward flows from the neighborhood of the two critical points (red and blue lines on the right plot correspondingly).

7. 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.

Acknowledgements

We acknowledge stimulating discussions with D. Blaschke and Yu.B. Ivanov.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Fedoryuk, M.V. (1976) The Asymptotics of the Fourier Transform of the Exponential Function of a Polynomial. Doklady Akademii Nauk, 227, 580-583.
[2]	Fedoryuk, M. (1989) Asymptotic Methods in Analysis. In: Gamkrelidze, R.V., Ed., Analysis I, Springer, Berlin, 83-191. https://doi.org/10.1007/978-3-642-61310-4_2
[3]	Witten, E. (2011) Analytic Continuation of Chern-Simons Theory. AMS/IP Studies in Advanced Mathematics, 50, 347.
[4]	Alexandru A., Başar, G. and Bedaque, P. (2016) Monte Carlo Algorithm for Simulating Fermions on Lefschetz Thimbles. Physical Review D, 93, Article ID: 014504. https://doi.org/10.1103/PhysRevD.93.014504
[5]	Alexandru, A., Başar, G., Bedaque, P.F. and Ridgway, G. (2017) Schwinger-Keldysh Formalism on the Lattice: A Faster Algorithm and Its Application to Field Theory. Physical Review D, 95, Article ID: 114501. https://doi.org/10.1103/PhysRevD.95.114501
[6]	Parisi, G., Klauder, J., Petersen, W., Ambjorn, J., Yang, S., Karsch, F. and Wyld, H. (1988) On Complex Probabilities. Stochastic Quantization, 131, 381.
[7]	Klauder, J.R. (1983) A Langevin Approach to Fermion and Quantum Spin Correlation Functions. Journal of Physics A: Mathematical and General, 16, L317. https://doi.org/10.1088/0305-4470/16/10/001
[8]	Klauder, J.R. (1984) Coherent-State Langevin Equations for Canonical Quantum Systems with Applications to the Quantized Hall Effect. Physical Review A, 29, 2036. https://doi.org/10.1103/PhysRevA.29.2036
[9]	Fujii, H., Honda, D., Kato, M., Kikukawa, Y., Komatsu, S. and Sano, T. (2013) Hybrid Monte Carlo on Lefschetz Thimbles—A Study of the Residual Sign Problem. Journal of High Energy Physics, 2013, Article No. 147. https://doi.org/10.1007/JHEP10(2013)147
[10]	Mukherjee, A., Cristoforetti, M. and Scorzato, L. (2013) Metropolis Monte Carlo Integration on the Lefschetz Thimble: Application to a One-Plaquette Model. Physical Review D, 88, Article ID: 051502. https://doi.org/10.1103/PhysRevD.88.051502
[11]	Di Renzo, F. and Eruzzi, G. (2015) Thimble Regularization at Work: From Toy Models to Chiral Random Matrix Theories. Physical Review D, 92, Article ID: 085030. https://doi.org/10.1103/PhysRevD.92.085030
[12]	Cristoforetti, M., Di Renzo, F., Scorzato, L., Collaboration, A., et al. (2012) New Approach to the Sign Problem in Quantum Field Theories: High Density QCD on a Lefschetz Thimble. Physical Review D, 86, Article ID: 074506. https://doi.org/10.1103/PhysRevD.86.074506
[13]	Weinberg, S. (1995) The Quantum Theory of Fields. Vol. 1: Foundations. Cambridge University Press, Cambridge.
[14]	Mou, Z.-G., Saffin, P.M., Tranberg, A. and Woodward, S. (2019) Real-Time Quantum Dynamics, Path Integrals and the Method of Thimbles. Journal of High Energy Physics, 2019, Article No. 94. https://doi.org/10.1007/JHEP06(2019)094
[15]	Berges, J., Borsanyi, S., Sexty, D. and Stamatescu, I.-O. (2007) Lattice Simulations of Real-Time Quantum Fields. Physical Review D, 75, Article ID: 045007. https://doi.org/10.1103/PhysRevD.75.045007
[16]	Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. and Teller, E. (1953) Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21, 1087. https://doi.org/10.1063/1.1699114
[17]	Hastings, W.K. (1970) Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57, 97-109. https://doi.org/10.1093/biomet/57.1.97