A subpolynomial-time algorithm for the free energy of one-dimensional quantum systems in the thermodynamic limit

We introduce a classical algorithm to approximate the free energy of local, translation-invariant, one-dimensional quantum systems in the thermodynamic limit of infinite chain size. While the ground state problem (i.e., the free energy at temperature $T = 0$) for these systems is expected to be computationally hard even for quantum computers, our algorithm runs for any fixed temperature $T>0$ in subpolynomial time, i.e., in time $O((\frac{1}{\varepsilon})^{c})$ for any constant $c>0$ where $\varepsilon$ is the additive approximation error. Previously, the best known algorithm had a runtime that is polynomial in $\frac{1}{\varepsilon}$. Our algorithm is also particularly simple as it reduces to the computation of the spectral radius of a linear map. This linear map has an interpretation as a noncommutative transfer matrix and has been studied previously to prove results on the analyticity of the free energy and the decay of correlations. We also show that the corresponding eigenvector of this map gives an approximation of the marginal of the Gibbs state and thereby allows for the computation of various thermodynamic properties of the quantum system.


Introduction
Multipartite quantum systems are described by a Hilbert space, which is a tensor product of the single-particle d-dimensional spaces. The behaviour of a quantum-many body system is described by a Hamiltonian which models the interaction between the different particles. Of particular interest are k-local Hamiltonians that can be written as a sum of terms acting nontrivially on at most k particles, with k being a constant. At thermal equilibrium, the system is described by the Gibbs state where β = 1/T is the inverse temperature, and Z β (H) = tr e −βH is the partition function. The free energy of the system at inverse temperature β is defined as At zero temperature, i.e., β = +∞, F β (H) becomes λ min (H), the ground energy of H. The problem of computing the ground energy for a given local Hamiltonian is known to be QMA-complete [1], and is the central problem in the area of Hamiltonian complexity [2]. This problem remains QMA-complete even if we restrict ourselves to 2-local Hamiltonians that are translation-invariant on a chain [3,4], i.e., H = i h i,i+1 where the operators h i,i+1 are given by some Hermitian operator h (the same one for every i) acting on particles i and i + 1. 1 In order to understand the physical properties of the system at nonzero temperature, it is crucial to understand not only the ground energy, but also the free energy function F β (H) as a function of β > 0 [5]. Indeed, computing F β (H) and its derivatives with respect to β and parameters of the Hamiltonian determines phase transitions and gives access to fundamental physical properties of the system in thermal equilibrium such as the internal energy, specific heat, or magnetic susceptibility [6].

Main result
In this paper, we focus on 2-local translation-invariant quantum systems on an infinite chain. As the free energy scales with the system size, in the thermodynamic limit of infinite systems we consider the free energy per particle f β (h). Note that f β (h) only depends on the finite matrix h of size d 2 × d 2 . Our objective is to design an algorithm to approximate f β (h) with a good scaling in terms of the target error ε and the local dimension d. As argued in recent works on Hamiltonian complexity in the thermodynamic limit [7,8], understanding the dependence of the complexity in terms of the desired precision for infinite systems is often closer to capturing the fundamental problems in many-body physics than understanding the dependence in the system size. Our main result is an algorithm that given as input h and a target error ε outputs an approximation of f β (h) and of the k-particle marginals of the Gibbs state.

Theorem 1.1. There is a deterministic algorithm that takes as input a Hermitian operator
h acting on C d ⊗ C d satisfying h ≤ 1 and ε ∈ (0, 1/e) and outputs an approximatioñ f β satisfying |f β − f β (h)| ≤ ε, where f β (h) is the free energy per particle of the infinite translation-invariant Hamiltonian on a chain defined by h. For any fixed β > 0, the running time of the algorithm is exp O log d log(1/ε) log log (1/ε) . Moreover, this algorithm can also compute an ε-approximation of the marginal of the Gibbs state on an interval of size k with the same running time for any fixed k.
Before describing the algorithm and proof method, we make some remarks and discuss related works.
Remarks We note that the temperature dependence of the running time is hidden in the O(.) notation as we are interested in the algorithm for fixed temperature. If we want to make the dependence on the inverse temperature β explicit, the running time takes the form exp O log d log(1/ε) log log(1/ε) exp(O(β)) , where O(.) only hides universal constants. 2 We note that, from the QMA-hardness results for the ground energy problem [3], it follows that we should not expect a polynomial dependence in β and 1/ε and as such an exponential scaling in β is unavoidable in our algorithm.
To appreciate the algorithm we use to prove Theorem 1.1, it is instructive to consider first a naive algorithm for this problem sometimes called exact diagonalization. The idea of this algorithm is to consider the Hamiltonian H [1,n] = n−1 i=1 h i,i+1 restricted to only n particles. Computing the free energy per particle f β,n of H [1,n] (for any value of β) can √ log(n/ε)) , where theÕ depends linearly on β. As with the algorithm of [12], it can be readily applied to the infinite translation-invariant chain by setting n = 1/ε. This leads to an algorithm that is quasi-linear in 1/ε for computing thermal expectation values for infinite translation-invariant chains. This algorithm can be used to compute the free energy using the argument in [17,Lemma 12] at the cost of an additional polynomial overhead in 1/ε, thus leading to an algorithm for computing the free energy for infinite translation-invariant chains that is polynomial in 1/ε.
In a different regime, classical algorithms have been designed in [17] to compute the free energy of dense Hamiltonians based on convex relaxations. These algorithms have a runtime that is exponential in 1/ε.
Besides this line of work on provably convergent algorithms to which our work shall also contribute, there are numerous algorithms that effectively address the problem despite having no convergence results or only in special cases. For the free energy problem this includes most notably Quantum Monte Carlo methods [18,19]. These probabilistic algorithms lack rigorous results on their runtime except for few special cases and are known to fail for Hamiltonians that have the so-called sign problem. Another example of effective algorithms for the ground state energy problem (β = +∞) in one dimension are tensor networks and the DMRG algorithm [20,21]. The convergence of a related algo-rithm to the ground state energy has been proven under the additional assumption that the Hamiltonian is gapped [22].
Proof technique Before giving an overview of the algorithm establishing Theorem 1.1, it is worth mentioning that the analogous classical problem has a very simple solution. In fact, using the technique of transfer matrices (see e.g. [23]), for any β the free energy per particle can be obtained from the eigenvalue of a simple d×d matrix and thus the problem reduces to standard numerical algorithms applied to some fixed matrix. This implies very efficient algorithms for any β including β = +∞.
However, the quantum case is significantly more complicated. This is illustrated for example by the fact that, when β = +∞, the problem is QMA-hard, and also that the simple Markov property for classical Gibbs states in one-dimension does not hold in the quantum setting. In his seminal work, Araki [24] proposed a quantum analogue of the transfer matrix, but it is a linear map between infinite-dimensional spaces. In our algorithm, we use a finite-dimensional approximation to this map. Our main technical result is to prove that the spectral radii of these finite-dimensional approximations converge superexponentially fast to e −βf β (h) . The algorithm (see Algorithm 2) is then simply to choose the finite-dimensional approximation parameter L for the transfer matrix as a function of the desired precision ε and then compute the spectral radius of the corresponding linear map.

