# Certified algorithms for equilibrium states
of local quantum Hamiltonians

Hamza Fawzi Omar Fawzi Samuel O. Scalet11footnotemark: 1 Department of Applied
Mathematics and Theoretical Physics, University of Cambridge, United
Kingdom.Univ Lyon, Inria, ENS Lyon, UCBL, LIP, Lyon, France

###### Abstract

We design algorithms for computing expectation values of observables in the
equilibrium states of local quantum Hamiltonians, both at zero and positive
temperature. The algorithms are based on hierarchies of convex relaxations
over the positive semidefinite cone and the matrix relative entropy cone, and
give certified and converging upper and lower bounds on the desired
expectation value. In the thermodynamic limit of infinite lattices, this shows
that expectation values of local observables can be approximated in finite
time, which contrasts with recent undecidability results about properties of
infinite quantum lattice systems. In addition, when the Hamiltonian is
commuting on a 2-dimensional lattice, we prove fast convergence of the
hierarchy at high temperature leading to a runtime guarantee for the algorithm
that is polynomial in the desired error.

##  1 Introduction

Condensed matter physics studies the physical properties of equilibrium states
of many-body systems. The central question is to determine properties of the
many-body state as a function of the interaction between the constituents of
the system. The topic of Hamiltonian complexity [Osb12, GHL+15] studies this
question from a complexity-theoretic point of view. The results in Hamiltonian
complexity suggest that efficient algorithms answering this question are
unlikely to exist. In fact, determining the ground energy is QMA-hard [KKR06],
even for translation-invariant systems [GI09]. Computing expectation values of
observables in the ground state is even harder than the ground energy [Amb14].
In the thermodynamic limit when the number of systems is taken to infinity,
finding good approximations to the energy can be hard for fixed Hamiltonians
[AI22, WC22] and the spectral gap is even uncomputable [CPGW15].

#### Certified algorithms

These complexity results put severe limitations on provably efficient and
correct classical and quantum algorithms for the fundamental questions in
condensed matter physics. However, these limitations only apply to algorithms
that have the required properties _for all_ valid instances of the problem.
Moreover, typically the instances showing hardness are highly contrived. In
order to avoid these limitations, we consider here _certified algorithms_
where we require the correctness condition for all instances, but relax the
“provable efficiency for all instances” condition. We say that an algorithm
for computing $\mathsf{p}^{*}(H)$ is _certified_ if on input $H$, it outputs a
pair of numbers $(\mathsf{p}^{\min}_{\ell},\mathsf{p}^{\max}_{\ell})$ such
that we have
$\mathsf{p}^{*}(H)\in[\mathsf{p}^{\min}_{\ell},\mathsf{p}^{\max}_{\ell}]$ for
all $H$ and all $\ell$. We further require that as $\ell\to\infty$, the
interval $[\mathsf{p}^{\min}_{\ell},\mathsf{p}^{\max}_{\ell}]$ converges to
the desired value $\mathsf{p}^{*}(H)$. Here, $\ell$ is a parameter that
governs the runtime of the algorithm such that fast convergence in $\ell$
corresponds to an efficient algorithm. The above cited complexity results show
e.g., that if $\mathsf{p}^{*}(H)$ is the ground energy of $H$, then for small
$\ell$, the interval $[\mathsf{p}^{\min}_{\ell},\mathsf{p}^{\max}_{\ell}]$ has
to be large for some $H$. However, for any input $H$ of interest, we can
always run the algorithm and the returned information will be correct and
certifiably so. As such, no a priori analysis of the convergence speed is
needed in order to obtain rigorous approximations of the quantity of interest.

#### Setup and main result

We now present the setup in an informal way. We refer to Section 2 for more
formal statements. We assume that our system is described by a local
Hamiltonian on the lattice $\mathbb{Z}^{D}$, which can be formally written as

$H=\sum_{X\subset\mathbb{Z}^{D}}h_{X}$

where $h_{X}$ is the interaction term acting only on the finite set of sites
$X\subset\mathbb{Z}^{D}$. We assume the Hamiltonian $H$ to be translation-
invariant, and finite range, meaning that $h_{X}$ is nonzero only when $X$ is
contained in a ball of fixed radius for a standard metric.111The relaxations
we define can be easily adapted for arbitrary local Hamiltonians on a finite
number of sites. The set of equilibrium states at temperature $T\geq 0$ is the
set of states $\rho$ that minimize $\operatorname{tr}(H\rho)-TS(\rho)$, where
$\operatorname{tr}(H\rho)$ is the energy, and $S(\rho)$ is the entropy. Note
that for infinite systems the thermal equilibrium state is not necessarily
unique, and this property lies at the heart of the existence of phase
transitions.

Given the description of such a system, our objective is to determine physical
properties of the corresponding equilibrium states. For example, consider a
local observable $O_{x}$ acting on a single lattice site $x\in\mathbb{Z}^{D}$,
e.g., the magnetization in the z-direction $O_{x}=Z_{x}$ for a spin-$1/2$
system. A physical property of interest is then the magnetization per site
which is given by the average expectation value of the observables $O_{x}$. If
the equilibrium state $\rho$ is translation-invariant then for any
$x,y\in\mathbb{Z}^{D}$, $\operatorname{tr}(\rho O_{x})=\operatorname{tr}(\rho
O_{y})$ is equal to the average value of $O$ on the system. Consider now the
set of values $\operatorname{tr}(\rho O_{x})$ that an equilibrium state at
temperature $T=1/\beta$ can take. Due to nonuniqueness of equilibrium states,
this set is in general an interval which we will denote by

$[\langle O\rangle^{\min}_{\beta},\langle
O\rangle^{\max}_{\beta}]:=\left\\{\operatorname{tr}(\rho O_{x}):\rho\text{ is
a translation-invariant equilibrium state at temperature
$T=1/\beta$}\right\\}\subset\mathbb{R}.$

In fact, a criterion to detect nonuniqueness of equilibrium states is $\langle
O\rangle^{\min}_{\beta}<\langle O\rangle^{\max}_{\beta}$. A prominent
(classical) example of this is the spontaneous magnetization of the 2D Ising
model below the critical temperature, whereby $O$ is the magnetization and
$\langle O\rangle^{\min}_{\beta}<0<\langle O\rangle^{\max}_{\beta}$ for all
$\beta$ above some critical value. We can now state our main result.

######  Theorem 1.1 (Certified algorithms for translation-invariant
expectation values, informal).

There is a family of finite and explicit convex optimization problems
$(\mathsf{P}^{\min}_{\ell})_{\ell\in\mathbb{N}}$ and
$(\mathsf{P}^{\max}_{\ell})_{\ell\in\mathbb{N}}$ (depending on $H$ and
$\beta\in[0,+\infty]$ and $O$) satisfying the following properties:

  * •

$(\mathsf{P}^{\min}_{\ell})$ and $(\mathsf{P}^{\max}_{\ell})$ have at most
$\exp(\mathcal{O}(\ell^{D}))$ variables and $\exp(\mathcal{O}(\ell^{D}))$
constraints

  * •

The optimal values $\mathsf{p}^{\min}_{\ell}$ and $\mathsf{p}^{\max}_{\ell}$
of $(\mathsf{P}^{\min}_{\ell})$ and $(\mathsf{P}^{\max}_{\ell})$ satisfy

$\mathsf{p}^{\min}_{\ell}\leq\langle O\rangle^{\min}_{\beta}\leq\langle
O\rangle^{\max}_{\beta}\leq\mathsf{p}^{\max}_{\ell}$

and $\mathsf{p}^{\min}_{\ell}\uparrow\langle O\rangle^{\min}_{\beta}$ and
$\mathsf{p}^{\max}_{\ell}\downarrow\langle O\rangle^{\max}_{\beta}$ as
$\ell\to\infty$.

Numerical results illustrating this algorithm can be found in Section 4. We
note that it is straightforward to adapt the hierarchies we develop to finite
systems and convergence will be exact when $\ell$ is taken to be the size of
the system.

#### Thermodynamic limit and decidability

We chose to focus on the thermodynamic limit due to its importance for
understanding phase transitions. In fact, phase transitions can be
characterized by nonanalyticities in the interval $[\langle
O\rangle^{\min}_{\beta},\langle O\rangle^{\max}_{\beta}]$ as a function of
$\beta$ or the parameters of the Hamiltonian. One difficulty we would like to
highlight regarding the thermodynamic limit concerns the definition of
equilibrium states at any given temperature $T\geq 0$. Indeed, it was actually
shown in [BCW21] that computing local observables on the ground state of
infinite systems is undecidable. This might seem to contradict our result.
This is not the case however, due to a subtle point in the definition of
ground-state observables in the infinite limit. The authors of [BCW21] define
those as limits of ground-state observables in finite systems, i.e., the
observable is computed for the ground-state of a truncated Hamiltonian with
open boundary condition $H_{\Lambda}=\sum_{X\subset\Lambda}h_{X}$, where
$\Lambda\subset\mathbb{Z}^{D}$ is finite, and afterwards the limit
$\Lambda\uparrow\mathbb{Z}^{D}$ is taken. The definition we use in our work,
which is standard in the operator algebraic framework, is more general and in
particular contains limits of ground-state observables of finite system
Hamiltonians with any choice of boundary conditions. It is simple to construct
a Hamiltonian where fixing the boundary condition while taking the limit
excludes natural ground states (see Remark 2.1 and also Remark 3.3 concerning
thermal states of the 2D Ising model). The comparison however illustrates that
even the existence of any convergent algorithm for observables in ground
states is far from obvious, which is what we achieve in Theorem 1.1. For the
operator algebraic definition of observables we define the following promise
problem.

######  Definition 1.2.

Equilibrium Observable Problem
---
Input: |  Local dimension $d$, local Hamiltonian term $h$ with rational coefficients, local observable $O$ with rational coefficients, rational numbers $a$ and $b$ with $b>a$, temperature $\beta\in[0,+\infty]$
Promise: |  Either $\langle O\rangle_{\beta}^{\min}>b$ or $\langle O\rangle_{\beta}^{\max}<a$
Question: |  Output YES if $\langle O\rangle_{\beta}^{\min}>b$, NO otherwise

Note, that while our definition of thermal states is more general, an
equivalent promise for the undecidable problem is guaranteed in [BCW21]. We
prove as a consequence of Theorem 1.1:

######  Theorem 1.3 (Decidability of translation-invariant ground state and
thermal observables).

The Equilibrium Observable Problem is decidable.

#### Main ingredients

As we are dealing with infinite systems, we adopt the operator algebraic
formulation of quantum mechanics [BR87, BR97]. In this setting, a state is
defined as a normalized positive linear functional
$\omega:\mathcal{A}\to\mathbb{C}$ on the algebra of observables
$\mathcal{A}$.222We omit the formal definition of $\mathcal{A}$ at this point,
and refer to Section 2 instead. At $T=0$, an equilibrium state or ground state
is defined via the condition

$\omega(a^{*}[H,a])\geq 0\;\;\forall a\in\mathcal{A}$ (1)

where $[H,a]$ is the commutator. Intuitively, the condition above expresses
the fact that the energy of $\omega$ has to increase under any local
perturbation. At positive temperature $T=1/\beta>0$, thermal equilibrium
states are defined via the Kubo-Martin-Schwinger (KMS) condition333Equation
(2) is the KMS condition for finite systems. See Section 2 for the rigorous
statement of the KMS condition for infinite systems.

$\omega(ba)=\omega(ae^{-\beta H}be^{\beta H})\quad\forall a,b\in\mathcal{A}.$
(2)

In the classical case, when the Hamiltonian is diagonal in a basis
$\\{|\sigma\rangle\\}_{\sigma}$, Equation (2) reduces to

$\omega(|\sigma\rangle\\!\langle\sigma|)=\omega(|\sigma^{\prime}\rangle\\!\langle\sigma^{\prime}|)\exp(-\beta(\langle\sigma|H|\sigma\rangle-\langle\sigma^{\prime}|H|\sigma^{\prime}\rangle)),$
(3)

for all $\sigma,\sigma^{\prime}$. When $\sigma^{\prime}$ is obtained from
$\sigma$ by flipping a single spin, these are known as the spin-flip equations
and are exploited in the bootstrap approach for the classical Ising model
[CGL+22].

To obtain Theorem 1.1, a natural idea is to relax the set of $\beta$-KMS
states and impose the condition (2) for a subset of observables $a,b$ lying in
some finite-dimensional space. However, the main obstruction one is faced with
is that even if $a,b$ are local observables the expression $ae^{-\beta
H}be^{\beta H}$ is in general not local, except for commuting Hamiltonians. As
such, even though (2) form a set of linear equations on $\omega$, they involve
the value of $\omega$ on nonlocal observables. We circumvent this issue by
using another characterization of thermal equilibrium states via so-called
Energy-Entropy Balance (EEB) inequalities [AS77]: this is an infinite set of
scalar convex inequalities, each indexed by a local observable
$a\in\mathcal{A}$, which carve out the set of $\beta$-KMS states. On the one
hand these inequalities are better suited than the $\beta$-KMS condition
because they preserve locality, i.e., they only require the value of the state
$\omega$ on a finite region around the support of $a$. On the other hand, a
drawback of these inequalities is that there are infinitely many of them, and
unlike the $\beta$-KMS condition their dependence on $a$ is nonlinear, and so
it is not sufficient to choose a basis of local observables. A key ingredient
to prove our theorem is to formulate a matrix generalization of such
inequalities using the matrix-valued relative entropy function

