Parallel tempering algorithm for integration over Lefschetz thimbles

TL;DR

This paper introduces a parallel tempering algorithm using flow time as an auxiliary variable to improve Monte Carlo sampling over Lefschetz thimbles, effectively addressing the sign problem in complex action integrations.

Contribution

The paper proposes a novel parallel tempering method employing flow time as a tempering parameter to enhance exploration across Lefschetz thimbles in Monte Carlo simulations.

Findings

Successfully applied to the (0+1)-dimensional Thirring model at finite density.

Reproduces analytic results for large flow times, confirming effectiveness.

Addresses the difficulty of exploring multimodal distributions in complex actions.

Abstract

The algorithm based on integration over Lefschetz thimbles is a promising method to resolve the sign problem for complex actions. However, this algorithm often meets a difficulty in actual Monte Carlo calculations because the configuration space is not easily explored due to the infinitely high potential barriers between different thimbles. In this paper, we propose to use the flow time of the antiholomorphic gradient flow as an auxiliary variable for the highly multimodal distribution. To illustrate this, we implement the parallel tempering method by taking the flow time as a tempering parameter. In this algorithm, we can take the maximum flow time to be sufficiently large such that the sign problem disappears there, and two separate modes are connected through configurations at small flow times. To exemplify that this algorithm does work, we investigate the (0+1)-dimensional massive Thirring model at finite density and show that our algorithm correctly reproduces the analytic results for large flow times such as T=2.

Equations

$$⟨\mathcal{O}(x)⟩=\frac{{\int }_{{\mathbb{R}}^N}dx{e}^{-S(x)}\mathcal{O}(x)}{{\int }_{{\mathbb{R}}^N}dx{e}^{-S(x)}}.$$

$$\frac{d{z}^i}{dt}$$

$$z_t:{\Sigma }_0\ni x\mapsto z_t(x)\equiv z(t;x)\in {\Sigma }_t.$$

$$⟨\mathcal{O}(x)⟩=\frac{{\int }_{{\Sigma }_t}dz{e}^{-S(z)}\mathcal{O}(z)}{{\int }_{{\Sigma }_t}dz{e}^{-S(z)}},$$

$$⟨\mathcal{O}(x)⟩$$

$$\displaystyle=\frac{\Bigl{\langle}e^{i\bigl{[}{\rm arg}\det\!J_{t}(x)-S_{I}(z_{t}(x))\bigr{]}}\mathcal{O}\bigl{(}z_{t}(x)\bigr{)}\Bigr{\rangle}_{S_{\rm eff}}}{\Bigl{\langle}e^{i\bigl{[}{\rm arg}\det\!J_{t}(x)-S_{I}(z_{t}(x))\bigr{]}}\Bigr{\rangle}_{S_{\rm eff}}}.$$

$$\displaystyle S_{\rm eff}(x;t)\equiv S_{R}\bigl{(}z_{t}(x)\bigr{)}-\ln\,\bigl{|}\det\!J_{t}(x)\bigr{|}.$$

$$\displaystyle\frac{d(J_{t}(x))_{ij}}{dt}=\frac{d}{dt}\Bigl{(}\frac{\partial z_{t}^{i}(x)}{\partial x^{j}}\Bigr{)}=\frac{\partial}{\partial x^{j}}\Bigl{[}\frac{\partial S}{\partial z^{i}}(z_{t}(x))\Bigr{]}^{\ast}=\Bigl{[}\frac{\partial^{2}S}{\partial z^{i}\partial z^{k}}(z_{t}(x))\,\frac{\partial z^{k}_{t}(x)}{\partial x^{j}}\Bigr{]}^{\ast}=\Bigl{[}H_{ik}(z_{t}(x))\,(J_{t})_{kj}(x)\Bigr{]}^{\ast}.$$

$$\frac{d{z}_t^i}{dt}$$

$$\frac{d{(J_t)}_{ij}}{dt}$$

$$\prod _{\alpha }{e}^{-S(x_{\alpha };{\lambda }_{\alpha })}.$$

