# Efficient thermalization and universal quantum computing
with quantum Gibbs samplers

Cambyse Rouzé<EMAIL_ADDRESS>Inria, Télécom Paris - LTCI, Institut
Polytechnique de Paris, 91120 Palaiseau, France  Daniel Stilck França
<EMAIL_ADDRESS>Univ Lyon, ENS Lyon, UCBL, CNRS, Inria, LIP,
F-69342, Lyon Cedex 07, France  Álvaro M. Alhambra<EMAIL_ADDRESS>Instituto de Física Téorica UAM/CSIC, C. Nicolás Cabrera 13-15, Cantoblanco,
28049 Madrid, Spain

###### Abstract

The preparation of thermal states of matter is a crucial task in quantum
simulation. In this work, we prove that a recently introduced, efficiently
implementable dissipative evolution thermalizes to the Gibbs state in time
scaling polynomially with system size at high enough temperatures, and for any
Hamiltonian that satisfies a Lieb-Robinson bound, such as local Hamiltonians
on a lattice. Furthermore, we show the efficient adiabatic preparation of the
associated purifications or “thermofield double” states. To the best of our
knowledge, these are the first results rigorously establishing the efficient
preparation of high temperature Gibbs states and their purifications. In the
low-temperature regime, we show that implementing this family of dissipative
evolutions for inverse temperatures logarithmic in the system’s size is
polynomially equivalent to standard quantum computation. On a technical level,
for high temperatures, our proof makes use of the mapping of the generator of
the evolution into a Hamiltonian, and then analysing it as perturbation of the
Hamiltonian corresponding to infinite temperature. For low temperature, we
instead perform a perturbation at zero temperature of the Laplace transform of
the energy observable at fixed runtime, and resort to circuit-to-Hamiltonian
mappings akin to the proof of universality of quantum adiabatic computing.
Taken together, our results show that a family of quasi-local dissipative
evolutions efficiently prepares a large class of quantum many-body states of
interest, and has the potential to mirror the success of classical Monte Carlo
methods for quantum many-body systems.

##  I Introduction

Markov Chain Monte Carlo (MCMC) methods are one of the go-to tools for
sampling from classical Gibbs states (GS) of spin systems [1]. They exhibit
efficiency in practice [2] and are also provably efficient in some cases, such
as at high temperatures [3]. These methods traditionally rely on “local
updates”, where a few spins are flipped based on local energy changes, a
feature crucial for both practical implementation and theoretical analysis.
However, extending such algorithms to quantum systems has been a formidable
challenge. Despite significant research efforts [4, 5, 6, 7, 8, 9, 10, 11,
12], general quantum algorithms capable of (quasi-)local updates that also
provably converge to a quantum GS have remained elusive, until the recent
breakthroughs of [13, 14].

The present work delves into the quantum algorithm proposed in [13, 14],
showing its ability to efficiently prepare quantum states of interest. We
investigate both the high-temperature and the low-temperature regimes. In the
high-temperature regime, we show that above a certain constant threshold
temperature, the Lindbladians constructed in [13] efficiently converge to the
GS for all Hamiltonians obeying a Lieb-Robinson bound, which includes local
Hamiltonians on a lattice [15]. To the best of our knowledge, this is the
first unconditional result that shows efficient preparation of high-
temperature quantum Gibbs states on a quantum computer for any lattice
dimension [5, 6, 10, 8]. Interestingly, our results also imply the efficient
adiabatic preparation (starting with the state at inverse temperature
$\beta=0$) of the associated purifications of thermofield doubles. The
preparation of such purified states is of independent interest [16, 17, 18],
and this scheme is the first provably efficient for a general class of models.

For low temperatures, we show that implementing the Lindbladians of [13] for
$\beta=\Omega(\log(n))$ already gives rise to a universal model for quantum
computing equivalent to the circuit model. To this end, we show that the
Lindbladians corresponding to circuit-to-Hamiltonian mappings efficiently
converge to a state with a polynomial overlap with the ground state
corresponding to the output of the circuit. This state can then be leveraged
to efficiently obtain the ground state through a series of local measurements.
These results are in analogy with those of [19], where the authors established
similar results for temperatures $\beta=\textrm{poly}(n)$ and polynomial
rounds of mid-circuit measurements, compared to one round at the end for our
construction. However, we may need to repeat the procedure a polynomial number
of times to succeed in preparing the state. Thus, our construction has an
exponentially better dependency on the temperature and does not require mid-
circuit measurements.

On a technical level, for high temperature, our analysis establishes a
spectral gap for the Lindbladians defined in [13]. This is done by viewing the
high-temperature generator as a perturbation of its infinite temperature
counterpart, and employing results on quasi-local perturbations of gapped
Hamiltonians [20]. At low temperature, we also follow a perturbative approach,
but now at $\beta=+\infty$. We are inspired by the classical results of [21]
and the recent quantum results of [19] to show that the Lindbladians of [13]
have the power to simulate geometrically local quantum circuits already at
moderate temperatures $\beta=\Omega(\log(n))$. Taken together, our findings
show that the algorithm of [13] is a powerful tool to prepare non-trivial
quantum many-body states, and to replicate the success of classical MCMC in
the quantum realm.

##  II Quantum Gibbs sampling

We start by recalling the generator of the quantum Monte Carlo algorithm
recently introduced in [13]. For $\beta>0$ and a Hamiltonian $H$ over the
lattice $\Lambda\equiv[0,L]^{D}$ with $n:=|\Lambda|$, we denote
$\gamma(\omega):=\operatorname{exp}(-{(\beta\omega+1)^{2}}/{2})$ and we
consider the Lindbladian

$\displaystyle\mathcal{L}_{\beta}(\rho)=-i[B,\rho]+\sum_{a\in\Lambda,\alpha\in[3]}\int_{-\infty}^{\infty}\gamma(\omega)\mathcal{D}^{a,\alpha}_{\omega}(\rho)\,d\omega$
(1)

where $\mathcal{D}^{a,\alpha}_{\omega}$ stands for the dissipative term
corresponding to the jump operators

$\displaystyle
A^{a,\alpha}(\omega):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{iHt}A^{a,\alpha}e^{-iHt}e^{-i\omega
t}\,f(t)\,dt\quad$

with
$f(t):=\operatorname{exp}(-{t^{2}}/{\beta^{2}})\,\sqrt{\beta^{-1}\sqrt{2/\pi}}\,$.
Above, $A^{a,1}\equiv\sigma_{x},\,A^{a,2}\equiv\sigma_{y}$ and
$A^{a,3}\equiv\sigma_{z}$ denote the $1$-local Pauli matrices on site $a$. The
coherent part is defined through the self-adjoint matrix

$\displaystyle
B=\sum_{a\in\Lambda,\alpha\in[3]}\,\int_{-\infty}^{\infty}\beta^{-2}b_{1}(t/\beta)e^{-iHt}$
$\displaystyle\qquad\times\int_{-\infty}^{\infty}b_{2}(t^{\prime}/\beta)e^{iHt^{\prime}}A^{a,\alpha}e^{-2iHt^{\prime}}A^{a,\alpha}e^{iHt^{\prime}}dt^{\prime}e^{iHt}dt$

for some rapidly decaying functions $b_{1},b_{2}$ with
$\|b_{1}\|_{L_{1}},\|b_{2}\|_{L_{1}}\leq 1$ whose exact expression we recall
in Appendix A.2. The Lindbladian constructed above has the main advantage of
being both quasi-local and reversible in the Heisenberg picture:

$\displaystyle\langle
X,\mathcal{L}^{\dagger}_{\beta}(Y)\rangle_{\sigma_{\beta}}=\langle\mathcal{L}^{\dagger}_{\beta}(X),Y\rangle_{\sigma_{\beta}}$

where $\langle
X,Y\rangle_{\sigma_{\beta}}:=\operatorname{tr}(X^{\dagger}\sigma_{\beta}^{\frac{1}{2}}Y\sigma_{\beta}^{\frac{1}{2}})$
denotes the KMS inner product, and where the Gibbs state
$\sigma_{\beta}\propto{\operatorname{exp}(-\beta H)}$ is the unique fixed
point of the evolution generated by $\mathcal{L}_{\beta}$.

##  III Spectral gap for high-temperature quantum Gibbs sampler

Next, we consider a local Hamiltonian $H$ on the lattice $\Lambda$, defined
through a function that maps any non-empty finite set $X\subset\Lambda$ to a
self-adjoint element $h_{X}$ supported in $X$ with
$\max_{X}\left\|h_{X}\right\|_{\infty}\leq h$. The Hamiltonian of the system
in a region $A\subset\Lambda$ is then defined as the sum over interactions
included in $A$:

$\displaystyle H_{A}:=\sum_{X\subseteq A}h_{X}\,.$

In addition, we assume that the Hamiltonian is such that no interaction
$h_{X}$ has support on more than $k$ sites, and each site $u$ appears on at
most $l$ non-zero operators $h_{X}$. We refer to such Hamiltonian as
$(k,l)$-local. As often, we will make the implicit assumption that $H$ is one
element of a family of Hamiltonians indexed by the lattice size, in which case
the constants $h,k$ and $l$ are assumed to be independent of $L$. We show that
for any $\beta\leq\beta^{*}=\mathcal{O}(1)$ the spectral gap of
$-\mathcal{L}_{\beta}^{\dagger}$, i.e. the smallest nonzero eigenvalue, is
constant.

######  Theorem 1.

In the notations of the previous paragraph, for any $(k,l)$-local Hamiltonian
$H$, there exists a constant $\beta^{*}>0$ independent of $L$ such that for
$\beta\leq\beta^{*}$, the spectral gap of $-{\mathcal{L}}^{\dagger}_{\beta}$,
is lower bounded by $\frac{1}{2\sqrt{2}e^{1/4}}$.

By a standard argument recallled in Appendix B, Theorem 1 implies the
following guarantee on the convergence of the Gibbs sampler:

######  Corollary 2.

For any $\epsilon>0$, the evolution generated by $\mathcal{L}_{\beta}$ gets
$\epsilon$-close to $\sigma_{\beta}$ in polynomial time, i.e. for any initial
state $\rho$,

$\displaystyle\|e^{t\mathcal{L}_{\beta}}(\rho)-\sigma_{\beta}\|_{1}\leq\epsilon\quad\text{
for all}\quad t=\Omega(\log(1/\epsilon)+|\Lambda|)\,.$ (2)

Furthermore, we can prepare an $\epsilon$-approximation of $\sigma_{\beta}$ on
a quantum computer with $\tilde{\mathcal{O}}(|\Lambda|)$ Hamiltonian
simulation time and $\tilde{\mathcal{O}}(|\Lambda|)$ two qubit gates, where
$\tilde{\mathcal{O}}$ omits terms poly-logarithmic in $\epsilon^{-1}$ and
$|\Lambda|$.

Given Eq. 2, the result on the efficient preparation on a quantum computer
immediately follows from [13, Theorem 1.2], where the authors show how to
efficiently simulate the Lindbladians of Eq. (1) on a quantum computer. In the
next paragraphs, we briefly explain the proof steps to obtain the bound on the
spectral gap, which are further detailed in Appendix B.

##### Mapping the Lindbladians to Hamiltonians:

for any $\beta>0$ we define as in [14]

$\displaystyle\widetilde{\mathcal{L}}_{\beta}(X):=\sigma_{\beta}^{-\frac{1}{4}}\mathcal{L}_{\beta}(\sigma_{\beta}^{\frac{1}{4}}X\sigma_{\beta}^{\frac{1}{4}})\sigma_{\beta}^{-\frac{1}{4}}\equiv\sum_{a,\alpha}\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}\,.$
(3)

Letting $\ket{\sqrt{\sigma_{\beta}}}$ be the vectorization of
$\sqrt{\sigma_{\beta}}$, we see that it is in the kernel of
$\widetilde{\mathcal{L}}_{\beta}$. Furthermore, as
$\widetilde{\mathcal{L}}_{\beta}$ is related to the original Lindbladian by a
similarity transformation, they both share the same spectrum, and
$\widetilde{\mathcal{L}}_{\beta}$ is self-adjoint with respect to the Hilbert-
Schmidt inner product for all $\beta>0$, by KMS symmetry of
$\mathcal{L}_{\beta}^{\dagger}$. Thus, we can interpret
$\widetilde{\mathcal{L}}_{\beta}$ as a Hamiltonian. Furthermore, its Hilbert
space is independent of $\sigma_{\beta}$.

##### Decomposing the generator:

we can write the generator as the following telescopic sum:

$\displaystyle\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}=\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,0}+\sum_{r=0}^{\infty}\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r+1}-\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r},$

where for $r>0$ and a site $a$,
$\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r+1}$ is defined as the generator
for the state with the reduced Hamiltonian $H_{B_{a}(r)}$ associated to the
ball $B_{a}(r)$ centered in site $a$ and of radius $r$. Thus,
$\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r+1}-\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r}$
is strictly supported on $B_{a}(r+1)$.

##### The $\beta\to 0$ limit:

by construction, for $\beta\to 0$, the $\widetilde{\mathcal{L}}_{\beta}$
converges to the generator of the $1-$local depolarizing channel. When
performing the mapping to Hamiltonians described in Eq. (3), we thus obtain a
Hamiltonian $\tilde{H}_{0}$ that is frustration-free, gapped, and satisfies
local topological quantum order (LTQO) as in [22, 20]. It is straightforward
to show these properties, so we leave the detailed derivations to [23].
Showing the convergence to the depolarizing semigroup requires slightly more
work. We sketch the argument here, leaving the details for [23]. In essence,
the choice of the functions $f,b_{1}$ and $b_{2}$ ensures that they become
delta functions in the limit $\beta\to 0$. It is not difficult to see that at
$\beta=0$, the matrices we integrate over are just proportional to the jump
operators $A^{a,i}$, which are those of the depolarizing channel.

##### Quasi-locality of the generator:

the next crucial step is to bound the generators in
$\widetilde{\mathcal{E}}_{\beta,r}^{a,\alpha}:=\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r+1}-\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r}$
as a function of $r$ and $\beta$. By using Lieb-Robinson bounds, we find that,
for small enough $\beta$

$\left\|\mathcal{E}_{\beta,r}^{a,\alpha}\right\|_{2\rightarrow
2}\leq\mathcal{O}\left(\frac{\beta^{r}}{\beta^{r}_{0}}\right)\,,$ (4)

for some constant $\beta_{0}$ depending on the locality constants $k,l$ of the
Hamiltonian $H$. This means that we can see the finite $\beta$ contributions
to $\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}$ as a quasi-local perturbation
whose strength vanishes with $\beta<\beta_{0}$ (see [23] for further details).
With this, we find that $\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}$ satisfies
all the assumptions from [20], from which we can guarantee the lower bound on
the gap in Theorem 1. The expression of $\beta^{*}$ then depends on
$\beta_{0}$, the Lieb-Robinson velocity of $H$, and the lattice dimension. The
proof extends straightforwardly to systems with exponentially decaying
interactions with similar Lieb-Robinson bounds [24], and it is likely to also
extend to systems with polynomially decaying interactions [25, 26]

### Comparison with previous work

There is an extensive literature on sampling from quantum GS on a quantum
computer. It includes numerous schemes with exponential runtimes for general
Hamiltonians and temperatures [27, 4, 28, 12, 29, 11], including a range of
variational algorithms [30, 31]. Others are aimed at mimicking the natural
thermalization process of open systems [32, 33, 34], for a fast covergence is
only guaranteed in restricted cases such as commuting Hamiltonians [6, 10]. To
the best of our knowledge, Corollary 2 together with the algorithm of [13] are
the first to give efficient (polynomial-time) preparation for general non-
commuting models, even with the restriction to high temperatures. Note that
Corollary 2 also works with exponentially decaying and noncommuting
interactions, and does not require any assumptions on the Gibbs state besides
the temperature. Previous polynomial-time algorithms are limited to 1D [35] or
commuting Hamiltonians [36]. A result with a comparable degree of generality
is that of [5]. However, there the authors need to assume conditions on the
decay of correlations of the Gibbs state. Although these are believed to hold
at high enough temperatures, they also require the application of gates acting
on $\log^{D}(n)$ qubits, where $D$ is the dimension of the lattice. Except for
$D=1$, standard compilation results would then give a quasi-polynomial
$e^{\mathcal{O}(\log^{D}(n))}$ overhead to compile such gates, so the
resulting algorithm is quasi-polynomial.

##  IV Adiabatic preparation of the purified Gibbs state

The proof of Theorem 1 involves a lower bound on the gap of the Lindbladian,
done by mapping this Lindbladian to a Hamiltonian with the same gap, but in a
doubled Hilbert space. By construction, the Hamiltonian
$\widetilde{\mathcal{L}}_{\beta}$ is quasi-local, and has the purified Gibbs
state as a fixed point. This means that $\widetilde{\mathcal{L}}_{\beta}$ is a
provable realization of Hamiltonians such as those suggested in e.g. [37, 16].

Additionally, the fact that we can control the gap starting from $\beta=0$
means that we can straightforwardly devise an adiabatic scheme to prepare the
purified Gibbs state. The adiabatic path consists of lowering the temperature
from $\beta=0$, so that our parametrized Hamiltonian is
$\mathscr{H}(s)\equiv\widetilde{\mathcal{L}}_{s\beta}$ with $s\in[0,1]$, and
$\beta\leq\beta^{*}$. This way, by initially preparing the GS of
$\mathscr{H}(0)$, which consists of a product of Bell pairs
$\bigotimes_{i}\ket{\Phi^{+}}_{i}$, we can adiabatically sweep towards the
purified Gibbs state (or _thermofield double_) $\ket{\sqrt{\sigma_{\beta}}}$.

######  Theorem 3.

In the notations of the previous paragraphs, and given $\beta^{*}$ as
introduced in Theorem 1, for any $\beta\in[0,\beta^{*}]$, the minimum time
required to adiabatically prepare a state $\epsilon$-close to the purified
Gibbs state at inverse temperature $\beta$ is

$\displaystyle T_{\operatorname{ad}}=\mathcal{O}\left(\frac{(\beta
n)^{3}}{\epsilon^{2}}\right)\,.$ (5)

The proof follows from the lower bound of the gap from Theorem 1, as well as
standard estimates on the runtime of the adiabatic algorithm [38]. To obtain
the estimate of Eq. 5, in [23] we give upper bounds to the first and second
derivatives of the Hamiltonian $\mathscr{H}(s)$ with respect to the parameter
$s$.

The generator $\widetilde{\mathcal{L}}_{\beta}$ thus provides a polynomial-
time scheme to prepare high temperature purified Gibbs states. Being able to
prepare these states is key in quantum simulations of certain models of
entangled black holes [39, 16, 18, 40], and their generation allows for the
measurement of OTOCs [17]. Because of this, many algorithms to prepare them
have been devised [41, 42, 30, 43], and even implemented in small experiments
[44]. However, these are largely limited to variational approaches, without
guarantees of efficiency or scalability. In contrast, the performance of our
adiabatic algorithm at high temperatures is at most cubic in system size and
quadratic in the inverse precision.