$D_{op}(A\|B)=A^{1/2}\log(A^{1/2}B^{-1}A^{1/2})A^{1/2}.$

We show that an infinite set of scalar Energy-Entropy Balance inequalities can
be compactly formulated by a single convex matrix inequality of the form
$D_{op}(A_{\omega}\|B_{\omega})\preceq\beta C_{\omega}$ for some suitable
matrices $A_{\omega},B_{\omega},C_{\omega}$ that depend linearly on $\omega$.
To the best of our knowledge, this is the first application of the matrix
relative entropy function in the design of convex relaxations in quantum
information. Optimization problems involving this function can be solved using
interior-point methods [FS23] or via semidefinite programming approximations
[FSP19].

#### Related work

A special case of the problem considered in this work is when the observable
$O$ is the energy. Certified algorithms do exist for the ground energy of
local Hamiltonians: one can combine semidefinite programming relaxations
[Maz06, PNA10, BP12, BH12] for lower bounds and variational methods such as
tensor networks [BC17, CPGSV21, Whi92] for upper bounds. Using such two-sided
bounds on the ground energy, the recent works [WSF+23, Han20] obtain bounds
for expectation values of local observables in the ground state, although no
convergence guarantees are given. The approach based on imposing the
additional constraint (1), which leads to convergence guarantees, was proposed
independently and concurrently in [AKV+23], where it was derived as a special
case of a method to strengthen semidefinite relaxations for noncommutative
polynomial optimization problems by incorporating optimality conditions as
additional constraints. These papers however do not discuss the case of
positive temperature since the equilibrium states are not characterized by a
noncommutative polynomial optimization problem which is linear in the state.
One could try to use lower bounds for the free energy coming from convex
relaxations [PH11], but it is not clear how to use such bounds for
observables. We note however that for classical systems, thermal observables
can be obtained via the bootstrap approach [CGL+22, PRV19], by directly
imposing the KMS conditions (3). The resulting hierarchies were shown to be
asymptotically convergent [CS23], even if only a weaker set of constraints are
imposed. For quantum systems, to the best of our knowledge, Theorem 1.1
provides the first certified algorithms for general observables in thermal
states.

Note that, in some restricted settings, such as 1D gapped systems [LVV15] at
zero temperature or 1D systems at positive temperature [Has06, MSVC15, KAA21,
FFS23b], there are provably efficient algorithms computing representations of
equilibrium states and thus expectation values, but such algorithms are
tailored to these settings. In addition, for arbitrary dimensions and high
temperatures provably efficient algorithms for computing the free energy exist
[KKB20, HMS20, MH21], which can be used to compute observables [BCGW22].

#### Quantifying convergence speed

Theorem 1.1 does not provide quantitative guarantees on the convergence speed.
As previously mentioned, one cannot hope to prove general fast convergence
guarantees as even the special case where the observable $O$ corresponds to
energy is unlikely to have a polynomial-time quantum algorithm, even when
$D=1$ [GI09].

However, one can expect provably fast convergence for some classes of
Hamiltonians. We illustrate this by showing exponential convergence in two
regimes for which it is known that no phase transitions can occur. In the
high-temperature regime and in one dimension at any nonzero temperature, the
set of equilibrium states reduces to a singleton, and for any local observable
$O$ we have $\langle O\rangle^{\min}_{\beta}=\langle
O\rangle^{\max}_{\beta}=\langle O\rangle_{\beta}$. Assuming a commuting local
Hamiltonian $H$, the theorem below shows exponential convergence to $\langle
O\rangle_{\beta}$ in the level $\ell$ of the convex optimization hierarchy.

######  Theorem 1.4 (Quantitative convergence rate, informal).

Assume that $D\leq 2$, and that $H$ is translation-invariant and _commuting_ ,
i.e., $[h_{X},h_{Y}]=0$ for all $X,Y\subset\mathbb{Z}^{D}$. For $D=1$ and
$\beta_{1}=\infty$ or for $D=2$ and some $\beta_{1}>0$, we have for all
$0\leq\beta<\beta_{1}$, and for any local observable $O$, $\langle
O\rangle^{\min}_{\beta}=\langle O\rangle^{\max}_{\beta}=\langle
O\rangle_{\beta}$ and, using the same notation as in Theorem 1.1

$\max\\{\mathsf{p}^{\max}_{\ell}-\langle O\rangle_{\beta},\langle
O\rangle_{\beta}-\mathsf{p}^{\min}_{\ell}\\}\leq c_{1}\|O\|e^{-c_{2}\ell}$

for some constants $c_{1},c_{2}>0$ depending on the dimension and the
interaction.

This result addresses an open problem raised in [CS23] about the speed of
semidefinite programming hierarchies for the (classical) Ising model, in
particular whether exponential convergence holds away from criticality.
Theorem 1.4 establishes such a statement for the more general class of
commuting local Hamiltonians.

#### Outlook

The preliminary numerical experiments in Section 4 demonstrate that the
approach proposed in this paper is very promising and could ultimately be used
in computing phase diagrams. Our result opens several avenues for future work:
First, one can exploit additional convex constraints that are known to hold
for marginals of equilibrium states. For example, one can add entropy-based
inequalities as proposed in [FFS23a]. Furthermore, other valid inequalities
tailored to specific models can be added to the convex relaxation such as
reflection positivity which was used in [CGL+22] for the Ising model. Second,
it would be very interesting to combine the methods from this paper with
variational methods (e.g., from tensor networks) to obtain more accurate
bounds, such as the recent work on the ground energy problem [KSDN22]. Another
natural question is whether one can use our algorithms within hybrid
classical-quantum algorithms for quantum simulation problems [ZWJ22]. On the
analytical side, it remains open whether convergence guarantees for other
classes of Hamiltonian can be proven. Promising candidates would be regimes in
which other classical algorithms are efficient such as thermal states in 1D
and at high-temperature (for general noncommuting Hamiltonians) as well as
gapped ground states in 1D.

##  2 Preliminaries

We consider spin systems on the infinite $D$-dimensional lattice
$\mathbb{Z}^{D}$ where $D\geq 1$ and we adopt the operator algebraic point of
view. We describe the setup briefly here, and we refer to [Nac06] or [BR97,
Section 6.2] for more details.

For any site $x\in\mathbb{Z}^{D}$, the local Hilbert space
$\mathcal{H}_{x}\simeq\mathbb{C}^{d}$ has dimension $d$. The Hilbert space
associated to a finite subset $X\subset\Lambda$ is the tensor product
$\mathcal{H}_{X}=\otimes_{x\in X}\mathcal{H}_{x}$, and we let
$\mathcal{A}_{X}=\mathcal{B}(\mathcal{H}_{X})$ be the algebra of observables
supported on $X$. If $x\in\mathbb{Z}^{D}$, we let
$\tau_{x}:\mathcal{A}_{X}\to\mathcal{A}_{X+x}$ be the translation operator by
$x$. Let

$\mathcal{A}_{\text{loc}}=\cup_{X\subset\mathbb{Z}^{D}}\mathcal{A}_{X}$ (4)

be the set of local observables, and
$\mathcal{A}=\overline{\mathcal{A}_{\text{loc}}}$ be its completion with
respect to the operator norm, i.e., $\mathcal{A}$ is the $C^{*}$-algebra of
quasi-local observables. A state is a linear functional
$\omega:\mathcal{A}\to\mathbb{C}$ such that $\omega(1)=1$,
$\omega(a^{*})=\overline{\omega(a)}$, and $\omega(a^{*}a)\geq 0$ for all
$a\in\mathcal{A}$.

#### Hamiltonians

Consider a local translation-invariant Hamiltonian $H$, which can be formally
written as

$H=\sum_{\begin{subarray}{c}X\subset\mathbb{Z}^{D}\end{subarray}}h_{X}$

where $h_{X}\in\mathcal{A}_{X}$ are the local interaction terms, and the
summation is over finite subsets $X$ only. We require that the Hamiltonian is
translation-invariant, i.e., $h_{X+x}=\tau_{x}(h_{X})$ for all
$X\subset\mathbb{Z}^{D}$ and $x\in\mathbb{Z}^{D}$. We assume in this paper
that the Hamiltonian $H$ is finite-range, i.e., that $h_{X}$ is nonzero only
when $X$ is contained in a ball of certain fixed radius $r$ around any lattice
site. If $\Lambda$ is a finite subset of $\mathbb{Z}^{D}$, we let

$H_{\Lambda}=\sum_{X\subset\Lambda}h_{X}\in\mathcal{A}_{\Lambda}$ (5)

and

$\tilde{H}_{\Lambda}=\sum_{X\cap\Lambda\neq\emptyset}h_{X}\in\mathcal{A}_{\bar{\Lambda}}$
(6)

where $\bar{\Lambda}=\Lambda\cup\partial^{ex}\Lambda$, and
$\partial^{ex}\Lambda$ is the external boundary of $\Lambda$, i.e.,
$\partial^{ex}\Lambda=\\{y\in\Lambda^{c}:\exists
Y\subset\mathbb{Z}^{D},h_{Y}\neq 0,y\in Y,\Lambda\cap Y\neq\emptyset\\}$. Note
that even though $H$ is not a well-defined element of $\mathcal{A}$, for any
$a\in\mathcal{A}_{\text{loc}}$, we formally write $[H,a]$ for the well-defined
element of $\mathcal{A}_{\text{loc}}$ that is equal to

$[H,a]=[\tilde{H}_{\Lambda},a],$

when $a\in\mathcal{A}_{\Lambda}$.

#### Ground states

A state $\omega:\mathcal{A}\to\mathbb{C}$ is called a _ground state_ of $H$ if

$\omega(a^{*}[H,a])\geq 0\quad\forall a\in\mathcal{A}_{\text{loc}}.$ (7)

In the case of finite systems, this condition is equivalent to saying that the
density matrix of $\omega$ is supported on the eigenspace of $H$ of minimal
eigenvalue. To see this, we note that for finite systems, states are described
by density matrices $\omega(a)=\operatorname{tr}[\rho a]$. Choosing the
operators $a$ in the above definition for a given eigenbasis of the
Hamiltonian $e_{i}$ as $e_{j}e_{i}^{*}$, we see that

$\rho_{ii}E_{j}-\rho_{ii}E_{i}\geq 0$

where we denote matrix elements of $\rho_{ij}$ in the energy eigenbasis and
the corresponding eigenvalues by $E_{i}$. This is fulfilled only if
$\rho_{ii}=0$ for all $i$ corresponding to eigenvectors that are not in the
ground-space sector. By writing

$\rho=\left(\begin{matrix}A&B\\\ B^{*}&C\end{matrix}\right)$

where the blocks correspond to the ground-space and its orthogonal complement,
we can see that $C$ has nonpositive diagonal entries, but is also positive
semidefinite and thereby $0$. By the Schur complement lemma we can then
conclude that $B=0$ and thereby $\rho$ is supported on the ground space.

######  Remark 2.1 (Relation to thermodynamic limit).

It is immediate to verify that the limit of ground states of finite
Hamiltonians $H_{\Lambda}$ of increasing size will automatically satisfy
condition (7). Actually this is even true for any choice of “boundary
condition” $\Delta_{\Lambda^{c}}$ acting outside $\Lambda$. Indeed, one can
easily show that if $\omega$ is such that

$\omega(a)=\lim_{\Lambda\uparrow\mathbb{Z}^{D}}\langle\Phi^{\Lambda}|a|\Phi^{\Lambda}\rangle$

where for each $\Lambda$, $\Phi^{\Lambda}$ is a ground state of the finite
Hamiltonian $H_{\Lambda}+\Delta_{\Lambda^{c}}$, then $\omega$ will satisfy
(7). Imposing such arbitrary boundary conditions while taking the limit may
actually be needed in order to capture the different ground states in the
thermodynamic limit. To illustrate this, consider the following one-
dimensional classical system with three states $\\{0,1,2\\}$ per site given by
the Hamiltonian

