# Thermal Area Law for Lattice Bosons

Marius Lemm<EMAIL_ADDRESS>Oliver Siebert  oliver.siebert@uni-
tuebingen.de Department of Mathematics, University of Tübingen, Auf der
Morgenstelle 10, 72076 Tübingen, Germany

###### Abstract

A physical system is said to satisfy a thermal area law if the mutual
information between two adjacent regions in the Gibbs state is controlled by
the area of their boundary. Thermal area laws have been derived for systems
with bounded local interactions such as quantum spin systems. However, for
lattice bosons these arguments break down because the interactions are
unbounded. We rigorously derive a thermal area law for a class of bosonic
Hamiltonians in any dimension which includes the paradigmatic Bose-Hubbard
model. The main idea to go beyond bounded interactions is to introduce a
quasi-free reference state with artificially decreased chemical potential by
means of a double Peierls-Bogoliubov estimate.

Quantum correlations, particularly entanglement, are central to the
foundations of quantum mechanics and have received increasing attention in
recent years due to the rise of quantum information theory and quantum
computing. They can be quantified and compared by associating suitable entropy
functionals to subsystems. The tools from quantum information theory provide
novel and useful ways to analyze quantum many-body systems. A prime example is
the so-called entanglement entropy of a region $A$, defined as the von Neumann
entropy of the reduced density matrix on $A$ which is obtained by taking the
partial trace over the complement of $A$.

It is expected that in typical quantum many-body systems with translation-
invariant short-ranged interactions the entanglement entropy of the ground
state satisfies an area law — meaning that it is bounded by a constant times
the boundary surface area of $A$ (as opposed to the trivial bound which would
entail the volume of $A$). The area law captures our physical intuition that
correlations are concentrated on short distances and therefore only occur
across the boundary cut. It is extremely useful in practice as it severely
restricts the admissible many-body states for approximating ground states
(i.e., quantum matter) and can thus serve to overcome the notorious curse of
dimensionality through the famous density matrix renormalization group (DMRG)
numerical algorithm [1] and its higher-dimensional extensions [2, 3, 4]. The
connection is clearest for 1D lattice systems, where a state satisfies an area
law if and only if it is representable as a matrix product state (MPS) with
fixed bond dimension independent of the system size [5, 6]. For detailed
reviews also covering the higher-dimensional situation, see [7, 8].

Area laws for the entanglement entropy as described above have their origins
in the holographic principle in the context of quantum gravity [9] and have
been numerically observed in a large number of many-body systems. They have
been derived for gapped 1D spin systems [10, 11, 12, 13, 14], for 1D quantum
states with finite correlation lengths [15, 16], gapped harmonic lattice
systems [17, 18, 19, 20], ground states in the same gapped phase as others
obeying an area law [21, 22], models whose Hamiltonian spectra satisfy related
conditions [23, 24], certain frustration-free spin systems [25], tree-graph
systems [26], models exhibiting local topological order [27] and high-
dimensional systems under additional assumptions such as frustration-freeness
[24, 28, 23, 29, 30, 31]. There has also been recent progress for certain
long-range interactions [32]. A general statement in higher dimensions remains
elusive and this problem is known as the area law conjecture.

The analog of the area law for Gibbs states is called the thermal area law.
Gibbs states are important for fundamental reasons, since realistic systems
are never at zero temperature, and for exciting novel applications such as
quantum machine learning [33, 34] and semidefinite programming solvers [35,
36]. The efficient simulability of Gibbs states via tensor networks hinges on
thermal area laws [37, 38, 39]. For Gibbs states, the entanglement between two
regions $A$ and $B$ is quantified by their mutual information

$I(A:B)=S(\rho_{A})+S(\rho_{B})-S(\rho_{AB}).$

Mutual information is the central tool for quantifying correlations in mixed
quantum states [40, 41, 37] and it has found applications diverse research
areas such as condensed-matter physics and general relativity. At zero
temperature, $I(A:B)$ reduces to the entanglement entropy. An area law at
positive temperature was derived in a seminal work of Wolf et. al [42] who
proved for local and bounded interactions that

$I(A:B)\leq C\beta|\partial A|$

whenever $A$ and $B$ are disjoint and together make up the entire lattice.
Very recently, the $\beta$-scaling result was improved to $\beta^{2/3}$ [37]
which matters for the experimentally relevant regime of low temperatures
($\beta\to\infty$). This dependence is not far from optimal, since Gottesman
and Hastings found a 1D model for which the scaling of the mutual information
is at least $\beta^{1/5}$ for large $\beta$ [43]. Moreover, a generalization
to various Rényi generalizations of the mutual information was given in [38]
and this has important applications to simulating and approximating Gibbs
states [39]. Further results of thermal area laws were established for free
fermions [44], for the entanglement negativity (instead of the mutual
information) [45], as a result of rapid mixing for dissipative quantum lattice
systems [46, 47] (with a logarithmic correction) and numerically for some spin
chains showing a $\log\beta$-dependence [48]. A current account of thermal
area laws and related phenomena is given in [41]. There has also been recent
progress in the experimental verification of area laws for both the
entanglement entropy and the mutual information at positive temperature by
means of ultra-cold atom simulators [49].

The applications of thermal areas laws mirror those of their ground state
counterparts. In particular, a thermal area law can be used for the important
problem of approximating the Gibbs state by matrix product operators
(MPO)—mixed-state cousins of the standard matrix product states—and their
higher-dimensional analogs [41, 37, 39, 50]. Just like their pure-state
counterpart MPO approximations allow for a compact representation of the state
and efficient numerical calculation. It is generally believed 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
[41]. Another direct way to use a thermal area law is to note that controlling
the mutual information automatically controls all standard correlation
functions. More precisely, the mutual information upper bounds all truncated
correlations $\operatorname{tr}(M_{A}\otimes
M_{B}\rho_{AB})-\operatorname{tr}(M_{A}\rho_{A})\operatorname{tr}(M_{B}\rho_{B})$
between arbitrary bounded observables $M_{A}$ and $M_{B}$; see [42] and
Corollary 3 below.

A critical limitation of the existing results is that they only hold for
bounded interactions bounded local Hilbert space dimension. This is naturally
the case for quantum spin systems and lattice fermions. However, for lattice
bosons as described, e.g., by the paradigmatic Bose-Hubbard Hamiltonian, the
interactions are unbounded and the standard arguments fail.

There has recently been a surge of interest in bosonic lattice systems and
related models for two main reasons: (i) The Hamiltonians can be
experimentally fine-tuned for cold atoms in optical lattices [51], which makes
them promising platform for quantum simulation and quantum engineering, see
also [52], and (ii) in many cases the standard techniques of quantum
information theory fail (including the derivations of the thermal area law)
because of the unbounded interactions.

Area laws for non-interacting bosons were considered in [19, 20, 53], but the
case of interacting bosons, including the paradigmatic Bose-Hubbard model
proved elusive to rigorous analysis. A partly numerical investigation was
given in [54] with a focus on the phase transition from Mott insulator to
superfluid by analogy with other symmetry breaking transitions [55, 56, 57].
Recently, Abrahamsen et al. rigorously proved an area law for gapped ground
states of 1D bosonic lattice Hamiltonians in [58] by using a truncation of the
local Hilbert spaces and a quantum number tail bound from [59]. Their work
only concerns the zero temperature case and leaves open the positive
temperature case.

Given all this background, it is a fundamental and also rather timely question
whether a thermal area law holds for interacting bosonic lattice gases. In
this work, we answer this question in the affirmative for a broad class of
interacting bosonic lattice gases.

That is, we rigorously derive the first thermal area law for a broad class of
bosonic Hamiltonians in any dimension including the paradigmatic Bose-Hubbard
model. In a nutshell, our result states that under natural assumptions on the
bosonic lattice gases (e.g., short-ranged hopping), we again have the bound