$$\displaystyle w_{\alpha}(x,x^{\prime})={\rm min}\Bigl{(}1,\,\frac{e^{-S(x^{\prime};\lambda_{\alpha})-S(x;\lambda_{\alpha+1})}}{e^{-S(x;\lambda_{\alpha})-S(x^{\prime},\lambda_{\alpha+1})}}\Bigr{)},$$

$$w_{\alpha }(x,x^{\prime })e^{-S(x;{\lambda }_{\alpha })-S(x^{\prime },{\lambda }_{\alpha +1})}=w_{\alpha }(x^{\prime },x)e^{-S(x^{\prime };{\lambda }_{\alpha })-S(x,{\lambda }_{\alpha +1})}.$$

$$\Delta S_{\mathrm{eff},\alpha }$$

$$\displaystyle=S_{R}(z^{\prime}_{t_{\alpha}})-\ln\,\bigl{|}\det\!J^{\prime}_{t_{\alpha}}\bigr{|}-S_{R}(z_{t_{\alpha}})+\ln\,\bigl{|}\det\!J_{t_{\alpha}}\bigr{|}$$

$$\displaystyle w_{\alpha}(x_{\alpha},x_{\alpha+1})\equiv{\rm min}\Bigl{(}1,\,\frac{e^{-S_{\rm eff}(x_{\alpha+1};\,t_{\alpha})-S_{\rm eff}(x_{\alpha};\,t_{\alpha+1})}}{e^{-S_{\rm eff}(x_{\alpha};\,t_{\alpha})-S_{\rm eff}(x_{\alpha+1};\,t_{\alpha+1)}}}\Bigr{)}.$$

$$\displaystyle\langle\mathcal{O}(x)\rangle\thickapprox\frac{\sum_{a}e^{i\bigl{[}{\rm arg}\det\!J_{T}^{(a)}-S_{I}(z_{T}^{(a)})\bigr{]}}\mathcal{O}\bigl{(}z_{T}^{(a)}\bigr{)}}{\sum_{a}e^{i\bigl{[}{\rm arg}\det\!J_{T}^{(a)}-S_{I}(z_{T}^{(a)})\bigr{]}}}.$$

$$Z={\int }_{\mathrm{PBC}}[d\varphi (\tau )]{\int }_{\mathrm{ABC}}[d\bar{\psi }(\tau )d\psi (\tau )]e^{-S[\varphi ,\bar{\psi },\psi ]},$$

$$\displaystyle S[\phi,\bar{\psi},\psi]=\int_{0}^{\beta}d\tau\,\Bigl{[}-\bar{\psi}\bigl{(}\gamma^{0}(\partial_{0}+i\phi+\mu)+m\bigr{)}\psi+\frac{1}{2g^{2}}\,\phi^{2}\Bigr{]},\quad\biggl{[}\gamma^{0}=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}\biggr{]},$$

$$Z=\int (dU)e^{-S(U,\bar{\chi },\chi )},$$

$$\displaystyle S(U,\bar{\chi},\chi)=-\sum_{n,\ell=1}^{N}\bar{\chi}_{n}\,D_{n\ell}(U)\,\chi_{\ell}+\frac{1}{2g^{2}a}\,\sum_{n=1}^{N}\Bigl{[}1-\frac{1}{2}(U_{n}+U_{n}^{-1})\Bigr{]}$$

$$Z$$

$$\displaystyle Z=\frac{e^{-N\alpha}}{2^{N-1}}\bigl{[}\cosh(N\mu)\,I_{1}^{N}(\alpha)+\rho_{+}\,I_{0}^{N}(\alpha)\bigr{]},$$

$$⟨\bar{\chi }\chi ⟩$$

$$\displaystyle S(U)\equiv\frac{1}{2g^{2}}\,\sum_{n}\Bigl{[}1-\frac{1}{2}\,(U_{n}+U_{n}^{-1})\Bigr{]}-\ln\det D(U)$$

$$\displaystyle\langle\bar{\chi}\chi\rangle=\Bigl{\langle}\frac{1}{N}\,{\rm tr\,}\Bigl{(}D^{-1}(U)\frac{\partial D(U)}{\partial m}\Bigr{)}\Bigr{\rangle}=\Bigl{\langle}\frac{1}{N}\,{\rm tr\,}D^{-1}(U)\Bigr{\rangle},$$