##  V Zero-temperature Gibbs sampling and universal quantum computing

Next, we consider the near-zero temperature regime, where we prove the
existence of Gibbs samplers that reach states with inverse polynomial overlap
with the ground states of $\mathsf{BQP}$-hard Hamiltonians at mildly large
$\beta=\Omega(\log(n))$ in polynomial time. This result, which we prove in
Appendix C, leads to a dissipative variant of the well-known equivalence
between the complexity class $\mathsf{BQP}$ of decision problems which can be
solved through the use of polynomially-sized uniform families of quantum
circuits, and the class $\mathsf{AdiabQP}$ of decision problems decidable by
means of unitary quasi-adiabatic evolutions [45, 38, 46, 47]. More precisely,

######  Definition 4 (The class $\mathsf{GibbsQP}$).

This is the class of decision problems decidable through the measurement of
arbitrary $5$-qubit observables at the output of polynomially sized uniform
families of Gibbs samplers associated with $(k,l)$-local Hamiltonians for some
constants $k,l$ independent of the instance size, inverse temperature
$\beta=\Omega(\log(n))$, initialized in the maximally mixed state, and run up
to time $T=\mathcal{O}(\operatorname{poly}(n))$, where $n$ denotes the
problem’s instance size. As commonly assumed, the precision of the
implementation of real parameters involved is inverse polynomially small in
the size of the problem instance.

######  Theorem 5.

$\mathsf{BQP}=\mathsf{AdiabQP}=\mathsf{GibbsQP}$.

Before providing some details about the proof of Theorem 5, which we postpone
to Appendix C.3, let us make some observations and comparison to the existing
literature on dissipative quantum computing. As previously mentioned, the
equivalence $\mathsf{BQP}=\mathsf{AdiabQP}$ was derived in the seminal paper
[45], and later refined in [47] by exhibiting a simpler circuit-to-Hamiltonian
(CTH) construction with a family of $2$-local, $2$-dimensional, qubit
Hamiltonians (see also previous refinements in [38, 46, 47]). That
$\mathsf{GibbsQP}\subseteq\mathsf{BQP}$ follows from the quantum algorithm
prescribed in [14, 13]. Here, we are proving the opposite inclusion that
$\mathsf{BQP}\subseteq\mathsf{GibbsQP}$ by making use of circuit-to-
Hamiltonian constructions akin to those appearing in the proofs of
$\mathsf{BQP}\subseteq\mathsf{AdiabQP}$.

##### Proof sketch

: Similar to the quantum adiabatic algorithm, our method consists of encoding
the output states of a uniform, polynomially sized family of quantum circuits
$C$ into the ground states of a uniform family of local Hamiltonians $H_{C}$.
However, while the adiabatic theorem requires a path between $H_{C}$ and a
trivial Hamiltonian with gap closing at most polynomially fast, we instead
need to show the rapid convergence towards states close to the ground state of
$H_{C}$ for a Gibbs sampling algorithm with a properly selected filter
function $\gamma(\omega)$. It turns out that the Gaussian filter introduced
before Eq. (1) does not work, as the Markov chain in the eigenbasis does not
necessarily converge to the ground state. However, it was conjectured in [13]
that better mixing properties might be attained by means of linear
combinations of Gaussian filters without loss of KMS symmetry. In particular,
the following filter that is reminiscent of Metropolis sampling can be
achieved that privileges transitions to lower energy levels:

$\displaystyle\gamma^{M}(\omega):=e^{-\beta\max\big{(}\omega+\frac{1}{2\beta},0\big{)}}\,.$

Consider the generator ${\mathcal{L}}_{\beta,H}$ constructed through jumps
$\\{A^{a}\\}_{a\in[m]}$ and Gaussian filter $\gamma$ replaced by the
Metropolis filter $\gamma^{M}$ and associated Hamiltonian $H$ on $n^{\prime}$
qubits, with constant minimum energy and Bohr frequency gaps, and
$M=\operatorname{poly}(n^{\prime})$ distinct eigenvalues. In Proposition 8, we
prove that it satisfies the following continuity bound:

$\displaystyle\|{\mathcal{L}}_{\infty,H}-{\mathcal{L}}_{\beta,H}\|_{1\to
1}=\mathcal{O}\Big{(}mM^{4}e^{-\Delta\,\beta}\Big{)}$ (6)

for some constant $\Delta>0$ depending on the spectrum of $H$. In words, the
generator ${\mathcal{L}}_{\infty,H}$ provides a good approximation to
${\mathcal{L}}_{\beta,H}$ as long as $\beta=\Omega(\log(n^{\prime}))$. Next in
Proposition 10, we show the following perturbation bound at fixed $\beta$:
given $H=H_{0}+V$ a perturbed Hamiltonian over $n^{\prime}$ qubits with
$\|V\|_{\infty}$ small enough,

$\displaystyle\|\mathcal{L}_{\beta,H}\\!-\\!\mathcal{L}_{\beta,H_{0}}\|_{1\to
1}\\!=\\!\mathcal{O}\Big{(}m\beta\|V\|_{\infty}\Big{(}\|H_{0}\|_{\infty}\\!-\\!\ln\|V\|_{\infty}\\!\Big{)}\Big{)}.$
(7)

We resort to CTH constructions that are frustration free and have ground state
energy $0$. It then suffices to measure all energy terms and obtain $0$ to
ensure we have prepared the ground state. For this procedure to succeed with
polynomially small probability, it suffices to obtain a state with at least a
polynomial overlap with the ground state. We aim to prove that, when
initialized in the maximally mixed state, the evolution generated by
$\mathcal{L}_{\beta,H}$ reaches rapidly states with at least inverse
polynomial overlap with $|\psi_{0}\rangle$. For this, it is enough to control
the Laplace transform of the random variable corresponding to the measurement
of $H$ in the evolved state $\rho_{t}:=e^{t\mathcal{L}_{\beta,H}}(\rho)$ given
input state $\rho$. Indeed, by a simple probabilistic trick known as Herbst
argument, denoting by $E_{1}$ the energy of the lowest excited states, for any
$\theta>0$

$\displaystyle 1-\langle\psi_{0}|\rho_{t}|\psi_{0}\rangle\leq e^{-\theta
E_{1}}\operatorname{tr}\left(\rho_{t}e^{\theta H}\right)\,.$

In the standard case where
$E_{1}=\Omega(\operatorname{poly}(n^{\prime})^{-1})$ and
$\|H_{0}\|_{\infty}=\mathcal{O}(\operatorname{poly}(n^{\prime}))$, combining
Eq.(6) and Eq.(7), one can replace the state $\rho_{t}$ on the upper bound
above by the state $\rho^{\prime}_{t}:=e^{t\mathcal{L}_{\infty,H_{0}}}(\rho)$,
and $e^{\theta H}$ by $e^{\theta H_{0}}$, at the cost of an inverse
polynomially small additive error, as long as
$t=\mathcal{O}(\operatorname{poly}(n^{\prime}))$ and
$\theta,\|V\|_{\infty}=\mathcal{O}(\operatorname{poly}(n^{\prime})^{-1})$.

We have just argued that, to show that the overlap between $\rho_{t}$ and
$|\psi_{0}\rangle$ is at least inverse polynomial in $n^{\prime}$, it suffices
to show that the Laplace transform
$\operatorname{tr}\left(\rho_{t}^{\prime}e^{\theta H_{0}}\right)$ is inverse
polynomially close to $1$. This is taken care of in Lemma 11. There, we show
that, as long as the Hamiltonian satisfies a Cheeger-like inequality
associated to a certain constant $C>0$, we have

$\displaystyle\operatorname{tr}\left(\rho_{t}^{\prime}e^{\theta
H_{0}}\right)\leq
1+e^{-\frac{C}{2}\theta\Delta_{E}e^{-\Delta_{E}\theta}t+\theta\|H_{0}\|_{\infty}}\,,$

where $\Delta_{E}\equiv\Delta_{E}(H_{0})$ denotes the minimal gap between any
two eigenvalues of $H_{0}$. We then show that the clock Hamiltonian satisfies
this inequality. Finally, we show a perturbed version of the inequality in
Prop. 12 to conclude that the same is true for a suitable CTH, i.e. we can
achieve $\Omega(\textrm{poly}(n)^{-1})$ overlap with the ground state of the
CTH in polynomial time.

### Comparison with previous work

The idea of solving $\mathsf{BQP}$-hard problems by engineering Lindbladians
can be traced back to [48, 49]. However, the Lindbladians constructed there
did not emulate the naturally occurring physical process of thermalization. In
contrast to these early works, the authors of the recent paper [19] developed
a thermal gradient descent algorithm inspired by natural physical processes
that generate low-energy states. This algorithm involves intertwining local
dissipation through Gibbs sampling, followed by measurements to confirm the
reduction in energy at the output of the dissipation. While this new model of
quantum computing is polynomially equivalent to the circuit model, the
measurements at each step necessitate a possibly prohibitive overhead in
practice. Moreover, the algorithm is guaranteed to work only at inverse
temperatures that scale as $\beta=\mathcal{O}(T^{19})$ in the size $T$ of the
circuit being emulated. In contrast, our approach is proved to work already at
$\beta=\mathcal{O}(\log(n))$ for circuits of polynomial size, and does not
require the active choice of locally dissipating depending on intermediate
energy measurements: we simply need to let the system thermalize globally with
a Hamiltonian corresponding to a CTH encoding for a polynomial amount of time.
We then measure a series of local tests corresponding to the terms of another
CTH mapping. If they all accept, which we show happens with inverse
polynomially large probability, we obtain the output of the circuit. Thus,
repeating this procedure polynomially many times suffices to reach the ground
state. In that sense, our strategy is arguably closer in spirit to the quantum
adiabatic algorithm. That said, and while we do not further investigate this
possibility in the present study, the dissipative nature of this procedure
compared to quantum adiabatic computing might provide a natural fault-
tolerance to errors. In contrast, we recall that the proof of a version of the
threshold theorem for adiabatic quantum computing remains a major open problem
to this date [50, 51]. Finally, our proof is arguably much more
straightforward than that of [19].

##  VI Conclusion

In this work, we establish the first results on the efficient dissipative
preparation of high-temperature Gibbs states and their purifications, and
showed that the class of Lindbladian evolutions of [13] is BQP complete
already at moderate temperatures. While the latter should be taken more as a
proof-of-principle result, we hope that further work can improve the
construction. Taken together, our results show that the Lindbladians of [13]
have the full potential to replicate the success of MCMC for quantum many-body
states. It would also be interesting to improve the results presented here by
showing a modified log-Sobolev inequality [10, 9, 8] for the generators
instead of the spectral gap, as this would provide convergence not only in
polynomial but even logarithmic time. Furthermore, implementing the
Lindbladian of [13] still requires sophisticated quantum circuits. It would be
interesting to see if in our regimes it is possible to implement them more
simply, perhaps even through a microscopic model of system-environment
interactions such as the one for Davies maps [52]. This would be particularly
interesting for near-term computation, as it is known that dissipative
preparation can be more robust than simulating unitary dynamics [53, 54, 55].
More generally, it would be interesting to investigate further how errors
affect these time evolutions, and for instance, whether the computational
model GibbsQP admits a threshold theorem.

##### Acknowledgements:

AMA acknowledges useful discussions with Tim Hsieh and David Gosset. CR would
like to thank Simone Warzel for fruitful discussions on Hamiltonian gaps, and
acknowledges financial support from the ANR project QTraj
(ANR-20-CE40-0024-01) of the French National Research Agency (ANR). DSF
acknowledges funding from the European Union under Grant Agreement 101080142
and the project EQUALITY. AMA acknowledges support from the Spanish Agencia
Estatal de Investigacion through the grants “IFT Centro de Excelencia Severo
Ochoa SEV-2016-0597" and “Ramón y Cajal RyC2021-031610-I”, financed by
MCIN/AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR.

## References

  * Levin and Peres [2017] D. A. Levin and Y. Peres,  _Markov chains and mixing times_ , Vol. 107 (American Mathematical Soc., 2017).
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of markov chain monte carlo_ (CRC press, 2011).
  * Martinelli [1999] F. Martinelli, Lectures on probability theory and statistics (Saint-Flour, 1997) 1717, 93 (1999).
  * Temme  _et al._ [2011] K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, Nature 471, 87–90 (2011).
  * Brandão and Kastoryano [2019] F. G. S. L. Brandão and M. J. Kastoryano, Communications in Mathematical Physics 365, 1 (2019).
  * Kastoryano and Brandão [2016] M. J. Kastoryano and F. G. S. L. Brandão, Communications in Mathematical Physics 344, 915 (2016).
  * Bardet  _et al._ [2021] I. Bardet, A. Capel, A. Lucia, D. Pérez-García, and C. Rouzé, Journal of Mathematical Physics 62, 061901 (2021), https://doi.org/10.1063/1.5142186 .
  * Ángela Capel  _et al._ [2021] Ángela Capel, C. Rouzé, and D. S. França, The modified logarithmic sobolev inequality for quantum spin systems: classical and commuting nearest neighbour interactions (2021), arXiv:2009.11817 [quant-ph] .
  * Bardet  _et al._ [2024] I. Bardet, Á. Capel, L. Gao, A. Lucia, D. Pérez-García, and C. Rouzé, Communications in Mathematical Physics 405, 42 (2024).
  * Bardet  _et al._ [2023] I. Bardet, Á . Capel, L. Gao, A. Lucia, D. Pérez-Garc\́text{id}a, and C. Rouzé, Physical Review Letters 130, 10.1103/physrevlett.130.060401 (2023).
  * Zhang  _et al._ [2023] D. Zhang, J. L. Bosse, and T. Cubitt, Dissipative quantum gibbs sampling (2023), arXiv:2304.04526 [quant-ph] .
  * Shtanko and Movassagh [2023] O. Shtanko and R. Movassagh, Preparing thermal states on noiseless and noisy programmable quantum processors (2023), arXiv:2112.14688 [quant-ph] .
  * Chen  _et al._ [2023a] C.-F. Chen, M. J. Kastoryano, and A. Gilyén, An efficient and exact noncommutative quantum gibbs sampler (2023a), arXiv:2311.09207 [quant-ph] .
  * Chen  _et al._ [2023b] C.-F. Chen, M. J. Kastoryano, F. G. S. L. Brandão, and A. Gilyén, Quantum thermal state preparation (2023b), arXiv:2303.18224 [quant-ph] .
  * Nachtergaele and Sims [2006] B. Nachtergaele and R. Sims, Communications in Mathematical Physics 265, 119–130 (2006).
  * Cottrell  _et al._ [2019] W. Cottrell, B. Freivogel, D. M. Hofman, and S. F. Lokhande, Journal of High Energy Physics 2019, 58 (2019), arXiv:1811.11528 [cond-mat, physics:hep-th, physics:quant-ph].
  * Sundar  _et al._ [2022] B. Sundar, A. Elben, L. K. Joshi, and T. V. Zache, New Journal of Physics 24, 023037 (2022), arXiv:2107.02196 [cond-mat, physics:nlin, physics:quant-ph].
  * Brown  _et al._ [2023] A. R. Brown, H. Gharibyan, S. Leichenauer, H. W. Lin, S. Nezami, G. Salton, L. Susskind, B. Swingle, and M. Walter, PRX Quantum 4, 010320 (2023), publisher: American Physical Society.
  * Chen  _et al._ [2023c] C.-F. Chen, H.-Y. Huang, J. Preskill, and L. Zhou, Local minima in quantum systems (2023c), arXiv:2309.16596 [quant-ph] .
  * Michalakis and Zwolak [2013] S. Michalakis and J. P. Zwolak, Communications in Mathematical Physics 322, 277 (2013).
  * Caputo  _et al._ [2011] P. Caputo, F. Martinelli, F. Simenhaus, and F. L. Toninelli, Communications on Pure and Applied Mathematics 64, 778 (2011).
  * Bravyi  _et al._ [2010] S. Bravyi, M. B. Hastings, and S. Michalakis, Journal of Mathematical Physics 51, 10.1063/1.3490195 (2010).
  * [23] See the Supplemental Material for the detailed proofs of the main results.
  * Hastings and Koma [2006] M. B. Hastings and T. Koma, Communications in Mathematical Physics 265, 781 (2006).
  * Tran  _et al._ [2021] M. C. Tran, A. Y. Guo, C. L. Baldwin, A. Ehrenberg, A. V. Gorshkov, and A. Lucas, Phys. Rev. Lett. 127, 160401 (2021).
  * [26] C.-F. (Anthony) Chen, A. Lucas, and C. Yin, Reports on Progress in Physics 86, 116001 (2023).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Physical Review Letters 103, 220502 (2009), arXiv:0905.2199 [quant-ph].
  * Chowdhury and Somma [2017] A. N. Chowdhury and R. D. Somma, Quantum Information and Computation 17, 10.26421/QIC17.1-2 (2017).
  * Holmes  _et al._ [2022] Z. Holmes, G. Muraleedharan, R. D. Somma, Y. Subasi, and B. Şahinoğlu, Quantum 6, 825 (2022).
  * Chowdhury  _et al._ [2020] A. N. Chowdhury, G. H. Low, and N. Wiebe, A Variational Quantum Algorithm for Preparing Quantum Gibbs States (2020), arXiv:2002.00055 [quant-ph].
  * Wang  _et al._ [2021] Y. Wang, G. Li, and X. Wang, Physical Review Applied 16, 054035 (2021), publisher: American Physical Society.
  * Chen and Brandão [2023] C.-F. Chen and F. G. S. L. Brandão, Fast thermalization from the eigenstate thermalization hypothesis (2023), arXiv:2112.07646 [quant-ph] .
  * Li and Wang [2023] X. Li and C. Wang, in  _50th International Colloquium on Automata, Languages, and Programming (ICALP 2023)_, Leibniz International Proceedings in Informatics (LIPIcs), Vol. 261, edited by K. Etessami, U. Feige, and G. Puppis (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2023) pp. 87:1–87:20.
  * Rall  _et al._ [2023] P. Rall, C. Wang, and P. Wocjan, Quantum 7, 1132 (2023).
  * Bilgin and Boixo [2010] E. Bilgin and S. Boixo, Physical Review Letters 105, 170405 (2010).
  * Ge  _et al._ [2016] Y. Ge, A. Molnár, and J. I. Cirac, Phys. Rev. Lett. 116, 080503 (2016).
  * Feiguin and Klich [2013] A. E. Feiguin and I. Klich, Hermitian and non-Hermitian thermal Hamiltonians (2013), arXiv:1308.0756 [cond-mat, physics:physics].
  * Ambainis and Regev [2006] A. Ambainis and O. Regev, An elementary proof of the quantum adiabatic theorem (2006), arXiv:quant-ph/0411152 [quant-ph] .
  * Maldacena and Qi [2018] J. Maldacena and X.-L. Qi, Eternal traversable wormhole (2018), arXiv:1804.00491 [cond-mat, physics:gr-qc, physics:hep-th].
  * Nezami  _et al._ [2023] S. Nezami, H. W. Lin, A. R. Brown, H. Gharibyan, S. Leichenauer, G. Salton, L. Susskind, B. Swingle, and M. Walter, PRX Quantum 4, 010321 (2023), arXiv:2102.01064 [hep-th, physics:quant-ph].
  * Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Physical Review Letters 123, 220502 (2019).
  * Martyn and Swingle [2019] J. Martyn and B. Swingle, Physical Review A 100, 032107 (2019).
  * Faílde  _et al._ [2023] D. Faílde, J. Santos-Suárez, D. A. Herrera-Martí, and J. Mas, Hamiltonian Forging of a Thermofield Double (2023).
  * Zhu  _et al._ [2020] D. Zhu, S. Johri, N. M. Linke, K. A. Landsman, N. H. Nguyen, C. H. Alderete, A. Y. Matsuura, T. H. Hsieh, and C. Monroe, Proceedings of the National Academy of Sciences 117, 25402 (2020).
  * Aharonov  _et al._ [2008] D. Aharonov, W. van Dam, J. Kempe, Z. Landau, S. Lloyd, and O. Regev, SIAM Review 50, 755 (2008), https://doi.org/10.1137/080734479 .
  * Kempe  _et al._ [2006] J. Kempe, A. Kitaev, and O. Regev, SIAM Journal on Computing 35, 1070 (2006), https://doi.org/10.1137/S0097539704445226 .
  * Oliveira and Terhal [2008] R. Oliveira and B. M. Terhal, Quantum Info. Comput. 8, 900–924 (2008).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nature Physics 5, 633–636 (2009).
  * Kraus  _et al._ [2008] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Childs  _et al._ [2001] A. M. Childs, E. Farhi, and J. Preskill, Phys. Rev. A 65, 012322 (2001).
  * Young  _et al._ [2013] K. C. Young, M. Sarovar, and R. Blume-Kohout, Phys. Rev. X 3, 041013 (2013).
  * Davies [1974] E. B. Davies, Communications in Mathematical Physics 39, 91 (1974).
  * Cubitt  _et al._ [2015] T. S. Cubitt, A. Lucia, S. Michalakis, and D. Perez-Garcia, Communications in Mathematical Physics 337, 1275–1315 (2015).
  * Borregaard  _et al._ [2021] J. Borregaard, M. Christandl, and D. Stilck França, npj Quantum Information 7, 10.1038/s41534-021-00363-9 (2021).
  * Kim and Swingle [2017] I. H. Kim and B. Swingle, Robust entanglement renormalization on a noisy quantum computer (2017), arXiv:1711.07500 [quant-ph] .
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, Communications in Mathematical Physics 28, 251 (1972).
  * Hastings [2004] M. B. Hastings, Phys. Rev. B 69, 104431 (2004).
  * Haah  _et al._ [2021] J. Haah, M. B. Hastings, R. Kothari, and G. H. Low, SIAM Journal on Computing 52, FOCS18 (2021).
  * Kitaev  _et al._ [2002] A. Y. Kitaev, A. Shen, and M. N. Vyalyi,  _Classical and quantum computation_ , 47 (American Mathematical Soc., 2002).
  * Note [1] In fact, the clock Hamiltonian constructed in [19] is slightly different from the one we consider here, however it is easy to verify that their ground states agree by a simple modification of the proof given in [19, Appendix J.1].
  * Note [2] In principle they only need to be local w.r.t. a fully connected graph, so that we can map any quantum circuit to this graph with a polynomial overhead.