$I(A:B)\leq C\beta|\partial A|$

for all $\beta\geq 1$. The precise result is Theorem 2 below. The
$\beta$-scaling is the same as that found by [42] and our proof uses the same
basic idea as a starting point, namely to use the Gibbs variational principle
to bound the mutual information by a difference of boundary energy
expectations (Lemma 1). However, the unbounded interactions then pose
technical difficulties which we overcome by introducing a quasi-free reference
state with artificially decreased chemical potential by means of a double
Peierls-Bogoliubov estimate. This reduces us to computations with quasi-free
states that can be completed with Wick’s rule. Further details are explained
below and in the Supplemental Material [60].

An immediate application of the result is a general bound on correlation
functions of Gibbs states, cf. Corollary 3. Further physical applications,
e.g., to the approximability of bosonic Gibbs states in the spirit of [41, 37]
are postponed to future work.

##  I Setup for infinite-dimensional Hilbert spaces

A technical point in the description of many-boson systems is that the local
Hilbert spaces are infinite-dimensional since the particle number is
unbounded. For this reason, we include this short preliminary section in which
we recall the elegant approach to thermal area laws via the Gibbs variational
principle by Wolf et al. [42] and note that it adapts straightforwardly to the
infinite-dimensional situation. These abstract results are then utilized for
the Bose-Hubbard model in the following section.

Let $\mathfrak{h}$ be the local Hilbert space for one site, $\Lambda$ a finite
set and $A,B\subseteq\Lambda$ such that $A\cap B=\emptyset$ and $A\sqcup
B=\Lambda$. We set

$\displaystyle\mathcal{H}_{A}=\bigotimes_{x\in
A}\mathfrak{h},\quad\mathcal{H}_{B}=\bigotimes_{x\in
B}\mathfrak{h},\quad\mathcal{H}_{AB}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}.$

We suppose that the Hamiltonian can be decomposed as (cf. Figure 1)

$\displaystyle H_{AB}=H_{A}\otimes\mathds{1}+\mathds{1}\otimes
H_{B}+H_{\partial}$ (1)

Figure 1: Periodic box for $d=2$ and $L=9$ decomposed into two regions $A$ and
$B$. The bonds connecting $A$ and $B$ make up the boundary Hamiltonian
$H_{\partial}$

For details about operator domains, which are relevant because our operators
are unbounded, see [60]. Let $\mathcal{T}_{1}^{+}(\mathcal{H})$ denote the set
of all density matrices on $\mathcal{H}$. Then we define the free energy as

$F^{\beta}_{AB}=\inf_{\rho\in\mathcal{T}_{1}^{+}(\mathcal{H}_{AB})}\mathcal{F}^{\beta}(\rho),$

where

$\displaystyle\mathcal{F}^{\beta}(\rho)=\operatorname{tr}(H_{AB}\rho)-\frac{S(\rho)}{\beta},\quad
S(\rho)=-\operatorname{tr}(\rho\ln\rho).$

The Gibbs state

$\rho_{AB}=\frac{e^{-\beta H_{AB}}}{Z_{AB}},\quad
Z_{AB}=\operatorname{tr}e^{-\beta H_{AB}},$

minimizes the free energy, i.e.,
$F^{\beta}_{AB}=\mathcal{F}^{\beta}(\rho_{AB})$. (See Proposition 5 for a
proof of this in the infinite-dimensional setting.) We reduce to the $A$,
respectively $B$ subsystem by taking partial traces of the Gibbs state

$\displaystyle\rho_{A}=\operatorname{tr}_{\mathcal{H}_{B}}(\rho_{AB}),\qquad\rho_{B}=\operatorname{tr}_{\mathcal{H}_{A}}(\rho_{AB}).$

The general bound from [42] straightforwardly extends to the infinite-
dimensional setting as follows.

######  Lemma 1 (Boundary energy controls mutual informat.)

Let $\beta>0$ and suppose that $\operatorname{tr}e^{-\beta H_{AB}}<\infty$.
Then

$I(A:B)\leq\beta\operatorname{tr}(H_{\partial}(\rho_{A}\otimes\rho_{B}-\rho_{AB})).$

###### Proof.

We start from

$\mathcal{F}^{\beta}(\rho_{AB})=F^{\beta}_{AB}\leq\mathcal{F}^{\beta}(\rho_{A}\otimes\rho_{B}).$

Using $S(\rho_{A}\otimes\rho_{B})=S(\rho_{A})+S(\rho_{B})$, we obtain

$I(A:B)\leq\beta\operatorname{tr}(H_{AB}(\rho_{A}\otimes\rho_{B}-\rho_{AB}))$

and the right-hand side equals
$\beta\operatorname{tr}(H_{\partial}(\rho_{A}\otimes\rho_{B}-\rho_{AB}))$ by
basic properties of the partial trace. ∎

##  II Main Result

Now we consider the Bose-Hubbard model in the framework of the previous
section. Let $\Lambda\equiv\Lambda_{L}$ denote a box of side length $L$ in the
$d$-dimensional lattice with periodic boundary conditions. For
$x,y\in\Lambda_{L}$ we write $x\sim y$ if $x$ and $y$ are nearest neighbors in
the periodized lattice. At each site lives a bosonic particle described by the
local Hilbert space $\mathfrak{h}=\ell^{2}({\mathbb{N}})$.

The total Hilbert space
$\mathcal{H}_{AB}=\otimes_{x\in\Lambda_{L}}\mathfrak{h}$ is isomorphic to the
Fock space $\mathfrak{F}(\ell^{2}(\Lambda_{L}))$. On it, we consider the Bose-
Hubbard Hamiltonian

$\displaystyle H_{AB}=-J\sum_{x\sim
y}a_{x}^{\dagger}a_{y}+\frac{U}{2}\sum_{x\in\Lambda_{L}}n_{x}(n_{x}-1)-\mu\mathcal{N},$
(2)

where $J\in\mathbb{R}$ represents the strength of the kinetic nearest-neighbor
hopping, $U>0$ the strength of the on-site repulsion and $\mu\in\mathbb{R}$
the chemical potential, and $\mathcal{N}=\sum_{x\in\Lambda_{L}}n_{x}$ is the
total number operator. The Hamiltonian is self-adjoint on a suitable domain
$\mathcal{D}(H_{AB})$; see e.g., [61, 60].

We are now ready to state the main result. We decompose the box into two
regions $A$ and $B$ as shown in Figure 1. We write $L_{A}\in\\{1,\ldots,L\\}$
for the width of region $A$.

######  Theorem 2 (Main result: thermal area law)

For all $\beta,U,\mu>0$, we have

$\displaystyle I(A:B)\leq c(J,U,\mu)\max\\{1,\beta\\}L^{d-1}.$

Since $L^{d-1}$ is the size of the boundary between $A$ and $B$, this proves
the thermal area law.

A few remarks are in order: (i) The repulsiveness assumption $U>0$ is
necessary as it ensures stability of the system. (ii) The assumption that
$\mu>0$ is standard, cf. the usual phase diagram in Figure 13(a) in [51].
Indeed, note that if one would take $\mu$ sufficiently negative, then the
system becomes effectively devoid of particles in the grand-canonical setting.
(iii) The constant $c(J,U,\mu)$ can be made explicit from (S5). (iv) The
maximum $\max\\{1,\beta\\}$ means that the bound behaves as $\beta$ for low
temperatures in accordance with [42]. We recall that some growth in $\beta$ is
strongly expected without any gap assumption [43]. For high temperature
$\beta<1$, we find a temperature independent lower bound which matches the
classical situation [42].

As explained in [42], a bound on the mutual information implies a bound for
the correlation of any pair of bounded observables and the same is true in the
bosonic setting. Let $M_{A},M_{B}$ be two bounded self-adjoint operators on
$\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively. We denote their
truncated correlation function as

