††thanks: Both authors contributed equally.††thanks: Both authors contributed
equally.

# Quantum advantage and stability to errors in analogue quantum simulators

Rahul Trivedi  Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1,
85748 Garching, Germany.  Munich Center for Quantum Science and Technology
(MCQST), Schellingstr. 4, D-80799 Munich, Germany.  Electrical and Computer
Engineering, University of Washington, Seattle, WA - 98195, USA  Adrian Franco
Rubio  Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, 85748
Garching, Germany.  Munich Center for Quantum Science and Technology (MCQST),
Schellingstr. 4, D-80799 Munich, Germany.  J. Ignacio Cirac  Max-Planck-
Institut für Quantenoptik, Hans-Kopfermann-Str. 1, 85748 Garching, Germany.
Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4,
D-80799 Munich, Germany.

###### Abstract

Several quantum hardware platforms, while being unable to perform fully fault-
tolerant quantum computation, can still be operated as analogue quantum
simulators for addressing many-body problems. However, due to the presence of
errors, it is not clear to what extent those devices can provide us with an
advantage with respect to classical computers. In this work we consider the
use of noisy analogue quantum simulators for computing physically relevant
properties of many-body systems both in equilibrium and undergoing dynamics.
We first formulate a system-size independent notion of stability against
extensive errors, which we prove for Gaussian fermion models, as well as for a
restricted class of spin systems. Remarkably, for the Gaussian fermion models,
our analysis shows the stability of critical (gapless) models at zero
temperature which have long-range correlations. Furthermore, we analyze how
this stability may lead to a quantum advantage, for the problem of computing
the thermodynamic limits of many-body models, in the presence of a constant
error rate and without any explicit error correction.

##  I Introduction

Quantum information processing systems hold the promise of solving a number of
problems in physics and computer science faster than their classical
counterparts [1, 2]. However, most quantum algorithms with theoretical
performance guarantees require a fault-tolerant quantum computer [3, 4, 5].
While in principle possible, implementing a fault tolerant quantum computer is
a technological challenge that could still take a long time to solve. This has
motivated several investigations trying to identify both quantum algorithms,
as well as physically relevant computational problems, that can be addressed
by quantum hardware in the near term and without any explicit error
correction.

Analog quantum simulators, wherein a target Hamiltonian is mimicked by an
experimentally controllable system, have shown some promise in solving
problems arising in many-body physics in the near term [6, 7, 8]. From a
theoretical standpoint, several previous works have developed the notion of
universal quantum simulators [9, 10, 11]. More practically, they offer several
distinct advantages in solving many-body problems. First, in many of those
problems, there is no need for a universal gate set to implement the target
Hamiltonian, and thus the requirements are much milder than those to build a
quantum computer. Furthermore, the simulation can be performed just by letting
the system evolve, without the need to Trotterize the evolution into a quantum
circuit, which typically incurs in a rapid proliferation of errors during the
algorithm [8, 12, 13, 14, 7]. Moreover, since one is typically interested in
the thermodynamic limit, this only requires obtaining the expectation values
of certain intensive quantities, like the energy or magnetization per lattice
site. Even if the quantum state of the simulator is orthogonal to the ideal
one, those quantities may still be correct, or have a small error, raising the
expectation that the degree of robustness required to be much lower than for
other problems in quantum computing.

These expectations make quantum simulators very promising in providing some
advantage with respect to classical computers when addressing typical quantum
many-body problems that appear in different branches of physics. However, the
standard notion of quantum advantage [15, 16, 17, 18, 19, 20], that deals with
the computational effort required to solve a problem as a function of the
system size, must be revisited since one may be interested in the
thermodynamic limit as opposed to a finite system [21, 22]. Additionally, the
fact that quantum simulators do not implement error correction casts doubts on
their applicability and has to be carefully taken into account. More
specifically, analog quantum simulators in most applications try to implement
a target Hamiltonian $H$ which is typically geometrically local i.e.
$H=\sum_{\alpha}h_{\alpha}$, with $\|{h_{\alpha}}\|\leq 1$. Each term
$h_{\alpha}$ is implemented by engineering local interactions between
different subsystems in the quantum simulator. In practice, however, each of
these local interactions are slightly different from their target.
Consequently, the quantum simulator implements a perturbed Hamiltonian
$H^{\prime}=H+\delta\sum_{\alpha}v_{\alpha}$, where $\|{v_{\alpha}}\|\leq 1$
and $\delta$ measures the “hardware error” incurred by the simulator locally.
While a well designed experimental setup can, in principle, achieve $\delta\ll
1$, the total error between $H^{\prime}$ and $H$ grows with $n$, the number of
qubits in the quantum simulator (i.e. the problem size) —
$\|{H^{\prime}-H}\|\sim\delta\text{poly}(n)$. Due to this proliferation of
errors, the error in the state of the quantum simulator, relative to the
target state, both for dynamics or thermal equilibrium, will only be small if
$\delta<O(1/\text{poly}(n))$ — this could potentially limit the applicability
to quantum simulators to only small-scale problems.

In this paper we show that many physically relevant problems are stable and
avoid this worst-case proliferation of errors, and thus can be addressed with
the help of quantum simulators without implementing any error correction. We
also propose a notion of quantum advantage, in the presence of errors, for
such problems, where the figure of merit is the computational time to obtain
an intensive quantity in the thermodynamic limit with a hardware-limited
precision.

In order to address the stability with respect to noise, we consider the
quantum simulation of both equilibrium and non-equilibrium problems. We first
study geometrically local Gaussian fermion models, and show that $k-$local
observables (or their weighted averages) are stable for the problem of
constant-time dynamics. Furthermore, when considering the ground state and
Gibbs state of Gaussian fermion models, we show that translationally invariant
$k-$local observables are stable without any restrictive assumptions on the
spectrum of the model — our results hold not only for gapped models, but also
for gapless models and without any assumption on the frustration freeness of
the model. Finally, we also consider spin systems with geometrically local
interactions, and show that well-known locality properties [23, 24, 25, 26,
27] of these models imply stability of local observables in constant-time
dynamics, ground states and Gibbs states but with more restrictive assumptions
on the model (e.g. the Hamiltonian being stably gapped for ground states, or
exhibiting exponential clustering of correlations for Gibbs states).
Additionally, for Gaussian fermion models, we numerically verify that the the
ground state problem is, in fact, stable to errors in both the gapped and
gapless cases, as is theoretically predicted by our stability analysis. Based
on these results, we hypothesize that many physically relevant, gapped or
gapless, models could potentially be within the reach of near-term analogue
quantum simulators.

Finally, we turn to the question of whether the class of problems that are
stable under errors on quantum simulators also provide some advantage over the
best known classical algorithms for solving the same problems. Here, instead
of framing the algorithm run-time as a function of system-size (i.e. the
number of spins/fermions) and the desired precision, we treat the system-size
as a truncation parameter introduced in approximating the more physically
relevant thermodynamic limit [21, 22]. This allows us to express the run-time
of the classical and quantum algorithms in terms of a single parameter — the
precision in the approximated thermodynamic limit. We argue that for many
problems for which we have evidence of stability to errors, we also expect an
advantage, with respect to precision, in using the quantum simulator — in the
absence of errors, for the problem of constant-time dynamics this advantage is
superpolynomial with respect to precision and exponential with respect to
evolution time in two or higher dimensional local Hamiltonians, and for ground
state computation it is exponential in desired precision. In the presence of
errors, there is a fundamental limit to the precision that can be achieved by
the quantum simulator which is determined by the hardware error — the
advantage of the quantum simulator can then be assessed by scalings of the
classical run-time required to achieve this hardware-limited precision with
the hardware error.

##  II Stability of analogue quantum simulators

Abstractly, an analogue quantum simulation can be considered to configure a
target Hamiltonian $H$ on $n-$spins and measure an observable $O$ on a state
$\rho_{H}$ associated with this Hamiltonian (e.g. a state obtained from an
initial product state after evolution for finite time under $H$, or the Gibbs
state/ground state of $H$). The target Hamiltonian can be expressed as a sum
of terms each of which correspond to interactions between groups of spins i.e.

$\displaystyle H=\sum_{\Lambda}h_{\Lambda}.$ (1)

In an experimentally realistic setting, in the process of configuring this
target Hamiltonian, an error can be incurred in each term i.e. the Hamiltonian
that is implemented on the hardware is instead given by

$\displaystyle H^{\prime}=\sum_{\Lambda}h_{\Lambda}^{\prime},$ (2)

where we are only guaranteed that
$\|{h_{\Lambda}-h_{\Lambda}^{\prime}}\|\leq\delta$ for all $\Lambda$ and for
some $\delta$ — $\delta$ can be considered to be a measure of the ‘hardware’
error incurred in the quantum simulator. We note that since the number of
terms in the Hamiltonian in Eqs. 1 and 2 will typically be
$\textnormal{poly}(n)$, we can only guarantee that
$\|{H-H^{\prime}}\|\leq\delta\textnormal{poly}(n)$ — consequently, for the
analogue quantum simulator to be accurate without error correction, we would
require the hardware error to scale as $1/\text{poly}(n)$. This would
typically make it intractable to obtain thermodynamic limits accurately with
analogue quantum simulators since the hardware error would have to be made
very small.

Figure 1: Schematic depiction of our error model for analogue quantum
simulator. A target Hamiltonian $H$, expressed as sum of Hamiltonian terms
modelling interactions between groups of spins, when implemented on an
analogue quantum simulator would have an error $\delta$ per term. The
parameter $\delta$ thus controls the accuracy of the hardware implementation.

An alternative viewpoint would be to ask if there are certain interesting
many-body problems for which a good estimate for the thermodynamic limit can
be produced with a hardware with constant errors. This motivates us to look
for ‘quantum simulation tasks’ which are stable to an error in each term of
the Hamiltonian, as made precise in the following definition.

######  Definition 1 (Stable quantum simulation task).

The quantum simulation task on $n$ spins of measuring an observable $O_{n}$ in
a state $\rho_{H_{n}}$ associated with a target Hamiltonian $H_{n}$ is said to
be stable if the corresponding state $\rho_{H_{n}^{\prime}}$ for the perturbed
Hamiltonian $H_{n}^{\prime}$ satisfies

$|{\textnormal{Tr}(O_{n}\rho_{H_{n}})-\textnormal{Tr}(O_{n}\rho_{H_{n}^{\prime}})}|\leq
f(\delta),$

for some $f$, independent of $n$, such that $f(\delta)\to 0$ as $\delta\to 0$.

Thus, if a quantum simulation task is stable as per this definition, we can
hope to be able to estimate the thermodynamic limit of the observable on a
quantum simulator with the hardware error determined entirely by the precision
desired in the thermodynamic limit. These problems would not require the
hardware error to be scaled down with system size even in the absence of error
correction, and we can consider them to be problems that analogue quantum
simulators can conceivable solve in the near term.

The stability of quantum simulation tasks for many-body physics problems,
while hard to prove rigorously, can be intuitively understood. For simplicity,
consider the specific setting in which both the target Hamiltonian $H$ and
$V=H^{\prime}-H$ are translationally invariant — in this case, if the
thermodynamic limit of an observable (e.g. a local observable) $O$ in a non-
equilibrium or equilibrium state associated with the translationally invariant
Hamiltonian $H(s)=H+sV$ exists and is a continuous function of $s$, then it
can be seen that it is indeed stable in the sense of definition 1. However,
this argument falls short of a full proof since in an actual experiment, even
for target Hamiltonians that are translationally invariant, the presence of
errors can make it translationally varying thus making it hard to define its
thermodynamic limit. In particular, disorder can induce Anderson or many-body
localization which can even result in local observables being unstable. In the
following sections, however, we show that for several many-body problems,
physically interesting observables and order parameters are stable even to
translationally varying errors in the Hamiltonian.

##  III Stability of Gaussian fermion models

To make further progress on whether many-body simulation tasks can be
considered to be stable in the sense defined above, we consider physically
relevant setting of free-fermion models. We will consider fermions arranged on
a $d-$dimensional lattice with $L$ sites in each direction
$\mathbb{Z}_{L}^{d}$, and at each site we have $D$ fermionic modes — we denote
by $c_{x}^{\alpha}$ for $x\in\mathbb{Z}_{L}^{d},\alpha\in\\{1,2\dots 2D\\}$
the Majorana operators associated with each site $x$. We consider a target
Hamiltonian with local interactions given by

$H=\sum_{\begin{subarray}{c}x,y\in\mathbb{Z}_{L}^{d}\\\ d(x,y)\leq
R\end{subarray}}\sum_{\alpha,\beta=1}^{2D}h^{\alpha,\beta}_{x,y}c_{x}^{\alpha}c_{y}^{\beta}.$
(3)

where $R$ is the range of interaction, and we assume that all interaction
strengths are bounded by a constant, $|h^{\alpha,\beta}_{x,y}|\leq J$. Due to
local hardware errors, the quantum simulator instead implements a perturbed
free-fermion Hamiltonian $H^{\prime}$,

$H^{\prime}=\sum_{\begin{subarray}{c}x,y\in\mathbb{Z}_{L}^{d}\\\ d(x,y)\leq
R\end{subarray}}\sum_{\alpha,\beta=1}^{2D}h^{\prime\alpha,\beta}_{x,y}c_{x}^{\alpha}c_{y}^{\beta},$
(4)

such that

$|h^{\alpha,\beta}_{x,y}-h^{\prime\alpha,\beta}_{x,y}|\leq\delta.$

Generically, we expect $\|{H-H^{\prime}}\|=\Theta(n)$, where $n=2DL^{d}$ is
the total number of fermions. However, in the following subsections, we show
the stability of the quantum simulation task of measuring several important
classes of observables, which include typical intensive order parameters, in
both dynamical and equilibrium states of these models.

###  III.1 Finite-time dynamics

We are interested in the expectation value of Gaussian observables $O$ which
are either $k-$local, i.e. they act on a set
$\mathcal{S}\subseteq\mathbb{Z}_{L}^{d}$ of $k$ sites

$\displaystyle
O_{0}=\sum_{x,y\in\mathcal{S}}\sum_{\alpha,\beta=1}^{2D}o_{x,y}^{\alpha,\beta}c_{x}^{\alpha}c_{y}^{\beta},$
(5)

or weighted averages of $k-$local Gaussian observables i.e. are of the form
$\sum_{i=1}^{M}w_{i}O_{i}$, where $O_{i}$ is of the form Eq. 5,
$\sum_{i=1}^{M}|{w_{i}}|=1$ and $M$ can possibly grow with $n\sim L^{d}$.

We consider an arbitrary Gaussian initial state, and let the target state
$\rho_{H}$ be the state obtained on evolving it with $H$ for time $t$. We show
in Appendix B that

######  Proposition 1.

The quantum simulation task of measuring $k-$local Gaussian observables, or
their weighted sums, for constant-time dynamics under a spatially local
Gaussian Hamiltonian is stable with $f(\delta)=O(t\delta)$.