##  Appendix A Basic notations and setup

Given a finite set $V$, we denote by $\mathcal{H}_{V}=\bigotimes_{v\in
V}\mathcal{H}_{v}$ the Hilbert space of $n=|V|$ qubits (_i.e._ ,
$\mathcal{H}_{v}\equiv\mathbb{C}^{2}$ for all $v\in V$) and by
$\mathcal{B}_{V}$ the algebra of linear operators on $\mathcal{H}_{V}$. The
trace on $\mathcal{B}_{V}$ is denoted by
$\operatorname{tr}\left(\cdot\right)$. $\mathcal{O}_{V}$ corresponds to the
space of self-adjoint linear operators on $\mathcal{H}_{V}$. $\mathcal{S}_{V}$
denotes the set of quantum states. For any subset $A\subseteq V$, we use the
standard notations $\mathcal{O}_{A},\mathcal{S}_{A}\ldots$ for the
corresponding objects defined on subsystem $A$. For systems made of $n$
qudits, we sometimes also use the notations $[n]$ when refering to the set
$\\{1,\dots,n\\}$. Given a state $\rho\in\mathcal{S}_{V}$, we denote by
$\rho_{A}$ its marginal onto the subsystem $A$. For any $X\in\mathcal{O}_{V}$,
we denote by $\|X\|_{1}:=\operatorname{tr}\left(|X|\right)$ its trace norm, by
$\|X\|_{\infty}:=\sup_{\||\psi\rangle\|\leq 1}\|X|\psi\rangle\|$ its operator
norm, and by $\|X\|_{2}:=\sqrt{\operatorname{tr}\left(|X|^{2}\right)}$ its
Hilbert-Schmidt norm. Following standard quantum information-theoretic
terminology, operators acting on susystems $A\subseteq V$ are often denoted by
$O_{A}$. The identity map is denoted by $I$. We denote by $O^{\dagger}$ the
adjoint of an operator $O\in\mathcal{B}_{V}$, and by $\Phi^{\dagger}$ the
Hilbert-Schmidt dual of a superoperator on $\mathcal{B}_{V}$.

###  A.1 Local Hamiltonians and Lieb-Robinson bounds

Next, we consider a quantum spin system on a $D$-dimensional finite lattice
$\Lambda\equiv[0,L]^{D}$, with $n:=|\Lambda|=(L+1)^{D}$, and a local
Hamiltonian $H$ on the lattice, defined through a function that maps any non-
empty finite set $X\subset\Lambda$ to a self-adjoint element $h_{X}$ supported
in $X$ with $\max_{X}\left\|h_{X}\right\|_{\infty}\leq h$. The Hamiltonian of
the system in a region $A\subset\Lambda$ is then defined as the sum over
interactions included in $A$:

$\displaystyle H_{A}:=\sum_{X\subseteq A}h_{X}\,.$

In addition, we assume that the Hamiltonian is such that no interaction
$h_{X}$ has support on more than $k$ sites, and each site $u$ appears on at
most $l$ non-zero operators $h_{X}$. We refer to such Hamiltonian as
$(k,l)$-local. As often, we will make the implicit assumption that $H$ is one
element of a family of Hamiltonians indexed by the lattice size, in which case
the constants $h,k$ and $l$ are assumed to be independent of $L$.

For such Hamiltonians, given the so-called Lieb-Robinson velocity
$J:=\max_{u\in\Lambda}\sum_{Z\ni u}|Z|\left\|h_{Z}\right\|_{\infty}$, it is
well-known that, for any operator $A^{u}$ supported on site $u\in\Lambda$ [56,
57, 15, 24, 58]:

$\left\|e^{-iHt}A^{u}e^{iHt}-e^{-iH_{B_{u}(r)}t}A^{u}e^{iH_{B_{u}(r)}t}\right\|_{\infty}\leq\left\|A^{u}\right\|_{\infty}\frac{(2J|t|)^{r}}{r!}\leq\left\|A^{u}\right\|_{\infty}\frac{(2hkl|t|)^{r}}{r!}\,,$
(8)

where $B_{u}(r)$ denotes the ball of radius $r$ with respect to the graph
distance centered at $u$, and where the second bound results from
$(k,l)$-locality.

###  A.2 Quantum Gibbs sampling with Gaussian filters

We now recall the generator of the quantum Monte Carlo algorithm introduced in
[13]: for $\beta>0$ and a Hamiltonian $H$, we denote
$\gamma(\omega):=\operatorname{exp}(-\frac{(\beta\omega+1)^{2}}{2})$ and we
consider the Lindbladian defined for all $\rho\in\mathcal{S}_{\Lambda}$ as

$\displaystyle\mathcal{L}_{\beta}(\rho)$
$\displaystyle=-i[B,\rho]+{\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\int_{-\infty}^{\infty}\gamma(\omega)\,\Big{(}A^{a,\alpha}(\omega)\rho
A^{a,\alpha}(\omega)^{\dagger}-\frac{1}{2}\big{\\{}A^{a,\alpha}(\omega)^{\dagger}A^{a,\alpha}(\omega),\rho\big{\\}}\Big{)}d\omega}\equiv\sum_{a}\sum_{\alpha}\mathcal{L}_{\beta}^{a,\alpha}(\rho)\,,$
(9)

with jump operators

$\displaystyle
A^{a,\alpha}(\omega):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{iHt}A^{a,\alpha}e^{-iHt}e^{-i\omega
t}\,f(t)\,dt\qquad\text{ where }\quad
f(t):=e^{-\frac{t^{2}}{\beta^{2}}}\,\sqrt{\beta^{-1}\sqrt{2/\pi}}\,.$

Above, $A^{a,1}\equiv\sigma_{x},A^{a,2}\equiv\sigma_{y}$ and
$A^{a,3}\equiv\sigma_{z}$ denote the $1$-local Pauli matrices on site $a$. The
coherent part is defined through the self-adjoint matrix

$\displaystyle
B:=\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}b_{1}(t)e^{-i\beta
Ht}\int_{-\infty}^{\infty}b_{2}(t^{\prime})e^{i\beta
Ht^{\prime}}A^{a,\alpha}e^{-2i\beta Ht^{\prime}}A^{a,\alpha}e^{i\beta
Ht^{\prime}}dt^{\prime}e^{i\beta
Ht}dt\equiv\sum_{a}\sum_{\alpha}B_{\beta}^{a,\alpha}\,,$

for some smooth, rapidly decaying functions $b_{1},b_{2}$ with
$\|b_{1}\|_{L_{1}},\|b_{2}\|_{L_{1}}\leq 1$ defined as

$\displaystyle
b_{1}(t):=2\sqrt{\pi}\,e^{\frac{1}{8}}\,\left(\frac{1}{\operatorname{cosh}(2\pi
t)}\ast_{t}\,\sin(-t)e^{-2t^{2}}\right)$ (10) $\displaystyle
b_{2}(t):=\frac{1}{2\pi}\,\sqrt{\frac{1}{\pi}}\,\operatorname{exp}\big{(}-4t^{2}-2it\big{)}\,.$
(11)

The Lindbladian constructed above has the main advantage of being both quasi-
local and reversible (i.e. in “detailed balance”) in the Heisenberg picture:

$\displaystyle\langle
X,\mathcal{L}^{\dagger}_{\beta}(Y)\rangle_{\sigma_{\beta}}=\langle\mathcal{L}^{\dagger}_{\beta}(X),Y\rangle_{\sigma_{\beta}}\,,$

where $\langle
X,Y\rangle_{\sigma_{\beta}}:=\operatorname{tr}\left(X^{\dagger}{\sigma_{\beta}}^{\frac{1}{2}}Y{\sigma_{\beta}}^{\frac{1}{2}}\right)$
denotes the KMS inner product, and where the Gibbs state

$\displaystyle\sigma_{\beta}=\frac{e^{-\beta
H}}{\operatorname{tr}\left(e^{-\beta H}\right)}$

is the unique invariant state of the evolution generated by
$\mathcal{L}_{\beta}$.

In what follows, we will make use of a well-known mapping of KMS-symmetric
Lindbladians into Hamiltonians: for any $\beta>0$, we define the operator

$\displaystyle\widetilde{\mathcal{L}}_{\beta}(X):=\sigma_{\beta}^{-\frac{1}{4}}\mathcal{L}_{\beta}(\sigma_{\beta}^{\frac{1}{4}}X\sigma_{\beta}^{\frac{1}{4}})\sigma_{\beta}^{-\frac{1}{4}}\equiv\sum_{a,\alpha}\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}\,.$
(12)

Clearly, since $\sigma_{\beta}$ is a steady state of the evolution, and unique
due to the choice of the jumps $A^{a,\alpha}$,
$\operatorname{Ker}(\widetilde{\mathcal{L}}_{\beta})=|\sqrt{\sigma_{\beta}}\rangle$.
Moreover, $\widetilde{\mathcal{L}}_{\beta}$ is self-adjoint with respect to
the Hilbert-Schmidt inner product for all $\beta>0$, by the KMS symmetry of
$\mathcal{L}_{\beta}^{\dagger}$. Our strategy in what follows will be to show
that $\widetilde{\mathcal{L}}_{\beta}$ can be expressed as a quasi-local
perturbation of a certain $1$-local gapped Hamiltonian (the limit
$\beta\rightarrow 0$) and to study the stability of the latter against such
perturbations.

##  Appendix B Constant gap at high temperature

In this appendix, we show our first main result, namely that the Lindbladians
$\mathcal{L}_{\beta}^{\dagger}$ are gapped at high enough temperature:

######  Theorem 1.

In the notations of Appendix A, for any $(k,l)$-local Hamiltonian $H$, there
exists a constant $\beta^{*}>0$ independent of $L$ such that for
$\beta<\beta^{*}$, the spectral gap of ${\mathcal{L}}^{\dagger}_{\beta}$, is
lower bounded by $\frac{1}{2\sqrt{2}e^{1/4}}$.

In general, proving that a dissipative evolution is gapped is not sufficient
to ensure the rapid convergence towards its fixed points, as the Liouvillian
might have exponentially large Jordan blocks. However the situation turns out
to be easier in the case of a Gibbs sampler.

######  Corollary 2.

In the setting of Theorem 1, for any $\beta<\beta^{*}$ and $\epsilon>0$, the
evolution generated by $\mathcal{L}_{\beta}$ gets $\epsilon$-close to
$\sigma_{\beta}$ in polynomial time. More precisely, for any
$\rho\in\mathcal{S}_{\Lambda}$,

$\displaystyle\|e^{t\mathcal{L}_{\beta}}(\rho)-\sigma_{\beta}\|_{1}\leq\epsilon\qquad\text{
for all}\qquad t=\Omega(\log(1/\epsilon)+|\Lambda|)\,.$

###### Proof.

The corollary above is a consequence of a very standard result, whose proof we
recall for the reader’s convenience. By equivalence of noncommutative weighted
$L_{p}$ norms and the KMS symmetry of $\mathcal{L}_{\beta}^{\dagger}$, the
trace distance can be controlled as follows: given
$X:=\sigma_{\beta}^{-{1}/{2}}\rho\sigma_{\beta}^{-{1}/{2}}$,

$\displaystyle\|e^{t\mathcal{L}_{\beta}}(\rho)-\sigma_{\beta}\|_{1}\leq\|e^{t\mathcal{L}^{\dagger}_{\beta}}(X)-I\|_{\sigma_{\beta}}\leq
e^{-\frac{t}{2\sqrt{2}e^{1/4}}}\|X-I\|_{\sigma_{\beta}}\leq
2\|\sigma_{\beta}^{-1}\|_{\infty}e^{-\frac{t}{2\sqrt{2}e^{1/4}}}$
$\displaystyle\leq
2^{|\Lambda|+1}\,e^{2\beta^{*}\|H\|_{\infty}-\frac{t}{2\sqrt{2}e^{1/4}}}$
$\displaystyle\leq
2^{|\Lambda|+1}\,e^{2\beta^{*}lh|\Lambda|-\frac{t}{2\sqrt{2}e^{1/4}}}\,,$

where
$\left\|Y\right\|_{\sigma_{\beta}}=\left\|\sigma_{\beta}^{1/4}Y\sigma_{\beta}^{1/4}\right\|_{2}$.
The result follows by setting the upper bound above to $\epsilon$ and solving
for $t$.

∎

The rest of this appendix is devoted to proving Theorem (1) and to applying it
to the adiabatic algorithm. In Section B.1, we recall some well-known gap
stability results for the class of frustration-free gapped Hamiltonians
satisfying a condition known as local topological quantum order. In Section
B.2, we prove that the maps $\widetilde{\mathcal{L}}_{\beta}$ satisfy the
conditions for stability of the gap recalled in Section B.1. Finally, we
analyze the efficiency of the adiabatic evolution associated to the generator
$\widetilde{\mathcal{L}}_{\beta}$ in Section B.3.

###  B.1 Stability of frustration-free gapped Hamiltonians

As previously mentioned, our main result at high temperatures builds on a
well-known result [22, 20] about the stability of certain gapped Hamitonians
under small quasi-local perturbations: let $\widetilde{H}_{0}$ be a
Hamiltonian defined on a $D$-dimensional lattice $\Lambda$ of side-length
$L+1$ with corresponding Hilbert space
$\widetilde{\mathcal{H}}_{\Lambda}=\bigotimes_{j\in\Lambda}\widetilde{\mathcal{H}}_{j}$,
which satisfies the following requirements:

  * 1.

It can be written as a sum of geometrically local terms
$\widetilde{H}_{0}=\sum_{u\in\Lambda}\widetilde{Q}_{u}$, where each
interaction $\widetilde{Q}_{u}$ is positive semidefinite and acts non-
trivially on a Hilbert space supported on the ball $B_{u}(1)$ of radius $1$
around site $u$;

  * 2.

It satisfies periodic boundary conditions;

  * 3.

For $\widetilde{P}_{0}$ the projector onto the groundspace of
$\widetilde{H}_{0}$, we have $\widetilde{H}_{0}\widetilde{P}_{0}=0$, i.e. the
ground state energy is fixed to $0$. Moreover,
$\widetilde{Q}_{u}\widetilde{P}_{0}=0$ for all $u\in\Lambda$ (frustration-
freeness);

  * 4.

$\widetilde{H}_{0}$ has a gap $\gamma_{L}\geq\gamma>0$ above the groundstate
subspace, for all $L\geq 2$, where $\gamma$ is a constant independent of the
system size $L$.

Given a ball $A\equiv B_{u}(r)$ with $r\leq L^{*}\leq L$ and an observable
$O_{A}$ supported on $A$, we set $A(\ell):=B_{u}(r+\ell)$ and denote

$\displaystyle
c_{\ell}(O_{A}):=\frac{\operatorname{tr}\left(\widetilde{P}_{A(\ell)}O_{A}\right)}{\operatorname{tr}\left(\widetilde{P}_{A(\ell)}\right)}\,,\qquad
c(O_{A}):=\frac{\operatorname{tr}\left(\widetilde{P}_{0}O_{A}\right)}{\operatorname{tr}\left(\widetilde{P}_{0}\right)}\,,$

