Thermal Area Law for Lattice Bosons

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.


Introduction
In quantum many-body systems with translationinvariant short-ranged interactions the entanglement entropy of the ground state typically 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]- [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 higherdimensional situation, see [7], [8].
The analog of the area law for Gibbs states is called the thermal area law.Gibbs states are important for fundamental reasons and their efficient simulability via tensor networks hinges on thermal area laws [34]- [36].For Gibbs states, the total correlations between two regions A and B are quantified by their mutual information [34], [37], [38] I(A : B) = S(ρ A ) + S(ρ B ) − S(ρ AB ), where S denotes the von Neumann entropy, ρ AB the Gibbs state of the full system, and ρ A , ρ B the reduced density matrices corresponding to A and B, respectively.At zero temperature, I(A : B) reduces to twice the entanglement entropy.An area law at positive temperature was derived in a seminal work of Wolf et.al [39] who proved for local and bounded interactions that with ∂ AB being the boundary region between A and B, and whenever A and B are disjoint and together make up the entire lattice.Very recently, the βscaling result was improved to β 2/3 [34] which matters for the experimentally relevant regime of low temperatures (β → ∞).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 β 1/5 for large β [40].Moreover, a generalization to various Rényi generalizations of the mutual information was given in [35] and this has important applications to simulating and approximating Gibbs states [36].Further results of thermal area laws were established for free fermions [41], for the entanglement negativity (instead of the mutual information) [42], as a result of rapid mixing for dissipative quantum lattice systems [43], [44] (with a logarithmic correction) and numerically for some spin chains showing a log β-dependence [45].A current account of thermal area laws and related phenomena is given in [38].
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 [46].Thermal area laws provide us with universal properties of Gibbs states independent of the system size.They are especially valuable for lower temperatures, while at high temperatures many analytic properties are available, e.g.exponential decay of correlations [47]- [52], the large deviation principle [53]- [55] or the approximate quantum Markov property [56], [57].Such characterizations are important for the computation of Gibbs states, in general an NP-hard problem [58], [59], which in turn is fundamental in novel applications such as quantum machine learning [60], [61], semidefinite programming solvers [62], [63] or the imaginary-time evolution in the framework of near-term quantum devices [64]- [70].
In particular, area laws can be used for the approximation of Gibbs state by matrix product operators (MPO) and their higher-dimensional analogs [34], [36], [38], [71].In special cases like one-dimensional quantum spin systems they gave rise to concrete algorithms for the efficient approximation by MPO, where the temperature behavior in the area law determines the maximal bond dimension [34].Such algorithms also proved to be useful for the representation of Gibbs states as a convex combination of MPS [72], or for the approximation of ground states under a low-energy-density assumption typically observed for gapped systems.It is generally believed that the scaling of the mutual information with β in the low temperature regime is related to the computational complexity of the ground space of the models [38].Another direct way to use a thermal area law is to note that controlling the mutual information automatically controls all standard correlation functions, see [39] and Appendix A.3 below.
A critical limitation of the existing results is that they only hold for bounded interactions and 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 three main reasons: (i) The Hamiltonians can be experimentally fine-tuned for cold atoms in optical lattices [73], which makes them promising platform for quantum simulation and quantum engineering, see also [74].(ii) Bosonic encoding can be used in quantum information processing which can provide multiple advantages over finite-dimensional discrete-variable (DV) codes [75], [76].(iii) 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], [77], but the case of interacting bosons, including the paradigmatic Bose-Hubbard model proved elusive to rigorous analysis.A partly numerical investigation was given in [78] with a focus on the phase transition from Mott insulator to superfluid by analogy with other symmetry breaking transitions [79]- [81].Recently, Abrahamsen et al. rigorously proved an area law for gapped ground states of 1D bosonic lattice Hamiltonians in [82] by using a truncation of the local Hilbert spaces and a quantum number tail bound from [83].Their work only concerns the zero temperature case and leaves open the positive temperature case.
In this work 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 The precise result is Theorem 2 below.The β-scaling is the same as that found by [39] 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 quasifree states that can be completed with Wick's rule.Further details are explained below and in the appendix.

