# Quantum many-body systems in thermal equilibrium

Álvaro M. Alhambra<EMAIL_ADDRESS>Max-Planck-Institut für
Quantenoptik, Hans-Kopfermann-Strasse 1, D-85748 Garching, Germany  Instituto
de Física Teórica UAM/CSIC, C/ Nicolás Cabrera 13-15, Cantoblanco, 28049
Madrid, Spain

###### Abstract

The thermal or equilibrium ensemble is one of the most ubiquitous states of
matter. For models comprised of many locally interacting quantum particles, it
describes a wide range of physical situations, relevant to condensed matter
physics, high energy physics, quantum chemistry and quantum computing, among
others. We give a pedagogical overview of some of the most important universal
features about the physics and complexity of these states, which have the
locality of the Hamiltonian at its core. We focus on mathematically rigorous
statements, many of them inspired by ideas and tools from quantum information
theory. These include bounds on their correlations, the form of the
subsystems, various statistical properties, and the performance of classical
and quantum algorithms. We also include a summary of a few of the most
important technical tools, as well as some self-contained proofs.

###### Contents

  1. I Introduction
    1. I.1 Scope and content
    2. I.2 Completely analytical interactions
  2. II Mathematical preliminaries and notation
    1. II.1 Operator norms
    2. II.2 Information-theoretic quantities
    3. II.3 Lattice notation
    4. II.4 Asymptotic notation
  3. III An overview of technical tools
    1. III.1 Connected cluster expansion
    2. III.2 Thermal locality estimates
    3. III.3 Quantum belief propagation
    4. III.4 Selected trace inequalities
  4. IV Correlations
    1. IV.1 Correlations between neighbouring regions: Thermal area law
    2. IV.2 Decay of long-range correlations
    3. IV.3 A refined correlation decay: Conditional mutual information
  5. V Locality of temperature
    1. V.1 Local indistinguishability
      1. V.1.1 Proof in 1D
    2. V.2 Hamiltonian of mean force
    3. V.3 Non-equilibrium thermodynamics with strong coupling
  6. VI Statistical properties
    1. VI.1 Measurement statistics and concentration bounds
      1. VI.1.1 Chernoff-Hoeffding bound
      2. VI.1.2 Large deviation bound
      3. VI.1.3 Berry-Esseen theorem
    2. VI.2 Equivalence of ensembles
  7. VII Algorithms and complexity of thermal states
    1. VII.1 An efficient classical algorithm for the 1D partition function
    2. VII.2 Hardness of approximating partition functions
    3. VII.3 Tensor network methods
    4. VII.4 Quantum algorithms for preparing thermal states
  8. VIII Commuting Hamiltonians
  9. IX Conclusions and open questions
    1. IX.1 List of open questions
  10. A Derivations of Gibbs states
    1. A.1 Gibbs weights from the ultra-weak coupling assumption
    2. A.2 Jaynes’ maximum entropy principle
  11. B Miscellaneous proofs

##  I Introduction

We are currently at the dawn of the age of synthetic quantum matter.
Increasingly better experiments on a variety of quantum platforms are
improving in size and controllability at unprecedented rates, aided by the
current impulse of quantum information science and technology. This gives very
good prospects to the exploration of the physics of complex quantum many body
systems. Our aspiration to better understand these systems is very well
motivated from a scientific perspective, but also potentially from the
industrial one: unlocking the potential of complex quantum systems may bring
surprising advances to the engineering of new materials or chemical compounds
in the future. It may also yield computational tools with unprecedented
capabilities for a still unknown range of applications.

Many of the most commonly studied materials and current experimental platforms
are described by an arrangement of quantum particles in some sort of lattice
configuration. Due to the spatial decay of electromagnetic forces, each of
these particles only interacts appreciably with their immediate vicinity,
which causes the interactions between them to be local.

In this tutorial, we focus on the properties of these important systems when
at thermal equilibrium, so that they are accurately described by the so-called
thermal or Gibbs state. We review and explain some of their most important
universal properties, covered from a mathematical perspective. That is, we
focus on statements that can be proven about states of the Gibbs form

$\rho_{\beta}=\frac{e^{-\beta H}}{Z}=\frac{1}{Z}\sum_{l}e^{-\beta
E_{l}}|E_{l}\rangle\langle E_{l}|,$ (1)

where $H=\sum_{l}E_{l}|E_{l}\rangle\langle E_{l}|$ is the Hamiltonian, $\beta$
is the inverse temperature and $Z\equiv\text{Tr}[e^{-\beta H}]$ is the
partition function. The Hamiltonian describes the interactions between the $N$
particles, which are restricted to short-ranged or local. A “local
Hamiltonian” is a Hermitian operator $H$ in the finite-dimensional Hilbert
space of $N$ $d$-dimensional particles $(\mathbb{C}^{d})^{\otimes N}$. It is
defined as a sum of terms

$H=\sum_{i}h_{i}\otimes\mathbb{I},$ (2)

each of which has support (i.e. acts non-trivially) on at most $k$ particles,
and bounded strength, such that

$\max_{i}||h_{i}||=h$ (3)

For a definition of the operator norm $||\cdot||$ see Sec. II.1 below.
Typically, the Hamiltonians are scaled so that $h=\mathcal{O}(1)$.

In what follows we just write the terms as $h_{i}$ for simplicity. These
constitute the individual interactions, which are typically arranged in a
lattice of a small dimension. A simple example is e.g. the transverse-field
Ising model in one dimension with open boundary conditions

$H_{\text{Ising}}=\sum_{j=1}^{N-1}\left(J\sigma^{X}_{j}\sigma^{X}_{j+1}+\Delta\sigma^{Z}_{j}\right)+\Delta\sigma^{Z}_{N}.$
(4)

Here, $k=2$ and the interactions are arranged on a 1D chain.

The idea of a local Hamiltonian is very general, and involves many different
models describing a wide range of situations, of interest for many fields of
physics, chemistry and computer science. The only thing they have in common is
the _locality_ of the interactions. We aim to understand mathematically how
this fact alone constrains both the physics and the computational complexity
when combined with thermal fluctuations.

The thermal states of these general local Hamiltonians appear in many
different contexts, and are interesting for a wide variety of reasons. Some of
the main ones are:

  * •

It is one of the most ubiquitous states of quantum matter: typical experiments
happen at finite temperature, where the quantum system at hand is weakly
coupled to some external radiation field or phonon bath, that drive it to the
thermal state. For completeness, we sketch the standard argument of how the
weak coupling assumption leads to states of the Gibbs form in App. A.1.

  * •

The thermal state is also important when studying not just systems with an
external bath, but also in the evolution of isolated quantum systems, even
when their global state is pure: in very generic cases, these end up being
“their own bath”, and the individual subsystems thermalize to the Gibbs
ensemble Gogolin and Eisert (2016); D’Alessio  _et al._ (2016).

  * •

From a general condensed matter/material science standpoint, we are very
interested in numerous questions about the physics at finite temperature: How
are conserved quantities (e.g. charge, energy) propagated in a state close to
equilibrium? How does the system respond to small or large perturbations away
from equilibrium?

  * •

Systems at thermal equilibrium (both quantum and classical) display
interesting phase transitions in certain (low) temperature regimes (e.g.
classical Ising model in 2D). It is thus relevant to study what are their
universal properties both in and away from the critical points.

  * •

They are also important from the point of view of quantum phases of matter and
topological order at zero temperature. It has been widely established that
thermal states of local models in dimension $D-1$ appear in the entanglement
spectrum of $D$\- dimensional ground states Li and Haldane (2008). As such,
understanding their structure should also help us in elucidating the low
energy behaviour of many interesting systems.

  * •

They very naturally appear in information theory and inference as the
distributions that best reproduce partial current knowledge of a system. This
is justified by Jaynes’ principle, which we explain in App. A.2.

  * •

These states are also important for computation. For instance, being able to
sample from the thermal distribution of local models is a typical subroutine
for certain classical and quantum algorithms Brandão and Svore (2017); Brandão
_et al._ (2019); van Apeldoorn  _et al._ (2020); G.S L. Brandão  _et al._
(2022). They are also a very naturally occurring data structure in both
classical and quantum machine learning Amin  _et al._ (2018); Kieferová and
Wiebe (2017); Biamonte  _et al._ (2017); Bairey  _et al._ (2019); Anshu  _et
al._ (2021a) (often under the name of _Boltzmann machines_).

  * •

It is known from quantum computational complexity that the low energy subspace
of local Hamiltonians is able to encode the solution to very hard
computational problems: finding the ground state energy is QMA complete Kempe
_et al._ (2005). Thus, it is widely believed that even a quantum computer
should not be able to do it in polynomial time. This then at least also
applies to the thermal state at very low temperature, and motivates the study
of how the complexity changes as the temperature rises Aharonov  _et al._
(2013); Hastings (2013).

There are many different specific aspects that one could explore, but here we
focus on the following, which we believe to be of particular importance:

  * •

The correlations between the particles at different parts of the lattice. In
particular, how does the structure of those correlations are structured in
relation to the geometry of the lattice.

  * •

The states of the subsystems that a thermal state can take, and how are they
related to few-particle Gibbs states.

  * •

The statistical physics properties of these systems at equilibrium, including
Jaynes’ principle, concentration bounds and equivalence of ensembles.

  * •

The efficiency of classical and quantum algorithms for the generation and
manipulation of thermal states, and the computation of expectation values and
partition functions.

More specifically, we focus on these topics for a broad family of Gibbs states
that can be understood as being away from phase transitions within the phase
space. It is for this region of the parameter space of Hamiltonians that the
mathematical results described here are typically more tractable and give
insightful results. We describe this further in Sec. I.2.

The general topic of this tutorial, and the particular results explained here,
are a small part of the exciting past, present and future efforts to
understand the physics and complexity of quantum many-body systems. We hope to
contribute to the understanding and cross-fertilization of the many different
angles that the quantum many-body problem can take. See also e.g. Kliesch and
Riera (2018); Gogolin and Eisert (2016) for previous references with partially
overlapping content.

###  I.1 Scope and content

Throughout this tutorial, we cover statements that have a precise mathematical
formulation, many of them motivated by a quantum information theoretic
perspective. This notably includes a short exposition of a few key technical
tools in Sec. III. These have not previously appeared together, but rather
separately explained in the literature with various levels of detail,
depending on the context and usage. We hope that this encourages new,
potentially unexpected, applications thereof.

Along the sections with the actual physical and computational results, we
write the proofs of some of the simpler or more important ones explained
throughout. This includes at least one main result per section, which should
serve as a pedagogical example. For the rest, some of which have more detailed
or involved derivations, we refer the reader to the original works cited along
the text. One of our main hopes is that after reading this tutorial even the
more technical works will be more easily accessible to a wider range of
researchers. Because of this, rather than the traditional Theorem-Proof
structure of most mathematical physics writing, we have chosen a more
streamlined style for the presentation which allows for more physical
explanations and intuitions of the steps. This will hopefully contribute to a
wider readability.

Most of what we describe are known results, with at most some small
improvements, with the proofs either being the same or slightly simplified
versions of previous ones. The relevant references are included, but this does
not mean that all of the previous relevant ones are listed here: we are
certainly missing to mention a very large body of work. This includes many
relevant papers on mathematical physics, but also a lot of important physics
literature that covers these topics from perspectives that are beyond our
scope: based on numerical methods, theory work on experimental
implementations, as well as all experimental results.

###  I.2 Completely analytical interactions

Before we proceed, we should put the content of this tutorial into context
more precisely. In the mathematical study of many body systems, one of the
biggest points of interest are phase transitions, including the study of order
parameters, symmetry breaking, and other very well established ideas which aim
at classifying the possible kinds of phase transitions. For instance, there
are important models, such as the paradigmatic Ising in two and three
dimensions, that have a very well understood phase transition at a given
temperature.

However, for many classes of models and most regions of their parameter space,
local Hamiltonians do not have a thermal phase transition. This “one phase
region” (as it is sometimes referred to Lenci and Rey-Bellet (2005)) likely
contains the “simplest” cases of thermal equilibrium, including those of non-
interacting gases. These situations are typically characterized by the
analiticity of the partition function and other closely related facts, such
as:

  * •

The convergence of the cluster expansion.

  * •

The localization of correlations in the lattice, and absence of long range
order.

  * •

The approximation of marginals with local Gibbs states, and the idea of
locality of temperature.

  * •

Concentration properties of local observables.

  * •

Efficiency of approximation, either with quantum or classical algorithms.

  * •

The existence and boundedness of log-Sobolev constants.

It is expected that many or all of these simplifying facts are equivalent, in
that a model that obeys one (such as the analiticity of the partition
function) will also display the other features. In the classical case, a large
number of conditions are known to be equivalent to the analiticity of the
partition function. The study of this problem was initiated by the seminal
work of Dobrushing and Shlosman Dobrushin and Shlosman (1987), aiming at
characterizing these “completely analytical interactions” in terms of many
diferent equivalent conditions (12 in the original article). In the quantum
case, much less is known about the equivalence of the analiticity of the
partition function with other physical facts, although some important steps
have been taken (see e.g. Harrow  _et al._ (2020); Capel  _et al._ (2020)).

The main aim of this tutorial is to cover results that show that Gibbs states
have simplifying features with respect to generic quantum states. Following
that spirit, most (although not all) of the results explained here apply to
this “phase” or universality class of Gibbs states, in which those simplifying
facts are expected to hold. In fact, every element of the list above is
individually considered in each of the sections below. Because of that, we do
not cover an important part of the literature where analytical results are
typically much harder to obtain. For instance, those studying the effect of
phase transitions in e.g. the simulability of Gibbs states, the types of
correlations that can arise, and others.

##  II Mathematical preliminaries and notation

###  II.1 Operator norms

A basic but very important mathematical tool in this context are the Schatten
$p$-norms for operators, as well as the different inequalities between them.
These norms are maps from the space of operators to $\mathbb{R}$, as
$M\rightarrow||M||_{p}$, that obey the following properties:

  * •

Homogeneous: If $\alpha$ is a scalar, $||\alpha M||_{p}=|\alpha|||M||_{p}$.

  * •

Positive: $||M||_{p}\geq 0$.

  * •

Definite: $||M||_{p}=0\iff M=0$.

  * •

Triangle inequality: $||M_{1}+M_{2}||_{p}\leq||M_{1}||_{p}+||M_{2}||_{p}$.

For a given operator $M$ with singular values $\\{\alpha_{l}^{M}\\}$ and
$p\in[1,\infty)$, they are defined as

$||M||_{p}\equiv\text{Tr}[|M|^{p}]^{\frac{1}{p}}=\left(\sum_{l}(\alpha_{l}^{M})^{p}\right)^{\frac{1}{p}}.$
(5)

The more important ones are the operator norm
$||M||\equiv||M||_{\infty}=\max_{l}|\alpha_{l}^{M}|$, the Hilbert-Schmidt
$2$-norm $||M||_{2}=\text{Tr}[MM^{\dagger}]^{1/2}$ and the 1-norm or trace
norm

$||M||_{1}=\max_{||P||\leq 1}\text{Tr}[MP].$ (6)

Thus $|\text{Tr}[M]|\leq||M||_{1}$, with equality for positive operators. For
quantum states, $\text{Tr}[\rho]=||\rho||_{1}=1$.

Typically we measure the “strength” of an observable with the operator norm,
and the closeness of two quantum states with the trace norm
$||\rho-\sigma||_{1}$, since it is related to the probability of
distinguishing them under measurements. The $2$-norm, on the other hand, is
often the easiest one to compute in practice. Also note the very important
Hölder’s inequality

$||M_{1}M_{2}||_{p}\leq||M_{1}||_{q_{1}}||M_{2}||_{q_{2}},$ (7)

which holds for $\frac{1}{p}=\frac{1}{q_{1}}+\frac{1}{q_{2}}$ (e.g.
$p=q_{1}=1,q_{2}=\infty$). A particularly useful corollary is the Cauchy-
Schwarz inequality, when $q_{1}=q_{2}=2$ and $p=1$,

$|\text{Tr}[M_{1}^{\dagger}M_{2}]|^{2}\leq\text{Tr}[M_{1}M_{1}^{\dagger}]\text{Tr}[M_{2}M_{2}^{\dagger}].$
(8)

###  II.2 Information-theoretic quantities

Let us define the von Neumann entropy of a quantum state $\rho$ 111The $\log$s
here are defined with base $e$.

$S(\rho)=-\text{Tr}[\rho\log(\rho)],$ (9)

which, roughly speaking, quantifies the uncertainty we have about the
particular state. It is bounded by $0\leq S(\rho)\leq\log d$. The lower bound
is obtained by choosing $\rho$ pure, and the upper bound by the identity
$\rho=\mathbb{I}/d$. Another important quantity is the Umegaki relative
entropy

$D(\rho|\sigma)=\text{Tr}[\rho(\log\rho-\log\sigma)],$ (10)

which is a measure of distinguishability of quantum states. It obeys Pinsker’s
inequality

$D(\rho|\sigma)\geq\frac{1}{2}||\rho-\sigma||_{1}^{2},$ (11)

which relates the relative entropy with the 1-norm. It is strictly positive
for $\rho\neq\sigma$, and vanishes otherwise. It is also closely related to
the non-equilibrium free energy

