# A Multilevel Blocking Approach to the Sign Problem in Real-Time Quantum Monte Carlo Simulations

## Abstract

We propose a novel approach toward the general solution of the sign problem in real-time path-integral simulations. Using a recursive multilevel blocking strategy, this method circumvents the sign problem by synthesizing the phase cancellations arising from interfering quantum paths on multiple timescales. A number of numerical examples with one or a few explicit degrees of freedom illustrate the usefulness of the method.

Path integrals  provide an elegant alternative to the operator formulation of quantum mechanics. Because they are easily adapted to many-body systems, quantum Monte Carlo (QMC) simulations based on path integrals can potentially yield exact results for the dynamics of condensed phase quantum systems. A number of early attempts to use QMC simulations for real-time dynamics  demonstrated their potential, but these studies also uncovered the ubiquitous “dynamical sign problem” — interference among quantum paths leads to large statistical noise that increases linearly with the number of possible paths, which in turn grows exponentially with the timescale of the problem. Consequently, real-time QMC simulations have been limited to problems of very short timescales. Several attempts to extend the timescale of real-time QMC simulations have appeared , all of which rely on a common idea — by using various filtering schemes, the phase cancellations can be numerically stabilized by damping out the non-stationary regions in configuration space. Such filtering methods were able to extend the timescale somewhat, but they were generally unable to reach timescales of interest in typical chemical systems.

Here we propose a novel approach based on a blocking strategy which may provide a general solution to the dynamical sign problem for very long times. The blocking strategy asserts that by sampling groups of paths called blocks, the sign problem is always reduced compared to sampling single paths . Though this blocking strategy seems simple, its practical implementation is cumbersome, especially when going out to long times. Because the number of paths grows exponentially with the timescale, the number of blocks also grows immensely. To cure this problem, we first define elementary blocks and then group them together into larger blocks. Blocks of different sizes are introduced on several timescales called levels. After taking care of the sign cancellations within all blocks on a finer level, the entire sign problem can then be transferred to the next coarser level. By recursively proceeding from the bottom level (shortest timescale) up to the top (longest timescale), the troublesome numerical instabilities associated with the sign problem can be systematically avoided. In a slightly different form, this multilevel blocking (MLB) algorithm has recently been applied to treat the “fermion sign problem” in many-fermion imaginary-time simulations . As most technical details of the MLB scheme can be found in , we only give a brief description of the algorithm below, stressing the main ideas and its differences from the fermion formulation.

For concreteness, we describe the MLB method for the calculation of an equilibrium time-correlation function,

$$A(0)B(t)=\frac{\mathrm{Tr}\left\{e^{(\beta \mathrm{}+it)H/\mathrm{}}Ae^{+itH/\mathrm{}}B\right\}}{\mathrm{Tr}\left\{e^{(\beta \mathrm{}+it)H/\mathrm{}}e^{+itH/\mathrm{}}\right\}}.$$
(1)   
With minor modifications, the MLB method also applies to other dynamical properties like the thermally symmetrized correlation function , $`C_s(t)=Z^1\mathrm{Tr}\{e^{(\beta \mathrm{}/2+it)H/\mathrm{}}Ae^{(\beta /\mathrm{}/2it)H/\mathrm{}}B\}`$, with $`Z`$ being the partition function. In terms of path integrals, the traces in (1) involve two quantum paths, one propagated backward in time for the duration $`t`$ and the other propagated in complex time for the duration $`ti\beta \mathrm{}`$. Discretizing each of the two paths into $`P`$ slices, the entire cyclic path has a total of $`2P`$ slices. A slice on the first half of them has length $`t/P`$, and on the second half $`(ti\beta \mathrm{})/P`$. We require $`P=2^L`$ which defines the total number of levels $`L`$. Denoting the quantum numbers (e.g., spin or position variables) at slice $`j`$ by $`𝒓_j`$, $`\{𝒓_1,\mathrm{}𝒓_{2P}\}`$ is a discrete representation of a path, and the correlation function (1) reads

$$\frac{𝑑𝒓_1\mathrm{}𝑑𝒓_{2P}B(𝒓_{2P})A(𝒓_P)_{j=1}^{2P}(𝒓_j,𝒓_{j+1})_0}{𝑑𝒓_1\mathrm{}𝑑𝒓_{2P}_{j=1}^{2P}(𝒓_j,𝒓_{j+1})_0},$$
(2)   
where the level-0 bond $`(𝒓_j,𝒓_{j+1})_0`$ is simply the short-time propagator between slices $`j`$ and $`j+1`$, and $`𝒓_{2P+1}=𝒓_1`$. A direct application of the QMC method would sample these paths using the modulus of the integrand in the denominator as the weight.