$\displaystyle h_{i,i+1}(s,s^{\prime})=\left\\{\begin{array}[]{cc}0&\text{ if
}s=s^{\prime}\in\\{0,1\\}\\\ 2&\text{ if
}(s,s^{\prime})\in\\{(0,1),(1,0),(2,2),(2,1),(0,2),(1,2)\\}\\\ -1&\text{ if
}s=2,s^{\prime}=0.\\\ \end{array}\right.$

For the finite Hamiltonian $H_{\Lambda}$ with
$\Lambda=\\{-\ell,\dots,\ell\\}$, the unique ground state is the sequence
$(2,0,\ldots,0)$ with energy $-1$. Importantly, note that for the observable
$O=|0\rangle\\!\langle 0|_{0}$ at site $0$, we have that for any $\ell\geq 1$,
its value on the ground state of $H_{\Lambda}$ is $1$. However, another valid
ground state of the infinite chain is the all-1 state. For this ground state
in the thermodynamic limit, the observable takes a value of $0$.

#### Thermal equilibrium states

Thermal equilibrium states at an inverse temperature $\beta=1/T\geq 0$ are
most commonly defined in terms of the Kubo-Martin-Schwinger (KMS) conditions.
For finite systems, these conditions can be simply written as

$\omega(ba)=\omega(ae^{-\beta H}be^{\beta H})\quad\forall
a,b\in\mathcal{A}_{\text{loc}}.$ (8)

It is easy to verify that the Gibbs state
$\omega(a)=\operatorname{tr}(e^{-\beta H}a)/\operatorname{tr}e^{-\beta H}$
satisfies the equality conditions above, and that it is also the unique such
state. To state the KMS-conditions for infinite systems, one first needs to
introduce the time evolution operator

$\alpha_{t}(a)=\lim_{\Lambda\uparrow\mathbb{Z}^{D}}e^{iH_{\Lambda}t}ae^{-iH_{\Lambda}t}\qquad\forall
a\in\mathcal{A}_{\text{loc}},t\in\mathbb{R}$

where the limit is taken over a suitable sequence of finite subsets $\Lambda$.
For the finite-range interaction we consider here, $\alpha_{t}$ exists and is
strongly continuous [BR97, Theorem 6.2.4]. For a given time-evolution
operator, there is a dense subset $\mathcal{A}_{\alpha}\subset\mathcal{A}$
such that $F(\alpha_{t}(a))$ is entire for all $a\in\mathcal{A}_{\alpha}$ and
bounded linear functionals $F$ [BR87, Proposition 2.5.22]. We say that a state
$\omega:\mathcal{A}\to\mathbb{C}$ is a $\beta$-KMS state if [BR97, Definition
5.3.1]

$\omega(ba)=\omega(a\alpha_{i\beta}(b))\quad\forall
a,b\in\mathcal{A}_{\alpha}.$ (9)

In fact, it is sufficient to ensure this condition for any dense
$\alpha$-invariant *-subalgebra of $\mathcal{A}_{\alpha}$. (Notice that for
finite-dimensional systems, the condition above is the same as (8).) The set
of translation-invariant $\beta$-KMS states will be denoted by
$\mathcal{G}_{\beta}$:

$\mathcal{G}_{\beta}=\left\\{\omega:\mathcal{A}\to\mathbb{C}:\omega(1)=1,\;\omega(a^{*}a)\geq
0\;\forall
a\in\mathcal{A}_{\text{loc}},\;\omega(a)=\omega(\tau_{x}(a))\;\forall
a\in\mathcal{A}_{\text{loc}},x\in\mathbb{Z}^{D},\;\text{ and
}\eqref{eq:kms}\right\\}.$ (10)

The set $\mathcal{G}_{\beta}$ is nonempty for all $\beta\geq 0$. This can be
shown by taking an appropriate limit of finite Gibbs states on a sequence of
$\Lambda\uparrow\mathbb{Z}^{D}$. A crucial aspect of infinite systems (i.e.,
systems in the thermodynamic limit) is that the set of thermal equilibrium
states $\mathcal{G}_{\beta}$ is not necessarily reduced to a singleton. The
existence of many $\beta$-KMS states signals the presence of different
thermodynamic phases.

##  3 Convex relaxations via energy-entropy balance inequalities

Our goal is to obtain a convex relaxation for the marginals of equilibrium
states $\omega\in\mathcal{G}_{\beta}$. A natural way to obtain a convex
relaxation is to impose only a subset of the $\beta$-KMS equality conditions
(9). The problem with the KMS linear equality constraints is that even if
$a,b$ are local observables, the observable $a\alpha_{i\beta}(b)$ is in
general not local (unless, say, the Hamiltonian is commuting).

###  3.1 Energy-entropy balance inequalities

The energy-entropy inequalities is a family of convex inequalities that
characterize $\beta$-KMS states and that only involve local observables of
$\omega$. They were first derived by Araki and Sewell in [AS77].

######  Theorem 3.1 (Energy-entropy balance (EEB) inequalities, see [BR97,
Theorem 5.3.15]).

Let $H$ be a finite-range Hamiltonian. Then $\omega:\mathcal{A}\to\mathbb{C}$
is a $\beta$-KMS state if, and only if, it satisfies444The function
$x\log(x/y)$ is defined as equal to 0 if $x=0$ and $+\infty$ if $x>0$ and
$y=0$.

$\displaystyle\omega([H,a])=0\quad\forall a\in\mathcal{A}_{\text{loc}},$ (11a)
$\displaystyle\omega(a^{*}a)\log\left(\frac{\omega(a^{*}a)}{\omega(aa^{*})}\right)\leq\beta\omega(a^{*}[H,a])\quad\forall
a\in\mathcal{A}_{\text{loc}}.$ (11b)

######  Remark 3.2.

The cited theorem requires the conditions to hold for all elements in the
domain of the generator $\delta$ of the time-evolution operator $\alpha_{t}$.
On $\mathcal{A}_{\text{loc}}$ this operator takes values $\delta(a)=i[H,a]$.
It is, however, defined on a larger domain given by the closure of this
operator, i.e., the commutator formula on $\mathcal{A}_{\text{loc}}$ is a core
for $\delta$. By the definition of the closure of an operator,
$a\in\mathcal{A}$ is said to be in the domain of $\delta$, i.e., $a\in
D(\delta)$ if there exists a sequence $a_{i}\in\mathcal{A}_{\text{loc}}$ such
that $a_{i}\to a$ and $\delta(a_{i})\to\delta(a)$. Due to the latter, it is
sufficient to impose the EEB inequalities on $\mathcal{A}_{\text{loc}}$ only,
despite the lack of continuity of $\delta$.

The convexity of the constraints (11b) in $\omega$ follow from convexity of
the function $(x,y)\mapsto x\log(x/y)$. Furthermore if $H$ is a local
Hamiltonian, and $a$ is a local operator, then the constraints (11) only
involve local observables of $\omega$. Note that for $\beta=+\infty$ (zero
temperature), the constraint (11b) reduces to the inequality (7) which
characterizes the ground states of $H$.

######  Remark 3.3 (Relation to thermodynamic limit).

Using conditions (11) one can show that any limit of finite Gibbs states
$\exp(-\beta H_{\Lambda})/\operatorname{tr}\exp(-\beta H_{\Lambda})$ as
$\Lambda\uparrow\mathbb{Z}^{D}$ will be a valid $\beta$-KMS state. Actually
the same is true for any sequence of finite Gibbs states of
$H_{\Lambda}+\Delta_{\Lambda^{c}}$ where $\Delta_{\Lambda^{c}}$ is an
arbitrary term acting only on the boundary of $\Lambda$.

Similarly to the example for ground-states in Remark 2.1, we demonstrate how
different thermal states can arise from different choices of boundary
conditions. Consider the 2D Ising model at some temperature $0<T<T_{c}$, and
let for any $\ell\in\mathbb{N}$, $\rho_{\ell,\beta}$ be the finite Gibbs state
on a region $\\{-\ell,\dots,\ell\\}^{2}$ with open boundary condition. Then by
the spin-flip symmetry of the Ising Hamiltonian, we know that
$\operatorname{tr}(Z_{x}\rho_{\ell,\beta})=0$ for any site $x$ in the region
$\\{-\ell,\dots,\ell\\}^{2}$. Thus this means that for any limit
$\omega_{\beta}$ of the finite Gibbs state with open boundary condition, we
will have $\omega_{\beta}(Z_{x})=0$. It is known however that for the 2D Ising
model below the critical temperature, there are two (extremal) infinite-volume
Gibbs states $\omega_{\beta}^{+}$ and $\omega_{\beta}^{-}$ such that
$\omega_{\beta}^{-}(Z_{x})<0<\omega_{\beta}^{+}(Z_{x})$, which can be reached
by taking limits of finite Hamiltonians with fixed boundary condition, i.e.,
aligned spins on the boundary [FV17].

###  3.2 Convex relaxations

The EEB inequalities above can be used to obtain rigorous lower and upper
bounds on the value of any observable on thermal states via convex
optimization. For any integer $\ell\geq 0$, let

$\Lambda_{\ell}=\\{-\ell,\ldots,\ell\\}^{D}\subset\mathbb{Z}^{D}$

be a finite sublattice of $\Lambda$, and let
$\mathcal{A}_{\ell}\subset\mathcal{A}$ be the subalgebra of observables
supported on $\Lambda_{\ell}$. We assume for simplicity that $H$ has only
nearest-neighbour interactions so that for $a\in\mathcal{A}_{\ell-1}$,
$[H,a]\in\mathcal{A}_{\ell}$. Let $O$ be a local observable whose support is
contained in $\Lambda_{\ell}$, and consider the optimization problem

$\displaystyle\underset{\tilde{\omega}:\mathcal{A}_{\ell}\to\mathbb{C}}{\text{min/max}}\quad$
$\displaystyle\tilde{\omega}(O)$ (12a) s.t. $\displaystyle\tilde{\omega}(1)=1$
(12b) $\displaystyle\tilde{\omega}(a^{*}a)\geq 0\quad\forall
a\in\mathcal{A}_{\ell}$ (12c)
$\displaystyle\tilde{\omega}(\tau_{x}(a))=\tilde{\omega}(a)\quad\forall
a\in\mathcal{A}_{\ell-1}\;\forall x\in\\{-1,0,+1\\}^{D}$ (12d)
$\displaystyle\tilde{\omega}([H,a])=0\;\forall a\in\mathcal{A}_{\ell-1}$ (12e)
$\displaystyle\tilde{\omega}(a^{*}a)\log\left(\frac{\tilde{\omega}(a^{*}a)}{\tilde{\omega}(aa^{*})}\right)\leq\beta\tilde{\omega}(a^{*}[H,a])\quad\forall
a\in\mathcal{A}_{\ell-1}.$ (12f)

The feasible set of this convex optimization problem is an _outer relaxation_
for the set of $\Lambda_{\ell}$-marginals of translation-invariant thermal
equilibrium states $\mathcal{G}_{\beta}$. This relaxation is obtained by
imposing the EEB inequalities on the finite-dimensional subspace
$\mathcal{A}_{\ell-1}$, in addition to the normalization, positivity, and
translation-invariance constraints (12b), (12c), (12d). Besides the linear
equations (12b), (12d), (12e), the constraint (12c) can be encoded as the
positivity of a Hermitian matrix of size
$|\Lambda_{\ell}|\times|\Lambda_{\ell}|$. However, the last set of constraints
(12f) form an infinite set of scalar inequalities, and it is unclear how to
formulate them in a compact way. Clearly one can sample certain observables
$a_{i}$ ($i=1,\ldots,m$) and impose the $m$ scalar inequalities, however, this
poses the question of how to choose these $a_{i}$’s.

#### Matrix EEB inequality

We introduce the following matrix version of the EEB inequality which allows
us to impose the EEB inequality for all observables $a$ in a finite-
dimensional subspace. The matrix EEB inequality makes use of the matrix
relative entropy function [FK89, FFS90] which is defined as

$D_{op}(A\|B)=A^{1/2}\log(A^{1/2}B^{-1}A^{1/2})A^{1/2}$

for any two $m\times m$ positive semidefinite matrices $A,B\succeq 0$ such
that $A\ll B$ (i.e., $\text{im}\,A\subset\text{im}\,B$). The function $D_{op}$
is the operator perspective [Eff09, ENG11] of the negative logarithm function
(which is operator convex) and satisfies the following properties, see
[FFS90]:

  * •

Homogeneity: $D_{op}(\lambda A\|\lambda B)=\lambda D_{op}(A\|B)$ for
$\lambda>0$

  * •

Operator convexity555For homogeneous functions, convexity is equivalent to
subadditivity:
$D_{op}(\sum_{k}A_{k}\|\sum_{k}B_{k})\preceq\sum_{k}D_{op}(A_{k}\|B_{k})$

  * •

Transformer inequality: For any rectangular matrix $K$,

$D_{op}(K^{*}AK\|K^{*}BK)\preceq K^{*}D_{op}(A\|B)K.$ (13)

Furthermore, the function $D_{op}$ is closed, in the sense that
$\\{(A,B,T):A\ll B\text{ and }D_{op}(A\|B)\preceq T\\}$ is a closed convex set
[FS23, Theorem B.1]. Convex optimization problems involving the matrix
relative entropy can be solved using interior-point methods [FS23] or via
semidefinite programming approximations [FSP19].

######  Theorem 3.4 (Matrix EEB inequality).

Let $H$ be a finite-range Hamiltonian and $\omega:\mathcal{A}\to\mathbb{C}$ be
a $\beta$-KMS state. Let $a_{1},\ldots,a_{m}\in\mathcal{A}_{\text{loc}}$ and
define the $m\times m$ matrices

$\displaystyle A_{ij}$ $\displaystyle=\omega(a_{i}^{*}a_{j})$ (14)
$\displaystyle B_{ij}$ $\displaystyle=\omega(a_{j}a_{i}^{*})$ $\displaystyle
C_{ij}$ $\displaystyle=\omega(a_{i}^{*}[H,a_{j}]).$

Then $A,B,C$ are Hermitian, $A,B$ are positive semidefinite with $A\ll B$ and
$D_{op}(A\|B)\preceq\beta C$.

###### Proof.

It is clear that $A,B$ are Hermitian and that they are positive semidefinite.
That $A\ll B$ follows from the fact that for a KMS state
$\omega(bb^{*})=0\implies\omega(b^{*}b)=0$ for any $b\in\mathcal{A}$. (This
can be seen e.g., from (11b)). To see that $C$ is Hermitian, note that since
$\omega$ is a KMS state we have $\omega([H,b])=0$ for any
$b\in\mathcal{A}_{\text{loc}}$ (see (11a)). Then