where $\widetilde{P}_{B}(\epsilon)$ is the projection onto the subspace of
eigenstates of $\widetilde{H}_{B}=\sum_{B_{v}(1)\subset B}\widetilde{Q}_{v}$
with energy at most $\epsilon\geq 0$ and we set
$\widetilde{P}_{B}(0)\equiv\widetilde{P}_{B}$. The Hamiltonian
$\widetilde{H}_{0}$ satisfies local topological order if, for each fixed
$1\leq\ell\leq L-r$,

$\displaystyle\Big{\|}P_{A(\ell)}O_{A}P_{A(\ell)}-c_{\ell}(O_{A})P_{A(\ell)}\Big{\|}\leq\|O_{A}\|\,\Delta_{0}(\ell)\,.$

where $\Delta_{0}(\ell)$ is a decaying function of $\ell$. Here $L^{*}$
denotes a cut-off parameter for topological order, which corresponds to the
code distance in the case of a stabilizer Hamiltonian. Finally, we say that
$\widetilde{H}_{0}$ is locally gapped with gap $\gamma(r)>0$ iff for each
$u\in\Lambda$ and $r\geq 0$,
$\widetilde{P}_{B_{u}(r)}(\gamma(r))=\widetilde{P}_{B_{u}(r)}$. When
$\gamma(r)$ decays at most polynomially in $r$, we say that
$\widetilde{H}_{0}$ satisfies the local gap condition.

Next, we consider a perturbation $\widetilde{V}$ with strength $\widetilde{J}$
and decay rate $f(r)$, with $f(r)\leq 1$, $r\geq 0$, that is,
$\widetilde{V}:=\sum_{u\in\Lambda}\sum_{r=0}^{L}\widetilde{V}_{u}(r)$, such
that $\widetilde{V}_{u}(r)$ has support on $B_{u}(r)$ and satisfies
$\|\widetilde{V}_{u}(r)\|\leq\widetilde{J}f(r)$ for some rapidly decaying
function $f(r)$ with $f(r)\leq 1$.

######  Theorem 6 (Stability of gap for frustration-free gapped Hamiltonians
[20]).

Let $\widetilde{H}_{0}$ be a frustration-free Hamiltonian as defined above
with spectral gap $\gamma$ satisfying the local topological quantum order and
local-gap conditions. For a perturbation $\widetilde{V}$ of strength
$\widetilde{J}$ and decay rate $f(r)$ there exist constants
$\widetilde{J}_{0}>0$ and $L_{0}\geq 2$ such that for
$\widetilde{J}\leq\widetilde{J}_{0}$ and $L\geq L_{0}$, the spectral gap of
$\widetilde{H}_{0}+\widetilde{V}$ is bounded from below by $\gamma/2$.

The constants $\widetilde{J}_{0}$ and $\widetilde{L}_{0}$ depend on the
unperturbed Hamiltonian $\widetilde{H}_{0}$ and the form of the decay $f(r)$.
More specifically, they depend on the Lieb-Robinson velocity of
$\widetilde{H}_{0}+\widetilde{V}$, the global $\gamma$ and local $\gamma(r)$
gaps, and the volume of a $D$-dimensional sphere if radius $1$ (so that
$\widetilde{J}_{0}$ decreases for higher dimensional lattices). Their specific
form can be determined from [20].

###  B.2 Quantum Monte Carlo algorithm as a gapped Hamiltonian: proof of
Theorem 1

We now show that the Hamiltonian $\widetilde{\mathcal{L}}_{\beta}$ defined in
(12) on the Hilbert space of Hilbert-Schmidt operators over
$(\mathbb{C}^{2})^{\otimes\Lambda}$ can be seen as the perturbation
$\widetilde{H}_{0}+\widetilde{V}$ of a Hamiltonian $\widetilde{H}_{0}$
(corresponding to $\beta\rightarrow 0$) that satisfies the conditions of
Theorem 6 for $\beta$ small enough, from which we infer that it has a finite
spectral gap. To do this, we study the convergence of
$\widetilde{\mathcal{L}}_{\beta}$ as $\beta\to 0$ in norm $\|.\|_{2\to 2}$. We
first consider the following telescopic sum decomposition of
$\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}$:

$\displaystyle\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}=\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,0}+\sum_{r=0}^{\infty}\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r+1}-\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r},$

where for any $r\in\mathbb{N}$, $\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r}$
is defined with respect to the jump operators $A^{a,\alpha}_{r}(\omega)$ as
well as the state $\sigma_{r}$ expressed analogously to $A^{a,\alpha}(\omega)$
and $\sigma_{\beta}$ up to replacing $H$ with the reduced Hamiltonian
$H_{B_{a}(r)}$. Therefore the perturbations
$\widetilde{\mathcal{E}}^{a,\alpha}_{\beta,r}:=\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r+1}-\widetilde{\mathcal{L}}^{a,\alpha}_{\beta,r}$
are supported on regions $B_{a}(r+1)$. The $r=0$ case corresponds to the
possible $1$-local, single qubit, terms in the Hamiltonian.

Considering the definition of the Lindbladian in Sec. A, and the respective
Hamiltonian in Eq. (12), we can decompose the maps
$\widetilde{\mathcal{L}}_{\beta,r}^{a,\alpha}$ into

$\displaystyle\widetilde{\mathcal{L}}_{\beta,r}^{a,\alpha}:=\widetilde{\mathcal{C}}_{\beta,r}^{a,\alpha}+\widetilde{\Psi}_{\beta,r}^{a,\alpha}+\widetilde{\mathcal{D}}_{\beta,r}^{a,\alpha}\,,$

where

$\displaystyle\widetilde{\mathcal{C}}_{\beta,r}^{a,\alpha}(X):=-i\Big{[}\Delta_{{\beta,r}}^{-\frac{1}{4}}(B^{a,\alpha}_{\beta,r})X-X\Delta_{\beta,r}^{\frac{1}{4}}(B^{a,\alpha}_{\beta,r})\Big{]}\,,$
$\displaystyle\widetilde{\mathcal{D}}_{\beta,r}^{a,\alpha}(X):=-\frac{1}{2}\int\Big{[}\Delta_{\beta,r}^{-\frac{1}{4}}\big{(}A^{a,\alpha}_{r}(\omega)^{\dagger}A^{a,\alpha}_{r}(\omega)\big{)}X+X\Delta_{\beta,r}^{\frac{1}{4}}\big{(}A^{a,\alpha}_{r}(\omega)^{\dagger}A^{a,\alpha}_{r}(\omega)\big{)}\Big{]}\gamma(\omega)d\omega\,,$
$\displaystyle\widetilde{\Psi}_{\beta,r}^{a,\alpha}(X):=\int\gamma(\omega)\,\Delta_{\beta,r}^{-\frac{1}{4}}(A^{a,\alpha}_{r}(\omega))X\Delta_{\beta,r}^{\frac{1}{4}}(A^{a,\alpha}_{r}(\omega)^{\dagger})\,d\omega\,.$

Above, we denoted by $\Delta_{\beta,r}^{z}(X)=\sigma_{r}^{z}X\sigma_{r}^{-z}$
the $z$-th power of the modular operator corresponding to the state
$\sigma_{r}\propto{\operatorname{exp}({-\beta{H_{B_{a}(r)}}})}$. We aim at
controlling the norm $\|\widetilde{\mathcal{E}}^{a,\alpha}_{\beta,r}\|_{2\to
2}$, showing that it decays sufficiently fast with $r$. From Proposition A.1
and Corollary A.5 in [13] we can see that

$\displaystyle\widetilde{\mathcal{C}}_{\beta,r}^{a,\alpha}(X)+\widetilde{\mathcal{D}}_{\beta,r}^{a,\alpha}(X)=\frac{1}{2}\left(N_{\beta,r}^{a,\alpha}X+XN_{\beta,r}^{a,\alpha}\right),$
(13)

where we define

$\displaystyle
N_{\beta,r}^{a,\alpha}=\int_{-\infty}^{\infty}n_{1}(t)e^{-i\beta
tH_{B_{a}(r)}}\left(\int_{-\infty}^{\infty}n_{2}(t^{\prime})A_{r}^{a,\alpha\,\dagger}(\beta
t^{\prime})A_{r}^{a,\alpha}(-\beta t^{\prime})dt^{\prime}\right)e^{i\beta
tH_{B_{a}(r)}}dt,$ (14)

with $A_{r}^{a,\alpha}(t)=e^{itH_{B_{a}(r)}}A^{a,\alpha}e^{-itH_{B_{a}(r)}}$
and

$\displaystyle
n_{1}(t):=\frac{\sqrt{\pi}}{2}\,\left(\frac{1}{\operatorname{cosh}(2\pi
t)}\ast_{t}\,e^{-2t^{2}}\right)$ (15) $\displaystyle
n_{2}(t):=8b_{2}(t)=\frac{4}{\pi^{\frac{3}{2}}}\,\operatorname{exp}\big{(}-4t^{2}-2it\big{)}\,.$
(16)

Notice that these are slightly different to the functions in (10) above due to
the mapping in Eq. (12). We can thus bound

$\displaystyle\left\|\widetilde{\mathcal{C}}_{\beta,r+1}^{a,\alpha}(X)+\widetilde{\mathcal{D}}_{\beta,r+1}^{a,\alpha}(X)-\widetilde{\mathcal{C}}_{\beta,r}^{a,\alpha}(X)+\widetilde{\mathcal{D}}_{\beta,r}^{a,\alpha}(X)\right\|_{2}\leq\left\|X\right\|_{2}\left\|N_{\beta,r+1}^{a,\alpha}-N_{\beta,r}^{a,\alpha}\right\|_{\infty},$

which is such that

$\displaystyle\left\|N_{\beta,r+1}^{a,\alpha}-N_{\beta,r}^{a,\alpha}\right\|_{\infty}\leq\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}n_{1}(t)|n_{2}(t^{\prime})|\left\|\mathcal{N}(t,t^{\prime})\right\|_{\infty}dt^{\prime}dt$

where, by the Lieb-Robinson bound recalled in Equation (8), since
$\|A^{a,\alpha}\|_{\infty}=1$,

$\displaystyle\left\|\mathcal{N}(t,t^{\prime})\right\|_{\infty}=\left\|A_{r+1}^{a,\alpha\,\dagger}(\beta(t^{\prime}-t))A_{r+1}^{a,\alpha}(\beta(-t^{\prime}-t))-A_{r}^{a,\alpha\,\dagger}(\beta(t^{\prime}-t))A_{r}^{a,\alpha}(\beta(-t-t^{\prime}))\right\|_{\infty}\leq
2\frac{\left(\beta J(|t|+|t^{\prime}|)\right)^{r}}{r!},$ (17)

so that, integrating over $t,t^{\prime}$ and upper bounding the integrals, we
find

$\displaystyle\left\|N_{\beta,r+1}^{a,\alpha}-N_{\beta,r}^{a,\alpha}\right\|_{\infty}\leq
2e^{\frac{3\pi^{2}}{4}}\left(\frac{\beta J}{4}\right)^{r}.$

Now we bound the transition part $\widetilde{\Psi}_{\beta,r}^{a,\alpha}(X)$.
As can be seen in [13] (Corollary A.3), we can write

$\displaystyle\widetilde{\Psi}_{\beta,r}^{a,\alpha}(X)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}h_{-}(t_{-})h_{+}(t_{+})A_{r}^{a,\alpha}(\beta(t_{+}-t_{-}))XA_{r}^{a,\alpha}(\beta(-t_{+}-t_{-}))^{\dagger}dt_{+}dt_{-},$

where

$\displaystyle h_{+}(t)=e^{-1/4-4t^{2}}\quad,\quad
h_{-}(t)=\frac{1}{\pi}e^{-2t^{2}}.$

Considering that, again by the Lieb-Robinson bound (8),

$\displaystyle\left\|A_{r+1}^{a,\alpha}(\beta(t_{+}-t_{-}))XA_{r+1}^{a,\alpha}(\beta(-t_{+}-t_{-}))^{\dagger}-A_{r}^{a,\alpha}(\beta(t_{+}-t_{-}))XA_{r}^{a,\alpha}(\beta(-t_{+}-t_{-}))^{\dagger}\right\|_{2}\leq\left\|X\right\|_{2}2\frac{\left(\beta
J\,(|t_{+}|+|t_{-}|)\right)^{r}}{r!},$

we can bound

$\displaystyle\left\|\widetilde{\Psi}_{\beta,r+1}^{a,\alpha}(X)-\widetilde{\Psi}_{\beta,r}^{a,\alpha}(X)\right\|_{2}\leq\left\|X\right\|_{2}\frac{\left(\frac{\sqrt{3}\beta
J}{4}\right)^{r}}{\Gamma(1+\frac{r}{2})}\,.$

Putting all the bounds together, we can write,

$\displaystyle\|\widetilde{\mathcal{E}}_{\beta,r}^{a,\alpha}\|_{2\to 2}\leq
4e^{\frac{3\pi^{2}}{4}}\left(\frac{\beta
J}{4}\right)^{r}+\frac{\left(\frac{\sqrt{3}\beta
J}{4}\right)^{r}}{\Gamma(1+\frac{r}{2})}\,.$ (18)

The next step is to determine the point we perturb around
$\widetilde{\mathcal{L}}^{a,\alpha}_{0,0}$. Notice that
$\lim_{\beta\rightarrow 0}N_{\beta,r}^{a,\alpha}=\lambda{I}$ for some constant
$\lambda$, and so $\lim_{\beta\rightarrow
0}\widetilde{\mathcal{C}}_{\beta,r}^{a,\alpha}(X)+\widetilde{\mathcal{D}}_{\beta,r}^{a,\alpha}(X)=\lambda
X$, which is the identity channel up to $\lambda$. Additionally,

$\displaystyle\lim_{\beta\rightarrow
0}\widetilde{\Psi}_{\beta,0}^{a,\alpha}(X)=\frac{1}{2\sqrt{2}e^{1/4}}A^{a,\alpha}XA^{a,\alpha\,\dagger}.$

This means that $\lambda=-\frac{1}{\sqrt{2}e^{1/4}}$. When summing over all
Pauli matrices $A^{a,\alpha}$, we obtain the generator of the local
depolarizing channel. Finally, we bound the remaining contribution, which is
the norm difference
$\|\widetilde{\mathcal{L}}_{\beta,0}^{a,\alpha}-\widetilde{\mathcal{L}}_{0,0}^{a,\alpha}\|_{2\to
2}$ between the $1$-qubit generator at $\beta$ and that at $\beta=0$. Note
that

$\displaystyle\left\|N_{\beta,0}^{a,\alpha}-\lambda{I}\right\|_{\infty}$
$\displaystyle\leq\int_{-\infty}^{\infty}n_{1}(t)\int_{-\infty}^{\infty}|n_{2}(t^{\prime})|\left\|A_{0}^{a,\alpha\,\dagger}(\beta
t^{\prime})A_{0}^{a,\alpha}(-\beta
t^{\prime})-{I}\right\|_{\infty}dt^{\prime}dt$ (19)
$\displaystyle\leq\int_{-\infty}^{\infty}n_{1}(t)\int_{-\infty}^{\infty}|n_{2}(t^{\prime})|\left\|\sum_{m,m^{\prime}}(-1)^{m^{\prime}}(\beta
t^{\prime})^{m+m^{\prime}}\frac{K_{m}(A_{0}^{a,\alpha\,\dagger})K_{m^{\prime}}(A_{0}^{a,\alpha})}{m!m^{\prime}!}-{I}\right\|_{\infty}dt^{\prime}dt$
$\displaystyle\leq\int_{-\infty}^{\infty}n_{1}(t)\int_{-\infty}^{\infty}|n_{2}(t^{\prime})|\sum_{m+m^{\prime}>0}(2\beta
ht^{\prime})^{m+m^{\prime}}\frac{1}{m!m^{\prime}!}dt^{\prime}dt$
$\displaystyle\leq\int_{-\infty}^{\infty}n_{1}(t)\int_{-\infty}^{\infty}|n_{2}(t^{\prime})|\left(e^{4\beta
ht^{\prime}}-1\right)dt^{\prime}dt$
$\displaystyle\leq\frac{2}{\pi}\left(e^{\beta^{2}h^{2}}-1\right)\times\int_{-\infty}^{\infty}n_{1}(t)dt=\sqrt{\pi}\left(e^{\beta^{2}h^{2}}-1\right),$

where we bounded
$\left\|K_{m}(A_{0}^{a,\alpha\,})\right\|_{\infty}\leq(2h)^{m}$ for the $m$-th
nested commutator $K_{m}$ with $1$-local terms. For the last part, we compare
$\widetilde{\Psi}_{\beta,0}^{a,\alpha}(X)$ with the depolarizing generator
$\widetilde{\Psi}_{0,0}^{a,\alpha}(X)=\frac{1}{2\sqrt{2}e^{1/4}}A^{a,\alpha}XA^{a,\alpha\,\dagger}$
in a similar way as Eq. (19), so that

$\displaystyle\left\|\widetilde{\Psi}_{\beta,0}^{a,\alpha}(X)-\widetilde{\Psi}_{0,0}^{a,\alpha}(X)\right\|_{2}$
$\displaystyle\leq\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}h_{-}(t_{-})h_{+}(t_{+})\left\|A_{0}^{a,\alpha}(\beta(t_{+}-t_{-}))XA_{0}^{a,\alpha}(\beta(-t_{+}-t_{-}))^{\dagger}-A_{0}^{a,\alpha}XA_{0}^{a,\alpha\,\dagger}\right\|_{2}dt_{+}dt_{-}$
$\displaystyle\leq\left\|X\right\|_{2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}h_{-}(t_{-})h_{+}(t_{+})\times$
$\displaystyle\quad\quad\quad\left(\left\|A_{0}^{a,\alpha}(\beta(t_{+}-t_{-}))-A_{0}^{a,\alpha}\right\|_{\infty}+\left\|A_{0}^{a,\alpha\,\dagger}(\beta(-t_{+}-t_{-}))-A_{0}^{a,\alpha\,\dagger}\right\|_{\infty}\right)dt_{+}dt_{-}$
$\displaystyle\leq\left\|X\right\|_{2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}2h_{-}(t_{-})h_{+}(t_{+})\left(e^{2\beta
h(t_{+}-t_{-})}+e^{2\beta h(-t_{+}-t_{-})}-2\right)dt_{+}dt_{-}$
$\displaystyle=\frac{e^{\frac{3}{4}\beta^{2}h^{2}}-1}{\sqrt{2}e^{1/4}}\,\|X\|_{2}\,.$