Setup for infinite-dimensional Hilbert spaces
A technical point in the description of many-boson systems is that the local Hilbert spaces are infinitedimensional 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. [39] 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 h be the local Hilbert space for one site, Λ a finite set and A, B ⊆ Λ such that A ∩ B = ∅ and A ⊔ B = Λ.We set We suppose that the Hamiltonian can be decomposed as For details about operator domains, which are relevant because our operators are unbounded, see Appendix A.1.Let T + 1 (H) denote the set of all density matrices on H. Then we define the free energy as where The Gibbs state minimizes the free energy, i.e., F β AB = F β (ρ AB ).(See Proposition 4 for a proof of this in the infinitedimensional setting.)We reduce to the A, respectively B subsystem by taking partial traces of the Gibbs state The general bound from [39] straightforwardly extends to the infinite-dimensional setting as follows.
Lemma 1 (Boundary energy controls mutual informat.)Let β > 0 and suppose that tr e −βH AB < ∞.Then Proof.We start from Using S(ρ A ⊗ ρ B ) = S(ρ A ) + S(ρ B ), we obtain and the right-hand side equals ) by basic properties of the partial trace.

Main Result
Now we consider the Bose-Hubbard model in the framework of the previous section.Let Λ ≡ Λ L denote a box of side length L in the d-dimensional lattice with periodic boundary conditions.For x, y ∈ Λ L we write x ∼ 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 h = ℓ 2 (N).
The total Hilbert space H AB = ⊗ x∈Λ L h is isomorphic to the Fock space F(ℓ 2 (Λ L )).On it, we consider the Bose-Hubbard Hamiltonian where J ∈ R represents the strength of the kinetic nearest-neighbor hopping, U > 0 the strength of the on-site repulsion and µ ∈ R the chemical potential, and N = x∈Λ L n x is the total number operator.The Hamiltonian is self-adjoint on a suitable domain D(H AB ); see e.g., [84] and Appendix A.1.
We are now ready to state the main result.We decompose the box into two regions A and B with boundary region as shown in Figure 1.
Theorem 2 (Main result: thermal area law) For all β, U, µ > 0, we have 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 µ > 0 is standard, cf. the usual phase diagram in Figure 13(a) in [73].Indeed, note that if one would take µ sufficiently negative, then the system becomes effectively devoid of particles in the grand-canonical setting.(iii) The constant c(J, U, µ) can be made explicit from (11).(iv) The maximum max{1, β} means that the bound behaves as β for low temperatures in accordance with [39].We recall that some growth in β is strongly expected without any gap assumption [40].For high temperature β < 1, we find a temperature independent lower bound which matches the classical situation [39].(v) As explained in [39], 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, see Appendix A.3.
We close the presentation by discussing several extensions of the result which can be obtained from the same methods.The proof can be extended to Hamiltonians of the form where M ∈ N, f ≥ 0 is a polynomially bounded function, growing faster than x M , and J ν1...ν l x1...x k y1...y k is uniformly bounded and finite-range.Moreover, the proof also works if we add to the original Hamiltonian density-density interactions of the form provided that U is sufficiently large, and J xy is uniformly bounded and finite-range.Furthermore, instead of finite-range interactions, we can consider hopping terms decaying at infinity fast enough, e.g., − x,y J xy a † x a y for J xy ≥ 0, x ∈ Z d , satisfying Finally, 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.
The translation-invariance of the underlying lattice is in fact necessary in our current proof.Dropping this assumption would require to control the number of bosons potentially accumulating on the boundaries of the whole system, see also the proof of Proposition 5.It is an interesting open problem to remove the translation-invariance assumption.Furthermore, as we control the H 0 term with the on-site interaction W , we also need the rather strong decay assumption on J xy .Therefore, an interesting problem is to develop an alternative approach allowing for long-range interactions with α arbitrary close to d.