$C_{ij}=\omega(a_{i}^{*}[H,a_{j}])=\omega(a_{i}^{*}(Ha_{j}-a_{j}H))=\omega(a_{i}^{*}Ha_{j}-Ha_{i}^{*}a_{j})=\omega([a_{i}^{*},H]a_{j})=\overline{\omega(a_{j}^{*}[H,a_{i}])}=\overline{C_{ji}}.$

We now prove the inequality $D_{op}(A\|B)\preceq\beta C$. We present the proof
here in the finite-dimensional case. We treat the general case in Appendix A.
We thus assume that $\mathcal{A}$ is a finite-dimensional matrix algebra. Let
$L:\mathcal{A}\to\mathcal{A}$ be the commutator operator with the Hamiltonian
$H$, i.e., $L(a)=[H,a]$. The map $L$ is self-adjoint with respect to the
Hilbert-Schmidt inner product. Observe that for any $a\in\mathcal{A}$,
$e^{-\beta L}(a)=e^{-\beta H}ae^{\beta H}$ (This can be easily seen e.g., by
checking that $a(\beta):=e^{-\beta H}ae^{\beta H}$ is the solution of the ODE
$a^{\prime}(\beta)=-L(a(\beta))$ with $a(0)=a$.) Let
$L=\sum_{k}\lambda_{k}P_{k}$ be a spectral decomposition of $L$. By the KMS
condition we have

$B_{ij}=\omega(a_{j}a_{i}^{*})=\omega(a_{i}^{*}e^{-\beta
L}(a_{j}))=\sum_{k}e^{-\beta\lambda_{k}}\omega(a_{i}^{*}P_{k}(a_{j}))=\sum_{k}e^{-\beta\lambda_{k}}A^{(k)}_{ij}$

where we let $A^{(k)}_{ij}=\omega(a_{i}^{*}P_{k}(a_{j}))$. Since
$\sum_{k}P_{k}=I$ and $\sum_{k}\lambda_{k}P_{k}=L$, we have

$\sum_{k}A^{(k)}=A\quad\text{ and }\quad\sum_{k}\lambda_{k}A^{(k)}=C.$

Then we have

$\displaystyle
D_{op}(A\|B)=D_{op}\left(\sum_{k}A^{(k)}\|\sum_{k}e^{-\beta\lambda_{k}}A^{(k)}\right)\preceq\sum_{k}D_{op}(A^{(k)}\|e^{-\beta\lambda_{k}}A^{(k)})=\sum_{k}\beta\lambda_{k}A^{(k)}=\beta
C$

as desired. ∎

#### A convex relaxation using the matrix relative entropy

We now consider the modified convex optimization problem

$\displaystyle\underset{\tilde{\omega}:\mathcal{A}_{\ell}\to\mathbb{C}}{\text{min/max}}\quad$
$\displaystyle\tilde{\omega}(O)$ (15a) s.t. $\displaystyle\tilde{\omega}(1)=1$
(15b) $\displaystyle\tilde{\omega}(a^{*}a)\geq 0\quad\forall
a\in\mathcal{A}_{\ell}$ (15c)
$\displaystyle\tilde{\omega}(\tau_{x}(a))=\tilde{\omega}(a)\quad\forall
a\in\mathcal{A}_{\ell-1}\;\forall x\in\\{-1,0,+1\\}^{D}$ (15d)
$\displaystyle\tilde{\omega}([H,a])=0\;\forall a\in\mathcal{A}_{\ell-1}$ (15e)
$\displaystyle
D_{op}\left(\,[\tilde{\omega}(a_{i}^{*}a_{j})]_{ij}\,,\,[\tilde{\omega}(a_{j}a_{i}^{*})]_{ij}\,\right)\;\preceq\;\beta[\tilde{\omega}(a_{i}^{*}[H,a_{j}])]_{ij}.$
(15f)

where in the last constraint, $\\{a_{i}\\}$ is a basis of
$\mathcal{A}_{\ell-1}$. One can show, using (13), that the infinite number of
scalar constraints (12f) are implied by the finite positive semidefinite
constraint (15f). We prove this in the next proposition.

######  Proposition 3.5.

Consider the optimization problem (12) and assume $\\{a_{i}\\}$
($i=1,\ldots,m$) is a basis of $\mathcal{A}_{\ell-1}$. Then the constraints
(12f) are all implied by the matrix inequality (15f).

###### Proof.

Let $a$ be an arbitrary element of $\mathcal{A}_{\ell-1}$, and let
$a=\sum_{i}x_{i}a_{i}$, with $x_{i}\in\mathbb{C}$, be its decomposition in the
basis $\\{a_{i}\\}$. Also let $A,B,C$ be the matrices defined in (14), so that
the inequality (15f) can be written as $D_{op}(A\|B)\preceq\beta C$. If we let
$x=(x_{1},\ldots,x_{m})\in\mathbb{C}^{m}$, then note that

$\displaystyle x^{*}Ax$
$\displaystyle=\sum_{ij}\overline{x_{i}}x_{j}\omega(a_{i}^{*}a_{j})=\omega(a^{*}a)$
$\displaystyle x^{*}Bx$
$\displaystyle=\sum_{ij}\overline{x_{i}}x_{j}\omega(a_{j}a_{i}^{*})=\omega(aa^{*})$
$\displaystyle x^{*}Cx$
$\displaystyle=\sum_{ij}\overline{x_{i}}x_{j}\omega(a_{i}^{*}[H,a_{j}])=\omega(a^{*}[H,a]).$

It follows from (13) that

$\omega(a^{*}a)\log\left(\frac{\omega(a^{*}a)}{\omega(aa^{*})}\right)=D_{op}(x^{*}Ax\|x^{*}Bx)\leq
x^{*}D_{op}(A\|B)x\leq\beta x^{*}Cx=\beta\omega(a^{*}[H,a])$

as desired. ∎

###  3.3 Convergence

Let $(\mathsf{P}^{\min}_{\ell})$ and $(\mathsf{P}^{\max}_{\ell})$ be
respectively the minimization and maximization problems (15), and let
$\mathsf{p}^{\min}_{\ell}\leq\mathsf{p}^{\max}_{\ell}$ be their optimal
values. The next theorem shows the asymptotic convergence of
$\mathsf{p}^{\min}_{\ell}$ and $\mathsf{p}^{\max}_{\ell}$ to $\langle
O\rangle^{\min}_{\beta}$ and $\langle O\rangle^{\max}_{\beta}$ respectively.

######  Theorem 3.6.

Let $H$ be a nearest-neighbour Hamiltonian on $\mathbb{Z}^{D}$,
$O\in\mathcal{A}_{\text{loc}}$ a local observable, and $\beta\in[0,+\infty]$.
Define

$[\langle O\rangle^{\min}_{\beta},\langle
O\rangle^{\max}_{\beta}]=\\{\omega(O):\omega\in\mathcal{G}_{\beta}\\}.$

Then $\mathsf{p}^{\min}_{\ell}\uparrow\langle O\rangle^{\min}_{\beta}$ and
$\mathsf{p}^{\max}_{\ell}\downarrow\langle O\rangle^{\max}_{\beta}$ as
$\ell\to\infty$.

###### Proof.

Note that $\mathcal{G}_{\beta}$ is convex and weak-* compact [BR97, Theorem
5.3.30] and so $\\{\omega(O):\omega\in\mathcal{G}_{\beta}\\}$ is a closed
interval.

If $\omega\in\mathcal{G}_{\beta}$, then for any $\ell$, its restriction to
$\mathcal{A}_{\ell}$ is feasible for the optimization problems (15). Thus this
means that $\mathsf{p}^{\min}_{\ell}\leq\langle
O\rangle^{\min}_{\beta}\leq\langle
O\rangle^{\max}_{\beta}\leq\mathsf{p}^{\max}_{\ell}$ for any $\ell$.
Furthermore, observe that $\mathsf{p}^{\min}_{\ell}$ is monotonic
nondecreasing, and $\mathsf{p}^{\max}_{\ell}$ is monotone nonincreasing.
Indeed assume $\tilde{\omega}_{\ell}:\mathcal{A}_{\ell}\to\mathbb{C}$ is
feasible for (15) and consider its restriction to $\mathcal{A}_{k}$ for
$k<\ell$. We want to show that it is feasible for the level $k$ of the
relaxation (15). It is immediate that the restriction satisfies the
constraints (15b)-(15e) at the level $k$ of the relaxation. To prove (15f) we
use the transformer inequality (13). Let $\\{b_{s}\\}$ be a basis of
$\mathcal{A}_{k-1}$ and note that we can express each $b_{s}$ as a linear
combination of the $a_{i}\in\mathcal{A}_{\ell-1}$, i.e.,
$b_{s}=\sum_{i}K_{si}a_{i}$ for some $K_{si}\in\mathbb{C}$. As such we have

$\tilde{\omega}(b_{s}^{*}b_{t})=\sum_{ij}\overline{K_{si}}K_{tj}\tilde{\omega}(a_{i}^{*}a_{j})$

so that we can write
$[\tilde{\omega}(b_{s}^{*}b_{t})]_{s,t}=\overline{K}[\tilde{\omega}(a_{i}^{*}a_{j})]_{i,j}\overline{K}^{*}$
and similarly for $[\tilde{\omega}(b_{t}b_{s}^{*})]_{s,t}$. Hence we get

$\displaystyle
D_{op}([\tilde{\omega}(b_{s}^{*}b_{t})]_{s,t}\;,\;[\tilde{\omega}(b_{t}b_{s}^{*})]_{s,t})$
$\displaystyle=D_{op}(\overline{K}[\tilde{\omega}(a_{i}^{*}a_{j})]_{i,j}\overline{K}^{*}\;,\overline{K}[\tilde{\omega}(a_{j}a_{i}^{*})]_{i,j}\overline{K}^{*})$
$\displaystyle\preceq\overline{K}D_{op}([\tilde{\omega}(a_{i}^{*}a_{j})]_{i,j}\;,\;[\tilde{\omega}(a_{j}a_{i}^{*})]_{i,j})\overline{K}^{*}$
$\displaystyle\preceq\beta\overline{K}[\tilde{\omega}(a_{i}^{*}[H,a_{j}])]_{i,j}\overline{K}^{*}$
$\displaystyle=[\tilde{\omega}(b_{s}^{*}[H,b_{t}])]_{s,t},$

as desired.

Since $(\mathsf{p}^{\min}_{\ell})$ and $(\mathsf{p}^{\max}_{\ell})$ are
monotonic they must have limits
$\mathsf{p}^{\min}_{\infty}:=\lim\mathsf{p}^{\min}_{\ell}$ and
$\mathsf{p}^{\max}_{\infty}:=\lim\mathsf{p}^{\max}_{\ell}$ respectively. We
want to show that $\mathsf{p}^{\min}_{\infty}=\langle O\rangle^{\min}_{\beta}$
and $\mathsf{p}^{\max}_{\infty}=\langle O\rangle^{\max}_{\beta}$.

For each $\ell$, let $\tilde{\omega}_{\ell}:\mathcal{A}_{\ell}\to\mathbb{C}$
be a feasible solution to (15) such that
$\mathsf{p}^{\min}_{\ell}\leq\tilde{\omega}_{\ell}(O)\leq\mathsf{p}^{\min}_{\ell}+1/\ell$
(if the minimum value of the optimization problem is attained, we simply take
$\tilde{\omega}$ the optimal solution which satisfies
$\tilde{\omega}_{\ell}(O)=\mathsf{p}^{\min}_{\ell}$). We can extend
$\tilde{\omega}_{\ell}$ to a state on the full $C^{*}$-algebra $\mathcal{A}$
[BR87, Proposition 2.3.24]. We now have a sequence of states
$\tilde{\omega}_{\ell}$ on $\mathcal{A}$. Since the set of states on
$\mathcal{A}$ is weak-* compact [BR87, Theorem 2.3.15] we can extract from
$\tilde{\omega}_{\ell}$ a subsequence $\tilde{\omega}_{n(\ell)}$ that weak-*
converges to a state $\tilde{\omega}$ on $\mathcal{A}$. We now show that
$\tilde{\omega}$ is a $\beta$-KMS translation-invariant state.

To prove translation-invariance, let $a\in\mathcal{A}_{\text{loc}}$ and
$x\in\mathbb{Z}^{D}$. Since $a$ is a local operator, we know that
$a,\tau_{x}(a)\in\mathcal{A}_{n(\ell)}$ for all large enough $\ell$. This
means that

$\tilde{\omega}(\tau_{x}(a))=\lim_{\ell}\tilde{\omega}_{n(\ell)}(\tau_{x}(a))=\lim_{\ell}\tilde{\omega}_{n(\ell)}(a)=\tilde{\omega}(a)$

where the middle inequality follows from the translation-invariant constraint
imposed in the optimization problem (15). This is true for all
$a\in\mathcal{A}_{\text{loc}}$ and $x\in\mathbb{Z}^{D}$, and so
$\tilde{\omega}$ is translation-invariant.