This shows that $\sum_{a,\alpha}\widetilde{\mathcal{L}}_{\beta,0}^{a,\alpha}$
converges to the generator of the local depolarizing semigroup with a
multiplicative constant $\lambda$ as $\beta\to 0$. Clearly the latter is
gapped with $\gamma=-\lambda=\frac{1}{\sqrt{2}e^{1/4}}$, frustration-free,
satisfies local topological order with sharp cutoff and is also locally
gapped. Therefore,
$\widetilde{H}_{0}:=\widetilde{\mathcal{L}}_{0,0}=\sum_{a,\alpha}\widetilde{\mathcal{L}}_{0,0}^{a,\alpha}$
satisfies the conditions of the unperturbed Hamiltonian in Theorem 6. Together
with the exponential decay of Eq. (18), this shows that for $\beta<4/J$, the
perturbed generator $\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}$ obeys all of
the assumptions of Theorem 6. In this case, both global and local gaps are
$\alpha$, the decay $f(r)$ is given by Eq. (18), and the existence of a finite
Lieb-Robinson velocity follows from the quasi-locality of
$\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}$ [24]. This implies that there
exists $\beta^{*}=\mathcal{O}(1)$ depending on the parameter
$\widetilde{J}_{0}$ introduced in Theorem 6 (which itself depends on the
lattice dimension and the Lieb-Robinson velocity of
$\widetilde{\mathcal{L}}^{a,\alpha}_{\beta}$), such that for
$\beta<\beta^{*}\leq 4/J$ the spectral gap of
$\widetilde{\mathcal{L}}_{\beta}$, and therefore that of
$\mathcal{L}_{\beta}$, is lower bounded by $\frac{1}{2\sqrt{2}e^{1/4}}$.

###  B.3 Adiabatic algorithm

We now make use of the gap proved in Theorem 1 to devise an adiabatic
evolution on the $n$-fold tensor product of Hilbert spaces
$\widetilde{\mathcal{H}}_{v}\simeq\mathbb{C}^{2}\otimes\mathbb{C}^{2}$
interpolating between $\mathscr{H}(0):=\widetilde{\mathcal{L}}_{0}$ and
$\mathscr{(}1):=\widetilde{\mathcal{L}}_{\beta}$. For this, we resort to the
following standard estimate on the performance of the adiabatic evolution:

######  Lemma 7.

[38] Let $\mathscr{H}(s)$ with $s\in[0,1]$ be an adiabatic path with minimum
gap $\gamma$. The minimum adiabatic time $T_{\operatorname{ad}}$ required to
achieve an error $\epsilon$ in the final state is bounded as

$T_{\operatorname{ad}}\leq\frac{10}{\epsilon^{2}}\max\left\\{\frac{\left\|\frac{d\mathscr{H}}{ds}\right\|_{\infty}^{3}}{\gamma^{4}},\frac{\left\|\frac{d\mathscr{H}}{ds}\right\|_{\infty}\left\|\frac{d^{2}\mathscr{H}}{ds^{2}}\right\|_{\infty}}{\gamma^{3}}\right\\}$
(20)

In our case, the adiabatic path is
$\mathscr{H}(s)=\widetilde{\mathcal{L}}_{\beta s}$, with $\beta_{s}:=s\beta$,
starting at $\beta_{0}=0$ and ending at the desired temperature
$\beta=\beta_{1}$ at $s=1$. Thus, we need to estimate the derivatives of the
Hamiltonian path $\\{\mathscr{H}_{s}\\}_{s\in[0,1]}$ w.r.t. parameter $s$,
done in the subsections below. In particular, in (26) and (27) we show that

$\displaystyle\max_{s\in[0,1]}\left\|\frac{d}{ds}{\mathscr{H}}(s)\right\|_{\infty}=\max_{\beta^{\prime}\in[0,\beta]}\beta^{\prime}\left\|\frac{d}{d\beta^{\prime}}\widetilde{\mathcal{L}}_{\beta^{\prime}}\right\|_{2\to
2}\leq
61\,\beta\,n\,\max_{a,\alpha}\left\|\left[H,A^{a,\alpha}\right]\right\|_{\infty}$
(21)
$\displaystyle\max_{s\in[0,1]}\left\|\frac{d^{2}}{ds^{2}}{\mathscr{H}}(s)\right\|_{\infty}=\max_{\beta^{\prime}\in[0,\beta]}(\beta^{\prime})^{2}\left\|\frac{d^{2}}{d\beta^{\prime
2}}\widetilde{\mathcal{L}}_{\beta^{\prime}}\right\|_{2\to 2}\leq
C\,\beta^{2}\,n\left(\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}+\left\|[H,A^{a,\alpha}]\right\|_{\infty}^{2}\right)\,,$
(22)

for some constant $C>0$. Morever, notice that for a $(k,l)$-local Hamiltonian,
with the notation of Sec. A,

$\displaystyle\max_{a,\alpha}\left\|\left[H,A^{a,\alpha}\right]\right\|_{\infty}$
$\displaystyle\leq 2hl$
$\displaystyle\max_{a,\alpha}\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}$
$\displaystyle\leq 4h^{2}l^{2}k,$

which are $\mathcal{O}$(1) constants. Combining these bounds with Theorem 1,
we obtain the bound on the runtime of the adiabatic evolution preparing the
purified Gibbs state corresponding to $\sigma_{\beta}$:

######  Theorem 3 (Run-time of the adiabatic evolution).

In the notations of the previous paragraphs, and given $\beta^{*}$ as
introduced in Theorem 1, for any $\beta\in[0,\beta^{*}]$, the minimum time
required to adiabatically prepare a state $\epsilon$-close to the purified
Gibbs state at inverse temperature $\beta$ is bounded as

$\displaystyle T_{\operatorname{ad}}\leq\frac{C^{\prime}kl^{2}(\beta
hn)^{3}}{\epsilon^{2}}\,,$ (23)

for some universal constant $C^{\prime}>0$.

####  B.3.1 First derivative

It remains to bound the first and second derivatives of
$\widetilde{\mathcal{L}}_{\beta}$ w.r.t. the parameter $\beta$. This
Hamiltonian can be written in terms of two contributions
$\widetilde{\mathcal{L}}_{\beta}=\mathscr{H}_{\text{cd}}+\mathscr{H}_{\text{tr}}$.
First, it is possible to bound at once the coherent term corresponding to
$\widetilde{\mathcal{C}}_{\beta}$, which is
$\Delta_{\beta}^{-\frac{1}{4}}(B)\otimes{I}+{I}\otimes\Delta_{\beta}^{-\frac{1}{4}}(B)^{*}$
in the purified picture, and the decay term, corresponding to
$\widetilde{\mathcal{D}}_{\beta}$. Similarly to Eq. (13), in Corollary III.5
from [13], it is shown that their total contribution to the effective
Hamiltonian is of the form
$\mathscr{H}_{\text{cd}}=\frac{1}{2}(N\otimes{I}+{I}\otimes N^{*})$ in the
purified picture, where $N^{*}$ stands for the complex conjugate of $N$ when
decomposed in the computational basis, and where we recall that

$\displaystyle
N=\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}n_{1}(t)e^{-i\beta
Ht}\int_{-\infty}^{\infty}n_{2}(t^{\prime})e^{i\beta
Ht^{\prime}}A^{a,\alpha}e^{-2i\beta Ht^{\prime}}A^{a,\alpha}e^{i\beta
Ht^{\prime}}dt^{\prime}e^{i\beta Ht}dt,$

with functions $n_{1}$ and $n_{2}$ defined in Eq. (15) and (16) respectively.
We see that

$\displaystyle\frac{d}{d\beta}N=\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}n_{1}(t)\int_{-\infty}^{\infty}n_{2}(t^{\prime})\mathcal{M}(t,t^{\prime})dt^{\prime}dt.$

where

$\displaystyle\mathcal{M}(t,t^{\prime})$ $\displaystyle=-it\left[H,e^{i\beta
H(t^{\prime}-t)}A^{a,\alpha}e^{-2i\beta Ht^{\prime}}A^{a,\alpha}e^{i\beta
H(t^{\prime}+t)}\right]$ (24) $\displaystyle+it^{\prime}e^{i\beta
H(t^{\prime}-t)}[H,A^{a,\alpha}]e^{-2i\beta Ht^{\prime}}A^{a,\alpha}e^{i\beta
H(t^{\prime}+t)}+it^{\prime}e^{i\beta H(t^{\prime}-t)}A^{a,\alpha}e^{-2i\beta
Ht^{\prime}}[H,A^{a,\alpha}]e^{i\beta H(t^{\prime}+t)}\,.$

After some algebra, we find that

$\displaystyle\left\|\mathcal{M}(t,t^{\prime})\right\|_{\infty}\leq
2\left\|[H,A^{a,\alpha}]\right\|_{\infty}\left(|t|+|t^{\prime}|\right),$

so that we can bound

$\displaystyle\left\|\frac{d}{d\beta}N\right\|_{\infty}$ $\displaystyle\leq
2\max_{a,\alpha}\left\|[H,A^{a,\alpha}]\right\|_{\infty}\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\left\|A^{a,\alpha}\right\|_{\infty}^{2}\int_{-\infty}^{\infty}|n_{1}(t)|\int_{-\infty}^{\infty}|n_{2}(t^{\prime})|\left(|t|+|t^{\prime}|\right)dt^{\prime}dt$
$\displaystyle=\frac{6e^{\frac{\pi^{2}}{2}}}{\sqrt{2}\pi^{2}}\,n\,\max_{a,\alpha}\left\|[H,A^{a,\alpha}]\right\|_{\infty}.$

Here, we used that $\|A^{a,\alpha}\|_{\infty}=1,$ and evaluated the integrals
explicitly.

Similarly, we upper bound the derivative of the transition term of the
Hamiltonian, corresponding to the contribution of $\widetilde{\Psi}_{\beta}$,
which we can write as

$\displaystyle\mathscr{H}_{\text{tr}}=\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}h_{-}(t_{-})h_{+}(t_{+})A^{a,\alpha}\left(\beta(t_{+}-t_{-})\right)\otimes
A^{a,\alpha\,*}\left(\beta(t_{+}-t_{-})\right)dt_{-}dt_{+},$

where we have re-scaled $t,t^{\prime}$ so that $h_{-}(t_{-}),h_{+}(t_{+})$ are
independent of $\beta$. We get

$\displaystyle\left\|\frac{d}{d\beta}\mathscr{H}_{\text{tr}}\right\|_{\infty}$
$\displaystyle\leq\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}2h_{-}(t_{-})h_{+}(t_{+})\left\|A^{a,\alpha}\right\|_{\infty}\left\|\frac{d}{d\beta}\Delta_{\beta}^{i(t_{+}-t_{-})}(A^{a,\alpha})\right\|_{\infty}$
$\displaystyle\leq
6\,n\,\max_{a,\alpha}\left\|\left[H,A^{a,\alpha}\right]\right\|_{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|t_{+}-t_{-}|h_{-}(t_{-})h_{+}(t_{+})dt_{-}dt_{+},$
$\displaystyle\leq\frac{3\left(\sqrt{2}+2\right)}{4\sqrt[4]{e}\sqrt{\pi}}\,n\,\max_{a,\alpha}\left\|\left[H,A^{a,\alpha}\right]\right\|_{\infty},$

where we used once again that $\left\|A^{a,\alpha}\right\|_{\infty}=1$ and

$\displaystyle\left\|\frac{d}{d\beta}\Delta_{\beta}^{i(t_{+}-t_{-})}A^{a,\alpha}\right\|_{\infty}$
$\displaystyle\leq|t_{+}-t_{-}|\left\|\left[H,A^{a,\alpha}\right]\right\|_{\infty}.$
(25)

Putting everything together, we have that

$\displaystyle\left\|\frac{d}{d\beta}\widetilde{\mathcal{L}}_{\beta}\right\|_{2\to
2}\leq
61\,n\,\max_{a,\alpha}\left\|\left[H,A^{a,\alpha}\right]\right\|_{\infty}.$
(26)

####  B.3.2 Second derivative

For the coherent and decaying parts, we have

$\displaystyle\frac{d^{2}}{d\beta^{2}}N=\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}n_{1}(t)\int_{-\infty}^{\infty}n_{2}(t^{\prime})\frac{d}{d\beta}\mathcal{M}(t,t^{\prime})dt^{\prime}dt,$

we can differentiate $\frac{d}{d\beta}\mathcal{M}(t,t^{\prime})$ from Eq. (24)
to obtain a somewhat lengthy expression of nested commutators. From this we
obtain

$\displaystyle\left\|\frac{d}{d\beta}\mathcal{M}(t,t^{\prime})\right\|_{\infty}\leq\left(\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}+\left\|[H,A^{a,\alpha}]\right\|_{\infty}^{2}\right)\left(c_{1}t^{2}+c_{2}t^{\prime
2}+c_{3}tt^{\prime})\right),$

for $c_{i}$ some $\mathcal{O}(1)$ constants. Proceeding as above, this means
that, for some constant $C_{\text{cd}}$,

$\displaystyle\left\|\frac{d^{2}}{d\beta^{2}}N\right\|_{\infty}\leq\frac{C_{\text{cd}}}{2}\left(\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}+\left\|[H,A^{a,\alpha}]\right\|_{\infty}^{2}\right).$

For the transition part, using the chain rule and the triangle inequality, we
have that

$\displaystyle\left\|\frac{d^{2}}{d\beta^{2}}\mathscr{H}_{\operatorname{tr}}\right\|_{\infty}$
$\displaystyle\leq\sum_{a\in\Lambda}\sum_{\alpha\in[3]}\,\int_{-\infty}^{\infty}2h_{-}(t_{-})h_{+}(t_{+})\left(\left\|\frac{d^{2}}{d\beta^{2}}\Delta_{\beta}^{i(t_{+}-t_{-})}A^{a,\alpha}\right\|_{\infty}+\left\|\frac{d}{d\beta}\Delta_{\beta}^{i(t_{+}-t_{-})}A^{a,\alpha}\right\|_{\infty}^{2}\right).$

The second term is bounded by (25). The first follows from the fact that

$\displaystyle\frac{d^{2}}{d\beta^{2}}\Delta_{\beta}^{i(t_{+}-t_{-})}A^{a,\alpha}=-(t_{+}-t_{-})^{2}[H,[H,\Delta_{\beta}^{i(t_{+}-t_{-})}A^{a,\alpha}]],$

so that
$\left\|\frac{d^{2}}{d\beta^{2}}\Delta_{\beta}^{i(t_{+}-t_{-})-\frac{1}{4}}A^{a,\alpha}\right\|_{\infty}\leq(t_{+}-t_{-})^{2}\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}$.
The contribution is bounded by

$\displaystyle\left\|\frac{d^{2}}{d\beta^{2}}\mathcal{H}_{\operatorname{tr}}\right\|_{\infty}$
$\displaystyle\leq
6n\,\max_{a,\alpha}\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}\int_{-\infty}^{\infty}h_{-}(t_{-})h_{+}(t_{+})(t_{+}-t_{-})^{2}dt_{-}dt_{+}$
$\displaystyle\leq\frac{3\left(3\pi\sqrt{2}+8\right)}{16\sqrt[4]{e}\pi}n\,\max_{a,\alpha}\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}.$

The numerical constant above is close to $1$, so that finally we bound the
second derivative as

$\displaystyle\left\|\frac{d^{2}}{d\beta^{2}}\widetilde{\mathcal{L}}_{\beta}\right\|_{2\to
2}\leq
n\,\left(\left\|[H,[H,A^{a,\alpha}]]\right\|_{\infty}+\left\|[H,A^{a,\alpha}]\right\|_{\infty}^{2}\right)\left(C_{\operatorname{cd}}+1\right).$
(27)

##  Appendix C Zero temperature quantum Gibbs sampling and the complexity
class $\mathsf{GibbsQP}$

In this appendix, we prove the equivalence of the classes $\mathsf{BQP}$ and
$\mathsf{GibbsQP}$, that is Theorem 5. In short, our proof relies on a
perturbative analysis of the generator of the Gibbs sampling evolution studied
in the previous sections around $T:=\frac{1}{\beta}=0$. The Hamiltonian we
consider is the original Kitaev circuit Hamiltonian already used in the proof
that $\mathsf{BQP}=\mathsf{AdiabQP}$ [45]. Interestingly, a similar argument
was already considered in the classical literature (see e.g. [21]). However,
we need a different filter $\gamma(\omega)$ than the Gaussian one used in the
previous section, since it leads to a direct sum of generators fixing the
projections onto any eigenspaces of the corresponding Hamiltonian in the low
temperature limit. In contrast, Metropolis and Glauber filters favour
transitions to lower energies and hence stabilise the ground states of $H$ in
the low temperature limit.

Here, we make use of the Metropolis-type filter already introduced in [13]:
given inverse temperature $\beta>0$ and tunable width $\sigma_{E}>0$ to be
chosen later,

$\displaystyle\gamma^{M}(\omega):=e^{-\beta\max\big{(}\omega+\frac{\beta\sigma_{E}^{2}}{2},0\big{)}}\,.$

Using the results of [13], we directly find

$\displaystyle\mathcal{L}^{\dagger}_{\beta,\sigma_{E}}(X)=i\sum_{\nu}\big{[}B_{\nu},X\big{]}+\sum_{a}\sum_{\nu_{1},\nu_{2}}\alpha_{\nu_{1},\nu_{2}}\left((A^{a}_{\nu_{1}})^{\dagger}XA_{\nu_{2}}^{a}-\frac{1}{2}\big{\\{}(A_{\nu_{1}}^{a})^{\dagger}A_{\nu_{2}}^{a},X\big{\\}}\right)\equiv
i[B,X]+\mathcal{D}_{\beta,\sigma_{E}}^{\dagger}(X)\,,$

where the jumps $\\{A^{a}\\}_{a\in[m]}$ with $\|A^{a}\|_{\infty}=1$ will be
specified later, $A^{a}_{\nu}:=\sum_{E_{i}-E_{j}=\nu}P_{i}A^{a}P_{j}$,

$\displaystyle B_{\nu}:=\sum_{a}\sum_{\begin{subarray}{c}\nu_{1},\nu_{2}\\\
\nu_{1}-\nu_{2}=\nu\end{subarray}}\,\frac{\operatorname{tanh}(-\beta(\nu_{1}-\nu_{2})/4)}{2i}\,\alpha_{\nu_{1},\nu_{2}}(A^{a}_{\nu_{2}})^{\dagger}A_{\nu_{1}}^{a}\equiv\sum_{a}\sum_{\begin{subarray}{c}\nu_{1},\nu_{2}\\\
\nu_{1}-\nu_{2}=\nu\end{subarray}}\,b_{\nu_{1},\nu_{2}}(A^{a}_{\nu_{2}})^{\dagger}A_{\nu_{1}}^{a}\,,$

and

$\displaystyle\alpha_{\nu_{1},\nu_{2}}=\frac{1}{4}e^{-\beta\frac{(\nu_{1}+\nu_{2})}{2}-\frac{(\nu_{1}-\nu_{2})^{2}}{8\sigma_{E}^{2}}}\,\left(1+\operatorname{Erf}\Big{(}\frac{\nu_{1}+\nu_{2}-\beta\sigma_{E}^{2}}{2\sqrt{2}\sigma_{E}}\Big{)}+e^{\frac{1}{2}\beta(\nu_{1}+\nu_{2})}\operatorname{Erfc}\Big{(}\frac{\nu_{1}+\nu_{2}+\beta\sigma_{E}^{2}}{2\sqrt{2}\sigma_{E}}\Big{)}\right)\,.$
(28)