$D(\rho|\rho_{\beta})=\beta\text{Tr}[\rho H]-S(\rho)+\log Z\equiv\beta
F_{\beta}(\rho)+\log Z,$ (12)

which also shows that the equilibrium free energy is
$F_{\beta}(\rho_{\beta})=-\beta^{-1}\log Z$. This distance measure also
naturally appears in the derivation of Jayne’s principle, as shown in App.
A.2.

From these quantities we can also define the quantum mutual information,
which, given a bipartite state $\rho^{AB}$ on subsystems $A$ and $B$, with
$\text{Tr}_{B}[\rho^{AB}]=\rho_{A},\text{Tr}_{A}[\rho^{AB}]=\rho^{B}$,
quantifies the correlations between $A$ and $B$ as

$\displaystyle I(A:B)_{\rho}$
$\displaystyle=S(\rho^{A})+S(\rho^{B})-S(\rho^{AB})$ (13)
$\displaystyle=D(\rho^{AB}|\rho^{A}\otimes\rho^{B}).$ (14)

In particular, it is zero if and only if $\rho^{AB}=\rho^{A}\otimes\rho^{B}$.
For all these three functions we can also define their corresponding Rényi
generalizations. See Khatri and Wilde (2020); Scalet  _et al._ (2021) for
details.

A further, perhaps more refined quantity is the _conditional mutual
information_ (CMI), defined as

$\displaystyle I(A:C|B)_{\rho}=$ $\displaystyle
S(\rho^{AB})+S(\rho^{BC})-S(\rho^{ABC})-S(\rho^{B})$ $\displaystyle=$
$\displaystyle I(A:BC)_{\rho}-I(A:B)_{\rho}.$ (15)

This perhaps less known quantity is behind many non-trivial statements in
quantum information theory (see Section 11.7 in Wilde (2017) for more
details). In a nutshell, it measures how much $A$ and $C$ share correlations
that are _not_ mediated by $B$. That is, if this quantity is small, most of
the correlations between $A$ and $C$ (which may be weak) are in reality
correlations between $A$ and $B$ and $B$ and $C$.

###  II.3 Lattice notation

In what follows we need some technical definitions regarding the properties of
the Hamiltonian and the lattice. The lattice is a hypergraph which we denote
by $\Lambda=\\{V,E\\}$ with vertex set $V$ and hyperedges $E$. To each vertex
we associate a Hilbert space of dimension $d$, $\mathbb{C}^{d}$. The number of
particles is $N=|V|$, and the number of hyperedges is $|E|$. The locality of
the Hamiltonian can be expressed by a parameter $\mathfrak{d}$, defined as the
as the largest number of hyperedges adjacent to any individual hyperedge.

We can separate the vertices into subregions, such as $V_{A}$, and we denote
with $\partial_{A}\in V_{A}$ the sites at the boundary of that region (that
is, with at least one hyperedge connecting to $\setminus\Lambda_{A}$), of
which there are $|\partial_{A}|$. For simplicity, we often refer to regions as
$A,B,..$ instead of $V_{A},V_{B},...$. We also need the notion of “distance”
between two regions, $\text{dist}(A,B)$, defined as the smallest number of
overlapping hyperedges that connect a vertex from $A$ with a vertex from $B$.

To define the Hamiltonian, we associate local interactions to hyperedges, such
that $H=\sum_{i\in E}h_{i}$. For an operator $h_{i}$, the set of vertices on
which it has non-trivial support is $\text{supp}(h_{i})$. We have already
specified that each $h_{i}$ is such that $|\text{supp}(h_{i})|\leq k$ (that
is, the hyperedges have size at most k), so that for constant $k$,
$N\propto|E|$. We also note that that $||h_{i}||\leq h$ and introduce the
following quantity

$J=\max_{x\in V}\sum_{i:x\in\text{supp}(h_{i})}||h_{i}||,$ (16)

that is, $J$ upper bounds the norm of the interactions that act on an
individual vertex.

###  II.4 Asymptotic notation

The so-called asymptotic or Bachmann–Landau notation succinctly describes the
asymptotic behaviour of a function when the argument grows large. It is
typically used when in a particular expression there are constant factors that
we are happy to omit, that are unnecessarily cumbersome, or when we only have
partial knowledge of the exact expression but know the asymptotic behaviour.
We say that , given functions $f(N),g(N)\geq 0$:

  * •

$f(N)=\mathcal{O}(g(N))$ if there are constants $M,N_{0}>0$ such that $\forall
N>N_{0}$, $f(N)\leq Mg(N)$.

  * •

$f(N)=\tilde{\mathcal{O}}(g(N))$ is similar to $\mathcal{O}(g(N))$ but with
possible additional poly-logarithmic factors, so that instead $\forall
N>N_{0}$, $f(N)\leq Mg(N)\text{polylog}(g(N))$.

  * •

$f(N)=o(g(N)$ if for every $\varepsilon>0$ there exists a $N_{0}>0$ such that
$\forall N>N_{0}$, $f(N)\leq\epsilon g(N)$.

  * •

$f(N)=\Omega(g(N))$ if there are constants $M,N_{0}>0$ such that $\forall
N>N_{0}$, $f(N)\geq Mg(N)$.

These are the most commonly used symbols of this notation, all of which appear
below.

##  III An overview of technical tools

When studying quantum thermal states from a mathematical point of view, what
we often need is some way of simplifying the operator $e^{-\beta H}$, in a way
that makes the particular problem at hand mathematically tractable. This is
usually achieved by expressing the relevant function of $e^{-\beta H}$ in
simpler terms. Potential issues that complicate this are:

  1. 1.

The exponential of a local operator is not a local operator, due to the high
order terms in the expansion, and could in principle be arbitrarily
complicated.

  2. 2.

The individual elements in the Hamiltonian Eq. (2) do not commute with each
other. Thus we cannot divide the exponential of the Hamiltonian into smaller
pieces by iterating simple identities like
$e^{-\beta(H_{1}+H_{2})}\stackrel{{\scriptstyle?}}{{=}}e^{-\beta
H_{1}}e^{-\beta H_{2}}$.

The locality of the Hamiltonian helps make these two problems often not as
serious as they could be in general situations. There is a number of tools to
deal with this, and we now describe some of the most relevant ones. Below, we
explain how the cluster expansion in Sec. III.1 helps with issue 1, while
there are at least two different techniques in Sec. III.2 and III.3 that help
us with issue 2.

###  III.1 Connected cluster expansion

This is a powerful set of ideas whose origins can be traced back to a wide set
of the classic (and classical) literature on mathematical physics and
statistical mechanics (see e.g. Ruelle (1969)) , initiated in Mayer and
Montroll (1941). It has traditionally been used to prove the analyticity of
the partition function at high temperatures and other regimes, so it serves as
an ideal tool to characterize the completely analytical interactions from Sec.
I.2. More recently, it has also been used to study the existence of
computationally efficient approximation schemes to it (see e.g. Mann and
Helmuth (2021)). The technique is flexible and general enough that it can also
cover objects beyond partition functions, such as characteristic functions and
other related quantities.

For simplicity we here focus on the high temperature expansion 222This type of
expansion is not limited to high temperatures, it can be used to expand around
any other parameter, such as a local magnetization. See e.g. Chapter 5 of
Friedli and Velenik (2017). . The starting point is the logarithm of the
partition function $\log Z\equiv\log\text{Tr}[e^{-\beta H}]$. Let us consider
its Taylor expansion around $\beta=0$

$\log Z=\sum_{m}\frac{\beta^{m}}{m!}K_{m}.$ (17)

One can then ask: what is the radius of convergence of this Taylor series?
More precisely, we would like to know whether there is some $\beta^{*}$
independent of the system size such that for $0\leq\beta<\beta^{*}$ we have
that:

  * •

The function $\log Z$ is analytic.

  * •

The $m$-th derivative at $\beta=0$ is such that

$\left|\frac{\text{d}^{m}\log Z}{\text{d}\beta^{m}}\right|=|K_{m}|\leq
C_{1}N(\beta/\beta^{*})^{m}m!,$ (18)

for some constant $C_{1}$.

  * •

The truncated Taylor series gives a good approximation as

$\left|\log Z-\sum_{m=0}^{M}\frac{\beta^{m}}{m!}K_{m}\right|\leq
C_{1}N\frac{(\beta/\beta^{*})^{(M+1)}}{1-(\beta/\beta^{*})}.$ (19)

There are various ways to narrow down the radius of convergence of this
series, but they all revolve around the idea of writing $\log Z$ in terms of
connected clusters.

A cluster is a multiset (that is, a set counting multiplicities) of
Hamiltonian terms $h_{i}$ (or alternatively, of hyperedges $\\{i\in E\\}$),
which can appear more than once. A given cluster W has size $|\textbf{W}|$
equal to the number of elements in the multiset (counting multiplicities
$\mu^{\textbf{W}}_{i}$, so that
$|\textbf{W}|=\sum_{\\{i\in\textbf{W}\\}}\mu^{\textbf{W}}_{i}$). Moreover, W
is connected if the hypergraph with hyperedges $i\in\textbf{W}$ is connected.
Let us define the set of all clusters of size at most $|\textbf{W}|=m$ with
$\mathcal{C}_{m}$, and the set of all connected clusters as $\mathcal{G}_{m}$.
For instance, $\mathcal{G}_{1}$ is the set of $\\{h_{i}\\}$, $\mathcal{G}_{2}$
are the pairs $\\{h_{i},h_{j}\\}$ provided that $i,j$ are adjacent or $i=j$
(in which case $\mu_{i}^{\textbf{W}}=2$). See Fig. 1 further illustrations of
a connected and a disconnected cluster.

Figure 1:  Illustration of the clusters, defined as a multiset of interaction
terms or hyperedges. In this example, the interaction is a graph on a square
lattice. The cluster on the left is connected $\textbf{W}\in\mathcal{G}_{m}$,
while the one on the right is disconnected. The thickness of the lines
represents the multiplicities $\mu_{i}^{\textbf{W}}$ of the edges, which may
appear any number of times in a cluster as long as
$\sum_{i}\mu_{i}^{\textbf{W}}=m$.

Now, let us define the Hamiltonian with auxiliary variables
$\\{\lambda_{i}\\}$ as $H(\lambda)=\sum_{i}\lambda_{i}h_{i}$. We use this to
introduce the cluster derivative

$\mathcal{D}_{\textbf{W}}=\left.\prod_{i\in\textbf{W}}\left(\frac{\partial}{\partial\lambda_{i}}\right)^{\mu^{\textbf{W}}_{i}}\right|_{\lambda=0}.$
(20)

Here, the subscript $\lambda=0$ means to set $\lambda_{i}=0$ for all $i$ after
taking the derivatives. We thus write

$\beta^{m}K_{m}=\sum_{\textbf{W}\in\mathcal{C}_{m}}\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta
H(\lambda)}].$ (21)

What we have done here is to simply write each moment of the Taylor series as
a sum of the contributions of all clusters W, without further specifying what
each contribution looks like.

We now prove the key simplification stemming from this expression. Let
$\textbf{W}\notin\mathcal{G}_{m}$, so that we have
$\textbf{W}=\textbf{W}_{1}\cup\textbf{W}_{2}$ where
$\textbf{W}_{1},\textbf{W}_{2}$ are non-overlapping clusters. This allowsus to
define $h_{\textbf{W}_{1}},h_{\textbf{W}_{2}}$ as the Hamiltonian terms in
those clusters, so that
$\text{supp}(h_{\textbf{W}_{1}})\cap\text{supp}(h)_{\textbf{W}_{2}}=\O$. We
then have

$\displaystyle\mathcal{D}_{\textbf{W}}$ $\displaystyle\log\text{Tr}[e^{-\beta
H(\lambda)}]=\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta(h_{\textbf{W}_{1}}(\lambda)+h_{\textbf{W}_{2}}(\lambda))}]$
$\displaystyle=\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta
h_{\textbf{W}_{1}}(\lambda)}]+\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta
h_{\textbf{W}_{2}}(\lambda)}]$ $\displaystyle=0.$ (22)

This means we can write the moments in terms of connected clusters only

$\beta^{m}K_{m}=\sum_{\textbf{W}\in\mathcal{G}_{m}}\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta
H(\lambda)}].$ (23)

This reduces the number of contributions to $K_{m}$ dramatically, and makes it
possible to estimate them. One way to show the convergence of the series (see
Kuwahara and Saito (2020a); Haah  _et al._ (2022); Wild and Alhambra (2023))
is to prove the following:

  * •

The number of connected clusters of size $m$ is bounded by $Nc_{1}^{m}$ for
some constant $c_{1}$ Malyshev (1980); Dobrushin (1996).

  * •

The size of each cluster derivative for a cluster of size $|\textbf{W}|=m$ is
at most

$\left|\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta
H(\lambda)}]\right|\leq(\beta c_{2})^{m}m!$ (24)

for some constant $c_{2}$ Haah  _et al._ (2022); Wild and Alhambra (2023).

The constants here can usually be taken to be simple functions of the lattice
parameters $\mathfrak{d},k,J,h$ and of some property of the interaction graph.
For instance, in Haah  _et al._ (2022); Wild and Alhambra (2023), it is shown
that $c_{1}=e\mathfrak{d}$, and that $c_{2}=2eh(\mathfrak{d}+1)$ (see Eq. (3)
for the definition of $h$). These facts together imply that $|K_{m}|\leq
2e^{2}h\mathfrak{d}(\mathfrak{d}+1)\equiv(\beta^{*})^{-1}\simeq\mathcal{O}(h)$,
so that the partition function is analytic within a disk in the complex plane
of radius $\beta<\beta^{*}$, and is also well approximated by its Taylor
series.

Beyond this argument for the convergence of the series, there are other more
general abstract methods for proving convergence of this type of quantity, in
terms of the so-called _polymer models_ Kotecký and Preiss (1986); Dobrushin
(1996); Fernández and Procacci (2007); Mann and Helmuth (2021). See Chapter 5
of Friedli and Velenik (2017) for an introduction.

So far we have only discussed convergence of the series. However, the cluster
expansion can be used to device efficient approximation schemes. The main idea
is to prove that the individual Taylor terms can be computed efficiently. This
requires two separate steps:

  * •

The set of all clusters of size $m$ can be enumerated in time
$\text{poly}(N)\times\exp(\mathcal{O}(m))$ Helmuth  _et al._ (2020); Haah  _et
al._ (2022)

  * •

Each cluster derivative can be computed exactly in time
$\text{poly}(N)\times\exp(\mathcal{O}(m))$ Björklund  _et al._ (2007); Mann
and Helmuth (2021); Haah  _et al._ (2022).

We can thus add all the contributions from the different derivatives to obtain
$K_{m}$ in time $\text{poly}(N)\times\exp(\mathcal{O}(m))$. This, together
with Eq. (19), implies that by calculating the Taylor series up to a degree
$M=\mathcal{O}\left(\log(N/\epsilon)\times\log(\beta^{*}/\beta)\right)$ there
exists an $\epsilon$-close additive approximation to $\log Z$ that can be
computed in time $\text{poly}(N,\epsilon^{-1})$.

We do not expect to be able to prove many general statements at all
temperatures, due to the presence of thermal phase transitions, and to the
fact that the ground state energy is computationally hard to estimate Kempe
_et al._ (2005). However, there are specific models in the literature for
which the convergence can be guaranteed for larger ranges of temperatures (see
e.g. Helmuth and Mann (2022) and references therein). See Sec. VII.2 for more
details.

This same technique also allows for e.g. the computation of expectation values
such as $\text{Tr}[h_{i}e^{-\beta H}/Z]$ by differentiating by an extra
$\lambda_{i}$ in the cluster derivative. It can be also applied to other
similar objects such as characteristic functions of the form
$\text{Tr}[e^{\alpha A}e^{-\beta H}/Z]$ for some local observable $A$, which
allows for the derivation of probability theory statements, as explained
further in Sec. VI.1.

###  III.2 Thermal locality estimates

We now show the first method to decompose the thermal state into a product of
smaller operators despite the non-commutativity, which is related to the
general idea of operator growth. Consider an operator $A$ with local support
on some small region on the lattice. For simplicity, let this region be such
that $|\text{supp}(A)|\leq k$.

An interesting quantity to study is the operator evolved in Euclidean or
imaginary time $\beta$ under the Hamiltonian $H$,

$A(i\beta)=e^{-\beta H}Ae^{\beta H}.$ (25)

This is in analogy with the Heisenberg-picture operator
$A(t)=e^{itH}Ae^{-itH}$, which can be understood in terms of the well-known
Lieb-Robinson bounds Lieb and Robinson (2004), that state that the support of
$A(t)$ is mostly confined to a linear lightcone. In many situations, one will
want to choose $A$ here to be one of the $h_{i}$ operators.

It then makes sense to ask the following question: what is the locality of the
Euclidean-evolved operator $A(i\beta)$? Perhaps surprisingly, this can be
dramatically different to the real-time case: there is no general linear
growth with the inverse temperature $\beta$, but a much wilder dependence on
it.

The main difference is that $e^{-itH}$ is a unitary matrix, while $e^{-\beta
H}$ is not. This means that results that exploit unitarity, such as the
aforementioned Lieb-Robinson bounds, do not apply straightforwardly. Our best
way forward seems then to analyze $A(i\beta)$ in terms of nested commutators

