Rapid mixing of path integral Monte Carlo for 1D stoquastic Hamiltonians

TL;DR

This paper provides a rigorous analysis of the mixing time for path integral Monte Carlo applied to 1D stoquastic Hamiltonians, justifying its effectiveness for certain disordered quantum spin systems at specific temperatures.

Contribution

It offers the first formal proof of rapid mixing for PIMC in 1D disordered stoquastic models, including long-range and nearest-neighbor interactions.

Findings

Proves PIMC mixes rapidly at inverse temperatures logarithmic in system size.

Uses canonical paths method to analyze convergence.

Justifies PIMC's use for approximating partition functions in these models.

Abstract

Path integral quantum Monte Carlo (PIMC) is a method for estimating thermal equilibrium properties of stoquastic quantum spin systems by sampling from a classical Gibbs distribution using Markov chain Monte Carlo. The PIMC method has been widely used to study the physics of materials and for simulated quantum annealing, but these successful applications are rarely accompanied by formal proofs that the Markov chains underlying PIMC rapidly converge to the desired equilibrium distribution. In this work we analyze the mixing time of PIMC for 1D stoquastic Hamiltonians, including disordered transverse Ising models (TIM) with long-range algebraically decaying interactions as well as disordered XY spin chains with nearest-neighbor interactions. By bounding the convergence time to the equilibrium distribution we rigorously justify the use of PIMC to approximate partition functions and expectations of observables for these models at inverse temperatures that scale at most logarithmically with the number of qubits. The mixing time analysis is based on the canonical paths method applied to the single-site Metropolis Markov chain for the Gibbs distribution of 2D classical spin models with couplings related to the interactions in the quantum Hamiltonian. Since the system has strongly nonisotropic couplings that grow with system size, it does not fall into the known cases where 2D classical spin models are known to mix rapidly.

Equations

$$\mathcal{H}=\sum _{j=1}^n{H}_{j,j+1},$$

$$⟨z|H_{j,j+1}|z^{\prime }⟩\le 0\text{for all}z,z^{\prime }\in \{-1,1{\}}^n\text{with}z\ne z^{\prime }.$$

$${\rho }_{\beta }:=\frac{e^{-\beta \mathcal{H}}}{{\mathcal{Z}}_{\beta }},{\mathcal{Z}}_{\beta }:=\mathrm{tr}\left(e^{-\beta \mathcal{H}}\right),$$

$$\mathcal{H}=-\sum _{j=1}^n{\Gamma }_j{\sigma }_j^x+\sum _{j=1}^n{K}_j^z{\sigma }_j^z+\sum _{j,k=1}^n{K}_{jk}^{zz}{\sigma }_j^z{\sigma }_k^z,\Gamma :=\underset{j=1,...,n}{\min }{\Gamma }_j>0,$$

$$K_j^z\in [-1,1],\text{ and }|K_{jk}^{zz}|\le |i-j{|}^{-(2+\xi )}\text{ for }\xi >0.$$

$$H=-\sum _{j=1}^n{\Gamma }_j{\sigma }_j^x-\sum _{j=1}^n{K}_j^{xx}{\sigma }_j^x{\sigma }_{j+1}^x-\sum _{j=1}^n{K}_j^{yy}{\sigma }_j^y{\sigma }_{j+1}^y+\sum _{j=1}^n\left(K_j^{zz}{\sigma }_j^z{\sigma }_{j+1}^z+K_j^z{\sigma }_j^z\right),$$

$$K_j^{zz},K_j^z\in [-1,1],K_j^{xx}\ge 0,K_j^{yy}\in [-K_j^{xx},K_j^{xx}]\text{and}\Gamma :=\underset{j=1,...,n}{\min }{\Gamma }_j>0.$$

$$|{\tilde{\mathcal{Z}}}_{\beta }-{\mathcal{Z}}_{\beta }|$$

$$|\tilde{\mathcal{Q}}-\mathrm{tr}[\mathcal{Q}{\rho }_{\beta }]|$$

$$\mathrm{poly}(n,e^{\beta },{\Gamma }^{-1},{\delta }_{\text{mult}}^{-1},\log \left({\delta }_{\text{add}}^{-1}\right),\log ({\delta }_{\text{fail}}^{-1}))$$

$$\mathrm{poly}(2^{\beta {\left(\log \frac{n}{{\delta }_{\text{mult}}}\right)}^2},\log (1/{\delta }_{\text{add}}),\log (1/{\delta }_{\text{fail}}))$$

$$E_{\beta }(\mathbf{z})=\sum _{(i,j)\in \Lambda }\left(\frac{K_{j,j+1}^{zz}}{L}{z}_{i,j}{z}_{i,j+1}+\frac{K_j^z}{L}{z}_{i,j}-{\beta }^{-1}{J}_i{z}_{i,j}{z}_{i+1,j}\right),z_{L+1,j}:=z_{1,j},$$

$$\sum _{j=1}^{n/2}\sum _{k=\frac{n}{2}+1}^n|j-k{|}^{-(2+\xi )}=\mathcal{O}({\xi }^{-2}+n^{-\xi }),$$

$$E_{\beta }(\mathbf{z}):=\sum _{(i,j)\in \Lambda }\left(\frac{K_j^{zz}}{L}{z}_{i,j}{z}_{i,j+1}+\frac{K_j^z}{L}{z}_{i,j}-{\beta }^{-1}{h}_{i,j}\left(z_{i,j},z_{i,j+1},z_{i+1,j},z_{i+1,j+1}\right)\right),$$

$$H\to H-\Gamma \sum _{i=1}^n{\sigma }_i^x$$

$$Z_{\beta ,L}:=\mathrm{tr}\left[{\left(e^{\frac{A}{2L}}{e}^{\frac{B}{L}}{e}^{\frac{A}{2L}}{e}^{\frac{C}{L}}\right)}^L\right].$$

$$e^{-\delta }{\mathcal{Z}}_{\beta }\le Z_{\beta ,L}\le e^{\delta }{\mathcal{Z}}_{\beta },$$

$$Z_{\beta ,L}=\mathrm{tr}\left[{\left(e^{\frac{A}{2L}}{e}^{\frac{B}{L}}{e}^{\frac{A}{2L}}\right)}^L\right]=\mathrm{tr}\left[{\left(e^{\frac{A}{L}}{e}^{\frac{B}{L}}\right)}^L\right].$$

$$Z_{\beta ,L}$$

$$=\sum _{z_1,\dots ,z_L}\prod _{i=1}^L⟨z_i|e^{\frac{A}{L}}{e}^{\frac{B}{L}}|z_{i+1}⟩,z_{L+1}:=z_1$$

$$=\sum _{z_1,\dots ,z_L}\prod _{i=1}^L{e}^{\frac{A(z_i)}{L}}⟨z_i|e^{\frac{B}{L}}|z_{i+1}⟩.$$

$$⟨z_i|e^{\frac{B}{L}}|z_{i+1}⟩=⟨z_i|e^{-\frac{\beta \Gamma }{L}{\sum }_{j=1}^n{\sigma }_j^x}|z_{i+1}⟩=\prod _{j=1}^n⟨z_{i,j}|e^{-\frac{\beta \Gamma }{L}{\sigma }_j^x}|z_{i+1,j}⟩.$$

$$A=-\beta \mathrm{diag}(H),B=-\beta \sum _{i\in \text{Odd}}\left[H_i-\mathrm{diag}(H_i)\right],C=-\beta \sum _{i\in \text{Even}}\left[H_i-\mathrm{diag}(H_i)\right],$$

$$Z_{\beta ,L}$$

$$=\sum _{z_1,\dots ,z_{2L}}\prod _{i\in \text{Odd}}^{2L-1}⟨z_i|e^{\frac{A}{2L}}{e}^{\frac{B}{L}}|z_{i+1}⟩⟨z_{i+1}|e^{\frac{A}{2L}}{e}^{\frac{C}{L}}|z_{i+2}⟩$$

$$=\sum _{z_1,\dots ,z_{2L}}{e}^{{\sum }_{i=1}^{2L}\frac{A(z_i)}{2L}}\prod _{i\in \text{Odd}}^{2L-1}⟨z_i|e^{\frac{B}{L}}|z_{i+1}⟩⟨z_{i+1}|e^{\frac{C}{L}}|z_{i+2}⟩,$$

