Tempered transitions between thimbles

TL;DR

This paper introduces a tempered transition algorithm to improve Monte Carlo sampling in quantum field theories with complex actions, effectively addressing the sign problem by enabling sampling across multiple modes.

Contribution

The paper proposes a novel tempered proposals algorithm that enhances sampling efficiency in complex action quantum field theories, building on the holomorphic gradient flow method.

Findings

Tempered proposals successfully sample the 0+1 dimensional Thirring model at finite density.

The method overcomes the limitations of standard Monte Carlo sampling in multi-modal distributions.

The approach demonstrates potential for broader applications in complex quantum systems.

Abstract

Quantum field theories with complex actions cannot be investigated using importance sampling due to the sign problem. One possible solution is to use the holomorphic gradient flow, a method we introduced related to the Lefschetz thimbles idea. In many cases the probability distribution generated by this method is multi-modal and standard Monte-Carlo sampling fails. We propose an algorithm that incorporates tempered proposals to solve this problem. We apply this algorithm to the 0+1 dimensional Thirring model at finite density for a parameter set where standard sampling fails and show that tempered proposals cure this problem.

Tables

Table 1. Table 1 : Results for the chiral condensate at μ = 1.0 𝜇 1.0 \mu=1.0 for the exact solution, a Monte Carlo trapped to the global minimum thimble and a Monte Carlo utilizing tempered transitions.

Method	$⟨\overline{\psi }\psi ⟩$
Exact	0.575
with tempered transition	0.602(20)
w/o tempered transitions	0.470(02)

Equations

$$⟨\mathcal{O}⟩=\frac{{\int }_{{\mathbb{R}}^N}d{\zeta }_i{e}^{-S[\zeta ]}\mathcal{O}[\zeta ]}{{\int }_{{\mathbb{R}}^N}d{\zeta }_i{e}^{-S[\zeta ]}}=\frac{{\int }_{\mathcal{M}}d\varphi e^{-S[\varphi ]}\mathcal{O}[\varphi ]}{{\int }_{\mathcal{M}}d{\varphi }_i{e}^{-S[\varphi ]}},$$

$$\frac{d{\varphi }_i}{dt}=\overline{\frac{\partial S}{\partial {\varphi }_i}},{\varphi }_i(0)={\zeta }_i,$$

$$\frac{d{S}_R}{dt}=\frac{1}{2}\left(\frac{\partial S}{\partial {\varphi }_i}\frac{d{\varphi }_i}{dt}+\overline{\frac{\partial S}{\partial {\varphi }_i}}\overline{\frac{d{\varphi }_i}{dt}}\right)={\left|\frac{\partial S}{\partial {\varphi }_i}\right|}^2>0,\frac{d{S}_I}{dt}=\frac{1}{2i}\left(\frac{\partial S}{\partial {\varphi }_i}\frac{d{\varphi }_i}{dt}-\overline{\frac{\partial S}{\partial {\varphi }_i}}\overline{\frac{d{\varphi }_i}{dt}}\right)=0$$

$$Z={\int }_{\mathcal{M}}d{\varphi }_i{e}^{-S[\varphi ]}={\int }_{{\mathbb{R}}^N}d{\zeta }_i\det \left(\frac{\partial {\varphi }_i}{\partial {\zeta }_j}\right)e^{-S[\varphi (\zeta )]}.$$

$$\frac{d{J}_{ij}}{dt}=\overline{H_{ik}}\overline{J_{kj}},H_{ij}\equiv \frac{{\partial }^{2}S}{\partial {\varphi }_i\partial {\varphi }_k},J_{ij}(0)={\delta }_{ij},$$

$$⟨\mathcal{O}⟩$$

$$p_i({\varphi }_i){\hat{T}}_{i+1}({\hat{\varphi }}_i\to {\hat{\varphi }}_{i+1})=p_i({\varphi }_{i+1}){\hat{T}}_{i+1}({\hat{\varphi }}_{i+1}\to {\hat{\varphi }}_i).$$

$$\mathcal{P}({\varphi }_0\to {\varphi }_0^{\prime })$$