$\displaystyle A(i\beta)$
$\displaystyle=\sum_{m=0}^{\infty}\frac{(-\beta)^{m}}{m!}[H,[H,....,[H,A]...]$
$\displaystyle\equiv\sum_{m=0}^{\infty}\beta^{m}C_{m}(A).$ (26)

It is easy to see that the $m$-th term in this expansion has support on a
connected region whose furthermost point is a distance $m$ away from $A$. The
question then becomes: how does this expansion in terms of $\beta$ converge?
We now discuss how the answer to this question (either high temperatures and
1D) appears to also be restricted to the one-phase region of completely
analytical interactions from Sec. I.2.

It can be shown that, for general interaction graphs, Ruelle (1969); Kuwahara
_et al._ (2016)

$||C_{m}(A)||\leq k||A||(2Jk)^{m},$ (27)

with $J,k$ as defined in Sec. II.3. This statement is very much related to the
bound on the number of connected clusters in Sec. III.1 above, since only
connected clusters contribute to the nested commutators. Eq. (27) implies that
as long as $\beta<(2Jk)$, the expansion can be controlled as a geometric
series, from which we obtain

$\displaystyle||A(i\beta)||\leq k||A||\frac{1}{1-2\beta Jk}$ (28)
$\displaystyle||A(i\beta)-\sum_{m=0}^{M}\beta^{m}C_{m}(A)||\leq
k||A||\frac{(2\beta Jk)^{M+1}}{1-2\beta Jk}.$ (29)

Given that the $m$-th nested commutator can have support on at most $k\times
m$ sites, the latter equation means that $A(i\beta)$ is, roughly speaking,
localized within the subset of vertices a distance at most $k\times m$ away
from $\text{supp}(A)$. It is known that one cannot extend this result to
temperatures lower than $\beta\simeq\mathcal{O}(1)$, since there exists a 2D
lattice in which the terms in the nested commutators in Eq. (26) add up
constructively, in a way that the norm of $A(i\beta)$ grows with system size,
and diverges as $N\rightarrow\infty$ Bouch (2015).

On the other hand, it has been known for some time Araki (1969) that when the
lattice is a one-dimensional chain, the nested commutators grow more slowly,
so that this type of convergence happens for all temperatures. For simplicity,
we show explicitly the result for $k=2$ combining Pérez-García and Pérez-
Hernández (2020) and Bouch (2015), which is

$\displaystyle||A(i\beta)||\leq||A||f(\beta,J)\exp(f(\beta,J))$ (30)
$\displaystyle||A(i\beta)-\sum_{m=0}^{M}\beta^{m}C_{m}(A)||\leq
15||A||e^{-(M+1)}$ $\displaystyle\quad\quad\quad\forall M>g(\beta,J).$ (31)

Here, we have defined $f(\beta,J)\equiv 16\beta J\exp(1+8\beta J)$ and
$g(\beta,J)\equiv\exp(240e^{2}\beta J)-1$ 333The numerical constants in these
functions are likely sub-optimal, but this is very rarely important for
applications.. The intuitive reason for these is that the geometric bound of
Eq. (27) can be improved in this case as Bouch (2015) (again, for $k=2$)

$||C_{m}(A)||\leq 15||A||\left(\frac{240eJ}{\log(m+1)}\right)^{m}.$ (32)

Notice that, because of the logarithm, the series in Eq. (26) is not
geometric, and converges for all $\beta$. For further explanations of these
points see also Avdoshkin and Dymarsky (2020).

So far, we have described how does $A(i\beta)$ approximate its Taylor
expansion. An alternative approximation commonly considered is to the operator
$e^{-\beta H_{\Lambda_{m}}}Ae^{\beta H_{\Lambda_{m}}}$, where
$H_{\Lambda_{m}}=\sum_{\text{supp}(h_{i})\in\Lambda_{m}}h_{i}$ and
$\Lambda_{m}$ is a subset of the hypergraph corresponding to some region
$V_{m}$ of the full lattice $\Lambda$. One can then consider how the norm

$||A(i\beta)-e^{-\beta H_{\Lambda_{m}}}Ae^{\beta H_{\Lambda_{m}}}||$ (33)

decays with $m$ in terms of how $\Lambda_{m}$ is defined (typically, some
hyper-sphere centered around $A$). The analysis and convergence turn out to be
almost the same as the one for the moments $C_{m}(A)$ above. The reason is
that the difference between $\sum_{m=0}^{M}\beta^{m}C_{m}(A)$ and $e^{-\beta
H_{\Lambda_{m}}}Ae^{\beta H_{\Lambda_{m}}}$ are essentially the higher order
terms in $\beta$ of the latter, which are also suppressed. See e.g. Pérez-
García and Pérez-Hernández (2020) for a detailed analysis of the 1D case or
e.g. Lemma 20 in Kuwahara  _et al._ (2021) for a proof in higher dimensions.

One of the main reasons why both of these approximations are interesting is
that they are related to the following propagator

$E_{A}\equiv e^{-\beta(H+A)}e^{\beta
H}=\mathcal{T}\left(e^{-\int_{0}^{\beta}A(s)\text{d}s}\right),$ (34)

where $A(s)=e^{-sH}Ae^{sH}$ and $\mathcal{T}$ denotes the usual time-ordered
integral. This is such that

$e^{-\beta(H+A)}=E_{A}e^{-\beta H}.$ (35)

This operator $E_{A}$ can be used, for instance, to decompose $e^{-\beta H}$
as a product of its parts by e.g. choosing $A$ as the Hamiltonian at the
boundary of two regions. This operator can be analyzed through a usual Dyson
series in terms of powers of $e^{-xH}Ae^{xH}$. Assuming $\beta<(2Jk)^{-1}$, it
can be shown that $E_{A}$ has bounded norm as it follows from Eq. (28) and
(34) that

$||E_{A}||\leq\exp\left(\int_{0}^{\beta}\text{d}s||e^{-sH}Ae^{sH}||\right)\leq\left(\frac{1}{1-2\beta
Jk}\right)^{\frac{||A||}{2\beta J}}.$ (36)

In Appendix B we also show that it is approximately localized in a similar way
as $e^{-xH}Ae^{xH}$ is. This means that there exist an operator $E_{A}(l)$
with support restricted to a distance at most $l$ away from $A$ such that for
$\beta<(2Jk)^{-1}$

$||E_{A}-E_{A}(l)||\leq\beta k||A||\frac{(2\beta Jk)^{l+1}}{(1-2\beta
Jk)^{\frac{||A||}{2\beta J}+1}}.$ (37)

Also, notice that if $[H,A]=0$, then $E_{A}=e^{-\beta A}$. With the right
choice of $H,A$, the operator $E_{A}$ can thus be thought of as a “transfer
operator”. Corresponding results also exists for 1D using Eq. (30) and (31).

Alternatively, one can also define the following operator

$E^{\prime}_{A}\equiv e^{-\beta(H+A)}e^{\beta H}e^{\beta A},$ (38)

with the difference that $H$ and $A$ are now treated on equal footing. In this
case, $E^{\prime}_{A}$ is just the multiplicative error term in the first
order Trotter product formula, which can be similarly analyzed through the
expansion of $A(i\beta)$ (see the thorough analysis of Trotter errors in
Childs  _et al._ (2021) for more details). These Trotter errors are most
commonly analyzed in the context of digital quantum simulation Lloyd (1996),
for which it is often convenient to go to higher orders in the decomposition.

We finish this subsection outlining a result in 1D related to this discussion,
which follows from bounds on the quantity in Eq. (33). It appeared first in
Araki (1969), and it features in Sections V.1.1 and VII.1. Let us define
$E^{l}_{A}=e^{-\beta(H_{l}+A)}e^{\beta H_{l}}$, where $H_{l}$ are the
interaction terms a distance at most $l$ away from $\text{supp}(A)$. It can be
shown that

$\displaystyle||E_{A}||$ $\displaystyle\leq C_{1}$ (39)
$\displaystyle||E_{A}-E^{l}_{A}||$ $\displaystyle\leq
C_{2}\frac{q^{l}}{(l+1)!},$ (40)

where $C_{1},C_{2}$ and $q>1$ are constants depending on $k,J,\beta$ which we
do not show explicitly for simplicity, although notice that $C_{1}$ will be
essentially the exponential of Eq. (30) . The proof is similar to that of
Appendix B, together with a bound on Eq. (33). We refer the reader to e.g.
Araki (1969); Pérez-García and Pérez-Hernández (2020) for further details.

The approximations $E_{A}(l)$ and $E_{A}^{l}$ to the operator $E_{A}$ are
important in that they allow us to decompose $e^{-\beta H}$ into a product of
smaller local operators despite the Hamiltonian being non-commuting. For
instance, they will be useful in the arguments of Sec. V.1.1.

###  III.3 Quantum belief propagation

An idea related to the previous locality estimates appeared first Hastings
(2007a), and has more recently featured in several results about Gibbs states
on lattices Kim (2012); Kato and Brandão (2019); Ejima and Ogata (2019);
Brandão and Kastoryano (2019); Harrow  _et al._ (2020); Kuwahara  _et al._
(2021); Anshu  _et al._ (2021a). It is a tool similar to that of the previous
Sec. III.2, in that it also allows us to decompose the thermal state as a
product of smaller, localized operators, which makes certain calculations more
tractable. The goal is to be able to divide the big operator $e^{-\beta H}$
into smaller pieces, that allow, for instance, to prove that the Gibbs state
can be sequentially generated, or that local perturbations only have effect in
the near vicinity.

This is part of a celebrated series of works including the decay of
correlations for gapped ground states Hastings and Koma (2006) or the area law
of entanglement in one dimension Hastings (2007b) which show how Lieb-Robinson
bounds (a dynamical statement) can be used to prove static properties about
ground and thermal states. The derivation here is a particular example of that
idea, but see Hastings (2010, 2021) for overviews that go beyond Gibbs states.

The goal is to construct a quasi-local operator $O_{A}^{m}$ (to be defined
below) with support near $\text{supp}(A)$ such that

$e^{-\beta(H+A)}\simeq O_{A}^{m}e^{-\beta H}(O_{A}^{m})^{\dagger}.$ (41)

Notice the difference with Eq. (35), where we only multiply $e^{-\beta H}$
with an operator from the left.

We start by considering the “perturbed” Hamiltonian $H(s)=H+sA$ and the
following derivative

$\frac{de^{-\beta H(s)}}{ds}=-\frac{\beta}{2}\left\\{e^{-\beta
H(s)},\Phi_{\beta}^{H(s)}(A)\right\\}$ (42)

where, if $H(s)=\sum_{i}E_{i}(s)|{i(s)}\rangle\langle{i(s)}|$ is the energy
eigenbasis,

$\Phi_{\beta}^{H(s)}(A)_{ij}=\langle
i(s)|A|j(s)\rangle\tilde{f}_{\beta}(E_{i}(s)-E_{j}(s)),$ (43)

where $\tilde{f}_{\beta}(\omega)=\frac{\tanh(\beta\omega/2)}{\beta\omega/2}$.
With
$f_{\beta}(t)=\frac{4}{\beta\pi}\log\left(\frac{e^{\pi|t|}{\beta}+1}{e^{\pi|t|}{\beta}-1}\right)$
the Fourier transform of $\tilde{f}_{\beta}(\omega)$ (see Appendix B of Anshu
_et al._ (2021a)), we can also write

$\Phi_{\beta}^{H(s)}(A)=\int_{-\infty}^{\infty}\text{d}tf_{\beta}(t)e^{-itH(s)}Ae^{itH(s)}.$
(44)

The proof leading to Eq. (42) that explains hte appearance of
$\tilde{f}_{\beta}(\omega)$ is shown in App. B.

Since $||e^{-itH(s)}Ae^{itH(s)}||=||A||$,
$||\Phi_{\beta}^{H(s)}(A)||\leq||A||$ by the triangle inequality and Eq.
(184). Moreover, it can also be approximated by a localized operator around
the support of $A$ by using Lieb-Robinson bounds Lieb and Robinson (2004). In
particular, when $H(s)$ is local, it can be shown that Barthel and Kliesch
(2012)

$\displaystyle||e^{-itH(s)}Ae^{itH(s)}$
$\displaystyle-e^{-itH_{\Lambda_{m}}(s)}Ae^{itH_{\Lambda_{m}}(s)}||$
$\displaystyle\leq m^{D-1}b||A||e^{c^{\prime}(vt-m)},$ (45)

where $H_{\Lambda_{m}}(s)$ is the restriction of $H(s)$ to the sum of local
terms that are at most a distance $m$ away from the support of $A$, $v$, $b$
and $c^{\prime}$ are constants, and $D$ is the dimension of the interaction
lattice. This should be reminiscent of the $H_{\Lambda_{m}}$ appearing in Eq.
(33), with the only difference being that now the evolution is for real times.

Eq. (45) allows us to establish that after time-evolving $A$ for a short time,
the support is still approximately localized, with a radius growing with time.
This also allows us to define $\Phi_{\beta}^{H_{\Lambda_{m}}(s)}(A)$, which is
close to the original $\Phi_{\beta}^{H(s)}(A)$ in the following sense

$\displaystyle\frac{1}{||A||}||\Phi_{\beta}^{H(s)}(A)-\Phi_{\beta}^{H_{\Lambda_{m}}(s)}(A)||$
(46) $\displaystyle\leq
m^{D-1}b\int_{-m/2v}^{m/2v}e^{c^{\prime}(v|t|-m)}f_{\beta}(t)\text{d}t+4\int_{m/2v}^{\infty}f_{\beta}(t)\text{d}t$
$\displaystyle\leq m^{D-1}be^{-c^{\prime}m/2}+\frac{1}{e^{\frac{\pi m}{2\beta
v}}-1},$

where in the first line, after the triangle inequality, we divided the
integral into two different ranges, and used Eq. (45) in the first range and
$||\Phi_{\beta}^{H(s)}(A)||\leq||A||$ in the second. The bounds on the
integrals were obtained using the properties of $f_{\beta}(t)$ from Eq. (184)
and (186). As a result, for $m$ large enough, the difference between the two
operators is exponentially decaying in $m$.

We can now integrate Eq. (42) between $s=0$ and $s=1$ to obtain

$e^{-\beta(H+A)}=O_{A}e^{-\beta H}O_{A}^{\dagger},$ (47)

where

$O_{A}=\mathcal{T}e^{-\frac{\beta}{2}\int_{0}^{1}\text{d}s\Phi_{\beta}^{H(s)}(A)}.$
(48)

Similarly to Eq. (34) above (see the analyisis of the operator $E_{A}$ in
Appendix B), this operator has a bounded norm, since

$\displaystyle||O_{A}||\leq
e^{\frac{\beta}{2}\int_{0}^{1}\text{d}s||\Phi_{\beta}^{H(s)}(A)||}=e^{\frac{\beta}{2}||A||},$
(49)

Additionally, it is also approximately localized exactly around the support of
$A$. Let us define the operator $O_{A}^{m}$ in the natural way

$O_{A}^{m}\equiv\mathcal{T}e^{-\frac{\beta}{2}\int_{0}^{1}\text{d}s\Phi_{\beta}^{H_{\Lambda_{m}}(s)}(A)}.$
(50)

We can use an argument analogous to that used in App. B to show that $O_{A}$
and $O_{A}^{m}$ are exponentially close, as, given Eq. (46), for large enough
$m$,

$\displaystyle||O_{A}-O_{A}^{m}||$
$\displaystyle\leq\frac{\beta||O_{A}||}{2}\int_{0}^{1}\text{d}s||\Phi_{\beta}^{H(s)}(A)-\Phi_{\beta}^{H_{\Lambda_{m}}(s)}(A)||$
(51) $\displaystyle\leq\beta||A||e^{\frac{\beta||A||}{2}-\Omega(m)}.$ (52)

These bounds can be compared to Eq. (39) and (40), which are of a very similar
nature. There are, however, two important differences between $O_{A}$ and
$E_{A}$ in Eq. (34) above:

  * •

Since it is based on the Lieb-Robinson bound, the operator $O_{A}$ is well-
behaved in all lattices and at all temperatures, in the sense that it has a
bounded norm and is approximately localized. This is as opposed to $E_{A}$,
which is likely a large operator in high dimensions and low temperatures. That
is, this result holds for all Gibbs states, irrespective of whether they are
completely analytical interactions or not.

  * •

On the other hand, to recover $e^{-\beta(H+A)}$ from $e^{-\beta H}$ we require
left and right multiplication with $O_{A},O_{A}^{\dagger}$, as opposed to Eq.
(35), which may be problematic in some applications. In particular, we should
not expect it to be a key ingredient in proving results that do not hold at
all temperatures such as the decay of correlations or the analyticity of the
partition function.

###  III.4 Selected trace inequalities

In the past two subsections we have explored how to analyze perturbations to
the Hamiltonian in the Gibbs operator $e^{-\beta(H+A)}$ via $E_{A}$ and
$O_{A}$. When considering traces, simpler identities hold. We exemplify this
with two with very elementary implications and proofs, which can be found in
(at least) Lenci and Rey-Bellet (2005). The first one is about the stability
of partition functions. Let $H_{1},H_{2}$ be Hermitian operators. Then we have

$\left|\log\text{Tr}[e^{H_{1}+H_{2}}]-\log\text{Tr}[e^{H_{1}}]\right|\leq||H_{2}||.$
(53)

The proof is just as follows:

$\displaystyle\left|\log\text{Tr}[e^{H_{1}+H_{2}}]-\log\text{Tr}[e^{H_{1}}]\right|$
$\displaystyle\quad=\left|\int_{0}^{1}\frac{\text{d}}{\text{d}t}\log\text{Tr}\left[e^{H_{1}+tH_{2}}\right]\text{d}t\right|$
$\displaystyle\quad\leq\int_{0}^{1}\left|\frac{\text{Tr}[H_{2}e^{H_{1}+tH_{2}}]}{\text{Tr}[e^{H_{1}+tH_{2}}]}\right|\text{d}t$
$\displaystyle\quad\leq||H_{2}||,$ (54)

where in the last inequality we have simply used Hölder’s inequality Eq. (7)
with $q_{1}=1,q_{2}=\infty$. If we take e.g. $H_{1}=-\beta H,H_{2}=-\beta A$,
this implies that changing the Hamiltonian by $A$ changes the log-partition
function at most by $\beta||A||$.

The second is a similar result that holds for expectation values of positive
operators. Let $H_{1},H_{2}$ be as before, and let $C>0$. Then

$\displaystyle\left|\log\text{Tr}[Ce^{H_{1}+H_{2}}]-\log\text{Tr}[Ce^{H_{1}}]\right|$
(55)
$\displaystyle\quad\leq\int_{0}^{1}\text{d}t\int_{-1/2}^{1/2}\text{d}s||e^{-s(H_{1}+tH_{2})}H_{2}e^{s(H_{1}+tH_{2})}||.$

The proof can be found in Appendix B. This norm can then be bounded with the
results from Sec. III.2, to scale as $\propto||H_{2}||$. The resulting
expression can for instance be used for analyzing characteristic functions of
observables $F$ by taking $C=e^{\alpha F}$ for $\alpha\in\mathbb{R}$ Lenci and
Rey-Bellet (2005).

More generally, in the practice of mathematical quantum physics, whether it is
from the many body, the QI, or any other perspective, many important proof
ingredients take the form of inequalities, either between operators, traces,
or norms (such as those already mentioned in Sec. II.1). There are too many to
give a reasonably complete overview here but we refer the reader to e.g.
Bhatia (2013); Carlen (2010).

##  IV Correlations

One of the more important questions when studying many body systems is: how
and how much are the different parts correlated? Intuitively, the stronger
these correlations, and the longer their range, the more complex a state is -
the reason being that we cannot think of the large system as a collection of
simpler, weakly correlated parts. The obvious extreme example is that of an
uncorrelated gas, in which the particles do not interact and have completely
independent properties.

For thermal states with local interactions, we can expect that locality will
make the state far from generic, in a way that constraints its complexity.
Intuitively, it should cause the correlations to be “localized”, meaning that
particles are only correlated with their vicinity as given by the lattice
geometry. For a rough intuition, consider the first terms of the Taylor series

$e^{-\beta
H}=\mathbb{I}-\beta\sum_{i}h_{i}+\frac{\beta^{2}}{2}\sum_{i,j}h_{i}h_{j}+...$
(56)

That is, at very high temperatures we approach the trivial uncorrelated state
$\propto\mathbb{I}$ and the leading order term includes only $k$-local
couplings, with only higher order terms coupling far away particles. We thus
expect that the correlations between particles will generally be weaker $i)$
the higher the temperature and $ii)$ the larger their distance on the
interaction graph. This is one of the main ways of understanding how the Gibbs
states of the “one phase region” of Sec. I.2 are very different from generic
states.