$$⟨\mathcal{O}(U)⟩\equiv \frac{\int (dU)e^{-S(U)}\mathcal{O}(U)}{\int (dU)e^{-S(U)}}.$$

Full text

1. 1 Introduction

2. 2 Algorithm

3. 2.1 Integration over Lefschetz thimbles (review)

4. 2.2 Implementation of parallel tempering

5. 3 Example

6. 4 Conclusion and outlook

1 Introduction

There are many cases where one needs to deal with complex actions. One important example for high energy/nuclear physics is quantum chromodynamics (QCD) at finite density (see Ref. [1] for a review on the recent developments). However, since a complex action does not give a real and positive Boltzmann weight, one cannot directly resort to the traditional Markov chain Monte Carlo methods to estimate correlation functions. The so-called reweighting algorithm (which absorbs the phase of the weight into observables) is highly ineffective when the imaginary part of the action becomes very large (such as in the thermodynamic limit), 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. The difficulty of numerical evaluation in such a situation is termed the sign problem.

There have been many proposals to circumvent the sign problem. One of the approaches that are currently under intense study is the use of integration over Lefschetz thimbles [2] (see also Refs. [3, 4, 5, 6, 7, 8, 9, 10]), which we will call the Lefschetz thimble algorithm hereafter. There, the original real-valued variable (say, $x=(x^{i})\in{\mathbb{R}}^{N}$) is complexified according to the antiholomorphic gradient flow (sometimes called the upward flow in the literature) $\dot{z}^{i}=\left[\partial_{i}S(z)\right]^{\ast}$. In the original Lefschetz thimble algorithm, as will be reviewed in Sect. 2.1, the flow time is taken to infinity to map the original configuration space to a union of Lefschetz thimbles. Since the imaginary part of the action is constant on each thimble, the sign problem disappears if the path integral is made over the Lefschetz thimbles. However, since two different thimbles are separated by an infinitely high potential barrier, one needs to invent some machinery to incorporate contributions from all relevant thimbles.

Recently, Alexandru et al. [11] was a very interesting proposal to consider configurations on a manifold that is obtained from the original configuration manifold by a finite amount of flow time. This algorithm was a great success in various models [11, 12], but reducing the amount of flow time may also reduce the effectiveness against the sign problem, and one does not know a priori whether the chosen flow time avoids both the sign problem and the multimodal problem simultaneously.

In this paper, as a versatile tool for Monte Carlo calculations of models with complex actions, we propose a Lefschetz thimble algorithm where the flow time is used as an auxiliary variable for the highly multimodal distribution. There can be various methods to realize this idea, and in this paper we implement the parallel tempering method because of its simplicity, by taking the flow time as a tempering parameter. There, we consider a set of manifolds corresponding to various flow times. Two separate modes at large flow times (where the sign problem no longer exists) are then connected by passing through configurations at small flow times (where the original sign problem exists but the multimodality is expected to be mild). Since the sample average is taken only with respect to the largest flow time, we need not worry about the sign problem at small flow times, although we take into account configurations there.

This paper is organized as follows. In Sect. 2 we first review the basics of the Lefschetz thimble algorithm based on Refs. [7, 11], and then implement the parallel tempering method in the algorithm by taking the flow time as a tempering parameter. In Sect. 3 we investigate the $(0+1)$-dimensional massive Thirring model at finite density, and show that our algorithm correctly reproduces the analytic results for large flow times such as $T=2$. Section 4 is devoted to the conclusion and outlook for future work.

2 Algorithm

2.1 Integration over Lefschetz thimbles (review)

We consider a real $N$-dimensional dynamical variable, $x=(x^{i})\in{\mathbb{R}}^{N}$, with action $S(x)$ that may take complex values for real-valued $x$. Our main concern is to evaluate the expectation values of functions of $x$:

[TABLE]

We assume that $|e^{-S(x)}|$ decreases rapidly enough in the limit $x\to\pm\infty$, and that $e^{-S(z)}$ and $\mathcal{O}(z)$ are entire functions when regarded as functions of $z=(z^{i})\in{\mathbb{C}}^{N}$. The integration region can then be changed to any other region $\Sigma$ in ${\mathbb{C}}^{N}$ as long as it is obtained as a continuous deformation of the original region with the boundary fixed at infinity. We consider as such an integration region the submanifold that is obtained from the following antiholomorphic gradient flow $z(t;x)$ with a flow time $t$:

[TABLE]

In fact, the flow defines a map from the original integration region $\Sigma_{0}\equiv{\mathbb{R}}^{N}$ to a real $N$-dimensional submanifold $\Sigma_{t}$ in ${\mathbb{C}}^{N}$:

[TABLE]

We thus see that (2.1) can be rewritten as

[TABLE]

which can be further rewritten as a reweighted integral over ${\mathbb{R}}^{N}$ [7] as111 One may also take as an integration region the tangent space to the critical point of the dominant thimble [7] if the integral can be well approximated by the integration around the critical point.

[TABLE]

Here, J_{t}(x)\equiv\bigl{(}\partial z_{t}^{i}(x)/\partial x^{j}\bigr{)} is the Jacobi matrix, S_{R}\bigl{(}z_{t}(x)\bigr{)} (S_{I}\bigl{(}z_{t}(x)\bigr{)}) is the real (imaginary) part of S\bigl{(}z_{t}(x)\bigr{)}, and $\langle\ast\rangle_{S_{\rm eff}}$ is the expectation value taken with respect to

[TABLE]

The Jacobi matrix is obtained by solving the combined differential equations with respect to $t$ [7]:222The second equation is obtained by differentiating the first equation with respect to $x^{j}\in{\mathbb{R}}$;

\displaystyle\frac{d(J_{t}(x))_{ij}}{dt}=\frac{d}{dt}\Bigl{(}\frac{\partial z_{t}^{i}(x)}{\partial x^{j}}\Bigr{)}=\frac{\partial}{\partial x^{j}}\Bigl{[}\frac{\partial S}{\partial z^{i}}(z_{t}(x))\Bigr{]}^{\ast}=\Bigl{[}\frac{\partial^{2}S}{\partial z^{i}\partial z^{k}}(z_{t}(x))\,\frac{\partial z^{k}_{t}(x)}{\partial x^{j}}\Bigr{]}^{\ast}=\Bigl{[}H_{ik}(z_{t}(x))\,(J_{t})_{kj}(x)\Bigr{]}^{\ast}.

[TABLE]

where $H_{ij}(z)\equiv\partial^{2}S(x)/\partial z^{i}\partial z^{j}$ is the Hesse matrix for the action $S(z)$.

The key point is that the right-hand side of (2.4) or (2.5) does not depend on the flow time $t$, so that one can set $t$ to an arbitrary value that is convenient for actual calculation. Note that under the flow the real part S_{R}\bigl{(}z_{t}(x)\bigr{)} does not decrease while the imaginary part S_{I}\bigl{(}z_{t}(x)\bigr{)} is kept constant, because (d/dt)S\bigl{(}z_{t}(x)\bigr{)}=\bigl{|}\partial_{z^{i}}S\bigl{(}z_{t}(x)\bigr{)}\bigr{|}^{2}\geq 0. In the original Lefschetz thimble algorithm, one takes the limit $t\to\infty$, in which $\Sigma_{t}$ approaches a union of connected components (Lefschetz thimbles), and the action has a constant imaginary part on each thimble.333Generically there is a single critical point $z_{\sigma}$ on each connected component $\mathcal{J}_{\sigma}$, and $\mathcal{J}_{\sigma}$ is obtained as the set of orbits flowing out of $z_{\sigma}$. The complementary submanifold to $\mathcal{J}_{\sigma}$ in ${\mathbb{C}}^{N}$ consists of orbits that flow into $z_{\sigma}$, and will be denoted by $\mathcal{K}_{\sigma}$. The integrations in (2.5) are dominated by points near the intersection of $K_{\sigma}$ and ${\mathbb{R}}^{N}$. In a generic situation the phase change coming from $J_{t}(x)$ is sufficiently mild, so that the Monte Carlo calculation for the expression (2.5) is free from sign problems. However, two different thimbles are also disconnected in the sense of Monte Carlo sampling because $S_{R}$ increases indefinitely near the boundary of each thimble. This multimodality of distribution makes the Monte Carlo calculation impractical, especially when contributions from more than one thimble are relevant to estimating expectation values.