$$\pi (\mathbf{z}):=\frac{1}{Z_{\beta ,L}}{e}^{{\sum }_{i=1}^L\frac{A(z_i)}{L}}\prod _{i=1}^L⟨z_i|e^{\frac{B_i}{L}}|z_{i+1}⟩.$$

$${\mathcal{Z}}_{\beta }(\zeta ):=\mathrm{tr}[e^{-\beta \mathcal{H}+\zeta \mathcal{O}}],$$

$${⟨\mathcal{O}⟩}_{\beta }=\frac{1}{{\mathcal{Z}}_{\beta }(0)}\frac{\partial {\mathcal{Z}}_{\beta }}{\partial \zeta }(0).$$

$$e^{-\delta }\le \frac{{\hat{\mathcal{Z}}}_{\beta }(0)}{{\mathcal{Z}}_{\beta }(0)},\frac{{\hat{\mathcal{Z}}}_{\beta }(\zeta )}{{\mathcal{Z}}_{\beta }(\zeta )}\le e^{\delta },$$

$$\widehat{{⟨\mathcal{O}⟩}_{\beta }}:=\frac{{\hat{\mathcal{Z}}}_{\beta }(\zeta )-{\hat{\mathcal{Z}}}_{\beta }(0)}{\zeta {\hat{\mathcal{Z}}}_{\beta }(0)}.$$

$$|\widehat{{⟨\mathcal{O}⟩}_{\beta }}-{⟨\mathcal{O}⟩}_{\beta }|\le 2\sqrt{\delta }\Vert \mathcal{O}\Vert$$

$$\left|{\mathcal{Z}}_{\beta }(\zeta )-{\mathcal{Z}}_{\beta }(0)-\zeta \frac{\partial {\mathcal{Z}}_{\beta }}{\partial \zeta }(0)\right|$$

$$\le {\zeta }^2\Vert \mathcal{O}{\Vert }^2\underset{{\zeta }^{\prime }\in [0,\zeta ]}{\max }{\mathcal{Z}}_{\beta }({\zeta }^{\prime })$$

Full text

In this work we analyze the mixing time of PIMC for 1D stoquastic Hamiltonians, including disordered transverse Ising models (TIM) with long-range algebraically decaying interactions as well as disordered XY spin chains with nearest-neighbor interactions. By bounding the convergence time to the equilibrium distribution we rigorously justify the use of PIMC to approximate partition functions and expectations of observables for these models at inverse temperatures that scale at most logarithmically with the number of qubits.

The mixing time analysis is based on the canonical paths method applied to the single-site Metropolis Markov chain for the Gibbs distribution of 2D classical spin models with couplings related to the interactions in the quantum Hamiltonian. Since the system has strongly nonisotropic couplings that grow with system size, it does not fall into the known cases where 2D classical spin models are known to mix rapidly.

1 Introduction

The application of Markov chain Monte Carlo methods to thermal states of a limited class of quantum systems was first proposed by Suzuki, Miyashita and Kuroda 40 years ago [43]. Their method applied to Hamiltonians with real and nonpositive off-diagonal matrix elements in a particular basis, a property that was called “being free of the sign problem” for many years until the term “stoquastic” was eventually adopted [12] to emphasize a connection with nonnegative matrices. Since then there have been a growing number of works on rigorous classical simulations for various stoquastic systems [11, 8, 15, 9]. Many quantum systems of physical and computational interest are stoquastic, including spinless particles moving on arbitrary graphs with position dependent interactions, as well as generalized transverse Ising models, which are particularly notable for their use in quantum annealing [31, 19, 1].

Statement of results.

A general $n$-qubit 1D stoquastic Hamiltonian with nearest-neighbor interactions has the form,

[TABLE]

where each $H_{j,j+1}$ acts non-trivially on (at most) the sites $j,j+1$ (with $n+1$ identified with 1) and each $\|H_{j,j+1}\|\leq 1$. Define the “computational basis” $\{\left|1\right\rangle,\left|-1\right\rangle\}$ for each qubit to be the eigenstates of the $\sigma^{z}$ operator. We assume $\mathcal{H}$ is stoquastic in the computational basis [12, 11], which means that in this basis each $H_{j,j+1}$ is a real symmetric matrix with nonpositive off-diagonal entries,

[TABLE]

For a given inverse temperature $\beta$ we consider the quantum partition function $\mathcal{Z}_{\beta}$ and the quantum Gibbs state $\rho_{\beta}$,

[TABLE]

and also the expectation value of observables $\mathcal{Q}$ in the thermal state, $\langle\mathcal{Q}\rangle_{\beta}:=\operatorname{tr}[\mathcal{Q}\rho_{\beta}]$. For the special class of systems called 1D generalized transverse Ising models (TIM) we go beyond (1) and allow systems with long-range interactions that decay sufficiently quickly (faster than inverse square) with distance,

[TABLE]

The form (5) simplifies both the presentation of the PIMC method and our proof that the mixing time of PIMC for these models is at most $\textrm{poly}(n,e^{\beta},\Gamma^{-1})$. We also show the mixing time of PIMC is at most $\operatorname{poly}(n,e^{\beta},\Gamma^{-1})$ for generalized TIM with $\sigma^{x}\sigma^{x}$ and $\sigma^{y}\sigma^{y}$ interactions,

[TABLE]

The mixing time bounds can be used to construct an FPRAS (fully polynomial-time randomized approximation scheme) for the partition function of models of the form (5) and (6) and inverse temperatures $\beta={\cal O}(\log n)$. In addition to the partition function there is a similar FPRAS for the expectation value of observables that are sparse and row-computable (in the computational basis) in the thermal state $\rho_{\beta}$ for $\beta={\cal O}(\log n)$. For general 1D stoquastic Hamiltonians (1) we prove that the PIMC mixing time is at most quasi-polynomial in the system size, which yields a quasi-polynomial time approximation algorithm for the partition function and for the approximation of observables.

Theorem 1**.**

There is an algorithm which takes as input an inverse temperature $\beta$ and an $n$-qubit 1D stoquastic Hamiltonian $H$ and outputs an estimate $\tilde{Z}_{\beta}$ of the quantum partition function $\mathcal{Z}_{\beta}$.

If also given an observable $\mathcal{Q}\succeq 0$ that is sparse and row-computable in the computational basis, the algorithm outputs an estimate $\tilde{\mathcal{Q}}$ of its expectation value in the state $\rho_{\beta}$.

These estimates satisfy

[TABLE]

with probability $\geq 1-\delta_{\text{fail}}$. The algorithm runs in time at most

[TABLE]

for models of the form (5) and (6), and in time at most

[TABLE]

for models of the form (1).

Relation to previous work.

A variety of quantum Monte Carlo methods have been devised for estimating equilibrium properties of stoquastic Hamiltonians, which can be broadly divided into diffusion methods [21, 40, 27] and path integral methods [43, 42, 36, 2, 39]. Diffusion Monte Carlo methods use Hamiltonian matrix elements to define transition rates for a substochastic random process, along with birth and death population dynamics or importance sampling to control statistical fluctuations. Examples of proven efficient convergence for diffusion Monte Carlo methods include simulating frustration-free stoquastic adiabatic computation [11] and simulating the “go with the winners” method [26]. Disadvantages of diffusion methods include quantum amplitude vs quantum probability obstructions for some variants [27, 13], as well as the need for importance sampling with trial wave functions that may in general not be available [8].

PIMC methods, in contrast, are based on Markov chains defined on an enlarged state space consisting of many coupled copies of a high-temperature classical system along a so-called “imaginary-time” direction [43]. PIMC has been used to give an FPRAS for the ferromagnetic case of TIM [8] on general interaction graphs by a reduction to a classic result on sampling the Gibbs distribution of classical ferromagnetic Ising models [28]. Also, the Glauber dynamics for an infinite dimensional variant of PIMC has been shown to have optimal temporal mixing for ferromagnetic TIM at high temperature [14] and for ferromagnetic TIM on regular trees [35]. Here “ferromagnetic” means of the form in (5) but with each $K_{j}^{z}=0$ and each $K_{j}^{zz}\leq 0$, i.e. lower energy is achieved by aligning adjacent spins. Since such interactions cannot be frustrated, these models are considered easier to analyze. Disadvantages of PIMC methods include the significant overhead created by adding the imaginary-time direction, and also the possibility of slow mixing due to topological obstructions [3, 23]. The general presence of obstructions to rapid mixing for the standard versions of both diffusion and path integral methods also motivates the continuing development of new algorithms for simulating stoquastic Hamiltonians [2, 9].