An important motivation for this is that, as we will see in later sections,
the situations in which the correlations are weaker or short range roughly
correspond to those in which we expect better algorithms for the description
of thermal states. This is perhaps most clearly the case in the context of
tensor network methods. We now proceed to describe (and even prove) the more
important ways in which these correlations are constrained.

###  IV.1 Correlations between neighbouring regions: Thermal area law

One of the more important statements about correlations in quantum many-body
systems is the area law. This roughly states that a measure of correlations
between two adjacent regions is upper-bounded by a number proportional to the
size of their mutual boundary.

Traditionally, this has been mostly studied in the context of ground states,
which are pure. There, the relevant measure of correlations is the
entanglement entropy, or some Rényi version of it. In that context, an area
law for the entanglement entropy is believed to hold for all ground states of
models with a gap Eisert  _et al._ (2010). This can be proven in 1D Hastings
(2007b); Arad  _et al._ (2013) and in some cases in 2D Anshu  _et al._
(2022a). The interest in it is largely due to its relation to other phenomena,
such as phase transitions Vidal  _et al._ (2003), the decay of long-range
correlations Brandão and Horodecki (2014) or the effectiveness of certain
tensor network algorithms Verstraete and Cirac (2006); Ge and Eisert (2016).

For thermal states, a very general area law can be shown to hold for systems
in any dimension, at all temperatures. We now give a short proof of this
statement, and then discuss its significance (see Wolf  _et al._ (2008) for
the original reference). In this case, since it is a mixed state, an
appropriate measure of correlations is the mutual information in Eq. (13).

Let us partition our interaction graph into two subsets of particles $A,B$,
with a thermal state $\rho_{\beta}^{AB}$. We start with the very simple
thermodynamic observation that the free energy $F$ from Eq. (12) of the
thermal state is lower than that of any other state (this follows from Eq.
(12)), and in particular

$F_{\beta}(\rho_{\beta}^{AB})\leq
F_{\beta}(\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}).$ (57)

Writing out the free energy explicitly as $F_{\beta}(\rho)=\text{Tr}[\rho
H]-\beta^{-1}S(\rho)$ and rearranging yields

$S(\rho_{\beta}^{A}\otimes\rho_{\beta}^{B})-S(\rho_{\beta}^{AB})\leq\beta\left(\text{Tr}[H\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}]-\text{Tr}[H\rho_{\beta}^{AB}]\right).$
(58)

Given that the entropy is additive $S(\rho\otimes\sigma)=S(\rho)+S(\sigma)$
notice that the LHS is exactly the mutual information
$I(A:B)_{\rho_{\beta}^{AB}}$ from Eq. (13). Since the Hamiltonian is local, we
can write it as

$H=H_{A}+H_{B}+H_{I},$ (59)

where $H_{A},H_{B}$ have support only on $A,B$ respectively, and $H_{I}$ is
the interaction between them (with support on both). By definition, the
expectation values of $H_{A}$ and $H_{B}$ coincide on both states
$\text{Tr}[(H_{A}+H_{B})\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}]=\text{Tr}[(H_{A}+H_{B})\rho_{\beta}^{AB}]$,
so that

$\displaystyle\beta\left(\text{Tr}[H\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}]-\text{Tr}[H\rho_{\beta}^{AB}]\right)$
(60)
$\displaystyle=\beta\left(\text{Tr}[H_{I}\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}]-\text{Tr}[H_{I}\rho_{\beta}^{AB}]\right).$

Now we can use a few of the operator inequalities from Section II.1 to obtain

$\displaystyle\text{Tr}[H_{I}\rho_{\beta}^{A}$
$\displaystyle\otimes\rho_{\beta}^{B}]-\text{Tr}[H_{I}\rho_{\beta}^{AB}]\leq||H_{I}(\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}-\rho_{\beta}^{AB})||_{1}$
$\displaystyle\leq||H_{I}||\times||\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}-\rho_{\beta}^{AB}||_{1}$
(61)
$\displaystyle\leq||H_{I}||\times(||\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}||_{1}+||\rho_{\beta}^{AB}||_{1})=2||H_{I}||.$

Putting Eq. (58) and (IV.1) together we have the final result

$I(A:B)_{\rho_{\beta}^{AB}}\leq 2\beta||H_{I}||.$ (62)

This is the _area law_ for the mutual information of a thermal state: it
implies that the strength of the correlations of systems $A,B$ cannot depend
on their size, but that it grows at most as their common boundary. For a local
Hamiltonian, we have that

$||H_{I}||\leq 2kh|\partial_{AB}|,$ (63)

where $\partial_{AB}=\partial_{A}\cup\partial_{B}$, $h$ is defined in Eq. (3)
and $k$ is the largest support of any $h_{i}$. Notice that with
$|\partial_{AB}|$ we do not mean the size of the boundary of systems $A,B$
together, but the number of elements of $\partial_{A}$ that are connected to
$\partial_{B}$ by hyperedges. We show this schematically in Fig. 2. This is to
be contrasted with the most general upper bound on the mutual information,
which is $I(A:B)\leq\min\\{\log(d_{A}),\log(d_{B})\\}$ (since $\log
d_{A}\propto|A|$ the largest possible scaling is a “volume law” instead).

.

Figure 2:  Under an area law, the correlations between regions A and B grow at
most as the size of their mutual boundary $\partial_{AB}$.

What this strongly suggests (although it does not quite prove) is that the
correlations between $A$ and $B$ are localized around the mutual boundary, and
that the bulks of $A$ and $B$ are mostly uncorrelated. That is, the only
relevant information about $A$ that $B$ contains is about the region of $A$
that is near their boundary.

This statement, as can be seen from the proof, holds for all temperatures and
all interaction graphs, which is likely as general as it can be , at least for
systems with finite-dimensional Hilbert spaces Lemm and Siebert (2022). The
drawback of that generality, however, is that it will be unable to signal
important phenomena that happens only at specific temperature ranges, such as
thermal phase transitions, or an efficient classical or quantum simulability.
That is, Eq. (62) does not narrow down the set of “completely analytical
interactions” for Sec. I.2 in any meaningful way. Other more specific versions
of the thermal area law in the literature may have more potential in this
regard Sherman  _et al._ (2016); Scalet  _et al._ (2021).

The temperature dependence of Eq. (62) can be improved to
$\tilde{\mathcal{O}}(\beta^{2/3})$ Kuwahara  _et al._ (2021). This can be
proven with a variety of techniques, including those of Sec. III.2 and Sec.
III.3, as well as methods originated in the study of ground states Arad  _et
al._ (2013). This dependence is not far from optimal, since there exists a 1D
model for which the scaling of the MI is at least $\mathcal{O}(\beta^{1/5})$
at low temperatures Gottesman and Hastings (2010).

Many important physical models have a very different temperature dependence,
such as $\log(\beta+1)$ Znidaric  _et al._ (2008); Bernigau  _et al._ (2015).
Classical systems, on the other hand, have an upper bound that is independent
of the temperature, as $I(A:B)_{\rho_{\beta}^{AB}}\leq|\partial_{AB}|\log d$
Wolf  _et al._ (2008). All these suggests that the scaling of the mutual
information with $\beta$ in the low temperature regime is related to the
computational complexity of the ground space of the models.

###  IV.2 Decay of long-range correlations

An important fact about thermal states is that often their spatially separated
parts are very weakly correlated. Let $C,D$ be regions such that their
distance is $\text{dist}(C,D)$ (see Fig. 3). We focus on measures of
correlations evaluated at the marginals on these regions
$\text{Tr}_{\setminus(CD)}[\rho_{\beta}]=\rho_{\beta}^{CD}$. For instance,
taking the mutual information, we expect that in general

$I(C:D)_{\rho_{\beta}^{CD}}\leq f\left(\text{dist}(C,D)\right),$ (64)

where $f$ is some rapidly decaying function. In fact, we expect that for
completely analytical interactions $f(l)\leq
K|\partial_{C}||\partial_{D}|e^{-l/\xi}$, where $K>0$ is some constant,
$\partial_{C,D}$ is the size of the boundary of each region, and $\xi$ is the
thermal _correlation length_ that depends on the temperature and other
parameters, but not on $l$ or the system’s size.

This has been proven in translation invariant 1D systems at all temperatures
Bluhm  _et al._ (2022). The main idea behind it is to use the locality
estimates from III.2, and in particular the properties of the operator $E_{A}$
in Eq. 35. That this also holds for high enough temperatures in all dimensions
follows the cluster expansion applied to the mutual information Kuwahara  _et
al._ (2021).

.

Figure 3: Regions $C,D$ in the lattice are separated by a distance
$\text{dist}(C,D)$. The mutual information between these two regions typically
decays exponentially with their distance.

A more commonly stated but weaker condition is the decay of the connected two-
point correlators. This usually takes the form

$\displaystyle\frac{\left|\text{Tr}[\rho_{\beta}M_{C}\otimes
M_{D}]-\text{Tr}[\rho_{\beta}M_{C}]\text{Tr}[\rho_{\beta}M_{D}]\right|}{||M_{C}||||M_{D}||}$
$\displaystyle\leq K|\partial_{C}||\partial_{D}|e^{-\text{dist}(C,D)/\xi},$
(65)

where here $M_{C}$ and $M_{D}$ have support on regions $C,D$, respectively.
That this is weaker than the decay of the mutual information follows from
Pinsker’s inequality applied to the marginal on regions $C,D$.

It is known that correlations decay exponentially at large temperatures for
arbitrary interaction graphs Kliesch  _et al._ (2014); Fröhlich and Ueltschi
(2015). This can be shown with the cluster expansion Ueltschi (2004). Instead
of showing the proof in full generality, we can already see a simple but
instructive case by noticing that, with the notation of Sec. III.1,

$\displaystyle\frac{d\log
Z(\lambda)}{d\lambda_{i}d\lambda_{j}}\Big{|}_{\lambda=1}=\beta^{2}(\langle
h_{i}h_{j}\rangle-\langle h_{i}\rangle\langle h_{j}\rangle).$ (66)

At the same time, when one differentiates over two variables
$\lambda_{i},\lambda_{j}$, the nonzero contributions come from clusters that
contain both

$\displaystyle\frac{d\log Z(\lambda)}{d\lambda_{i}d\lambda_{j}}$ (67)
$\displaystyle=\frac{d}{d\lambda_{i}d\lambda_{j}}\left(\sum_{m}\frac{1}{m!}\sum_{\begin{subarray}{c}\textbf{W}\in\mathcal{G}_{m}\\\
i,j\in\textbf{W}\end{subarray}}\mathcal{D}_{\textbf{W}}\log\text{Tr}[e^{-\beta
H(\lambda)}]\right).$

However, connected clusters such that $i,j\in\textbf{W}$ belong to
$\mathcal{G}_{m}$ with $m\geq\text{dist}(i,j)$. This means that the lowest
moment that appears in the correlation function is $K_{m}$ with
$m=\text{dist}(i,j)$. If Eq. 18 holds, then the correlation function will
decay at least as $(\beta/\beta^{*})^{\text{dist}(i,j)}$, mirroring (65). A
similar argument also holds for arbitrary few body-observables $h_{C},h_{D}$
if one consider an appropriate cluster expansion of the perturbed Hamiltonian
$H+h_{C}+h_{D}$. See e.g. Ueltschi (2004); Fröhlich and Ueltschi (2015) for
more general results.

The connection between this type of correlation decay and the analiticity of
the partition function is very well understood in the classical case, where
they are known to be equivalent Dobrushin and Shlosman (1987). In the quantum
case, it is only known that a condition stronger than the analiticity of $\log
Z$ (called _analiticity after measurement_) implies decay of correlations. See
Harrow  _et al._ (2020) for details.

Intuitively, both these properties are related to the absence of phase thermal
phase transitions: at those critical points, the correlation function diverges
and the partition function becomes non-analytic. Since there are known phase
transitions at finite temperature (e.g. 2D classical Ising model), the
exponential decay does not hold for all thermal states at all temperatures. As
such, this can typically be thought of as a characteristic property of the set
of completely analytical interactions from Sec. I.2.

The decay of correlations is an important fact: it shows that the different
parts of the system behave almost completely independently. A state with this
property should then share many large-scale features with an uncorrelated gas,
in which the particles are not interacting at all. This has as a wealth of
related physical consequences. For instance, it is associated with basic
statistical physics facts covered in Sec. VI, in particular the validity of
the central limit theorem and related results on concentration properties of
thermal states Brandão and Cramer (2015); Anshu (2016) and the phenomenon of
equivalence of ensembles Brandão and Cramer (2015); Tasaki (2018); Kuwahara
and Saito (2020b). It also features in the proof of _local
indistinguishabiliy_ in Sec. V.1.1.

###  IV.3 A refined correlation decay: Conditional mutual information

In Sec. IV.1, we mentioned that the area law itself does not quite imply that
the correlations in a system are localized, in the sense that a particular
subsystem is only appreciably correlated with its vicinity. There is, however,
a significantly stronger statement about correlations that does imply it in a
clear way. 444The ideas described here can be seen as the quantum analogues of
a much stronger statement that holds for classical probability distributions:
the Hammersley-Clifford theorem Hammersley and Clifford (1971).

This is the property of being an approximate _quantum Markov state_ Hayden
_et al._ (2004a), which is defined in terms of the decay of the CMI in Eq.
(15). Let us consider three regions $A,B,C$ such that $B$ shields $A$ from
$C$. A simple example is given in Fig. 4, or in Fig. 5 for 1D.

.

Figure 4: In this configuration, the region $B$ shields $A$ from $C$, such
that the minimum distance between $A$ and $C$ is given by the shortest path
from $A$ to $C$ through $B$.

.

Figure 5: In this chain, the middle region $B$ shields $A$ from $C$, and their
distance is related to the size of $B$.