A very interesting proposal made in Refs. [11] is to use a finite amount $t$, which is chosen to be large enough to avoid the sign problem but also not too large in order to enable exploration in the configuration space. However, one does not know a priori whether the adopted value of $t$ is actually free from the two obstacles (the sign and multimodal problems) simultaneously. We will show that one can solve both simultaneously if we implement the parallel tempering method in their algorithm with the flow time as a tempering parameter.

2.2 Implementation of parallel tempering

The basic idea of the parallel tempering algorithm [13, 14, 15] is the following. Suppose that we want to estimate expectation values with action $S(x;\lambda)$, where $x\in{\mathbb{R}}^{N}$ is a dynamical variable and $\lambda$ is the parameter (such as the temperature) that we want to use for a Monte Carlo calculation. The point is that even when the distribution is multimodal for the original $\lambda$ (e.g., when $\lambda$ represents a very low temperature), the multimodality can be made mild if one takes another value $\tilde{\lambda}$ (e.g., $\tilde{\lambda}$ corresponding to a very high temperature). So, if the configuration space is enlarged such that the parameter can change gradually between $\lambda$ and $\tilde{\lambda}$, two separate configurations for the original $\lambda$ will be connected by passing through configurations at parameters near $\tilde{\lambda}$. The parallel tempering algorithm enables the move of configurations among different $\lambda$ by enlarging the configuration space from ${\mathbb{R}}^{N}=\{x\}$ to the set of $A+1$ replicas, $({\mathbb{R}}^{N})^{A+1}=\{(x_{0},x_{1},\ldots,x_{A})\}$. We there assign $\lambda_{\alpha}$ to replica $\alpha$ ($\alpha=0,1,\ldots,A$), such that $\lambda_{0}=\tilde{\lambda}$ and $\lambda_{A}=\lambda$ and that $\lambda_{\alpha}$ and $\lambda_{\alpha+1}$ are sufficiently close to each other.444The computational cost required for the parallel tempering method can be roughly estimated to be proportional to $(A+1)/X$, where $X$ is the minimum of the acceptance rates for all pairs of adjacent replica (see Step 3 below). In this paper, we will set the parameters $A$ and ${\lambda_{\alpha}}$ such that the minimum acceptance rate is well above $50\%$ (see Fig. 5). We set up an irreducible, aperiodic Markov chain for the enlarged configuration space such that the probability distribution for $(x_{0},x_{1},\ldots,x_{A})$ eventually approaches the equilibrium distribution proportional to

[TABLE]

We finally take sample averages only with respect to a sample taken from $\alpha=A$. The simplest algorithm to realize this idea555 Of course, there can be many variations on this algorithm. is to swap two configurations of two adjacent replicas $\alpha$ and $\alpha+1$, (i.e., to update the configuration $(x_{\alpha}=x,\,x_{\alpha+1}=x^{\prime})$ to $(x_{\alpha}=x^{\prime},\,x_{\alpha+1}=x)$) with the probability

[TABLE]

which obviously satisfies the detailed balance condition

[TABLE]

Our proposal is to take the flow time $t$ as such a tempering parameter. The basic algorithm is then as follows.666To make discussions simple, we only take the flow time as a tempering parameter. The algorithm can be readily extended such that other parameters are included as extra tempering parameters.

• •

Step 0. Fix the maximum flow time $T$, which should be sufficiently large such that the sign problem disappears there, and pick up flow times $\{t_{\alpha}\}$ from the interval $[0,T]$ with $t_{0}=0<t_{1}<\cdots<t_{A}=T$. The values of $A$ and $t_{\alpha}$ are determined manually or adaptively to optimize the acceptance rate in Step 3 below.

• •

Step 1. Choose an initial value $x_{\alpha}\in{\mathbb{R}}^{N}$ for each replica $\alpha$, and numerically solve the differential equations (2.7) and (2.8) to obtain the triplet $(x_{\alpha},z_{t_{\alpha}},J_{t_{\alpha}})$.

• •