Algorithm 2:
Algorithm for computing the free energy per particle. The constant C is a number that can be obtained from our proofs. See Section 3 for more precise definitions.
Parameters: Inverse temperature β, universal constant C Input: d local dimension, Hamiltonian term h ∈ C d 2 ×d 2 such that h ≤ 1, error ε Output:f β approximation to the free energy f β (h) 1 L ← log(1/ε) exp(C(β + 1))/ log(log(1/ε)); /* parameter for approximation */ /* matrix representation of a linear map from C d L−1 ×d L−1 to itself: */ 2 L * L (·) ← tr L e −βH [1,L] /2 e βH [2,L] /2 (1 ⊗ ·)e βH [2,L] /2 e −βH [1,L] /2 ; 3 r L ← spectral radius of L * L ; 4f β ← −(log r L )/β ; To analyse the algorithm we make extensive use of Araki's expansionals [24] to show that the marginal ρ L on the first L − 1 sites of the infinite Gibbs state, is an approximate eigenvector of the finite-dimensional map L * L , i.e., that L * L (ρ L ) − e −βf β (h) ρ L 1 decays superexponentially fast in L. By using variational expressions of the spectral radius of positive maps (so called Collatz-Wielandt formula), this allows us to show that the spectral radius of L * L is superexponentially close to e −βf β (h) . We note that standard perturbation bounds for eigenvalues of non-normal operators have a very bad dependence on dimension, and are thus not usable here, see e.g., [25,Chapter VIII]. To prove that the corresponding eigenvector of L * L is close to ρ L , we establish a quantitative primitivity condition for L * L , i.e., we prove that a sufficiently high power of L * L maps nonzero positive semidefinite operators to positive definite ones. Using tools from the Perron-Frobenius theory of positive operators -more precisely the Hilbert projective metric-, this allows us to show that ρ L is superexponentially close to the eigenvector of L * L associated to its spectral radius.
Let us also mention that the above techniques based on [24] are specific to one dimension. In particular, the results in there are related to the absence of thermal phase transitions, which do occur in higher dimensions so an extension of this approach to higher dimension is not possible. See also [26,27] for hardness results of the classical partition function on bounded degree graphs.

Numerical implementation
We also implement our algorithm and run it on a Hamiltonian for which the free energy function is known exactly. We observe very small errors (machine precision) already for moderate choices of L. We also observe that for this example the scaling of the error with inverse temperature is better than the worst-case estimates derived theoretically.
Organization The paper is structured as follows, we first give in Section 2 the mathematical definitions in which the problem is formulated and review some technical results needed in our proof. We then present in Section 3 the statements of our main results including error bounds for the computed free energy and Gibbs state as well as a bound on the runtime of the corresponding algorithm. Sections 4 to 6 contain the proofs of our main results. In Section 7 we present our numerical results.

Preliminaries
We start by introducing some notation for the description of quantum many-body systems and recap results from [24] needed for the proofs in this manuscript. . If I ⊂ J are two intervals in Z, we will often identify elements of A I as elements of A J by tensoring with the identity on J \ I. The adjoint map is the partial trace over the sites in J \ I, tr J\I : A J → A I .