Since this quantifies how many of the correlations between $A$ and $C$ are not
mediated through $B$, we thus expect that it becomes small as the size of $B$
grows, and $A,C$ are further apart. This is perhaps the strongest sense in
which correlations can be localized.

This is studied in one dimensional systems in Kato and Brandão (2019). By
choosing $A,B,C$ to be adjacent regions of the chain (see Fig. 5), Kato and
Brandão (2019) shows that

$I(A:C|B)_{\rho_{\beta}}\leq c_{1}|B|e^{-c_{2}\sqrt{|B|}}.$ (68)

It is expected that the decay is $e^{-\Omega(|B|)}$ as opposed to Eq. (68),
which may be important for certain applications Kim (2017); Kato and Brandão
(2019); Brandão and Kastoryano (2019). The key technique is the quantum belief
propagation from Sec. III.3, but the locality bounds from Sec. III.2 are also
sufficient. The idea is to use those results to define a completely positive
map corresponding to a particular POVM outcome, that can be used in a “measure
until success” strategy. This result, however, relies on the exponential decay
of correlations in 1D, and also on bound on the correlation length of the form
$\xi\leq e^{\mathcal{O}(\beta)}$, which is currently not known. Interestingly,
Kato and Brandão (2019) also shows a converse statement: any state with a
sufficiently fast decaying CMI approximates the thermal state of some local
Hamiltonian.

In larger dimensions, the work Kuwahara  _et al._ (2020) shows that a non-
commutative analogue of the cluster expansion (in which there are no traces or
expectation values, but rather operators) suffices to study this problem. In
particular, if the cluster expansion corresponding to the object
$\log\text{Tr}_{\setminus A}[e^{-\beta H}]$ converges, exponentially, one
obtains an exponential decay of the form

$I(A:C|B)_{\rho_{\beta}}\leq k_{1}\min\\{|\partial A|,|\partial
C|\\}\left(\frac{\beta}{\beta^{*}_{\text{NC}}}\right)^{-k_{2}\times\text{dist}(A,C)},$
(69)

which only works for high enough temperatures
$\beta<\beta^{*}_{\text{NC}}=\mathcal{O}(1)$. This expansion is more involved
than the one described in Sec. III.1 in that the individual terms of the
expansion may not commute (since no trace is being taken when expanding
$\log\text{Tr}_{\setminus A}[e^{-\beta H}]$). Thus, the constant
$\beta^{*}_{\text{NC}}$ need not be the same as the $\beta^{*}$.

The significance of a fast decay of the CMI is highlighted by the idea of the
Petz map Petz (1986); Hayden  _et al._ (2004b). An important result in this
regard states that, given a tripartite state $\rho=\rho^{ABC}$, there exists a
CPTP map $\mathcal{N}(\cdot)_{B\rightarrow BC}$ (that is, acting on $B$, and
with output on BC) such that Fawzi and Renner (2015); Junge  _et al._ (2018)

$I(A:C|B)_{\rho}\geq 2||\rho^{ABC}-\mathcal{N}(\rho^{AB})_{B\rightarrow
BC}||_{1}.$ (70)

The map $\mathcal{N}$ usually goes under the name of _recovery map_. See
Sutter (2018) for an overview.

A fast decay of the CMI thus guarantees that the Gibbs state on $A,B,C$ can be
reconstructed from $\rho_{AB}$ by acting locally on $B$ (and importantly, not
on $A$), such that $\mathcal{I}_{A}\otimes\mathcal{R}_{B\rightarrow
BC}(\rho^{AB})\simeq\rho^{ABC}$, with $\mathcal{R}_{B\rightarrow BC}$ some CP
map taking only $B$ as input. This gives a way of sequentially preparing the
whole thermal state from its smaller components, which can potentially be used
e.g. for quantum algorithms (see Sec. VII).

##  V Locality of temperature

In the previous section we focused on the correlations between different
parts. Now, we move the spotlight to features of individual subsystems. That
is, if we divide the system into $A$ and its compliment $\setminus A$, what
does $\text{Tr}_{\setminus A}[\rho_{\beta}]$ look like? In the rest of the
section we drop the subscript $\beta$ for simplicity of notation.

Consider first the trivial case: if the particles are non-interacting, it
holds that the marginal on $A$ is the thermal state of $H_{A}$ which is the
Hamiltonian that acts only on subsystem $A$. That is

$\rho^{A}\equiv\text{Tr}_{\setminus A}[\rho]=\frac{e^{-\beta H_{A}}}{Z_{A}}.$
(71)

Now, how does Eq. (71) change when we introduce local (and potentially strong)
interactions? Can we identify the state of a subsystem with some thermal
state? How different is it from $\frac{e^{-\beta H_{A}}}{Z_{A}}$? This general
question sometimes goes under the name of _locality of temperature_ Kliesch
_et al._ (2014).

There are (to the author’s knowledge) two different but related answers to
this: the idea of _local indistinguishability_ and also the notion of
_Hamiltonian of mean force_. We now explain both of them, elaborate on their
significance for thermodynamics, and also give a proof of the simplest
instance of the first (in 1D).

###  V.1 Local indistinguishability

Given the above discussion on the decay of correlations, we expect that the
state of a local subsystem will not depend much on the parts that are far away
enough from it. A possible way to phrase this is that the local marginal
$\rho_{A}$ is indistinguishable from the marginal of a much smaller thermal
state, with a Hamiltonian that acts only in the vicinity of $A$. We now make
this intuition precise.

Let us refer to partitions into $ABC$ such as those in Fig. 4 or Fig. 5, and
write the Hamiltonian with the following terms:

$H=H_{A}+H_{AB}+H_{B}+H_{BC}+H_{C}.$ (72)

We now have the full thermal state $\rho$, as well as a thermal state
supported on $A,B$ defined as

$\rho_{0}^{AB}=\frac{e^{-\beta(H_{A}+H_{B}+H_{AB})}}{Z_{AB}},$ (73)

that is, without the terms in $H$ that have support in $C$. One can also think
of this as the marginal of the thermal state
$\rho_{0}^{AB}\otimes\rho_{0}^{C}\equiv e^{-\beta(H_{AB}+H_{C})}/Z_{AB}Z_{C}$
in which we have removed the interactions $H_{BC}$ between $AB$ and $C$.
Notice that $\rho_{0}^{AB}\neq\rho^{AB}$ due to the presence of $H_{BC}$. This
is, however, just a small local term.

The main idea is that if $B$ is large enough, these two states are almost
indistinguishable on $A$. Let us assume that the connected correlations from
Eq. (65) decay with function $f(\text{dist}(C,D))$. Then, the following upper
bound holds for some constant $K>0$ Brandão and Kastoryano (2019)

$\displaystyle||$
$\displaystyle\text{Tr}_{BC}[\rho]-\text{Tr}_{B}[\rho_{0}^{AB}]||_{1}$ (74)
$\displaystyle\leq
K|\partial_{C}|\left(f(\text{dist}(A,C))+e^{-\Omega(\text{dist}(A,C))}\right).$

The first term in the RHS comes from the decay of correlations assumption. The
second comes from using the QBP technique in Sec. III.3. The exponential decay
of this quantity thus holds whenever both the correlations decay fast enough,
and Lieb-Robinson bounds hold. An alternative proof for high temperatures
using the cluster expansion can also be found in Kliesch  _et al._ (2014).

A straightforward consequence of this is that we do not need to know the whole
state to compute local quantities. If we care about some kind of local order
parameter, or want to compute currents or else between some part and its
surroundings, we can calculate them without having to diagonalize a huge
matrix of size $\exp(N)$, but rather just focus on a much smaller region. This
is particularly useful in translation-invariant systems.

####  V.1.1 Proof in 1D

We now show the full proof of this result in the case of one dimension. The
more general one, however, is essentially the same and can be found in Brandão
and Kastoryano (2019). It uses previously mentioned results, and shares some
steps and ideas that appear in other fundamental questions including the proof
of the absence of phase transitions in 1D Araki (1969); Harrow  _et al._
(2020) or of decay of correlations Araki (1969); Bluhm  _et al._ (2022). It
will also be a key ingredient in the algorithm of Sec. VII.1.

We focus on the restricted setting of a chain, that we divide into three parts
$A,B,C$, such that $B$ is in the middle and $A$ is a small subsystem at the
end of the chain, as in Fig. 6.

.

Figure 6: Choice of regions for the proof in Sec. V.1.1, and depiction of the
distance $l$ on which the operator $E_{BC}^{l}$ acts within $B$.

The aim is a small upper bound on

$||\text{Tr}_{BC}[\rho]-\text{Tr}_{B}[\rho_{0}^{AB}]||_{1}=\max_{||M_{A}||\leq
1}|\text{Tr}[M_{A}(\rho-\rho_{0}^{AB}\otimes\rho_{0}^{C})]|,$ (75)

where $M_{A}$ has support on $A$ only, and the equality comes from the
definition of the $1$-norm. Now, let us define the following two operators

  * •

$E_{BC}=e^{\beta H}e^{-\beta(H-H_{BC})}$,

  * •

$E^{l}_{BC}=e^{\beta(H^{l}_{C}+H^{l}_{B}+H_{BC})}e^{-\beta(H^{l}_{B}+H^{l}_{C})}$,
where $H^{l}_{B}$ and $H^{l}_{C}$ are the terms of $H_{B}$ and $H_{C}$ that
are a distance at most $l$ from the boundary terms $H_{BC}$.

The parameter $l$ is free, so we can choose to our convenience. We refer now
to the result from Araki (1969) in Eq. (39) and (40), from which it follows
that

$\displaystyle||E_{BC}||$ $\displaystyle\leq C_{1}$ (76)
$\displaystyle||E_{BC}-E^{l}_{BC}||$ $\displaystyle\leq
C_{2}\frac{q^{1+l}}{(l+1)!}$ (77)

That is, the operator $E_{BC}$ has bounded norm and, since we can approximate
it by $E^{l}_{BC}$ with some $l<\text{dist}(A,C)$, its support on region $A$
is super-exponentially suppressed in $l$ (due to the factorial, which always
dominates over $q^{l}$). In what follows, we choose $l=|B|/2$. Notice that by
definition $\rho_{0}^{AB}\otimes\rho_{0}^{C}=\frac{Z}{Z_{AB}Z_{C}}\rho
E_{BC}$.

Let us define $M_{A}^{*}$ to be the operator that optimizes the RHS of Eq.
(75). With the triangle inequality we can write

$\displaystyle|\text{Tr}[M^{*}_{A}(\rho-\rho_{0}^{AB}\otimes\rho_{0}^{C})]|$
(78)
$\displaystyle\leq\left|\text{Tr}[M^{*}_{A}(\rho-\frac{Z}{Z_{AB}Z_{C}}\rho
E^{l}_{BC})]\right|$
$\displaystyle+\left|\text{Tr}[M^{*}_{A}(\frac{Z}{Z_{AB}Z_{C}}\rho
E^{l}_{BC}-\rho_{0}^{AB}\otimes\rho_{0}^{C})]\right|.$

Let us now upper-bound these two terms independently. The second can be
bounded with Eq. (77) and Hölder’s inequality applied twice.

$\displaystyle\left|\text{Tr}[M^{*}_{A}\frac{Z}{Z_{AB}Z_{C}}\rho
E^{l}_{BC}-\rho_{0}^{AB}\otimes\rho_{0}^{C})]\right|$ (79)
$\displaystyle=\left|\text{Tr}[M^{*}_{A}\frac{Z}{Z_{AB}Z_{C}}\rho(E^{l}_{BC}-E_{BC})]\right|$
(80)
$\displaystyle\leq\frac{Z}{Z_{AB}Z_{C}}||M^{*}_{A}||||\rho||_{1}||E^{l}_{BC}-E_{BC}||$
(81) $\displaystyle\leq\frac{Z}{Z_{AB}Z_{C}}\times
C_{2}\frac{q^{1+l}}{(1+l)!}.$ (82)

Given Eq. (53), $\max\\{\frac{Z}{Z_{AB}Z_{C}},\frac{Z_{AB}Z_{C}}{Z}\\}\leq
e^{\beta||H_{BC}||}$ , which is a constant that only depends on $\beta,k,J$.
Thus, this second term is super-exponentially suppressed in $|B|$.

For the first term, we require the decay of correlations property Eq. (65)
(which holds in 1D under the assumption of translation invariance). Since
$l=|B|/2$,

$\displaystyle\left|\text{Tr}[M^{*}_{A}\rho
E^{l}_{BC}]-\text{Tr}[M^{*}_{A}\rho]\text{Tr}[\rho E^{l}_{BC}]\right|$ (83)
$\displaystyle\leq Ke^{-\frac{|B|}{2\xi}}||E^{l}_{BC}||\leq
2KC_{1}e^{-\frac{|B|}{2\xi}},$ (84)

where for the last inequality we used
$||E^{l}_{BC}||\leq||E_{BC}||+||E^{l}_{BC}-E_{BC}||\leq 2C_{1}$, which holds
for sufficiently large $l$. We can now write

$\displaystyle\left|\text{Tr}[M^{*}_{A}(\rho-\frac{Z}{Z_{AB}Z_{C}}\rho)E^{l}_{BC}]\right|$
(85)
$\displaystyle\leq\left|\text{Tr}[M^{*}_{A}\rho]-\frac{Z}{Z_{AB}Z_{C}}\text{Tr}[M^{*}_{A}\rho]\text{Tr}[\rho
E^{l}_{BC}]\right|+2KC_{1}e^{-\frac{|B|}{2\xi}}$
$\displaystyle\leq\left(1-\frac{Z}{Z_{AB}Z_{C}}\text{Tr}[\rho
E^{l}_{BC}]\right)+2KC_{1}e^{-\frac{|B|}{2\xi}},$

where we used the triangle inequality in the first line, and Hölder’s
inequality $\text{Tr}[M^{*}_{A}\rho]\leq||M^{*}_{A}||\leq 1$ to get to the
second. Finally, we can use Eq. (77) again after another application of
Hölder’s inequality

$\displaystyle|\text{Tr}[\rho E^{l}_{BC}]-\text{Tr}[\rho
E_{BC}]|\leq||E_{BC}-E^{l}_{BC}||\leq C_{2}\frac{q^{1+l}}{(1+l)!},$ (86)

and since $\text{Tr}[\rho E_{BC}]=\frac{Z}{Z_{AB}Z_{C}}\leq
e^{\beta||H_{BC}||}$ we obtain

$\displaystyle\left|\text{Tr}[M^{*}_{A}(\rho-\frac{Z}{Z_{AB}Z_{C}}\rho
E^{l}_{BC})]\right|$ (87) $\displaystyle\leq
C_{2}e^{\beta||H_{BC}||}\frac{q^{1+l}}{(1+l)!}+2KC_{1}e^{-\frac{|B|}{2\xi}}.$
(88)

This finishes the proof. Putting everything together, we see that we have
upper-bounded our target quantity in Eq. (75) by a small number related to the
error term in the decay of correlations and Araki’s result. Without writing
the constants explicitly, and just on the leading exponential error, the final
result is

$||\text{Tr}_{BC}[\rho]-\text{Tr}_{B}[\rho_{0}^{AB}]||_{1}\leq
e^{-\Omega(|B|)},$ (89)

where $\Omega(x)$ is defined in Sec. II.4.

For simplicity, we have only dealt with the case of a 1D chain, where $A$ is
at the end of it. To generalize the proof, one just needs to define an
analogous partition $ABC$ in higher dimensions (see Fig. 4 for 2D) and then
remove all the different interaction terms from $H_{BC}$ one by one. Here, we
have done it with the operator $E_{BC}$, but this can also be done with the
(suitably defined) QBP operator $O_{A}$ from Sec. III.3, and the result Eq.
(74) is essentially unchanged.

###  V.2 Hamiltonian of mean force

The state $\rho_{A}$ is obviously the thermal state of _some_ Hamiltonian on
$A$, since we can always define

$\tilde{H}_{A}\equiv\beta^{-1}\log\text{Tr}_{\setminus A}[e^{-\beta H}],$ (90)

which is in general different from $H_{A}$. This is the so-called Hamiltonian
of mean force Miller (2018). Notice that it can be defined up to some additive
constant chosen at will.

How does this Hamiltonian compare to the “bare” Hamiltonian $H_{A}$, which
disregards the interactions between $A$ and the rest? That is, we would like
to understand the norm and locality of the operator
$\Phi_{A}\equiv\tilde{H}_{A}-H_{A}$. This turns out to be a difficult problem,
very much related to both the decay of mutual information and of the
conditional mutual information from Sec. IV.

One potential result is as follows. Since the interactions are local, it makes
sense that, if the size of $A$ is much larger than the number of nearest
neighbours $k$, most of the weight of $\Phi_{A}$ is localized around its
boundary with the rest of the system, of size $|\partial A|$.

If the size of $A$ is much larger than the number of nearest neightbours $k$,
most of the weight of $\Phi_{A}$ is localized around its boundary with the
rest of the system, of size $|\partial A|$. The precise question is: can we
approximate $\Phi_{A}$ with another operator $\Phi_{A}^{l}$ that only has
support on sites a distance $l$ away from the boundary? Using a non-
commutative cluster expansion, as in Sec. IV.3, Theorem 2 in Kuwahara  _et
al._ (2020) shows that, for any temperature $\beta$ above a threshold one
$\beta^{*}_{\text{NC}}>0$, one can define a $\Phi_{A}^{l}$ such that