We now show that $\tilde{\omega}$ is a $\beta$-KMS state by showing that the
conditions (11) hold. Let $a\in\mathcal{A}_{\text{loc}}$. Since $a$ is a local
observable we have $a\in\mathcal{A}_{n(\ell)-1}$ for all large enough $\ell$.
Hence this implies that

$\displaystyle\tilde{\omega}_{n(\ell)}([H,a])=0$
$\displaystyle\tilde{\omega}_{n(\ell)}(a^{*}a)\log\left(\frac{\tilde{\omega}_{n(\ell)}(a^{*}a)}{\tilde{\omega}_{n(\ell)}(aa^{*})}\right)\leq\beta\tilde{\omega}_{n(\ell)}(a^{*}[H,a])$

for all large enough $\ell$. By taking the limit $\ell\to\infty$ this shows
that $\tilde{\omega}$ satisfies the conditions (11) and is thus a $\beta$-KMS
state. Finally note that

$\tilde{\omega}(O)=\lim_{\ell}\tilde{\omega}_{n(\ell)}(O)\leq\lim_{\ell}\mathsf{p}^{\min}_{n(\ell)}+1/n(\ell)=\mathsf{p}^{\min}_{\infty}.$

Since $\mathsf{p}^{\min}_{\infty}\leq\langle
O\rangle^{\min}_{\beta}\leq\tilde{\omega}(O)$ we necessarily get equality
$\mathsf{p}^{\min}_{\infty}=\langle O\rangle^{\min}_{\beta}$. The same
argument can be used to show that $\mathsf{p}^{\max}_{\infty}=\langle
O\rangle^{\max}_{\beta}$. ∎

###  3.4 Decidability

It has been shown in [BCW21] that computing the phase diagram of a 2D nearest
neighbour translation-invariant Hamiltonian in the case $\beta=\infty$ is
undecidable. This is shown for both a definition of phase in terms of gapped
and gapless Hamiltonians as well as in terms of an order parameter. It is
interesting to compare the latter with our result. More precisely, the
construction of the authors gives a Hamiltonian (described by a nearest
neighbour interaction $h(n)$) for every input $n$ to a universal Turing
machine and a local observable $O$ such that

  * •

If the Turing machine halts on input $n$ then for any sequence of ground
states $|\psi_{L}\rangle$ of the Hamiltonian on an $L\times L$ square with
open boundary conditions
$\lim_{L\to\infty}\langle\psi_{L}|O|\psi_{L}\rangle=0$.

  * •

If the Turing machine does not halt on input $n$, then for any sequence of
ground states $|\psi_{L}\rangle$ of the Hamiltonian on an $L\times L$ square
with open boundary conditions
$\lim_{L\to\infty}\langle\psi_{L}|O|\psi_{L}\rangle=1$.

The limit values in this construction however crucially depend on the choice
of boundary conditions. While the boundary condition is not mentioned
explicitly as it is always kept open, in fact a technique from [GI09] is used
to encode a boundary term corresponding to an energy shift of a Hamiltonian
into the translation-invariant interaction term. Before this shift the
construction yields a nonnegative ground-state energy in the nonhalting case
and a negative one in the halting case. Adding this energy shift and adding a
subspace $\\{|0\rangle\\}$ with a trivial product state with zero energy to
the Hamiltonian, the ground state in the nonhalting case will be given by this
product state. The observable $O=|0\rangle\\!\langle 0|$ is then given by a
projector onto this subspace.

The definition of ground-states in our work is different and in fact the set
of equilibrium states $\mathcal{G}_{\infty}$ includes the limits of ground-
states of Hamiltonians for any boundary condition, see Remark 2.1. In this
version, as opposed to the one above, our hierarchies do prove the
decidability of the equilibrium observable problem as defined in Definition
1.2.

######  Theorem 3.7.

The Equilibrium Observable Problem is decidable.

###### Proof.

Let us start with the case $\beta=\infty$. The algorithm consists of iterating
over $\ell=2,3,\ldots$ and alternating between the minimization and
maximization of the convex relaxations (15). Note that for the ground-state
Eq. (15f) reduces to
$0\preceq\;\beta[\tilde{\omega}(a_{i}^{*}[H,a_{j}])]_{ij}$, such that the
problem becomes an SDP in rational coefficients that can be solved exactly.
More specifically, we use quantifier elimination [Ana98, BPR06] to decide
whether the hierarchy augmented by $\tilde{\omega}(O)\geq a+(b-a)/2$ or
augmented by $\tilde{\omega}(O)\leq a+(b-a)/2$ is feasible. The algorithm
stops as soon as one of them becomes infeasible and outputs NO or YES
respectively. Due to Theorem 1.1 this is guaranteed to happen at sufficiently
large $\ell$.

In the case $\beta<\infty$, we face additional difficulties due to the matrix
logarithm, as our convex relaxation is no longer an SDP. We will use SDP
approximations of the logarithm, however, those have convergence guarantees
that are only uniform on intervals that are bounded away from 0. A way to deal
with this issue is to strengthen our convex relaxation by an a priori lower
bound on the moment matrices derived from an analysis of the Hamiltonian. This
is proven in the following Lemma by adapting a result from [BLMT23] to the
infinite system moment problem setting.

######  Lemma 3.8.

For a given local translation-invariant Hamiltonian $H$ on $\mathbb{Z}^{D}$
and temperature $\beta$, there exists a KMS-state
$\omega\in\mathcal{G}_{\beta}$ such that the following is true: for any finite
set $\Lambda\subset\mathbb{Z}^{D}$, there exists an explicit constant
$C_{\Lambda}$ depending on $\beta$, the locality and norm of the Hamiltonian,
and $\Lambda$, such that the marginal $\rho_{\Lambda}\in\mathcal{A}_{\Lambda}$
of $\omega$, defined by

$\operatorname{tr}[\rho_{\Lambda}a]=\omega(a)\qquad\forall
a\in\mathcal{A}_{\Lambda}$

is lower bounded

$\rho_{\Lambda}\succeq C_{\Lambda}.$

###### Proof.

Without loss of generality we assume that the Hamiltonian is defined on qubits
since general qudit systems can be embedded by encoding each qudit in a
Hilbert space $\mathbb{C}^{2^{m^{D}}}$ such that $d\leq 2^{m^{D}}$,
corresponding to a hypercube of spins. This results in a Hamiltonian with
higher but still finite range.

In [BLMT23, Corollary 2.14] it is shown that for a given qubit Hamiltonian
$H_{V}$ on a finite set $V$ and for a fixed $\Lambda\subset V$ there exist
explicit constants $c_{\Lambda},c^{\prime}_{\Lambda}$ such that the following
is true: for any Hermitian observable $X=\sum_{i}x_{i}a_{i}$ supported on
$\Lambda\subset V$ (here, $x_{i}\in\mathbb{R}$ and $a_{i}$ is the orthonormal
Pauli basis supported on $\Lambda$) we have

$\operatorname{tr}[X^{2}\tilde{\rho}_{V}]\geq\exp(-c_{\Lambda}-\beta
c^{\prime}_{\Lambda})\max x_{i}^{2},$

where $\tilde{\rho}_{V}=e^{-\beta H_{V}}/\operatorname{tr}e^{-\beta H_{V}}$ is
the Gibbs state at inverse temperature $\beta$ for the given finite
Hamiltonian $H_{V}$. Crucially, the constants
$c_{\Lambda},c^{\prime}_{\Lambda}$ do not depend on $|V|$, but only on the
locality of the Hamiltonian and the degree of its so-called dual interaction
graph, which are fixed.666The dual interaction graph has nodes for each
Hamiltonian term and edges between terms that overlap If we take the
weak*-limit of Gibbs states $\tilde{\rho}_{V}$ on increasing system sizes
$V\uparrow\mathbb{Z}^{D}$ (taking an appropriate subsequence to ensure
convergence, see [BR97, Proposition 6.2.15]) we get a KMS-state $\omega$ which
satisfies

$\omega\left(X^{2}\right)\geq\exp(-c_{\Lambda}-\beta
c^{\prime}_{\Lambda})\max_{i}x_{i}^{2}$

for all $X\in\mathcal{A}_{\Lambda}$, where the $x_{i}$ are the coefficients of
$X$ in the Pauli basis of $\mathcal{A}_{\Lambda}$.

Now let $\rho_{\Lambda}$ be the marginal of $\omega$ on $\Lambda$, and let $X$
be the projector onto a one-dimensional eigenspace of $\rho_{\Lambda}$ with
lowest eigenvalue. Note that the $x_{i}$ are real as $X$ is a Hermitian
operator and that

$\max_{i}x_{i}^{2}\geq\frac{1}{4^{|\Lambda|}}\|(x_{i})_{i}\|_{2}^{2}=\frac{1}{4^{|\Lambda|}}\|X\|_{F}^{2}=\frac{1}{4^{|\Lambda|}}.$

Thereby, we conclude

$\rho_{\Lambda}\succeq\operatorname{tr}\left[\rho_{\Lambda}X^{2}\right]=\omega\left(X^{2}\right)\geq\frac{1}{4^{|\Lambda|}}\exp(-c_{\Lambda}-\beta
c^{\prime}_{\Lambda})=:C_{\Lambda}.$

∎

We now choose the basis of $\mathcal{A}_{\Lambda}$ in our convex relaxation
(15) to be orthonormal. We see that the above Lemma proves a lower bound on
the lowest eigenvalue of the moment matrix $\omega$ that corresponds to the
KMS-state given in the Lemma:

$\min_{\|x\|_{2}^{2}=1}\langle
x|[\tilde{\omega}(a_{i}^{*}a_{j})]_{ij}|x\rangle=\min_{\|x\|_{2}^{2}=1}\sum_{ij}\overline{x_{i}}x_{j}\omega(a_{i}^{*}a_{j})=\operatorname{tr}[\rho
X^{*}X]\geq C_{\Lambda_{\ell}}\operatorname{tr}[X^{*}X]=C_{\Lambda_{\ell}}$
(16)

and similarly for $[\tilde{\omega}(a_{j}a_{i}^{*})]_{ij}$. We use this to
strengthen our convex relaxation (15) by adding the contraints

ij $\displaystyle\succeq C_{\Lambda_{\ell}}^{\prime}>0,$ (17)
$\displaystyle[\tilde{\omega}(a_{j}a_{i}^{*})]_{ij}$ $\displaystyle\succeq
C_{\Lambda_{\ell}}^{\prime}>0$

where $C^{\prime}_{\Lambda_{\ell}}$ is a rational lower bound to
$C_{\Lambda_{\ell}}$. Note that Lemma 3.8 need not hold for KMS states that
arise from different boundary conditions as those might come with different
constants $c$, $c^{\prime}$. In the promise problem setting this is not
relevant though as the relaxation still contains at least the open boundary
Gibbs state and thereby the respective SDP problem is always feasible.
Therefore, the optimal values of the modified problems are two sequences
$\mathsf{p}^{\min}_{\ell}\leq\mathsf{p}^{{}^{\prime}\min}_{\ell}\leq\mathsf{p}^{{}^{\prime}\max}_{\ell}\leq\mathsf{p}^{\max}_{\ell}$.
Note that monotonicity and thereby convergence are no longer guaranteed for
$\mathsf{p}^{{}^{\prime}\min}_{\ell},\mathsf{p}^{{}^{\prime}\max}_{\ell}$.
However, we only need that the optimal values exceed the decision threshold,
which is guaranteed by convergence of
$\mathsf{p}^{\min}_{\ell},\mathsf{p}^{\max}_{\ell}$ and the above inequality.

After this modification, since the moment matrix is bounded away from zero and
upper bounded by 1, the eigenvalues in the argument of the logarithm for the
operator relative entropy are contained in an interval $[a,1/a]$ for some
$a>0$. We now introduce a semidefinite approximation of $D_{op}$ from [FSP19].
For any integer $m\geq 1$ let $h_{m}(x)=2^{m}\left(x^{1/2^{m}}-1\right)$ which
satisfy $\log(x)\leq h_{m+1}(x)\leq h_{m}(x)$ and note that $h_{m}\to\log$
uniformly on any compact interval in $(0,\infty)$ since $0\leq
h_{m}(x)-\log(x)\leq 2^{-m}[(x-1)^{2}+(x^{-1}-1)^{2}]$ for all $x>0$. Consider
the _noncommutative perspective_ [Eff09, ENG11] of $-h_{m}$ which we denote by
$D_{op}^{[m]}$:

$D_{op}^{[m]}(X,Y):=P_{-h_{m}}(X,Y)=-X^{1/2}h_{m}(X^{-1/2}YX^{-1/2})X^{1/2}.$

Then $D_{op}^{[m]}$ is jointly convex in $(X,Y)$, has a finite semidefinite
programming formulation [FSP19], and satisfies $D_{op}^{[m]}(X\|Y)\preceq
D_{op}(X\|Y)$ since $-h_{m}\leq-\log$. By changing the function $D_{op}$ in
the entropy constraint (15f) by $D_{op}^{[m]}$ we get a semidefinite program
$(\mathsf{P}^{{}^{\prime}\min}_{\ell,m})$ with rational coefficients whose
optimal value satisfies
$\mathsf{p}^{{}^{\prime}\min}_{\ell,m}\leq\mathsf{p}^{{}^{\prime}\min}_{\ell}$
and analogously for the upper bounds.