$\mathcal{C}(A,B)=\operatorname{tr}(M_{A}\otimes
M_{B}\rho_{AB})-\operatorname{tr}(M_{A}\rho_{A})\operatorname{tr}(M_{B}\rho_{B}).$

######  Corollary 3 (Control on truncated correlations)

For all $\beta,U,\mu>0$, we have

$\displaystyle\frac{\mathcal{C}(A,B)^{2}}{2\left\|M_{A}\right\|^{2}\left\|M_{B}\right\|^{2}}\leq
2c(J,U,\mu)\max\\{1,\beta\\}L^{d-1},$

with the same constant as in Theorem 2.

###### Proof of Corollary 3.

This follows from the quantum Pinsker inequality
$I(A:B)\geq\frac{1}{2}\left\|\rho_{AB}-\rho_{A}\otimes\rho_{B}\right\|_{1}^{2}$
and $\left\|X\right\|_{1}\geq\operatorname{tr}(XY)/\left\|Y\right\|$. We
mention that the standard proof of the quantum Pinkser inequality, cf. [62,
Theorem 1.15] extends to infinite dimensions since the data processing
inequality holds for trace-class operators [63]. ∎

We close the presentation by discussing several extensions of the result which
can be obtained from the same methods. First, the proof extends to on-site
interaction terms $\sum_{x}f(n_{x})$ for any function $f\geq 0$ which is
polynomially bounded and grows at least linearly at infinity.

Second, it is possible to replace the nearest-neighbor hopping term
$-J\sum_{x\sim y}a_{x}^{\dagger}a_{y}$ by a more general hopping. For example,
we can treat higher order hopping if the on-site repulsion is sufficiently
strong. More precisely, a quite general choice of hopping would be

$\displaystyle\sum_{k=1}^{M}\sum_{l=1}^{k}$
$\displaystyle\sum_{\begin{subarray}{c}\\{x_{1},\ldots,x_{l},\\\
y_{1},\ldots,y_{l}\\}\subset\Lambda_{L}\end{subarray}}\sum_{\nu_{1}+\ldots+\nu_{l}=k}J^{\nu_{1}\ldots\nu_{l}}_{x_{1}\ldots
x_{l}y_{1}\ldots y_{l}}$
$\displaystyle\times(a^{\dagger}_{x_{1}})^{\nu_{1}}\ldots(a^{\dagger}_{x_{l}})^{\nu_{l}}a_{y_{1}}^{\nu_{1}}\ldots
a_{y_{l}}^{\nu_{l}},$

where $M\in{\mathbb{N}}$ such that the function $f$ (of the generalized on-
site interaction) is growing faster than $x^{M}$, and
$J^{\nu_{1}\ldots\nu_{l}}_{x_{1}\ldots x_{k}y_{1}\ldots y_{k}}$ is uniformly
bounded and finite-range in the sense that it vanishes if the distance between
two points is greater than some fixed constant. Furthermore, we can also
consider a sufficiently at infinity decaying hopping term, e.g.,
$-\sum_{x,y}J_{xy}a_{x}^{\dagger}a_{y}$ for $J_{x}\geq 0$,
$x\in{\mathbb{Z}}^{d}$, satisfying $J_{xy}\leq(1+|x-y|)^{-\alpha}$ for
sufficiently large $\alpha>d$.

Third, the underlying lattice structure can be easily modified as well,
though, the constant will be less explicit since it depends on the spectrum of
the graph Laplacian. In summary, the thermal area law can be proved for an
entire class of bosonic lattice gases in any dimension.

##  III Sketch of proof of Theorem 2

The detailed proof of Theorem 2 is given in the supplemental material [60].
Here we give a sketch of the main ideas.

We decompose the Hamiltonian as

$H_{AB}=H_{0}+W+(\mu-2dJ)\mathcal{N}$

where $H_{0}=-J\sum_{x\sim y}a_{x}^{\dagger}a_{y}+2dJ\mathcal{N}$ is the
shifted kinetic term and $W=\frac{U}{2}\sum_{x\in X}n_{x}(n_{x}-1)$ is the on-
site interaction.

To use Lemma 1, we again decompose the Hamiltonian as
$H_{AB}=H_{A}+H_{B}+H_{\partial}$ where we define the subsystem Hamiltonians
with open boundary conditions along the cut. More precisely, for
$X\in\\{A,B\\}$, we set

$\displaystyle H_{X}=-J\sum_{\begin{subarray}{c}x\sim y,\\\ x,y\in
X\end{subarray}}a_{x}^{\dagger}a_{y}+\frac{U}{2}\sum_{x\in
X}n_{x}(n_{x}-1)-\mu\sum_{x\in X}n_{x}$

and $H_{\partial}=-J\sum_{\begin{subarray}{c}x\sim y,\\\ x\in A,y\in
B\end{subarray}}(a_{x}^{\dagger}a_{y}+a_{y}^{\dagger}a_{x})$.

Notice that the boundary Hamiltonian $H_{\partial}$ contains $L^{d-1}$ terms,
so the right-hand side of Lemma 1 seems to exhibits the desired scaling in $L$
and $\beta$. The main challenge is that the hopping terms
$a_{x}^{\dagger}a_{y}+a_{y}^{\dagger}a_{x}$ between $A$ and $B$ regions are
unbounded in contrast to the cases of spin systems or lattice fermions.

Our first idea is that since expectations with respect to the full Gibbs state
are rather difficult to handle, we aim for the expectation in a quasifree
state which can be computed via Wick’s rule. To this end, we want to remove
the $W$ term in $e^{-\beta H_{AB}}$. The technical tool to rigorously
implement such a shift in the operator exponent will be a double application
of the Peierls-Bogoliubov inequality [64], which has a long history in the
study of quantum many-body systems and quantum information theory. Applying it
twice, we obtain

$\displaystyle\frac{\operatorname{tr}(Pe^{K})}{\operatorname{tr}(e^{K})}\leq\frac{\operatorname{tr}(Pe^{K+P})}{\operatorname{tr}(e^{K+P})}.$
(3)

We use this bound with $K=-\beta H_{AB}$ and
$P=\beta(W-(\mu+\gamma-2dJ)\mathcal{N})$ so that the new effective Hamiltonian
is

$K+P=-\beta(H_{0}+(\gamma+2dJ)\mathcal{N})$

and $e^{-\beta(H_{0}+(\gamma+2dJ)\mathcal{N})}$ is a trace-class quasifree
state, so that expectations can be calculated via Wick’s rule. Here we
introduced a parameter $\gamma>0$ large enough in order to make
$e^{-\beta(H_{0}+(\gamma+2dJ)\mathcal{N})}$ normalizable. This can be
interpreted as an artificial decrease of the chemical potential which is
introduced to stabilize the system by balancing the loss of the repulsion $W$
from the exponential. (Indeed, note that without help from $\gamma$, we would
get $e^{-\beta(H_{0}-(\mu+2dJ)\mathcal{N})}$ which has infinite trace for
$\mu>-2dJ$.)

The final expression that we arrive at via (3) can then be evaluated using
Wick’s theorem. Subsequently, by means of the corresponding one-particle
density operator, we obtain rather explicit expressions which amount to a
Riemann sum of the density of Planck’s law. This density decays for large
$\beta$ in an integrable way. In particular, we see that the resulting
expression as well as the error terms are bounded for $\beta\geq 1$ and
Theorem 2 follows. For the details, see [60].

##  IV Conclusions

We presented a rigorous proof for a thermal area law for the Bose-Hubbard
model and related bosonic lattice gases. This result closes a gap in the
recently growing literature about the quantum information theory of lattice
bosons and provides a positive-temperature counterpart to the area law for
gapped bosonic ground states (in 1D) [58].