To our knowledge, the simulations in this work are the first provably efficient Monte Carlo methods for a large class of generalized TIM systems, which differ from ferromagnetic TIM by allowing couplings of both signs between neighboring spins, as well as local fields of both signs. A generalized TIM in 2D or higher dimensions can encode NP-complete constraint satisfaction problems, which is the basis for the use of these models in quantum annealing [31, 19, 1]. Besides their relative simplicity from the standpoint of experimental implementation, generalized TIM on interaction graphs of degree 3 are also known to be universal for stoquastic quantum annealing [10, 17]. The application of Monte Carlo methods to stoquastic quantum annealing is called simulated quantum annealing, and a significant open question is whether some variant of simulated quantum annealing can efficiently simulate stoquastic quantum annealing. This open question has motivated several recent proposals for non-stoquastic quantum annealing [16, 25, 46, 41]. Some recent progress has been made by analyzing specific models for which both quantum annealing and simulated quantum annealing can be exponentially faster than classical simulated annealing [20, 38, 32, 37, 6, 7, 15, 30].

Another class of simulation methods for 1D quantum systems are those based on matrix product states [45]. These methods do not require the Hamiltonian to be stoquastic and have been applied in rigorous classical simulations of systems with limited entanglement [22, 33, 4]. Their runtime scales exponentially with spectral gap, while our runtime bound for PIMC scales exponentially with the inverse temperature, which is qualitatively similar. These kinds of simulations are not expected to obsolete PIMC in general, however, as the latter is used in practice to simulate highly entangled stoquastic systems. Furthermore, the version of PIMC that we consider is closely related to the practical implementations of PIMC, and so this work can be seen as a rigorous analysis of a heuristic approximation algorithm that was originally developed and applied with great success by the physics community.

Organization of the remaining sections.

Section 2 reviews the PIMC method and sketches the proof of our mixing time analysis. Section 3.1 defines the Suzuki-Trotter approximation to the partition function which is used to map 1D stoquastic Hamiltonians onto 2D classical spin systems. Section 3.2 and Section 3.3 describe the use of samples from the Gibbs distribution of the 2D classical system to approximate quantum observables and the quantum partition function. Section 3.4 proves a concentration theorem for PIMC that justifies the analysis of a restricted state space in later sections. Section 4.1 reviews the canonical paths method that we use to analyze the PIMC mixing time. This method is then applied to generalized TIM in Section 4.3, generalized TIM with XX interactions in Section 4.4, and general 1D stoquastic Hamiltonians in Section 4.5.

2 Overview

In the PIMC method the Hamiltonian (1) is related by Suzuki’s quantum-to-classical mapping [42, 44] to a system of $\{\pm 1\}$ classical spins on a 2D lattice $\Lambda=\{1,\ldots,L\}\times\{1,\ldots,n\}$. The new site index $\{1,\ldots,L\}$ is known in physics as the “imaginary-time” direction, and it is often useful to think of the set of spin configurations $\Omega:=\{-1,1\}^{n\times L}$ as consisting of “time slices” along the $L$ direction, i.e. for $\mathbf{z}\in\Omega$ we write $\mathbf{z}:=(z_{1},\ldots,z_{L})$ with each $z_{i}\in\{-1,1\}^{n}$. The mapping introduces a local energy function $E_{\beta}:\Omega\rightarrow\mathbb{R}$ on the classical spin configurations in such a way that an estimate of the classical partition function $Z:=\sum_{\mathbf{z}\in\Omega}e^{-\beta E_{\beta}(\mathbf{z})}$ can be used to estimate the quantum partition function $\mathcal{Z}$, and so that samples from the Gibbs distribution $\pi(\mathbf{z}):=e^{-\beta E_{\beta}(\mathbf{z})}/Z$ can be used to estimate expectation values in the quantum Gibbs state. For the nearest-neighbor version of the transverse Ising models in (5) the energy function $E_{\beta}$ is

[TABLE]

where $z_{i,j}$ denotes the spin in configuration $\mathbf{z}$ at the site $(i,j)\in\Lambda$, and $J_{i}:=\frac{1}{2}\log\coth\left(\frac{\beta\Gamma_{i}}{L}\right)$ (see Section 3.1). The couplings along the imaginary-time direction are ferromagnetic, but the couplings along the spatial direction can have varying signs as well as local fields. Our choices of parameters are such that $L=\operatorname{poly}(n)$, and so the general case considers $J_{i}$ with magnitude up to $\mathcal{O}(\log n)$. Therefore the couplings of the 2D model are highly anisotropic, as in figure 1.

The PIMC method proceeds by attempting to sample from $\pi$ and approximate $Z$ using the Markov chain Monte Carlo method with single-site Metropolis transition probabilities. With this in mind, the primary question is to determine the mixing time $\tau$ for this Markov chain, which is the minimum number of steps of this Markov chain needed to sample from a distribution that is close to $\pi$. We will use the method of canonical paths that was originally introduced by Jerrum and Sinclair (reviewed in [29, 34]) in order to bound the mixing time. The idea of the method is that for every pair of states $\mathbf{x},\mathbf{y}\in\Omega$ the Markov chain $P$ needs to route an amount of “traffic” equal to $\pi(\mathbf{x})\pi(\mathbf{y})$ along some path $\gamma_{\mathbf{x},\mathbf{y}}$ which connects $\mathbf{x}$ to $\mathbf{y}$ using the transitions of the Markov chain. Define the congestion across a particular transition $e=(\mathbf{z},\mathbf{z}^{\prime})$ to be the ratio of the total traffic through it (i.e. $\sum_{\mathbf{x},\mathbf{y}:e\in\gamma_{\mathbf{x},\mathbf{y}}}\pi(\mathbf{x})\pi(\mathbf{y})$) to the amount of probability flow across it at equilibrium, which equals $\pi(\mathbf{z})P(\mathbf{z},\mathbf{z}^{\prime})$. The mixing time can then be upper bounded in terms of the maximum congestion through any transition, using standard arguments which we review in Section 4.1.

Our problem now reduces to finding good sets of paths between every pair of configurations $\mathbf{x},\mathbf{y}\in\Omega$. We use a set of paths inspired by those used to analyze the (unweighted) random walk on the boolean hypercube [29]. Specifically we will replace the entries of $\mathbf{x}$ with entries of $\mathbf{y}$ one at a time. A crucial requirement is that the energies of the intermediate states do not get too large, as this would create a large congestion. For the model (9) the paths first update the classical spins corresponding to the site of qubit 1 in consecutive imaginary-time order, and then proceed to do the same for the sites of qubits 2 through $n$. For each qubit $j\in\{1,\ldots,n\}$ the line of sites $\{(1,j),(2,j),\ldots,(L,j)\}$ is called the worldline of the $j$-th qubit, and so these canonical paths have the form of updating each qubit worldline consecutively, always finishing one worldline before starting the next.

An example of an intermediate step along this path is depicted in Fig. 2. Notice that intermediate steps have only introduced a constant number of broken bonds between spins that are neighbors in the imaginary time direction, each contributing an amount of energy that is ${\cal O}(\beta^{-1}\log n)$, and at most $L$ broken bonds between spins that are neighbors in the space direction, each contributing energy ${\cal O}(1/L)$. The relatively low energy of the intermediate configurations along the path is the key to the formal proof of rapid mixing for PIMC applied to 1D generalized TIM in Section 4.3. For 1D generalized TIM with long-range interactions, the energy of the intermediate configuration $\mathbf{z}$ will include contributions from spins that are arbitrarily far away from the cut. Using the assumption on the rate of decay of the magnitude of these interactions, the total contribution to the energy is upper bounded by

[TABLE]

where the bound is proven by viewing the summation as a right Riemann sum of the corresponding integral, which yields an upper bound because the integrand is a monotonically decreasing function. Since we regard $\xi$ as a constant this contribution to the energy is $\mathcal{O}(1)$.

For Hamiltonians of the form (6) the 2-local off-diagonal terms give rise to a classical energy function with interactions involving four spins at a time,

[TABLE]