Furthermore using the constraint (17), the uniform convergence of $h_{m}$ to
$\log$, and the continuity of $D_{op}$ on positive definite matrices, we know
that
$\mathsf{p}^{{}^{\prime}\min}_{\ell,m}\uparrow\mathsf{p}^{{}^{\prime}\min}_{\ell}$
and similarly
$\mathsf{p}^{{}^{\prime}\max}_{\ell,m}\downarrow\mathsf{p}^{{}^{\prime}\max}_{\ell}$
as $m\uparrow\infty$. By choosing an exhaustive sequence $(\ell_{i},m_{i})$ of
$\mathbb{N}^{+}\times\mathbb{N}^{+}$ and deciding the two feasibility problems
as in the ground-state case we obtain an algorithm that decides the problem in
finite time.

∎

###  3.5 Commuting Hamiltonians and the high-temperature regime

In this section we focus on commuting Hamiltonians $H$ which satisfy
$[h_{X},h_{Y}]=0$ for any $X,Y\subset\mathbb{Z}^{D}$. For such Hamiltonians
one can define a convex relaxation for the set of thermal equilibrium states
directly from the KMS equality condition. Indeed, in this case if
$b\in\mathcal{A}_{\Lambda}\subset\mathcal{A}_{\text{loc}}$ where $\Lambda$ is
finite, then

$\alpha_{t}(b)=\lim_{\Lambda^{\prime}\uparrow\mathbb{Z}^{D}}e^{itH_{\Lambda^{\prime}}}be^{-itH_{\Lambda^{\prime}}}=e^{it\tilde{H}_{\Lambda}}be^{-it\tilde{H}_{\Lambda}}$
(18)

where $\tilde{H}_{\Lambda}$ is as defined in (6). Clearly, in this case
$\alpha_{t}(b)$ can be extended to an entire function for all
$b\in\mathcal{A}_{\text{loc}}$. As such, the $\beta$-KMS condition for a state
$\omega$ can be written as

$\omega(ba)=\omega(ae^{-\beta\tilde{H}_{\Lambda}}be^{\beta\tilde{H}_{\Lambda}})\qquad\forall
a\in\mathcal{A}_{\bar{\Lambda}},\;b\in\mathcal{A}_{\Lambda}$ (19)

for all $\Lambda\subset\mathbb{Z}^{D}$ finite. (Recall that
$\bar{\Lambda}=\Lambda\cup\partial^{ex}\Lambda$ where $\partial^{ex}\Lambda$
is the external boundary of $\Lambda$.) A consequence of these KMS
conditions777This is because $\alpha_{t}(h_{X})=h_{X}$, since the Hamiltonian
is commuting. is that

$\omega([a,h_{X}])=0\;\;\forall a\in\mathcal{A}_{\bar{\Lambda}},\;\forall
X\subset\bar{\Lambda}$ (20)

for all $\Lambda\subset\mathbb{Z}^{D}$ finite. Given a fixed
$\Lambda\subset\mathbb{Z}^{D}$ finite, let888Note that (19) for a fixed
$\Lambda$ only implies (20) for $X\subset\Lambda$ (and not $\bar{\Lambda}$).
That’s why we impose (20) explicitly.

$\displaystyle\mathcal{G}_{\Lambda,\beta}^{\text{comm}}=\Bigl{\\{}\omega:\mathcal{A}_{\bar{\Lambda}}\to\mathbb{C}\text{
s.t. }$ $\displaystyle\omega(1)=1,\omega(a^{*}a)\geq 0\;\;\forall
a\in\mathcal{A}_{\bar{\Lambda}},$ (21)
$\displaystyle\omega(a)=\omega(\tau_{x}(a))\;\forall
a\in\mathcal{A}_{\Lambda},x\in\\{-1,0,1\\}^{D},\text{ and
}\eqref{eq:localKMScomm}\text{ and }\eqref{eq:localKMScomm2}\Bigr{\\}}.$

Since (19) is bilinear in $(a,b)$ it only needs to be imposed on a basis
$\\{a_{i}\\}$ of $\mathcal{A}_{\bar{\Lambda}}$ and $\\{b_{j}\\}$ of
$\mathcal{A}_{\Lambda}$. As such, the convex set (21) is the feasible set of a
semidefinite program. It is an outer relaxation for the set of
$\bar{\Lambda}$-marginals of $\beta$-KMS states.

Given an observable $O\in\mathcal{A}_{\Lambda}$, one can consider the convex
optimization problems

$\min/\max\;\\{\tilde{\omega}(O)\text{ s.t.
}\tilde{\omega}\in\mathcal{G}_{\Lambda,\beta}^{\text{comm}}\\}.$ (22)

The solution of this minimization (resp. maximization) problem is a _lower
bound_ on $\langle O\rangle^{\min}_{\beta}$ (resp. _upper bound_ on $\langle
O\rangle^{\max}_{\beta}$).

The asymptotic convergence of the hierarchy (22) for a sequence of increasing
finite lattices $\Lambda_{\ell}$ to $\langle O\rangle^{\min}_{\beta}$ and
$\langle O\rangle^{\max}_{\beta}$ can be proved in exactly the same way as in
Theorem 3.6. For this hierarchy however we are also able to give a
quantitative rate of convergence in the high temperature regime.

######  Theorem 3.9.

Assume that the Hamiltonian is translation-invariant, commuting and that the
lattice dimension $D\leq 2$. For $D=1$ and $\beta_{1}=\infty$ or for $D=2$ and
some $\beta_{1}>0$, we have for all $0\leq\beta<\beta_{1}$, and for any local
observable $O$ the following:

  * •

$\langle O\rangle^{\min}_{\beta}=\langle O\rangle^{\max}_{\beta}=\langle
O\rangle_{\beta}$

  * •

If we let $\mathsf{p}^{\min}_{\ell}$ and $\mathsf{p}^{\max}_{\ell}$ be the
optimal values of the convex relaxations (22) with
$\Lambda=\Lambda_{\ell}=\\{-\ell,\ldots,\ell\\}^{D}$, then
$\mathsf{p}^{\min}_{\ell}\leq\langle
O\rangle_{\beta}\leq\mathsf{p}^{\max}_{\ell}$, and furthermore

$\max\left\\{\langle
O\rangle_{\beta}-\mathsf{p}^{\min}_{\ell},\mathsf{p}^{\max}_{\ell}-\langle
O\rangle_{\beta}\right\\}\leq c_{1}\|O\|e^{-c_{2}\ell}.$

for some constants $c_{1},c_{2}>0$ depending on the interaction, temperature
and the support of $O$.

###### Proof.

We start with an important lemma showing that states in
$\mathcal{G}_{\Lambda,\beta}^{\text{comm}}$ can be expressed as a product of
the local (i.e., finite) Gibbs state on $\Lambda$ and a boundary term on
$\partial^{ex}\Lambda$, after a suitable perturbation that eliminates the
surface interaction terms. The statement is closely related to [BR97, Prop.
6.2.17]. The difference is that the state $\omega$ in the lemma below is not
assumed to be a “full” KMS state, but just a state in
$\mathcal{G}_{\Lambda,\beta}^{\text{comm}}$.

######  Lemma 3.10.

Assume $H$ is a commuting Hamiltonian. Let $\Lambda\subset\mathbb{Z}^{D}$ be a
finite set, and let $\omega\in\mathcal{G}_{\Lambda,\beta}^{\text{comm}}$. Let

$W=\sum_{\begin{subarray}{c}X\cap\Lambda\neq\emptyset\\\
X\cap\Lambda^{c}\neq\emptyset\end{subarray}}h_{X}\in\mathcal{A}_{\bar{\Lambda}}$

be the surface interaction term. Let

$\varphi^{\Lambda}(x)=\frac{\operatorname{tr}(e^{-\beta
H_{\Lambda}}x)}{\operatorname{tr}e^{-\beta H_{\Lambda}}}.$

There exists a state
$\hat{\omega}:\mathcal{A}_{\partial^{ex}\Lambda}\to\mathbb{C}$ on the boundary
such that using the product state
$\omega^{W}:\mathcal{A}_{\bar{\Lambda}}\to\mathbb{C}$ defined by

$\omega^{W}(xc)=\varphi^{\Lambda}(x)\hat{\omega}(c)$

we can express $\omega$ as a perturbation of the local Gibbs state with a
boundary condition

$\omega(O)=\frac{\omega^{W}(Oe^{-\beta W})}{\omega^{W}(e^{-\beta W})}.$

###### Proof.

We define

$\omega^{W}(a)=\frac{\omega(e^{\beta W}a)}{\omega(e^{\beta
W})}=\frac{\omega(ae^{\beta W})}{\omega(e^{\beta W})},$

where the equality follows using the $\beta$-KMS condition and the fact that
$H$ is commuting. We first prove that it has the claimed form, i.e., is
product with the factor on $\Lambda$ being $\varphi^{\Lambda}$. Let
$a,b\in\mathcal{A}_{\Lambda}$ and $c\in\mathcal{A}_{\partial^{ex}\Lambda}$.
Then condition (19) tells us that

$\displaystyle\omega^{W}(bac)$ $\displaystyle=\frac{1}{\omega(e^{\beta
W})}\omega(bace^{\beta W})$ (23) $\displaystyle=\frac{1}{\omega(e^{\beta
W})}\omega(ace^{\beta
W}e^{-\beta\tilde{H}_{\Lambda}}be^{\beta\tilde{H}_{\Lambda}})$
$\displaystyle=\frac{1}{\omega(e^{\beta W})}\omega(ace^{-\beta
H_{\Lambda}}be^{\beta H_{\Lambda}}e^{\beta W})$
$\displaystyle=\omega^{W}(ace^{-\beta H_{\Lambda}}be^{\beta H_{\Lambda}})$
$\displaystyle=\omega^{W}(ae^{-\beta H_{\Lambda}}be^{\beta H_{\Lambda}}c)$

where in the last line we used the fact that
$c\in\mathcal{A}_{\partial^{ex}\Lambda}$ and so it commutes with $e^{-\beta
H_{\Lambda}}be^{\beta H_{\Lambda}}\in\mathcal{A}_{\Lambda}$.

From (23) we get that for any fixed $c\in\mathcal{A}_{\partial^{ex}\Lambda}$,
the state $x\in\mathcal{A}_{\Lambda}\mapsto\omega^{W}(xc)/\omega^{W}(c)$
satisfies the KMS condition for the finite-dimensional Hamiltonian
$H_{\Lambda}$. Since this is a finite system, this means that necessarily for
any $x\in\mathcal{A}_{\Lambda}$,

$\frac{\omega^{W}(xc)}{\omega^{W}(c)}=\frac{\operatorname{tr}(e^{-H_{\Lambda}}x)}{\operatorname{tr}e^{-H_{\Lambda}}}.$

This proves that $\omega^{W}$ is product and has the marginal on $\Lambda$ as
claimed. We conclude by computing by

$\omega(O)=\frac{\omega(O)}{\omega(e^{\beta W})}\frac{\omega(e^{\beta
W})}{\omega(I)}=\frac{\omega^{W}(Oe^{-\beta W})}{\omega^{W}(e^{-\beta W})},$

using the definition of $\omega^{W}$ twice. This proves Lemma 3.10. ∎

We are now ready to prove exponential convergence of the convex optimization
hierarchy to $\langle O\rangle_{\beta}$ for $0\leq\beta<\beta_{1}$. We start
with the 1D case. Our main tool will be the decay of correlations and local
indistinguishability of Gibbs states for one-dimensional chains [BCPH22].
Without loss of generality let us assume a 2-local interaction; the general
case follows from blocking the sites first. Let
$\omega\in\mathcal{G}_{\Lambda_{\ell},\beta}^{\text{comm}}$ and consider the
perturbation $\omega^{W}$ from Lemma 3.10. The surface interaction term $W$
has support on the sites $-(\ell+1),-\ell,\ell,\ell+1$ and $\|\exp(\pm\beta
W)\|\leq\exp(2\beta\|h\|)$ is bounded by a constant. We first apply the decay
of correlations to bound the effect of the perturbation. From Lemma 3.10 we
have

$\omega(O)=\frac{\omega^{W}(Oe^{-\beta W})}{\omega^{W}(e^{-\beta W})}.$

Let us introduce the operator
$Y=(id_{\Lambda}\otimes\hat{\omega}_{\partial^{ex}\Lambda})(e^{-\beta
W})\in\mathcal{A}_{\Lambda}$, which satisfies $\|Y\|\leq\exp(2\beta\|h\|)$. We
further decompose $Y=Y_{L}Y_{R}$ where $Y_{L}$, $Y_{R}$ are supported on
$-\ell$ and $\ell$ respectively. We can write

$\omega(O)=\frac{\varphi^{\Lambda}(OY)}{\varphi^{\Lambda}(Y)}$

and thereby

$\displaystyle|\omega(O)-\varphi^{\Lambda}(O)|$
$\displaystyle\leq\frac{|\varphi^{\Lambda}(OY)-\varphi^{\Lambda}(OY_{L})\varphi^{\Lambda}(Y_{R})|}{\varphi^{\Lambda}(Y)}$
$\displaystyle\quad+\frac{|\varphi^{\Lambda}(OY_{L})\varphi^{\Lambda}(Y_{R})-\varphi^{\Lambda}(Y_{L})\varphi^{\Lambda}(Y_{R})\varphi^{\Lambda}(O)|}{\varphi^{\Lambda}(Y)}$
$\displaystyle\quad+\frac{|\varphi^{\Lambda}(O)\varphi^{\Lambda}(Y_{L})\varphi^{\Lambda}(Y_{R})-\varphi^{\Lambda}(O)\varphi^{\Lambda}(Y)|}{\varphi^{\Lambda}(Y)}$
$\displaystyle\leq\|O\|K^{\prime}\exp(-\alpha\ell)$