$||\Phi_{A}-\Phi_{A}^{l}||\leq\frac{e}{2\beta}(\beta/\beta^{*}_{\text{NC}})^{l/k}|\partial
A|.$ (91)

That is, $\Phi_{A}$ can be exponentially well approximated with an operator
localized around the boundary. See Fig. 7 for an illustration.

.

Figure 7: Illustration of the regions of the Hamiltonian of mean force. The
correction term $\Phi_{A}$ is exponentially well approximated by
$\Phi_{A}^{l}$, which has support on the region $\partial A$ only.

A similar result is expected to hold in 1D, but this is a so far open problem.
See Trushechkin  _et al._ (2022) for a recent overview on this topic for a
different set of models, and its implications.

###  V.3 Non-equilibrium thermodynamics with strong coupling

The ideas of this section may help understand how a priori complex non-
equilibrium thermodynamic processes may be tractable in practice. In many
situations of interest, the starting point at $t=0$ is an equilibrium state
with time-dependent Hamiltonian

$H(t)=H_{S}(t)+H_{B}+H_{I}.$ (92)

This includes a system Hamiltonian $H_{S}(t)$, a bath Hamiltonian $H_{B}$, and
an interaction $H_{I}$ between them 555This interaction might also depend on
time, which here we do not consider for simplicity.. The driving of $H_{S}(t)$
for $t>0$ takes the system away from equilibrium, and different thermodynamic
quantities can then be studied.

Textbook thermodynamics are typically centered around macroscopic systems such
as gases, where one can take the weak coupling limit $H_{I}\ll H_{S},H_{B}$.
In many body quantum systems, such as the ones considered here, this limit may
no longer apply, which creates a number of difficulties for thermodynamic
considerations (see e.g. Talkner and Hänggi (2020); Miller (2018); Strasberg
and Esposito (2020)). These, however, can be dealt with if one considers the
effective state on the system

$\rho^{S}(t)=\frac{\text{Tr}_{B}[e^{-\beta H(t)}]}{\text{Tr}_{SB}[e^{-\beta
H(t)}]}\equiv\frac{e^{-\beta\tilde{H}_{S}(t)}}{\tilde{Z}_{S}(t)},$ (93)

where it is convenient to define the Hamiltonian of mean force
$\tilde{H}_{S}(t)$ such that
$\tilde{Z}_{S}(t)\equiv\frac{\text{Tr}_{SB}[e^{-\beta
H(t)}]}{\text{Tr}_{B}[e^{-\beta H_{B}}]}$. This way, for instance, the non-
equilibrium free energy (and thus the second law) is defined as

$F_{S}(t)=\text{Tr}_{S}[\tilde{H}_{S}(t)\rho^{S}(t)]-\frac{1}{\beta}S(\rho^{S}(t)).$
(94)

See e.g. Strasberg and Esposito (2020) for details. These quantities may be
difficult to calculate. While $\tilde{H}_{S}(t)$ may be inferred from the
system alone, $\tilde{Z}_{S}(t)$ may depend on the system but also potentially
on its relation to the whole bath.

Consider the rather general situation in which the system to be a small part
of a lattice Hamiltonian of a Gibbs state with short-ranged correlations, and
the bath to be the rest of the lattice. In that case, the situation simplifies
dramatically. First, the discussion from Sec. V.2 suggests that typically the
Hamiltonian of mean force $\tilde{H}_{S}$ does not necessarily depend on the
whole bath, but has some corrections which only depend on the region around
the boundary between $S$ and $B$.

Then, for the effective partition function $\tilde{Z}_{S}$, the same fact
follows from local indistinguishability. To show this, notice that

$\displaystyle\frac{\tilde{Z}_{S}(t)}{Z_{S}(t)}\equiv$
$\displaystyle\frac{\text{Tr}_{SB}[e^{-\beta H(t)}]}{\text{Tr}_{SB}[e^{-\beta
H_{S}(t)}\otimes e^{-\beta H_{B}}]}$ (95) $\displaystyle=$
$\displaystyle\text{Tr}_{SB}[O^{\dagger}_{I}O_{I}\frac{e^{-\beta
H_{S}(t)}\otimes e^{-\beta H_{B}}}{\text{Tr}_{SB}[e^{-\beta H_{S}(t)}\otimes
e^{-\beta H_{B}}]}]$ $\displaystyle=$
$\displaystyle\text{Tr}_{SB}[(O^{l}_{I})^{\dagger}O^{l}_{I}\frac{e^{-\beta
H_{S}(t)}\otimes e^{-\beta H_{B}}}{\text{Tr}_{SB}[e^{-\beta H_{S}(t)}\otimes
e^{-\beta H_{B}}]}]$ $\displaystyle+e^{\mathcal{O}(\beta||H_{I}(t)||-l)}$ (96)

where $O_{I}$ is the belief propagation operator from Eq. (48) with $A=H_{I}$
and $O^{l}_{I}$ is its local approximation in Eq. (50). Defining $B_{l}$ to be
the region of the bath that is at most a distance $2l$ away from $S$, we have
that

$\displaystyle\frac{\tilde{Z}_{S}(t)}{Z_{S}(t)}=$
$\displaystyle\text{Tr}_{SB_{l}}[(O^{l}_{S})^{\dagger}O^{l}_{S}\frac{e^{-\beta
H_{S}(t)}\otimes e^{-\beta H_{B_{l}}}}{\text{Tr}_{SB_{l}}[e^{-\beta
H_{S}(t)}\otimes e^{-\beta H_{B_{l}}}]}]$
$\displaystyle+e^{\mathcal{O}(\beta||H_{I}||-l)}.$ (97)

This shows that the effective partition function can be approximated to
multiplicative error $\epsilon$ by computing the expectation value of
$O^{\dagger}_{I}O_{I}$ on the system and a region of the bath a distance
$l=\log(\epsilon^{-1})+\mathcal{O}(\beta||H_{I}||)$ away from the small
system. In $D$ dimensions, assuming $||H_{I}||=\mathcal{O}(1)$, the
computational cost of exact diagonalization is
$e^{\mathcal{O}(l)}=e^{\mathcal{O}(\log^{D}\epsilon^{-1})}$, independent of
the size of the bath.

##  VI Statistical properties

We now explain and prove some important statistical features of thermal
states. These are central statements of the field of statistical physics and
characterize the ensembles involved: the thermal (or canonical) and the
microcanonical, as well as the grand canonical or others, when relevant. In
contrast to the results of other sections, all those shown here (as well as
their proofs) apply equally to classical models.

###  VI.1 Measurement statistics and concentration bounds

In Sec. IV we saw how in many instances of thermal states, in particular for
those in the “one phase region”, the different subsystems tend not to have
strong correlations. This has a number of consequences, and we now explore an
important one that shows that their large-scale statistical properties
resemble those of non-interacting/statistically independent systems. These are
_concentration bounds_ , akin to the (perhaps more widely known) central limit
theorem.

The setting is as follows: let us consider a $k$-local observable
$A=\sum_{j}A_{j}$, such that $A_{j}$ has support on at most $k$ sites. The
best example is the energy, but also other properties like magnetization
$\sum_{j}\sigma_{j}^{Z}$.

The expectation value of any such observable can be thought of as a
macroscopic property of the system (such as the average magnetization of the
material). While we expect that there will be thermal fluctuations around that
average value, our intuition from thermodynamics tells us that any such large-
scale property should have a definite value, almost free of fluctuations. This
is due to one of the most basic ideas from probability theory: the measurement
statistics of sums of independent random variables greatly concentrate around
the average. The main conclusion is that if we measure an observable $A$ on a
thermal state, the outcome will be very close to the average $\langle
A\rangle_{\beta}$ with overwhelmingly large probability. That is, the
distribution

$P_{A,\beta}(x)=\text{Tr}[\rho\delta(x-A)].$ (98)

which is the probability of obtaining outcome $x$ when measuring $A$, is
highly peaked around the average $\langle
A\rangle_{\beta}=\text{Tr}[\rho_{\beta}A]$. This has important implications
for the validity of thermodynamic descriptions of these systems, in that
averaged macroscopic quantities characterize the large system of many
particles whose properties we do not know with any certainty.

In the theory of probability, there are various types of results
characterizing distributions comprised of many independent (or close to
independent) variables. Their proofs most often involve constraining the
characteristic function $\langle e^{\lambda A}\rangle_{\beta}$, where
$\lambda$ may be real or imaginary. We now describe some of them.

####  VI.1.1 Chernoff-Hoeffding bound

This is a concentration bound that reads

$\displaystyle P_{A,\beta}(|x-\langle A\rangle_{\beta}|>\delta)\leq
2\exp\left(-\frac{\delta^{2}}{4c\bar{A}}\right),$ (99)

where $\bar{A}\equiv\sum_{j}||A_{j}||$. Thus if $\delta^{2}\gg c\bar{A}$, the
probability of measuring $A$ to be away from $\langle A\rangle_{\beta}$ by at
least $\delta$ is exponentially small. The most common proof technique is via
a bound on the characteristic function of the form

$\log\langle e^{\tau(A-\langle A\rangle_{\beta})}\rangle_{\beta}\leq
c\tau^{2}\bar{A},$ (100)

for some $\mathcal{O}(1)$ constant $c>0$ (which might depend on $\beta$ and
the parameters of the Hamiltonian) and a wide enough range of $\tau$. From
this it follows that

$\displaystyle P_{A,\beta}(x-\langle A\rangle_{\beta}>\delta)=\int_{x-\langle
A\rangle_{\beta}>\delta}\text{Tr}[\rho\delta(x-A)]$ (101)
$\displaystyle=\int_{x-\langle A\rangle_{\beta}>\delta}\text{Tr}[\rho
e^{\tau(A-\langle A\rangle_{\beta})}e^{-\tau(A-\langle
A\rangle_{\beta})}\delta(x-A)]$ (102) $\displaystyle\leq
e^{-\tau\delta}\text{Tr}[\rho_{\beta}e^{\tau(A-\langle
A\rangle_{\beta})}]\leq\exp\left(-\tau\delta+c\tau^{2}\bar{A}\right).$ (103)

One can follow the same steps for the range $\langle
A\rangle_{\beta}-x>\delta$. Then, choosing $\tau=\delta/(2c\bar{A})$ yields
Eq. (99).

This result can be very easily shown for independent random variables or for
independent spins. For interacting spins, Eq. (100) was shown in Kuwahara and
Saito (2020a) with the cluster expansion technique from Sec. III.1, which
holds for all dimensions and all temperatures $\beta<\beta^{*}$. To see this,
notice that a bound of the form of Eq. (100) follows from proving the
convergence of the expansion of $\log\langle e^{\tau(A-\langle
A\rangle_{\beta})}\rangle_{\beta}$ to second order. The main result from Anshu
(2016) proves a slightly weaker version of Eq. (99) with a different
technique, only assuming the decay of correlations from Sec. IV.2.

####  VI.1.2 Large deviation bound

A related important type of concentration bound is given by large deviation
theory. This is the branch of probability theory concerned with understanding
the likelihood of very rare events, and has a long history as one of the most
important mathematical frameworks for studying statistical physics. For
instance it gives a way of describing the equilibrium properties of large
ensembles (as is also the case here), or for predicting the long-time
behaviour of non-equilibrium processes such as Brownian motion. See Touchette
(2009) for an excellent overview of the main results and their consequences
for classical systems.

The basic idea is that given any set of measurement outcomes $\mathcal{A}$, we
would like to identify whether there always exists a _rate function_
$I_{\mathcal{A}}$ such that

$\lim_{N\rightarrow\infty}-\frac{\log
P_{A,\beta}(x\in\mathcal{A})}{N}=I_{\mathcal{A}}.$ (104)

If this is the case, the dominant behaviour of $P_{A,\beta}(x\in\mathcal{A})$
is essentially a decaying exponential $P_{A,\beta}(x\in\mathcal{A})\simeq
e^{-NI_{\mathcal{A}}+o(N)}$, unless $I_{\mathcal{A}}=0$. This means that, in
the thermodynamic limit, the measurement statistics of $A$ are extremely
peaked around the points where the rate function vanishes $I_{\mathcal{A}}=0$.

This is slightly stronger than the Chernoff-Hoeffding inequality, in that it
can in principle give an exact expression of the probability distribution for
large enough $N$. However, we do not always know how large an $N$ is “enough”,
and for finite $N$, it often does not give an expression as explicit as Eq.
(99).

Again, the proof strategy most often involves the characteristic function. In
particular, the Gärtner-Ellis theorem states that a sufficient condition is
that the function

$g(\tau)=\lim_{N\rightarrow\infty}\frac{\log\langle e^{\tau
A}\rangle_{\beta}}{N}$ (105)

exists and is differentiable. This has been shown using the cluster expansion
in Netočný and Redig (2004) for $1$-local observables, and upper bounds on the
rate for general observables have been shown using the locality estimates from
Sec. III.2 in Lenci and Rey-Bellet (2005); Hiai  _et al._ (2007). The full
large deviation principle was shown in 1D in Ogata (2010). An alternative
proof can be found in Ogata and Rey-Bellet (2011).

####  VI.1.3 Berry-Esseen theorem

Another interesting probability theory result is the Berry-Esseen theorem
Berry (1941); Esseen (1942); Brandão and Cramer (2015), which can be thought
of as a refinement of the central limit theorem for a finite sample size
(which in this case is the system size $N$). Let us define the cumulative
distribution function

$\displaystyle F(x)=\int_{-\infty}^{x}P_{A,\beta}(x)\text{d}x$ (106)

as well as the equivalent for a Gaussian with the same average and variance