The $h_{i,j}$ terms (which are implicitly defined by the Gibbs distribution (22)) correspond to strong ferromagnetic interactions between neighboring pairs of spins in neighboring time slices, and violating them increases the energy by an amount that is $\mathcal{O}(\log n)$. If $z_{i,j}\neq z_{i+1,j}$ and $z_{i,j+1}\neq z_{i+1,j+1}$ then we say there is a 2-local jump between time slice $i$ and time slice $i+1$, because the resulting energy contribution from the 4-body interaction will depend on the magnitude of the 2-local off-diagonal quantum terms i.e. $K^{xx}_{j}$ and $K^{yy}_{j}$.

The canonical paths we have described so far will lead to intermediate configurations that potentially create or destroy a non-constant number of 2-local jumps across the domain wall pictured in figure 2. If $K^{xx}_{j}$ is larger than $K_{j}^{x}$, then configurations $\mathbf{x}$ and $\mathbf{y}$ which contain a non-constant number of 2-local jumps will result in the canonical paths described so far overloading some edges with large congestion. To overcome this we include additional path segments along the canonical paths that remove 2-local jumps from the worldlines near the domain wall in figure 2. This greatly increases the number of paths passing through each edge, but we are still able to obtain a $\operatorname{poly}(n,e^{\beta})$ upper bound on the congestion by noticing that most of these paths have endpoints with very small stationary probabilities.

For Hamiltonians of the form (1) it is possible that the off-diagonal terms in (1) do not allow for the PIMC Markov chain with single-site transitions to be ergodic. In this situation some configurations $\mathbf{z}\in\Omega$ will have $\pi(\mathbf{z})=0$, which creates a new obstacle to the canonical paths we have described so far. To restore ergodicity under single-site updates we add a fictitious transverse field to the Hamiltonian with a weak coupling that goes to zero with system size:

[TABLE]

with $\Gamma\leq\delta_{\text{mult}}/n$ so that the subsequent estimate for the partition function will not be further than the allowed error tolerance away from its intended value.

Once the fictitious field is added it once again makes sense to attempt to use the same canonical paths that were applied to systems of the form (6). But now the fact that the off-diagonal matrix elements of the Hamiltonian are no longer symmetric under the global operation of flipping all spins means that removing a single (1-local or 2-local) jump from a given worldline can decrease the stationary probability of that configuration by a factor scaling exponentially in the number of jumps remaining in that worldline. To overcome this problem we show in Section 3.4 that typical configurations distributed according to $\pi$ will have at most $\mathcal{O}(\beta\log n)$ jumps in any particular worldline. This means that the worst-case intermediate configurations just described may have a stationary probability reduced by $n^{O(\beta\log n)}$, and this factor bottlenecks the mixing time estimate. Using these facts we describe canonical paths in Section 4.5 within a restricted state space $\Omega^{*}$ in which every configuration has ${\cal O}(\beta\log n)$ jumps per worldline, and by either restricting the PIMC Markov chain to $\Omega^{*}$ or applying the leaky random walk conductance bound [15] we obtain a bound on the mixing time within the subset $\Omega^{*}$ that establishes theorem 1.

3 Path integral Monte Carlo

Sections 3.1 through 3.3 rigorously present the underpinnings of the standard PIMC method: the quantum-to-classical mapping, and its application to approximating expectation values of observables as well as the partition function. Section 3.4 proves a general concentration of measure bound which shows that the number of jumps along any worldline in PIMC is $\mathcal{O}(\beta\log n)$ with high probability. Note that this result applies to general stoquastic Hamiltonians (not just the 1D case). This result is used to analyze a restricted state space in Section 4.5.

3.1 Quantum-to-classical mapping

The first step of the quantum-to-classical mapping is to split $H$ as in (1) into a sum of commuting Hamiltonians. For each local term $H_{i}$ we construct local operators $A_{i},B_{i},C_{i}$ such that $-\beta H_{i}=A_{i}+B_{i}+C_{i}$ and $[A_{i},A_{j}]=[B_{i},B_{j}]=[C_{i},C_{j}]=0$ for all $i,j$. Let $A=\sum_{i}A_{i}$, $B=\sum_{i}B_{i}$, $C=\sum_{i}C_{i}$, and hereafter let $A_{i}$ match the diagonal part of $H_{i}$ in the computational basis, while $B_{i}$ and $C_{i}$ contain the off-diagonal terms of $H_{i}$ when $i$ is odd or even, respectively. For positive integers $L$ define the Suzuki-Trotter approximation of the partition function $Z=\operatorname{tr}\left[e^{A+B+C}\right]$,

[TABLE]

The fact that $\lim_{L\rightarrow\infty}Z_{\beta,L}=\mathcal{Z}_{\beta}$ is the content of the 1875 Lie product formula, while for finite $L$ the Suzuki-Trotter approximation is related to the quantum partition function (3) by

[TABLE]

for some $\delta$ that is ${\cal O}\left(\|H\|^{3}/L^{3}\right)$. This error bound follows from two facts about the matrix exponential: a variant of the usual Suzuki-Trotter expansion, and a continuity bound. Both appear in [8] for transverse Ising models, and the proof for three commuting layers $A,B,C$ is similar and so we omit the details.

Before considering the general case of $H$ as in (1) it is useful to consider the simplified case of (12) for generalized TIM of the form (5) by taking $A=-\beta\operatorname{diag}(H)$, $B=-\beta\left(H-\operatorname{diag}(H)\right)$, and $C=0$ so that the Suzuki-Trotter approximation reduces to

[TABLE]

For transverse Ising models, the next step is to expand the Suzuki-Trotter approximation (14) to the partition function by inserting several complete sets of states,

[TABLE]

where in the last step we have used the fact that $A=-\beta\operatorname{diag}(H)$ is diagonal in the computational basis and defined $A(z_{i}):=\langle z_{i}|A|z_{i}\rangle$. The Gibbs distribution of the classical 2D spin system is defined by equating $e^{-\beta E_{\beta}(z_{1},\ldots,z_{L})}$ with the summand of (17). The closed form (9) is obtained applying the identity $e^{\omega\sigma^{x}}=\cosh(\omega)I+\sinh(\omega)\sigma^{x}$ to

[TABLE]

For a general 1D stoquastic Hamiltonian of the form in (1) we will apply (12) with

[TABLE]

Expanding (12) with $2L$ complete sets of states,

[TABLE]

where periodic boundary conditions $z_{2L+1}:=z_{1}$ are taken in (20). At this point the essential difference from the transverse Ising case is clear: the off-diagonal part of the Hamiltonian is split into two commuting layers, these layers alternate along the imaginary-time direction, and this causes the coupling between classical spins at neighboring time slices to alternate as well. For ease of notation define $B_{i}$ to be $2B$ if $i$ is odd, and $2C$ if $i$ is even, and rescale $L\rightarrow L/2$ to obtain a distribution $\pi$ from (21),

[TABLE]

3.2 Approximating observables

In this section we consider an observable ${\cal O}$ for which we would like to estimate $\left\langle{\cal O}\right\rangle_{\beta}:=\operatorname{tr}[{\cal O}\rho_{\beta}]$. We justify this by demonstrating closeness of $\left\langle{\cal O}\right\rangle_{\beta}$ to a value $\langle O\rangle_{\beta,L}$ that is defined in terms of the Suzuki-Trotter approximation in equation (34). The expectation $\left\langle{\cal O}\right\rangle_{\beta}$ can be expressed as the derivative of a partition function by defining

[TABLE]

and using the fact that

[TABLE]

The next proposition shows that $\langle{\cal O}\rangle_{\beta}$ can be approximated by a finite difference.

Proposition 2**.**

Suppose that $\hat{\mathcal{Z}}_{\beta}(0)$ and $\hat{\mathcal{Z}}_{\beta}(\zeta)$ satisfy

[TABLE]

for $\zeta=\sqrt{3\delta}/\|{\cal O}\|$. Then the estimator

[TABLE]

satisfies the bound

[TABLE]

Proof.

We use first Taylor’s theorem, then the bound $\operatorname{tr}AB\leq\|A\|\|B\|_{1}$ to obtain:

[TABLE]

Dividing both sides of (27) by $\zeta\mathcal{Z}_{\beta}(0)$ we have

[TABLE]

Then

[TABLE]

In the last step we have used the assumption that $\delta\leq 1/21$ to simplify the numerical constants. A similar argument shows that $\widehat{\left\langle{\cal O}\right\rangle_{\beta}}-{\left\langle{\cal O}\right\rangle_{\beta}}\geq-2\sqrt{\delta}\|{\cal O}\|$, thus completing the proof. ∎