The idea of the proof is based on the general idea in [42] together with a
Peierls-Bogoliubov argument which artificially decreases the chemical
potential in order to get a trace-class free reference state. The method is
highly robust and extends to many other bosonic lattice gases and any lattice
dimension.

Natural follow-up problems include the approximability of bosonic thermal
states by generalized matrix product operators in the spirit of [50]. This is
related to area laws for the generalized Rényi entropy [6], so it would be
interesting to extend the present results to some Rényi generalizations of the
mutual information (as in [38]), see also [43].

On the one hand, our result is of fundamental nature in the quantum
information theory of lattice systems. On the other hand, it paves the way for
future information-theoretic studies of lattice bosons such as the ones
described above. The goal is to unlock the full potential of these
experimentally finely tunable systems for modern applications such as quantum
machine learning [33, 34] and semidefinite programming solvers [35, 36].

## Acknowledgments

The authors are grateful to Álvaro Alhambra and Ángela Capel Cuevas for useful
comments on a draft version of the manuscript.

## References

  * White [1993] S. R. White, Density-matrix algorithms for quantum renormalization groups, Physical Review B 48, 10345 (1993).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions, arXiv preprint cond-mat/0407066 (2004).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group: a short introduction, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369, 2643 (2011).
  * Stoudenmire and White [2012] E. M. Stoudenmire and S. R. White, Studying two-dimensional systems with the density matrix renormalization group, Annu. Rev. Condens. Matter Phys. 3, 111 (2012).
  * Dalzell and Brandão [2019] A. M. Dalzell and F. G. Brandão, Locally accurate MPS approximations for ground states of one-dimensional gapped local Hamiltonians, Quantum 3, 187 (2019).
  * Verstraete and Cirac [2006] F. Verstraete and J. I. Cirac, Matrix product states represent ground states faithfully, Physical Review B 73, 094423 (2006).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Reviews of Modern Physics 82, 277 (2010).
  * Ge and Eisert [2016] Y. Ge and J. Eisert, Area laws and efficient descriptions of quantum many-body states, New Journal of Physics 18, 083026 (2016).
  * Bousso [2002] R. Bousso, The holographic principle, Reviews of Modern Physics 74, 825 (2002).
  * Hastings [2007a] M. B. Hastings, An area law for one-dimensional quantum systems, Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007a).
  * Arad  _et al._ [2012] I. Arad, Z. Landau, and U. Vazirani, Improved one-dimensional area law for frustration-free systems, Physical Review B 85, 195145 (2012).
  * Arad  _et al._ [2013] I. Arad, A. Kitaev, Z. Landau, and U. Vazirani, An area law and sub-exponential algorithm for 1D systems, arXiv preprint arXiv:1301.1162 (2013).
  * Arad  _et al._ [2017] I. Arad, Z. Landau, U. Vazirani, and T. Vidick, Rigorous rg algorithms and area laws for low energy eigenstates in 1D, Communications in Mathematical Physics 356, 65 (2017).
  * Huang [2014] Y. Huang, Area law in one dimension: Degenerate ground states and Renyi entanglement entropy, arXiv preprint arXiv:1403.0327 (2014).
  * Brandão and Horodecki [2013] F. G. Brandão and M. Horodecki, An area law for entanglement from exponential decay of correlations, Nature Physics 9, 721 (2013).
  * Cho [2018] J. Cho, Realistic area-law bound on entanglement from exponentially decaying correlations, Physical Review X 8, 031009 (2018).
  * Audenaert  _et al._ [2002] K. Audenaert, J. Eisert, M. Plenio, and R. Werner, Entanglement properties of the harmonic chain, Physical Review A 66, 042327 (2002).
  * Plenio  _et al._ [2005] M. B. Plenio, J. Eisert, J. Dreissig, and M. Cramer, Entropy, entanglement, and area: analytical results for harmonic lattice systems, Physical Review Letters 94, 060503 (2005).
  * Cramer and Eisert [2006] M. Cramer and J. Eisert, Correlations, spectral gap and entanglement in harmonic quantum systems on generic lattices, New Journal of Physics 8, 71 (2006).
  * Cramer  _et al._ [2006] M. Cramer, J. Eisert, M. B. Plenio, and J. Dreissig, Entanglement-area law for general bosonic harmonic lattice systems, Physical Review A 73, 012309 (2006).
  * Van Acoleyen  _et al._ [2013] K. Van Acoleyen, M. Mariën, and F. Verstraete, Entanglement rates and area laws, Physical Review Letters 111, 170501 (2013).
  * Mariën  _et al._ [2016] M. Mariën, K. M. Audenaert, K. Van Acoleyen, and F. Verstraete, Entanglement rates and the stability of the area law for the entanglement entropy, Communications in Mathematical Physics 346, 35 (2016).
  * Hastings [2007b] M. B. Hastings, Entropy and entanglement in quantum ground states, Physical Review B 76, 035114 (2007b).
  * Masanes [2009] L. Masanes, Area law for the entropy of low-energy states, Physical Review A 80, 052104 (2009).
  * de Beaudrap  _et al._ [2010] N. de Beaudrap, M. Ohliger, T. Osborne, and J. Eisert, Solving frustration-free spin systems, Physical Review Letters 105, 060504 (2010).
  * Abrahamsen [2019] N. Abrahamsen, A polynomial-time algorithm for ground states of spin trees, arXiv preprint arXiv:1907.04862 (2019).
  * Michalakis [2012] S. Michalakis, Stability of the area law for the entropy of entanglement, arXiv preprint arXiv:1206.6900 (2012).
  * Brandao and Cramer [2015] F. G. Brandao and M. Cramer, Entanglement area law from specific heat capacity, Physical Review B 92, 115134 (2015).
  * Cho [2014] J. Cho, Sufficient condition for entanglement area laws in thermodynamically gapped spin systems, Physical Review Letters 113, 197204 (2014).
  * Anshu  _et al._ [2020] A. Anshu, I. Arad, and D. Gosset, Entanglement subvolume law for 2D frustration-free spin systems, in  _Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing_ (2020) pp. 868–874.
  * Anshu  _et al._ [2022] A. Anshu, I. Arad, and D. Gosset, An area law for 2d frustration-free spin systems, in  _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_ (2022) pp. 12–18.
  * Kuwahara and Saito [2020] T. Kuwahara and K. Saito, Area law of noncritical ground states in 1D long-range interacting systems, Nature Communications 11, 1 (2020).
  * Amin  _et al._ [2018] M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, Quantum boltzmann machine, Physical Review X 8, 021050 (2018).
  * Anshu  _et al._ [2021] A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, Sample-efficient learning of interacting quantum systems, Nature Physics 17, 931 (2021).
  * Brandao and Svore [2017] F. G. Brandao and K. M. Svore, Quantum speed-ups for solving semidefinite programs, in  _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ (IEEE, 2017) pp. 415–426.
  * Van Apeldoorn  _et al._ [2020] J. Van Apeldoorn, A. Gilyén, S. Gribling, and R. de Wolf, Quantum sdp-solvers: Better upper and lower bounds, Quantum 4, 230 (2020).
  * Kuwahara  _et al._ [2021] T. Kuwahara, Á. M. Alhambra, and A. Anshu, Improved thermal area law and quasilinear time algorithm for quantum Gibbs states, Physical Review X 11, 011047 (2021).
  * Scalet  _et al._ [2021] S. O. Scalet, Á. M. Alhambra, G. Styliaris, and J. I. Cirac, Computable Rényi mutual information: Area laws and correlations, Quantum 5, 541 (2021).
  * Alhambra and Cirac [2021] A. M. Alhambra and J. I. Cirac, Locally accurate tensor networks for thermal states and time evolution, PRX Quantum 2, 040331 (2021).
  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang, Quantum computation and quantum information (2002).
  * Alhambra [2022] Á. M. Alhambra, Quantum many-body systems in thermal equilibrium, arXiv preprint arXiv:2204.08349 (2022).
  * Wolf  _et al._ [2008] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Area laws in quantum systems: mutual information and correlations, Physical Review Letters 100, 070502 (2008).
  * Gottesman and Hastings [2010] D. Gottesman and M. B. Hastings, Entanglement versus gap for one-dimensional spin systems, New Journal of Physics 12, 025002 (2010).
  * Bernigau  _et al._ [2015] H. Bernigau, M. J. Kastoryano, and J. Eisert, Mutual information area laws for thermal free fermions, Journal of Statistical Mechanics: Theory and Experiment 2015, P02008 (2015).
  * Sherman  _et al._ [2016] N. E. Sherman, T. Devakul, M. B. Hastings, and R. R. Singh, Nonzero-temperature entanglement negativity of quantum spin models: Area law, linked cluster expansions, and sudden death, Physical Review E 93, 022128 (2016).
  * Kastoryano and Eisert [2013] M. J. Kastoryano and J. Eisert, Rapid mixing implies exponential decay of correlations, Journal of Mathematical Physics 54, 102201 (2013).
  * Brandao  _et al._ [2015] F. G. Brandao, T. S. Cubitt, A. Lucia, S. Michalakis, and D. Perez-Garcia, Area law for fixed points of rapidly mixing dissipative quantum systems, Journal of Mathematical Physics 56, 102202 (2015).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. Prosen, and I. Pižorn, Complexity of thermal states in quantum spin chains, Physical Review A 78, 022103 (2008).
  * Tajik  _et al._ [2022] M. Tajik, I. Kukuljan, S. Sotiriadis, B. Rauer, T. Schweigler, F. Cataldini, J. Sabino, F. Møller, P. Schüttelkopf, S.-C. Ji, _et al._ , Experimental verification of the area law of mutual information in quantum field theory, arXiv preprint arXiv:2206.10563 (2022).
  * Guth Jarkovský  _et al._ [2020] J. c. v. Guth Jarkovský, A. Molnár, N. Schuch, and J. I. Cirac, Efficient description of many-body systems with matrix product density operators, PRX Quantum 1, 010304 (2020).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Reviews of Modern Physics 80, 885 (2008).
  * Childs  _et al._ [2014] A. M. Childs, D. Gosset, and Z. Webb, The bose-hubbard model is qma-complete, in  _International Colloquium on Automata, Languages, and Programming_ (Springer, 2014) pp. 308–319.
  * Cramer  _et al._ [2007] M. Cramer, J. Eisert, and M. Plenio, Statistics dependence of the entanglement entropy, Physical Review Letters 98, 220603 (2007).
  * Alba  _et al._ [2013] V. Alba, M. Haque, and A. M. Läuchli, Entanglement spectrum of the two-dimensional bose-hubbard model, Physical Review Letters 110, 260403 (2013).
  * Metlitski and Grover [2011] M. A. Metlitski and T. Grover, Entanglement entropy of systems with spontaneously broken continuous symmetry, arXiv preprint arXiv:1112.5166 (2011).
  * Kallin  _et al._ [2011] A. B. Kallin, M. B. Hastings, R. G. Melko, and R. R. Singh, Anomalies in the entanglement properties of the square-lattice Heisenberg model, Physical Review B 84, 165134 (2011).
  * Song  _et al._ [2011] H. F. Song, N. Laflorencie, S. Rachel, and K. Le Hur, Entanglement entropy of the two-dimensional heisenberg antiferromagnet, Physical Review B 83, 224410 (2011).
  * Abrahamsen  _et al._ [2022] N. Abrahamsen, Y. Su, Y. Tong, and N. Wiebe, Entanglement area law for 1D gauge theories and bosonic systems, arXiv preprint arXiv:2203.16012 (2022).
  * Tong  _et al._ [2021] Y. Tong, V. V. Albert, J. R. McClean, J. Preskill, and Y. Su, Provably accurate simulation of gauge theories and bosonic systems, arXiv preprint arXiv:2110.06942 (2021).
  * [60] M. Lemm and O. Siebert, Supplemental material: Thermal area law for lattice bosons.
  * Faupin  _et al._ [2022] J. Faupin, M. Lemm, and I. M. Sigal, Maximal speed for macroscopic particle transport in the bose-hubbard model, Phys. Rev. Lett. 128, 150602 (2022).
  * Ohya and Petz [2004] M. Ohya and D. Petz,  _Quantum entropy and its use_ (Springer Science & Business Media, 2004).
  * Lindblad [1975] G. Lindblad, Completely positive maps and entropy inequalities, Communications in Mathematical Physics 40, 147 (1975).
  * Ruelle [1999] D. Ruelle,  _Statistical mechanics: Rigorous results_ (World Scientific, 1999).
  * Carlen [2010] E. Carlen, Trace inequalities and quantum entropy: an introductory course, in  _Entropy and the quantum_, Contemp. Math., Vol. 529 (Amer. Math. Soc., Providence, RI, 2010) pp. 73–140.
  * Ruskai [1972] M. B. Ruskai, Inequalities for traces on von neumann algebras, Communications in Mathematical Physics 26, 280 (1972).
  * Ruelle [1969] D. Ruelle,  _Statistical Mechanics: Rigorous Results_ (New York: W.A. Benjamin, 1969).
  * Bratteli and Robinson [2003] O. Bratteli and D. Robinson,  _Operator Algebras and Quantum Statistical Mechanics: Equilibrium States. Models in Quantum Statistical Mechanics_, Theoretical and Mathematical Physics (Springer Berlin Heidelberg, 2003).