Step 2. For each $\alpha$, construct a Metropolis process to update the value of $x$. Explicitly, we take a value $x^{\prime}_{\alpha}$ from $x_{\alpha}$ using a symmetric proposal distribution, and recalculate the triplet $(x^{\prime}_{\alpha},z^{\prime}_{t_{\alpha}},J^{\prime}_{t_{\alpha}})$ using $x^{\prime}_{\alpha}$ as the initial value. We then update $x_{\alpha}$ to $x^{\prime}_{\alpha}$ with the probability ${\rm min}(1,e^{-\Delta S_{{\rm eff},\alpha}})$, where

[TABLE]

(recall that S_{\rm eff}(x;t)=S_{R}\bigl{(}z_{t}(x)\bigr{)}-\ln\,\bigl{|}\det\!J_{t}(x)\bigr{|}, Eq. (2.6)). We repeat the process sufficiently many times such that local equilibrium is realized for each $\alpha$.

• •

Step 3. Starting from $\alpha=0$ through $\alpha=A-1$, swap the values of $x$ between two adjacent replicas $\alpha$ and $\alpha+1$ with the probability

[TABLE]

• •

Step 4. After repeating Steps 2 and 3 sufficiently many times, get a triplet $(x_{A},z_{t_{A}=T},J_{t_{A}=T})$ from $\alpha=A$ as an element of a sample.

• •

Step 5. After repeating Steps 2 to 4, we obtain a sequence of triplets $\{(x_{A}^{(a)},z_{t_{A}=T}^{(a)},J_{t_{A}=T}^{(a)})\}\}$, which we use to estimate the expectation value:

[TABLE]

Note that the action with $t_{0}=0$ is the original action for which the sign problem exists. However, as can be seen from (2.14), the sample average is taken only with respect to the action with $t_{A}=T$, so that we need not worry about the original sign problem, although we include configurations near $t_{A}=0$. Also, $t_{A}=T$ can be taken to be sufficiently large such that the sign problem disappears if we complement intermediate flow times sufficiently (with larger $A$). Thus, this simple algorithm solves the two obstacles simultaneously: the original sign problem at $t_{0}=0$ is resolved at $t_{A}=T$ while the multimodal problem at $t_{A}=T$ is resolved by passing through configurations near $t_{0}=0$.777We have implicitly assumed that the action at $t_{0}=0$ does not cause a multimodality in the configuration space. If this is not the case, we further introduce other parameters (such as the coefficient of the action) as extra tempering parameters or prepare flow times $\{t_{\alpha}\}$ such that $t_{0}<0$.

3 Example

In this section we investigate the $(0+1)$-dimensional massive Thirring model at finite density [16, 17, 10] to exemplify that the algorithm given in the previous section does work.

The $(0+1)$-dimensional massive Thirring model is defined from the standard $(1+1)$-dimensional massive Thirring model by dimensional reduction. With an auxiliary field $\phi(\tau)$, the continuum representation of the grand partition function $Z={\rm tr\,}e^{-\beta(H-\mu Q)}$ is given by the path integral

[TABLE]

where the Euclidean action takes the form

[TABLE]

and the $\phi$ integral (the $\psi$, $\bar{\psi}$ integral) obeys the periodic (antiperiodic) boundary condition. In realizing the model on the lattice, we discretize the Euclidean time as $\tau=na$ $(n=1,\ldots,N)$ with $\beta=Na$ ($N$: even), and follow the prescription of Ref. [17], where $\phi(\tau)=\phi_{n}$ is treated as a $U(1)$ gauge potential and is combined with the chemical potential $\mu$ to become a link variable of a complexified gauge group, $e^{(i\phi(\tau)+\mu)a}=e^{i\phi_{n}a}e^{\mu a}\equiv U_{n}\,e^{\mu a}$, as proposed in Ref. [18]. Then, by using a staggered fermion formulation, the grand partition function is given by

[TABLE]

where $(dU)=\prod_{n}(dU_{n})\equiv\prod_{n}[d(\phi_{n}a)/2\pi]$ (with $\phi_{n}a\in(-\pi,\pi]$) and

[TABLE]