Thus, this generator is similar to the one in Sec. A.2, but with different
coefficients for the jump operators, as per the Metropolis or Glauber filters.

###  C.1 Perturbation bounds

Choosing $\sigma_{E}=\frac{1}{\beta}$ and taking the limit $\beta\to\infty$,
one can see that $\alpha_{\nu_{1},\nu_{2}}\to\frac{1}{2}$ whenever
$\nu_{1}=\nu_{2}<0$, and
$\alpha_{0,0}=\frac{1}{2}\operatorname{Erfc}(1/2\sqrt{2})$ for all $\beta$.
Therefore,

$\displaystyle\mathcal{L}^{\dagger}_{\infty,0}(X)=\frac{1}{2}\sum_{a,\nu<0}\,\left((A^{a}_{\nu})^{\dagger}XA^{a}_{\nu}-\frac{1}{2}\big{\\{}(A^{a}_{\nu})^{\dagger}A^{a}_{\nu},X\big{\\}}\right)+\frac{1}{2}\operatorname{Erfc}(1/2\sqrt{2})\sum_{a}\left((A^{a}_{0})^{\dagger}XA^{a}_{0}-\frac{1}{2}\big{\\{}(A^{a}_{0})^{\dagger}A^{a}_{0},X\big{\\}}\right)\,.$
(29)

Our first technical result is a control of the rate at which
$\mathcal{L}_{\beta}^{\dagger}\equiv{\mathcal{L}}^{\dagger}_{\beta,\frac{1}{\beta}}$
converges to
${\mathcal{L}}_{\infty}^{\dagger}\equiv\mathcal{L}^{\dagger}_{\infty,0}$ as
$\beta\to\infty$. In what follows, we will also denote the generator
$\mathcal{L}^{\dagger}_{\beta}$ by $\mathcal{L}^{\dagger}_{\beta,H}$ in order
to keep track of the Hamiltonian it depends on.

######  Proposition 8.

For any Hamiltonian $H:=\sum_{i}E_{i}P_{i}\in\mathcal{O}_{V}$, $|V|=n$, with
$M$ distinct eigenvalues and minimum gap between any two energies
$\Delta_{E}:=\min_{j,k}\\{|E_{j}-E_{k}|\\}$ and between any two Bohr
frequencies $\Delta_{\nu}:=\min_{j,k}\\{|\nu_{j}-\nu_{k}|\\}$, and
$\beta>(2\Delta_{E})^{-1}$,

$\displaystyle\|{\mathcal{L}}^{\dagger}_{\infty,H}-{\mathcal{L}}^{\dagger}_{\beta,H}\|_{\infty\to\infty}\leq
3mM^{4}\max\left\\{\,e^{-\frac{-\beta^{2}\Delta_{\nu}^{2}}{8}},\,\frac{4\beta\Delta_{E}}{2\beta\Delta_{E}-1}e^{-\frac{1}{4}\beta\Delta_{E}}\right\\}\,.$
(30)

###### Proof.

From the discussion at the beginning of the section, it suffices to control
the rate of convergence of the coefficients $\alpha_{\nu_{1},\nu_{2}}$ and
$b_{\nu_{1},\nu_{2}}$ as $\beta\to\infty$ for all Bohr frequencies
$\nu_{1},\nu_{2}$. First of all, we control the convergence of the
coefficients $\alpha_{\nu_{1},\nu_{2}}$: whenever $\nu_{1}=\nu_{2}>0$, we get

$\displaystyle\alpha_{\nu_{1},\nu_{1}}$
$\displaystyle=\frac{1}{4}e^{-\beta{\nu_{1}}}\,\left(1+\operatorname{Erf}\Big{(}\frac{2\beta\nu_{1}-1}{2\sqrt{2}}\Big{)}+e^{\beta\nu_{1}}\operatorname{Erfc}\Big{(}\frac{2\beta\nu_{1}+1}{2\sqrt{2}}\Big{)}\right)$
$\displaystyle\leq\frac{1}{4}\,\left(2e^{-\beta\nu_{1}}+2e^{-\frac{1}{8}(1+2\beta\nu_{1})^{2}}\right)$
$\displaystyle\leq 2e^{-\frac{1}{2}\beta\Delta_{E}}\,,$

where we used Chernoff bound in the first inequality: for all $x>0$

$\displaystyle\operatorname{Erfc}(x)\leq\frac{e^{-x^{2}}}{x\sqrt{\pi}}\,.$

In contrast, whenever $\nu_{1}=\nu_{2}<0$, we have that, whenever
$2\beta\Delta_{E}>1$

$\displaystyle\left|\alpha_{\nu_{1},\nu_{1}}-\frac{1}{2}\right|=\left|\frac{1}{4}\,\left(-1-\operatorname{Erf}\Big{(}\frac{1}{2\sqrt{2}}(1+2\beta\nu_{1})\Big{)}+e^{-\beta\nu_{1}}\operatorname{Erfc}\Big{(}\frac{1}{2\sqrt{2}}(1-2\beta\nu_{1})\Big{)}\right)\right|\leq
e^{-\frac{1}{4}\beta\Delta_{E}}\,\frac{2\beta\Delta_{E}}{2\beta\Delta_{E}-1}\,.$

Next, whenever $\nu_{1}\neq\nu_{2}$, we have that

$\displaystyle\alpha_{\nu_{1},\nu_{2}}\leq
e^{-\frac{\beta^{2}(\nu_{1}-\nu_{2})^{2}}{8}}\leq
e^{-\frac{\beta^{2}\Delta_{\nu}^{2}}{8}}\,.$

Finally, we need to control $b_{\nu_{1},\nu_{2}}$. Since whenever
$\nu_{1}=\nu_{2}$, $b_{\nu_{1},\nu_{2}}=0$, it suffices to control the limit
of $b_{\nu_{1},\nu_{2}}$ to $0$ for $\nu_{1}\neq\nu_{2}$. We get

$\displaystyle\Big{|}b_{\nu_{1},\nu_{2}}\Big{|}\leq\frac{1}{2}\,\Big{|}\alpha_{\nu_{1},\nu_{2}}\Big{|}\leq\frac{1}{2}\,e^{-\frac{\beta^{2}\Delta_{\nu}^{2}}{8}}\,.$

Using all the bounds derived above, we finally have that

$\displaystyle\|{\mathcal{L}}^{\dagger}_{\beta,H}-{\mathcal{L}}_{\infty,H}^{\dagger}\|_{\infty\to\infty}\leq
3mM^{4}\max\left\\{\,e^{-\frac{-\beta^{2}\Delta_{\nu}^{2}}{8}},\,\frac{4\beta\Delta_{E}}{2\beta\Delta_{E}-1}e^{-\frac{1}{4}\beta\Delta_{E}}\right\\}$

as claimed. ∎

In what follows, we perform a perturbative analysis of the generators at low
temperature. For this, it will be more convenient to work in the time
representation. In [13], it was proven that $B_{\beta,H}:=\sum_{\nu}\,B_{\nu}$
is such that $\|B_{\beta,H}-B^{\eta}_{\beta,H}\|_{\infty}\leq
m\eta\beta\|H\|_{\infty}$, with $0<\eta<1$,

$\displaystyle
B^{\eta}_{\beta,H}:=\sum_{a}\,\int_{-\infty}^{\infty}\,b_{1}(t)\,e^{-i\beta
Ht}\left(\int_{-\infty}^{\infty}\,b^{\eta}_{2}(t^{\prime})\,A_{H}^{a\dagger}(\beta
t^{\prime})A_{H}^{a}(-\beta
t^{\prime})dt^{\prime}+\frac{1}{16\sqrt{2}\,\pi}(A^{a})^{\dagger}A^{a}\right)e^{i\beta
Ht}dt$ (31)

where $O_{H}(t):=e^{itH}Oe^{-itH}$, and with functions

$\displaystyle
b_{1}(t):=2\sqrt{\pi}\,e^{\frac{1}{8}}\,\left(\frac{1}{\operatorname{cosh}(2\pi
t)}\ast\sin(-t)\operatorname{exp}(-2t^{2})\right)\qquad\text{ with
}\quad\|b_{1}\|_{L_{1}}<1$ $\displaystyle
b^{\eta}_{2}(t)=\frac{1}{4\sqrt{2}\,\pi}\,\frac{\operatorname{exp}(-2t^{2}-it)+\mathsf{1}(|t|\leq\eta)i(2t+i)}{t(2t+i)}\qquad\text{
with
}\quad\|b^{\eta}_{2}\|_{L_{1}}<\frac{1}{5}+\frac{1}{2\sqrt{2}\pi}\ln(1/\eta)\,.$

Moreover, it is easy to verify that

$\displaystyle\|t\mapsto|t|\,b_{1}(t)\|_{L_{1}}\leq
2\sqrt{\pi}e^{\frac{1}{8}}\,\int_{-\infty}^{\infty}\,|t|\,e^{-2t^{2}}\,dt=\sqrt{\pi}\,e^{\frac{1}{8}}\,,$
$\displaystyle\|t\mapsto|t|\,b^{\eta}_{2}(t)\|_{L_{1}}\leq 1+\eta\,.$ (32)

Using similar ideas, the dissipative part can be re-expressed as follows:

######  Lemma 9.

In the notations of the previous paragraph we denote, for any $0<\eta<1$,

$\displaystyle\mathcal{D}^{\eta,\dagger}_{\beta,H}(X):=\sum_{a}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{ie^{-\frac{x^{2}}{2}}\upsilon^{\eta}(y)}{4\pi\sqrt{\pi/2}}\,F^{a}(X,\beta(x+y)/2,\beta(x-y)/2)dxdy+\sum_{a}\int_{-\infty}^{\infty}\frac{e^{-\frac{2t^{2}}{\beta^{2}}}}{\beta
2\sqrt{\pi/2}}\,F^{a}(X,t,t)dt$ (33)

where

$\displaystyle\upsilon^{\eta}(y):=\frac{e^{-\frac{y^{2}}{2}+\frac{iy}{2}}+\mathsf{1}(|y|\leq\eta)i(i+y)}{y(i+y)}\,\qquad\text{
with }\quad\|\upsilon^{\eta}\|_{L_{1}}\leq C+2\ln(1/\eta)$

for some universal constant $C$, and

$\displaystyle
F^{a}(X,t,t^{\prime}):={A_{H}^{a\dagger}(t)XA_{H}^{a}(t^{\prime})-\frac{1}{2}\Big{\\{}A_{H}^{a\dagger}(t)A_{H}^{a}(t^{\prime}),X\Big{\\}}}\,.$

Then the map
$\mathcal{D}_{\beta,H}:=\sum_{a}\int\gamma^{M}(\omega)\mathcal{D}^{a}_{\omega}\,d\omega$,
with
$\mathcal{D}_{\omega}^{a}(X):=A^{a}(\omega)XA^{a}(\omega)^{\dagger}-\frac{1}{2}\big{\\{}A^{a}(\omega)^{\dagger}A^{a}(\omega),X\big{\\}}$,
satisfies

$\displaystyle\left\|\mathcal{D}^{\eta,\dagger}_{\beta,H}-\mathcal{D}^{\dagger}_{\beta,H}\right\|_{\infty\to\infty}\leq
m\beta\eta C^{\prime}\|H\|_{\infty}$

for some universal constant $C^{\prime}$.

###### Proof.

We start by re-expressing the map $\mathcal{D}_{\beta,H}^{\dagger}$ as
follows:

$\displaystyle\mathcal{D}_{\beta,H}^{\dagger}(X)$
$\displaystyle=\sum_{a}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\,\frac{\gamma^{M}(\omega)e^{-\frac{t^{2}+(t^{\prime})^{2}}{\beta^{2}}}}{\beta
2\pi\sqrt{\pi/2}}\,\left(A_{H}^{a\dagger}(t)XA_{H}^{a}(t^{\prime})-\frac{1}{2}\Big{\\{}A_{H}^{a\dagger}(t)A_{H}^{a}(t^{\prime}),X\Big{\\}}\right)\,e^{-i\omega(t-t^{\prime})}\,d\omega\,dt\,dt^{\prime}$
$\displaystyle=\sum_{a}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\,\frac{\hat{\gamma}^{M}(t-t^{\prime})e^{-\frac{t^{2}+(t^{\prime})^{2}}{\beta^{2}}}}{\beta
2\pi\sqrt{\pi/2}}\,{\left(A_{H}^{a\dagger}(t)XA_{H}^{a}(t^{\prime})-\frac{1}{2}\Big{\\{}A_{H}^{a\dagger}(t)A_{H}^{a}(t^{\prime}),X\Big{\\}}\right)}\,dt\,dt^{\prime}$
$\displaystyle=\sum_{a}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{ie^{-\frac{t^{2}+(t^{\prime})^{2}}{\beta^{2}}+\frac{i(t-t^{\prime})}{2\beta}}}{2\pi\sqrt{\pi/2}(\beta+i(t-t^{\prime}))(t-t^{\prime})}\,F^{a}(X,t,t^{\prime})dtdt^{\prime}+\sum_{a}\int_{-\infty}^{\infty}\frac{e^{-\frac{2t^{2}}{\beta^{2}}}}{\beta
2\sqrt{\pi/2}}\,F^{a}(X,t,t)dt$
$\displaystyle=\sum_{a}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{ie^{-\frac{u^{2}+v^{2}}{2\beta^{2}}+\frac{iv}{2\beta}}}{4\pi\beta^{2}\sqrt{\pi/2}(1+i\frac{v}{\beta})\frac{v}{\beta}}\,F^{a}(X,(u+v)/2,(u-v)/2)dudv+\sum_{a}\int_{-\infty}^{\infty}\frac{e^{-\frac{2t^{2}}{\beta^{2}}}}{\beta
2\sqrt{\pi/2}}\,F^{a}(X,t,t)dt$
$\displaystyle=\sum_{a}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{ie^{-\frac{x^{2}+y^{2}}{2}+\frac{iy}{2}}}{4\pi\sqrt{\pi/2}(1+i{y})y}\,F^{a}(X,\beta(x+y)/2,\beta(x-y)/2)dxdy+\sum_{a}\int_{-\infty}^{\infty}\frac{e^{-\frac{2t^{2}}{\beta^{2}}}}{\beta
2\sqrt{\pi/2}}\,F^{a}(X,t,t)dt,$

where we used that $\hat{\gamma}^{M}(t)\equiv\int\gamma^{M}(\omega)e^{-i\omega
t}d\omega=\frac{i\beta e^{\frac{it}{2\beta}}}{(\beta+it)t}+\pi\delta(t)$ in
the sense of distributions, made the change of variables $u=t+t^{\prime}$,
$v=t-t^{\prime}$ in the penultimate line, and $x=u/\beta$, $y=v/\beta$ in the
last one. Next, given $0<\eta<1$ we consider the function

$\displaystyle\upsilon^{\eta}(y):=\frac{e^{-\frac{y^{2}}{2}+\frac{iy}{2}}+\mathsf{1}(|y|\leq\eta)i(i+y)}{y(i+y)}\,.$

It is easy to see that

$\displaystyle\upsilon^{\eta}(y)=\upsilon^{1}(y)-\frac{i\,\mathsf{1}(\eta<|y|\leq
1)}{y}\,.$

By the triangle inequality combined with Hölder’s inequality, we find that

$\displaystyle\|\upsilon^{\eta}\|_{L_{1}}$
$\displaystyle\leq\|y\mapsto(1+y^{2})^{-1}\|_{L_{2}}\,\|y\mapsto(1+y^{2})\upsilon^{1}(y)\|_{L_{2}}+\left\|y\mapsto\frac{\mathsf{1}(\eta<|y|\leq
1)}{y}\right\|_{L_{1}}$
$\displaystyle\leq\sqrt{\frac{\pi}{2}}\,C+2\ln(1/\eta)$

where in the last line we have used that

$\displaystyle
C^{2}\equiv\int_{-\infty}^{\infty}(1+y^{2})^{2}|\upsilon^{1}(y)|^{2}\,dy<\infty$

We are now ready to control the norm difference
$\left\|\mathcal{D}^{\eta,\dagger}_{\beta,H}-\mathcal{D}^{\dagger}_{\beta,H}\right\|_{\infty\to\infty}$.
By definition, we have that

$\displaystyle(\mathcal{D}^{\eta,\dagger}_{\beta,H}-\mathcal{D}^{\dagger}_{\beta,H})(X)=\sum_{a}\int_{-\infty}^{\infty}\frac{-e^{-\frac{x^{2}}{2}}}{4\pi\sqrt{\pi/2}}\left(\int_{-\eta}^{\eta}\frac{1}{y}\,F^{a}(X,\beta(x+y)/2,\beta(x-y)/2)\,dy\,\right)dx\,.\
$ (34)

Let us now focus on the integrand in $y$ above. We have

$\displaystyle F^{a}(X,\beta(x+y)/2,\beta(x-y)/2)$
$\displaystyle=e^{iH\frac{\beta}{2}(x+y)}A^{a\dagger}e^{-iH\frac{\beta}{2}(x+y)}Xe^{iH\frac{\beta}{2}(x-y)}A^{a}e^{-iH\frac{\beta}{2}(x-y)}$
$\displaystyle-\frac{1}{2}\Big{\\{}e^{iH\frac{\beta}{2}(x+y)}A^{a\dagger}e^{-iH\frac{\beta}{2}(x+y)}e^{iH\frac{\beta}{2}(x-y)}A^{a}e^{-iH\frac{\beta}{2}(x-y)},X\Big{\\}}$
$\displaystyle\equiv A_{1}+A_{2}\,.$

Let us first focus on the first term $A_{1}$: decomposing the $y$-dependent
pieces into a sum of cosinus and sinus, denoting $u:=\beta x/2$ and $v:=\beta
y/2$, and writing $c(x):=\cos(x)$ and $s(x)=\sin(x)$ we get

$\displaystyle
A_{1}=e^{iHu}\left(c(Hv)+is(Hv)\right)A^{a\dagger}e^{-iHu}\left(c(Hv)-is(Hv)\right)Xe^{iHu}\left(c(Hv)-is(Hv)\right)A^{a}e^{-iHu}\left(c(Hv)+is(Hv)\right)\,.$

From here, it is clear that this product can be developed into a sum of $16$
terms, the last of which is

$\displaystyle
C(v):=e^{iHu}c(Hv)A^{a\dagger}e^{-iHu}c(Hv)Xe^{iHu}c(Hv)A^{a}e^{-iHu}c(Hv)$