where to each term we can apply the uniform clustering from [BCPH22, Theorem
6.2], which extends the decay of correlations already known from Araki’s
seminal paper [Ara69] to the case of Gibbs states on finite intervals. The
constant $K^{\prime}$ here comes from the interaction dependent constant from
the cited result combined with the (constant) norms of $Y^{-1}$, $Y_{L}$,
$Y_{R}$ and the locality of $O$ reducing the distance between the supports of
$O$ and $Y$. The constant $\alpha$ depends on the interaction.

We conclude by [PGPH23, Theorem 4.16 (ii)] (see also [Ara69]), which proves
that expectations in the global Gibbs state are approximated by local ones,
with an error decaying exponentially with the distance to the boundary, i.e.,

$\displaystyle|\varphi^{\Lambda}(O)-\langle
O\rangle_{\beta}|\leq\|O\|K^{\prime\prime}\exp(-\delta\ell).$

Again, $\delta$ and $K^{\prime\prime}$ depend on the interaction and the
latter in addition on the locality of $O$. This proves Theorem 3.9 in the case
$D=1$, for all finite $\beta\geq 0$.

In higher dimension we expect the theorem to break down for arbitrary
temperatures due to the existence of phase transitions. Nevertheless, at
sufficiently high temperatures, the two main ingredients, decay of
correlations and local indistinguishability still hold. A technical problem we
are facing is however the (super)exponential growth of $e^{\beta W}$ since the
support and norm of $W$ scale with $\ell^{D-1}$, where $\ell$ is the distance
between $O$ and the boundary. In 2D, however, we can still trade off this
exponential growth against the decay rates by choosing sufficiently high
temperatures. We closely follow the proof of the 1D case. Again we define
$Y=(id_{\Lambda}\otimes\hat{\omega}_{\partial^{ex}\Lambda})(e^{-\beta
W})\in\mathcal{A}_{\Lambda}$. This time, however, we can only bound
$\|Y\|,\|Y^{-1}\|\leq\exp(\beta\ell r\|h\|)$, where $r$ is a constant
depending on the locality of the Hamiltonian. For the decay of correlations we
can resort to the high temperature result in [KGK+14, Theorem 2]. Below a
critical inverse temperature this gives

$\displaystyle|\omega(O)-\varphi^{\Lambda}(O)|$
$\displaystyle=\frac{|\varphi^{\Lambda}(OY)-\varphi^{\Lambda}(O)\varphi^{\Lambda}(Y)|}{\varphi^{\Lambda}(Y)}$
$\displaystyle\leq\|O\|\|Y\|\|Y^{-1}\|C\exp(-\gamma\ell)$
$\displaystyle\leq\|O\|C\exp((2\beta r\|h\|-\gamma)\ell)$

where $C$ and $\gamma$ depend on the interaction and temperature. Note that
$\gamma$ grows unboundedly for $\beta\to 0$. By choosing $\beta$ sufficiently
small such that $2\beta r\|h\|<\gamma$ this implies exponential decay.

Finally, local indistinguishability is proven in [KGK+14, Corollary 2] again
above the same critical temperature and with the same decay rate. The
formulation proves convergence of marginals in 1-norm distance but by choosing
the region to be the support of $O$ the result for observables follows:

$|\varphi^{\Lambda}(O)-\langle
O\rangle_{\beta}|\leq\|O\|K^{\prime\prime}\exp(-\gamma\ell)$

∎

######  Remark 3.11.

We remark that in the proof above, we do not make explicit use of the
uniqueness of the $\beta$-KMS state for $0\leq\beta<\beta_{1}$. This
uniqueness actually follows from the proof since it holds for any choice of
KMS-state and any observable. The used results on decay of correlations based
on cluster expansions are in fact the same techniques used to prove the
uniqueness of the KMS-state at sufficiently high temperature.

##  4 Numerical experiments

We consider the one-dimensional transverse field Ising model

$H=-\sum_{i\in\mathbb{Z}}Z_{i}Z_{i+1}+gX_{i}$ (24)

where $X_{i}$ and $Z_{i}$ are the usual Pauli matrices at site $i$ and $g\geq
0$. The model is exactly solvable [Pfe70], and the ground state exhibits a
phase transition at $g=1$. For $g<1$, the ground space is degenerate, and for
$g>1$, the ground state is unique. The model is gapped for $g\neq 1$, and
gapless at $g=1$. There is no phase transition with respect to the temperature
as this is a one-dimensional model.

Let $\mathcal{G}_{g,\beta}$ be the set of translation-invariant equilibrium
states for the parameter $g$ of the Hamiltonian, and at inverse temperature
$\beta$. For $g<1$, $|\mathcal{G}_{g,\infty}|>1$, while for $g>1$ or
$\beta<\infty$ we have $|\mathcal{G}_{g,\beta}|=1$.

We focus on two observables, namely the $Z$-magnetization per site, i.e.,
$O=Z_{0}$, and the nearest-neighbour correlation function $O=Z_{0}Z_{1}$.

#### Magnetization

When $g=0$, the model is a simple one-dimensional Ising model, and there are
two translation-invariant ground states, namely the all-up and all-down
states. This means that at $g=0$ the magnetization can take any value between
$-1$ and $+1$. For general $g\geq 0$, it has been shown that [Pfe70, Eq.
(3.12a)]

$\\{\omega(Z_{0}):\omega\in\mathcal{G}_{g,\infty}\\}=[-M_{z}(g),M_{z}(g)]$

where

$M_{z}(g)=\begin{cases}(1-g^{2})^{1/8}&\text{if }g<1\\\ 0&\text{if }g\geq
1.\end{cases}$

Figure 1 (left) shows the true values of the magnetization as a function of
$g$, along with the bounds obtained from the convex relaxation (15). The
relaxations were solved respectively with a 4-site and 5-site region
$\Lambda$. (Note that the semidefinite program for a $L$-site region involves
a positivity constraint of a matrix of size $4^{L-1}$.) We see as expected
that the bounds with $L=5$ are tighter than the bounds for $L=4$. We also
observe that the bounds seem very accurate away from the critical point $g=1$.

In Figure 1 (right) we plot the lower and upper bounds on the magnetization at
positive temperature $\beta=1$. Note that for positive temperature, there is a
unique thermal equilibrium state and the average magnetization is 0 for all
$g\geq 0$.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$-1$$-0.5$$0$$0.5$$1$$g$Magnetization
$\langle Z_{0}\rangle$ at $\beta=\infty$$L=4$$L=5$

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$-1$$-0.5$$0$$0.5$$1$$g$Magnetization
$\langle Z_{0}\rangle$ at $\beta=1$$L=4$$L=5$

Figure 1: Upper and lower bounds on magnetization $\langle Z_{0}\rangle$ for
thermal equilibrium states ($\beta=\infty$ and $\beta=1$) of the quantum Ising
model (24), as a function of $g$. The shaded region is the interval
corresponding to $L=5$.

#### Spin-spin correlation function

Next, we consider the nearest-neighbour correlation function $O=Z_{0}Z_{1}$.
The exact value for all choices of $g$ and $\beta\in[0,+\infty]$ is given by
[Sac11, Eq. (10.58)]:999Here the interval reduces to a single point, even when
the set of equilibrium states is degenerate for $g<1$ and $\beta=\infty$.

$\langle Z_{0}Z_{1}\rangle_{g,\beta}=\frac{1}{\pi}\int_{0}^{\pi}\frac{1+g\cos
k}{\sqrt{1+2g\cos k+g^{2}}}\tanh(J\beta\sqrt{1+2g\cos k+g^{2}})dk.$

Figure 2 compares the bounds obtained from the convex relaxations with the
true value. For this observable we see a bigger gap with the true value,
however this gap seems to decay quickly as $g\gg 1$ away from the critical
region. We also notice that at $\beta=\infty$ (ground state), the upper bounds
are much closer to the true value than the lower bounds, especially for $g<1$.

$0.5$$1$$1.5$$0$$0.5$$1$$g$Correlation $\langle Z_{0}Z_{1}\rangle$ at
$\beta=\infty$$L=4$$L=5$

$0.5$$1$$1.5$$0.2$$0.4$$0.6$$0.8$$1$$g$Correlation $\langle Z_{0}Z_{1}\rangle$
at $\beta=1$$L=4$$L=5$

Figure 2: Upper and lower bounds on correlation $\langle Z_{0}Z_{1}\rangle$
for thermal equilibrium states ($\beta=\infty$ and $\beta=1$) of the quantum
Ising model (24), as a function of $g$. The shaded region is the interval
corresponding to $L=5$.

For $\beta<\infty$ the optimization problems involving the matrix relative
entropy were solved using the semidefinite programming approximations [FSP19].
We used the Mosek interior-point solver [ApS22] for all the computations,
except for the case $L=5,\beta<\infty$ where we used the first-order solver
SCS [OCPB16].

## Acknowledgments

We thank Daniel Stilck França for helpful discussions. HF acknowledges funding
from UK Research and Innovation (UKRI) under the UK government’s Horizon
Europe funding guarantee EP/X032051/1. OF acknowledges funding by the European
Research Council (ERC Grant AlgoQIP, Agreement No. 851716) as well as by the
European Union’s Horizon 2020 within the QuantERA II Programme under Grant
VERIqTAS Agreement No 101017733. SOS acknowledges support from the UK
Engineering and Physical Sciences Research Council (EPSRC) under grant number
EP/W524141/1.

##  Appendix A Proof of Theorem 3.4

We adopt the same notations and approach as in the proof of [BR97, Theorem
5.3.15]. In [BR97, page 88] it is shown that for any
$a\in\mathcal{A}_{\text{loc}}$ there are positive Radon measures
$\mu_{a},\nu_{a}$ such that

$\displaystyle\mu_{a}(1)$ $\displaystyle=\omega(a^{*}a)$ (25)
$\displaystyle\mu_{a}(-p)$
$\displaystyle=-\int_{-\infty}^{\infty}pd\mu_{a}(p)=\omega(a^{*}[H,a])$
$\displaystyle\nu_{a}(1)$ $\displaystyle=\omega(aa^{*})$ $\displaystyle
d\nu_{a}(p)$ $\displaystyle=e^{p\beta}d\mu_{a}(p).$

The measures $\mu_{a}$ and $\nu_{a}$ are defined by the implicit relations

$\mu_{a}(\hat{f})=(2\pi)^{-1/2}\int_{-\infty}^{\infty}dtf(t)\omega(a^{*}\alpha_{t}(a))\quad\text{
and
}\quad\nu_{a}(\hat{f})=(2\pi)^{-1/2}\int_{-\infty}^{\infty}dtf(t)\omega(\alpha_{t}(a)a^{*})$

where $\alpha_{t}$ is the time evolution operator, and $\hat{f}$ is the
Fourier transform of $f$,

$\hat{f}(p)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}f(z)e^{-ipz}dz,\qquad
f(z)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}\hat{f}(p)e^{ipz}dp.$

Let $\\{a_{1},\ldots,a_{m}\\}$ be a finite subset of
$\mathcal{A}_{\text{loc}}$. By noting that $\mu_{a}$ and $\nu_{a}$ are
quadratic in $a$, we can define the following positive $m\times m$ Hermitian
matrix-valued measures $\mu$ and $\nu$ by

$[\mu(\hat{f})]_{kl}=(2\pi)^{-1/2}\int_{-\infty}^{\infty}dtf(t)\omega(a_{k}^{*}\alpha_{t}(a_{l})),\quad\text{
and
}\quad[\nu(\hat{f})]_{kl}=(2\pi)^{-1/2}\int_{-\infty}^{\infty}dtf(t)\omega(\alpha_{t}(a_{l})a_{k}^{*}),\quad(1\leq
k,l\leq m).$

(In the finite-dimensional case, if $L=\sum_{j}\lambda_{j}P_{j}$ is a spectral
decomposition of $L=[H,\cdot]$, we have
$\mu_{k,l}=\sum_{j}\delta_{-\lambda_{j}}\omega(a_{k}^{*}P_{j}(a_{l}))$ and
$\nu_{k,l}=\sum_{j}\delta_{-\lambda_{j}}\omega(P_{j}(a_{l})a_{k}^{*})$.)

The matrix-valued measure $\mu$ is related to the measures $\mu_{a}$ by the
following identity

$x^{*}\mu(f)x=\sum_{kl}\overline{x_{k}}x_{l}[\mu(f)]_{kl}=\mu_{\sum_{k}x_{k}a_{k}}(f)$

for all $x\in\mathbb{C}^{m}$ and all integrable functions
$f:\mathbb{R}\to\mathbb{R}$.

Using (25) we see that $\mu(1)=A$. Indeed for any $x\in\mathbb{C}^{m}$ we have

$x^{*}\mu(1)x=\mu_{\sum_{k}x_{k}a_{k}}(1)=\omega\left(\left(\sum_{k}x_{k}a_{k}\right)^{*}\left(\sum_{k}x_{k}a_{k}\right)\right)=\sum_{k,l}\overline{x_{k}}x_{l}\omega(a_{k}^{*}a_{l})=x^{*}Ax.$