We first assign all slices along the discretized path to different levels $`\mathrm{}=0,\mathrm{},L`$ (see Figure 1). Each slice $`j=1,\mathrm{},2P`$ belongs to a unique level $`\mathrm{}`$, such that $`j=(2k+1)2^{\mathrm{}}`$ and $`k`$ is a nonnegative integer. For instance, slices $`j=1,3,5,\mathrm{}`$ belong to level $`\mathrm{}=0`$, slices $`j=2,6,10,\mathrm{}`$ to $`\mathrm{}=1`$, etc. The MLB algorithm starts by sampling only configurations which are allowed to vary on slices associated with the finest level $`\mathrm{}=0`$, using the weight $`𝒫_0=|(𝒓_1,𝒓_2)_0\mathrm{}(𝒓_{2P},𝒓_1)_0|`$. The short-time level-$`0`$ bonds are then employed to synthesize longer-time level-$`1`$ bonds that connect the even-$`j`$ slices. Subsequently the level-1 bonds are used to synthesize level-$`2`$ bonds, and so on. In this way the MLB algorithm moves recursively from the finest level ($`\mathrm{}=0`$) up to increasingly coarser levels until $`\mathrm{}=L`$, where the measurement is done using $`𝒓_{2P}`$ and $`𝒓_P`$.

Generating a MC trajectory containing $`K`$ samples for each slice on level $`\mathrm{}=0`$ and storing these samples, we compute the level-1 bonds according to

$`(𝒓_j,𝒓_{j+2})_1`$ $`=`$ $`\mathrm{\Phi }[(𝒓_j,𝒓_{j+1})_0(𝒓_{j+1},𝒓_{j+2})_0]_{𝒫_0}`$ (3)   
$`=`$ $`K^1{\displaystyle \underset{𝒓_{j+1}}{}}\mathrm{\Phi }[(𝒓_j,𝒓_{j+1})_0(𝒓_{j+1},𝒓_{j+2})_0],`$ (4)   
where the summation $`_{𝒓_{j+1}}`$ extends over the $`K`$ samples, and $`\mathrm{\Phi }[z]=e^{i\mathrm{arg}(z)}`$ denotes the phase. For a complete solution of the sign problem, the sample number $`K`$ has to be sufficiently large . The bonds (3) contain crucial information about the sign cancellations on the previous level $`\mathrm{}=0`$. Their benefit becomes clear when rewriting the integrand of the denominator in (2) as

$`𝒫_0\times (𝒓_2,𝒓_4)_1\mathrm{}(𝒓_{2P2},𝒓_{2P})_1(𝒓_{2P},𝒓_2)_1.`$

Comparing this to (2), we notice that the entire sign problem has been transferred to the next coarser level.

In the next step, the sampling is carried out on level $`\mathrm{}=1`$ in order to compute the next-level bonds, using the weight $`𝒫_0𝒫_1`$ with $`𝒫_1=|(𝒓_2,𝒓_4)_1\mathrm{}(𝒓_{2P},𝒓_2)_1|`$. Generating a sequence of $`K`$ samples for each slice on level $`\mathrm{}=1`$, and storing these samples, we then calculate the level-2 bonds in analogy with (3),

$`(𝒓_j,𝒓_{j+4})_2=\mathrm{\Phi }\left[(𝒓_j,𝒓_{j+2})_1(𝒓_{j+2},𝒓_{j+4})_1\right]_{𝒫_0𝒫_1},`$

and iterate the process up to the top level by employing analogously defined level-$`\mathrm{}`$ bonds. The correlation function (1) can then be computed from

$$\frac{B(𝒓_{2P})A(𝒓_P)\mathrm{\Phi }[(𝒓_P,𝒓_{2P})_L(𝒓_{2P},𝒓_P)_L]_𝒫}{\mathrm{\Phi }[(𝒓_P,𝒓_{2P})_L(𝒓_{2P},𝒓_P)_L]_𝒫},$$
(5)   
with the positive definite MC weight $`𝒫=𝒫_0𝒫_1\mathrm{}𝒫_L`$. The denominator in (5) gives the average phase and indicates to what extent the sign problem has been solved. Under the direct QMC method, the average phase decays exponentially with $`t`$ and is typically close to zero. With the MLB algorithm, however, the average phase remains close to unity even for long times, with a CPU time requirement $`t`$. The price to pay for the stability of the algorithm is the increased memory requirement $`K^2`$ associated with having to store the sampled configurations.