We point out that the dependence of the error between the observable in
perturbed and unperturbed models on $t$ is independent of the dimensionality
of the lattice $d$ — this result is thus stronger than what would be expected
simply from locality, wherein the error would be expected to grow as
$t\times(\text{Number of sites in the light cone at time }t)\propto t^{d+1}$ —
we revisit this in section IV.

###  III.2 Equilibrium

We next study the stability properties of the ground state as well as the
Gibbs state. We again consider the target Hamiltonian described while studying
the stability of constant-time dynamics, but make an additional physically
reasonable assumption on the density of modes of the target Hamiltonian $H$.

######  Assumption 1.

We assume that the number of eigenfrequencies $n_{\eta}$ of $H$, which are
eigenvalues of the matrix $h_{x,y}^{\alpha,\beta}$ defining the target
Hamiltonian, lying in the interval $[-\eta,\eta]$ satisfy the upper bound

$\displaystyle n_{\eta}\leq nf_{h}(\eta)+\kappa(\eta,n),$ (6)

where $n=DL^{d}$ is the number of fermionic modes, $f_{h}(\eta)$ is some
function which $\to 0$ as $\eta\to 0$ and $\kappa(\eta,n)$ is $o(n)$ for any
fixed $\eta$.

Physically, this assumption formalizes the expectation that any frequency
interval would have an extensive number of eigenmodes lying within it — we
expect this to be true for most physically relevant models. In particular, for
translationally invariant models, such an estimate would be true with
$f_{h}(\eta)$ being governed by the derivatives of the dispersion relation of
the free-fermion bounds near $\omega=0$.

The observables we consider while analyzing ground states are translationally
invariant Gaussian observables generated $k-$locally, i.e., if $O_{0}$ is a
$k-$local observable of the form of Eq. 5, then we consider observables of the
form

$O=\frac{1}{n}\sum_{x\in\mathbb{Z}_{L}^{d}}\tau_{x}(O_{0}),$

where $\tau_{x}$ is the observable $O_{0}$ translated by $x$. Considering the
target state $\rho_{H}$ to be the ground state of the free-fermion
Hamiltonian, we then obtain the following proposition (proved in appendix C)

######  Proposition 2.

The quantum simulation task of measuring translationally invariant Gaussian
observables generated $k-$locally, in the ground state of a spatially local
free-fermion model whose density of modes satisfies Eq. 6 is stable with
$f(\delta)=O(\sqrt{\delta})+f_{h}(O(\delta^{1/4}))$.

We note that the stability result above holds very generally, with only a mild
assumption on the density of modes of the free-fermion model. In particular,
it holds for free-fermion models which are not gapped, i.e. the energy
separation between the ground state and the first excited state decreases with
$n$, as long as the observable under consideration is translationally
invariant.

The translational invariance of the observables considered here is key to
obtaining the stability result for ground states — translationally varying
observables need not be stable, even if they are intensive and spatially
local. A simple counter-example here is that of Anderson localization —
consider $H$ to be a 1D translationally invariant tight-binding model i.e.
$H=\sum_{i=1}^{n-1}(a_{i+1}^{\dagger}a_{i}+a_{i}^{\dagger}a_{i+1})$, with
errors $\sum_{i=1}^{n}\delta_{i}a_{i}^{\dagger}a_{i}$ where $v_{i}$ is chosen
uniformly at random between $[-\delta,\delta]$. In the absence of error, the
ground state of $H$ is completely delocalized across the spin-chain. In the
presence of errors, no matter how small, this model is known to be localized.
Now, for every $\delta_{1},\delta_{2}\dots\delta_{n}$ and , consider the
translationally varying intensive observable
$O_{\delta_{1},\delta_{2}\dots\delta_{n}}$ given by the average particle
numbers on $\Theta(1/\delta)$ sites around the site where the ground state is
localized. This observable, when measured in the delocalized ground state of
the unperturbed Hamiltonian $H$, yields an expected value of $0$ as
$n\to\infty$. On the other hand, in the ground state of the perturbed
localized model it will yield an expected value of $\Theta(1)$. Thus, not all
translationally varying observables can be stable, even if we restrict
ourselves to free-fermion models, with the observables being intensive and
spatially local.

Finally, we consider the Gibbs state, $\rho_{H}=e^{-\beta
H}/\text{Tr}(e^{-\beta H})$ where the inverse-temperature $\beta$ is a
constant independent of $n$, corresponding to the spatially local free-fermion
Hamiltonian above, and again study the stability of translationally invariant
Gaussian observables that are generated $k-$locally. We show in appendix D
that

######  Proposition 3.

The quantum simulation task of measuring translationally invariant Gaussian
observables generated $k-$locally, in the Gibbs state at inverse-temperature
$\beta$ of a spatially local free-fermion is stable with
$f(\delta)=O(\beta\sqrt{\delta})$.

We point out that, in contrast to the corresponding result for ground state,
this stability result corresponding to the Gibbs state does not rely on an
assumption on the density of modes of the target Hamiltonian. However,
$f(\delta)$ grows with $\beta$, so this result does not directly imply the
stability of the ground state since $\beta$ would in general have to be
increased with $n$ for the Gibbs state to approximate the ground state.

As an illustrative example to corroborate proposition 2, we numerically study
the ground state of the free fermionic Su–Schrieffer–Heeger (SSH) model on $n$
fermiions with periodic boundary condtion:

$H_{\text{SSH}}[J]\equiv\sum_{i=1}^{n}{t_{i}a_{i}^{\dagger}a^{\phantom{\dagger}}_{i+1}}+\text{H.c},~{}~{}~{}~{}~{}t_{i}=\begin{cases}1&i\text{
odd,}\\\ J&i\text{ even.}\end{cases}$ (7)

where $a_{n+1}\equiv a_{n}$. This model displays a (topological) phase
transition at $J=1$, where the gap closes as $1/n$, and is gapped otherwise.

Figure 2: Numerical study of the error in intensive translationally varying
observable in the SSH model. The observable that we study here is
$O=H_{\text{SSH}}[J]/n$, where $H_{\text{SSH}}[J]$ is the Hamiltonian of the
ideal SSH model (Eq. 7) (a) The error in the expected value of the observable
$O$ in the ground state between the perturbed and unperturbed Hamiltonians, as
a function of $\delta$, the hardware-error, and the number of sites $n$. For
both gapped ($J=0.5,1.5$) and gapless ($J=1.0$) cases, we see that the error
in $O$ becomes independent of $n$ as $N\to\infty$. (b) Numerically extracted
error between the perturbed and unperturbed models for $n\to\infty$ as a
function of $\delta$, and its fit with $\delta^{2}$. (c) The error between the
perturbed and unperturbed model as a function of $J$ — for the same hardware
error $\delta$, this error peaks at $J=1$ which is also the point at which the
gap in the unperturbed model closes. All the errors are computed by averaging
over 500 random instances of perturbed models.

We consider measuring the energy density $H_{\text{SSH}}[J]/n$ of the
unperturbed Hamiltonian in the presence and absence of errors. Figure 2(a)
shows impact of changing system size on this energy density — we see that for
both gapped ($J=0.5,1.5$) and gapless ($J=1.0$) cases, the errors in the
energy density becomes independent of $n$ as $n\to\infty$, verifying the
expectation in proposition 2. Furthermore, we show the error in the energy
density for large $n$ as a function of $\delta$ in Fig. 2(b) and see that,
consistent with proposition 2, this error $\to 0$ as $\delta\to 0$. Finally,
Fig. 2(c) shows this error as a function of $J$ — we see the error peak near
$J=1$ (i.e. the point where the gap in the Hamiltonian closes), and that it is
smaller for values of $J$ where the model is gapped.

##  IV Stability in quantum spin systems

While for free fermion models, we could prove tight stability results with
minimal assumptions on the model, looser stability results hold for quantum
spin systems under more restrictive assumptions on their many-body spectrum.
In this section, we outline several quantum simulation tasks whose stability
is a consequence of several locality results that have already been
established in the many-body literature [23, 24, 25, 28, 29, 30, 26, 27]. Our
emphasis here is to cross-examine these results from the standpoint of quantum
simulation, and discuss their implications on noisy quantum simulators.

We restrict ourselves to spin systems that are spatially local — $n$ spins are
arranged on a lattice $\mathfrak{L}\subset\mathbb{Z}^{d}$ in $d-$dimensions,
and the Hamiltonian is of the form

$H=\sum_{x\in\mathfrak{L}}h_{x},$

where $h_{x}$ only acts on spins within a distance $R$ from $x$, and
$\|{h_{x}}\|\leq J$ for all $x$. We consider different physically relevant
observables, and states (either ground states, Gibbs states or states
generated by constant-time dynamics) that are associated with $H$, and study
the stability of the quantum simulation problem of computing a local
observable in these states.

###  IV.1 Finite time dynamics

Consider first the setting where an initial state $(\ket{0}\bra{0})^{\otimes
n}$ is evolved under the Hamiltonian $H$ for a time $t$ that is independent of
$n$ — we are thus interested in $\rho_{H}=e^{-iHt}(\ket{0}\bra{0})^{\otimes
n}e^{iHt}$. We consider observables $O$ that are either local (i.e. they only
act non-trivially on an $n-$independent subset of spins), or of the form

$\displaystyle O=O_{1}O_{2}\dots O_{k},$ (8)

where $O_{1},O_{2}\dots O_{k}$ are local and $k$ is independent of $n$, or

$O=\sum_{i=1}^{M}w_{i}O_{i},$

where $\sum_{i=1}^{M}|{w_{i}}|=1$, $O_{i}$ are of the form of Eq. 8 and $M$
can possibly grow with $n$. For these observables, the stability of this
quantum simulation task can be stated:

######  Proposition 4.

The quantum simulation task of measuring $k-$local observables, or their
weighted averages, for constant-time dynamics under a spatially local
Hamiltonian is stable with
$f(\delta)=O(t^{d+1}\delta)+O(t\delta\log^{d}(\delta^{-1}))$.

The proof of this result, provided in appendix E is straightforward and uses
the Lieb-Robinson’s bounds [23, 24] to approximate the Heisenberg picture
evolution of local observables with that corresponding to the Hamiltonian
truncated within their light cones, and then uses a perturbation theory bound
on this truncated Hamiltonian.

Note also that for large $t$, the error between the target observable and the
observable measured on the quantum simulator grows as $t^{d+1}$ — this is
looser than the corresponding result in Gaussian fermion models (proposition
1), where the error grows only as $t$. Furthermore, since this error bound
becomes loose with $t$, it prevents us from using this result to understand
the stability of quantum simulation tasks which are aimed to studying the
ground state properties of many-body Hamiltonians, and use the adiabatic
algorithms [7] that evolve a Hamiltonian for $t\sim\text{poly}(n)$. To address
these problems, we consider the stability of the ground state and Gibbs states
of local Hamiltonians using a different approach.

###  IV.2 Equilibrium

We next study the stability of the task of simulating the ground state and
Gibbs states of $H$. We first consider the problem of measuring $k-$local
observables, and their weighted averages, in the ground state. We assume that
$H$ is gapped i.e. the energy difference between the ground state and the
first excited state is $\geq\Delta$ independent of $n$. Furthermore, we also
assume that the Hamiltonian remains gapped in the presence of errors i.e. any
Hamiltonian $H^{\prime}$ such that

$\displaystyle H^{\prime}=\sum_{x\in\mathfrak{L}}h_{x}^{\prime}\text{ with
}\|{h_{x}-h_{x}^{\prime}}\|\leq\delta\ \text{for all }x,$ (9)

is gapped with gap $\Delta$ — we refer to such a target Hamiltonian $H$ to be
_stably gapped_ with gap $\Delta$. We point out that the stability of the gap
in the presence of errors or perturbations has only been shown for certain
frustration free models with local topological order [31, 32, 33, 34],
although we posit it as a reasonable physical assumption. The stability of
this quantum simulation task is a direct consequence of the spectral flow
method, or quasi-adiabatic continuation, proposed by Hastings and co-workers
[25, 26] which shows that map taking ground state of $H$ to the ground state
of $H^{\prime}$ can be well-approximated by a local unitary. We thus obtain
the following proposition and we include a proof of this in appendix F.

######  Proposition 5.

The quantum simulation task of measuring $k$-local observables, or their
weighted averages, in the ground state of stably gapped spatially local
Hamiltonians is stable with $f(\delta)=O(\delta)$.

We point out, again, that the choice of observables here is crucial to having
a stable quantum simulation problem — even for stably gapped Hamiltonians,
non-local observables would in general not yield a stable quantum simulation
task. A counter-example here would be the unperturbed Hamiltonian
$H=\sum_{i=1}^{n}Z_{i}$, and a perturbed Hamiltonian
$H^{\prime}=\sum_{i=1}^{n}Z_{i}+\delta X_{i}$. It can be seen that both these
models are gapped for small $\delta$ — the gap of $H$ is 2, and gap of
$H^{\prime}$ is $2\sqrt{1+\delta^{2}}$. However, the non-local observable
$O=(\ket{0}\bra{0})^{\otimes n}$, when measured in the ground state of $H$ is
1 and when measured in the ground state of $H^{\prime}$ evaluates to
$(1+\delta^{2})^{-n/2}$ — thus, at any $\delta>0$, the error in this
observable $\to 1$ as $n\to\infty$ showing that it is unstable in the sense of
definition 1.

We next consider the Gibbs state of $H$ at some temperature $\beta$
independent of $n$, and assume that the Gibbs state has an exponential
clustering of correlation [27] i.e. for any two observables $A$, $B$ separated
by distance $l$,

$\left|{\langle A\otimes B\rangle-\langle A\rangle\langle
B\rangle}\right|\leq\|{A}\|\|{B}\|O(e^{-c_{2}l}),$

for some constant $c_{2}$. Furthermore, as in the case of ground states, we
assume that this exponential clustering of correlations is stable i.e. for any
perturbation $H^{\prime}$ of $H$ of the form of Eq. 9, this exponential
clustering property holds. In this case, we again obtain that the problem of
measuring $1-$local observables and their weighted averages is stable.

######  Proposition 6.

The quantum simulation task of measuring 1-local observables, or their
weighted averages, in the Gibbs state of spatially local Hamiltonians with
stable exponential clustering of correlations is stable with
$f(\delta)=O(\log^{1-1/d}(1/\delta)e^{-\Omega(\log^{1/d}(1/\delta))}\big{)}$.

##  V Discussion: Quantum advantage in the thermodynamic limit

In many problems in physics, the quantities of interest are the value of
certain intensive observables in the thermodynamic limit i.e. when the system
size $n\to\infty$. Thus, typical complexity arguments behind quantum advantage
that rely on the scaling of the computational effort with the system size have
to be revised and adapted to this limit. One way of doing that is to
concentrate on the computational complexity of obtaining the thermodynamic
limit as a function of a prescribed error, $\varepsilon$ [21, 22]. A quantum
advantage can then be obtained if the computational time on a quantum computer
required to reach a precision of $\varepsilon$ scales more favorably with it
than with a classical computer. For instance, in a quantum computer this
scaling can be $T_{\rm Q}={\rm poly}(1/\varepsilon)$ while in a classical one
$T_{\rm cl}={\rm exp}({\log^{2}(1/\varepsilon}))$ or $T_{\rm cl}={\rm
exp}({O(1/\varepsilon)})$, in which case we will respectively have a
superpolynomial or exponential quantum advantage. Furthermore, even this
notion of advantage relies on error-free computation — for quantum simulation
in the presence of errors, a new notion of quantum advantage needs to be
introduced.