By (25), (13), and the triangle inequality, the expectation $\langle{\cal O}\rangle_{\beta}$ can be approximated by

[TABLE]

where $Z_{\beta,L}(\zeta)$ is a Suzuki-Trotter approximation to (23). If ${\cal O}$ is diagonal use the form in Section 3.1 with $A\rightarrow A+\zeta{\cal O}$,

[TABLE]

Evaluating the derivative at $\zeta=0$,

[TABLE]

where the final expression follows by the symmetry of $\pi$ under cyclic permutations of the time slices. For off-diagonal ${\cal O}$,

[TABLE]

Evaluating the derivative at $\zeta=0$,

[TABLE]

where the last step follows by symmetry of $\pi$ under even cyclic permutations of the time slices. Since the summand of (42) is efficiently computable and has bounded magnitude we can estimate it using a standard Monte Carlo procedure (see lemma 3) once one has an efficient meanns of obtaining independent samples from $\pi$.

3.3 Approximating partition functions

In Section 4 we will prove that for any tolerance $\delta>0$ the PIMC method can efficiently output an element of $\Omega$ sampled from a distribution $\tilde{\pi}$ which is within $|\tilde{\pi}-\pi|\leq\delta$. Here we explain how these samples can be used to approximate the partition function with an application of the standard telescoping product technique [28, 5]111The techniques used in this section to estimate the partition function (or equivalently the free energy) are necessary because the entropy term in the free energy cannot be directly estimated as an observable. If one is only interested in the estimation of specific observables then it suffices to use the end result(42).. We express the partition function $Z(\beta)$ as a product of terms involving the partition function at higher temperatures,

[TABLE]

where the temperature schedule $\beta_{0}<\beta_{1}<\ldots<\beta_{k}$ with $\beta_{0}=0$ and $\beta_{k}=\beta$ is chosen so that each ratio in the telescoping product is bounded. Defining $w_{\beta_{i}}(\mathbf{z}):=Z(\beta_{i})\pi_{\beta_{i}}(\mathbf{z})$, note that

[TABLE]

At infinite temperature the quantum partition function is $\mathcal{Z}(0)=2^{n}$. Further, a uniform schedule with $k=\mathcal{O}(\beta\|H\|\log(\beta\|H\|))$ ensures that

[TABLE]

Note that in infinite temperature limit $\beta=0$ the classical effective systems (9) and (10) have infinite coupling strength between spins in the imaginary-time direction, and so the Markov chain we analyze in Section 4 will not be ergodic. This is ununsual from the perspective of classical spin systems, but it happens in this case because of the temperature dependence of the couplings. To deal with this the first expectation $\left\langle w_{\beta_{1}}/w_{\beta_{0}}\right\rangle_{\pi_{\beta_{0}}}$ in the telescoping product can be computed using the fact that $\pi_{\beta_{0}}$ is equivalent to the uniform distribution on the subset $\Omega_{\textnormal{classical}}\subseteq\Omega$ consisting of configurations $(z_{1},\ldots,z_{L})$ which have $z_{i}=z_{j}$ for all $i,j=1,\ldots,L$.

Finally, in order to bound the number of samples needed to compute the expectation values in the telescoping product we will make use of the Hoeffding bound.

Lemma 3** ([24]).**

Let $X_{1},\ldots,X_{t}$ be independent random variables satisfying $|X_{i}|\leq 1$ and $\operatorname*{\mathbb{E}}[X_{i}]=\bar{X}$. Then

[TABLE]

Define $M(z)=e^{-2\delta}w_{\beta_{i}}(z)/w_{\beta_{i+1}}(z)$, so that $|M(z)|\leq 1$. Observe that $|\langle M\rangle_{\pi}-\langle M\rangle_{\tilde{\pi}}|\leq\|\pi-\tilde{\pi}\|_{1}\|M\|\leq\delta$. We will estimate $\langle M\rangle_{\tilde{\pi}}$ by drawing $t$ samples from $\tilde{\pi}$, which we call $\mathbf{z}^{(1)},\ldots,\mathbf{z}^{(t)}$. Our estimator will simply be the sample mean, namely

[TABLE]

By Lemma 3 we have

[TABLE]

We conclude that taking $t=O(1/\delta^{2})$ samples yields an $O(\delta)$ approximation to $\langle M\rangle_{\tilde{\pi}}$ with high probability. Therefore a polynomial number of samples will suffice to estimate each term of the telescoping product.

3.4 Restricting the PIMC state space

In this section we justify restricting PIMC to a subset of the state space that limits the number of jumps in any particular worldline to be ${\cal O}(\beta\log n)$ (this restricted state space will be used for the mixing time analysis in Section 4.5). For a configuration $\mathbf{z}=(z_{1},...,z_{L})$ the number of jumps in the worldline of the $j$-th qubit can be expressed in terms of the Pauli operator $X_{j}$,

[TABLE]

The strategy is to bound the $k$-th moment of the random variable $d_{j}$ with respect to the QMC stationary distribution $\pi$ and use the fact that

[TABLE]

These moments satisfy

[TABLE]

where $A,B_{i}$ are defined in Section 3.1. To merge the two products in the expression above into a simpler trace expression, we want to show the following inequality for any neighboring time slices $z_{i}$ and $z_{i+1}$,

[TABLE]

where we assume the local terms are normalized so that for some $\Gamma>0$

[TABLE]

If $\langle z_{i}|X_{j}|z_{i+1}\rangle=0$ then (50) is satisfied. Otherwise if $\langle z_{i}|X_{j}|z_{i+1}\rangle=1$ then the numerator on the RHS of (50) satisfies

[TABLE]

To upper bound the denominator of the RHS of (50) , the fact that $B_{i}$ is a commuting Hamiltonian implies

[TABLE]

where in the last step we have used $\|H_{j}\|\leq 1$, ${\frac{\beta}{L}}\leq\ln(2)$ and the convexity of $e^{x}$. Using $\ln(2)/2\geq 1/4$ we have established (50).

To apply this result to (49) we will use the permutation symmetry of the $i_{1},...,i_{k}$ variables to restrict to a sum over ordered $k$-tuples $(i_{1},...,i_{k})$ which satisfy $i_{1}\leq i_{2}\leq\ldots i_{k}$.

[TABLE]

Defining $t_{r}=(i_{r}-i_{r-1})\beta/L$ (as well as $t_{1}=r_{1}\beta/L,t_{k}=(L-i_{k})\beta/L$) this expression can be seen as the trace of a product of operators,

[TABLE]

where

[TABLE]

and $L_{i}$,$R_{t_{i}}$ can each denote one of several products of operators depending on $t_{i}$, which is determined by where the $X_{j}$ operators disrupt the sequence,

[TABLE]

The purpose of this decomposition is that each $M_{t_{i}}$ is a PSD operator raised to an positive integer power, and the rest of the Matrices $L_{t_{i}},R_{t_{i}},X_{j}$ are all either PSD operators or are explicit products of 2,3, or 4 PSD operators. Now the traces in (54) can be bounded using classic inequality due to Ky Fan [18, Theorem 1]. Let $P_{1}...P_{N}$ satisfy $P_{i}\succeq 0$ for each $i$, and let $\sigma(P_{i})$ be a diagonal matrix with the eigenvalues of $P_{i}$ along the diagonal in descending order. Then

[TABLE]

For our purposes define the notation $\sigma(R_{t_{i}})$ to be the product of diagonal matrices of eigenvalues of the operators that make up $R_{t_{i}}$ e.g. if $R_{t_{i}}=e^{A/L}e^{B_{\textrm{even}}/L}$ then $\sigma(R_{t_{i}})=\sigma\left(e^{A/L}\right)\sigma\left(e^{B_{\textrm{even}}/L}\right)$. Furthermore we may take $B_{\textrm{even}},B_{\textrm{odd}}\succeq 0$ without loss of generality (by shifting the overall Hamiltonian by a multiple of the identity), so that $\mathbbm{1}\succeq\sigma(R_{t_{i}}),\sigma(L_{t_{i}})$ for each $t_{i}$. Therefore each trace in (54) satisfies

[TABLE]

where the fact that $\mathbbm{1}\succeq\sigma(R_{t_{i}}),\sigma(L_{t_{i}})$ is used to go from (55) to (56). Now for $k\ll L$ we have $Z_{\beta,L}\approx Z_{\beta,L-2k}$, so this implies that $\mathbbm{E}_{\pi}d^{k}_{j}\leq k!\beta^{k}$. Therefore taking $k=c\log n$ for some $c>0$ and using $k!\leq(k/e)^{k}$ yields