$$\mathcal{A}({\varphi }_0\to {\varphi }_0^{\prime })=\min \left(1,\frac{p_1({\varphi }_0)}{p_0({\varphi }_0)}\cdots \frac{p_n({\varphi }_{n-1})}{p_{n-1}({\varphi }_{n-1})}\frac{p_{n-1}({\varphi }_{n-1}^{\prime })}{p_n({\varphi }_{n-1}^{\prime })}\cdots \frac{p_0({\varphi }_0^{\prime })}{p_1({\varphi }_0^{\prime })}\right).$$

$$p({\varphi }_0)\mathcal{P}({\varphi }_0\to {\varphi }_0^{\prime })\mathcal{A}({\varphi }_0\to {\varphi }_0^{\prime })=p({\varphi }_0^{\prime })\mathcal{P}({\varphi }_0^{\prime }\to {\varphi }_0)\mathcal{A}({\varphi }_0^{\prime }\to {\varphi }_0).$$

$$Z=\left[\prod _{t=1}^N{\int }_0^{2\pi }\frac{d{\varphi }_t}{2\pi }\right]\det D{e}^{-\frac{1}{2g^2}{\sum }_{t=1}^N(1-\cos {\varphi }_t)}\equiv \left[\prod _{t=1}^N{\int }_0^{2\pi }\frac{d{\varphi }_t}{2\pi }\right]e^{-S[\varphi ]},$$

$$S[\varphi ]$$

$$D_{t,t^{\prime }}$$

$$⟨\bar{\psi }\psi ⟩=\frac{1}{N}\frac{\partial }{\partial m}\log Z=\frac{(1+m^2{)}^{-1/2}{I}_0^N(\frac{1}{2g^2})\sinh (N{\sinh }^{-1}(m))}{I_1^N(\frac{1}{2g^2})\cosh (N\mu )+I_0^N(\alpha )\cosh (N{\sinh }^{-1}(m))}$$

Full text

I Introduction

Non-perturbative results from quantum field theories can be obtained using stochastic sampling, as long as the path integral can be represented as a sum over real and positive contributions. This is no longer the case when one considers phenomenologically interesting problems regarding systems at finite density, like QCD at non-zero baryon density, or questions related to real-time dynamics. In these cases, the integrand is complex and direct Monte-Carlo sampling cannot be applied. The standard work-around, reweighting, that samples according to a positive measure and corrects for the difference by introducing a complex fluctuating phase in the observables, fails due to phase oscillations; this is the infamous sign problem. Recently a possible solution was proposed by Cristoforetti et al. Cristoforetti:2012su : We start by embedding the integration domain of the path integral in a complex space. Using the analytical properties of the integrand we then deform the integration manifold without changing the integral value. The proposal by Cristoforetti was to use the a manifold that corresponds to a union of Lefschetz thimbles. This thimble decomposition is always possible and it has the advantage that the integrand’s complex phase on each thimble is constant. When the integral is dominated by the contribution of one thimble, the sign problem is effectively solved Cristoforetti:2013wha ; Mukherjee:2013aga ; Fujii:2013sra ; Alexandru:2015xva ; Alexandru:2016san .

In cases where more than one thimble contributes significantly to the integral, the problem is significantly harder. Sampling algorithms have to be able not only to sample each thimble but also be able to dynamically transition between thimbles in order to properly take into account their relative contribution. Additionally, significant analytical work is needed to identify all the thimbles contributing to the integral. In an earlier paper Alexandru:2015sua we proposed a method that sidesteps this task. The idea is to use a class of manifolds generated by the holomorphic gradient flow and parametrized by the flow time $T_{\text{flow}}$ interpolating smoothly between the original integration domain ($T_{\text{flow}}=0$) and the thimble decomposition ($T_{\text{flow}}=+\infty$). While the value of the integral on all these manifolds is the same, the phase fluctuations become progressively smaller as we increase $T_{\text{flow}}$, allowing us to use reweighting. On the other hand as $T_{\text{flow}}$ increases the probability distribution becomes multi-modal with growing potential barriers between different modes. For algorithms that rely on small-step updates, such distributions are difficult to sample since the transitions rate between modes becomes very small. For some systems there are values of $T_{\text{flow}}$ for which the sign problem is mild and the transition rate between modes is good such that standard small-step algorithms can be used Alexandru:2016gsd ; Alexandru:2016ejd .