Supplemental Material:
Thermal Area Law for Lattice Bosons

##  V Preliminaries

###  V.1 Operator domains

For the general setup, we require the following statements about operator
domains. We assume that $H_{A}$ and $H_{B}$ are densely defined symmetric
operators on $\mathcal{D}(A)$ and $\mathcal{D}(B)$, respectively and the
boundary Hamiltonian is defined on $H_{\partial}$ on
$\mathcal{D}(A)\otimes\mathcal{D}(B)$. We assume that $H_{AB}$ is essentially
self-adjoint on $\mathcal{D}(A)\otimes\mathcal{D}(B)$ which is then the
appropriate domain for the equality

$\displaystyle H_{AB}=H_{A}\otimes\mathds{1}+\mathds{1}\otimes
H_{B}+H_{\partial}$ (S1)

Concerning the Bose-Hubbard model, in the total Fock space
$\mathfrak{F}(\ell^{2}(\Lambda_{L}))$, we consider the dense domain

$\mathfrak{F}_{\operatorname{fin}}(\ell^{2}(\Lambda_{L}))=\\{\psi\in\mathfrak{F}(\ell^{2}(\Lambda_{L})):\exists
n_{0}\in{\mathbb{N}}:\forall n\geq n_{0}:\psi_{n}=0\\}.$

The operator $H_{AB}$ is self-adjoint on the largest domain
$\mathcal{D}(H_{AB})$, where it can be defined, cf. [61]. It then follows from
$[H_{AB},\mathcal{N}]=0$ that $H_{AB}$ is indeed essentially self-adjoint on
$\mathfrak{F}_{\operatorname{fin}}(\ell^{2}(\Lambda_{L}))$. Moreover, for
subsystems $X\in\\{A,B\\}$, we always set
$\mathcal{D}(X)=\mathfrak{F}_{\operatorname{fin}}(\mathcal{H}_{X})$. The above
properties can then be verified.

###  V.2 Trace inequalities in infinite dimensions

######  Lemma 4 (Peierls-Bogoliubov inequality)