Setup and basic notations
Let h ∈ C d 2 ×d 2 be the two-body Hamiltonian acting on two sites 3 , and let h i,i+1 ∈ A [i,i+1] the corresponding Hamiltonian when acting on sites {i, i + 1}. If a ≤ b, we let For an inverse temperature β > 0 and any integer N > 1, the partition function is defined as Z β,N (h) = tr exp −βH [1,N ] .
The infinite translation-invariant free energy per particle is defined as: This limit is known to exist, see e.g., [28,Theorem 15.5]. By appropriately rescaling h, one can assume without loss of generality that β = 1; indeed, it is immediate to check that f β (h) = (1/β)f 1 (βh). In the rest of the paper we will thus assume that β = 1, and omit the subscript β from the definition f β (h). We will also suppress the dependence on h in the notation, as it will be clear from the context, and just write To keep track of the complexity of our algorithm in the inverse temperature, we make the following convenient definitions: Definition 2.1. We say that a constant G ≥ 0 is bounded in Hamiltonian norm (BHN) if there exist some numerical constants C 1 , C 2 > 0 such that We say that a function ε(L) ≥ 0 is superexponentially decaying if there exist some numerical constants C 1 , C 2 , C 3 such that It is immediate to check that products of BHN -constants have property BHN and that the product of a BHN -constant and a superexponentially decaying function is superexponentially decaying.
In the classical or commuting case, the operators E L are independent of L and act nontrivially only on particles 1 and 2, and the linear map L * L plays the role of a transfer matrix. For the general quantum case, it is a highly nontrivial result of [24] that the operator E L is close to an operator with support only on the first sites, that it is bounded uniformly in L, and that it converges for L → ∞. The map L * L can then be thought of as a finite-dimensional approximation of an infinite-dimensional "noncommutative transfer matrix" introduced by Araki [24]. 4 The key quantity of interest to us will be the spectral radius of L * L , which we will denote by r L r L = max{|λ| : λ is an eigenvalue of L * L }.
The following lemma collects the results from [24] that we will need for the analysis of our algorithm. Note that, for convenience, we use the formulation of [29] which mostly considers exponentially decaying Hamiltonians instead of strictly local ones but still recovers the following results when restricting to local interactions. The dependency of the constants on h might not be immediately clear from the formulation in [29], so we discuss how to obtain the claimed dependencies in Appendix A.

Lemma 2.2 ([29, Proposition 4.2]).
There exists a BHN -constant G and a superexponentially decaying function ε(L) such that for any L ≥ 2 From the above it also follows immediately (by renaming constants) that

Gibbs states
We also introduce marginals of thermal states on finite systems and in the thermodynamic limit. For m > L ≥ 2, these are defined as In [24] (see also [29,Lemma 4.15]) it was shown that the limit lim m→∞ ρ L,m exists, which we denote by ρ L :

Description of the algorithm and overview of the analysis
In this section, we present our main convergence results and the resulting runtime of the proposed algorithm. Given a 2-body Hamiltonian h ∈ C d 2 ×d 2 , our first main result shows that the sequence − log r L (where r L is the spectral radius of the map L * L introduced in (5)) converges superexponentially fast to the free energy per site f = f (h) of the infinite chain, namely: In order to compute other properties of the Gibbs state such as expectation values of local observables or correlation functions, it is useful to also have a description of the Gibbs state. It turns out that the eigenvector corresponding to the spectral radius of L * L approximates the marginal ρ L of the Gibbs state defined in (8). This is the object of the next theorem (recall the definition of a BHN -constant in Def. 2.1): There exist BHN -constants G, G and a superexponentially decaying function ε(L) such that the following is true.
Note that the above theorem proves convergence to the one-sided Gibbs state, i.e., the Gibbs state of a chain that extends to infinity only in one direction. We focus on this case to simplify the presentation. The case of the chain that is infinite in both directions can be handled in the exact same way by using the two-sided version of the noncommutative transfer operator as defined in [24]. Alternatively, the one-sided version of the algorithm can be used in a black box way to obtain the two-sided marginal by recasting the Hamiltonian as we explain in Appendix D. We also note that one can compute the marginals of the Gibbs state by computing certain directional derivatives of the free energy function. This is described in more detail in Appendix B and it directly yields the marginals for the two-sided chain, as the free energy is the same in both cases. It is, however, inefficient to use this method to compute marginals of many particles.
Theorems 3.1 and 3.2 imply that the free energy and the k-particle marginals of the Gibbs state can be approximated up to any given error by choosing an appropriate L, and computing the spectral radius of L * L and the corresponding eigenvector. By making this choice and its dependence on the problem input explicit, we obtain the following corollary.
and outputs an approximationf such that Furthermore, for some numerical constants C 1 , . . . , C 5 > 0 and a given system size k, there exists an algorithm that takes the same inputs as above, runs in time at most and outputs an approximation v ∈ A [1,k] While we have made both the dependence on ε and the dependence on h explicit, we want to point out that our algorithm should be considered for some fixed bound on h (recall that we have assumed β = 1, for general β, this value is the inverse temperature times the norm of the Hamiltonian). As pointed out in the introduction, the ground state problem for the system we consider is computationally hard. If an algorithm was able to solve the free energy problem efficiently in 1/ε and h , it could also approximate the ground state energy efficiently by taking h → ∞. We formulate and prove such a reduction in Appendix C, i.e., we establish the QMA EXP -hardness of the infinite translation-invariant free energy problem with the temperature as an additional problem input. This shows that, unless QMA EXP = EXP, no algorithm can have a running time of the form exp(polylog( h , 1/ε)).
We note that the temperature dependence in the second part of the corollary is worse than in the first part 5 , which could be an artefact of our proof method (see Remark 5.6). In fact, it is also possible to compute thermodynamic quantities as derivatives of the free energy with respect to parameters in the Hamiltonian. This allows us to just use the first part of our proposed algorithm to access other properties of the Gibbs state. The details of this approach are worked out in detail in Appendix B for 2-local observables. In this case, this slightly improves the temperature dependence, as it removes the second term in (10), while still having the same dependence as the first term. The reason for this is that the error bounds for approximations of the derivative involve the second derivative of the free energy, which can only be bounded by invoking results from [24] as it is related to the analyticity of the free energy. We note however that this method is not practical to compute the average of k-local observables for moderately large k > 2, as it requires to block the system first which results in a sytem with local dimension d k/2 and interaction strength up to (k/2) h . On the other hand, our algorithm allows us to get L-local marginals for free. In addition, as illustrated in the numerical results Section 7, our algorithm also allows us to compute entropic quantities on the marginals such as the (conditional) mutual information.
4 Convergence speed of the spectral radius of L * L (proof of Theorem 3.1) In this section, we prove Theorem 3.1, showing that − log r L converges superexponentially fast in L to the free energy f , where r L is the spectral radius of the finite-dimensional linear map L * L defined in (5). The key to our proof is to show that ρ L , the marginal of the Gibbs state on [1, L − 1] (see (7)), is an approximate eigenvector of L * L , in the sense that where ε(L) is a superexponentially decaying function in L. This will be proven in Corollary 4.4. We start by giving an alternative expression for the free energy. While the limit of the normalized log-partition function exists, this does not immediately imply that the unnormalized partition function grows in approximately constant steps. It can, however, be shown in the case of a translation-invariant interaction, using locality results of Gibbs states: Lemma 4.1. The free energy defined in (3) also admits the following expression: Moreover Proof. Using the continuity of the logarithm, we have to prove the convergence of Thus for N, N > M we have Let ε > 0. By Lemma 2.2, there is a large enough M such that for any N, N > M the last two terms have magnitude ≤ ε. Furthermore, since the sequence ρ M +1,N converges as N → ∞ (see (8)), there is a large enough N 0 such that for any N, Together with the boundedness of E † M +1 E M +1 from Lemma 2.2, this tells us that the sequence (Z N +1 /Z N ) is Cauchy and thus converges.
The moreover part follows from The second technical lemma we need gives a bound on the operator norm of ρ L,m , the marginal of the finite Gibbs state defined in (7), and its inverse.
Proof. We consider the operator Taking the partial trace (which is a positive map), we get, with γ = tr e −H [L,m] > 0, and recalling thatρ L,m = tr [L,m] e −H [1,m] , We then divide by Z m which normalizesρ L,m Using translation invariance, the factor γ/Z m can be bounded as follows The next lemma is crucial, and shows that ρ L is an approximate eigenvector of the linear map L * L with eigenvalue e −f , for the trace distance.