For systems where the sign problem only becomes manageable when $T_{\text{flow}}$ is large and the probability distribution is multi-modal a possible solution is to use the method of tempered transitions. The basic idea is to use a set of small-steps to build a large mode-to-mode move Neal1996 . The sequence of steps is constructed using a set of guiding distributions that overlap well sequentially but gradually lower and raise the potential barriers between modes. In this paper we discuss a proposal where the guiding probabilities are generated by changing the $T_{\text{flow}}$ and apply it to the Thirring model in $0+1$ dimensions, a system whose path integral decomposition requires multiple thimbles Alexandru:2015xva . The plan of the paper is the following: in Section II we review the relevant details for the holomorphic flow and thimble decomposition, in Section III we review the details of the tempered transitions algorithm as it applies to our problem, and in Section IV we review the relevant details for the $0+1$ Thirring model and present the results of our simulations.

II Holomorphic Gradient Flow

Here we show how to deform the domain of path integration in order to ameliorate the sign problem Cristoforetti:2013wha ; Alexandru:2015sua . The starting point is Cauchy’s theorem, which allows one to deform the domain of path integration ($\mathbb{R}^{N}$) into a submanifold ${\cal M}$ of complex space ($\mathbb{C}^{N}\approx\mathbb{R}^{2N}$) without changing the value of the path integral:

[TABLE]

where $\zeta_{i}$, $i=1,\ldots,N$ are real field variables. The sign problem arises because $S[\zeta]$ is not real, leading to rapid phase oscillations in the path integral. The goal is to find a manifold, ${\cal M}$, such that Cauchy’s theorem applies and $S[\phi]$ does not oscillate as rapidly for $\phi\in{\cal M}$ as it does for $\zeta\in\mathbb{R}^{N}$. One way to construct such a manifold is to identify every field configuration in the original integration domain ($\mathbb{R}^{N}$) as an initial condition for the following set of first order differential equations known as the holomorphic gradient flow equations:

[TABLE]

integrated up to a fixed “flow time” $T_{\text{flow}}$. Here the bar on the RHS of Eq. 2 denotes complex conjugation. These equations map a particular field configuration, $\zeta\in\mathbb{R}^{N}$ to a point $\phi(T_{\text{flow}})\in\mathbb{C}^{N}$. We will call this motion the “flow”. The map defined by the flow $\zeta\mapsto\phi(T_{\text{flow}})$ is one-to-one as Eq. 2 is first order. Therefore flowing $\mathbb{R}^{N}$ generates a N real-dimensional manifold in ${\cal M}\subset\mathbb{C}^{N}$ (i.e. an $N$ real-dimensional manifold, ${\cal M}$, embedded in $\mathbb{C}^{N}\simeq\mathbb{R}^{2N}$).

Having constructed ${\cal M}$ we now to establish that (i) Cauchy’s theorem applies on ${\cal M}$ $($so the path integral on ${\cal M}$ is equal to the path integral on \mathbb{R}^{N}$$) and (ii) the phase oscillations on ${\cal M}$ are milder than the phase oscillations on $\mathbb{R}^{N}$, which leads to a milder sign problem. First observe an important property of the flow equations: the real part of the action, $S_{R}$, increases monotonically along a flow trajectory, whereas the imaginary part, $S_{I}$, stays constant:

[TABLE]

where in the second equalities we used Eq. 2. For (i) to hold, it must be that points do not cross any singularities throughout the continuous deformation process, especially singularities that might arise when the field variables go to infinity. Indeed no singularities are crossed because $S_{R}$ increases monotonically with the flow Eq. 3, which ensures that the modulus of the integrand, $|e^{-S}|=e^{-S_{R}}$ decreases monotonically, and therefore remains bounded from above and damps the integral exponentially as the fields approach infinity 111We assumed that the path integral was convergent to begin with. This is indeed the case when the lattice spacing is finite. The standard renormalization procedure has to be followed to approach the continuum limit..