Let $\mathcal{N}\geq 0$ be self-adjoint operator with purely discrete
spectrum, let $\Pi_{N}:=\mathds{1}_{\mathcal{N}\leq N}$, $N\in{\mathbb{N}}$
and assume that $\Pi_{N}\mathcal{H}$ is finite-dimensional for all
$N\in{\mathbb{N}}$. Let $(\mathcal{D}(K),K)$ be a self-adjoint and
$(\mathcal{D}(P),P)$ be a symmetric operator such that $[K,\mathcal{N}]=0$ and
$[P,\mathcal{N}]=0$, and assume that $K+P$ is self-adjoint. Then we have

$\displaystyle\frac{\operatorname{tr}(Pe^{K})}{\operatorname{tr}(e^{K})}\leq\log\left(\frac{\operatorname{tr}(e^{K+P})}{\operatorname{tr}(e^{K})}\right),$
(S2)

provided that all traces exist. In particular, we obtain (3), i.e.,

$\displaystyle\frac{\operatorname{tr}(Pe^{K})}{\operatorname{tr}(e^{K})}\leq\frac{\operatorname{tr}(Pe^{K+P})}{\operatorname{tr}(e^{K+P})}.$
(S3)

###### Proof of Lemma 4.

The inequality (S2) is well-known for matrices, see for example [65, (2.14)]
(or more generally if $P$ is bounded [66]). Therefore, we have

$\displaystyle\frac{\operatorname{tr}(Pe^{K}\Pi_{N})}{\operatorname{tr}(e^{K}\Pi_{N})}\leq\log\left(\frac{\operatorname{tr}(e^{K+P}\Pi_{N})}{\operatorname{tr}(e^{K}\Pi_{N})}\right),$

where $\Pi_{N}=\mathds{1}_{\mathcal{N}\leq N}$. Taking the limit $N\to\infty$
yields the desired result. Finally, (3) follows from (S2) by means of

$\log\left(\frac{\operatorname{tr}(e^{K+P})}{\operatorname{tr}(e^{K})}\right)=-\log\left(\frac{\operatorname{tr}(e^{K})}{\operatorname{tr}(e^{K+P})}\right)\leq\frac{\operatorname{tr}(Pe^{K+P})}{\operatorname{tr}(e^{K+P})}.$

This proves Lemma 4. ∎

######  Proposition 5 (Gibbs variational principle)

Let $(H,\mathcal{D}(H))$ be a self-adjoint operator such that
$\operatorname{tr}e^{-\beta H}<\infty$ for all $\beta>0$. Let
$\mathcal{F}^{\beta}(\rho):=\left(\operatorname{tr}(H\rho)-\frac{S(\rho)}{\beta}\right)$.
Then

$\inf_{\rho\in\mathcal{T}_{1}^{+}(\mathcal{H})}\mathcal{F}^{\beta}(\rho)=\mathcal{F}^{\beta}(\rho_{H}),$

where $\rho_{H}=e^{-\beta H}/\operatorname{tr}e^{-\beta H}$.

###### Proof of Proposition 5.

We have

$\operatorname{tr}(K\ln K-K\ln P)\geq\operatorname{tr}(K-P)$

for all positive self-adjoint trace class operators $K,P$ [67]. This yields
[68, Section 5.3.1]

$\displaystyle\mathcal{F}^{\beta}(\rho)$
$\displaystyle=\beta^{-1}\operatorname{tr}(\rho\ln\rho-\rho\ln\rho_{H})-\beta^{-1}\ln\operatorname{tr}e^{-\beta
H}\geq-\beta^{-1}\ln\operatorname{tr}e^{-\beta
H}=\mathcal{F}^{\beta}(\rho_{H}).\qed$

##  VI Proof of Theorem 2

To begin, we notice that we may assume without loss of generality that $J>0$.
Indeed, if $J=0$, the Gibbs state is a product state (Mott insulator) and the
claim is trivial and if $J<0$ we can employ the unitary transformation
$a_{x}\to-a_{x}$ at every second lattice site to reduce to the case $J>0$.

###  VI.1 Step 1: Controlling boundary energy by particle number

In this section, we prove the following bound

###### Proposition 6

We have

$I(A:B)\leq\frac{8}{L}\beta J\operatorname{tr}(\mathcal{N}\rho_{AB}).$

###### Proof of Proposition 6.

Using Lemma 1 and the operator Cauchy-Schwarz inequality
$\pm(a_{x}^{\dagger}a_{y}+a_{y}^{\dagger}a_{x})\leq n_{x}+n_{y}$ we get

$\displaystyle I(A:B)\leq\beta J\sum_{\begin{subarray}{c}x\sim y,\\\ x\in
A,y\in
B\end{subarray}}\operatorname{tr}((n_{x}+n_{y})(\rho_{A}\otimes\rho_{B}+\rho_{AB})).$

Since, for $x\in A$,
$\operatorname{tr}(n_{x}\rho_{A}\otimes\rho_{B})=\operatorname{tr}(n_{x}\rho_{A})=\operatorname{tr}(n_{x}\rho_{AB})$
and similarly for $n_{y}$ and $\rho_{B}$, we obtain

$\displaystyle I$ $\displaystyle(A:B)\leq 2\beta
J\sum_{\begin{subarray}{c}x\sim y,\\\ x\in A,y\in
B\end{subarray}}\operatorname{tr}((n_{x}+n_{y})\rho_{AB}).$ (S4)

Observe that $H_{AB}$ is translation-invariant, i.e., for all
$x\in\Lambda_{L}$,

$\mathcal{T}_{x}^{*}H_{AB}\mathcal{T}_{x}=H_{AB},$

where $\mathcal{T}_{x}$ denotes the unitary translation operator by $x$,
$(\mathcal{T}_{x}\psi)(y)=\psi(y+x\mod L)$. This implies
$\mathcal{T}_{x}^{*}\rho_{AB}\mathcal{T}_{x}=\rho_{AB}$ and therefore,
$\operatorname{tr}(n_{x}\rho_{AB})=\operatorname{tr}(n_{y}\rho_{AB})$ for all
$x,y$. The summation in (S4) is over $2L^{d-1}$ terms, so

$\displaystyle\sum_{x\sim y,\leavevmode\nobreak\ x\in A,y\in
B}\operatorname{tr}((n_{x}+n_{y})\rho_{AB})=4L^{d-1}\operatorname{tr}(n_{x_{0}}\rho_{AB})=4L^{d-1}\frac{\operatorname{tr}(\mathcal{N}\rho_{AB})}{L^{d}},$

where $x_{0}\in\Lambda$ is some arbitrary element. This proves Proposition 6.
∎

###  VI.2 Step 2: Removing the interaction from the exponential

In the following we write
$\left<P\right>_{K}:=\operatorname{tr}(Pe^{-K})/\operatorname{tr}(e^{-K})$ for
operators $K$ and $P$ and we will also drop the subscript $AB$ and write
$H=H_{AB}$.

###### Proposition 7

Let $C_{0}=\frac{1}{4}\frac{U}{\gamma+\mu}$. For all $\gamma>0$, we have

$\displaystyle\left<\mathcal{N}\right>_{\beta H}\leq
2\bigg{(}\frac{C_{0}}{U/2}\left<W-(\mu+\gamma-2dJ)\mathcal{N}\right>_{\beta(H_{0}+\gamma\mathcal{N})}+\frac{C_{0}}{4}(1+C_{0}^{-1})^{2}L^{d}\bigg{)}.$

###### Proof of Proposition 7.

Let $C>0$. We have $n\leq Cn(n-1)+\frac{C}{4}(1+C^{-1})^{2}$ for all
$n\in{\mathbb{N}}$. Thus, we obtain