and all the others containing at least one $\sin(Hv)$ factor. By parity,
symmetry of the interval $[-\eta,\eta]$ around $0$ and since the integral over
$y$ above is even due to the $1/v$ dependence, we get that
$\int_{-\eta}^{\eta}C(v)/y=0$. In contrast, all the other $15$ terms are of
order $o(v)$ in the limit $v\to 0$, which implies that $y\mapsto S(x,y)$ is
integrable. A very similar argument for the term $A_{2}$ hence leads to the
conclusion that

$\displaystyle\left\|\mathcal{D}^{\eta,\dagger}_{\beta,H}-\mathcal{D}^{\dagger}_{\beta,H}\right\|_{\infty\to\infty}=m\beta
C^{\prime}\eta\|H\|_{\infty}\,.$ (35)

for some universal constant $C^{\prime}$. ∎

We end these preliminaries showing how the Lindblad generators are stable to
perturbations in the underlying Hamiltonian.

######  Proposition 10.

Let $H=H_{0}+V$ be a perturbed Hamiltonian over $n$ qubits, and let
${\mathcal{L}}^{\dagger}_{\beta,H_{0}}$, resp.
${\mathcal{L}}^{\dagger}_{\beta,H}$ be the above generators corresponding to
Hamiltonians $H_{0}$ and $H$. Then for any $0<\eta<1$

$\displaystyle\|{\mathcal{L}}^{\dagger}_{\beta,H}-{\mathcal{L}}^{\dagger}_{\beta,H_{0}}\|_{\infty\to\infty}\leq
m\beta\Big{(}C_{1}\eta(\|H_{0}\|_{\infty}+\|V\|_{\infty})+C_{2}\|V\|_{\infty}(1+\ln(1/\eta))\Big{)}\,,$

for some positive universal constants $C_{1}$ and $C_{2}$.

###### Proof.

From (31) and (C.1), it suffices to control the operator norms
$\|A^{a}_{H}(t)-A^{a}_{H_{0}}(t)\|_{\infty}$ for any $t\in\mathbb{R}$. From
obvious continuity bounds, we get that
$\|A^{a}_{H}(t)-A^{a}_{H_{0}}(t)\|_{\infty}\leq 2|t|\,\|V\|_{\infty}$.
Similarly,
$\|e^{-iHt}A^{a}_{H}(t^{\prime})e^{iHt}-e^{-iH_{0}t}A^{a}_{H_{0}}(t^{\prime})e^{iH_{0}t}\|_{\infty}\leq
2(|t|+|t^{\prime}|)\|V\|_{\infty}$. Therefore, together with the bounds (32),

$\displaystyle\|B_{\beta,H}-B_{\beta,H_{0}}\|_{\infty}$
$\displaystyle\leq\|B_{\beta,H}-B^{\eta}_{\beta,H}\|_{\infty}+\|B^{\eta}_{\beta,H}-B^{\eta}_{\beta,H_{0}}\|_{\infty}+\|B^{\eta}_{\beta,H_{0}}-B_{\beta,H_{0}}\|_{\infty}$
$\displaystyle\leq
m\eta\beta(\|H\|_{\infty}+\|H_{0}\|_{\infty})+4\beta\|V\|_{\infty}\sum_{a}\int_{-\infty}^{\infty}|b_{1}(t)|\left(\int_{-\infty}^{\infty}\,|b^{\eta}_{2}(t^{\prime})|\,(|t|+|t^{\prime}|)dt^{\prime}+|t|\right)dt$
$\displaystyle\leq
m\eta\beta(\|H\|_{\infty}+\|H_{0}\|_{\infty})+4\beta\,\|V\|_{\infty}m\,\left(1+\eta+\sqrt{\pi}e^{\frac{1}{8}}\Big{(}\frac{6}{5}+\frac{1}{2\sqrt{2}\pi}\ln(1/\eta)\Big{)}\right)\,.$

Similarly, we have that

$\displaystyle\big{\|}\mathcal{D}^{\dagger}_{\beta,H}-\mathcal{D}^{\dagger}_{\beta,H_{0}}\big{\|}_{\infty\to\infty}$
$\displaystyle\leq\big{\|}\mathcal{D}_{\beta,H}^{\dagger}-\mathcal{D}^{\eta\dagger}_{\beta,H}\big{\|}_{\infty\to\infty}+\big{\|}\mathcal{D}^{\eta\dagger}_{\beta,H}-\mathcal{D}^{\eta\dagger}_{\beta,H_{0}}\big{\|}_{\infty\to\infty}+\big{\|}\mathcal{D}^{\eta\dagger}_{\beta,H_{0}}-\mathcal{D}^{\dagger}_{\beta,H_{0}}\big{\|}_{\infty\to\infty}$
$\displaystyle\leq m\beta\eta
C^{\prime}(\|H\|_{\infty}+\|H_{0}\|_{\infty})+m\beta
C^{\prime\prime}\|V\|_{\infty}(1+\ln(1/\eta))$

The result follows by combining the above bounds. ∎

###  C.2 Controlling the time evolved Laplace transform

Here, we propose a simple proof of convergence of Gibbs samplers to states
that are close to the ground state of Hamiltonians encoding poly-sized quantum
circuits. The strategy that we follow is by a control of the Laplace transform
as a function of time: given $\theta,t\geq 0$ and an initial state $\rho$,

$\displaystyle\mathfrak{L}_{\rho,H,\beta}(\theta,t):=\operatorname{tr}\left(e^{t\mathcal{L}_{\beta,H}}(\rho)\,e^{\theta
H}\right)=\operatorname{tr}\left(e^{t\mathcal{L}_{\beta,H}}(\rho)P_{0}e^{\theta
H}\right)+\operatorname{tr}\left(e^{t\mathcal{L}_{\beta,H}}(\rho)P_{0}^{\perp}e^{\theta
H}\right)\,,$

where $P_{0}\equiv P_{0}(H)$, resp. $P^{\perp}_{0}\equiv P_{0}^{\perp}(H)$
denote the projection onto the ground space of $H$ and its orthogonal. First,
we consider an unperturbed Hamiltonian $H=H_{0}\geq 0$ with $P_{0}H_{0}=0$ to
be specified later and $\beta=\infty$. For simplicity, we also assume that the
initial state $\rho$ takes the form $\rho:=\sum_{i}\rho_{i}P_{i}(H_{0})$,
where $H_{0}=\sum_{i\neq 0}E_{i}(H_{0})P_{i}(H_{0})$ denotes the spectral
decomposition of $H_{0}$. The energies are ordered so that $0=E_{0}\leq
E_{1}\leq\dots\leq\|H_{0}\|_{\infty}$ (for instance, we can start the process
in the maximally mixed state, e.g. at $\beta=0$). Clearly,
$e^{t\mathcal{L}_{\infty,H_{0}}}(\rho)$ keeps this form for all $t\geq 0$ by
(29). Moreover, it is not difficult to see that
$\operatorname{tr}\left(e^{t\mathcal{L}_{\infty,H_{0}}}(\rho)P_{0}e^{\theta
H_{0}}\right)\leq 1$ for all $t\geq 0$. Hence,

$\displaystyle\mathfrak{L}_{\rho,H_{0},\infty}(\theta,t)\leq
1+\operatorname{tr}\left(e^{t\mathcal{L}_{\infty,H_{0}}}(\rho)P_{0}^{\perp}e^{\theta
H}\right)\,.$

Next, we show that
$\operatorname{tr}\left(e^{t\mathcal{L}_{\infty,H_{0}}}(\rho)P_{0}^{\perp}e^{\theta
H}\right)$ goes to $0$ quickly in $t$ for all $\theta>0$.

######  Lemma 11.

In the notations of the previous paragraph, assume that $A^{a}=A^{a\dagger}$
for all $a\in[m]$ and that the following Cheeger-type bound holds: there is a
constant $C>0$ such that, for all $i\neq 0$,

$\displaystyle\sum_{a}\operatorname{tr}\left(P_{i}(H_{0})A^{a}P_{i-1}(H_{0})A^{a}\right)\geq
C\operatorname{tr}\left(P_{i}(H_{0})\right)\,.$ (36)

Then for any initial state $\rho=\sum_{i}\rho_{i}P_{i}(H_{0})$, $t,\theta>0$,

$\displaystyle\mathfrak{L}_{\rho,H_{0},\infty}(\theta,t)\leq
1+e^{-\frac{C}{2}\theta\Delta_{E}e^{-\Delta_{E}\theta}t+\theta\|H_{0}\|_{\infty}}\,,$

where $\Delta_{E}\equiv\Delta_{E}(H_{0})$ denotes the minimal gap between any
two eigenvalues of $H_{0}$.

###### Proof.

We show an inequality of the form

$\displaystyle\operatorname{tr}\left(X\mathcal{L}_{\infty,H_{0}}^{\dagger}(P_{0}^{\perp}e^{\theta
H_{0}})\right)\leq-C^{\prime}\,\operatorname{tr}\left(XP_{0}^{\perp}e^{\theta
H_{0}}\right)$

for all $X=\sum_{i}x_{i}P_{i}\geq 0$ with $P_{i}\equiv P_{i}(H_{0})$, to
conclude that
$\operatorname{tr}\left(P_{0}^{\perp}e^{t\mathcal{L}_{\infty,H_{0}}}(\rho)e^{\theta
H_{0}}\right)\leq e^{-C^{\prime}t}\operatorname{tr}\left(\rho
P_{0}^{\perp}e^{\theta H_{0}}\right)\leq
e^{-C^{\prime}t+\theta\|H_{0}\|_{\infty}}$. So let’s try to do precisely that.
We have

$\displaystyle\operatorname{tr}\left(X\mathcal{L}_{\infty,H_{0}}^{\dagger}(P_{0}^{\perp}e^{\theta
H_{0}})\right)$
$\displaystyle\leq\frac{1}{4}\sum_{a,\nu<0}\,\big{(}e^{\theta\nu}-1\big{)}\operatorname{tr}\left(\big{\\{}X,(A^{a}_{\nu})^{\dagger}A^{a}_{\nu}\big{\\}}P_{0}^{\perp}e^{\theta
H_{0}}\right)$
$\displaystyle=\frac{1}{2}\sum_{a,\nu<0}\left(e^{\theta\nu}-1\right)\,\sum_{i\neq
0}\,x_{i}e^{\theta
E_{i}}\,\operatorname{tr}\left(P_{i}(A^{a}_{\nu})^{\dagger}A^{a}_{\nu}\right)$
$\displaystyle\leq\frac{1}{2}(e^{-\theta\Delta_{E}}-1)\sum_{i\neq
0}x_{i}e^{\theta
E_{i}}\,\sum_{a}\operatorname{tr}\left(P_{i}A^{a}P_{i-1}A^{a}\right)$
$\displaystyle\leq\frac{C}{2}(e^{-\theta\Delta_{E}}-1)\sum_{i\neq
0}x_{i}e^{\theta E_{i}}\,\operatorname{tr}\left(P_{i}\right)$
$\displaystyle\leq-\frac{C}{2}\theta\Delta_{E}e^{-\Delta_{E}\theta}\operatorname{tr}\left(XP_{0}^{\perp}e^{\theta
H_{0}}\right)$

where in the penultimate line we made use of the assumption in the statement
of the proposition. The result follows. ∎

We now bound the effect of perturbing the Hamiltonian on the Laplace
transform:

######  Proposition 12.

In the notations of Lemma 11 and Propositions 8 and 10, let $H=H_{0}+V$. Then:

$\displaystyle\mathfrak{L}_{\rho,H,\beta}(\theta,t)\leq
1+e^{-\frac{C}{2}\theta\Delta_{E}e^{-\Delta_{E}\theta}t+\theta\|H_{0}\|_{\infty}}+|t|\,3mM^{4}e^{\theta\|H\|_{\infty}}\max\left\\{\,e^{\frac{-\beta^{2}\Delta_{\nu}^{2}}{8}},\,\frac{4\beta\Delta_{E}}{2\beta\Delta_{E}-1}e^{-\frac{1}{4}\beta\Delta_{E}}\right\\}$
$\displaystyle\qquad+|t|m\beta
e^{\theta\|H\|_{\infty}}\Big{(}C_{1}\eta(\|H_{0}\|_{\infty}+\|V\|_{\infty})+C_{2}\|V\|_{\infty}(1+\ln(1/\eta))\Big{)}\,.$

for some universal constants $C_{1},C_{2}>0$, where $M$ is the number of
distinct eigenvalues of $H_{0}$, $\Delta_{E}$ the minimum gap between any two
energies of $H_{0}$, and $\Delta_{\nu}$ the minimum gap between any two Bohr
frequencies of $H_{0}$.

###### Proof.

By Proposition 8, we get that

$\displaystyle\mathfrak{L}_{\rho,H_{0},\beta}(\theta,t)$
$\displaystyle\leq\mathfrak{L}_{\rho,H_{0},\infty}(\theta,t)+\left\|e^{t\mathcal{L}_{\infty,H_{0}}}-e^{t\mathcal{L}_{\beta,H_{0}}}\right\|_{\infty\to\infty}\,e^{\theta\|H\|_{\infty}}$
$\displaystyle\leq
1+e^{-\frac{C}{2}\theta\Delta_{E}e^{-\Delta_{E}\theta}t+\theta\|H_{0}\|_{\infty}}+|t|\,3mM^{4}e^{\theta\|H\|_{\infty}}\max\left\\{\,e^{\frac{-\beta^{2}\Delta_{\nu}^{2}}{8}},\,\frac{4\beta\Delta_{E}}{2\beta\Delta_{E}-1}e^{-\frac{1}{4}\beta\Delta_{E}}\right\\}\,,$

where in the second ineqiality, we made use of the bound derived in Lemma 11.
Using the perturbation bound found in Proposition 10, we get

$\displaystyle\mathfrak{L}_{\rho,H,\beta}(\theta,t)$
$\displaystyle\leq\mathfrak{L}_{\rho,H_{0},\beta}(\theta,t)+\left\|e^{t\mathcal{L}_{\beta,H_{0}}}-e^{t\mathcal{L}_{\beta,H}}\right\|_{\infty\to\infty}\,e^{\theta\|H\|_{\infty}}$
$\displaystyle\leq
1+e^{-\frac{C}{2}\theta\Delta_{E}e^{-\Delta_{E}\theta}t+\theta\|H_{0}\|_{\infty}}+|t|\,3mM^{4}e^{\theta\|H\|_{\infty}}\max\left\\{\,e^{\frac{-\beta^{2}\Delta_{\nu}^{2}}{8}},\,\frac{4\beta\Delta_{E}}{2\beta\Delta_{E}-1}e^{-\frac{1}{4}\beta\Delta_{E}}\right\\}$
$\displaystyle\qquad+|t|m\beta
e^{\theta\|H\|_{\infty}}\Big{(}C_{1}\eta(\|H_{0}\|_{\infty}+\|V\|_{\infty})+C_{2}\|V\|_{\infty}(1+\ln(1/\eta))\Big{)}\,.$

∎

Making some additional assumptions on the spectrum of $H_{0}$ we immediately
obtain a bound on the overlap of a time-evolved state and the ground state:

######  Corollary 13.

In the same setting as Lemma 11 and Proposition 12, assume further that the
spectrum of $H_{0}$ is such that $M=\mathcal{O}(n^{a})$ for some exponent
$a>0$, that $H$ has a unique ground state $|\psi_{0}\rangle$ of energy zero,
and that $\Delta_{E}(H_{0}),\Delta_{\nu}(H_{0})=\Omega(1)$, norm
$\|H_{0}\|_{\infty}=\mathcal{O}(n)$, and $C=\Omega(n^{-b})$ for some positive
constant $b>0$ and $m=\mathcal{O}(n^{c})$ for $c>0$. Then for
$t=\Omega(n^{1+\epsilon})$, $k=2(1+2\epsilon+4a+c)$,
$\beta\geq\max\\{\Delta_{\nu}^{-1},\Delta_{E}^{-1}\\}k\log(n)$ and
$\|V\|_{\infty}=\mathcal{O}(n^{-(2+3\epsilon+4a+c)})$ for any $\epsilon>1$, we
have that, upon starting the quantum Monte Carlo algorithm in the maximally
mixed state $\tau$:

$\displaystyle\langle\psi_{0}|e^{t\mathcal{L}_{H,\beta}}(\tau)|\psi_{0}\rangle=\Omega(n^{-2\epsilon})\,.$
(37)

###### Proof.

From the conditions in the statement of the Corollary and by Proposition 12,
assuming that $\|V\|_{\infty}=o(1)$ and
$C^{\prime}\log(n)\geq\beta\geq\max\\{\Delta_{\nu}^{-1},\Delta_{E}^{-1}\\}k\log(n)$
for some constants $k,C^{\prime}>0$ that we will determine later, there are
positive constants $C_{i}$ such that

$\displaystyle\mathfrak{L}_{\rho,H,\beta}(\theta,t)\leq 1+C_{1}\,e^{\theta
C_{2}n}\left(e^{-\theta
C_{3}e^{-\Delta_{E}\theta}t}+tn^{4a+c}\beta\,\left\\{\frac{e^{-k\log(n)/2}}{2k-1}+\|V\|_{\infty}\,n\right\\}\right)\,.$

Choosing $\theta=(C_{2}n)^{-1}$, the above bound simplifies to

$\displaystyle\mathfrak{L}_{\rho,H,\beta}(\theta,t)\leq
1+C_{1}\,e\left(e^{-\frac{C_{3}}{C_{2}n}e^{-\frac{\Delta_{E}}{C_{2}n}}t}+tn^{4a+c}C^{\prime}\log(n)\,\left\\{\frac{n^{-\frac{k}{2}}}{2k-1}+\|V\|_{\infty}\,n\right\\}\right)\,.$

Next, choosing $t=\Omega(n^{1+\epsilon})$, $k=2(1+2\epsilon+4a+c)$ and
$\|V\|_{\infty}=\mathcal{O}(n^{-(2+3\epsilon+4a+c)})$ for any $\epsilon>0$, we
get that for some $K>0$

$\displaystyle\mathfrak{L}_{\rho,H,\beta}(\theta,t)\leq 1+Kn^{-\epsilon}\,.$

In words, the Laplace transform $\mathfrak{L}_{\rho,H,\beta}(\theta,t)$ can be
taken inverse polynomially close to $1$ using the bound from Prop. 12. Now,
choosing $\epsilon>1$ such that $e^{-\theta E_{1}}\leq 1-Kn^{-\epsilon}$
eventually, we get by a standard Markov inequality that

$\displaystyle
1-\langle\psi_{0}|e^{t\mathcal{L}_{H,\beta}}(\tau)|\psi_{0}\rangle$
$\displaystyle\leq e^{-\theta
E_{1}}\,\mathfrak{L}_{\rho,H,\beta}(\theta,t)\leq e^{-\theta
E_{1}}\Big{(}1+Kn^{-\epsilon}\Big{)}\leq\Big{(}1-Kn^{-\epsilon}\Big{)}\Big{(}1+Kn^{-\epsilon}\Big{)}=1-K^{2}n^{-2\epsilon}\,.$
(38)