For property (ii), consider the limit of large flow time, $T_{\text{flow}}$. In this limit, almost all the original field configurations in $\mathbb{R}^{N}$ will flow into configurations with very large $S_{R}$ and will contribute practically nothing to the path integral due to the exponential damping factor $e^{-S_{R}}$. Therefore in the large flow limit, the main support of the path integral will come from fields that flow to fixed points of Eq. 2 which are critical points of the action $\partial S/\partial\phi_{i}=0$ (i.e. classical solutions to the complexified equations of motion). Consider a point in $\mathbb{R}^{N}$ that flows to a critical point. Then an infinitesimal neighborhood around this point will flow to an $N$-real dimensional manifold ${\cal J}$, attached to that critical point. Since $S_{I}$ remains unchanged with the flow and the variation of the $S_{I}$ in an infinitesimal neighborhood is infinitesimal, $S_{I}$ will be approximately constant on ${\cal J}$. In the limit $T_{\text{flow}}\to\infty$, $S_{I}$ will be exactly constant. This $N$ dimensional manifold ${\cal J}$ attached to a critical point, over which $S_{I}$ is constant is known as a “Lefschetz thimble” $($the multi-dimensional generalization of the stationary phase/steepest descent contour familiar from complex analysis$)$. In the limit $T_{\text{flow}}\to\infty$, ${\cal M}$ will be a particular combination of thimbles, and $S_{I}$ will be piecewise constant on ${\cal M}$. For sufficiently large, but finite $T_{\text{flow}}$, $S_{I}$ will not exactly be piecewise constant, but will be approximately piecewise constant. Consequently, by tuning $T_{\text{flow}}$ we can continuously soften the severity of the phase oscillations caused by $e^{iS_{I}}$. In other words, the manifolds defined by different $T_{\text{flow}}$ interpolate between $\mathbb{R}^{N}$ ($T_{\text{flow}}=0$) where $S_{I}$ varies rapidly and the associated thimble decomposition of the path integral ($T_{\text{flow}}\to\infty$) where $S_{I}$ is piecewise constant.

It is desirable, from the point of view of the sign problem, to integrate over a highly flowed manifold because the only regions with appreciable support have nearly constant $S_{I}$. However there are costs to this procedure. First, it is numerically expensive to do so. More fundamentally, however, regions of support on the parameterization manifold $($where the Monte Carlo takes place$)$ of a highly flowed surface are separated in field space by extended regions with very large $S_{R}$. Consequently a Monte Carlo simulation must sample from a multi-modal distribution. The situation is worse when there are fermions involved, because in these cases, thimbles have boundaries: $N-1$ dimensional sub-manifolds where the fermion determinant vanishes and $S_{R}\rightarrow\infty$. As an illustration of this phenomenon we shown in Fig. 1 the action on a submanifold of ${\cal M}$ for various flow times. One way to avoid this situation is, instead of approaching the $T_{\text{flow}}\to\infty$ limit, to tune $T_{\text{flow}}$ to a finite value where the phase oscillations are under control yet the action barriers are not so high that the Monte Carlo gets trapped. In a variety of examples this is indeed possible Alexandru:2016gsd ; Schmidt:2017gvu ; Alexandru:2016san . However, depending on the model and, in particular at large volumes, it may not be possible to solve the sign problem and the isolated minima problem simultaneously by tuning the flow time. Thus we come to the central point of this paper: how to tackle multimodal distributions created by the holomorphic gradient flow.

Before concluding this section we outline the Monte Carlo computation that we will employ later on. Using the fact that the flow defines a one-to-one mapping between the initial field, $\phi_{i}(0)=\zeta_{i}$ and the flowed field $\phi_{i}(T_{\text{flow}})$ we can parameterize the path integral over ${\cal M}$ using real variables $\zeta_{i}$

[TABLE]

Notice that this a re-parameterization of ${\cal M}$ and is distinct from the contour deformation from $\mathbb{R}^{N}$ to ${\cal M}$ that we discussed earlier. Parameterizing ${\cal M}$ with real fields allows us to perform the Metropolis updates on $\mathbb{R}^{N}$. The Jacobian, $J_{ij}=\partial\phi_{i}/\partial\zeta_{j}$, associated with this change of variables also satisfies a flow equation

[TABLE]

which transports the local tangent space at $\zeta_{i}$ to the flowed point $\phi_{i}(T_{\text{flow}})$ along the flow trajectory followed by $\zeta_{i}$. The determinant of this Jacobian is a complex number which we combine with the action to define an effective action, $S_{\text{eff}}[\zeta]=S[\phi(\zeta)]-\log\det J$. In our Monte Carlo simulations, the configurations are sampled according to the real part of the effective action, $\operatorname{Re}S_{\text{eff}}[\zeta]=S_{R}[\phi(\zeta)]-\log|\det J|$. The phase of the Jacobian, along with the phase of the action, $\varphi(\zeta)\equiv\operatorname{Im}S_{\text{eff}}[\zeta]=S_{I}[\phi(\zeta)]-\mathop{\hbox{Im}}\det J$, is included via “reweighing”:

[TABLE]

Finally, even though we are performing a path integral over $\mathbb{R}^{N}$ due to the way we parameterize ${\cal M}$, the action and all the operators are still evaluated on ${\cal M}$. In particular the fluctuations of $S_{I}[\phi(\zeta)]$ which enter in are drastically milder than those over the original domain given by $S_{I}[\zeta]$.

III Tempered transitions

The method of tempered transitions was designed to perform Monte Carlo calculations in situations where the desired probability distribution is multimodal, that is, has more than one, well separated peaks Neal1996 . Multimodal distributions are challenging for Monte Carlo methods because, with most algorithms, the Monte Carlo chain ends up being trapped in one of the modes for a very large number of steps making it nearly impossible to sample properly.

We now quickly describe the method of tempered transitions in general terms, before applying it to a specific model. Suppose the distribution of interest is $p(\phi)$. Then a standard importance sampling technique is to make symmetric proposals $\phi\rightarrow\phi^{\prime}$ in the sampling space and accept such proposals with probability $\text{min}\{1,\frac{p(\phi^{\prime})}{p(\phi)}\}$. In order to achieve a reasonable acceptance rate, the proposed $\phi^{\prime}$ is chosen close to $\phi$. For a multimodal distribution this leads to the trapping alluded above. In the tempered transitions methods one makes a more sophisticated $($and computationally expensive$)$ proposal that has a fair likelihood of being on another mode and also of being accepted. This is achieved by considering a sequence of $n+1$ progressively flatter probability distributions $p_{i}(\phi)$, with $p_{0}(\phi)=p(\phi)$ being the desired distribution to sample from. For each of these distributions $p_{i}(\phi)$ consider a transition probability $\hat{T}_{i+1}(\hat{\phi}_{i}\rightarrow\hat{\phi}_{i+1})$ satisfying detailed balance with respect to $p_{i}(\phi)$:

[TABLE]

The transition probabilities $\hat{T}_{i}(\hat{\phi}_{i-1}\rightarrow\hat{\phi}_{i})$ can be chosen to be, for instance, Metropolis steps.

The proposed configuration $\phi^{\prime}$ is obtained from $\phi=\phi_{0}$ by evolving $\phi$ through a series of updates with the transition probabilities $\hat{T}_{i}(\hat{\phi}_{i}\rightarrow\hat{\phi}_{i+1})$ all the way up to $\hat{T}_{n}(\hat{\phi}_{n-1}\rightarrow\hat{\phi}_{n})$ and then again, in reverse order, down to $\hat{T}_{1}(\hat{\phi}_{1}^{\prime}\rightarrow\hat{\phi}_{0}^{\prime})$ (see Fig. 2). More precisely, the probability of proposing $\phi^{\prime}=\phi_{0}^{\prime}$ starting from $\phi=\phi_{0}$ along the path $\phi_{0}\rightarrow\phi_{1}\rightarrow...\rightarrow{\phi_{1}}^{\prime}\rightarrow{\phi_{0}}^{\prime}$ is given by:

[TABLE]

The final configuration $\phi_{0}^{\prime}$ is the proposed configuration that is accepted or not according to the acceptance probability

[TABLE]

It is straightforward to verify that the combination of a tempered transition and the accept/reject step satisfy detailed balance and thus samples the true distribution $p(\phi)$:

[TABLE]

While Eq. 10 guarantees the correctness of the method, its usefulness relies on its ability to generate proposals transitioning between modes with a high probability of acceptance. A heuristic discussion of how to choose the intermediate probabilities $p_{i}(\phi)$ is presented in Neal1996 . In summary, it is necessary that the probabilities $p_{i}(\phi)$ at the “top” of the ladder in Fig. 2 do not have well separated modes and that the subsequent distribution probabilities ($p_{i}(\phi)$ and $p_{i+1}(\phi)$) be close enough so the Monte Carlo chain along up and down the ladder in Eq. 8 can gently guide the configuration from one mode to another.