In this section we discuss potential advantage in quantum simulations for
dynamics and ground state problems. First, we will consider the ideal (and
unrealistic) situation where there are no errors. For dynamical problems we
will show that a quantum advantage can be reached for the problem of computing
the thermodynamic limit in the sense described above. For ground state
problems, we will adopt a different perspective than the one analyzed in Refs.
[21, 22], where they consider the most general scenario and even a quantum
computer may require an exponential time in terms of the required error.
Instead, we will consider promise problems where the many-body Hamiltonian
fulfills physically reasonable assumptions and discuss potential quantum
advantage. Most of our discussion is based on technical results related to
locality estimates and analysis of adiabatic quantum algorithms that can
already be found in literature — we focus on describing the implications of
these results for quantum simulation.

Next, we will consider the more realistic scenario where the quantum simulator
is subject to errors as analyzed in the previous sections. There, we will
argue that quantum advantage has to be based on the classical computational
time required to obtain the value of the intensive observable in the
thermodynamic limit with the same error as with the quantum simulator. Within
this framework, we will rely on the robustness results of the previous
sections, to argue that quantum advantage can persist for stable quantum
simulation tasks even in the presence of errors.

###  V.1 Quantum advantage for ideal quantum simulators

Let us consider a many-body model defined as a family of Hamiltonians
$\\{H_{n}\\}_{n\in\mathbb{N}}$ and observables $\\{O_{n}\\}_{n\in\mathbb{N}}$,
where $H_{n},O_{n}$ act on $n-$spins. We are interested in the expected value
of $O_{n}$ in a many body quantum state generated by the Hamiltonian $H_{n}$,
$\rho_{H_{n}}$. We furthermore assume that the models and observables under
consideration have a well-defined thermodynamic limit i.e.

$O^{*}:=\lim_{n\to\infty}\text{Tr}(\rho_{H_{n}}O_{n})$ (10)

exists.

Examples of physically relevant $\rho_{H_{n}}$ would include the states
obtained on evolving an initial state of $n$ spins after some time $t$, the
Gibbs state corresponding to $H_{n}$ or the ground state of $H_{n}$. We will
mostly focus on families of local Hamiltonians

$\displaystyle H_{n}=\sum_{x\in\mathbb{Z}_{L_{n}}^{d}}h_{x},$ (11)

where $h_{x}$ acts on spins in a cube of unit length with $x$ being its lower
left corner and satisfies $\|{h_{x}}\|\leq J$ for some $J>0$. Typically, one
considers translationally invariant models, where $h_{x}$ is just some $h_{0}$
translated to the position $x$. We will also only consider here observables
$O_{n}$ are then acting on a single or few lattice sites.

Now, given a precision $\varepsilon$, we can then choose $n(\varepsilon)$ such
that

$|{O^{*}-\textnormal{Tr}(O_{n}\rho_{H_{n}})}|\leq\varepsilon,$ (12)

i.e. approximate the thermodynamic limit by a finite-size problem. The run-
time of a quantum simulation or a classical simulation for the finite-size
problem can thus be expressed in terms of the precision $\varepsilon$ demanded
in the thermodynamic limit. This allows us to then compare the scalings of the
run-time of these algorithms with the precision $\varepsilon$, and declare an
algorithm _to have an advantage in precision_ compared to others depending on
their respective scaling.

We note that there are classical algorithms which directly operate in the
thermodynamic limit [35, 36]. We will not consider them here, since for ground
state problem in dimensions higher than two, it is not possible to give
rigorous scalings, and for dynamical ones they give the same scaling as the
ones considered here.

Figure 3: An erroneous quantum simulator can obtain the thermodynamic limit of
the perturbed model to a precision $\varepsilon^{\prime}$ in time
$\text{poly}(1/\varepsilon^{\prime})$ — this thermodynamic limit, however, can
have an error $f(\delta)$ from the target thermodynamic limit in the presence
of hardware error $\delta$.

####  V.1.1 Finite time quantum dynamics

We consider first an initial product state $\ket{0}^{\otimes n}$, and for
$t>0$, we take $\rho_{H_{n}}=e^{-iH_{n}t}\big{(}\ket{0}\bra{0}\big{)}^{\otimes
n}e^{iH_{n}t}$. The observable of interest is a _local observable_ $O$.

The existence of the thermodynamic limit is obtained directly using the Lieb-
Robinson bounds [23, 24], and it also characterizes the error between the
thermodynamic limit and its finite-size approximation.

######  Lemma 1.

For the problem of computing a local observable $O^{*}$ after evolving
$\ket{0}^{\otimes n}$ for finite-time $t$ with respect to a nearest-neighbour
Hamiltonian exists and fulfills Eq. 12 for
$n=\Theta(\log^{d}(1/\varepsilon)+t^{d})$.

On a quantum simulator, one would evolve
$n=\Theta(\log^{d}(1/\varepsilon)+t^{d})$ qubits for a time $t$ and measure
the observable. The procedure would be repeated $M$ times, so that the error
in the observable would decrease as $1/M^{2}$. In order to obtain an error
$\varepsilon$, the computation time would then be $O(t/\varepsilon^{2})$

On a classical computer, in general, we would have to compute the quantum
state on $O(\log^{d}(1/\varepsilon)+t^{d})$ spins. Using exact methods, this
would require a run-time that scales as
$2^{O(\log^{d}(\varepsilon^{-1}))+t^{d}}$. For a fixed time $t$, in two or
higher dimension, this run-time would be superpolynomial in $1/\varepsilon$
and exponential in $t$. One could also use other classical methods that obtain
a better scaling with the required precision. For instance, in Ref. [37] a
method based on cluster expansion was analyzed and which the computational
time scales as poly$(n)$, although super-exponentially with time. Compared to
such a classical algorithm, the quantum advantage is super-exponential with
respect to evolution time $t$.

####  V.1.2 Ground state

Consider next the problem of estimating local observables in the ground state
of many-body Hamiltonians in the thermodynamic limit. The convergence rate of
a finite-size approximation of a local observable to its thermodynamic limit
for the ground state problem is expected to depend on whether the model is
gapped (and hence the ground state has exponentially decaying correlations
[28, 29]) or gapless. While it is generally hard to rigorously characterize
the rate of convergence of a finite-size approximation to the thermodynamic
limit for ground states, it is physically reasonable to consider the following
cases:

  * •

Logarithmic Convergence: $n=\log^{d}(1/\varepsilon)$, where $d$ is the lattice
dimension.

  * •

Power-law Convergence: $n={\rm poly}(1/\varepsilon)$.

The first case is expected to hold for gapped models. In fact, under the local
topological quantum order condition [31, 34, 33], one can easily show that
this is the case. The second case is expected to hold for critical (gapless)
models. For instance, for the Gaussian fermionic Hamiltonians analyzed in the
previous section this is the case, under very general conditions for the Fermi
surface.

Furthermore, to ensure that there is a quantum algorithm that reaches the
ground state we will assume that $H_{n}$ is adiabatically connected to a
family of Hamiltonians $H_{n}^{0}$ with efficiently preparable ground states
with a minimal gap, $\Delta_{n}$ along the adiabatic path, fulfilling

  * •

Power-law Gap: $\Delta_{n}\geq 1/{\rm poly}(n)$.

This ensures that using the adiabatic algorithm one can reach the ground state
within an error $\varepsilon$ in a time $T_{\rm Q}={\rm
poly}(n,1/\varepsilon)$, or $T_{\rm Q}={\rm poly}(1/\varepsilon)$ if framed
entirely in terms of the precision of the thermodynamic limit111Note that we
could have also considered a constant gap, in which case, at least under
certain further assumptions on the Hamiltonian, it is provably possible to
reach the ground state in $T_{\rm Q}={\rm polylog}(n,1/\varepsilon)$ [38]..
This scaling of the gap is expected to hold for critical systems as well.

We expect that, generically, classical algorithms to compute a ground state
observable would be no better than exactly diagonalizing a finite-size
Hamiltonian. Thus, under the Power-Law Gap and Logarithmic convergence
conditions, a classical computer would require time
$\exp(O(\log^{d}(1/\varepsilon)))$ and thus one would obtain a superpolynomial
quantum advantage using a quantum simulator. Instead, under the Power-Law
Convergence conditions, a classical computer would require time
$\exp(O(\text{poly}(1/\varepsilon)))$ thus yielding an exponential quantum
advantage using a quantum simulator.

###  V.2 Quantum advantage for noisy quantum simulators

As argued above, in the absence of errors, we expect a quantum simulator to
provide an advantage over classical algorithms for computing thermodynamic
limits. The presence of errors, however, sets a limit on the precision that
can be obtained by quantum simulators. In the error models analyzed in
previous sections, such a limit is determined by the hardware error $\delta$.
In time $\textnormal{poly}(1/\varepsilon^{\prime})$, the quantum simulator is
expected to compute the thermodynamic limit of the perturbed model to a
precision $\varepsilon^{\prime}$ — the precision of the target thermodynamic
limit obtained is thus upper bounded by

$\varepsilon\leq O(\text{max}(\varepsilon^{\prime},f(\delta))),$

where $f(\delta)$ is given in Definition 1. From this, we immediately see that
in the presence of hardware errors, we do not gain any benefit in precision by
running the quantum simulator beyond a time needed to obtain
$\varepsilon^{\prime}=f(\delta)$ with which we can expect to compute the
target thermodynamic limit to a precision of $O(f(\delta))$. Based on the
stability analysis of the previous subsections, we generically expect
$f(\delta)=\text{poly}(\delta)$ for most stable many-body simulation tasks,
and thus to be able to obtain the thermodynamic limit to a precision of
$O(\text{poly}(\delta))$, determined entirely by the hardware error $\delta$,
in quantum-simulation time $O(\text{poly}(1/\delta))$.

A numerical illustration of this analysis is shown in Fig. 4 — here, we use
the adiabatic quantum algorithm to find the energy density observable in the
ground state of the critical SSH model (i.e. Eq. 7 with $J=1$). Figure 4(a)
shows the convergence of the energy density observable, in the absence of
errors, to its thermodynamic limit — we see that a power-law convergence is
obtained, as physically expected for gapless models. In Fig. 4(b), we use a
system-size that yields a precision of $O(f(\delta))$, as determined by the
stability bounds on the ground-state of this model, and simulate an adiabatic
algorithm to find the ground state in the presence of hardware error. We see
that, in the presence of errors, the accuracy in the achieved precision is
fundamentally limited by the hardware precision $\delta$ — Fig. 4 shows the
run-time of the adiabatic algorithm as a function of this hardware-limited
precision. We see that this run-time scales polynomially with $1/\varepsilon$,
where $\varepsilon$ is the hardware-limited precision that is achieved by the
adiabatic algorithm.

Figure 4: Numerical study of quantum adiabatic algorithm in the presence of
error. We consider using the adiabatic algorithm to find the energy density
observable for the critical SSH model in the thermodynamic limit (TL). (a)
Convergence of the energy density to the thermodynamic limit as $n\to\infty$ —
the scaling of $\varepsilon$ with $n$ reveals a power-law scaling that is
expected for gapless models. (b) The adiabatic algorithm in the presence of
hardware errors — the quantity being plotted is the error of the noisy
adiabatic algorithm from the thermodynamic limit of the noiseless model. The
precision achieved by the adiabatic algorithm is fundamentally limited by
hardware errors. (c) The final precision ($\varepsilon$) achieved by an
adiabatic algorithm in the presence of errors as a function of the adiabatic
algorithm run-time, confirming that $T\sim\text{poly}(1/\varepsilon)$ as
expected from our analysis. Thus, on decreasing the hardware error, the error
achievable by the noisy quantum simulator decreases and the run-time of the
quantum algorithm increases at-most polynomially.

To define a notion of advantage in the presence of noise, we can now compare
the classical and quantum run-times needed to achieve this hardware-limited
precision. Assuming $f(\delta)=\text{poly}(\delta)$, it follows from the
discussion in the previous subsection that we would need classical run-times
that are either superpolynomial or exponential in $\text{poly}(1/\delta)$ to
achieve the precision that can be achieved by quantum simulators in time
$\text{poly}(1/\delta)$. If $\delta$ is decreased by a constant factor, then
the run-time of the quantum simulator will only increase at-most polynomially
with this factor, while the run-time of the classical simulator will increase
by a super-polynomial or exponential factor.

More specifically, in the case of dynamics, we can evolve the noisy quantum-
simulator for time $t=O(\delta^{-1/(d+2)})$ and, from proposition 4, we obtain
an error in local observable that $\to 0$ as $\text{poly}(\delta)$. All known
classical algorithms, based either on exact methods or cluster expansion, will
have a run-time which is $\exp(O(\text{poly}(1/\delta)))$ thus one would
obtain an exponential quantum advantage in the sense of the previous
paragraph. Similary, for the ground state, depending on whether we have
logarithmic or power-law convergence, one would obtain a superpolynomial or an
exponential quantum advantage. We summarize our expectation of noisy quantum
advantage for the quantum simulation task of dynamics and ground states in the
table 1. This table is based on both the stability results provided in the
previous two sections, as well as the scaling conjectured above for the run-
times of classical algorithms.