with D_{n\ell}(U)\equiv(1/2)\bigl{(}U_{n}e^{\mu a}\delta_{n+1,\ell}-U_{n-1}^{-1}e^{-\mu a}\delta_{n-1,\ell}-U_{N}e^{\mu a}\delta_{n,N}\delta_{\ell,1}+U_{N}^{-1}e^{-\mu a}\delta_{n,1}\delta_{\ell,N}\bigr{)}+ma\,\delta_{n\ell}. We henceforth set $a=1$ and treat $m$, $\mu$, and $g^{2}$ as dimensionless parameters.

After carrying out the fermion integration, we obtain

[TABLE]

which can be calculated analytically [17] as

[TABLE]

where $\alpha\equiv 1/(2g^{2})$, $2\rho_{\pm}\equiv(\sqrt{m^{2}+1}+m)^{N}\pm(\sqrt{m^{2}+1}-m)^{N}$, and $I_{n}(\alpha)$ is the modified Bessel function of the first kind of order $n$. Using this expression, we obtain the analytic form of the chiral condensate as

[TABLE]

We now numerically evaluate the chiral condensate by using the complex action

[TABLE]

as

[TABLE]

where

[TABLE]

It is easy to check that $[\det D(U;\mu)]^{\ast}=\det D(U;-\mu)$, so that the second term of (3.8) is complex-valued for real $\mu$, and the sign problem should arise when $N$ is large.

We follow the steps given in Sect. 2.2 with the complexification of the variables $\phi_{n}\in(-\pi,\pi]$. We first set the largest flow time to $T=2$, and prepare flow times $\{t_{\alpha}\}$ $(\alpha=0,1,2,\ldots,20)$ with equal separations as $t_{0}=0,\,t_{1}=0.1,\,t_{2}=0.2,\,\ldots,t_{20}=2.0$.888The fermion determinant at flow time $t=0$ is given by $\det D=(1/2^{N})\,[\zeta+\zeta^{-1}+(\sqrt{2}+1)^{N}+(\sqrt{2}-1)^{N})]$ where $\zeta\equiv e^{i\sum_{n}\phi_{n}+N\mu}$. Note that $\det D$ vanishes at $\sum_{n}\phi_{n}=\pi$ $(\mbox{mod}~{}2\pi)$ when $\mu=\log(\sqrt{2}+1)\sim 0.881$, which gives a multimodal distribution even at $t=0$. However, this can be handled without introducing another tempering parameter, because each mode has a rather wide distribution (with not-too-small values near the boundary) and thus the whole configuration space can be easily explored by setting the interval of proposal distribution to be a large value, as in this paper. Figure 1 shows $\phi(t)\equiv(1/N)\sum_{n}\phi_{n}(t)$ at flow times $t=t_{\alpha}$ with the $\phi_{n}(0)$ set to the same value.

As Step 1 of our algorithm, we make a cold start ($\phi_{n}=0$) for every replica $\alpha$, and numerically solve the differential equations (2.7) and (2.8) with the adaptive 4th-order Runge–Kutta method to obtain the triplet $(x_{\alpha},z_{t_{\alpha}},J_{t_{\alpha}})$. We repeat the Metropolis process twenty times in Step 2,999 As a proposal distribution we use the uniform distribution within the interval $[-\epsilon,\epsilon]$, where $\epsilon$ is chosen randomly from $\{1,\,10^{-1},\,10^{-2},\ldots,\,10^{-[2t_{\alpha}+1]}\}$ ($[k]$ is the floor of $k$). which is followed by a single sequence of swapping of Step 3. We then repeat Steps 2 and 3 ten times (as Step 4). With the first 5 data points discarded as initial sweeps, we estimate correlation functions with $10^{4}$ data points. Figure 2 shows the absolute value of the denominator in (2.14) (divided by the sample size) as a function of $\mu$ with the other parameters set to $N=8$, $m=1$, $g^{2}=1/2$.

The blue points are the result for the flow time $T=0$. They correspond to the usual reweighting calculus, and show that the sign problem actually exists for $\mu\gtrsim 1.0$. The green (red) points are the result for the flow time $T=2$ without (with) the parallel tempering implemented. The results show that the sign problem disappears at $T=2$.

Figure 3 shows the chiral condensate $\langle\bar{\chi}\chi\rangle$ as a function of $\mu$.