Lemma 4.3.
There is a superexponentially decaying function ε(L) such that Proof. For any m > L, we have, with By Lemma 2.2 we have X − 1 ≤ ε(L) for some superexponentially decaying function ε(L) and X ≤ G 2 for some BHN -constant G. Let Z m = tr e −H [1,m] so that ρ L,m = ρ L,m /Z m . Let S be an arbitrary Hermitian operator in A [1,L−1] . Then, interpreting as necessary S ∈ A [1,m+1] by tensoring with the identity and denoting the Hilbert-Schmidt inner product as Note that (13)) and tr ρ L,m+1 = 1, the second term in the last line of (15) has a magnitude which is at most ε (L) S , where ε (L) is a superexponentially decaying function. Thus this shows that for any S, we have Letting m → ∞ gives us (14) as desired.
The following corollary shows that ρ L is an approximate eigenvector of L * L in the positive semidefinite order.

Corollary 4.4. There is a superexponentially decaying function ε(L) such that
Proof. From (14) we get, The above can be rewritten as Using Lemma 4.2 to bound ρ −1 L and e f ≤ e h , we obtain that ε(L) := e f ε 1 (L) ρ −1 L is still superexponentially decaying, and we have We can now finish the proof of Theorem 3.1.
Proof of Theorem 3.1. Since L L is a positive map with the same spectral radius as L * L , there exists w L ≥ 0 such that L L (w L ) = r L w L [31, Theorem 2.5]. Then, we have Since w L , ρ L > 0 (because w L ≥ 0 and ρ L > 0), we get from the above that r L ≤ e −f (1 + ε(L)). In a similar way, we get r L ≥ e −f (1 − ε(L)), and so as desired.