Problem |  Stability |  Quantum Run-Time w.r.t hardware-limited precision |  Classical Run-Time w.r.t. hardware-limited precision |  Conjectured Advantage with respect to hardware errors
---|---|---|---|---
Constant time dynamics for $t=\Theta(\delta^{-1/(d+1)})$ |  Provable with $f(\delta)=\tilde{O}(\text{poly}(\delta))$ |  $\tilde{O}(\text{poly}(1/\delta))$ |  $\exp(O(\text{poly}(1/\delta)))$ |  Exponential
Ground states of stably gapped Hamiltonians with logarithmic convergence and power-law gap |  Provable with $f(\delta)=O(\delta)$ |  $O(\text{poly}(1/\delta))$ |  $\exp(O(\log^{d}(\delta^{-1}))$ |  Superpolynomial
Ground states of gapless Hamiltonians with power-law convergence and power-law gap |  Conjectured from free-fermions with $f(\delta)=O(\text{poly}(\delta))$ |  $O(\text{poly}(1/\delta))$ |  $\exp(O(1/\delta))$ |  Exponential
Table 1: Summary of stability, hardware-limited precision, classical and
quantum run-times for thermodynamic limits in many-body problems with respect
to the hardware error $\delta$, together with their conjectured quantum
advantage in the presence of noise.

##  VI Conclusion

We have considered both the stability and quantum advantage of using near-term
analogue quantum simulators for thermodynamic limits of many-body problems in
physics. Based on both existing theoretical results in many-body literature,
and new technical results for free-fermion models, we argue that many
physically relevant many-body problems are stable to a constant rate of error
on the quantum hardware being used to solve them and thus are accessible in
near-term experiments. We also hypothesize that these algorithms have an
advantage, with respect to the obtained precision, in computing thermodynamic
limits of many-body problems. Our formulation and results provides some
evidence for near-term analogue quantum simulators being useful for solving
many-body problems.

Extending the stability results for gapless models to the case of quantum
spins, or non-Gaussian fermionic systems is an immediate open problem
suggested by our work. While previous work by Hastings [25] indicates that,
under some assumption on the density of states of the many-body model, such a
stability result could hold for gapless spin systems, it would be interesting
to see if restricting observables to being translationally invariant could
help improve these results. Furthermore, we have assumed a coherent error
model for analyzing stability of quantum simulation task, future directions
could include expanding these results to incoherent error models. Finally, any
progress on a rigorous analysis of impact of errors on adiabatic quantum
algorithms would help lay current experimental efforts on a strong theoretical
foundation.

###### Acknowledgements.

We acknowledge useful discussions from Dorit Aharonov and Álvaro M. Alhambra.
This research is funded by the German Federal Ministry of Education and
Research (BMBF) through EQUAHUMO (Grant No. 13N16066) within the funding
program quantum technologies - from basic research to market and by the Munich
Quantum Valley (MQV), which is supported by the Bavarian state government with
funds from the Hightech Agenda Bayern Plus. R.T. also acknowledges funding
from the Max-Planck Harvard Research Center for Quantum Optics (MPHQ)
postdoctoral fellowship. A.F.R. is supported by the Alexander von Humboldt
Foundation through a Postdoctoral Fellowship.

## References

  * Jozsa [2001] R. Jozsa, Quantum factoring, discrete logarithms, and the hidden subgroup problem, Computing in science & engineering 3, 34 (2001).
  * Montanaro [2016] A. Montanaro, Quantum algorithms: an overview, npj Quantum Information 2, 1 (2016).
  * Aharonov  _et al._ [1996] D. Aharonov, M. Ben-Or, R. Impagliazzo, and N. Nisan, Limitations of noisy reversible computation, arXiv preprint quant-ph/9611028 (1996).
  * Aharonov and Ben-Or [1997] D. Aharonov and M. Ben-Or, Fault-tolerant quantum computation with constant error, in  _Proceedings of the twenty-ninth annual ACM symposium on Theory of computing_ (1997) pp. 176–188.
  * Aharonov  _et al._ [2006] D. Aharonov, A. Kitaev, and J. Preskill, Fault-tolerant quantum computation with long-range correlated noise, Physical review letters 96, 050504 (2006).
  * Altman  _et al._ [2021] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksson, K.-M. C. Fu, _et al._ , Quantum simulators: Architectures and opportunities, PRX Quantum 2, 017003 (2021).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).
  * Cirac and Zoller [2012] J. I. Cirac and P. Zoller, Goals and opportunities in quantum simulation, Nature physics 8, 264 (2012).
  * Kraus  _et al._ [2007] C. Kraus, M. M. Wolf, and J. I. Cirac, Quantum simulations under translational symmetr, Phys. Rev. A 75, 022303 (2007).
  * Cubitt  _et al._ [2018] T. S. Cubitt, A. Montanaro, and S. Piddock, Universal quantum hamiltonians, Proceedings of the National Academy of Sciences 115, 9497 (2018).
  * Zhou and Aharonov [2021] L. Zhou and D. Aharonov, Strongly universal hamiltonian simulators, arXiv preprint arXiv:2102.02991 (2021).
  * Stilck França and Garcia-Patron [2021] D. Stilck França and R. Garcia-Patron, Limitations of optimization algorithms on noisy quantum devices, Nature Physics 17, 1221 (2021).
  * De Palma  _et al._ [2022] G. De Palma, M. Marvian, C. Rouzé, and D. S. França, Limitations of variational quantum algorithms: a quantum optimal transport approach, arXiv preprint arXiv:2204.03455 (2022).
  * González-García  _et al._ [2022] G. González-García, R. Trivedi, and J. I. Cirac, Error propagation in nisq devices for solving classical optimization problems, arXiv preprint arXiv:2203.15632 (2022).
  * Bouland  _et al._ [2019] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani, On the complexity and verification of quantum random circuit sampling, Nature Physics 15, 159 (2019).
  * Aaronson and Arkhipov [2011] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, in  _Proceedings of the forty-third annual ACM symposium on Theory of computing_ (2011) pp. 333–342.
  * Hamilton  _et al._ [2017] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian boson sampling, Physical review letters 119, 170501 (2017).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Zhu  _et al._ [2022] Q. Zhu, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, M. Gong, _et al._ , Quantum computational advantage via 60-qubit 24-cycle random circuit sampling, Science bulletin 67, 240 (2022).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, _et al._ , Strong quantum computational advantage using a superconducting quantum processor, Physical review letters 127, 180501 (2021).
  * Aharonov and Irani [2022] D. Aharonov and S. Irani, Hamiltonian complexity in the thermodynamic limit, in  _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_ (2022) pp. 750–763.
  * Watson and Cubitt [2022] J. D. Watson and T. S. Cubitt, Computational complexity of the ground state energy density problem, in  _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_ (2022) pp. 764–775.
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems, in  _Statistical mechanics_ (Springer, 1972) pp. 425–431.
  * Hastings [2004a] M. B. Hastings, Lieb-schultz-mattis in higher dimensions, Physical review b 69, 104431 (2004a).
  * Hastings and Wen [2005] M. B. Hastings and X.-G. Wen, Quasiadiabatic continuation of quantum states: The stability of topological ground-state degeneracy and emergent gauge invariance, Physical review b 72, 045141 (2005).
  * Bachmann  _et al._ [2012] S. Bachmann, S. Michalakis, B. Nachtergaele, and R. Sims, Automorphic equivalence within gapped phases of quantum lattice systems, Communications in Mathematical Physics 309, 835 (2012).
  * Brandão and Kastoryano [2019] F. G. Brandão and M. J. Kastoryano, Finite correlation length implies efficient preparation of quantum thermal states, Communications in Mathematical Physics 365, 1 (2019).
  * Hastings [2004b] M. B. Hastings, Locality in quantum and markov dynamics on lattices and networks, Physical review letters 93, 140402 (2004b).
  * Hastings and Koma [2006] M. B. Hastings and T. Koma, Spectral gap and exponential decay of correlations, Communications in mathematical physics 265, 781 (2006).
  * Nachtergaele and Sims [2006] B. Nachtergaele and R. Sims, Lieb-robinson bounds and the exponential clustering theorem, Communications in mathematical physics 265, 119 (2006).
  * Bravyi  _et al._ [2010] S. Bravyi, M. B. Hastings, and S. Michalakis, Topological quantum order: stability under local perturbations, Journal of mathematical physics 51, 093512 (2010).
  * Bravyi and Hastings [2011] S. Bravyi and M. B. Hastings, A short proof of stability of topological order under local perturbations, Communications in mathematical physics 307, 609 (2011).
  * Michalakis and Zwolak [2013] S. Michalakis and J. P. Zwolak, Stability of frustration-free hamiltonians, Communications in Mathematical Physics 322, 277 (2013).
  * Cirac  _et al._ [2013] J. I. Cirac, S. Michalakis, D. Pérez-García, and N. Schuch, Robustness in projected entangled pair states, Physical Review B 88, 115108 (2013).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Physical review letters 69, 2863 (1992).
  * Jordan  _et al._ [2008] J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Classical simulation of infinite-size quantum lattice systems in two spatial dimensions, Physical review letters 101, 250602 (2008).
  * Wild and Alhambra [2022] D. S. Wild and A. M. Alhambra, Classical simulation of short-time quantum dynamics, arXiv preprint arXiv:2210.11490 (2022).
  * Ge  _et al._ [2016] Y. Ge, A. Molnár, and J. I. Cirac, Rapid adiabatic preparation of injective projected entangled pair states and gibbs states, Physical review letters 116, 080503 (2016).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-robinson bounds and the generation of correlations and topological quantum order, Physical review letters 97, 050401 (2006).
  * Hastings [2007] M. B. Hastings, Quantum belief propagation: An algorithm for thermal quantum systems, Physical Review B 76, 201102 (2007).

##  Appendix A Notational preliminaries

The following is a list of notations that we will use for the mathematical
proofs in the following appendices. Here $v$ denotes a vector, and $A$ a
matrix.

  * •

$\|v\|_{p}\equiv\left(\sum_{i}{v_{i}^{p}}\right)^{\frac{1}{p}}$ denotes the
$p$-norm of the vector $v$ ($p=2$ is the Euclidean norm, for which we drop the
subscript, $\|{\cdot}\|\equiv\|{\cdot}\|_{2}$). In the $p\to\infty$ limit, it
becomes the max norm, ${\|v\|_{\infty}\equiv\max_{i}{|v_{i}|}}$.

  * •

$\|{A}\|_{\textnormal{op}}$ is the operator norm of $A$, i.e. the norm induced
by the Euclidean vector norm, or the $\infty$-Schatten norm.
$\|{A}\|_{\textnormal{op},1}$ denotes the trace norm, i.e. the 1-Schatten norm
of $A$.

  * •

$\text{vec}(A)$ denotes the vectorization of $A$, i.e. the vector whose
components are the matrix elements of $A$.

  * •

Unless otherwise mentioned, $\|{v}\|$, where $v$ is a vector, will denote
their $\ell^{2}$ norm and $\|{O}\|$, where $O$ is an operator, will be its
operator norm.

In the proofs relating to gaussian fermion models, for a hermitian operator
$O$ expressed as a quadratic form over the Majorana operators
$c_{x}^{\alpha}$,

$O=\sum_{x,y\in\mathbb{Z}^{d}_{L}}\sum_{\alpha,\beta=1}^{2D}o_{x,y}^{\alpha,\beta}c_{x}^{\alpha}c_{y}^{\beta},$

we will denote by $\tilde{O}$ the matrix of coefficients
$o_{x,y}^{\alpha,\beta}$, with the indices $(x,\alpha)$ corresponding to the
rows and $(y,\beta)$ corresponding to the columns. We will assume, without
loss of generality and unless otherwise mentioned, that $\tilde{O}$ is a
hermitian matrix with purely imaginary matrix elements.

##  Appendix B Proof of proposition 1 (Dynamics of free-fermion models)

We will need the following lemma on matrix norms:

######  Lemma 2.

Let $M$ be an $n\times n$ matrix such that (i) $|M_{ij}|\leq\delta,\forall
i,j$, and (ii) any row or column of $M$ has at most $r$ nonzero elements. Then
$\|{M}\|_{\textnormal{op}}\leq r\delta$.

###### Proof.

Let $v=(v_{i})$ be a vector, and denote $\mathcal{R}_{i}\equiv\\{j|M_{ij}\neq
0\\}$, $\mathcal{C}_{j}\equiv\\{i|M_{ij}\neq 0\\}$. Then,

$\displaystyle\|Mv\|^{2}$
$\displaystyle=\sum_{i=1}^{n}{\bigg{|}\sum_{j\in\mathcal{R}_{i}}{A_{ij}v_{j}}\bigg{|}^{2}}$
$\displaystyle\leq\sum_{i=1}^{n}{\left(\sum_{j\in\mathcal{R}_{i}}{|A_{ij}|^{2}}\sum_{k\in\mathcal{R}_{i}}{|v_{k}|^{2}}\right)}$
$\displaystyle\leq
r\delta^{2}\sum_{i=1}^{n}\sum_{k\in\mathcal{R}_{i}}{|v_{k}|^{2}}=r\delta^{2}\sum_{k=1}^{n}\sum_{i\in\mathcal{C}_{k}}{|v_{k}|^{2}}$
$\displaystyle\leq r^{2}\delta^{2}\|v\|^{2}\implies\|M\|_{\textnormal{op}}\leq
r\delta.$

∎

######  Lemma 3.

Given bounded Hermitian operators $H$ and $H^{\prime}$, for any bounded
operator $O$

$\left\|{e^{iH^{\prime}t}Oe^{-iH^{\prime}t}-e^{iHt}Oe^{-iHt}}\right\|_{\textnormal{op}}\leq
2\|{O}\|_{\textnormal{op}}\|{H-H^{\prime}}\|_{\textnormal{op}}t.$

###### Proof.

Consider the operator $\tilde{O}(t)\equiv
e^{-iHt}e^{iH^{\prime}t}Oe^{-iH^{\prime}t}e^{iHt}$. Note that

$\frac{d}{dt}\tilde{O}(t)=i\bigg{[}e^{-iHt}(H^{\prime}-H)e^{iH^{\prime}t}Oe^{-iH^{\prime}t}e^{iHt}-e^{-iHt}e^{iH^{\prime}t}Oe^{-iH^{\prime}t}(H^{\prime}-H)e^{iHt}\bigg{]},$

and consequently,

$\left\|{\frac{d}{dt}\tilde{O}(t)}\right\|\leq 2\|{H-H^{\prime}}\|\|{O}\|.$

We then immediately obtain that

$\left\|{e^{iH^{\prime}t}Oe^{-iH^{\prime}t}-e^{iHt}Oe^{-iHt}}\right\|\leq\int_{0}^{t}\left\|{\frac{d}{ds}\tilde{O}(s)}\right\|ds\leq
2\|{O}\|\|{H-H^{\prime}}\|t.$

∎

###### Proof (of proposition 1).

Let $\tilde{H},\tilde{O}$ be the $2DL^{d}\times 2DL^{d}$ matrices with
elements $h^{\alpha,\beta}_{x,y},o^{\alpha,\beta}_{x,y}$ (see Eqs. (3)-(5)),
and let $\Gamma_{0}$ be the correlation matrix in the Majorana basis
corresponding to the initial state $\rho_{0}$,

$(\Gamma_{0})^{\alpha,\beta}_{x,y}\equiv\frac{1}{2}\operatorname{Tr}{(\rho_{0}[c^{\alpha}_{x},c^{\beta}_{y}])}.$
(13)

After evolving for time $t$, the expectation values of the unperturbed and
perturbed observable are given by

$\langle
O(t)\rangle_{H}=\text{Tr}\big{(}\Gamma_{0}e^{i\tilde{H}t}\tilde{O}e^{-i\tilde{H}t}\big{)}\
\text{and }\langle
O(t)\rangle_{H^{\prime}}=\text{Tr}\big{(}\Gamma_{0}e^{i\tilde{H}^{\prime}t}\tilde{O}e^{-i\tilde{H}^{\prime}t}\big{)}.$

Using
$|\operatorname{Tr}{A}|\leq\operatorname{Tr}{|A|}=\|A\|_{\textnormal{op},1}$,
and Hölder’s inequality,
$\|AB\|_{\textnormal{op},1}\leq\|A\|_{\textnormal{op},1}\|B\|_{\textnormal{op}}$,
we have

$\displaystyle|{\langle O(t)\rangle_{H}-\langle O(t)\rangle_{H^{\prime}}}|$
$\displaystyle\leq\|{\tilde{O}}\|_{\text{op},1}\|{e^{-i\tilde{H}t}\Gamma_{0}e^{i\tilde{H}t}-e^{-i\tilde{H}^{\prime}t}\Gamma_{0}e^{i\tilde{H}^{\prime}t}}\|\leq
2t\|{\tilde{O}}\|_{\textnormal{op},1}\|{\Gamma_{0}}\|\|{\tilde{H}-\tilde{H}^{\prime}}\|,$