$\displaystyle\left<\mathcal{N}\right>_{\beta
H}=\left<\sum_{x\in\Lambda_{L}}n_{x}\right>_{\beta
H}\leq\sum_{x\in\Lambda_{L}}\left<Cn_{x}(n_{x}-1)+\frac{C}{4}(1+C^{-1})^{2}\right>_{\beta
H}=\frac{C}{U/2}\left<W\right>_{\beta H}+\frac{C}{4}(1+C^{-1})^{2}L^{d}.$

This leads to

$\displaystyle\left<\mathcal{N}\right>_{\beta
H}\left(1-\frac{C}{U/2}(\mu+\gamma)\right)$
$\displaystyle\leq\frac{C}{U/2}\left<W-(\mu+\gamma)\mathcal{N}\right>_{\beta
H}+\frac{C}{4}(1+C^{-1})^{2}L^{d}$
$\displaystyle\leq\frac{C}{U/2}\left<W-(\mu+\gamma-2dJ)\mathcal{N}\right>_{\beta
H}+\frac{C}{4}(1+C^{-1})^{2}L^{d}.$

We now apply twice the Peierls-Bogoliubov inequality, i.e., (3) in Lemma 4
with $P=\beta(W-(\mu+\gamma-2dJ)\mathcal{N})$ and $K=-\beta H$, and get

$\displaystyle\left<W-(\mu+\gamma-2dJ)\mathcal{N}\right>_{\beta
H}\leq\left<W-(\mu+\gamma-2dJ)\mathcal{N}\right>_{\beta(H_{0}+\gamma\mathcal{N})}.$

Using this in the previous bound and setting $C=C_{0}$ yields the desired
estimate. This proves Proposition 7. ∎

###  VI.3 Step 3: Calculation for quasi-free states

The estimate of the on-site interaction in the free reference states leads to
an estimate of the particle number on a specific site via the one-particle
density matrix. Here we get a Riemann sum of the density function of Planck’s
law.

We set

$f(\gamma,\beta,J):=\int_{[0,\frac{1}{2}]^{d}}\left(e^{4J\beta\sum_{j=1}^{d}\sin^{2}(\pi
x_{j})+\beta\gamma}-1\right)^{-1}{\mathrm{d}}x,$

and denote the error terms by

$\displaystyle\epsilon_{1}$ $\displaystyle=\frac{(L+1)^{d}}{L^{d}}-1,$
$\displaystyle\epsilon_{2}$
$\displaystyle=\frac{(e^{\beta\gamma}-1)^{-1}}{L^{d}}\frac{(2(2L+1))^{d-1}}{2L-1}.$

We shall use the following two lemmas.

###### Lemma 8

For all $\gamma>0$ and $x\in\Lambda_{L}$,

$\displaystyle\left<a_{x}^{\dagger}a_{x}\right>_{\beta(H_{0}+\gamma\mathcal{N})}\leq$
$\displaystyle 2^{d}(1+\epsilon_{1})f(\gamma,\beta,J)+\epsilon_{2}.$

###### Lemma 9

For all $\gamma>0$ and $\alpha\in(0,1)$,

$f(\gamma,\beta,J)\leq\frac{1}{2^{d}}\frac{e^{-\alpha\beta\gamma}}{(1-\alpha)\beta\gamma}.$

We postpone the proofs of these lemmas for now and show how they imply they
main result.

###### Proof of Theorem 2.

Wick’s theorem for quasi-free states [68] yields

$\displaystyle\left<n_{x}(n_{x}-1)\right>_{{\beta(H_{0}+\gamma\mathcal{N})}}=\left<(a_{x}^{\dagger}a_{x})^{2}\right>_{\beta(H_{0}+\gamma\mathcal{N})}=2\left<a_{x}^{\dagger}a_{x}\right>_{\beta(H_{0}+\gamma\mathcal{N})}^{2}.$

By Lemma 8,

$\displaystyle\left<W\right>_{\beta(H_{0}+\gamma\mathcal{N})}=2L^{d}\left<a_{x}^{\dagger}a_{x}\right>_{\beta(H_{0}+\gamma\mathcal{N})}^{2}\leq
2^{2d+2}L^{d}(1+\epsilon_{1})^{2}f(\gamma,\beta,J)^{2}+4L^{d}\epsilon_{2}^{2}.$

Using this and
$\left<-(\mu+\gamma-2dJ)\mathcal{N}\right>_{\beta(H_{0}+\gamma\mathcal{N})}\leq
0$ in the upper bound of Proposition 7 , we obtain that for all
$\gamma>2dJ-\mu$,

$\displaystyle\left<\mathcal{N}\right>_{\beta H}\leq
2L^{d}\bigg{(}\frac{C_{0}}{U/2}(2^{2d+2}(1+\epsilon_{1})^{2}f(\gamma,\beta,J)^{2}+4\epsilon_{2}^{2})+\frac{C_{0}}{4}(1+C_{0}^{-1})^{2}\bigg{)}.$

Next we use that
$f(\gamma,\beta,J)\leq\frac{1}{2^{d}}\frac{e^{-\alpha\beta\gamma}}{(1-\alpha)\beta\gamma}$
for all $\alpha\in(0,1)$, cf. Lemma 9. The upper bound for
$\left<\mathcal{N}\right>_{\beta H}$ then becomes

$\displaystyle
2L^{d}\bigg{(}\frac{1}{2}\frac{1}{\gamma+\mu}\left(4(1+\epsilon_{1})^{2}\frac{e^{-2\alpha\beta\gamma}}{(1-\alpha)^{2}\beta^{2}\gamma^{2}}+4\epsilon_{2}^{2}\right)+\frac{1}{16}\frac{U}{\gamma+\mu}+\frac{\gamma+\mu}{U}+\frac{1}{2}\bigg{)}.$

Estimating $\epsilon_{1}\leq 1$ and $\epsilon_{2}\leq
6^{d-1}(e^{\beta\gamma}-1)^{-1}$, we arrive at the bound

$\displaystyle
2L^{d}\bigg{(}\frac{1}{2(\gamma+\mu)}\bigg{(}16\frac{e^{-2\alpha\beta\gamma}}{(1-\alpha)^{2}\beta^{2}\gamma^{2}}+4\cdot
6^{2(d-1)}(e^{\beta\gamma}-1)^{-2}\bigg{)}+\frac{1}{16}\frac{U}{\gamma+\mu}+\frac{\gamma+\mu}{U}+\frac{1}{2}\bigg{)}.$

Finally, choosing $\gamma=\max\\{1/\beta,2dJ+1\\}$ (such that $\beta\gamma\geq
1$), considering the limit $\alpha\to 0$ and in view of Proposition 6, the
theorem holds with

$\displaystyle c(J,U,\mu)=16J\bigg{(}\frac{16+4\cdot
6^{2(d-1)}(e-1)^{-2}}{2(2dJ+1+\mu)}+\frac{1}{16}\frac{U}{2dJ+1+\mu}+\frac{2dJ+1+\mu}{U}+\frac{1}{2}\bigg{)}.$
(S5)

Distinguishing the cases for $\beta$ completes the proof of Theorem 2. ∎

###  VI.4 Laplace eigenvalues and eigenvectors for periodic boundary
conditions

The proofs proof of Lemma 8 requires information on the spectral theory of the
one-body graph Laplacian, for which introduce notation here.

Consider the discrete Laplacian $-\Delta\geq 0$ on the chain
$\\{0,1,\ldots,L\\}$, $L\in{\mathbb{N}}$, with periodic boundary conditions,
i.e., assume that the nodes $0$ and $L$ are connected and every vector on
$\ell^{2}(\\{0,1,\ldots,L\\})$ satisfies $v(0)=v(L)$. So it suffices to use
vectors $v\in\ell^{2}(\\{1,\ldots,L\\})$.

The $L$ eigenvalues $\lambda_{i}$ and eigenvectors $v_{i}$, $i=1,\ldots,L$ of
$-\Delta$ in this situation are given by