Remark 4.5.
What we have implicitly used in the proof of Theorem 3.1 above (Equation (17)) is the following well-known variational characterization of the spectral radius of a positive map T , known as the Collatz-Wielandt formula: It holds for any irreducible positive operator T and we have denoted r(T ) its spectral radius, see e.g., [32]. In our setting, we do not need to establish irreducibility as we are simply looking for a bound, and because the "candidate" eigenvector ρ L is positive definite. However to get that ρ L is close to the eigenvector of L * L associated to its spectral radius, we do establish irreducibility (and even the stronger property of primitivity) of L * L for sufficiently large L in the next section, see Lemma 5.5. In this section, we prove that the eigenvector of L * L associated to its spectral radius is superexponentially close in trace distance to ρ L , the marginal of the Gibbs state on [1, L−1]. The key to obtain this result is to prove that the linear maps L * L are primitive, i.e., that for a sufficiently large integer k, (L * L ) k = L * L • · · · • L * L maps any positive semidefinite operator, to a positive definite one. It is well-known in the Perron-Frobenius theory for positive operators [32], that primitivity implies that the spectral radius is a nondegenerate eigenvalue. The next theorem gives a quantitative version of this fact, and shows that a strengthening of the primitivity assumption, implies that approximate eigenvectors (in the sense of (11)) are close to the eigenvector associated to the spectral radius. Whereas such a result seems quite natural, we were not able to find it in the literature, so we include here a statement for general positive operators. Assume there exists an integer k such that the following is true: for any ψ ∈ E, T k (ψψ * ) = T • · · · • T (ψψ * ) > 0, and more precisely Then the eigenvector v ∈ B + (E) of T with tr(v) = 1 associated to its spectral radius is positive definite and unique. Furthermore, if x > 0 satisfies for some constant r > 0 and tr x = 1, then x − v 1 ≤ 2kCε.
In subsection 5.1, we prove Theorem 5.1 using as a main tool the Hilbert projective metric. In subsection 5.2 we show that the maps L * L satisfy the assumption (18) and use this to complete the proof of Theorem 3.2.

Proof of Theorem 5.1
We start by recalling the definition of the Hilbert projective metric on B + (E) \ {0}. We adopt notations from [33].
Given x, y ∈ B + (E), x, y = 0 we let Note that inf(x/y) = 1/ sup(y/x). The Hilbert metric is defined by It is easy to check that the following properties are satisfied: • Quasi-convexity (see e.g., [34, Lemma 6.2]): y) for any x, y > 0. For primitive maps, one can show that T is in fact a contraction for the Hilbert metric.
Then for any x, y ≥ 0, It is easy to check that if T is primitive and satisfies assumption (18), then ∆(T ) < ∞. Indeed, we can easily prove Proof. Since d H is quasi-convex in each of its arguments, one can restrict x and y in the definition of ∆(T ) to be rank-one. Furthermore, we have by the triangle inequality Note that d H (a, 1) = log(λ max (a)/λ min (a)), and so equation (18) says that d H (T x, 1) ≤ log(C) for any rank-one positive operator x. This proves the claim.
We are now ready to prove Theorem 5.1.

Proof of Theorem 5.1. By iterating (19) k times we get
This tells us that d H (T k x, x) ≤ 2k log(1 + ε). Applying Prop. 5.4 with T k we get that where in the last inequality we used Prop. 5.3, which asserts that ∆(T k ) ≤ 2 log C, which implies 1 − tanh(∆(T k )/4) ≥ 2/(C + 1) ≥ 1/C. Inequality (22) shows that x and v are close to each other in the Hilbert projective metric. It remains to prove a similar bound with the trace distance. We assume henceforth that tr x = tr v = 1.

Proof of Theorem 3.2
To apply Theorem 5.1 to our map L * L , we need to show that condition (18) holds. In particular, we need a lower bound on the lowest eigenvalue of the iterated channel applied to some rank-1 projector. We allow for a number of channel iterations that scales linearly with L and a lower bound that decays exponentially in L. There exist BHN -constants G, G , G , G , such that for L ≥ exp(G ) and for all ψ ∈ H [1,L−1] with ψ = 1 we have

Lemma 5.5.
Remark 5.6. We believe that Lemma 5.5 can be improved. In particular, it is natural to think that L * L is primitive for all L ≥ 2. In fact, we conjecture that for any L ≥ 2, and any ψ ∈ H [1,L−1] such that ψ = 1, (L * L ) L−1 (ψψ * ) ≥ e −O(L) 1. If this conjecture is true, then this will improve the dependence on h in the complexity estimate (10).
Our proof of Lemma 5.5 will require some additional results from [24] concerning the limit of the maps L * L and their adjoint, which we recall now.
Limit of L L Note that the adjoint of the map L * L is given by The key additional fact that we will require from [24] is that if L is the limit of the maps (L L ), then the iterations L n converge, as n → ∞ to a "rank-one" operator (property (iii) below). More precisely, we have 7 : We are now ready to prove Lemma 5.5. 6 Indeed, if −B ≤ A ≤ B then necessarily A 1 ≤ tr B. This follows by noting that if A = λ λP λ is a spectral decomposition of A, then − tr(BP λ ) ≤ tr(AP λ ) ≤ tr(BP λ ) and so A 1 = λ | tr(AP λ )| ≤ λ tr(BP λ ) = tr(B). 7 We use the notation A N for the closure, with respect to the operator norm, of the union of all local algebras A [1,n] , where operators in A [1,n] are embedded into A [1,m] for m > n by tensoring with the identity.
Proof of Lemma 5.5. Throughout the proof G i are BHN -constants. As a first step we note that we can equivalently prove the statement for the adjoint channel: Now using L and g from Lemma 5.7 and introducing an integer M that we choose later, where the choice of M will ensure that the term in brackets is positive. In the last step we applied Lemma 4.2 and Lemma 5.7 (iv) to lower bound the first term. The second term is due to Lemma 5.7 (iii). The last term is bounded by a superexponentially decaying function ε(L), which follows from the following telescope sum: We where C 3 is a numerical constant. We now impose that to bound the second term in brackets by 1 4 e −2 h L and similarly to ensure the same bound for the last term. We combine the above bounds into a new BHN -constant G such that L ≥ exp(G ) implies all the above inequalities. Under this condition we conclude with G = 2G 1 and G = 3 h G.