In our implementation of tempered transitions, we modulate the probability distributions up the ladder by adjusting the amount of flow we subject the parameterization manifold to. It can be seen from Fig. 1 that at $T_{\text{flow}}=0.0$, the action barriers on ${\cal M}$ are mild while for large $T_{\text{flow}}$ the action barriers are high. Therefore, we choose $p_{0}(\phi)$ to be the probability distribution generated with a large enough $T_{\text{flow}}$ to tame the sign problem and we choose $p_{n}(\phi)$ to be the probability distribution generated when $T_{\text{flow}}=0.0$ where fields are mobile. There much latitude in how to choose the probability distributions between $p_{0}(\phi)$ and $p_{n}(\phi)$, but not all choices perform the same. We choose to interpolate between a $T_{\text{flow}}=T_{\text{max}}$ and $T_{\text{flow}}=0.0$ linearly. Whatever the optimal choice may be, it is clear that $p_{i}(\phi)$ and $p_{i+1}(\phi)$ should be similar in some sense for such a choice. This is because if all $p_{i}(\phi)$ are the same then the tempered proposal is guaranteed to be accepted.

IV Results

We now apply the method of tempered transitions on $0+1$ Thirring model at finite density. For non-zero chemical potential the determinant of the Dirac matrix is complex and the theory suffers from a sign problem. We discretize the Euclidean time direction using staggered fermions. The lattice partition function for this theory is

[TABLE]

where the effective action and the Dirac matrix are explicitly given by

[TABLE]

Here $N$ is an even number equal to to the number of lattice sites and all dimensionful quantities are measured in units of the lattice spacing $a$ which we set to one. This discretized model is exactly solvable Pawlowski:2013pje ; an observable of interest is the chiral condensate, which is given by:

[TABLE]

where $I_{n}(z)$ denotes the modified Bessel function of the first kind of order $n$. We will use this exact solution to compare against.

As an example we take an $N=16$ lattice with $m=1.0$ and $g^{2}={1\over 6}$ over a range of chemical potentials. First, we reduce the sign problem for these parameters by shifting the domain of integration from the real hyperplane to the tangent plane of the thimble fixed to the global minimum of the action Alexandru:2015sua where the phase fluctuations are smaller than on the real plane. The remaining phase fluctuations can be tamed by flowing the tangent plane. Calculations performed with flow times in the range $0.2\lesssim T_{\text{flow}}\lesssim 0.4$ shows that at least three thimbles contribute significantly (see Fig. 3). For larger values of flow time, calculations performed with the Metropolis algorithm remain trapped in the region with $\langle\phi\rangle\approx 0$ for a very large number of steps (we have followed the Monte Carlo chain up to $10\times 10^{6}$ steps without seeing a transition to other thimbles). For some of the parameters explored, in particular when the chiral condensate drops quickly around the mass of the fermion $\psi$, computing observables when trapped to the global minimal thimble yields statistically incorrect results, see Fig. 4. The need to integrate over several thimbles at points of sharp variation in thermodynamic functions has been stressed before in fermionic models Tanizaki:2016cou .

As the discrepancy between the trapped numerical result and the exact result is largest at $\mu=1.0$, we restrict our further analysis to this value of the chemical potential. We apply the method of tempered transitions using a ladder with 1600 steps in each direction. At the bottom of the ladder the flow is 0.5 where the fields are trapped to a single thimble, and at the top the flow is 0.0 were the Monte Carlo chain explores all thimbles. We interpolate linearly in the flow as sub-transitions are applied. At each step of the ladder, we use a transition operator $T_{k}(\phi,\phi^{\prime})$ which applies 10 standard metropolis proposals at the flow time $t_{k}$ to the configuration $\phi$. The reason for using such an operator is to allow the newly relaxed distribution to equilibrate before relaxing it further. We find that this procedure reduces the “free energy” cost of a tempered transition and increases the likelihood of a tempered transition to be accepted.

There are two time scales in the equilibration of the Monte Carlo chain: a fast one for the equilibration within one mode and a slow one for the equilibration between modes. These separated time scales exist for the simple reason that only small proposals are required to explore an individual thimble while a large proposal is required to transition. This is a generic property of theories with multiple thimbles. In light of this observation, we construct our Monte Carlo as follows: between two tempered transitions we perform 1000 standard metropolis steps. This allows all relevant thimbles to be explored and the space within each thimble to be explored. Moreover such a division is necessary because tempered transitions are expensive compared to normal transitions. In spite of this cost, tempered proposals induce transitions at a substantially higher rate than standard Metropolis proposals. In this study, a tempered transition is composed of $10\times 3200=3.2\times 10^{4}$ individual metropolis steps. We find that roughly every $10^{\text{th}}$ proposal between the main thimble and the shoulder thimbles is accepted, so it takes the comuptational effort of $\sim 3\times 10^{5}$ metropolis steps to induce a transition. Therefore tempered proposals are at least $1-2$ orders of magnitude more efficient than standard metropolis proposals in this study.