Sketch of proof of Theorem 2
The detailed proof of Theorem 2 is given in Appendix B.Here we give a sketch of the main ideas.
We decompose the Hamiltonian as where H 0 = −J x∼y a † x a y + 2dJN is the shifted kinetic term and W = U x a y + a † y a x ).Notice that the boundary Hamiltonian H ∂ contains at most d |∂ AB | many summands, so the right-hand side of Lemma 1 seems to exhibits the desired scaling in L and β.The main challenge is that the hopping terms a † x a y + a † y a x between A and B 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 −β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 [85, (2.14)], which has a long history in the study of quantum many-body systems and quantum information theory.Applying it twice, we obtain We use this bound with K = −βH AB and P = β(W − (µ+γ −2dJ)N ) so that the new effective Hamiltonian is and e −β(H0+(γ+2dJ)N ) is a trace-class quasifree state, so that expectations can be calculated via Wick's rule.Here we introduced a parameter γ > 0 large enough in order to make e −β(H0+(γ+2dJ)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 γ, we would get e −β(H0−(µ+2dJ)N ) which has infinite trace for µ > −2dJ.) The final expression that we arrive at via (4) 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 β in an integrable way.In particular, we see that the resulting expression as well as the error terms are bounded for β ≥ 1 and Theorem 2 follows.For the details, see Appendix B.
For the discussed generalizations (3), one has to bound H 0 by W and then remove the generalized hopping term H 0 instead of W in the exponent by means of the Peierls-Bogoliubov argument.In the end, one obtains a trace just involving number operators, which can be easily computed as well.

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 positivetemperature counterpart to the area law for gapped bosonic ground states (in 1D) [82].
The idea of the proof is based on the general idea in [39] 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 [71].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 [35]), see also [40].Furthermore, in light of recent progress of Lieb-Robinson bounds for bosons [84], [86]- [91] one could explore the applicability of such bounds and their imaginary time counterparts in the context of thermal area laws, see also [34], [44] for connections between imaginary time Lieb-Robinson bounds and thermal area laws.
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 [60], [61] and semidefinite programming solvers [62], [63].

A.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 D(A) and D(B), respectively and the boundary Hamiltonian is defined on H ∂ on D(A) ⊗ D(B).We assume that H AB is essentially self-adjoint on D(A) ⊗ D(B) which is then the appropriate domain for the equality Concerning the Bose-Hubbard model, in the total Fock space F(ℓ 2 (Λ L )), we consider the dense domain The operator H AB is self-adjoint on the largest domain D(H AB ), where it can be defined, cf.[84].It then follows from [H AB , N ] = 0 that H AB is indeed essentially self-adjoint on F fin (ℓ 2 (Λ L )).Moreover, for subsystems X ∈ {A, B}, we always set D(X) = F fin (H X ).The above properties can then be verified.

A.2 Trace inequalities in infinite dimensions
Lemma 3 (Peierls-Bogoliubov inequality) Let N ≥ 0 be self-adjoint operator with purely discrete spectrum, let Π N := 1 N ≤N , N ∈ N and assume that Π N H is finite-dimensional for all N ∈ N. Let (D(K), K) be a self-adjoint and (D(P ), P ) be a symmetric operator such that [K, N ] = 0, [P, N ] = 0, K + P is self-adjoint, and e K , P e K , e K+P are all trace-class.Then we have In particular, if P e K+P is trace-class as well, we obtain (4), i.e., tr(P e K ) tr(e K ) ≤ tr(P e K+P ) tr(e K+P ) .
Proof of Lemma 3. The inequality ( 6) is well-known for matrices, see for example [85, (2.14)] (or more generally if P is bounded [92]).Therefore, we have where Π N = 1 N ≤N .Taking the limit N → ∞ yields the desired result.Finally, (4) follows from ( 6) by means of log tr(e K+P ) tr(e K ) = − log tr(e K ) tr(e K+P ) ≤ tr(P e K+P ) tr(e K+P ) .