Proof of Theorem 3.2.
We consider again BHN -constants G i . We know from Lemma 5.7 that L * L (Q) ≤ G 1 Q 1, and so this implies that for the iterated map (L * L ) G 2 L (Q) ≤ G G 2 L 1 Q 1. Together with Lemma 5.5 this implies that for any pure state Q = ψψ * , Furthermore, from Corollary 4.4 we know that there is a superexponentially decaying function ε(L) such that Applying Theorem 5.1, we get that for for some BHN -constant G and the superexponentially decaying function ε(L).

Proof of Corollary 3.3
Proof of Corollary 3.3. We start with the algorithm for the free energy. The argument is completely analogous for the Gibbs state. Theorem 3.1 tells us that the error incurred by estimating the free energy via the spectral radius of L * L is at most where G is a BHN -constant and C is of the form C 3 h from the theorem. We need to find a value of L for which the error above is guaranteed to be less than the desired error ε.
As the expression in (24) cannot be inverted analytically to get L in terms of ε, we use Stirling's approximation to upper bound (24). Thus a sufficient condition on our L is that it satisfies: Note that we replaced C by another constant C = C 3 ( h + 1) assuming without loss of generality C 3 > 1 to account for the factor e in the Stirling approximation we used: n! > (n/e) n (see [23,Lemma B.3]).
Representing the concave function in the exponent on the right-hand side as a minimum over linear functions, one can show that this is implied if the condition is fulfilled for some γ. Using γ = log(log(1/ε))/2, one can check that the value L = log(1/ε)(2 + 2 exp(2C − 1)) + 2 log(G) log(log(1/ε)) will work. The algorithm is now given by the computation of the spectral radius of a square matrix of size d 2(L−1) × d 2(L−1) , which can be done in time polynomial in its size. Making the constants explicit again, using the assumption ε < 1/e to include the log(G) in the prefactor of log(1/ε), and absorbing all numbers into the new constant C 1 > 0 the resulting runtime then reads exp log(d) log(1/ε) exp (C 1 ( h + 1)) log(log(1/ε)) .
For the Gibbs state approximation we simply replace C by a BHN -constant. By the exact same calculation we obtain for some positive constants C 1 , In addition, we have to choose L ≥ exp(G) for some BHN -constant G for Theorem 3.2 to hold and furthermore L has to be at least as big as the size of the marginal that we want to obtain, which is another input to our algorithm.

Numerical Results
We implemented our algorithm and tested it on a dimerized XY model which is exactly solvable [35,36]. For N + 1 spin-1/2 particles, the Hamiltonian is given by where the spin operators are given by the Pauli matrices S x,y = σ x,y /2. In the case γ = 1, this model is translation-invariant and is given by the interaction term h = − β 4 (σ x ⊗ σ x + σ y ⊗σ y ). For γ = 1 we do not have translation-invariance in the above sense, but the model still becomes translation-invariant by blocking particles, i.e., it is a translation-invariant Hamiltonian with local dimension d = 4 and the interaction term Note also that the system is gapless for γ = 1 and gapped otherwise [36]. As discussed in the introduction, for gapped systems there exists an efficient algorithm to approximate the ground state energy, while the problem is QMA-hard in general. We show results for the error of the log-partition function βf β (h) and its dependence on temperature and the parameter L in Figure 1 for γ = 1. The value of the log-partition function itself as well as the expectation value of the energy is depicted in Figure 2 for γ = 2. We compute the energy for the two-sided infinite version of the system using the procedure described in Appendix D, but also obtain the same results by using a two-sided version of the map L * L . In Figure 4, we show the decay of the mutual information between two particles depending on their distance and the same for the conditional mutual information where the conditioning is on the rest of the marginal obtained from the computation (see [9] for an analytic result on this decay at high temperature). These were computed using the eigenvector of L * L . Note that they cannot be obtained efficiently from the derivative method (see Theorem 3.2) as the derivative method detailed in Appendix B is not efficient for obtaining marginals of many particles. While the theoretical bounds on L for a given error derived in this paper are impractically large, we observe very accurate results for moderate choices of the parameter. Also the error for the free energy seems to grow less than exponential with β for the chosen example as opposed to our worst-case estimate, which is doubly exponential. In Figure 3, we show the β-dependency of the error in the energy calculation. The dependency seems to be similar to the one for the free energy calculation. This is in contrast with the worse error dependence in the estimates we proved for the k-particle marginals (see Corollary 3.3). All computations were done on a laptop computer, where computation of a single datapoint for d = 2, L = 11 takes about 13 s.     Figure 4: Decay of the conditional mutual information (a) and the mutual information (b) for the two-sided chain and γ = 1 between two particles with distance d(·, ·). The conditioning is on the remaining particles in the 8-site marginal. We observe a growth of the decay rate with temperature.
A A note on the constants in [29] We restated multiple results from [29], which form our main technical tools in this work. While in the reference, constants are mostly left with an implicit dependence on h , we aim to make this explicit in our bounds. In this appendix we review some of these constants to show that they are bounded is given in the main text. Also note that the paper mostly considers exponentially decaying interactions and mostly just exponential instead of super-exponential decays. We can however obtain the better results from there by restricting to local Hamiltonians again. We list the following constants for a 2-local interaction Φ.
The following constant appears in several results cited in our main text and is a BHNconstant due to the bound in [29,Section 4.2]: To obtain the superexponential decay of G l we start by upper bounding Ω * k (2). To that end a generally useful formula is the Lagrange remainder bound for the tail of the exponential series: Note that the terms in the outer summation in the definition of Ω * k (2) are zero whenever n < k/2 as for those one of the factors has to be zero. Furthermore, the number of vectors in N n with entries no larger than 2 is bounded by 3 n and we obtain Plugging this into the expression for G l yields which gives us the superexponential-decay of this quantity. We finally need a bound on the constants K 2 and δ 2 from [29, Theorem 4.12] which relies on [29,Theorem 4.11]. In the later theorem we have log (2) and thereby N = 3r = 3 log(a) log (2) . K 2 is bounded by which is a BHN -constant if 1 + l≥1 G l 2 l is as well. This is the case because, looking at the l-dependent term in G l which is BHN and the coefficient in G l is BHN as well.
We now consider the constants in the proof of [29,Theorem 4.12] (note that there will be another constant K 2 distinct from the one above). We have K = K 2 = 4C 2 |L r | 1,2 (1 + |h| 1,2 ). C 2 is BHN as shown above. For the superoperator norm we apply [29,Corollary 4.9] to bound where the first sum is BHN due to the decay of the G l as before and the second sum is upper bounded by G 4 |Q| 1,2 ≤ G 4 making the whole expression BHN .
A bound for |h| 1,2 can be deduced from [29,Theorem 4.10], which implies membership of h in certain "quasilocal" sets. Namely, equation (42) implies h l ≤ 2G 3 G l , which we can sum to obtain which is now BHN as the sum converges as before. This closes the argument for K.
The constant δ is given as with C 2 and N = 3 log(2) log(a) = 3 log(2) log(C 2 ) We estimate using the concavity of the logarithm and its derivative so finally δ ≥ 1 6C 2 log(C 2 )/ log (3) and the denominator is bounded by a BHN -constant.