In other words, we can ensure that the state at time $t$ has an inverse
polynomial overlap with the ground space of the perturbed Hamiltonian. ∎

###  C.3 Mapping quantum circuits to Gibbs states with circuit-to-Hamiltonian
constructions

We now show how to use the results of the last Section to conclude that the
class of dissipative evolutions we considered so far, when enlarged with the
ability to perform arbitrary few-qubit measurements, provides a model of
quantum computation that is polynomially equivalent to the circuit model for
$\beta=\Omega(\log(n))$. To achieve this, we resort to a slightly modified
circuit-to-Hamiltonian mapping, similar to the ones in [45, 59, 19]. Consider
an $n$-qubit circuit $C=U_{T}\dots U_{1}$ made of $T$ geometrically local
unitary gates $U_{1},\dots,U_{T}$, each operating on one or two neighbouring
qubits on a $2D$ lattice. For simplicity and without loss of generality, we
will assume that $T$ is divisible by $4$ in what follows. Next, we define the
following Hamiltonian on $n+T$ qubits:

$\displaystyle
H_{C,\lambda}:=H_{\operatorname{clock}}+\lambda(H_{\operatorname{in}}+H_{\operatorname{prop}})\,,$
(39)

where

$\displaystyle
H_{\operatorname{clock}}:=\sum_{t=1}^{T-1}\,I\otimes|01\rangle\langle
01|_{t,t+1}\equiv H_{0}$ $\displaystyle
H_{\operatorname{in}}:=\sum_{j=1}^{n}|1\rangle\langle
1|_{j}\otimes|10\rangle\langle 10|_{t_{j}-1,t_{j}}$ $\displaystyle
H_{\operatorname{prop}}:=\frac{1}{2}\,\sum_{t=1}^{T}H_{\operatorname{prop}}(t)\,,$

for some coupling $\lambda>0$ to be fixed later, and with

$\displaystyle
H_{\operatorname{prop}}(t):=I-h_{t}\big{(}U_{t}\otimes|110\rangle\langle
100|_{t-1,t,t+1}+U_{t}^{\dagger}\otimes|100\rangle\langle
110|_{t-1,t,t+1}\big{)}\qquad\text{ for each }1<t<T\,,$ $\displaystyle
H_{\operatorname{prop}}(1):=I-h_{1}(U_{1}\otimes|10\rangle\langle
00|_{1,2}+U_{1}^{\dagger}\otimes|00\rangle\langle 10|_{1,2})$ $\displaystyle
H_{\operatorname{prop}}(T):=I-h_{T}(U_{T}\otimes|11\rangle\langle
10|_{T-1,T}+U_{T}^{\dagger}\otimes|10\rangle\langle 11|_{T-1,T})\,,$

with coupling constants $h_{t}:=\sqrt{t(T-t+1)}$. Above, each of the times
$t_{j}$ in the definition of $H_{\operatorname{in}}$ corresponds to the first
times qubit $j$ is acted on. Clearly, the Hamiltonian $H_{0}$ has gap
$E_{1}=1$. Moreover, it was shown in [19] that $H_{C}$ has as unique ground
state 111In fact, the clock Hamiltonian constructed in [19] is slightly
different from the one we consider here, however it is easy to verify that
their ground states agree by a simple modification of the proof given in [19,
Appendix J.1].

$\displaystyle|\eta^{\prime}\rangle:=\sum_{t=0}^{T}\frac{1}{2^{T/2}}\sqrt{\binom{T}{t}}\,U_{t}\dots
U_{1}|0\rangle\otimes|0^{t}1^{T-t}\rangle\,,$

where $X_{k}$ denotes the Pauli $X$ operator acting on qubit $k$, whereas
$Z_{k}$ denotes the Pauli $Z$ operator acting on qubit $k$. At this point, we
choose the jump operators $A^{a}$ within the following set

$\displaystyle\big{\\{}A^{a}\big{\\}}_{a\in[2T-1]}:=\big{\\{}I\otimes
X_{t}\big{\\}}_{t=1}^{T}\cup\big{\\{}X_{t}\otimes
X_{t+1}\big{\\}}_{t=1}^{T-1}\,.$

The next Lemma is crucial to establish the constant $C$ of Lemma 11 when
$H_{0}=H_{\operatorname{clock}}$:

######  Lemma 14.

Let $T$ be divisible by $4$ and $s\in\\{0,1\\}^{T}$. Then:

  1. 1.

If $\bra{s}H_{\operatorname{clock}}\ket{s}<\frac{T}{4}$, then there is $i$
s.t. $X_{i}\ket{s}=\ket{x}$ and
$\bra{x}H_{\operatorname{clock}}\ket{x}=\bra{s}H_{\operatorname{clock}}\ket{s}+1$.

  2. 2.

If $\bra{s}H_{\operatorname{clock}}\ket{s}\geq\frac{T}{4}$, then there is $i$
s.t. either $X_{i}X_{i+1}\ket{s}=\ket{x}$ or $X_{i}\ket{s}=\ket{x}$, and
$\bra{x}H_{\operatorname{clock}}\ket{x}=\bra{s}H_{\operatorname{clock}}\ket{s}-1$.

###### Proof.

Let us show the first case. We claim that if
$\bra{s}H_{\operatorname{clock}}\ket{s}<\frac{T}{4}$ then $s$ must contain
$111$ or $000$ as a substring. Once we show this, the claim immediately
follows, as we can just flip the middle bit to increase the energy of the
clock by $1$. Now assume that the string does not contain $111$ or $000$. This
implies that at most after every four steps of the clock a transition from $0$
to $1$ needs to happen. Thus, the resulting energy would be at least $T/4$, a
contradiction. Let us now consider the second case. Let $k$ be the size of the
smallest block of consecutive $1s$ in the string. Assume $k\geq 4$. For each
such block we also need the string to contain at least one zero to have an
energy increase for the clock. But then each transition corresponds to a
substring of length at least $5$, which means that the total energy is at most
$\tfrac{T}{5}$. Thus, we conclude that $k\leq 3$, i.e. we have a substring of
the form $01^{k}0$ with $1\leq k\leq 3$. If $k=1,2$ we can then flip these $k$
$1s$ and erase one transition, decreasing the energy by one. If $k=3$, we
instead flip a $0$ between two strings of $1^{\prime}s$. ∎

Next, we use Lemma 14 in order to prove the assumption of Lemma 11 in the case
where $H_{0}=H_{\operatorname{clock}}$

######  Lemma 15.

Denoting the spectral decomposition
$H_{\operatorname{clock}}:=\sum_{i=1}^{T/2}iP_{i}$, for all $i\in[T/2]$, we
have

$\displaystyle\sum_{a}\operatorname{tr}\left(P_{i}A^{a}P_{i-1}A^{a}\right)\geq\min\left\\{1,\frac{6}{(T-1)^{2}}\right\\}\operatorname{tr}\left(P_{i}\right)\,.$
(40)

In other words, assumption 36 of Lemma 11 holds with
$C=\min\left\\{1,\,\frac{6}{(T-1)^{2}}\right\\}$.

###### Proof.

First, we argue that the dimension of the eigenspace associated to eigenvalue
$i$ of $H_{\operatorname{clock}}$ is of dimension

$\displaystyle\operatorname{tr}\left(P_{i}\right)=2^{n}\times\binom{T+1}{2i+1}\,.$

Clearly, the factor $2^{n}$ arises from the degeneracy provided by the
original system on which $H_{\operatorname{clock}}$ does not act. So without
loss of generality, we can assume that $n=0$. Then, we start by considering
the case $i=0$. Clearly, the ground space of $H_{\operatorname{clock}}$ is
spanned by the eigenstates $0^{T}$ and $1^{T}$ (there are
$2=2\times\binom{T-1}{0}$ of them), as well as any string of the form
$1^{t}0^{T-t}$ for $t\in[T-1]$ (there are $T-1=\binom{T-1}{1}$ of them).
Therefore, in total, the ground space is spanned by
$\operatorname{tr}\left(P_{0}\right)=2\binom{T-1}{0}+\binom{T-1}{1}=\binom{T+1}{1}$
eigenstates. Next we consider $i=1$: this is spanned by $\binom{T-1}{1}$
bitstrings of the form $0^{t}1^{T-t}$, $\binom{T-1}{2}$ bitstrings of the form
$0^{r}1^{s}0^{T-r-s}$, $\binom{T-1}{2}$ bitstrings of the form
$1^{r}0^{s}1^{T-r-s}$, and $\binom{T-1}{3}$ bitstrings of the form
$1^{r}0^{s}1^{t}0^{T-r-s-t}$, for a total of
$\operatorname{tr}\left(P_{1}\right)=\binom{T-1}{1}+2\binom{T-1}{2}+\binom{T-1}{3}=\binom{T+1}{3}$.
The general case follows by a simple iteration. From this, it is easy to
verify that for all $i\neq 0$,

$\displaystyle\operatorname{tr}\left(P_{i-1}\right)\geq\frac{6}{(T-1)^{2}}\operatorname{tr}\left(P_{i}\right).$

Next, we consider $i<T/4+1$. By Lemma 14, for any eigenstate
$|s\rangle=P_{i-1}|s\rangle$ there is a state
$|x_{s}\rangle=P_{i}|x_{s}\rangle$ and $t\in[T]$ such that $\langle
s|X_{t}|x_{s}\rangle=1$. Therefore

$\displaystyle\sum_{a}\operatorname{tr}\left(P_{i}A^{a}P_{i-1}A^{a}\right)\geq\sum_{|s\rangle=P_{i-1}|s\rangle}|\langle
s|X_{t}|x_{s}\rangle|^{2}=\operatorname{tr}\left(P_{i-1}\right)\geq\frac{6}{(T-1)^{2}}\,\operatorname{tr}\left(P_{i}\right)\,.$

Finally, we consider $i\geq T/4+1$. Again, by Lemma 14, for any
$\ket{s}=P_{i}\ket{s}$, there is $\ket{x_{s}}=P_{i-1}\ket{x_{s}}$ and $a$ such
that either $\langle s|A^{a}|x_{s}\rangle=1$ with $A^{a}=X_{t}$ or
$A^{a}=X_{t}X_{t-1}$. Therefore

$\displaystyle\sum_{a}\operatorname{tr}\left(P_{i}A^{a}P_{i-1}A^{a}\right)\geq\sum_{|s\rangle=P_{i}|s\rangle}|\langle
s|X_{t}|x_{s}\rangle|^{2}=\operatorname{tr}\left(P_{i}\right)\,.$

The result then follows.

∎

Combining Lemma 16 with Lemma 15, we directly get the inverse polynomial
overlap of states evolved with our Gibbs sampler during a polynomial time with
the history state $|\eta^{\prime}\rangle$:

######  Corollary 16.

For $t=\Omega((n+T)^{1+\epsilon})$, $\beta=\Omega(\log(n))$ and
$\|V\|_{\infty}=\mathcal{O}((n+T)^{-(7+3\epsilon)})$ for any $\epsilon>1$, we
have that, upon starting the quantum Monte Carlo algorithm in the maximally
mixed state $\tau$:

$\displaystyle\langle\psi_{0}|e^{t\mathcal{L}_{H,\beta}}(\tau)|\psi_{0}\rangle=\Omega((n+T)^{-2\epsilon})\,.$
(41)

To conclude, we show how to obtain the ground state by measuring a
frustration-free circuit-to-Hamiltonian mapping. For instance, we can consider
the same construction as before, but changing the $H_{\operatorname{prop}}(t)$
terms to

$\displaystyle H_{\operatorname{prop},FF}(t):=I\otimes(|100\rangle\langle
100|_{t-1,t,t+1}+|110\rangle\langle
110|_{t-1,t,t+1})-\big{(}U_{t}\otimes|110\rangle\langle
100|_{t-1,t,t+1}+U_{t}^{\dagger}\otimes|100\rangle\langle
110|_{t-1,t,t+1}\big{)},$

for $1<t<T$ and

$\displaystyle H_{\operatorname{prop},FF}(1):=I\otimes(|00\rangle\langle
00|_{1,2}+|10\rangle\langle 10|_{1,2})-(U_{1}\otimes|10\rangle\langle
00|_{1,2}+U_{1}^{\dagger}\otimes|00\rangle\langle 10|_{1,2})$ $\displaystyle
H_{\operatorname{prop},FF}(T):=I\otimes(|10\rangle\langle
10|_{T-1,T}+|11\rangle\langle 11|_{T-1,T})-(U_{T}\otimes|11\rangle\langle
10|_{T-1,T}+U_{T}^{\dagger}\otimes|10\rangle\langle 11|_{T-1,T})\,,$

Let us define, in analogy to what we had before:

$\displaystyle
H_{CFF,\lambda}:=H_{\operatorname{clock}}+\lambda(H_{\operatorname{in}}+H_{\operatorname{prop},FF})\,.$
(42)

We recall two important properties of this Hamiltonian:

######  Proposition 17.

Let $H_{CFF,\lambda}$ be defined as in Eq. (42). Then for all $\lambda>0$,
$H_{CFF,\lambda}$ is frustration-free and has the history state
$\ket{\eta}=\frac{1}{\sqrt{L+1}}\sum\limits_{l=0}^{L}\ket{\gamma_{l}}$ with
$\ket{\gamma_{l}}=(U_{t}\dots
U_{1})|0^{n}\rangle\otimes|0^{t}1^{T-t}\rangle\,$ as its ground state with
corresponding energy $0$. Furthermore, we have
$|\langle\eta^{\prime}|\eta\rangle|^{2}\geq\frac{1}{T+1}$.

###### Proof.

We refer to [45, Claim 3.3] for a proof of the first claim. In essence, the
authors show that $H_{CFF,\lambda}$ is a sum of positive semidefinite terms
and that $\ket{\eta}$ has energy $0$, from which it immediately follows that
$\ket{\eta}$ is the ground state. Furthermore, it is easy to see that:

$\displaystyle|\langle\eta^{\prime}|\eta\rangle|^{2}\geq\frac{1}{T+1}\sum_{t=0}^{T}\frac{\binom{T}{t}}{2^{T}}=\frac{1}{T+1}\,.$
(43)

∎

To show that a complexity class contains BQP, it suffices to show that we can
prepare $\ket{\eta}$ within that class in polynomial time for
$T=\textrm{poly}(n)$. However, we recall the following standard argument that
shows that it actually suffices to prepare a state with polynomial overlap
with $\ket{\eta}$ combined with the ability to measure the terms of
$H_{CFF,\lambda}$ to ensure that we can obtain $\ket{\eta}$ with a polynomial
overhead.

######  Proposition 18.

Let $H=\sum_{i=1}^{m}P_{i}$ with projections $P_{i}$ be a frustration-free
Hamiltonian with unique ground state $\ket{\psi}$ and $\rho$ a quantum state
such that $\langle\psi|\rho|\psi\rangle\geq p$. Suppose we measure the
sequence of POVMs $\\{E_{i}^{0}=P_{i},E_{i}^{1}=I-P_{i}\\}$ for $1\leq i\leq
m$. Then, conditioned on observing the sequence of outcomes $1^{m}$, the final
state is $\ket{\psi}$. Furthermore, we observe $1^{m}$ with probability at
least $p$.

###### Proof.

The proof is somewhat standard, but we repeat some details here for
convenience. The fact that if we measure $1^{n}$ we have prepared the state
follows from the fact that, as the Hamiltonian is frustration-free and the
ground state is unique,
$\bigcap\limits_{i=1}^{m}\textrm{kern}(P_{i})=\textrm{span}(\ket{\psi})$. Now
for the second claim, we note that we can write the probability of observing
the sequence $1^{m}$ as:

$\displaystyle\operatorname{tr}\left(E_{1}^{1}\rho
E_{1}^{1}\right)\times\frac{\operatorname{tr}\left(E_{2}^{1}E_{1}^{1}\rho
E_{1}^{1}E_{2}^{1}\right)}{\operatorname{tr}\left(E_{1}^{1}\rho
E_{1}^{1}\right)}\times\cdots\frac{\operatorname{tr}\left(\left[\prod_{i=1}^{m}E_{i}^{1}\right]^{\dagger}\rho\left[\prod_{i=1}^{m}E_{i}^{1}\right]\right)}{\operatorname{tr}\left(\left[\prod_{i=1}^{m-1}E_{i}^{1}\right]^{\dagger}\rho\left[\prod_{i=1}^{m-1}E_{i}^{1}\right]\right)}$
$\displaystyle=\operatorname{tr}\left(\left[\prod_{i=1}^{m}E_{i}^{1}\right]^{\dagger}\rho\left[\prod_{i=1}^{m}E_{i}^{1}\right]\right)$
$\displaystyle=\langle\psi|\rho|\psi\rangle\geq p\,,$

where we again used that $\ket{\psi}$ is the unique state in the kernel. ∎

We now have gathered all the technical ingredients required to discuss our
scheme for universal quantum computation at $\beta=\Omega(\log(n))$
temperature, GibbsQP. We assume that for $5$-local Hamiltonians 222in
principle they only need to be local w.r.t. a fully connected graph, so that
we can map any quantum circuit to this graph with a polynomial overhead we can
implement the Lindbladians $\mathcal{L}_{\beta,H_{C,\lambda}}$ with the
Metropolis filter function, $\beta=\Omega(\log(n))$, for polynomial times and
up to diamond norm that is inversely polynomially small. Furthermore, we
assume we can implement arbitrary $5$-qubit projective measurements.

######  Theorem 5.

$\mathsf{BQP}=\mathsf{GibbsQP}$.

###### Proof.

The inclusion of $\mathsf{GibbsQP}$ in $\mathsf{BQP}$ follows from the work of
[13], where they showed how to efficiently simulate these Lindbladian
evolutions for polynomial times and inverse temperatures. To see the other
direction, given any quantum circuit $C$ on $n$ qubits with a polynomial
number of gates $T$, let us consider the $5$ local Hamiltonian
$H_{CFF,\lambda}$ we obtain from Kitaev’s circuit to Hamiltonian construction
discussed before for $\lambda=\mathcal{O}((n+T)^{-(8+3\epsilon)})$. According
to the definition of $\mathsf{GibbsQP}$ we can simulate
$e^{t\mathcal{L}_{\beta,H_{C,\lambda}}}$ up to an error that is inversely
polynomially small, and we can initiate the system to the maximally mixed
state.

According to Cor. 16 we can obtain a state with a polynomial overlap with
$\ket{\eta^{\prime}}$ by evolving for a polynomial time. It then follows from
Prop. 17 that this state also has an overlap with $\ket{\eta}$ which is at
least polynomially small. Thus, combining it with Prop. 18, we conclude that
we have at least an inverse polynomially small probability of accepting all
the terms when we measure $H_{CFF,\lambda}$. Thus, by repeating this procedure
polynomially many times we prepare the state $\ket{\eta}$ with a high
probability of success. ∎