where the second inequality follows from lemma 3. It only remains to prove
that the last expression is independent of lattice-size $L$. Observe that
since $O$ only acts on $k$ sites, $\tilde{O}$ has at most $2kD$ nonzero
eigenvalues, thus

$\|{\tilde{O}}\|_{\textnormal{op},1}\leq{2kD}\|{\tilde{O}}\|.$

Furthermore, for any correlation matrix, $\|{\Gamma_{0}}\|\leq 1$. Finally, by
assumption, $\tilde{H}-\tilde{H}^{\prime}$ is a matrix that satisfies the
conditions of lemma 2 (its elements are bounded by $\delta$, and due to
$H,H^{\prime}$ being local, it has at most $2DR^{d}$ nonzero elements per row
and column), hence

$\|{\tilde{H}-\tilde{H}^{\prime}}\|\leq 2DR^{d}\delta$ (14)

and the proposition follows. ∎

##  Appendix C Proof of proposition 2 (Ground states of local gaussian
fermionic models)

In this appendix we will prove the stability of the expectation value of a
translationally invariant, k-locally generated Gaussian observable on the
ground state of a quadratic Hamiltonian. We first provide a lemma that uses
the translation invariance of a local observable to provide an error bound.

######  Lemma 4.

Consider a quadratic operator $O$ which is translationally invariant and
expressible as

$O=\frac{1}{n}\sum_{x\in\mathbb{Z}_{L}^{d}}\tau_{x}(O_{0}),$

where $n=L^{d}$ is the number of sites in $\mathbb{Z}_{L}^{d}$, $O_{0}$ is a
quadratic operator with a support on at most $k$ sites and $\tau_{x}$ is a
super-operator that translates an operator by $x$, then for any quadratic
operator $A_{0}$,

$\left|{\textnormal{Tr}(O^{\dagger}A_{0})}\right|\leq\frac{4D^{2}k}{n}\|{\tilde{O}_{0}}\|\|{\tilde{A}_{0}}\|_{\textnormal{op},1}.$

###### Proof.

We note that

$\textnormal{Tr}(\tilde{O}^{\dagger}\tilde{A}_{0})=\frac{1}{n}\sum_{x\in\mathbb{Z}_{L}^{d}}\text{Tr}\big{(}\tilde{O}_{0}\tau_{x}^{\dagger}(\tilde{A}_{0}))=\frac{1}{n}\sum_{x\in\mathbb{Z}^{d}}\text{Tr}\big{(}\tilde{O}_{0}\tau_{-x}(\tilde{A}_{0})\big{)}.$

Define $A\equiv n^{-1}\sum_{x\in\mathbb{Z}_{L}^{d}}\tau_{-x}(A_{0})$. We note
that $A$ is translationally invariant on the underlying lattice —
consequently, if $F$ is the $n\times n$ Fourier transform matrix, then
$F_{D}=F\otimes I_{2D}$ block diagonalizes $\tilde{A}$ i.e.
$F_{D}\tilde{A}F_{D}^{\dagger}$ will be a block diagonal matrix with $n$
blocks of size $2D\times 2D$. Then we use

$\left|{\operatorname{Tr}(\tilde{O}^{\dagger}\tilde{A}_{0})}\right|=\left|{\operatorname{Tr}(F_{D}^{\dagger}\tilde{O}^{\dagger}F^{\phantom{\dagger}}_{D}F_{D}^{\dagger}\tilde{A}F^{\phantom{\dagger}}_{D})}\right|\leq\left\|{\text{vec}\big{(}F_{D}^{\dagger}\tilde{O}_{0}F^{\phantom{\dagger}}_{D}\big{)}}\right\|_{\infty}\left\|{\text{vec}(F_{D}^{\dagger}\tilde{A}F^{\phantom{\dagger}}_{D})}\right\|_{1}$

where we have applied Hölder’s inequality to the norms of the vectorized
matrices. Now we bound each of the factors in the right hand side. Since the
operator $O_{0}$ has support only $k$ sites, $\tilde{O}_{0}$ only has
$2Dk\times 2Dk$ non-zero elements. Suppose that $\Pi_{O_{0}}$ is a diagonal
matrix with $1$s on the entries that correspond to non-zero elements of
$\tilde{O}_{0}$ — it then follows that
$\tilde{O}_{0}=\Pi_{O_{0}}\tilde{O}_{0}\Pi_{O_{0}}$. We further note that if
$f_{i}$ is the $i^{\text{th}}$ column of $F_{D}$ then

$\left\|{\text{vec}\big{(}F_{D}^{\dagger}\tilde{O}_{0}F_{D}\big{)}}\right\|_{\infty}=\sup_{i,j}{\left|{f_{i}^{\dagger}\Pi_{O_{0}}\tilde{O}_{0}\Pi_{O_{0}}f_{j}}\right|}\leq\|{\tilde{O}_{0}}\|\sup_{i,j}\|{\Pi_{O_{0}}f_{i}}\|\|{\Pi_{O_{0}}f_{j}}\|=\frac{2Dk}{n}\|{\tilde{O}_{0}}\|,$

where we have used that each entry of $F_{D}$ has magnitude $1/\sqrt{n}$ since
it is the Fourier transform matrix. Next, since
$F_{D}^{\dagger}\tilde{A}F_{D}$ is block diagonal with $N$ $2D\times 2D$
blocks, and labelling by $A_{1},A_{2}\dots A_{N}$ these blocks, we obtain that

$\left\|{\text{vec}(F_{D}^{\dagger}\tilde{A}F_{D})}\right\|_{1}=\sum_{i=1}^{N}\left\|{\text{vec}(A_{i})}\right\|_{1}\leq
2D\sum_{i=1}^{N}\left\|{A_{i}}\right\|_{\text{op},1}=2D\left\|{F_{D}^{\dagger}\tilde{A}F_{D}}\right\|_{\text{op},1}=2D\left\|{\tilde{A}}\right\|_{\text{op},1}.$

where we have used $\|{\text{vec}(M)}\|_{1}\leq n\|{M}\|_{\textnormal{op},1}$
for an $n\times n$ matrix222To see this, let
$\sigma_{ij}\equiv\text{sign}(M_{ji})$. Then
${\|{\sigma}\|_{\textnormal{op}}\leq n\|{\text{vec}(\sigma)}\|_{\infty}=n}$,
and $\|{\text{vec}(M)}\|_{1}=\operatorname{Tr}(\sigma
M)\leq\|{\sigma}\|_{\textnormal{op}}\|{M}\|_{\textnormal{op},1}\leq
n\|{M}\|_{\textnormal{op},1}$.. Finally, since
$A=\sum_{x\in\mathbb{Z}_{L}^{d}}\tau_{-x}(A_{0})/N$, it follows that
$\|{\tilde{A}}\|_{\text{op},1}\leq\|{\tilde{A}_{0}}\|_{\text{op},1}$.
Combining the above estimates, the lemma statement follows. ∎

The correlation matrix $\Gamma$ of the ground state of a quadratic Hamiltonian
$H$ with matrix of coefficients $\tilde{H}$ (see Eq. (3)) is given by

$\Gamma=\text{sign}(\tilde{H}),$

where $\text{sign}(x)=x/|{x}|$ for $x\neq 0$ and $0$ for $x=0$333The reader
may be familiar with the equivalent formulation in terms of complex fermions,
where the function to be applied to the Hamiltonian matrix to obtain the
correlation matrix of the ground state is of the Heaviside type, such that it
populates negative energy states and depopulates positive energy states. The
function of the sign function in the language of Majorana fermions is exactly
analogous.. The sign function applied on a matrix is to be understood as an
operator function i.e. as a function acting on the eigenvalues of the argument
while keeping the eigenvectors unchanged. Our proof will rely on a Fourier
series approximation to the sign function. Within the interval $(-\pi,\pi)$,
we will investigate the approximation of $\text{sign}(x)$ with
$\text{sign}_{M}(x)$, where

$\text{sign}_{M}(x)=\sum_{n=-M}^{M}c_{n}e^{inx}\ \text{where}\
c_{n}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\text{sign}(x)e^{-inx}dx.$

To analyze the error between $\text{sign}_{M}(x)$ and $\text{sign}(x)$, it is
convenient to express $\text{sign}_{M}(x)$ in terms of the Dirichlet kernel,

$\text{sign}_{M}(x)\equiv\int_{-\pi}^{\pi}D_{M}(x-y)\text{sign}(y)dy,$

where

$D_{M}(x)\equiv\frac{1}{2\pi}\sum_{n=-M}^{M}e^{-inx}=\frac{1}{2\pi}\frac{\sin[(M+1/2)x]}{\sin(x/2)}.$

Below, we provide two technical lemmas about the $\text{sign}_{M}$ function —
one that quantifies the approximation error between it and the exact sign
function, and the next that quantifies the maximum value of the
$\text{sign}_{M}$ function. Both of these lemmas will be used for the
perturbation theory analysis of the free-fermion ground state problem.

Figure 5: (left) Truncated Fourier series approximation
$\operatorname{sign}_{M}(x)$ to the $\operatorname{sign}(x)$ function, used in
the proof of proposition 2. (right) Truncated Fourier series approximation
$t_{M}(x)$ to the $\tanh(\beta x)$ function (for $\beta=1$), used in the proof
of proposition 3.

######  Lemma 5.

For all $\eta\leq|{x}|\leq\pi-\eta$ and $M>0$,

$\big{|}\textnormal{sign}(x)-\textnormal{sign}_{M}(x)\big{|}\leq\frac{1}{M}+\frac{1}{M\eta}.$

###### Proof.

We first consider $x\in[\eta,\pi]$. We note that

$\displaystyle\text{sign}_{M}(x)=\int_{0}^{\pi}D_{M}(x-y)dy-\int_{0}^{\pi}D_{M}(x+y)dy.$

Now, since $\int_{-\pi}^{\pi}D_{M}(y)dy=1$, we obtain that

$\int_{0}^{\pi}D_{M}(x-y)dy=1-\int_{0}^{\pi}D_{M}(x+y)dy,$

and thus

$|{\text{sign}_{M}(x)-\text{sign}(x)}|=2\left|{\int_{0}^{\pi}D_{M}(x+y)dy}\right|.$

Next, we apply integration by parts to obtain

$\displaystyle\int_{0}^{\pi}D_{M}(x+y)dy=\frac{1}{\pi(2M+1)}\bigg{(}\frac{\cos((M+1/2)(\pi+x))}{\cos(x/2)}+\frac{\cos((M+1/2)(\pi+x))}{\sin(x/2)}-$
$\displaystyle\frac{1}{2}\int_{0}^{\pi}\frac{\cos((M+1/2)(x+y))\cos((x+y)/2)}{\sin^{2}((x+y)/2)}dy\bigg{)},$

and therefore

$\displaystyle\left|{\int_{0}^{\pi}D_{M}(x+y)dy}\right|\leq\frac{1}{\pi(2M+1)}\bigg{(}\frac{1}{|{\cos(x/2)}|}+\frac{1}{|{\sin(x/2)}|}+\frac{1}{2}\int_{0}^{\pi}\frac{|{\cos((x+y)/2)}|dy}{\sin^{2}((x+y)/2)}\bigg{)}$
$\displaystyle\qquad\qquad\leq\frac{2}{\pi(2M+1)}\bigg{(}\frac{1}{|{\cos(x/2)}|}+\frac{1}{|{\sin(x/2)}|}-1\bigg{)},$

where in the last step we have used the integral

$\frac{1}{2}\int_{x}^{\pi+x}\frac{|{\cos(y/2)}|}{{\sin^{2}(y/2)}}dy=\frac{1}{\sin(x/2)}+\frac{1}{\cos(x/2)}-2.$

Now, for $x\in(\eta,\pi/2)$, $|{\cos(x/2)}|\geq 1/\sqrt{2}$ and
$|{\sin(x/2)}|\geq x/\pi\geq\eta/\pi$. Therefore, we obtain that

$\left|{\int_{0}^{\pi}D_{M}(x+y)dy}\right|\leq\frac{2}{\pi(2M+1)}\bigg{(}\sqrt{2}+\frac{\pi}{\eta}-1\bigg{)}.$

While this bound is true for $x\in[\eta,\pi/2]$, we note that both
$\text{sign},\text{sign}_{M}$ satisfy $f(x)=f(\pi-x)$ for $x\in[0,\pi]$ and
consequently this bound also holds for $x\in[\pi/2,\pi-\eta]$. Finally, since
for both $\text{sign},\text{sign}_{M}$, $f(x)=-f(-x)$, it follows that this
bound holds for $[-\pi+\eta,-\eta]\cup[\eta,\pi-\eta]$. ∎

######  Lemma 6.

For all $x\in[-\pi,\pi]$,

$|{\textnormal{sign}_{M}(x)}|\leq 5.$

###### Proof.

This proof is an adaptation of the standard technique based on Riemann
integration that is used to treat Gibbs phenomena in Fourier analysis. We
repurpose that technique to provide error bounds as a function of $M$ instead
of just concentrating on the asymptotic limit $M\to\infty$. Again, we only
consider $x\in[0,\pi/2]$, and extend the bound on $|{\text{sign}_{M}(x)}|$ to
the remaining interval by symmetry. We divide the interval $[0,\pi/2]$ into
$[0,\alpha_{0}/M]\cup[\alpha_{0}/M,\pi/2]$, where $\alpha_{0}$ is a constant
that we pick later.

Consider first $x\in[\alpha_{0}/M,\pi/2]$. An application of lemma 5 yields

$|{\text{sign}_{M}(x)}|\leq 1+\frac{2}{\pi(2M+1)}\bigg{(}\sqrt{2}-1+\frac{\pi
M}{\alpha_{0}}\bigg{)}.$

For large $M$, this bound scales as $\sim 1/\alpha_{0}$ and thus does not
allow us to provide an upper bound on $\text{sign}_{M}(x)$ for $x$ close to
$0$. For this, we use the representation of $\text{sign}_{M}(x)$ as a Fourier
series which approximates a Riemann integral of $\sin(t)/t$. Consider
$x\in[0,\alpha_{0}/M)$ and let $\alpha=xM$ ($\alpha\leq\alpha_{0}$). Note that

$\text{sign}_{M}(x)=\frac{2}{\pi}\sum_{k\in[1,M]|k\text{ is
odd}}\frac{2}{k}\sin\bigg{(}\frac{k\alpha}{M}\bigg{)}$

To bound the term in the summation, we observe that it is an approximation of
the Riemann integral of $\sin(\alpha x)/x$ in the interval $[0,1]$. In
particular, since $\sup_{x\in\mathbb{R}}|{(\sin x/x)^{\prime}}|\leq 2$,
Taylor’s theorem yields that