B Reduction of Local Observables to Free Energy
In this section we show how one can compute the expectation value of a 2-local observable P 1,2 ∈ A [1,2] in the thermodynamic limit, by computing the free energy of an appropriately perturbed Hamiltonian. More precisely, the quantity we are interested in is where we assume P ≤ 1. Our discussion is to some extent similar to the argument in [17,Lemma 11] and adapts it to the infinite translation-invariant case. The main idea is that a thermal expectation value can be written as a derivative of the free energy (or the partition function) with respect to a parameter introduced in the Hamiltonian. The derivative can then be approximated by a finite difference with an error that can be bounded by the second derivative. The case of infinite systems presents a particular challenge compared to the finite case, specifically because the free energy can be nonanalytic in the thermodynamic limit, which gives rise to phase transitions. One important consequence of Araki's result [24,Lemma 9.3] however, is that this does not happen in one-dimensional finite-range systems. We refine this result and provide a quantitative version, i.e., we establish a bound on the second derivative of the free energy.
This then allows us to quantitatively bound the error caused by a finite difference approximation of the derivative. Combining this approach with our Algorithm 2 to compute the free energy, we prove the following result.
Lemma B.1. For some numerical constants C 1 , C 2 > 0, there exists an algorithm that takes as input the local dimension of a quantum system d, its 2-local translation invariant Hamiltonian h, a 2-local normalized observable P , an additive error ε, runs in time at most and outputs an approximationμ such that |µ −μ| ≤ ε, where µ is defined in (26).
Note that we restrict to 2-local observables as our algorithm for the free energy (Algorithm 2) only applies to 2-local Hamiltonians supported on two neighbouring sites. The more general case of k-local observables (again supported on k contiguous sites) follows by first blocking the system. We note however that the resulting algorithm will have a bad dependence on k in its runtime, as the local dimension of the blocked Hamiltonian will be d k/2 and its interaction strength can be of the order of k h /2.