The other parameters are again set to $N=8$, $m=1$, $g^{2}=1/2$. The dotted line represents the analytic result (3.7). The blue points are the result for the flow time $T=0$ and exhibit large statistical errors, reflecting the sign problem. The green points are the result for the flow time $T=2$ without the parallel tempering implemented. They have small statistical errors, but exhibit statistically significant discrepancies from the analytic result. This should be attributed, as discussed in detail in Ref. [11], to the fact that the dominant contributions come only from a single thimble for such large $T$. The red points are the result for the flow time $T=2$ now with the parallel tempering implemented. They show a good agreement with the analytic result, which implies that contributions from various thimbles are correctly taken into account through the parallel tempering. This can be confirmed by Figs. 4 and 5. Figure 4 exhibits the histogram of $\phi=(1/N)\sum_{n}\phi_{n}$ at $T=2$.

We see that the configurations are concentrated on a single thimble if the parallel tempering is not implemented (left), while they are spread out on various thimbles if the parallel tempering is implemented (right). Figure 5 shows the average acceptance rates for the swaps between replicas $\alpha$ and $\alpha+1$ in Step 3.

We see that the swapping is carried out very well because the average acceptance rate is more than $50\%$ for all pairs $(\alpha,\alpha+1)$.

4 Conclusion and outlook

In this paper, as a versatile tool for Monte Carlo calculations of models with complex actions, we have proposed a Lefschetz thimble algorithm where the flow time is used as an auxiliary variable for the highly multimodal distribution. In particular, we implemented the parallel tempering method by taking the flow time $t$ as a tempering parameter. There, we prepare flow times $\{t_{\alpha}\}$ such that $t_{0}=0$ and $t_{A}=T$. The largest flow time $T$ can be taken to be sufficiently large such that the sign problem disappears there. Although the algorithm includes configurations at $t_{0}=0$, the original sign problem at $t_{0}=0$ does not enter the calculation because the sample average is taken only with respect to $t_{A}=T$, where the sign problem disappears. We have investigated the $(0+1)$-dimensional massive Thirring model at finite density to exemplify that the algorithm does work, and showed that contributions from multi thimbles are correctly taken into account even for such a large flow time as $T=2$.

We should investigate to what extent this algorithm is actually versatile. One interesting class of models, for which we can readily test our algorithm before applying it to QCD at finite density, is that of various types of large-$N$ random matrix models with complex actions. In fact, if the free energy is calculated by an integration over matrices themselves (not over their eigenvalues), the classical solutions and the corresponding Lefschetz thimbles do not have a useful meaning because the quantum corrections are of the same order as the leading term. It thus provides us with a good test of versatility to check whether correct results are obtained for such models where the thimble structure or its usefulness is not clear. As a related model, the numerical study of the triangle–hinge model [19, 20, 21] should also be interesting. The model is a sort of matrix model that generates 3D random volumes as a collection of triangles and hinges. In order to restrict the resulting configurations to tetrahedral decompositions, one needs to introduce a special form of interaction [19], which makes the action complex-valued (M. Fukuma, S. Sugishita, and N. Umeda, manuscript in preparation). A numerical study was made for a simplified model (with no restriction to tetrahedral decompositions), and the existence of a third-order phase transition is confirmed (manuscript in preparation). It is thus interesting to see whether the phase transition still exists when the restriction is imposed.

Besides the Lefschetz thimble algorithm, the complex Langevin algorithm [22, 23] is also under intense study as a promising method to solve the sign problem. Recently, a very interesting proposal was made by Bloch [24] (see also Ref. [25]) to evaluate correlation functions by reweighting complex Langevin trajectories using such parameters that satisfy known validity conditions [26, 27] to be free from wrong convergence problems [26, 28, 27, 29, 30]. It should be interesting to compare the extent of versatility between the reweighted complex Langevin algorithm and our parallel tempering algorithm with the flow time as a tempering parameter.

A study along these lines is now in progress and will be reported elsewhere.

Acknowledgements

The authors thank Hikaru Kawai, Jun Nishimura, Sotaro Sugishita, and Asato Tsuchiya for useful discussions. This work was partially supported by MEXT (Grant No. 16K05321).