A.3 Truncated correlations bound
Let M A , M B be two bounded self-adjoint operators on H A and H B , respectively.We denote their truncated correlation function as For all β, U, µ > 0, we have with the same constant as in Theorem 2. This follows from the quantum Pinsker inequality I(A : B) ≥ and ∥X∥ 1 ≥ tr(XY )/ ∥Y ∥.We mention that the standard proof of the quantum Pinkser inequality, cf.[95,Theorem 1.15] extends to infinite dimensions since the data processing inequality holds for trace-class operators [96].
Note however that (8) is trivial unless |∂ AB | stays bounded in the infinite-volume limit, e.g. if A or B are kept at fixed size, as its left-hand side is always bounded by one.Therefore, in order to find more regimes where this is useful, it is an interesting open question if our main estimate in Theorem 2 can be improved for β → 0, as it is the case for spin systems [34], [39].

B 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 → −a x at every second lattice site to reduce to the case J > 0.

B.1 Step 1: Controlling boundary energy by particle number
In this section, we prove the following bound

Proposition 5
We have Proof of Proposition 5. Using Lemma 1 and the operator Cauchy-Schwarz inequality ±(a † x a y + a † y a x ) ≤ n x + n y we get Since, for x ∈ A, tr(n x ρ A ⊗ ρ B ) = tr(n x ρ A ) = tr(n x ρ AB ) and similarly for n y and ρ B , we obtain Observe that H AB is translation-invariant, i.e., for all x ∈ Λ L , where T x denotes the unitary translation operator by x, (T x ψ)(y) = ψ(y + x mod L).This implies T * x ρ AB T x = ρ AB and therefore, tr(n x ρ AB ) = tr(n y ρ AB ) for all x, y.The summation in ( 9) is over |∂ AB | many terms, so where x 0 ∈ Λ is some arbitrary element.This proves Proposition 5.

B.2
Step 2: Removing the interaction from the exponential In the following we write ⟨P ⟩ K := tr(P e −K )/ tr(e −K ) for operators K and P and we will also drop the subscript AB and write H = H AB .

B.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 and denote the error terms by We shall use the following two lemmas.

Lemma 7
For all γ > 0 and

Lemma 8
For all γ > 0, We postpone the proofs of these lemmas for now and show how they imply they main result.
Remark 9. Notice that the error terms ϵ 1 and ϵ 2 converge to zero as L → ∞.It is of separate interest whether the remaining term 2 d f (γ, β, J) actually captures the correct asymptotic behavior of a † x a x β(H0+γN ) , cf.Lemma 7.
The L eigenvalues λ i and eigenvectors v i , i = 1, . . ., L of −∆ in this situation are given by and B.5 Proofs of Lemmas 7 and 8 In this section, we give the still outstanding proofs of Lemmas 7 and 8.
Proof of Lemma 7. Let x 0 = (1, . . ., 1) ∈ Λ L .By the formula for the one-particle density matrix, cf.e.g.[94,Prop. 5.2.28] and by translation-invariance we find where the sum over the i 1 , . . ., i k is defined to be one if k = 0.This can be proven like the binomial theorem, more precisely, Then taking α → 0 proves Lemma 8.

Figure 1 :
Figure 1: Periodic box for d = 2 and L = 20 decomposed into two regions A and B. The green region shows the boundary region ∂AB.The bonds connecting A and B make up the boundary Hamiltonian H ∂ .

2
x∈X n x (n x − 1) is the onsite interaction.To use Lemma 1, we again decompose the Hamiltonian asH AB = H A + H B + H ∂where we define the subsystem Hamiltonians with open boundary conditions along the cut.More precisely, for X ∈ {A, B}, we set where ρ H = e −βH / tr e −βH .Proof of Proposition 4. We have tr(K ln K − K ln P ) ≥ tr(K − P ) for all positive self-adjoint trace class operators K, P [93, Prop.2.5.3].This yields [94, Section 5.3.1]