B.1 Preliminaries
For the proof of Lemma B.1, we will need some additional results from [29] and [14] concerning the quasi-locality of the time evolution operator, and decay of correlations. We introduce the necessary notations here.
We also need the following notation for the thermal state of the finite system on [a, b]: The results needed for our proof are summarized in the following Lemma.
The algorithm in the proof of Lemma B.1 is based on computing the derivative of the free energy function with respect to a parameter in a perturbed Hamiltonian. The perturbed translation-invariant Hamiltonian we consider is defined by We denote the corresponding partition function Z N (ε) = tr e −H [1,N ] (ε) and free energy per site f N (ε) = − 1 N log Z N (ε), as well as the limit We also make use of the perturbed thermal state associated to the finite system on the interval .
For convenience, we letρ ε,N =ρ ε, [1,N ] . The first lemma shows that the expectation value µ defined in (26) is equal to the derivative of the free energy for the perturbed Hamiltonian f (ε) defined in (27). Lemma B.3. For any 2-local observable P ∈ A [1,2] we have where f (ε) is the free energy per site of the perturbed Hamiltonian, as defined in (27).
Proof. We assume without loss of generality that P = 1. It is easy to verify that for any N , Define µ(ε) = lim N →∞ tr[P 1,2ρε,[−N,N ] ]. Lemma B.2 (i) tells us that for any i ∈ [1, N ], for some constants C, δ > 0. It thus follows that where in the last step M is any integer in [1, N/2]. By taking M = N/2 we see that To prove that d dε f N → d dε f , we need the convergence to be uniform on all compact sets. While the constants C and δ do depend on the interaction strength h + ε, an upper bound on these constants grows monotonically with h + ε. Therefore, we just choose the constants for interaction strength h + ε to prove uniform convergence for ε ∈ [0, ε ]. The convergence of the derivatives then follows from a standard result in analysis [37,Theorem 7.17].
Before we proceed we need a basic result from calculus. We were not able to find a proof in the literature for this specific setting so we include a proof for completeness. Note that the statement does not require convergence of the derivatives g N .
Lemma B.4. Let g N : (a, b) → R be a sequence of continuously differentiable functions converging to a continuously differentiable function g. If the derivatives of the elements are uniformly bounded g N (x) ≤ C, then the same bound holds for the limit g (x) ≤ C.
where the last line can be derived using the Duhamel formula (compare also [17]). We rewrite the second derivative in a slightly more compact form as where we used the time evolution operator with respect to the perturbed Hamiltonian We fix j ∈ [1, N − 1] and s ∈ [0, 1] and focus on the sum tr P j,j+1 Γ is,ε [1,N ] (P k,k+1 )ρ ε,N − tr[P j,j+1ρε,N ] tr[P k,k+1ρε,N ].
The idea is that the imaginary-time evolved operator Γ is,ε [1,N ] (P k,k+1 ) is approximately supported in a region close to {k, k + 1} and that it is thereby approximately uncorrelated from P j,j+1 if |j − k| is large.
We used the quasi-locality to bound the first and third term after the first inequality (note that tr[P k,k+1ρε,N ] = tr[Γ is,ε [1,N ] (P k,k+1 )ρ ε,N ] as the time evolution operator commutes with

C Reducing the Ground Energy Problem to the Free Energy
We reduce the ground-state problem to the free energy problem to establish QMA EXPhardness of the following problem.
Definition C.1. The Free Energy for Infinite Translation-Invariant Hamiltonian FE-ITIH is defined as follows Problem parameter: Three polynomialsp,q,r and d the dimension of a particle Problem input: N specified in binary, β ≥ 0 and the matrix h, each specified in binary with at most logr(N ) + 1 bits Promise: The free energy density of the infinite system defined by h, i.e., f β (h) = lim N →∞ − 1 βN log tr exp(−βH [1,N ] ) and lies outside the interval [1/p(N ), 1/p(N ) + 1/q(N )] Output: Determine if free energy density is at most 1/p(N ) or at least 1/p(N ) + 1/q(N ) Lemma C.2. There exist parametersp,q,r,d such that the above problem is QMA EXPhard.
Proof. We start by giving a quantitative convergence bound for the free energy to the ground-state energy for β → ∞.
Using the well-known variational formula for the free energy, we can write the free energy density for a finite system as The existence of the limit lim N →∞ f β,N shows the existence of the limit ground-state energy e 0 = lim N →∞ e 0,N . As a result, we have We now take a QMA EXP -hard problem (p, q, r, d) ITIH as defined in [3, Theorem 2.7]. We set the parameters for the FE-ITIH problem:p = p,q = 2 · q,r = r + 2(log d)q and we use the same dimension parameter d. Given an instance (N, h) for the ITIH problem, we will consider the instance (N, β, h) where we let β = 2q(N ) log d. Note that β ≤r(N ) and thus can be specified using at most logr(N ) + 1 bits. We now check that the promise is satisfied: we know that e 0 is either ≤ 1 p(N ) or at least 1 p(N ) + 1 q(N ) . Now for the free energy density, as f β ≤ e 0 , we have either f β ≤ 1 p(N ) or f β ≥ 1 p(N ) + 1 q(N ) − log d β = 1 p(N ) + 1 2q(N ) . It then follows that for the instance we constructed, e 0 ≤ 1 p(N ) if and only if f β ≤ 1 which proves the desired result.
We can interpret this hardness result by saying that we cannot expect to have an algorithm that has a dependence of the form exp(polylog(β, 1/ε)), unless QMA EXP = EXP.

D Recasting a two-sided chain to a one-sided chain
In this appendix, we explain how the Gibbs state of an infinite two-sided chain can be obtained from the Gibbs state of an infinite one-sided chain for a suitably modified Hamiltonian.
We consider a local dimension d and 2-particle Hamiltonian h ∈ B(C d ⊗ C d ) as given. We construct a system where the local Hilbert space corresponds to two copies of the original one, i.e., H = C d ⊗ C d with local dimension d = d 2 . We denote the sites of the system by indices with subscript i u , i d referring to either of the two copies for particle i respectively. The Hamiltonian for the one-sided chain h ∈ B(C d ⊗ C d ) is then chosen as h 1,2 = h 2u,1u + h 1 d ,2 d + h 1u,1 d − h 2u,2 d , see Figure 5 for an illustration.
In the sum over Hamiltonian terms on all sites, the third and last term occur each once for each particle and thereby cancel. The only exception is the first particle, where only the positive term appears. The construction can be thought of as a chain winding at the first site with both infinite ends going in the same direction. The thermal state on this construction is then equivalent to the two-sided infinite thermal state when using the order of sites as ρ iu,··· ,2u,1u,1 d ,2 d ,3 d ,··· ,i d .