The results of a Monte Carlo using 2000 tempered transitions is shown in Fig. 5. We find that the combination of tempered transitions with 1600 ladder steps and regular Metropolis steps is sufficient to transition regularly between thimbles. In addition, we find that the inclusion of the neighboring thimbles reproduces the exact result up to statistical errors as can be seen in Table 1. Roughly $10\%$ of tempered transitions from the central thimble to the shoulder thimbles were accepted. Note however that this low acceptance rate is not due to poor proposals but because the shoulder thimbles carry roughly $10\%$ of the weight of the path integral. In fact, the tempered transitions make proposals to the next-to-nearest shoulder thimbles frequently. This can be seen in Fig. 6 where the tempered proposals are plotted as a function of Monte Carlo step. These transitions to the next-to-nearest shoulder thimbles are, however, not accepted because these thimbles contribute very little to the path integral. It is expected, then, that in a problem where many thimbles carry substantial weight, tempered transitions may provide a natural means to explore many thimbles without any a priori knowledge of their location.

V Conclusion

The holomorphic gradient flow approach to the computation of path integrals with a sign problem frequently leads to multimodal probability distributions. Each of the modes are related to a thimble contributing to the integral, as described by Picard-Lefschetz theory. This poses a challenge to numerical computations as Monte Carlo chains tend to get “stuck” in one of the modes for exponentially long times. We applied the method of tempered transitions to this problem. We take the tempering parameter, which controls the steepness of the probability landscape in each ladder of the tempering process, to be the flow time Fukuma:2017fjq by which the real space is transported by the holomorphic gradient flow. We demonstrated in a simple example that this procedure is feasible and that it allows for the proper sampling of the field space even in circumstances where a simpler Metropolis algorithm fails. It was found that a combination of tempered steps interspersed with regular Metropolis steps was the most efficient choice.

The method is not without drawbacks. The main one is the need, even in the simple model considered here, of a very large number of ladders steps during a tempered proposal and the associated large computational cost. This is the main difficulty that limited us in this paper to fairly small toy models. Our experience in scaling up the number of degrees of freedom is that the number of ladder steps required scales roughly linearly with the number of degrees of freedom. By itself this is not such a steep scaling but it should be kept in mind that other, steeper increases in cost are caused by the computation of the jacobian, increase in required flow time and, in models with more spacetime dimensions, increase on the size of the Dirac matrix (in the $0+1$ dimensional model discussed here the determinant has a closed form). This is not to say that significant improvements are not possible by adjusting some parameters in our simulations. For instance, there is tremendous latitude in choosing the values of the intermediate flows. Most applications of the tempered transitions method use the temperature as the parameter changing along the tempered ladder, leading to an exponential flattening of the probabilities distributions Neal1996 . A similar exponential flattening of probabilities can be obtained by choosing the intermediate flows to be equally spaced along the ladder. It could well be that a different choice is significantly better.

As we finished the present paper a study of the same model we discuss using a different tempering method appeared Fukuma:2017fjq . A direct comparison of the efficacy of the two methods is hindered by some small differences between the two calculations that are unrelated to the tempering method. For instance, in Fukuma:2017fjq , the manifolds of integration are obtained by flowing the real manifold, not the main tangent manifold as is done here. It is clear, however, that either method will be extremely costly when applied to realistic field theories.

Acknowledgements.

A.A. is supported in part by the National Science Foundation CAREER grant PHY-1151648 and by U.S. DOE Grant No. DE-FG02-95ER40907. A.A. gratefully acknowledges the hospitality of the Physics Department at the University of Maryland where part of this work was carried out. G.B. is supported by the U.S. Department of Energy under Contract No. DE-FG02-01ER41195. P.B. and N.C.W. are supported by U.S. Department of Energy under Contract No. DE-FG02-93ER-40762.