[TABLE]

4 Proofs of rapid mixing

Having established the correctness of the PIMC method in section 3 we now procede to analyze the mixing time of the associated Markov chains. The mixing time analysis is divided into cases. In section 4.1 we review the main techniques that are common to each case in our analysis. Section 4.3 contains the mixing analysis for transverse Ising models, while sections 4.4 and 4.5 analyze Hamiltonians of the form (6) and (1) , respectively.

4.1 Markov chains, mixing times, and canonical paths

Let $P$ be the transition matrix of a reversible Markov chain defined on a state space $\Omega$ with stationary distribution $\pi$. Let $\mathbf{z}$ be the state of the chain at time $t=0$, and define $P^{t}(\mathbf{z},\cdot)$ to be the distribution of the state of the chain at time t. The variation distance of the chain at time $t$ starting from the initial state $\mathbf{z}$ at time $t=0$ is

[TABLE]

Define $\tau(\epsilon)$ to be the worst-case time needed to be within variation distance $\epsilon$ of the stationary distribution,

[TABLE]

Define the mixing time to be $\tau_{\textrm{mix}}:=\tau(1/4)$, and note that $\tau(\epsilon)=\tau_{\textrm{mix}}\log\epsilon^{-1}$ [34].

To bound the mixing time of $P$ we use the method of canonical paths [29]. A path from $x\in\Omega$ to $y\in\Omega$ is a sequence $\gamma_{x,y}=(v_{1},\ldots,v_{k})$ of states with $v_{1}=x,v_{k}=y$ and $P(v_{i},v_{i}+1)>0$ for $i=1,\ldots,k-1$. A set $\mathcal{P}=\{\gamma_{x,y}\}$ containing a path that connects every pair of states is called a set of canonical paths. If $\ell=\max_{\gamma\in\mathcal{P}}|\gamma|$ is the maximum length of any path in $\mathcal{P}$ and $\pi_{\min}:=\min_{\mathbf{z}\in\Omega}\pi(\mathbf{z})$, then the mixing time $\tau_{\textnormal{mix}}$ is $\mathcal{O}(R\;\ell\ln\pi_{\min}^{-1})$, where the congestion $R$ is defined as

[TABLE]

In order to compute the congestion we will use a standard technique called encoding. For an edge $e=\left(v,v^{\prime}\right)$, let $\mathcal{P}(e)$ be the set of paths in $\mathcal{P}$ which pass through $e$. An encoding is an assignment of an injective function $\eta_{e}:\mathcal{P}(e)\rightarrow\Omega$ to every edge $e\in E(\Omega)$. More generally, it can be useful to consider injective encodings $\eta_{e}:\mathcal{P}(e)\rightarrow\Omega\times\Theta$ that map at most $G=|\Theta|$ paths through the $e$ to each state in $\Omega$ [28]. Let $\eta_{e}(\gamma_{x,y})=(\eta_{e}(x,y),\theta_{e}(x,y))$, and suppose there is some $M$ such that

[TABLE]

then the congestion satisfies

[TABLE]

where $P_{\min}:=\min_{(\mathbf{z},\mathbf{z}^{\prime})\in E(\Omega)}P(\mathbf{z},\mathbf{z}^{\prime})$, (62) follows from property (60), and (63) follows from the injectivity property of the encoding.

4.2 PIMC transition probabilities

To sample from the classical Gibbs distribution (22) the PIMC method uses a Markov chain which chooses a site in the classical lattice $\Lambda$ uniformly at random and proposes to flip the bit at that site to make a transition $\mathbf{z}\rightarrow\mathbf{z^{\prime}}$ with an acceptance probability $P(\mathbf{z},\mathbf{z^{\prime}})$ given by the Metropolis rule [34],

[TABLE]

As long as the state space is a connected graph, the transition matrix $P$ with transitions satisfying (64) will converge to the stationary distribution $\pi$.

For the sake of lower bounding the minimum transition probability $P_{\textrm{min}}$ from the previous section it therefore suffices to lower bound $\pi(\mathbf{z^{\prime}})/\pi(\mathbf{z})$ where $\mathbf{z}$ and $\mathbf{z^{\prime}}$ differ only on a single site. Suppose $\mathbf{z}$ differs from $\mathbf{z^{\prime}}$ at position $(i,j)$, then after canceling common factors in (22) we have

[TABLE]

Depending on the off-diagonal matrix elements in $B_{i}$, the numerator of (65) can potentially be zero. This means that the state space $\Omega$ may not be connected by single bit flips, because the stationary distribution does not have support on all of $\Omega$. This can be avoided in general by adding a 1-local transverse field field with sufficiently small magnitude $\Gamma$ to avoid disturbing the resulting approximation. With this we can lower bound the ratio (65) using the fact that $B_{i}$ is a nonnegative matrix, which implies

[TABLE]

The denominator of (65) is upper bounded by 1, and (66) implies that the numerator is at least $(\beta\Gamma/L)^{2}$, therefore

[TABLE]

4.3 Transverse Ising models

In this section we will describe a set of canonical paths which will be used to show rapid mixing of PIMC applied to Hamiltonians of the form (5). Let $\mathbf{x},\mathbf{y}\in\Omega$ be spin configurations with $\mathbf{x}=(x_{1,1},\ldots,x_{L,1},\ldots,x_{1,n}\ldots,x_{L,n})$ and $\mathbf{y}=(y_{1,1},\ldots,y_{L,1},\ldots,y_{1,n}\ldots,y_{L,n})$.

It will be convenient to choose an ordering on the space $[L]\times[n]$. We order bits first by the second coordinate and then the first, so that $(r,l)\preceq(s,m)$ if $l\leq m$ or if $r\leq s$ and $l=m$. For any $1\leq r\leq L$ and $1\leq m\leq n$ define

[TABLE]

The canonical path $\gamma_{\mathbf{x},\mathbf{y}}$ from $\mathbf{x}$ to $\mathbf{y}$ consists of simply replacing each bit of $\mathbf{x}$ with the corresponding bit of $\mathbf{y}$ following the above ordering. This can be thought of as $n$ steps, each of which consists of $L$ substeps that correspond to transitions of the Markov chain. The $r$-th substep of the $m$-th step consists of flipping the bit in the $r$-th Trotter slice of the $m$-th qubit if $x_{r,m}\neq y_{r,m}$ or leaving it alone if $x_{r,m}=y_{r,m}$. In the former case, this contributes an edge $(\mathbf{v},\mathbf{v^{\prime}})\in E(\Omega)$ with

[TABLE]

In the latter case (i.e. $x_{r,m}=y_{r,m}$), this step does not add an edge to the path. Thus the number of edges in $\gamma_{\mathbf{x},\mathbf{y}}$ equals the number of positions on which $\mathbf{x}$ and $\mathbf{y}$ disagree, which is always $\leq nL$. In Section 4.1 we defined an encoding of a path to be an injective map $\eta_{e}$ from the set of paths through $e$ to the state space $\Omega$. The encoding we’ll use for this path is

[TABLE]

where $\cdot$ to denote elementwise multiplication. To see that $\mathbf{x},\mathbf{y}$ are uniquely determined by $\mathbf{v}$ together with $\boldsymbol{\eta}=\eta_{(\mathbf{v},\mathbf{v^{\prime}})}(\gamma_{\mathbf{x},\mathbf{y}})$, observe that $\mathbf{v}\cdot\boldsymbol{\eta}=\mathbf{x}\cdot\mathbf{y}$, which specifies the bits that are flipped along the path $\gamma_{\mathbf{x},\mathbf{y}}$. Moreover, the edge $(\mathbf{v},\mathbf{v}^{\prime})$ specifies the location $(r,m)$ of the bit being flipped in the transition from $\mathbf{v}$ to $\mathbf{v}^{\prime}$. Let $(\mathbf{v}\cdot\boldsymbol{\eta})_{\leq(r,m)}$ denote the vector which equals $\mathbf{v}\cdot\boldsymbol{\eta}=\mathbf{x}\cdot\mathbf{y}$ on the coordinates in $[\leq(r,m)]$ and equals 1 on the coordinates in $[>(r,m)]$; define $(\mathbf{v}\cdot\boldsymbol{\eta})_{>(r,m)}$ analogously. Then we can recover $\mathbf{x}$ by calculating $\mathbf{v}^{\prime}\cdot(\mathbf{v}\cdot\boldsymbol{\eta})_{\leq(r,m)}$, and can recover $\mathbf{y}$ by calculating $\mathbf{v}^{\prime}\cdot(\mathbf{v}\cdot\boldsymbol{\eta})_{>(r,m)}$.