Since this is true for all $x\in\mathbb{C}^{m}$, and $\mu(1)-A$ is Hermitian,
it follows that $\mu(1)=A$. One shows in a similar way that $\nu(1)=B$ and
$\mu(-p)=-\int pd\mu(p)=C$. Furthermore, the last identity in (25) imply that
$d\nu(p)=e^{p\beta}d\mu(p)$.

We now prove the inequality $D_{op}(A\|B)\preceq\beta C$. Let $\gamma$ be the
positive real-valued measure defined by
$d\gamma(p)=\operatorname{tr}[d\mu(p)]$. Since the trace is faithful, the
measure $\mu$ has a matrix-valued density with respect to $\gamma$, i.e., we
can write $d\mu(p)=\tilde{A}(p)d\gamma(p)$ where
$\tilde{A}:\mathbb{R}\to\mathbf{H}^{m}_{+}$ (see e.g., [FZ07, Section 4]).
Then we can write

$\displaystyle D_{op}(A\|B)$ $\displaystyle=D_{op}\left(\int 1d\mu(p)\|\int
1d\nu(p)\right)$ $\displaystyle=D_{op}\left(\int\tilde{A}(p)d\gamma(p)\|\int
e^{p\beta}\tilde{A}(p)d\gamma(p)\right)$ $\displaystyle\preceq\int
D_{op}(\tilde{A}(p)\|e^{p\beta}\tilde{A}(p))d\gamma(p)$
$\displaystyle=\beta\int-p\tilde{A}(p)d\mu(p)=\beta\int-pd\mu(p)=C.$

## References

  * [AI22] Dorit Aharonov and Sandy Irani.  Hamiltonian complexity in the thermodynamic limit.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 750–763, 2022.
  * [AKV+23] Mateus Araújo, Igor Klep, Tamás Vértesi, Andrew J. P. Garner, and Miguel Navascués.  Karush-Kuhn-Tucker conditions for non-commutative optimization problems.  preprint, 2023.
  * [Amb14] Andris Ambainis.  On physical problems that are slightly more difficult than QMA.  In 2014 IEEE 29th Conference on Computational Complexity (CCC), pages 32–43. IEEE, 2014.
  * [Ana98] Hirokazu Anai.  On solving semidefinite programming by quantifier elimination.  In Proceedings of the 1998 American Control Conference. ACC (IEEE Cat. No.98CH36207), volume 5, pages 2814–2818 vol.5, 1998.
  * [ApS22] MOSEK ApS.  The MOSEK optimization toolbox for MATLAB manual. Version 10.0., 2022.
  * [Ara69] Huzihiro Araki.  Gibbs states of a one dimensional quantum lattice.  Communications in Mathematical Physics, 14(2):120–157, 1969.
  * [AS77] Huzihiro Araki and Geoffrey L Sewell.  KMS conditions and local thermodynamical stability of quantum lattice systems.  Communications in Mathematical Physics, 52:103–109, 1977.
  * [BC17] Jacob C. Bridgeman and Christopher T. Chubb.  Hand-waving and interpretive dance: an introductory course on tensor networks.  Journal of Physics A: Mathematical and Theoretical, 50(22):223001, 2017.
  * [BCGW22] Sergey Bravyi, Anirban Chowdhury, David Gosset, and Pawel Wocjan.  Quantum hamiltonian complexity in thermal equilibrium.  Nature Physics, 18(11):1367–1370, 2022.
  * [BCPH22] Andreas Bluhm, Ángela Capel, and Antonio Pérez-Hernández.  Exponential decay of mutual information for gibbs states of local hamiltonians.  Quantum, 6:650, 2022.
  * [BCW21] Johannes Bausch, Toby S. Cubitt, and James D. Watson.  Uncomputability of phase diagrams.  Nature Communications, 12(1):452, 2021.
  * [BH12] Thomas Barthel and Robert Hübener.  Solving condensed-matter ground-state problems by semidefinite relaxations.  Phys. Rev. Lett., 108:200404, 2012.
  * [BLMT23] Ainesh Bakshi, Allen Liu, Ankur Moitra, and Ewin Tang.  Learning quantum hamiltonians at any temperature in polynomial time.  arXiv preprint arXiv:2310.02243, 2023.
  * [BP12] Tillmann Baumgratz and Martin B. Plenio.  Lower bounds for ground states of condensed matter systems.  New Journal of Physics, 14(2):023027, 2012.
  * [BPR06] Saugata Basu, Richard Pollack, and Marie-Françoise Roy.  Algorithms in Real Algebraic Geometry.  Springer Berlin Heidelberg, 2006.
  * [BR87] Ola Bratteli and Derek W. Robinson.  Operator Algebras and Quantum Statistical Mechanics 1.  Springer, 2nd edition, 1987.
  * [BR97] Ola Bratteli and Derek W. Robinson.  Operator Algebras and Quantum Statistical Mechanics 2.  Springer, 2nd edition, 1997.
  * [CGL+22] Minjae Cho, Barak Gabai, Ying-Hsuan Lin, Victor A Rodriguez, Joshua Sandor, and Xi Yin.  Bootstrapping the Ising Model on the Lattice.  arXiv preprint arXiv:2206.12538, 2022.
  * [CPGSV21] J. Ignacio Cirac, David Pérez-García, Norbert Schuch, and Frank Verstraete.  Matrix product states and projected entangled pair states: Concepts, symmetries, theorems.  Rev. Mod. Phys., 93:045003, Dec 2021.
  * [CPGW15] Toby S. Cubitt, David Pérez-García, and Michael M. Wolf.  Undecidability of the spectral gap.  Nature, 528(7581):207–211, 2015.
  * [CS23] Minjae Cho and Xin Sun.  Bootstrap, Markov Chain Monte Carlo, and LP/SDP Hierarchy for the Lattice Ising Model.  arXiv preprint arXiv:2309.01016, 2023.
  * [Eff09] Edward G. Effros.  A matrix convexity approach to some celebrated quantum inequalities.  Proc. Natl. Acad. Sci. USA, 106(4):1006–1008, 2009.
  * [ENG11] Ali Ebadian, Ismail Nikoufar, and Madjid Eshaghi Gordji.  Perspectives of matrix convex functions.  Proc. Natl. Acad. Sci. USA, 108(18):7313–7314, 2011.
  * [FFS90] Jun I. Fujii, Masatoshi Fujii, and Yuki Seo.  An extension of the Kubo-Ando theory: solidarities.  Math. Japon, 35(2):387–396, 1990.
  * [FFS23a] Hamza Fawzi, Omar Fawzi, and Samuel O Scalet.  Entropy constraints for ground energy optimization.  arXiv preprint arXiv:2305.06855, 2023.
  * [FFS23b] Hamza Fawzi, Omar Fawzi, and Samuel O. Scalet.  A subpolynomial-time algorithm for the free energy of one-dimensional quantum systems in the thermodynamic limit.  Quantum, 7:1011, 2023.
  * [FK89] Jun I. Fujii and Eizaburo Kamei.  Relative operator entropy in noncommutative information theory.  Math. Japon, 34:341–348, 1989.
  * [FS23] Hamza Fawzi and James Saunderson.  Optimal self-concordant barriers for quantum relative entropies.  SIAM Journal on Optimization, 33(4):2858–2884, 2023.
  * [FSP19] Hamza Fawzi, James Saunderson, and Pablo A Parrilo.  Semidefinite approximations of the matrix logarithm.  Foundations of Computational Mathematics, 19:259–296, 2019.
  * [FV17] Sacha Friedli and Yvan Velenik.  Statistical mechanics of lattice systems: a concrete mathematical introduction.  Cambridge University Press, 2017.
  * [FZ07] Douglas R. Farenick and Fei Zhou.  Jensen’s inequality relative to matrix-valued measures.  Journal of mathematical analysis and applications, 327(2):919–929, 2007.
  * [GHL+15] Sevag Gharibian, Yichen Huang, Zeph Landau, Seung Woo Shin, et al.  Quantum hamiltonian complexity.  Foundations and Trends® in Theoretical Computer Science, 10(3):159–282, 2015.
  * [GI09] Daniel Gottesman and Sandy Irani.  The quantum and classical complexity of translationally invariant tiling and Hamiltonian problems.  In 2009 50th Annual IEEE Symposium on Foundations of Computer Science, pages 95–104. IEEE, 2009.
  * [Han20] Xizhi Han.  Quantum many-body bootstrap.  arXiv preprint arXiv:2006.06002, 2020.
  * [Has06] Matthew B. Hastings.  Solving gapped hamiltonians locally.  Physical Review B, 73(8), feb 2006.
  * [HMS20] Aram W. Harrow, Saeed Mehraban, and Mehdi Soleimanifar.  Classical algorithms, correlation decay, and complex zeros of partition functions of quantum many-body systems.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing. ACM, jun 2020.
  * [KAA21] Tomotaka Kuwahara, Álvaro M. Alhambra, and Anurag Anshu.  Improved thermal area law and quasilinear time algorithm for quantum gibbs states.  Physical Review X, 11(1), mar 2021.
  * [KGK+14] Martin Kliesch, Christian Gogolin, Micahel J. Kastoryano, Arnau Riera, and Jens Eisert.  Locality of temperature.  Physical Review X, 4(3), jul 2014.
  * [KKB20] Tomotaka Kuwahara, Kohtaro Kato, and Fernando GSL Brandão.  Clustering of conditional mutual information for quantum gibbs states above a threshold temperature.  Physical review letters, 124(22):220601, 2020.
  * [KKR06] Julia Kempe, Alexei Kitaev, and Oded Regev.  The complexity of the local hamiltonian problem.  SIAM Journal on Computing, 35(5):1070–1097, 2006.
  * [KSDN22] Ilya Kull, Norbert Schuch, Ben Dive, and Miguel Navascués.  Lower bounding ground-state energies of local hamiltonians through the renormalization group.  arXiv preprint arXiv:2212.03014, 2022.
  * [LVV15] Zeph Landau, Umesh Vazirani, and Thomas Vidick.  A polynomial time algorithm for the ground state of one-dimensional gapped local hamiltonians.  Nature Physics, 11(7):566–569, 2015.
  * [Maz06] David A. Mazziotti.  Quantum chemistry without wave functions: Two-electron reduced density matrices.  Accounts of Chemical Research, 39(3):207–215, 2006.  PMID: 16548509.
  * [MH21] Ryan L. Mann and Tyler Helmuth.  Efficient algorithms for approximating quantum partition functions.  Journal of Mathematical Physics, 62(2), feb 2021.
  * [MSVC15] Andras Molnar, Norbert Schuch, Frank Verstraete, and J Ignacio Cirac.  Approximating gibbs states of local hamiltonians efficiently with projected entangled pair states.  Physical Review B, 91(4):045138, 2015.
  * [Nac06] Bruno Nachtergaele.  Quantum spin systems.  In Encyclopedia of Mathematical Physics, pages 295–301. Academic Press, Oxford, 2006.
  * [OCPB16] Brendan O’donoghue, Eric Chu, Neal Parikh, and Stephen Boyd.  Conic optimization via operator splitting and homogeneous self-dual embedding.  Journal of Optimization Theory and Applications, 169:1042–1068, 2016.
  * [Osb12] Tobias J. Osborne.  Hamiltonian complexity.  Reports on progress in physics, 75(2):022001, 2012.
  * [Pfe70] Pierre Pfeuty.  The one-dimensional Ising model with a transverse field.  Annals of Physics, 57(1):79–90, 1970.
  * [PGPH23] David Pérez-García and Antonio Pérez-Hernández.  Locality estimes for complex time evolution in 1d.  Communications in Mathematical Physics, 399:929–970, 2023.
  * [PH11] David Poulin and Matthew B. Hastings.  Markov entropy decomposition: A variational dual for quantum belief propagation.  Physical Review Letters, 106(8), feb 2011.
  * [PNA10] Stefano Pironio, Miguel Navascués, and Antonio Acín.  Convergent relaxations of polynomial optimization problems with noncommuting variables.  SIAM Journal on Optimization, 20(5):2157–2180, jan 2010.
  * [PRV19] David Poland, Slava Rychkov, and Alessandro Vichi.  The conformal bootstrap: theory, numerical techniques, and applications.  Reviews of Modern Physics, 91(1):015002, 2019.
  * [Sac11] Subir Sachdev.  Quantum Phase Transitions.  Cambridge University Press, 2nd edition, 2011.
  * [WC22] James D. Watson and Toby S. Cubitt.  Computational complexity of the ground state energy density problem.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 764–775, 2022.
  * [Whi92] Steven R. White.  Density matrix formulation for quantum renormalization groups.  Phys. Rev. Lett., 69:2863–2866, Nov 1992.
  * [WSF+23] Jie Wang, Jacopo Surace, Irénée Frérot, Benoît Legat, Marc-Olivier Renou, Victor Magron, and Antonio Acín.  Certifying ground-state properties of many-body systems.  arXiv preprint arXiv:2310.05844, 2023.
  * [ZWJ22] Ruizhe Zhang, Guoming Wang, and Peter Johnson.  Computing ground state properties with early fault-tolerant quantum computers.  Quantum, 6:761, 2022.