$\displaystyle
G(x)=\frac{1}{\sigma_{A}\sqrt{2\pi}}\int_{-\infty}^{x}e^{\frac{-(y-\langle
A\rangle_{\beta}}{2\sigma^{2}_{A}}},$ (107)

where $\sigma^{2}_{A}=\langle A^{2}\rangle_{\beta}-\langle
A\rangle_{\beta}^{2}$. The distance between these two functions is bounded by
Esseen’s inequality Feller (1991), which states that, for all $T>0$

$\displaystyle\Delta$ $\displaystyle\equiv\max_{x}|F(x)-G(x)|$ (108)
$\displaystyle\leq\frac{18}{\sqrt{2\pi^{3}}T}+\frac{1}{\pi}\int_{0}^{T}\frac{\left|e^{-\frac{t^{2}}{2}}-\langle
e^{it\frac{A-\langle
A\rangle_{\beta}}{\sigma_{A}}}\rangle\right|}{t}\text{d}t.$ (109)

That is, the right hand side is small if the characteristic function inside
the integral is close to a Gaussian for rather long times $t$.

This can be shown by bounding the logarithm of the characteristic function
$\langle e^{it\frac{A-\langle
A\rangle_{\beta}}{\sigma_{A}}}\rangle=\text{Tr}[e^{it\frac{A-\langle
A\rangle_{\beta}}{\sigma_{A}}}\frac{e^{-\beta H}}{Z}]$ with the cluster
expansion from Sec. III.1. Assuming $h=\mathcal{O}(1)$, for short times
$t/\sigma_{A}\leq t^{*}$, with $t^{*}$ some $\mathcal{O}(1)$ constant, it can
be shown that it is close to the second order Taylor expansion as

$\displaystyle\left|\log\langle e^{it\frac{A-\langle
A\rangle_{\beta}}{\sigma_{A}}}\rangle-\frac{-t^{2}}{2}\right|\leq\mathcal{O}\left(\frac{Nt^{3}}{\sigma^{3}}\right),$
(110)

so that

$\displaystyle\langle e^{it\frac{A-\langle
A\rangle_{\beta}}{\sigma_{A}}}\rangle=e^{-\frac{t^{2}}{2}+\mathcal{O}\left(\frac{Nt^{3}}{\sigma^{3}}\right)}.$
(111)

To prove this, see for instance Theorem 13 in Wild and Alhambra (2023). This
allows us to bound the integral in Eq. (109) choosing
$T=\frac{t^{*}\sigma_{A}}{2}$, to achieve

$\Delta\leq\mathcal{O}\left(\frac{1}{\sigma_{A}}+\frac{N}{\sigma_{A}^{3}}\right),$
(112)

which, considering that $\sigma_{A}=\Omega(\sqrt{N})$, means that
$\Delta=\mathcal{O}(N^{-1/2})$. This means that the cumulative functions
$F(x)$ and $G(x)$ become increasingly similar with system size, which shows
that the probability $P_{A,\beta}(x)$ approaches a Gaussian in the
thermodynamic limit.

A different proof starting from the assumption of decay of correlations, can
be found in Brandão  _et al._ (2015).

###  VI.2 Equivalence of ensembles

We now prove an important statement in the study of statistical physics, which
goes back all the way to Boltzmann and Gibbs. In large systems, the average
macroscopic properties of both the thermal or canonical state, and of the
microcanonical ensemble, are essentially the same. This means that both
canonical and ergodic averages coincide in the thermodynamic limit, and shows
that the particular ensemble used for calculations does not necessarily
matter.

There are various similar statements in the literature Lima (1971, 1972);
Müller  _et al._ (2015); Brandão and Cramer (2015); Tasaki (2018); Kuwahara
and Saito (2020b, a), but the proof that we now show follows that of Tasaki
(2018); Kuwahara and Saito (2020b, a) and relies on the Chernoff-Höffding
bound. Let us define the extensive observable $A=\sum_{j}A_{j}$ (such as e.g.
the total magnetization $\sum_{j}^{N}\sigma_{j}^{Z}$) with thermal/canonical
average $\langle A\rangle_{\beta}$ which for simplicity we will set to
$\langle A\rangle_{\beta}=0$, while the microcanonical average is

$\langle
A\rangle_{E,\Delta}=\frac{1}{D_{N}(E,\Delta)}\sum_{E_{j}\in(E-\Delta,E)}\langle
E_{j}|A|E_{j}\rangle,$ (113)

where $E$ is the energy and $\Delta$ the width of the microcanonical window
(which might depend on $N$), and $|E_{j}\rangle$ is the energy eigenstate of
energy $E_{j}$. $D_{N}(E,\Delta)$ is a normalization constant counting the
number of eigenstates within the window. This motivates the following
probability distribution

$\displaystyle
P_{E,\Delta}(x)=\frac{1}{D_{N}(E,\Delta)}\sum_{E_{j}\in(E-\Delta,E)}\delta(x-\langle
E_{j}|A|E_{j}\rangle),$ (114)

which gives the probability of measuring $x=\langle E_{j}|A|E_{j}\rangle$ when
sampling eigenstates from the microcanonical ensemble.

First, we need to determine what is the energy that corresponds to temperature
$\beta$ and thus characterizes the microcanonical ensemble. Given the
temperature $\beta$, the microcanonical energy $E_{0}$ is such that

$E_{0}(\Delta,\beta)\equiv\text{argmax}_{E}D_{N}(E,\Delta)e^{-\beta E}.$ (115)

Assuming that the width is significantly different than the energy scales of
the system, $\Delta\ll\langle H\rangle_{\beta}$ (which is most typically the
case), this roughly implies that $E_{0}$ is the energy of the microstates
$\\{|E_{j}\rangle\\}$ that have the dominant weight in the canonical ensemble
(when the density of states is weighted by the factor $e^{-\beta E}$). We have
written the dependence on $\beta,\Delta$ explicitly in Eq. (115), but let us
now drop them for simplicity of notation.

We start by upper bounding the $m$-th (even) moment of $P_{E_{0},\Delta}(x)$

$\displaystyle\int_{-\infty}^{\infty}x^{m}P_{E_{0},\Delta}(x)$ (116)
$\displaystyle=\frac{1}{D_{N}(E_{0},\Delta)}\sum_{E_{j}\in(E_{0}-\Delta,E_{0})}|\langle
E_{j}|A|E_{j}\rangle|^{m}$ $\displaystyle\leq$
$\displaystyle\frac{1}{D_{N}(E_{0},\Delta)}\sum_{E_{j}\in(E_{0}-\Delta,E_{0})}|\langle
E_{j}|A^{m}|E_{j}\rangle|=\langle A^{m}\rangle_{E_{0},\Delta},$

where we we used the convexity of $x^{m}$ with $m$ even. The bound can easily
be expressed in terms of a canonical average as, since $A^{m}$ is positive,

$\displaystyle\langle
A^{m}\rangle_{E_{0},\Delta}=\frac{1}{D_{N}(E_{0},\Delta)}\sum_{E_{j}\in(E_{0}-\Delta,E_{0})}\langle
E_{j}|A^{m}|E_{j}\rangle$ (117) $\displaystyle\leq\frac{e^{\beta
E_{0}}}{D_{N}(E_{0},\Delta)}\sum_{E_{j}\in(E_{0}-\Delta,E_{0})}e^{-\beta
E_{j}}\langle E_{j}|A^{m}|E_{j}\rangle$ (118) $\displaystyle\leq\frac{e^{\beta
E_{0}}}{D_{N}(E_{0},\Delta)}\sum_{E_{j}\in(-\infty,\infty)}e^{-\beta
E_{j}}\langle E_{j}|A^{m}|E_{j}\rangle$ (119) $\displaystyle=\frac{Ze^{\beta
E_{0}}}{D_{N}(E_{0},\Delta)}\langle A^{m}\rangle_{\beta}.$ (120)

The factor $\frac{Ze^{\beta E_{0}}}{D_{N}(E_{0},\Delta)}$ can now be upper
bounded using the definition of the microcanonical ensemble and the
concentration bound. Let us define the following modified partition function
$\tilde{Z}\equiv\sum_{|E_{j}-E_{0}|\leq\delta}e^{-\beta E_{j}}$. If we also
set $\delta=KN^{1/2}$ with $K=\mathcal{O}\left(1\right)$ it follows from Eq.
(99) that

$\displaystyle\frac{\tilde{Z}}{Z}$ $\displaystyle=1-P_{H,\beta}(|x-\langle
H\rangle_{\beta}|\geq\delta)$ (121) $\displaystyle\geq
1-2\exp\left(-\frac{\delta^{2}}{4cJN}\right)\geq 1/2.$ (122)

Now divide the energy range in the sum in equal parts of width
$\Delta^{*}\equiv\min\\{\Delta,\beta^{-1}\\}$, such that the largest energy of
each interval is $E_{\nu}$, so that $E_{\nu+1}=E_{\nu}+\Delta^{*}$ and

$\displaystyle\tilde{Z}$
$\displaystyle\leq\sum_{\begin{subarray}{c}\nu\in\mathbb{Z}\\\
|E_{\nu}-E_{0}|\leq\Delta^{*}+\delta\end{subarray}}D_{N}(E_{\nu},\Delta^{*})e^{-\beta(E_{\nu}-\Delta^{*})}$
(123) $\displaystyle\leq
e^{\beta\Delta^{*}}\left(\frac{2\delta}{\Delta^{*}}+2\right)\max_{\nu}D_{N}(E_{\nu},\Delta^{*})e^{-\beta
E_{\nu}}$ (124)
$\displaystyle\leq\frac{1}{2}K^{\prime}\frac{N^{1/2}}{\Delta^{*}}D_{N}(E_{0},\Delta)e^{-\beta
E_{0}},$ (125)

with $K^{\prime}=\mathcal{O}(1)$, where the last inequality follows from the
fact that $D_{N}(E_{0},\Delta)$ is monotonic on $\Delta$. We thus have
$\int_{-\infty}^{\infty}x^{m}P_{E_{0},\Delta}(x)\leq
K^{\prime}\frac{N^{1/2}}{\Delta^{*}}\langle A^{m}\rangle_{\beta}$. To finish
this part of the proof we bound $\langle A^{m}\rangle_{\beta}$. It was shown
in Kuwahara and Saito (2020b, a) that the concentration inequality Eq. (99)
implies that

$\langle
A^{m}\rangle_{\beta}\leq\left(4c\bar{A}\right)^{m/2}\left(\frac{m}{2}\right)!.$
(126)

For completeness, we reproduce the proof in Appendix B. We are now in a
position to bound the tail of $P_{E,\Delta}(x)$ as

$\displaystyle P_{E,\Delta}(x\geq
x_{0})=\int_{x_{0}}^{\infty}P_{E,\Delta}(x)\text{d}x\leq\frac{1}{x_{0}^{m}}\int_{-\infty}^{\infty}x^{m}P_{E_{0},\Delta}(x)$
$\displaystyle\leq
K^{\prime}\frac{N^{1/2}}{\Delta^{*}}\left(\frac{4c\bar{A}}{x_{0}^{2}}\right)^{m/2}\left(\frac{m}{2}\right)!\leq
K^{\prime}\frac{N^{1/2}}{\Delta^{*}}\left(\frac{4mc\bar{A}}{x_{0}^{2}}\right)^{m/2}.$
(127)

Thus, choosing $m=\lfloor{\frac{x_{0}^{2}}{4ce\bar{A}}}\rfloor$. and repeating
for $x\leq-x_{0}$, leads to (let us now bring back the average $\langle
A\rangle_{\beta}$ explicitly, previously taken to be zero)

$P_{E,\Delta}(|x-\langle A\rangle_{\beta}|\geq x_{0})\leq
2eK^{\prime}\frac{N^{1/2}}{\Delta^{*}}\exp\left(-\frac{x_{0}^{2}}{8ce\bar{A}}\right).$
(128)

We are almost done. We now bound the difference between canonical and
microcanonical as

$\displaystyle|\langle A\rangle_{E,\Delta}-\langle
A\rangle_{\beta}|\leq\sum_{j}\frac{|\langle E_{j}|A|E_{j}\rangle-\langle
A\rangle_{\beta}|}{D_{E_{0},\Delta}}$ (129)
$\displaystyle\leq\int_{|x|\leq\bar{A}}P_{E,\Delta}(x-\langle
A\rangle_{\beta})(x-\langle A\rangle_{\beta})\text{d}x$ (130)
$\displaystyle\leq x_{0}+2\bar{A}P_{E,\Delta}(|x-\langle A\rangle_{\beta}|\geq
x_{0}),$ (131)

and so choosing
$x_{0}=\sqrt{8ce\bar{A}}\log(4\bar{A}eK^{\prime}\frac{N^{1/2}}{\Delta^{*}})$,
the fact that $\bar{A}\propto N$ yields, for some constant $K^{\prime\prime}$,

$\frac{1}{N}|\langle A\rangle_{E,\Delta}-\langle
A\rangle_{\beta}|\leq\frac{K^{\prime\prime}\log\frac{N^{3/2}}{\Delta^{*}}}{N^{1/2}},$
(132)

so that the difference vanishes in the thermodynamic limit. Notice that
$\Delta^{*}\equiv\min\\{\Delta,\beta^{-1}\\}$, so that in principle even
rather low temperatures and very small (up to exponentially small)
microcanonical windows are allowed. This is the final result. It states that
average properties are essentially the same, provided that the average energy
$E_{0}$ is determined by Eq. (115), and that the width $\Delta$ is not too
small. The fact that it can be up to exponentially small in system size is
rather strong, and related to weak statements of the eigenstate thermalization
hypothesis (see Mori (2016); Kuwahara and Saito (2020b)).

##  VII Algorithms and complexity of thermal states

When addressing specific problems in many-body physics, we would most often
like to understand whether they are fundamentally complex or not, in the
precise sense established by theoretical computer science. This can typically
done in two complementary ways:

  * •

By showing that there exists an algorithm with a provable performance and run-
time. Additionally, it is interesting if the algorithm can be explicitly
constructed, and implemented in practice.

  * •

By establishing that a problem, or a set of them, belong to or are complete or
hard for a certain complexity class.

This applies to both classical and quantum computation, and their respective
complexity classes.

Problems related to quantum thermal states can also be studied under this
light. The relevant ones include most notably the estimation of the partition
function, or the generation of either approximations to the thermal states (in
quantum computers) or their classical representations (in classical
computers).

As an illustrative example of what can be proven, we start with a simple
explicit algorithm that approximates the quantum partition functions in 1D in
polynomial time Kuwahara and Saito (2018). We then briefly review some other
important known results about the hardness of approximating partition
functions. The rest of the section includes an explanation of the current best
tensor network results, which are provably efficient in a wide range of
situations, and a short review of quantum algorithms for preparing thermal
states.

###  VII.1 An efficient classical algorithm for the 1D partition function

Using some of the results from the previous sections, we now show that,
assuming that $h,\beta=\mathcal{O}(1)$, and that exponential decay of
correlations holds, we can efficiently approximate the partition function in
1D. This is done with an algorithm with runtime
$\text{poly}(N,\varepsilon^{-1})$ that outputs $Z^{\prime}$, where

$|\log Z^{\prime}-\log Z|\leq\mathcal{O}\left(\varepsilon\right).$ (133)

This section follows the result and proof strategy from Kuwahara and Saito
(2018), with some minor modifications.

In one dimension, let us consider the partial Hamiltonian
$H_{j}=\sum_{i=1}^{j-1}h_{i}$, which includes the first $j-1$ interaction
terms as counted from the left, starting from the left-most $h_{1}$. Then,
define the partial partition function

$\displaystyle Z_{i}$ $\displaystyle=\text{Tr}[e^{-\beta(H_{i}+h_{i})}]$ (134)
$\displaystyle=\text{Tr}[O_{h_{i}}e^{-\beta(H_{i})}O_{h_{i}}^{\dagger}]\equiv\text{Tr}[e^{-\beta(H_{i})}A_{i}],$
(135)

where $O_{h_{i}}$ is the quantum belief propagation from Sec. III.3 and
$A_{i}=O_{h_{i}}^{\dagger}O_{h_{i}}$. Now, rewriting Eq. (134) notice the
simple iterative relation

$Z_{i}=Z_{i-1}\text{Tr}[\rho_{i}A_{i}],$ (136)

where $\rho_{i}=e^{-\beta H_{i}}/Z_{i-1}$. Thus we can write

$Z=d^{N}\prod_{i=1}^{|E|}\text{Tr}[\rho_{i}A_{i}],$ (137)

where $Z\equiv Z_{|E|}$ and $d^{N}=Z_{0}$. The key now is to use results from
Sec. III.3 to approximate $A_{i}$, and local indistinguishability from Sec.
V.1.1. Let $A_{i}^{l}\equiv(O^{l}_{h_{i}})^{\dagger}O^{l}_{h_{i}}$, so that

$\displaystyle||A_{i}-A_{i}^{l}||$
$\displaystyle=||A_{i}-O_{h_{i}}^{\dagger}O_{h_{i}}^{l}+O_{h_{i}}^{\dagger}O_{h_{i}}^{l}-A_{i}^{l}||$
(138) $\displaystyle\leq 2||O_{h_{i}}|||O_{h_{i}}-O_{h_{i}}^{l}||$ (139)
$\displaystyle\leq e^{\mathcal{O}(\beta h)}e^{-\Omega(l)},$ (140)

where in the first line we used the triangle inequality and in the second we
used both Eq. (49) and (51). Now, let us label by $\Lambda_{l^{*}}$ to be the
rightmost part of the chain of length $l^{*}$, with vertex set $V_{l^{*}}$ and
in which $H_{i+1}$ has support. Choose $l^{*}\in\mathbb{R}$ so that
$A_{i}^{l^{*}}$ has support in the right side of $V_{2l^{*}}$ and define
$\rho_{i}^{(l_{*})}=e^{-\beta H_{\Lambda_{l_{*}}}}/\text{Tr}{e^{-\beta
H_{\Lambda_{l_{*}}}}}$, where $H_{\Lambda_{l_{*}}}=\sum_{\text{supp}(h_{i})\in
V_{l^{*}}}h_{i}$.

The expectation value can be approximated as

$\displaystyle\left|\text{Tr}[\rho_{i}A_{i}]-\text{Tr}[\rho_{i}^{(2l_{*})}A^{l_{*}}_{i}]\right|$
(141)
$\displaystyle\leq\left|\text{Tr}[\rho_{i}A_{i}]-\text{Tr}[\rho_{i}A^{l_{*}}_{i}]\right|+\left|\text{Tr}[\rho_{i}A^{l_{*}}_{i}]-\text{Tr}[\rho_{i}^{(2l_{*})}A^{l_{*}}_{i}]\right|$
$\displaystyle\leq||A_{i}-A_{i}^{l^{*}}||+||A_{i}^{l^{*}}||\times||\text{Tr}_{\setminus
V_{l^{*}}}[\rho_{i}]-\text{Tr}_{V_{2l^{*}}\setminus
V_{l^{*}}}[\rho_{i}^{(2l_{*})}]||_{1}.$

This follows from the triangle inequality. The partial trace $\setminus
V_{l^{*}}$ is over the support of $\rho_{i}$ excluding vertices $V_{l^{*}}$.
Eq. (141) now has a form that we can upper bound. Since
$||A_{i}^{l^{*}}||\leq||A_{i}-A_{i}^{l^{*}}||+||A_{i}||\leq
e^{\mathcal{O}(\beta h)}$, we can use Eq. (140) to bound the first term, and
Eq. (89) with $|B|=l^{*}$ to bound the second. With these, we conclude that
there exists constants $c_{1},c_{2}$ depending on all the constants involved
(i.e. $\beta,h,J,k,c^{\prime},v$) such that

$\left|\text{Tr}[\rho_{i}A_{i}]-\text{Tr}[\rho_{i}^{(2l_{*})}A^{l_{*}}_{i}]\right|\leq
c_{1}e^{-c_{2}l^{*}}.$ (142)

The key feature of $\text{Tr}[\rho_{i}^{(2l_{*})}A^{l_{*}}_{i}]$ is that it is
an expectation value of an operator whose form we known explicitly, as per Eq.
(50), evaluated in a thermal state of size $2l^{*}$. This can be computed
exactly (or rather, with a subleading error) in a time
$\text{exp}\left(\mathcal{O}(l^{*})\right)$. Let us now choose a precision
$\varepsilon/N$ in Eq. (142), so that
$l^{*}=\mathcal{O}\left(\log{N/\varepsilon^{-1}}\right)$. This way, we have

$\displaystyle Z^{\prime}\equiv
d^{N}\prod_{i=1}^{|E|}\text{Tr}[\rho_{i}^{(2l_{*})}A^{l_{*}}_{i}]$ (143)
$\displaystyle=Z\prod_{i=1}^{|E|}\left(1+\frac{\varepsilon}{N\text{Tr}[\rho_{i}A_{i}]}\right)=Z\left(1+\mathcal{O}\left(\varepsilon\right)\right).$
(144)

The last equation comes from the fact that $N\propto|E|$ and that all
eigenvalues of $A_{i}$ are $\mathcal{O}(1)$, as per the definition in Eq.
(48). The algorithm thus consists of exactly calculating the numbers
$\\{\text{Tr}[\rho_{i}^{(2l_{*})}A^{l_{*}}_{i}]\\}$ exactly, and then
multiplying them, so that

$\left|\log Z^{\prime}-\log Z\right|\leq\mathcal{O}\left(\varepsilon\right).$
(145)

Since there are $|E|\propto N$ terms in $Z^{\prime}$, and each takes time
$\text{poly}(N\times\varepsilon^{-1})$, the final runtime is
$\text{poly}(N,\varepsilon^{-1})$, as desired. See also Fawzi  _et al._ (2023)
for a related result in the translation invariant setting.

###  VII.2 Hardness of approximating partition functions

In the previous section we have seen how the partition function can be
approximated in 1D in the sense of Eq. (133) as long as the temperature is
$\beta=\mathcal{O}(1)$. Moreover, through the cluster expansion we briefly
explained in Sec. III.1 how it can be approximated for any local model as long
as $\beta<\beta^{*}$, where $\beta^{*}$ is some fixed constant independent of
system size.

On the other hand, in the limit of $\beta\rightarrow\infty$, the log-partition
function equals the energy of the ground state. For classical models,
approximating this to a certain precision is an NP-complete problem. For local
quantum Hamiltonians, it is QMA hard. This means that there should be no
efficient classical or quantum algorithm to approximate log-partition
functions for low enough temperatures, both for classical and quantum models.
In fact, it is known that the classical problem is only slightly harder than
NP Stockmeyer (1983) 666More specifically, it is in the class
$\text{BPP}^{\text{NP}}$, see Bravyi  _et al._ (2022)., and that it is at
least $\\#P$ hard if complex interactions are allowed Goldberg and Guo (2017).
For the quantum case, the exact complexity class to which this belongs or is
complete for is not yet clear (see Bravyi  _et al._ (2022) for more details
and results).

There is still the expectation that for certain classes of interesting models
we can still compute the partition function efficiently, even with classical
algorithms and at very low temperatures. One notable example are quantum Monte
Carlo methods Bravyi and Gosset (2017); Crosson and Harrow (2021); Crosson and
Slezak (2020), which are restricted to Hamiltonians without the so-called
“sign problem” (or _stoquastic_ Bravyi  _et al._ (2008)). Other results cover
different specific kinds of models Bravyi  _et al._ (2022); Helmuth and Mann
(2022). Quantum algorithms for approximating general partition functions also
exists Poulin and Wocjan (2009); Bravyi  _et al._ (2022); Cade and Montanaro
(2018); Chowdhury  _et al._ (2021), but often come with exponential run-times.

Another relevant angle of this problem is the connection of efficient
algorithms to the idea of completely analytical interactions from Sec. I.2, as
well as the physics of phase transitions. The intuition is that a physical
phase transition in the system may come together with a _computational_ phase
transition in which approximating $\log Z$ becomes fundamentally harder. Along
these lines it has been shown that in quite a general setting Harrow  _et al._
(2020), that the analiticity of the log-partition function implies the
existence of an efficient algorithm. This can be understood in terms of the
setting of Sec. III.1: as long as the Taylor expansion converges well, we can
compute the individual coefficients of Eq. (17) efficiently.

###  VII.3 Tensor network methods

Tensor network (TN) techniques are perhaps the most successful way of
classically computing physical properties of quantum systems in 1D, and
sometimes 2D and often come with rigorous theoretical guarantees. This
includes most notably the regime of low energy physics Verstraete and Cirac
(2006); Hastings (2007b); Landau  _et al._ (2015); Arad  _et al._ (2017);
Huang (2015) and, as we now review, that of finite temperature too. See e.g.
Bridgeman and Chubb (2017); Eisert (2013); Orús (2019); Cirac  _et al._ (2021)
for introductory texts to this topics.

The main aim is to obtain a TN representation of an operator $M_{D}$ such that
$||e^{-\beta H}-M_{D}||_{1}\leq\varepsilon Z$, which then allows us to compute
all thermal expectation values up to error $\varepsilon$ as per Eq. (6). The
index $D$ labels the bond dimension which, roughly speaking, quantifies the
complexity of representing $M_{D}$. A TN of bond dimension $D$ requires a
memory $\propto N\times D^{2}$ to be stored. Intuitively, the approximation
operator $M_{D}$ should be made out of a sum or low-depth product of operators
with smaller support i.e. of size at most $\propto\log D$. This is graphically
described for 1D in Fig. 8.

.

Figure 8: Schematically, the way to prove that a 1D thermal state is a tensor
network is by decomposing it as a product of smaller operators. It then
follows from standard methods that the bond dimension of the tensor network
representation is related to the size of those operators.

One possible way to do this is by adding piece by piece from right to left,
aided by the results from Sec. III.2. To do this, first define segments of the
chain of length $l$ such that the Hamiltonian up to segment $j-1$ is defined
by $H_{j}$. This allows us to define the operator

$\Psi_{j}=e^{\beta H_{j}}e^{-\beta H_{j+1}},$ (146)

such that $e^{-\beta H}=e^{-\beta h_{1}}\prod_{j=1}\Psi_{j}$. Each operator
$\Psi_{j}$ can then be approximated by a localized operator $\Psi_{j}^{l}$
with support in a region of length $l$ as in Eq. (40), exponentially well in
$l$. This can take the form

$\Psi^{l}_{j}=e^{\beta h_{j}}e^{-\beta h_{j+1}},$ (147)

There are $\mathcal{O}(N)$ of those operators, and the error $e^{-\Omega(l)}$
of each approximation can be shown to contribute additively. Thus, choosing
$l\propto\log{N/\varepsilon}$ gives the desired $\epsilon$-good approximation
to $e^{-\beta H}$. By construction, the product of operators resembles that of
Fig. 8. In that case, the bond dimension can be straightforwardly assumed to
be $D\leq e^{\mathcal{O}(l)}=\text{poly}(N,\varepsilon^{-1})$, which is
already computationally efficient and likely close to optimal.

A non-trivial improvement to this can be found in Kuwahara  _et al._ (2021).
Roughly speaking, one can define an operator $\tilde{\Psi}_{j}^{l}$ in which
the exponential functions of Eq. (147) are approximated by their Taylor
series. In that case, we can put forward results about the bond dimension
required to represent polynomials of Hamiltonians Arad  _et al._ (2013). This
leads to an improvement of the bond dimension to $D\leq
e^{\tilde{\mathcal{O}}(\sqrt{l})}=\text{exp}\left(\tilde{\mathcal{O}}(\sqrt{\log(N/\varepsilon)})\right)$.
This is sub-linear in system size, much more computationally efficient.

In higher dimensions, the best-known method is a variation of the cluster
expansion proposed in Hastings (2006). There, the expansion is treated in a
slightly different way, to approximate the exponential $e^{-\beta H}$ rather
than the log-partition function as in Sec. III.1. Instead of counting the
number of individual clusters of at most size $m$, one has to consider
arbitrary products $\prod_{i\in W}h_{i}\equiv h(W)$ of terms $h_{i}$ from a
multiset $W$ of size $|W|$, that can be divided into connected clusters. Let
us label the multisets $W=\\{h_{i}\\}$ in which the biggest cluster has size
$M$ to be $C_{M}$. That this is consistent with the cluster expansion can be
seen from taking the exponential of Eq. (17) given the expression in terms of
clusters of the powers in Eq. (21).

The following result was proven in detail in Kliesch  _et al._ (2014). It
reads

$||e^{-\beta H}-\sum_{W\in C_{M}}\frac{(-\beta)^{|W|}}{|w|!}h(W)||_{1}\leq
Z\left(e^{N\frac{b(\beta)^{M}}{1-b(\beta)}}-1\right),$ (148)

where $b(\beta)<1$ for all $\beta<\beta^{*}=\mathcal{O}(1)$. In Molnar  _et
al._ (2015) the sum over clusters on the RHS of Eq. (148) was shown to be a
tensor network (in fact, a so-called PEPO) of bond dimension
$e^{\mathcal{O}(M)}$. Thus, by setting the RHS to be $\varepsilon$, we achieve
a TN approximation to $e^{-\beta H}$ with bond dimension
$\text{poly}(N,\varepsilon^{-1})$. This only holds for inverse temperatures
below $\beta^{*}$, but the result can be extended to arbitrary temperatures
simply by taking powers of the operator. This means the bond dimension grows
as $D\leq\exp\left(\mathcal{O}\left(\beta\log{\frac{\beta
N}{\varepsilon}}\right)\right)$ (see Molnar  _et al._ (2015) for the details).
This scheme has recently been numerically implemented in practice Vanhecke
_et al._ (2021).

These results might seem surprising, since they show that there are in
principle efficient TN representations for all dimensions and all temperatures
$\beta=\mathcal{O}(1)$. This contradicts the intuition (justified by numerical
works Czarnik  _et al._ (2012, 2019); Czarnik and Corboz (2019); Czarnik  _et
al._ (2021); Wietek  _et al._ (2021); Ran  _et al._ (2013)) that, at phase
transitions, when long-range correlations are present, such efficient schemes
should not exist. The caveat, however, is that in dimensions higher than one,
a TN representation is not enough to be able to extract numerical data
efficiently. This is because the contraction of TN can be a computationally
demanding task by itself Schuch  _et al._ (2007); Haferkamp  _et al._ (2020).
In fact, what we expect is that the ability to reliably contract a higher
dimensional tensor network is related to facts such as local
indistinguishability Verstraete and Cirac ; Cirac  _et al._ (2013), which
allows us to obtain reliable results by contracting suitably smaller regions.

Finally, let us note that by using the local indistinguishability from Sec.
V.1 (or even without it in 1D Huang (2021)) it can be shown that a much
smaller bond dimension is needed to simulate local properties Alhambra and
Cirac (2021) .

###  VII.4 Quantum algorithms for preparing thermal states

One of the most promising applications of quantum computers is the generation
of exotic states of matter in complex many-body models. The expectation is
that this should allows us to discover a potentially wide variety of physics,
and also serve as a subroutine in certain quantum algorithms, such as those
performing optimization tasks.

Because of this, a question that has been very much explored lately is that of
how to prepare thermal states of local Hamiltonians with a quantum computer.
This could be either a fully fledged fault-tolerant one or one more suitable
for the so-called NISQ (Noisy Intermediate Scale Quantum) devices. In the
following we review some of the currently existing ones, and also explain the
ideas that highlighting the complexity of the problem. We mostly focus on
those that have some provable performance guarantees. There are many others we
will not cover (such as e.g. Motta  _et al._ (2020); Lu  _et al._ (2021); Cohn
_et al._ (2020) and others), which often rely on some level of heuristic
arguments. These also include approaches such as variational algorithms Farhi
_et al._ (2014); Chowdhury  _et al._ (2020); Wang  _et al._ (2021); Cerezo
_et al._ (2021) or those based on quantum versions of metropolis sampling
Terhal and DiVincenzo (2000); Temme  _et al._ (2011); Yung and Aspuru-Guzik
(2012). These may nonetheless be more efficient in many physically relevant
settings.

#### General considerations

We have very strong evidence pointing that preparing thermal states is, in its
most general setting, not an easy task. The results on QMA hardness of the
local Hamiltonian problem Kempe  _et al._ (2005) show that there are
vanishingly small temperatures (scaling quickly with system size) at which the
preparation of $\rho_{\beta}$ is QMA complete. That is, not even a quantum
computer can do it efficiently Kitaev  _et al._ (2002). This is the case even
for 1D systems Aharonov  _et al._ (2009).

There are also reasons to believe that the problem is not easy even at
slightly higher temperatures. For instance, it was recently shown Anshu  _et
al._ (2022b) (following a famous conjecture Hastings (2013)) that there are
local models for which preparing states below a certain energy density
(including low temperature Gibbs states) requires a circuit of depth at least
$\log N$. This property, however, does not apply to lattices Aharonov  _et
al._ (2013).

In that sense, we expect that large interesting classes of models and
temperature ranges will have efficient algorithms. The locality of the model,
and some of its consequences from the previous sections, should often simplify
this task.

#### Algorithms based on purifications

These algorithms work for general Hamiltonians, and are often designed to be
run in a fully fault tolerant quantum computer, capable of applying any
quantum circuit without large errors. They aim to construct the following
state

$|\rho_{\beta}\rangle=\frac{1}{\sqrt{Z}}\sum_{l}e^{-\beta
E_{l}/2}|E_{l}\rangle_{A}|l\rangle_{\mathcal{A}},$ (149)

where the second subsystem $\mathcal{A}$ is made of auxiliary particles with
an orthogonal basis $\\{|l\rangle\\}$ such that, upon tracing out, yield
$\text{Tr}_{\mathcal{A}}[|\rho_{\beta}\rangle\langle\rho_{\beta}|]=\rho_{\beta}$.
The first stage of the algorithm involves preparing a state $|\psi\rangle$
with $|\rho_{\beta}\rangle$ as a component such that

$|\psi\rangle=\frac{1}{\mathcal{N}}|\rho_{\beta}\rangle|0\rangle_{R}+....\quad,$
(150)

where we have included a possible additional register $R$.

The first works proposing this scheme Poulin and Wocjan (2009); Chiang and
Wocjan (2010); Bilgin and Boixo (2010) instead apply the phase estimation
algorithm Kitaev (1995). To a good approximation, this algorithm acts as
follows 777This output is stricly speaking only approximate, but the error of
this approximation can be dealt with using additional registers (see Poulin
and Wocjan (2009))

$\displaystyle
U_{PE}(|E_{l}\rangle_{A}|0\rangle_{\mathcal{A}})=|E_{l}\rangle_{A}|l\rangle_{\mathcal{A}},$
(151)

that is, it “measures” the energy of system $A$ into the register
$\mathcal{A}$. Inputting a uniform superposition
$\frac{1}{d^{N/2}}\sum_{l}|E_{l}\rangle$ yields

$\displaystyle
U_{PE}(\frac{1}{d^{N/2}}\sum_{l}|E_{l}\rangle|0\rangle_{\mathcal{A}})=\frac{1}{d^{N/2}}\sum_{l}|E_{l}\rangle|l\rangle_{\mathcal{A}}.$
(152)

Now we add an additional qubit register on the state $|0\rangle_{R}$ and
rotate it to $|\theta\rangle=\cos\theta|0\rangle+\sin\theta|0\rangle$ by an
angle $\theta(E_{l},\beta)=\arccos(e^{-\frac{\beta E_{l}}{2}})$ conditioned on
the system, obtaining

$\displaystyle\frac{1}{d^{N/2}}\sum_{l}|E_{l}\rangle|l\rangle_{\mathcal{A}}|\theta(E_{l},\beta)\rangle_{R}=\sqrt{\frac{Z}{d^{N}}}|\rho_{\beta}\rangle|0\rangle_{R}+....\quad$
(153)

which is the target state with $\mathcal{N}=\sqrt{\frac{d^{N}}{Z}}$. Other
approaches use more recent quantum simulation ideas, such as the technique
based on sums of unitaries Chowdhury and Somma (2017), which leads to a better
dependence on the approximation error $\epsilon$ in many cases of interest.

Finally, to obtain $|\rho_{\beta}\rangle$ with high precision, one must then
apply amplitude amplification of the state $|\psi\rangle$, to output the
component of Eq. (153) corresponding to the register state $|0\rangle_{R}$.
The gate complexity of this, however, grows linearly in $\mathcal{N}$, which
sets the leading (almost) exponential gate cost of the algorithm.

There already exists improvements to this type of scheme. In one dimension,
one can instead implement this same algorithm connecting subsequent segments
of the chain, which can reduce the gate complexity to a polynomial $\sim
N^{\mathcal{O}(\beta)}$ Bilgin and Boixo (2010). Also, recent progress shows
that the phase estimation and amplitude amplification steps in these schemes
can instead be replaced by random circuits with post-selection Shtanko and
Movassagh (2021), making them more amenable to current technologies. For
commuting Hamiltonians, a purification in the form of a tensor network state
(a PEPS) can be very efficiently prepared through an adiabatic algorithm Ge
_et al._ (2016). See also the recent Holmes  _et al._ (2022), which produces a
purification of a thermal state $\propto e^{-\beta H_{1}}$ starting from that
of another Hamiltonian $H_{0}$, and is efficient when $||H_{0}-H_{1}||$ is not
too large. A potentially efficient scheme along these lines is the one
presented in Motta  _et al._ (2020).

#### Efficient algorithms from physical features

Perhaps the main caveat of most of the aforementioned algorithms is that they
are constructed for very general Hamiltonians: they do not always make a very
clear use of the physical features that we expect could simplify the problem,
such as locality or any one of its consequences. However, we expect that there
exists efficient schemes to prepare Gibbs states of local model that belong to
the class of “completely analytical interactions”, aided by facts such as
decay of correlations.

This is the case for the proposal in Brandão and Kastoryano (2019), whose
efficiency depends on two such factors: the speed of decay of CMI from Sec.
IV.3, and the error in the local indistinguishability from Sec. V.1. The
algorithm uses iterations of the recovery map that appeared in Eq. (70), which
are guaranteed to yield a low error if the CMI decays quickly enough. The main