$\displaystyle\lambda_{2k+1}$
$\displaystyle=4\sin^{2}\left(\frac{k\pi}{L}\right),\leavevmode\nobreak\
k=0,\ldots,\left\lfloor\frac{L-1}{2}\right\rfloor,$
$\displaystyle\lambda_{2k}$
$\displaystyle=4\sin^{2}\left(\frac{k\pi}{L}\right),\leavevmode\nobreak\
k=1,\ldots,\left\lfloor\frac{L}{2}\right\rfloor,$

and

$\displaystyle v_{1}(i)$ $\displaystyle=L^{-1/2},$ $\displaystyle v_{L}(i)$
$\displaystyle=L^{-1/2}(-1)^{i}\text{ if }L\text{ even},$ $\displaystyle
v_{2k+1}(i)$ $\displaystyle=\sqrt{2/L}\cos\left(\frac{\pi
k}{L}(2i-1)\right),\leavevmode\nobreak\
k=1,\ldots,\left\lfloor\frac{L-1}{2}\right\rfloor,$ $\displaystyle v_{2k}(i)$
$\displaystyle=\sqrt{2/L}\sin\left(\frac{\pi
k}{L}(2i-1)\right),\leavevmode\nobreak\
k=1,\ldots,\left\lfloor\frac{L-1}{2}\right\rfloor.$

###  VI.5 Proofs of Lemmas 8 and 9

In this section, we give the still outstanding proofs of Lemmas 8 and 9.

###### Proof of Lemma 8.

Let $x_{0}=(1,\ldots,1)\in\Lambda_{L}$. By the formula for the one-particle
density matrix, cf. e.g. [68, Prop. 5.2.28] and by translation-invariance we
find

$\displaystyle\left<a_{x}^{\dagger}a_{x}\right>_{\beta(H_{0}+\gamma\mathcal{N})}$
$\displaystyle=\left<\delta_{x},e^{-\beta\gamma}e^{-\beta(-J\Delta)}(\mathds{1}-e^{-\beta\gamma}e^{-\beta(-J\Delta)})^{-1}\delta_{x}\right>$
$\displaystyle=\left<\delta_{x_{0}},(e^{\beta(-J\Delta+\gamma)}-\mathds{1})^{-1}\delta_{x_{0}}\right>$
$\displaystyle=\sum_{i_{1},\ldots,i_{d}=1}^{L}(e^{\beta(J\sum_{j=1}^{d}\lambda_{i_{j}}+\gamma)}-1)^{-1}\prod_{j=1}^{d}\left|v_{i_{j}}(1)\right|^{2}.$

Let $\ell=\lceil\frac{L}{2}\rceil$. Then we have for any $c>0$

$\displaystyle\sum_{i=1}^{L}(e^{\beta(J\lambda_{i}+c)}-1)^{-1}\left|v_{i}\right|^{2}$
$\displaystyle=\frac{1}{L}(e^{\beta
c}-1)^{-1}+\sum_{i=1}^{\ell}(e^{\beta(J\lambda_{2i}+c)}-1)^{-1}\bigg{(}\left|v_{2i}(1)\right|^{2}+\mathds{1}_{L=2\ell+1}\left|v_{2i+1}(1)\right|^{2}\bigg{)}$
$\displaystyle\leq\frac{1}{L}(e^{\beta
c}-1)^{-1}+\frac{2}{L}\sum_{i=1}^{\ell}(e^{\beta(J\lambda_{2i}+c)}-1)^{-1}.$

Then, by induction over $d$ we find as an upper bound

$\displaystyle\left<a_{x}^{\dagger}a_{x}\right>_{\beta(H_{0}+\gamma\mathcal{N})}\leq$
$\displaystyle\frac{1}{L^{d}}(e^{\beta\gamma}-1)^{-1}+\frac{1}{L^{d}}\sum_{k=1}^{d-1}4^{k}\sum_{i_{1},\ldots,i_{k}=1}^{\ell}(e^{\beta(J\sum_{j=1}^{k}\lambda_{2i_{j}}+\gamma)}-1)^{-1}$
(S6)
$\displaystyle+\frac{2^{d}}{L^{d}}\sum_{i_{1},\ldots,i_{d}=1}^{\ell}(e^{\beta(J\sum_{j=1}^{d}\lambda_{2i_{j}}+\gamma)}-1)^{-1}.$
(S7)

Using $2\ell-1\leq L\leq 2\ell$ the terms in (S6) can be estimated by

$\displaystyle\frac{(e^{\beta\gamma}-1)^{-1}}{L^{d}}\sum_{k=0}^{d-1}(4\ell)^{k}$
$\displaystyle=\frac{(e^{\beta\gamma}-1)^{-1}}{L^{d}}\frac{(4\ell)^{d-1}-1}{4\ell-1}\leq\epsilon_{2}.$

The last term (S7) represents a Riemann sum of a function, decreasing in each
argument, which is evaluated at the minimal points of the hypercubes
$[(i_{1}-1)\pi/2,i_{1}\pi/2]\times\ldots\times[(i_{d}-1)\pi/2,i_{d}\pi/2]$,
i.e.,

$\displaystyle\frac{2^{d}}{L^{d}}\sum_{i_{1},\ldots,i_{d}=1}^{\ell}(e^{\beta(J\sum_{j=1}^{d}\lambda_{2i_{j}}+\gamma)}-1)^{-1}$
$\displaystyle=\frac{(2\ell)^{d}}{L^{d}\ell^{d}}\sum_{i_{1},\ldots,i_{d}=1}^{\ell}(e^{\beta(4J\sum_{j=1}^{d}\sin^{2}(i_{j}\pi/L)+\gamma)}-1)^{-1}$
$\displaystyle\leq\frac{(L+1)^{d}}{L^{d}\ell^{d}}\sum_{i_{1},\ldots,i_{d}=1}^{\ell}(e^{\beta(4J\sum_{j=1}^{d}\sin^{2}(i_{j}\pi/(2\ell))+\gamma)}-1)^{-1}$
$\displaystyle\leq\frac{(L+1)^{d}}{L^{d}}\int_{[0,1]^{d}}(e^{\beta(4J\sum_{j=1}^{d}\sin^{2}(x_{j}\pi/2)+\gamma)}-1)^{-1}{\mathrm{d}}x$
$\displaystyle=2^{d}\frac{(L+1)^{d}}{L^{d}}f(\gamma,\beta,J).\qed$

###### Proof of Lemma 9.

We have

$\displaystyle\int_{[0,\frac{1}{2}]^{d}}(e^{4J\beta\sum_{j=1}^{d}\sin^{2}(\pi
x_{j})+\beta\gamma}-1)^{-1}{\mathrm{d}}x$
$\displaystyle=e^{-\alpha\beta\gamma}\int_{[0,\frac{1}{2}]^{d}}(e^{4J\beta\sum_{j=1}^{d}\sin^{2}(\pi
x_{j})+(1-\alpha)\beta\gamma}-e^{-\alpha\beta\gamma})^{-1}{\mathrm{d}}x$
$\displaystyle\leq
e^{-\alpha\beta\gamma}\int_{[0,\frac{1}{2}]^{d}}(e^{4J\beta\sum_{j=1}^{d}\sin^{2}(\pi
x_{j})+(1-\alpha)\beta\gamma}-1)^{-1}{\mathrm{d}}x$ $\displaystyle\leq
e^{-\alpha\beta\gamma}\int_{[0,\frac{1}{2}]^{d}}((1-\alpha)\beta\gamma)^{-1}{\mathrm{d}}x=\frac{1}{2^{d}}\frac{e^{-\alpha\beta\gamma}}{(1-\alpha)\beta\gamma}.$

This proves Lemma 9 ∎