For any subset of the lattice sites $A\subseteq\Lambda$, define the restriction $E_{\beta}^{A}$ of the classical energy function $E_{\beta}$ by restricting the sum in (9) to sites $(i,j)\in A$, so that for all $\mathbf{z}\in\Omega$ we have

[TABLE]

Now we compute

[TABLE]

where

[TABLE]

We also need

[TABLE]

where

[TABLE]

and $B_{Y}$ is defined similarly. The value of $M$ needed to satisfy (60) for this encoding can now be determined by

[TABLE]

To upper bound the exponent in the worst-case define $J:=\max_{m=1,\ldots,n}J_{m}$, then

[TABLE]

To upper bound the sum we use the assumption $K^{zz}_{jk}\leq|j-k|^{-(2+\xi)}$ and regard it as a lower Riemann sum for an integral,

[TABLE]

and since $\xi>0$ it suffices in (63) to take $M=\mathcal{O}\left(e^{8J+4\beta}\right)$. If the 1D quantum system has periodic boundary conditions connecting qubit 1 to qubit $n$, then the same analysis would show that $M=\mathcal{O}\left(e^{8(J+\beta)}\right)$ suffices for (60). Therefore by (61)-(63) the congestion $R$ satisfies

[TABLE]

the factor of $2nL$ comes from the chain being lazy and the choice of random single-site updates. Note that $\max_{(\mathbf{z},\mathbf{z^{\prime}})}\left(\pi(\mathbf{z^{\prime}}/\pi(\mathbf{z})\right)$ is ${\cal O}(L^{2}\Gamma^{-2})$ from (67) and also $\ell:=\max_{\mathbf{x},\mathbf{y}}|\gamma_{\mathbf{x},\mathbf{y}}|=nL$ and $\ln\pi_{\min}^{-1}=\ln(Z\cdot e^{nLJ+n\beta})\leq nL(J+\beta+\ln 2)$, using $Z\leq 2^{nL}$. In Section 3.1 we show that $L=\operatorname{poly}(n,\delta_{\textrm{mult}}^{-1})$ suffices to approximate the partition function with the desired precision. Therefore the mixing time bounds in Section 4.1 imply that the PIMC method generates samples as was assumed in Section 3.3, which completes the proof of theorem 1 for models of the form (5).

4.4 Polynomial time mixing in theorem 1

In this section we will prove rapid mixing of PIMC for Hamiltonians of the form (6) by describing paths between arbitrary states $\mathbf{x},\mathbf{y}\in\Omega$. It will be convenient to represent states in the form $\mathbf{z}=[\bar{z}_{1},\ldots,\bar{z}_{n}]$ with $\bar{z}_{j}:=(z_{1,j},\ldots,z_{L,j})$ denoting the spin values along the $j$-th worldline in configuration $\mathbf{z}\in\Omega$. We say that a “double jump” occurs in configuration $\mathbf{z}$ at position $(i,j)$ if $\left(1+z_{i,j}z_{i+1,j}\right)\left(1+z_{i,j+1}z_{i+1,j+1}\right)\neq 0$. For any consecutive pair of worldlines $j,j+1$, define the number of double jumps $d(\bar{z}_{j},\bar{z}_{j+1})$,

[TABLE]

For any configuration $\mathbf{z}$ and worldline $j$ we will define a path segment to a new configuration $[\bar{z}_{1},\ldots,\bar{z}_{j-2},\bar{z}^{\prime}_{j-1},\bar{z}^{\prime\prime}_{j},\bar{z}^{\prime}_{j+1},\bar{z}_{j+2},\ldots,\bar{z}_{n}]$ with all of the double jumps removed from worldline $j$,

[TABLE]

The double prime notation indicates that the $j$-th worldline does not contain any double jumps, while the single prime indicates a worldline that shares no jumps with its double primed neighboring worldline.

There are many possible paths to choose from for removing double jumps in the $i$-th worldline. Due to the form of the Hamitonian (6), changing the placement of a double jump along the imaginary-time direction has no effect on the energy of the effective classical system, except possibly for the part which arises from the diagonal part of the Hamiltonian. Similar reasoning holds for the ordering of the single and double jumps; changing the ordering can only change the effective classical energy through the diagonal part of the Hamiltonian. Similarly, changing the imaginary-time position of the double jumps in $\bar{x}_{i-1},\bar{x}_{i}$ will have no effect on the contributions to the effective classical energy that comes from double jumps occuring in $\bar{x}_{i-2},\bar{x}_{i-1}$. Finally, note that changing the position of a double jump along the imaginary-time direction can only affect the classical energy through the diagonal part of the Hamiltonian, reducing the probability of the configuration by at most $e^{-4\beta}$ (the factor of 4 is due to the potential change from the diagonal part of the Hamiltonian when 3 worldlines are changed arbitrarily).

From the above considerations it follows that the imaginary-time position of the double jumps in a particular worldline can be changed freely with the stationary weight of the resulting configuration decreased by at most a factor of $e^{-4\beta}$ (the factor of 4 bounds the potential change from the diagonal part of the Hamiltonian when 3 worldlines are changed arbitrarily). Furthermore, since double jumps arising from the term $\sigma^{x}_{i}\sigma^{x}_{i+1}$ are simply bit-flips, two such jumps will cancel out (“annihilate”) if they coincide. This allows for paths that remove these double jumps in pairs simply by moving them around, with the additional guarantee that this will not lower the stationary weight by more than $e^{-4\beta}$. Therefore the path from $[\bar{z}_{1},\ldots,\bar{z}_{n}]$ to $[\bar{z}_{1},\ldots,\bar{z}_{j-2},\bar{z}^{\prime}_{j-1},\bar{z}^{\prime\prime}_{j},\bar{z}^{\prime}_{j+1},\bar{z}_{j+2},\ldots,\bar{z}_{n}]$ is defined to remove double jumps in pairs by shifting the position of the double jump at the least imaginary-time position along the imaginary-time direction until it encounters the next double jump and the pair is annihilated.

The canonical path from $\mathbf{x}$ to $\mathbf{y}$ proceeds through the worldlines in order $1,\ldots,n$, prepping each pair of consecutive worldlines to remove double jumps, then updating the $j$-th worldline from $\bar{x}^{\prime\prime}_{j}$ to $\bar{y}^{\prime\prime}_{j}$. The full update of the first worldline proceeds as,

[TABLE]

The purpose of this series of steps (each of which consists of many single-site updates) is to avoid passing through intermediate configurations with too many additional double jumps occuring across the 1D domain wall formed near the worldlines which are being updated. The worst case intermediate configurations along the paths will have the form,

[TABLE]

This system of paths is highly degenerate (a particular configuration such as (81) will appear in many paths). To bound the congestion we will use an encoding function $\eta_{(\mathbf{z},\mathbf{z}^{\prime})}:\mathcal{P}(\mathbf{z},\mathbf{z}^{\prime})\rightarrow\Omega\times\Theta$. If $\eta_{(\mathbf{z},\mathbf{z}^{\prime})}(\gamma_{\mathbf{x}\mathbf{y}})=(\mathbf{\eta},\theta)$, then

[TABLE]

and $\theta$ contains all of the information needed to reconstruct $\bar{x}_{j-1},\bar{x}_{j},\bar{x}_{j+1},\bar{y}_{j-1},\bar{y}_{j},\bar{y}_{j+1}$ from $\bar{x}^{\prime}_{j-1},\bar{x}^{\prime\prime}_{j}$, $\bar{x}^{\prime}_{j+1},\bar{y}^{\prime}_{j-1},\bar{y}^{\prime\prime}_{j},\bar{y}^{\prime}_{j+1}$. This lost information must specify the location of all of the jumps that have been removed from these 6 worldlines. This means that the encoding sends exponentially paths to each state, and so to obtain a polynomial time mixing bound we have to slightly generalize (63) to account for the fact that the stationary probability on the endpoints of most of these paths is sufficiently small.

For any $\mathbf{v}\in\Omega$ define $d_{j}(\mathbf{v}):=d(\bar{v}_{j-1},\bar{v}_{j})+d(\bar{v}_{j},\bar{v}_{j+1})$ (which is just the total number of double jumps in the $j$-th worldline of $\mathbf{v}$). Let $q$ be the number of additional double jumps in the $j$-th worldlines of $\mathbf{x}$ and $\mathbf{y}$ which are not present in the $j$-th worldlines of $\mathbf{z}$ and $\boldsymbol{\eta}$,

[TABLE]

Define $\mathcal{P}^{q}_{\mathbf{z},\mathbf{z}^{\prime}}$ to be the subset of canonical paths passing through $(\mathbf{z},\mathbf{z}^{\prime})$ with endpoints $\mathbf{x}$ and $\mathbf{y}$ that satisfy (83). With these definitions the congestion (59) becomes

[TABLE]

Let $\Theta=\cup_{q=1}^{6L}\Theta_{q}$ where $G(q)=|\Theta_{q}|={{6L}\choose{q}}$ is the number of ways to restore $q$ jumps to the undetermined worldlines. Now suppose we find $M(q)$ such that

[TABLE]

for all $\gamma_{\mathbf{x}\mathbf{y}}\in\mathcal{P}^{q}_{\mathbf{z},\mathbf{z}^{\prime}}$, then this implies

[TABLE]

To satisfy (85) we need $M(q)$ to satisfy

[TABLE]

where (90) to (91) follows because the the terms in the classical energy function which correspond to off-diagonal parts of the Hamiltonian involving qubits $j-2$ and $j+2$ are unaffected by the potentially changed spin values in the primed $j-1$ and $j+1$ worldlines. The factor of $e^{4\beta}$ comes from upper bounding the change in the classical energy function due to the diagonal terms of the Hamiltonian between qubits $j-2,j-1$ and $j+1,j+2$. The effect of the diagonal part of the Hamiltonian on the ratio expression remaining in (91) can be bounded by a factor of $e^{8\beta}$. The additional single jumps introduced near the domain wall can increase the ratio by at most $L^{12}$ (using (67), and noting that the worst case occurs when 2 new single jumps are introduced in each of 3 worldlines, in each of $\mathbf{z}$ and $\mathbf{\eta}$). Each of the double jumps in the numerator which are not present in the denominator decreases the ratio by at least $(2\beta/L)$. Therefore we may take $M(q)=e^{12\beta}(2\beta)^{q}L^{-q}$ and (88) becomes

[TABLE]

Using the fact that $P_{\min}^{-1}$ is $\operatorname{poly}(n,\beta,\Gamma^{-1})$ from (67), and the fact that the maximum length of the canonical paths used in this section is $\mathcal{O}(n^{2}L^{2})$, this shows that $\tau_{\textnormal{mix}}$ is $\operatorname{poly}(n,e^{\beta},\Gamma^{-1})$ as required to complete the proof of theorem 1.

4.5 Quasipolynomial time mixing in theorem 1

In this section we also work with a restricted state space $\Omega^{*}$, which is defined as the set of configurations with fewer than $2B:=4\beta\log(n)$ jumps in each of the worldlines,

[TABLE]

In Section 3.4 we show that $Z^{*}=\sum_{\mathbf{z}\in\Omega^{*}}e^{-\beta E_{\beta}(\mathbf{z})}$ satisfies $|Z-Z^{*}|\leq\delta Z$. The PIMC dynamics can be restricted to $\Omega^{*}$ by rejecting moves which would leave $\Omega^{*}$, or alternatively the canonical paths in this section bound the mixing time within $\Omega^{*}$ by the leaky random walk conductance bound in [15]. To show rapid mixing of PIMC within $\Omega^{*}$ we assign a canonical path $\gamma_{\mathbf{xy}}$ to every $\mathbf{x},\mathbf{y}\in\Omega^{*}$. The simple case for these paths occurs when $\mathbf{x},\mathbf{y}$ are in a subset $\Omega^{*}_{\textrm{inner}}\subseteq\Omega^{*}$ which only allows half as many jumps per worldline as $\Omega^{*}$ does,

[TABLE]

If $\mathbf{x},\mathbf{y}\in\Omega^{*}_{\textrm{inner}}$ then we use the same paths which were used to show rapid mixing for TIM in Section 4.3. Various intermediate points $\mathbf{v}$ along this path $\gamma_{\mathbf{x},\mathbf{y}}$ may leave $\Omega^{*}_{\textrm{inner}}$ but will remain in $\Omega^{*}$ since

[TABLE]

may have as many as $B$ jumps in the first $r$ time slices of $\bar{y}_{m}$ and $B$ jumps in the $L-r$ latter time slices of $\bar{x}_{m}$, so $\bar{v}_{m}$ will have at most $2B$ jumps in total and so $\mathbf{v}\in\Omega^{*}$.

For points $\mathbf{x},\mathbf{y}\in\Omega^{*}$ which are not in $\Omega^{*}_{\textrm{inner}}$, we must not produce intermediate points with too many jumps in any worldline. To accomplish this the paths will have a “clean up” step similar to the one which was used in the previous section, to reduce the number of 2-local jumps in a worldline. Consider an intermediate point of the path which has so far updated the worldlines of the first $m$ qubits (which could be $m=0$), $\mathbf{v}=[\bar{y}_{1},\ldots,\bar{y}_{m},\bar{x}_{m+1},\ldots,\bar{x}_{n}]$. If $\bar{x}_{m+1}$ has more than $B$ double jumps then we remove these in imaginary-time order. Since we can no longer “annihilate” these transitions in pairs as was done in the previous section, we instead dismantle them into 1-local jumps (relying on the fictitious transverse field) and then annihilate those against one another. The effect of this process on the stationary weight of the configuration will be analyzed after we have given the encoding function for these paths. The net result of this process results is a worldline $\bar{x}^{\prime}_{m+1}$ with fewer than $B$ jumps. Similarly let $\bar{y}^{\prime}_{m+1}$ be defined from $\bar{y}_{m+1}$ with all but the first $B$ jumps removed. Next we update $\bar{x}^{\prime}_{m+1}$ to $\bar{y}^{\prime}_{m+1}$ in the usual way, and since $\bar{x}^{\prime}_{m+1}$ and $\bar{y}^{\prime}_{m+1}$ each have fewer than $B$ jumps, the intermediate configurations along this path will have less than $2B$ jumps in every worldline and therefore remain in $\Omega^{*}$.

The encoding $\eta:\mathcal{P}_{(\mathbf{z},\mathbf{z}^{\prime})}\rightarrow\Omega^{*}\times\Theta$ for these paths is again defined by (82). This time the maximum degeneracy of the encoding is $\binom{4L}{B}$, since there are at most $4L$ positions to which at most $B$ 2-local jumps can be restored in the undetermined worldlines. Unfortunately, this time we cannot have a inequality analogous to (85), because the additional double jumps present in $\mathbf{x}$ and $\mathbf{y}$ can completely change the bit-values in wordlines $j-1$ and $j+1$. Unlike the case in the previous section these altered spin values can reduce the classical energy drastically because of the $h$ terms acting on worldlines $j-2,j-1$ and $j+1,j+2$. The stationary weight of $\mathbf{z}$ and $\boldsymbol{\eta}$ can be lowered from $\mathbf{x}$ and $\mathbf{y}$ by a factor which scales exponentially in the number of double jumps that $\mathbf{x}$ and $\mathbf{y}$ have in worldlines $j-2,j-1$ and $j+1,j+2$. Therefore,

[TABLE]

The fictitious transverse field has magnitude $\delta_{\textrm{mult}}/n$, so the congestion satisfies

[TABLE]

from which it follows that the mixing time $\tau_{\textrm{mix}}$ is $\mathcal{O}(2^{\beta(\log n\delta_{\text{mult}}^{-1})^{2}})$ and we obtain the overall PIMC runtime as claimed in theorem 1.

Acknowledgments

Part of this work was completed while EC was funded by the Institute for Quantum Information and Matter, an NSF Physics Frontiers Center (NSF Grant PHY-1125565) with support of the Gordon and Betty Moore Foundation (GBMF-12500028). AWH was funded by NSF grants CCF-1452616, CCF-1629809, CCF-1729369 and ARO contract W911NF-17-1-0433. The research is based upon work partially supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.