$\left|{\sum_{k\in[1,M]|k\text{ is
odd}}\frac{2}{k}\sin\bigg{(}\frac{k\alpha}{M}\bigg{)}-\int_{0}^{1}\frac{\sin\alpha
x}{x}dx}\right|\leq\frac{4\alpha^{2}}{M}\leq\frac{4\alpha_{0}^{2}}{M}.$

Finally, we note that

$\int_{0}^{1}\frac{\sin\alpha x}{x}dx\leq\alpha\leq\alpha_{0}.$

Thus, we obtain that for $x\in[0,\alpha_{0}/M)$,

$|{\text{sign}_{M}(x)}|\leq\alpha_{0}+\frac{4\alpha_{0}^{2}}{M}.$

Thus, for the entire interval $[0,\pi/2]$, we obtain that

$|{\text{sign}_{M}(x)}|\leq\text{max}\bigg{(}\alpha_{0}+\frac{4\alpha_{0}^{2}}{M},1+\frac{2}{\pi(2M+1)}\bigg{(}\sqrt{2}-1+\frac{\pi
M}{\alpha_{0}}\bigg{)}\bigg{)}.$

Since this holds for any $\alpha_{0}$, we choose $\alpha_{0}=1$. We then
obtain that

$|{\text{sign}_{M}(x)}|\leq\text{max}\bigg{(}1+\frac{4}{M},1+\frac{2}{\pi(2M+1)}\big{(}\sqrt{2}-1+{\pi
M}\big{)}\bigg{)}\leq 5\text{ for }M\geq 1.$

∎

###### Proof (of proposition 2).

The expectation value of the observable $O$ in the ground state of the
Hamiltonian $H$ is given by

$\langle O\rangle_{H}=\text{Tr}\big{(}\tilde{O}\
\text{sign}(\tilde{H})\big{)},\ \langle
O\rangle_{H^{\prime}}=\text{Tr}\big{(}\tilde{O}\
\text{sign}(\tilde{H^{\prime}})\big{)}.$

Without loss of generality, we will assume that $\tilde{H},\tilde{H}^{\prime}$
are normalized so that
$\|{\tilde{H}}\|,\|{\tilde{H}^{\prime}}\|\leq\frac{\pi}{2}$. This way all the
eigenfrequencies lie in the interval $[-\frac{\pi}{2},\frac{\pi}{2}]$. Note
that lemma 2 guarantees that this can be done with a constant normalization
factor, i.e. one that does not depend on the system size, and does not change
the ground state (note however that $\delta$ and $f_{h}$ would have to be
rescaled accordingly). Now, from lemma 4, it follows that

$\left|{\langle O\rangle_{H}-\langle
O\rangle_{H^{\prime}}}\right|\leq\frac{4D^{2}k}{n}\|{\tilde{O}_{0}}\|\|{\text{sign}(\tilde{H})-\text{sign}(\tilde{H}^{\prime})}\|_{\text{op},1}.$

Furthermore,

$\displaystyle\|{\text{sign}(\tilde{H})-\text{sign}(\tilde{H}^{\prime})}\|_{\text{op},1}\leq$
$\displaystyle\quad\quad\|{\text{sign}(\tilde{H})-\text{sign}_{M}(\tilde{H})}\|_{\text{op},1}+\|{\text{sign}(\tilde{H}^{\prime})-\text{sign}_{M}(\tilde{H}^{\prime})}\|_{\text{op},1}+\|{\text{sign}_{M}(\tilde{H})-\text{sign}_{M}(\tilde{H^{\prime}})}\|_{\text{op},1}.$

We bound each term on the right hand side separately. Consider
$\|{\text{sign}(\tilde{H})-\text{sign}_{M}(\tilde{H})}\|_{\text{op},1}$ —
denoting by $\lambda_{i}$ the eigenvalues of $\tilde{H}$ and for any $\eta>0$,
we can express it as

$\|{\text{sign}(\tilde{H})-\text{sign}_{M}(\tilde{H})}\|_{\text{op},1}=\sum_{i|\lambda_{i}\in[-\eta,\eta]}|{\text{sign}(\lambda_{i})-\text{sign}_{M}(\lambda_{i})}|+\sum_{i|\lambda_{i}\notin[-\eta,\eta]}|{\text{sign}(\lambda_{i})-\text{sign}_{M}(\lambda_{i})}|.$

The motivation behind splitting the error into these two terms is that, within
the interval $[-\frac{\pi}{2},\frac{\pi}{2}]$, the approximation of
$\operatorname{sign}{(\lambda)}$ by $\operatorname{sign}_{M}{(\lambda)}$ is
only good outside the neighbourhood of $0$ (see Fig. 5) — consequently, we
treat the eigenvalues of $\tilde{H}$ which lie within $\eta$ radius of $0$
separately from the rest. It now follows that from assumption 1 and lemma 6
that

$\sum_{i|\lambda_{i}\in[-\eta,\eta]}|{\text{sign}(\lambda_{i})-\text{sign}_{M}(\lambda_{i})}|\leq
6nf_{h}(\eta)+6\kappa(\eta,n).$

Furthermore, from lemma 5,

$\sum_{i|\lambda_{i}\notin[-\eta,\eta]}|{\text{sign}(\lambda_{i})-\text{sign}_{M}(\lambda_{i})}|\leq\frac{n}{M}\bigg{(}1+\frac{1}{\eta}\bigg{)}.$

Therefore, we obtain that

$\frac{1}{n}\|{\text{sign}(\tilde{H})-\text{sign}_{M}(\tilde{H})}\|_{\text{op},1}\leq
6f_{h}(\eta)+6\frac{\kappa(\eta,n)}{n}+\frac{1}{M}\bigg{(}1+\frac{1}{\eta}\bigg{)}.$

We can similarly analyze
$\|{\text{sign}(\tilde{H}^{\prime})-\text{sign}(\tilde{H}^{\prime})}\|_{\text{op},1}$.
Denote by $\lambda^{\prime}_{i}$ the eigenvalues of $\tilde{H}^{\prime}$ — it
follows from Weyl’s theorem that
$|{\lambda_{i}-\lambda_{i}^{\prime}}|\leq\|{\tilde{H}-\tilde{H^{\prime}}}\|_{\text{op}}\leq
2DR^{d}\delta$ (see Eq. (14)). Consequently, for sufficiently small, but
$\Theta(1)$, $\delta$, we obtain that

$\sum_{i|\lambda_{i}^{\prime}\in[-\eta,\eta]}|{\text{sign}(\lambda_{i}^{\prime})-\text{sign}_{M}(\lambda_{i}^{\prime})}|\leq
6nf_{h}(\eta+2DR^{d}\delta)+6\kappa(\eta+2DR^{d}\delta,n),$

and

$\sum_{i|\lambda_{i}^{\prime}\notin[-\eta,\eta]}|{\text{sign}(\lambda_{i}^{\prime})-\text{sign}_{M}(\lambda_{i}^{\prime})}|\leq\frac{n}{M}\bigg{(}1+\frac{1}{\eta}\bigg{)}.$

Finally, we consider
$\|{\text{sign}_{M}(H)-\text{sign}_{M}(H^{\prime})}\|_{\text{op},1}\leq
n\|{\text{sign}_{M}(H)-\text{sign}_{M}(H^{\prime})}\|$. Now, denoting by
$\\{c_{m}\\}_{m\in\mathbb{Z}}$ the Fourier series components of sign function,
then

$\|{\text{sign}_{M}(H)-\text{sign}_{M}(H^{\prime})}\|\leq\sum_{m=-M}^{M}|{c_{m}}|\|{e^{im\tilde{H}}-e^{im\tilde{H}^{\prime}}}\|\leq\sum_{m=-M}^{M}|{mc_{m}}|\|{\tilde{H}-\tilde{H}^{\prime}}\|.$

Using the explicit expression for $c_{m}$, we can immediately conclude that
$|{mc_{m}}|=2/\pi$ when $m$ is odd, and $0$ when $m$ is even. Therefore, we
obtain that

$\|{\text{sign}_{M}(H)-\text{sign}_{M}(H^{\prime})}\|\leq\frac{2(M+1)}{\pi}\|{\tilde{H}-\tilde{H}^{\prime}}\|\leq\frac{2(M+1)}{\pi}2DR^{d}\delta.$

Combining all of these estimates, we obtain that

$\frac{1}{n}\|{\text{sign}(\tilde{H})-\text{sign}(\tilde{H}^{\prime})}\|_{\text{op},1}\leq
cM\delta+\frac{2}{M}\bigg{(}1+\frac{1}{\eta}\bigg{)}+6\big{(}f_{h}(\eta)+f_{h}(\eta+c^{\prime}\delta)\big{)}+6\bigg{(}\frac{\kappa(\eta,n)}{n}+\frac{\kappa(\eta+c^{\prime}\delta,n)}{n}\bigg{)}.$

with $c,c^{\prime}$ constants. Since this is valid for any $\eta$ and $M$,
choosing $M=\delta^{-1/2}$ and $\eta=\delta^{1/4}$, we obtain the proposition.
∎

##  Appendix D Proof of proposition 3 (Gibbs state of free-fermion models)

The correlation matrix of a thermal state of a quadratic Hamiltonian can be
written in terms of the coefficient matrix $H$ of the latter as
$\Gamma=\tanh(\beta H)$. Note that the $\beta\to\infty$ limit yields the sign
function, which was used in the previous appendix to compute the ground state
correlation matrix. Indeed, the reasoning here will be similar to that of
appendix C, replacing the sign function with the hyperbolic tangent. The next
couple of lemmas discuss the Fourier series approximation of $\tanh{\beta x}$,
defined as

${t}_{M}(x)\equiv\sum_{n=-M}^{M}c_{n}e^{inx},\ \text{where
}c_{n}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\tanh{\beta x}e^{-inx}dx.$

######  Lemma 7.

For $M\geq 1$, and $x\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right]$,

$|{t_{M}(x)-\tanh{\beta x}}|\leq\frac{q(\beta)}{M},$

where $q(\beta)\equiv
12\pi^{2}\beta^{3}+2\pi^{2}\beta^{2}+\left(2+\frac{\pi^{2}}{2}\right)\beta+\left(\frac{4\sqrt{2}}{\pi}+\frac{\pi^{2}}{2}\right)=O(\beta^{3})$.

###### Proof.

We fix the value of $\beta$ and let $t(x)$ be the $2\pi$-periodic extension of
$\tanh{\beta x}$

$t(x)\equiv\tanh{\beta(x-2n\pi)},\qquad x-2n\pi\in[-\pi,\pi],\quad
n\in\mathbb{Z}$ (15)

Once again, it will be convenient to represent $t_{M}(x)$ in terms of the
Dirichlet kernel $D_{M}$. We note,

$t_{M}(x)=\int_{-\pi}^{\pi}{D_{M}(x-y)t(y)dy}=\int_{-\pi}^{\pi}{D_{M}(y)t(x-y)dy}=\int_{-\pi}^{\pi}{D_{M}(y)t(x+y)dy},$

and therefore, using that the Dirichlet kernel is normalized, we write

$t(x)-t_{M}(x)=\dfrac{1}{2}\int_{-\pi}^{\pi}{D_{M}(y)\left(2t(x)-t(x-y)-t(x+y)\right)dy}=\int_{0}^{\pi}{D_{M}(y)f_{x}(y)dy},$

where in the last step we have defined $f_{x}(y)\equiv 2t(x)-t(x-y)-t(x+y)$.
In the integration interval $[0,\pi]$, $f_{x}(y)$ is piecewise smooth with a
single jump discontinuity at $y=\pi-x$. We thus split the integral into the
two intervals $[0,\pi-x]$ and $[\pi-x,\pi]$ and apply integration by parts in
each of them. For the first one,

$\displaystyle\int_{0}^{\pi-x}{D_{M}(y)f_{x}(y)dy}=-\dfrac{1}{\pi}\dfrac{\cos{\left(\left(M+\frac{1}{2}\right)y\right)}}{2M+1}\dfrac{f_{x}(y)}{\sin\frac{y}{2}}\Bigg{|}_{y=0}^{\pi-x}+\dfrac{1}{(2M+1)\pi}\int_{0}^{\pi-x}{g_{x}(y)\dfrac{\cos{\left(\left(M+\frac{1}{2}\right)y\right)}}{\sin^{2}\frac{y}{2}}dy}$

where $g_{x}(y)\equiv
2\sin\frac{y}{2}f_{x}^{\prime}(y)-\cos{\frac{y}{2}}f_{x}(y)$. To bound this
expression, we will use the following properties of the functions
$f_{x}(y),g_{x}(y)$ on the interval $[0,\pi-x]$, where they are smooth:

$\displaystyle f_{x}(0)=f^{\prime}_{x}(0)=0,|f_{x}(y)|\leq
4,|f^{\prime}_{x}(y)|\leq 2\beta,|f^{\prime\prime}_{x}(y)|\leq
4\beta^{2},|f^{\prime\prime\prime}_{x}(y)|\leq 12\beta^{3},$ $\displaystyle
g_{x}(0)=g^{\prime}_{x}(0)=0,|g_{x}^{\prime\prime}(y)|\leq
24\beta^{3}+4\beta^{2}+\beta+1.$

These bounds follow from direct computation, and in the case of $g_{x}(y)$
they are easiest to see when expressed in terms of $f_{x}(y)$. They imply (via
Taylor’s theorem with second order remainder) that

$|g_{x}(y)|\leq(24\beta^{3}+4\beta^{2}+\beta+1)\dfrac{y^{2}}{2}$

which together with $\sin^{2}(y)\geq\frac{y^{2}}{\pi^{2}}$ will allow us to
bound the integral. Putting it all together, we have

$\left|{\int_{0}^{\pi-x}{D_{M}(y)f_{x}(y)dy}}\right|\leq\dfrac{4\sqrt{2}}{(2M+1)\pi}+\dfrac{\pi^{2}}{(2M+1)}(24\beta^{3}+4\beta^{2}+\beta+1)$

Now we proceed on to the second interval $y\in[\pi-x,\pi]$ and similarly
integrate by parts,

$\displaystyle\int_{\pi-x}^{\pi}{D_{M}(y)f_{x}(y)dy}=-\dfrac{1}{\pi}\dfrac{\cos{\left(\left(M+\frac{1}{2}\right)y\right)}}{2M+1}\dfrac{f_{x}(y)}{\sin\frac{y}{2}}\Bigg{|}_{y=\pi-x}^{\pi}+\dfrac{1}{(2M+1)\pi}\int_{\pi-x}^{\pi}{g_{x}(y)\dfrac{\cos{\left(\left(M+\frac{1}{2}\right)y\right)}}{\sin^{2}\frac{y}{2}}dy}.$

Now the bound on $g_{x}(y)$ from Taylor’s theorem no longer holds, due to the
discontinuity, but since $y=0$ is not in the integration interval, we can just
use the constant bound $|g(x)|\leq 4\beta+4$ to obtain

$\left|{\int_{\pi-x}^{\pi}{D_{M}(y)f_{x}(y)dy}}\right|\leq\dfrac{4\sqrt{2}}{(2M+1)\pi}+\dfrac{4}{(2M+1)}(\beta+1),$

and putting everything together the lemma follows. ∎

######  Lemma 8.