Now we illustrate the practical usefulness of the MLB method by several numerical examples. In each of these examples, we compute a time-correlation function. The average phase is larger than 0.6 for all data sets shown below. The decay in the average phase with $`t`$ is a result of the finiteness of $`K`$ . Choosing a larger $`K`$ allows for a larger average phase out to longer time at the cost of increased computer memory and CPU time. Each data point in even the most intensive calculation took no more than a few hours on an IBM RS 6000/590.

A. Harmonic oscillator. For $`H=p^2/2m+m\omega ^2x^2/2`$, the real and imaginary parts of $`x(0)x(t)`$ oscillate in time due to vibrational coherence. In dimensionless units $`m=\omega =1`$, the oscillation period is $`2\pi `$. Figure 2(a) shows MLB results for $`C(t)=\mathrm{Re}x(0)x(t)`$. With $`P=32`$ for the maximum time $`t=26`$, $`K=200`$ samples were used for sampling the coarser bonds. Within error bars, the data coincide with the exact result and the algorithm produces stable results free of the sign problem. Without MLB, the signal-to-noise ratio was practically zero for $`t>2`$.

B. Two-level system. For a symmetric two-state system, $`H=\frac{1}{2}\mathrm{\Delta }\sigma _x`$, the dynamics is controlled by tunneling. The spin correlation function $`\sigma _z(0)\sigma _z(t)`$ exhibits oscillations indicative of quantum coherence. Figure 2(b) shows MLB results for $`C(t)=\mathrm{Re}\sigma _z(0)\sigma _z(t)`$, Putting $`\mathrm{\Delta }=1`$, the tunneling oscillations have a period of $`2\pi `$. With $`P=64`$ for the maximum time $`t=64`$, only $`K=100`$ samples were used for sampling the coarser bonds. The data agree well with the exact result. Again the simulation is stable and free of the sign problem. Without MBL, the simulation failed for $`t>4`$.

C. Double-well potential. Next, we examine a double-well system with the quartic potential $`V(x)=x^2+\frac{1}{4}x^4`$. At low temperatures, interwell transfer occurs through tunneling motions on top of intrawell vibrations. These two effects combine to produce nontrivial structures in the position correlation function. At high temperatures, interwell transfer can also occur by classical barrier crossings. Figure 2(c) shows MLB results for $`C(t)=\mathrm{Re}x(0)x(t)`$. The slow oscillation corresponds to interwell tunneling, with a period of approximately 16. The higher-frequency motions are characteristic of intrawell oscillations. In this simulation, $`K=300`$ samples were used. The data reproduce the exact result well, capturing all the fine features of the oscillations. Again the calculation is stable and free of the sign problem, whereas a direct simulation failed for $`t>3`$.

D. Multidimensional tunneling system. As a final example, we consider a problem with three degrees of freedom, in which a particle in a double-well potential is bilinearly coupled to two harmonic oscillators. The quartic potential in the last example is used for the double-well, and the harmonic potential in the first example is used for both oscillators. The coupling constant between each oscillator and the tunneling particle is $`\alpha =1/2`$ in dimensionless units. For this example, we computed the correlation function $`C_s(t)`$ for the position operator of the tunneling particle. Direct application of MC sampling to $`C_s(t)`$ has generally been found unstable for $`t>\beta \mathrm{}/2`$ . In contrast, employing only moderate values of $`K=400`$ to 900, the MLB calculations allowed us to go up to $`t=10\beta \mathrm{}`$. (Notice that this $`K`$ is no larger than three, i.e. the number of dimensions, times the $`K`$ needed for one-dimensional systems.) Figure 2(d) shows MLB results for $`C_s^{}(t)=\mathrm{Re}C_s(t)`$. For the coupled system, the position correlations have lost the coherent oscillations and instead decay monotonically with time. Coupling to the medium clearly damps the coherence and tends to localize the tunneling particle.

The data presented here demonstrate that the MLB method holds substantial promise toward an exact and stable simulation method for real-time quantum dynamics computations of many-dimensional systems up to timescales of practical interest. Instead of the exponentially vanishing signal-to-noise ratio in a ordinary application of the Monte Carlo method to real-time path integral problems, the MLB method has a CPU requirement that grows only linearly with time. Moreover, the data we have so far seem to suggest that the memory requirement $`K`$ also grows only linearly with the dimensionality of the system.

This research has been supported by the National Science Foundation under grants CHE-9257094 and CHE-9528121, by the Sloan Foundation, the Dreyfus Foundation, and by the Volkswagen-Stiftung.