If $\\{c_{n}\\}_{n\in\mathbb{Z}}$ are the Fourier series coefficients of
$\tanh{\beta x}$ in the interval $[-\pi,\pi]$, then for $M\geq 1$

$\sum_{n=-M}^{M}|{nc_{n}}|\leq 2M(\beta+1).$

###### Proof.

This follows by a straightforward manipulation of $c_{n}$ — note that
$c_{0}=0$, and for $n\neq 0$, we obtain from integration by parts that

$c_{n}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\tanh{\beta x}\
e^{-inx}dx=\frac{1}{2\pi}\bigg{(}\frac{2i}{n}\tanh\beta\pi\
e^{-in\pi}+\frac{\beta}{in}\int_{-\pi}^{\pi}\frac{e^{-inx}}{\text{cosh}^{2}\beta
x}dx\bigg{)}.$

Consequently,

$|{c_{n}}|\leq\frac{1}{2\pi}\bigg{(}\frac{2}{n}+\frac{2\pi\beta}{n}\bigg{)}\leq\frac{\beta+1}{n}.$

From this bound, the lemma follows. ∎

###### Proof (of proposition 3).

We bound the error between $\langle O\rangle_{H,\beta}$ and $\langle
O\rangle_{H^{\prime},\beta}$ using the same procedure as for the ground state
(see appendix C) — the proof simplifies significantly because $\tanh{\beta x}$
does not have a discontinuity near $x=0$ (unlike the sign function). From
lemma 4 it follows that

$|{\langle O\rangle_{H,\beta}-\langle
O\rangle_{H^{\prime},\beta}}|\leq\frac{4D^{2}k}{n}\|{\tilde{O}_{0}}\|\|{\tanh{\beta\tilde{H}}-\tanh{\beta\tilde{H}^{\prime}}}\|_{\text{op},1}.$

We again split

$\displaystyle\|{\tanh{\beta\tilde{H}}-\tanh{\beta\tilde{H}^{\prime}}}\|_{\text{op},1}\leq$
$\displaystyle\qquad\|{\tanh{\beta\tilde{H}}-t_{M}(\tilde{H})}\|_{\text{op},1}+\|{\tanh{\beta\tilde{H}^{\prime}}-t_{M}(\tilde{H}^{\prime})}\|_{\text{op},1}+\|{t_{M}(\tilde{H}^{\prime})-t_{M}(\tilde{H})}\|_{\text{op},1}$

We will assume once again that $\|{H}\|,\|{H^{\prime}}\|\leq\frac{\pi}{2}$, so
that from lemma 7, it follows that

$\|{\tanh{\beta\tilde{H}}-t_{M}(\tilde{H})}\|_{\text{op},1},\|{\tanh{\beta\tilde{H}^{\prime}}-t_{M}(\tilde{H}^{\prime})}\|_{\text{op},1}\leq\frac{nq(\beta)}{M},$

and

$\|{t_{M}(\tilde{H})-t_{M}(\tilde{H}^{\prime})}\|_{\text{op},1}\leq
n\|{t_{M}(\tilde{H})-t_{M}(\tilde{H}^{\prime})}\|_{\text{op}}\leq
n\sum_{m=-M}^{M}|{c_{m}}|\|{e^{im\tilde{H}}-e^{im\tilde{H}^{\prime}}}\|_{\text{op}}.$

Furthermore, from lemmas 3 and 2 we have
$\|{e^{im\tilde{H}}-e^{im\tilde{H}^{\prime}}}\|_{\text{op}}\leq mc\delta$,
where $c=2DR^{d}$. Thus, from lemma 8, it follows that

$\|{t_{M}(\tilde{H})-t_{M}(\tilde{H}^{\prime})}\|_{\text{op},1}\leq
2nM(\beta+1)c\delta.$

Thus, we obtain that for any $M>1$,

$|{\langle O\rangle_{H,\beta}-\langle O\rangle_{H^{\prime},\beta}}|\leq
4D^{2}k\|{O_{0}}\|_{\text{op}}\bigg{(}\frac{2q(\beta)}{M}+2(\beta+1)cM\delta\bigg{)}.$

choosing $M=\sqrt{q(\beta)/c(\beta+1)\delta}$, we obtain the result. ∎

##  Appendix E Proof of proposition 4 (Dynamics of locally interacting spin
systems)

We will need the following two lemmas:

######  Lemma 9 (Lieb-Robinson bounds, Ref. [24, 39]).

Given a Hamiltonian defined on a lattice $\mathbb{Z}^{d}_{L}$,
$H=\sum_{\alpha\in\mathbb{Z}_{L}^{d}s}h_{\alpha}$ such that
$\|{h_{\alpha}}\|\leq J$ and $h_{\alpha}$ acts on sites $\alpha^{\prime}$ such
that $d(\alpha,\alpha^{\prime})\leq R$, then there exist positive constants
$\mu,v$ that depend only on the lattice such that

  1. (a)

For any two operators $A,B$ with support $S_{A},S_{B}$ and $l=d(S_{A},S_{B})$,

$\|{[A(t),B]}\|\leq|{S_{A}}|\|{A}\|\|{B}\|e^{-\mu
l}\big{(}e^{vJ|{t}|}-1\big{)},$

where $A(t)=e^{iHt}Ae^{-iHt}$.

  2. (b)

For any local operator $O$ with support $S_{O}$, and for $l>0$,

$\|{O(t)-O_{l}(t)}\|\leq\|{O}\|\left|{S_{O}}\right|e^{-\mu
l}\big{(}e^{vJt}-1\big{)},$

where $O_{l}(t)=e^{iH_{l}t}Oe^{-iH_{l}t}$ with
$H_{l}=H-\sum_{\alpha|d(S_{h^{\alpha}},S_{O})\geq l}h_{\alpha}$ being the
restriction of the Hamiltonian to a region within distance $l$ of $S_{O}$.

###### Proof (of proposition 4).

We consider the Hamiltonians $H$ and $H^{\prime}$ represented as

$H=\sum_{\alpha\in\mathbb{Z}_{L}^{d}}h_{\alpha},\qquad
H^{\prime}=\sum_{\alpha\in\mathbb{Z}_{L}^{d}}h_{\alpha}^{\prime}.$

We have, by assumption,
$\|{h_{\alpha}-h_{\alpha}^{\prime}}\|_{\text{op}}\leq\varepsilon$. Thus if
$\varepsilon<J$, then for all $\alpha\in\mathbb{Z}_{L}^{d}$,
$\|{h_{\alpha}}\|,\|{h_{\alpha}^{\prime}}\|\leq 2J$. Now,

$|{\langle O\rangle_{H,t}-\langle
O\rangle_{H^{\prime},t}}|\leq\sum_{K\in\\{H,H^{\prime}\\}}\left\|{e^{iKt}Oe^{-iKt}-e^{iK_{l}t}Oe^{-iK_{l}t}}\right\|+\left\|{e^{iH_{l}t}Oe^{-iH_{l}t}-e^{iH_{l}^{\prime}t}Oe^{-iH_{l}^{\prime}t}}\right\|.$

We note that by simply counting the number of sites in the support of the
truncated Hamiltonian,

$\|{H_{l}-H_{l}^{\prime}}\|\leq(R_{O}+2R+2l)^{d}\varepsilon,$

where $R_{O}=\text{diam}(S_{O})$. From lemma 3, it then follows that

$\left\|{e^{iH_{l}t}Oe^{-iH_{l}t}-e^{iH_{l}^{\prime}t}Oe^{-iH_{l}^{\prime}t}}\right\|\leq
2\varepsilon t\|{O}\|(R_{O}+2R+2l)^{d}.$

Using this together with the Lieb-Robinson bounds (lemma 9), we obtain that

$|{\langle O\rangle_{H,t}-\langle O\rangle_{H^{\prime},t}}|\leq
2\|{O}\||{S_{O}}|e^{-\mu l}\big{(}e^{2vJt}-1\big{)}+2\varepsilon
t\|{O}\|(R_{O}+2R+2l)^{d}.$

Choosing $l=2vJt/\mu+\log(2J/\varepsilon)/\mu$, we obtain that

$|{\langle O\rangle_{H,t}-\langle
O\rangle_{H^{\prime},t}}|\leq\|{O}\|\bigg{[}\frac{\varepsilon}{J}|{S_{O}}|\big{(}1-e^{-2vJt}\big{)}+2\varepsilon
t\bigg{(}R_{O}+2R+\frac{4vJt}{\mu}+\frac{2}{\mu}\log\bigg{(}\frac{2J}{\varepsilon}\bigg{)}\bigg{)}^{d}\bigg{]},$

which proves the proposition statement. ∎

##  Appendix F Proof of proposition 5 (Ground states of gapped local
Hamiltonians)

We will apply the formalism developed in Ref. for spectral flows for families
of gapped Hamiltonians. We are interested in a target spatially local
Hamiltonian $H$, expressed as

$H=\sum_{x\in\mathfrak{L}}h_{x},$

where $h_{x}$ acts only on spins with a distance $R$ of $x\in\mathfrak{L}$,
and $\|{h_{x}}\|\leq 1$. The implemented Hamiltonian $H^{\prime}$ is assumed
to have a similar form,

$H^{\prime}=\sum_{x\in\mathfrak{L}}\big{(}h_{x}+v_{x}\big{)},$

where $\|{v_{x}}\|\leq\varepsilon$ for all $x\in\mathfrak{L}$. We assume that
$H$ is stably gapped with gap $\Delta$ i.e. any $H^{\prime}$ of the above form
has an energy gap between the ground state and the first excited state that is
larger than $\Delta$. We consider the family of Hamiltonians, $H_{s}$, for
$s\in[0,1]$, defined by

$H_{s}=H+s(H^{\prime}-H)=\sum_{x\in\mathfrak{L}}h_{x}+sv_{x},$

and note that the assumption of being stably gapped is equivalent to $H_{s}$
being gapped, with the gap being larger than $\Delta$, for all $s\in[0,1]$.
Now, the spectral flow method allows us to construct a unitary $U(s)$ that
relates the ground state $\ket{G_{s=0}}$ of $H_{s=0}=H$ to the ground state
$\ket{G_{s}}$ of $H_{s}$ as provided in the following lemma.

######  Lemma 10 (From Ref. [26]).

Consider the unitary $U(s)$ obtained from

$\frac{d}{ds}U(s)=iD(s)U(s)\ \text{where
}D(s)=\int_{-\infty}^{\infty}W_{\Delta}(t)e^{-itH_{s}}(H^{\prime}-H)e^{itH_{s}}dt,$

where $W_{\Delta}\in L^{1}(\mathbb{R})$ is a real valued odd function which
satisfies

  1. (a)

$|{W_{\Delta}(t)}|$ is bounded and satisfies

$\displaystyle\|{W_{\Delta}}\|_{\infty}=\sup_{t\in\mathbb{R}}|{W_{\Delta}(t)}|=\frac{1}{2}.$
(16)

  2. (b)

For $t>0$, the function $I_{\Delta}(t)=\int_{t}^{\infty}|{W_{\Delta}(s)}|ds$
satisfies

$\displaystyle I_{\Delta}(t)\leq G(\Delta t),$ (17)

where $G(x)$ falls off faster than any polynomial as $x\to\infty$.

Then, $\ket{G_{s}}=U(s)\ket{G_{s=0}}$, where $\ket{G_{s}}$ is the ground state
of $H(s)$.

_Proof (of proposition 5)_. Using this result, we can straightforwardly show
the stability of the quantum simulation task of computing a local observable
in the ground state of $H$. To see this, we note that

$|{\bra{G_{0}}O\ket{G_{0}}-\bra{G_{s}}O\ket{G_{s}}}|=\left|{\bra{G_{0}}\bigg{(}O-U^{\dagger}(s)OU(s)\bigg{)}\ket{G_{0}}}\right|\leq\|{O-U^{\dagger}(s)OU(s)}\|\leq\int_{0}^{s}\|{[O,D(s^{\prime})]}\|ds^{\prime}.$

It then remains to bound $\|{[O,D(s^{\prime})]}\|$ — we can do this by
following lemma 4.7 in Ref. [26], and we reproduce this below — we start by
noting that

$\|{[A,D(s^{\prime})]}\|\leq\sum_{x\in\mathfrak{L}}\left\|{\int_{-\infty}^{\infty}W_{\gamma}(t)[A,e^{itH_{s}}v_{x},e^{-itH_{s}}]dt}\right\|.$

For each term in this summation, we further split the integral and bound it as

$\left\|{\int_{-\infty}^{\infty}W_{\Delta}(t)[A,e^{itH_{s}}v_{x},e^{-itH_{s}}]dt}\right\|\leq\left\|{\int_{|{t}|\leq
T_{x}}W_{\Delta}(t)[A,e^{itH_{s}}v_{x}e^{-itH_{s}}]dt}\right\|+\left\|{\int_{|{t}|>T_{x}}W_{\Delta}(t)[A,e^{itH_{s}}v_{x}e^{-itH_{s}}]dt}\right\|.$

For the first term, which only concerns with $|{t}|\leq T_{x}$, we use the
Lieb Robinson’s bound (lemma 9) and Eq. 16 to obtain

$\left\|{\int_{|{t}|\leq
T_{x}}W_{\Delta}(t)[O,e^{itH_{s}}v_{x}e^{-itH_{s}}]dt}\right\|\leq\|{O}\|\|{v_{x}}\||{S_{O}}|e^{-\mu
d(S_{O},S_{v_{x}})}\int_{0}^{T_{x}}(e^{vt}-1)dt\leq\|{O}\|\|{v_{x}}\||{S_{O}}|\frac{e^{-\mu
d(S_{O},S_{v_{x}})}e^{vT_{x}}}{v}.$

For the second term for $|{t}|\geq T_{x}$, we use Eq. 17 together with the
fact that $W_{\Delta}$ is an odd function and the simple bound
$\|{[O,e^{itH_{s}}v_{x}e^{-itH_{s}}]}\|\leq 2\|{O}\|\|{v_{x}}\|$ to obtain
that

$\left\|{\int_{|{t}|\geq
T_{x}}W_{\Delta}(t)[O,e^{itH_{s}}v_{x}e^{-itH_{s}}]dt}\right\|\leq
2\|{O}\|\|{v_{x}}\|\int_{|{t}|\geq T_{x}}|{W_{\Delta}(t)}|\leq
2\|{O}\|\|{v_{x}}\|G(\Delta T_{x}),$

Note that $T_{x}$ can be arbitrary in the above two estimates — choosing
$T_{x}=\mu d(S_{O},S_{v_{x}})/2v$, we obtain that

$\left\|{\int_{-\infty}^{\infty}W_{\Delta}(t)[O,e^{itH_{s}}v_{x},e^{-itH_{s}}]dt}\right\|\leq\|{O}\|\|{v_{x}}\|\bigg{[}\frac{|{S_{O}}|}{v}e^{-\mu
d(S_{O},S_{v_{x}})/2}+2G\bigg{(}\frac{\Delta\mu}{2v}d(S_{O},S_{v_{x}})\bigg{)}\bigg{]},$

and therefore, for all $s^{\prime}\in[0,s]$, we obtain that bound

$\|{[O,D(s^{\prime})]}\|\leq\|{O}\|\varepsilon\sum_{x\in\mathfrak{L}}\bigg{[}\frac{|{S_{O}}|}{v}e^{-\mu
d(S_{O},S_{v_{x}})/2}+2G\bigg{(}\frac{\Delta\mu}{2v}d(S_{O},S_{v_{x}})\bigg{)}\bigg{]}.$

Noting that the summand in the above expression decreases faster than any
polynomial in $d(S_{O},S_{v_{x}})$, we see that it will be upper bounded by a
constant independent of the size of the lattice $\mathfrak{L}$, thus
independent of $n$. This proves the proposition. $\hfill\square$

##  Appendix G Proof of proposition 6 (Gibbs state with exponential clustering
of correlations)

We begin by presenting a proof of a standard bound on the perturbation of
Gibb’s states of a Hamiltonian. This can be found in , and we reproduce it
here for the convenience of the reader.

######  Lemma 11.

Suppose $H$ and $V$ are two Hermitian bounded operators, then for any
$\beta\geq 0$:

  1. (a)

The partition functions satisfy

$\textnormal{Tr}(e^{-\beta(H+V)})\leq\textnormal{Tr}(e^{-\beta
H})e^{\beta\|{V}\|}$

  2. (b)

For any $O\succeq 0$, it follows that

$\textnormal{Tr}\big{(}Oe^{-\beta(H+V)}\big{)}\leq\textnormal{Tr}\big{(}Oe^{-\beta
H}\big{)}\exp\big{(}e^{\beta(\|{H}\|+\|{V}\|)}\|{V}\|\big{)}.$

###### Proof.

(a) We will use the Duhamel’s formula, which states that for any
differentiable bounded operator $F(t)$,

$\frac{d}{dt}e^{F(t)}=\int_{0}^{1}e^{(1-u)F(t)}\frac{dF(t)}{dt}e^{uF(t)}dt.$

Defining $H(s)=H+sV$, we note from the Duhamel’s formula that

$\left|{\frac{d}{ds}\text{Tr}\big{(}e^{-\beta
H(s)}\big{)}}\right|=\beta\left|{\int_{0}^{1}\text{Tr}\big{(}e^{-(1-u)\beta
H(s)}Ve^{-u\beta
H(s)}ds\big{)}du}\right|=\beta\left|{\int_{0}^{1}\text{Tr}\big{(}e^{-\beta
H(s)}V\big{)}du}\right|\leq\beta\|{V}\|_{\text{op}}\|{e^{-\beta
H(s)}}\|_{\text{op},1},$

where we have used the Holder’s inequality in the last step. Noting that
$\|{e^{-\beta H(s)}}\|_{\text{op},1}=\text{Tr}(e^{-\beta H(s)})$, we obtain
that

$\left|{\frac{d}{ds}\text{Tr}\big{(}e^{-\beta
H(s)}\big{)}}\right|\leq\beta\|{V}\|_{\text{op}}\text{Tr}\big{(}e^{-\beta
H(s)}\big{)}$

Therefore,

$|{\log\text{Tr}(e^{-\beta(H+V)})-\log\text{Tr}(e^{-\beta
H})}|\leq\int_{0}^{1}\left|{\frac{d}{ds}\log\text{Tr}(e^{-\beta
H(s)})}\right|ds=\int_{0}^{1}\left|{\frac{1}{\text{Tr}(e^{-\beta
H(s)})}\frac{d}{ds}\text{Tr}(e^{-\beta
H(s)})}\right|ds\leq\beta\|{V}\|_{\text{op}}.$

Thus, we obtain the lemma statement.

(b) We again use the Duhamel’s formula to obtain

$\displaystyle\left|{\frac{d}{ds}\text{Tr}(Oe^{-\beta
H(s)})}\right|=\beta\left|{\int_{0}^{1}\text{Tr}\big{(}Oe^{-(1-u)\beta
H(s)}Ve^{-u\beta
H(s)}\big{)}du}\right|\leq\beta\int_{0}^{1}\left|{\text{Tr}\big{(}Oe^{-(1-u)\beta
H(s)}Ve^{-u\beta H(s)}\big{)}}\right|du$

We note that

$\displaystyle\left|{\text{Tr}\big{(}Oe^{-(1-u)\beta H(s)}Ve^{-u\beta
H(s)}\big{)}}\right|$ $\displaystyle=\left|{\text{Tr}\big{(}e^{-\beta
H(s)/2}Oe^{-\beta H(s)/2}e^{-(1/2-u)\beta H(s)}Ve^{-(u-1/2)\beta
H(s)}\big{)}}\right|,$ $\displaystyle\leq\text{Tr}\big{(}e^{-\beta
H(s)/2}Oe^{-\beta H(s)/2}\big{)}\|{e^{-(1/2-u)\beta H(s)}Ve^{-(u-1/2)\beta
H(s)}}\|,$ $\displaystyle\leq\text{Tr}\big{(}Oe^{-\beta
H(s)}\big{)}\|{V}\|e^{\beta(\|{H}\|+s\|{V}\|)}.$

where we have used the fact that since $u\in[0,1]$, $|{u-1/2}|\leq 1/2$.
Therefore, we obtain that

$\left|{\frac{d}{ds}\log\textnormal{Tr}\big{(}Oe^{-\beta
H(s)}\big{)}}\right|\leq\|{V}\|e^{\beta\|{H}\|}e^{\beta
s\|{V}\|}\implies\left|{\log\text{Tr}(Oe^{-\beta
H})-\log\text{Tr}(Oe^{-\beta(H+V)})}\right|\leq\|{V}\|e^{\beta(\|{H}\|+\|{V}\|)}.$

This estimate yields the lemma statement. ∎

######  Lemma 12.

Given bounded hermitian operators $H$ and $V$, and any bounded hermitian
operator (observable) O,

$\left|{\textnormal{Tr}\bigg{(}\frac{Oe^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}-\textnormal{Tr}\bigg{(}\frac{Oe^{-\beta(H+V)}}{Z_{H+V}(\beta)}\bigg{)}}\right|\leq
2\|{O}\|_{\textnormal{op}}\bigg{(}\left|{\exp(e^{\beta(\|{H}\|+\|{V}\|)}\|{V}\|)-1}\right|+\bigg{|}{\exp(\beta\|{V}\|)-1}\bigg{|}\bigg{)}$

###### Proof.

This is a straightforward application of lemma 11. We denote by
$H^{\prime}=H+V$. For simplicity, we will analyze the operator
$O^{\prime}=O+\|{O}\|I$ — $O^{\prime}\succeq 0$ and

$\text{Tr}\bigg{(}\frac{O^{\prime}e^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}-\text{Tr}\bigg{(}\frac{O^{\prime}e^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}=\text{Tr}\bigg{(}\frac{Oe^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}-\text{Tr}\bigg{(}\frac{Oe^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}.$

We begin by noting that

$\left|{\textnormal{Tr}\bigg{(}\frac{O^{\prime}e^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}-\textnormal{Tr}\bigg{(}\frac{O^{\prime}e^{-\beta
H^{\prime}}}{Z_{H^{\prime}}(\beta)}\bigg{)}}\right|\leq\frac{\textnormal{Tr}(O^{\prime}e^{-\beta
H})}{Z_{H}(\beta)}\left|{1-\frac{\textnormal{Tr}(O^{\prime}e^{-\beta
H^{\prime}})}{\textnormal{Tr}(O^{\prime}e^{-\beta
H})}}\right|+\frac{\text{Tr}(O^{\prime}e^{-\beta
H^{\prime}})}{Z_{H^{\prime}}(\beta)}\left|{1-\frac{Z_{H^{\prime}}(\beta)}{Z_{H}(\beta)}}\right|$

Noting that

$\frac{\textnormal{Tr}(O^{\prime}e^{-\beta
H})}{Z_{H}(\beta)},\frac{\textnormal{Tr}(O^{\prime}e^{-\beta
H^{\prime}})}{Z_{H^{\prime}}(\beta)}\leq\|{O^{\prime}}\|_{\text{op}}\leq
2\|{O}\|_{\text{op}},$

we obtain that

$\left|{\textnormal{Tr}\bigg{(}\frac{O^{\prime}e^{-\beta
H}}{Z_{H}(\beta)}\bigg{)}-\textnormal{Tr}\bigg{(}\frac{O^{\prime}e^{-\beta
H^{\prime}}}{Z_{H^{\prime}}(\beta)}\bigg{)}}\right|\leq
2\|{O}\|_{\text{op}}\bigg{(}\left|{\exp(e^{\beta(\|{H}\|+\|{V}\|)}\|{V}\|)-1}\right|+\bigg{|}{\exp(\beta\|{V}\|)-1}\bigg{|}\bigg{)}$

∎

We next need the notion of exponentially-clustered correlations in a Gibb’s
state — which we reproduce below from Ref. [27]. We will consider Hamiltonians
on $\mathfrak{L}\subset\mathbb{Z}^{d}$ expressed as

$H=\sum_{x\in\mathfrak{L}}h_{x},$

where $h_{x}$ acts only on spins within a distance $R$ of $x\in\mathfrak{L}$.
We will denote by $\text{supp}(h_{x})\subseteq\mathfrak{L}$ the support of
$h_{x}$. Given $X\subseteq\mathfrak{L}$, we denote by $H_{X}$ the operator

$H_{X}=\sum_{x|\text{supp}(h_{x})\subseteq X}h_{x},$

i.e. $H_{X}$ is the Hamiltonian $H$ obtained on restricting $H$ to the set
$X$.

######  Definition 2.

A local Hamiltonian $H$ is said to have exponential clustering of correlations
at inverse-temperature $\beta$ if $\exists c_{1},c_{2}>0$ such that for all
operators $X\subset\mathfrak{L}$ and $A,B$ with
$\textnormal{supp}(A),\textnormal{supp}(B)\subset X$ with
$d(\textnormal{supp}(A),\textnormal{supp}(B))\geq l$,

$\left|{\textnormal{Tr}\big{(}A\otimes
B\sigma_{X}(\beta)\big{)}-\textnormal{Tr}\big{(}A\sigma_{X}(\beta)\big{)}\textnormal{Tr}\big{(}B\sigma_{X}(\beta)\big{)}}\right|\leq
c_{2}\|{A}\|\|{B}\|e^{-c_{1}l},$

where $\sigma_{X}(\beta)=e^{-\beta H_{X}}/\textnormal{Tr}[e^{-\beta H_{X}}]$
is the Gibb’s state corresponding to $H_{X}$ at inverse-temperature $\beta$.

An important property of Hamiltonians with exponential clustering of
correlations, which relies on quantum belief propagation [40] and is proved in
Ref. [27], is that local observables can be estimated locally.

######  Lemma 13 (From Ref. [27]).

Suppose $H$ is a local Hamiltonian on a finite lattice
$\mathfrak{L}\subset\mathbb{Z}^{d}$ with exponential clustering of
correlations at inverse temperature $\beta$. If
$\mathfrak{L}=\mathcal{A}\cup\mathcal{B}\cup\mathcal{C}$ such that
$\text{dist}(\mathcal{A},\mathcal{C})\geq l$, then $\exists
c_{1}^{\prime},c_{2}^{\prime}$ such that

$\left\|{\textnormal{Tr}_{\mathcal{B},\mathcal{C}}(\sigma_{\mathfrak{L}}(\beta))-\textnormal{Tr}_{\mathcal{B}}(\sigma_{\mathcal{A}\cup\mathcal{B}}(\beta))}\right\|_{\textnormal{tr}}\leq|{\partial\mathcal{C}}|c_{2}^{\prime}e^{-c_{1}^{\prime}l},$

where $\sigma_{X}(\beta)$ is the Gibb’s state corresponding to $H_{X}$ and
$\partial\mathcal{C}$ is the boundary between $B,C$.

_Proof (of proposition 6)_. We assume that both $H$ and $H^{\prime}$ have
exponential clustering of correlations and satisfy lemma 13. Suppose $O$ is a
local observable with support $S_{O}$ and consider $\mathcal{B}$ to be a
region around $S_{O}$ and $\mathcal{C}$ be the remainder of the lattice. We
also assume that $d(\mathcal{C},S_{O})\geq l$, for some $l$ to be chosen
later. We denote by $\sigma_{l}(\beta)$ and $\sigma_{l}^{\prime}(\beta)$ the
Gibb’s state, at inverse temperature $\beta$, corresponding to
$H_{S_{O}\cup\mathcal{B}}$ and $H^{\prime}_{S_{O}\cup\mathcal{B}}$
respectively, and by $\sigma(\beta),\sigma^{\prime}(\beta)$ the Gibb’s state
corresponding to $H$ and $H^{\prime}$. Now, from lemma 13 it follows that

$|{\text{Tr}(O\sigma(\beta))-\text{Tr}(O\sigma_{l}(\beta))}|,|{\text{Tr}(O\sigma^{\prime}(\beta))-\text{Tr}(O\sigma_{l}^{\prime}(\beta))}|\leq\|{O}\|d(2l+R_{O})^{d-1}c_{2}^{\prime}e^{-c_{1}^{\prime}l},$

where $R_{O}=\text{diam}(S_{O})$ and we have used that
$|{\partial\mathcal{C}}|\leq
d\times\textnormal{diam}(S_{O}\cup\mathcal{B})^{d-1}\leq d(2l+R_{O})^{d-1}$.
Furthermore, lemma 12 can be used to bound
$|{\text{Tr}(\sigma_{l}(\beta))-\text{Tr}(O\sigma_{l}^{\prime}(\beta))}|$. We
note that

$\|{H_{S_{O}\cup\mathcal{B}}-H_{S_{O}\cup\mathcal{B}}^{\prime}}\|\leq\delta(2l+R_{O})^{d}\text{
and }\|{H_{S_{O}\cup\mathcal{B}}}\|\leq(2l+R_{O})^{d}.$

Therefore,

$|{\text{Tr}(\sigma_{l}(\beta))-\text{Tr}(O\sigma_{l}^{\prime}(\beta))}|\leq
2\|{O}\|O(e^{\beta(2l+R_{O})^{d}}(2l+R_{O})^{d}\delta).$

Thus, from the triangle inequality we obtain the bound that

$|{\text{Tr}(\sigma(\beta))-\text{Tr}(O\sigma^{\prime}(\beta))}|\leq\|{O}\|\bigg{[}O\bigg{(}(2l+R_{O})^{d-1}e^{-c_{1}^{\prime}l}\bigg{)}+O\bigg{(}e^{\beta(2l+R_{O})^{d}}(2l+R_{O})^{d}\delta\bigg{)}\bigg{]}.$

Choosing $2l+R_{O}=\Theta(\log^{1/d}(1/\sqrt{\delta}))$, we obtain that

$\left|{\text{Tr}(O\sigma(\beta))-\text{Tr}(O\sigma^{\prime}(\beta))}\right|\leq\|{O}\|O(\log^{1-1/d}(1/\delta)e^{-\Omega(\log^{1/d}(1/\delta))}\big{)},$

which proves the lemma statement. $\hfill\square$