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

Hamza Fawzi<fnref n="1" />, Omar Fawzi<fnref n="2" />, and Samuel O. Scalet<fnref n="1" />

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((1/\epsilon)^c)$ for any constant $c > 0$ where $\epsilon$ is the additive approximation error. Previously, the best known algorithm had a runtime that is polynomial in $1/\epsilon$. 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.

## 1 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

$$ \rho = e^{-\beta H} / Z_{\beta}(H) \qquad (1) $$

where $\beta = 1/T$ is the inverse temperature, and $Z_{\beta}(H) = \text{tr}[e^{-\beta H}]$ is the partition function. The free energy of the system at inverse temperature $\beta$ is defined as

$$ F_{\beta}(H) = -\frac{1}{\beta} \log Z_{\beta}(H). \qquad (2) $$

At zero temperature, i.e., $\beta = +\infty$, $F_{\beta}(H)$ becomes $\lambda_{\min}(H)$, the ground energy of $H$. The problem of computing the ground energy for a given local Hamiltonian is known to be

<footnote n="1">¹Department of Applied Mathematics and Theoretical Physics, University of Cambridge, United Kingdom</footnote>

<footnote n="2">²Univ Lyon, Inria, ENS Lyon, UCBL, LIP, Lyon, France</footnote>
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 = \sum_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$.<fnref n="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 \times d^2$. Our objective is to design an algorithm to approximate $f_β(h)$ with a good scaling in terms of the target error $\varepsilon$ 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 $\varepsilon$ 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 $\mathbb{C}^d \otimes \mathbb{C}^d$ satisfying $\|h\| \le 1$ and $\varepsilon \in (0, 1/e)$ and outputs an approximation $\tilde{f}_β$ satisfying $|\tilde{f}_β - f_β(h)| \le \varepsilon$, 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 \frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)}))$. Moreover, this algorithm can also compute an $\varepsilon$-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(\cdot)$ 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 \frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)}) \exp(O(\beta)))$, where $O(\cdot)$ only hides universal constants.<fnref n="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]} = \sum_{i=1}^{n-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

<footnote n="1">¹Technically, because of the choice of specification of input in [3], the problem is complete for a scaled version of QMA called QMA<sub>EXP</sub>.</footnote>

<footnote n="2">²The dependence is worse when computing marginals of the Gibbs state, see Corollary 3.3.</footnote>
be done in time polynomial in $d^n$ by explicitly writing the $d^n \times d^n$ matrix $H_{[1,n]}$. The sequence $f_{\beta,n}$ does converge to $f_{\beta}$ as $n \to \infty$, but the convergence is in general slow with an error decaying as $1/n$ due to the missing interaction term at the boundary. As a result, for a desired precision $\epsilon$, we obtain a runtime which is exponential in $1/\epsilon$. In order to obtain the subpolynomial dependence on $1/\epsilon$ in Theorem 1.1, we need to develop a more refined algorithm.

**Related work** In recent years, there have been multiple works about computing the free energy (or equivalently the partition function) at finite inverse temperature for a given Hamiltonian. In particular, for local Hamiltonians on an arbitrary bounded-degree graph, algorithms have been developed in [9, 10, 11] with performance guarantees when the inverse temperature $\beta$ is below some critical inverse temperature. The runtime of these algorithms is polynomial or quasi-polynomial in the number of particles and in $1/\epsilon$. These works rely on the so called cluster expansion, which, at its core, is a Taylor expansion of the partition function at $\beta = 0$. Truncating this expansion at a certain order and bounding the remainder terms allows for the approximation of the free energy for $\beta$ small enough. Indeed, the sum of remainder terms no longer converges if $\beta$ is too large which introduces a critical inverse temperature above which such algorithms do not have convergence guarantees.

However, a different method was used in [12] to obtain an algorithm for all temperatures for one-dimensional finite quantum systems. This algorithm combines several results from the analysis of 1D-systems: quantum belief propagation [13], together with a locality result about Gibbs states which was only recently proven in general [14]. For a system of $n$ particles in 1D, the running time of the algorithm is $n(1/\epsilon)^{O(1)}$, where $O(1)$ is a constant that depends exponentially on $\beta$. This algorithm can readily be applied to the infinite translation-invariant chain by setting $n = 1/\epsilon$ and this leads to an algorithm that is polynomial in $1/\epsilon$. For its implementation, the algorithm involves several choices of length scales to ensure convergence and numerical integrations to obtain the operators from the belief propagation.

Another algorithm, based on tensor networks, for computing thermal expectation values for a system of $n$ particles was designed in [15, 16] and has a runtime with an improved dependence on $\epsilon$, namely $n e^{\tilde{O}(\sqrt{\log(n/\epsilon)})}$, where the $\tilde{O}$ depends linearly on $\beta$. As with the algorithm of [12], it can be readily applied to the infinite translation-invariant chain by setting $n = 1/\epsilon$. This leads to an algorithm that is quasi-linear in $1/\epsilon$ 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/\epsilon$, thus leading to an algorithm for computing the free energy for infinite translation-invariant chains that is polynomial in $1/\epsilon$.

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/\epsilon$.

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 ($\beta = +\infty$) 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 $\beta$ the free energy per particle can be obtained from the eigenvalue of a simple $d \times d$ matrix and thus the problem reduces to standard numerical algorithms applied to some fixed matrix. This implies very efficient algorithms for any $\beta$ including $\beta = +\infty$.

However, the quantum case is significantly more complicated. This is illustrated for example by the fact that, when $\beta = +\infty$, 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^{-\beta f_{\beta}(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 $\epsilon$ 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 $\beta$, universal constant $C$

**Input:** $d$ local dimension, Hamiltonian term $h \in \mathbb{C}^{d^2 \times d^2}$ such that $\|h\| \le 1$, error $\epsilon$

**Output:** $\tilde{f}_{\beta}$ approximation to the free energy $f_{\beta}(h)$

1. $L \leftarrow \log(1/\epsilon) \exp(C(\beta+1))/\log(\log(1/\epsilon));$ /* parameter for approximation */

/* matrix representation of a linear map from $\mathbb{C}^{d^{L-1} \times d^{L-1}}$ to itself: */

2. $\mathcal{L}_L^*(\cdot) \leftarrow \text{tr}_L \left( e^{-\beta H_{[1,L]}/2} e^{\beta H_{[2,L]}/2} (\mathbf{1} \otimes \cdot) e^{\beta H_{[2,L]}/2} e^{-\beta H_{[1,L]}/2} \right)$;

3. $r_L \leftarrow$ spectral radius of $\mathcal{L}_L^*$;

4. $\tilde{f}_{\beta} \leftarrow -(\log r_L)/\beta$;

To analyse the algorithm we make extensive use of Araki's expansionals [24] to show that the marginal $\rho_L$ on the first $L-1$ sites of the infinite Gibbs state, is an approximate eigenvector of the finite-dimensional map $\mathcal{L}_L^*$, i.e., that $\|\mathcal{L}_L^*(\rho_L) - e^{-\beta f_{\beta}(h)}\rho_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 $\mathcal{L}_L^*$ is superexponentially close to $e^{-\beta f_{\beta}(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 $\mathcal{L}_L^*$ is close to $\rho_L$, we establish a quantitative *primitivity* condition for $\mathcal{L}_L^*$, i.e., we prove that a sufficiently high power of $\mathcal{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 $\rho_L$ is superexponentially close to the eigenvector of $\mathcal{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.

## 2 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.

### 2.1 Setup and basic notations

We consider a quantum 1D chain on the half-line $[1, \infty)$. Let $d \ge 2$ be the local dimension, $\mathcal{H}_i \cong \mathbb{C}^d$ the local Hilbert space, and $\mathcal{H}_{[a,b]} = \bigotimes_{i \in [a,b]} \mathcal{H}_i$ be the Hilbert space associated to sites $[a, b]$, where $b \ge a$ are integers. We let $\mathcal{A}_{[a,b]} = \mathcal{B}(\mathcal{H}_{[a,b]})$ be the associated algebra of operators acting on $\mathcal{H}_{[a,b]}$. If $I \subset J$ are two intervals in $\mathbb{Z}$, we will often identify elements of $\mathcal{A}_I$ as elements of $\mathcal{A}_J$ by tensoring with the identity on $J \setminus I$. The adjoint map is the partial trace over the sites in $J \setminus I$, $\text{tr}_{J \setminus I}: \mathcal{A}_J \to \mathcal{A}_I$.

Let $h \in \mathbb{C}^{d^2 \times d^2}$ be the two-body Hamiltonian acting on two sites<fnref n="3" />, and let $h_{i,i+1} \in \mathcal{A}_{[i,i+1]}$ the corresponding Hamiltonian when acting on sites $\{i, i+1\}$. If $a \le b$, we let

$$H_{[a,b]} = \sum_{i=a}^{b-1} h_{i,i+1} \in \mathcal{A}_{[a,b]}.$$

For an inverse temperature $\beta > 0$ and any integer $N > 1$, the partition function is defined as

$$Z_{\beta,N}(h) = \text{tr} \exp(-\beta H_{[1,N]}).$$

The infinite translation-invariant free energy per particle is defined as:

$$f_{\beta}(h) = \lim_{N \to \infty} -\frac{1}{\beta N} \log Z_{\beta,N}(h). \qquad (3)$$

<footnote n="3">³Note that without loss of generality we restrict to 2-local Hamiltonians for simplicity. A more general $r$-local Hamiltonian can be reduced by blocking $r/2$ sites into a single local Hilbert space of dimension $d^{(r/2)}$ and collecting the Hamiltonian terms. This recast Hamiltonian can then be used as an input to the algorithm.</footnote>
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 $\beta = 1$; indeed, it is immediate to check that $f_{\beta}(h) = (1/\beta)f_1(\beta h)$. In the rest of the paper we will thus assume that $\beta = 1$, and omit the subscript $\beta$ from the definition $f_{\beta}(h)$. We will also suppress the dependence on $h$ in the notation, as it will be clear from the context, and just write

$$Z_N = Z_{1,N}(h) \quad \text{and} \quad f = f_1(h).$$

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 $\mathcal{G} \ge 0$ is bounded in Hamiltonian norm (BHN) if there exist some numerical constants $C_1, C_2 > 0$ such that

$$\mathcal{G} \le \exp(C_1 \exp(C_2 \|h\|)).$$

We say that a function $\varepsilon(L) \ge 0$ is superexponentially decaying if there exist some numerical constants $C_1, C_2, C_3$ such that

$$\varepsilon(L) \le \exp(C_1 \exp(C_2 \|h\|)) \frac{\exp(C_3 \|h\| L)}{(L/2)!}.$$

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.

## 2.2 The map $\mathcal{L}_L^*$

For an integer $L \ge 2$, define

$$E_L = e^{-H_{[1,L]}/2} e^{H_{[2,L]}/2} \in \mathcal{A}_{[1,L]} \quad (4)$$

and let $\mathcal{L}_L^* : \mathcal{A}_{[1,L-1]} \to \mathcal{A}_{[1,L-1]}$ be the positive linear map defined by

$$\mathcal{L}_L^*(Q) = \text{tr}_L(E_L(\mathbf{1} \otimes Q)E_L^\dagger). \quad (5)$$

In the classical or commuting case, the operators $E_L$ are independent of $L$ and act non-trivially only on particles 1 and 2, and the linear map $\mathcal{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 \to \infty$. The map $\mathcal{L}_L^*$ can then be thought of as a finite-dimensional approximation of an infinite-dimensional “noncommutative transfer matrix” introduced by Araki [24].<fnref n="4" /> The key quantity of interest to us will be the spectral radius of $\mathcal{L}_L^*$, which we will denote by $r_L$

$$r_L = \max\{|\lambda| : \lambda \text{ is an eigenvalue of } \mathcal{L}_L^*\}. \quad (6)$$

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.

<footnote n="4">⁴By extending the maps $\mathcal{L}_L^*$ to an infinite-dimensional setting a limit of these maps can be defined. In [24, 29] it has been shown that its eigenvector and eigenvalue are a thermal state and the free energy of the system respectively, regaining the meaning of transfer matrices as in the classical case.</footnote>
**Lemma 2.2** ([29, Proposition 4.2]). There exists a BHN-constant $\mathcal{G}$ and a superexponentially decaying function $\varepsilon(L)$ such that for any $L \ge 2$

$$ (i) \quad \|E_L\|, \|E_L^{-1}\| \le \mathcal{G} $$

$$ (ii) \text{ for } L \le M \text{ we have } \|E_L - E_M\|, \|E_L^{-1} - E_M^{-1}\| \le \varepsilon(L) $$

From the above it also follows immediately (by renaming constants) that

$$ (iii) \text{ for } L \le M \text{ we have } \|E_L E_M^{-1} - \mathbf{1}\| \le \varepsilon(L) $$

## 2.3 Gibbs states

We also introduce marginals of thermal states on finite systems and in the thermodynamic limit. For $m > L \ge 2$, these are defined as

$$ \begin{aligned} \bar{\rho}_{L,m} &= \text{tr}_{[L,m]} e^{-H_{[1,m]}} \in \mathcal{A}_{[1,L-1]} \\ \rho_{L,m} &= \bar{\rho}_{L,m} / \text{tr} \bar{\rho}_{L,m}. \end{aligned} \quad (7) $$

In [24] (see also [29, Lemma 4.15]) it was shown that the limit $\lim_{m \to \infty} \rho_{L,m}$ exists, which we denote by $\rho_L$:

$$ \rho_L = \lim_{m \to \infty} \rho_{L,m}. \quad (8) $$

# 3 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 \in \mathbb{C}^{d^2 \times d^2}$, our first main result shows that the sequence $-\log r_L$ (where $r_L$ is the spectral radius of the map $\mathcal{L}_L^*$ introduced in (5)) converges superexponentially fast to the free energy per site $f = f(h)$ of the infinite chain, namely:

**Theorem 3.1.** There is a superexponentially decaying function $\varepsilon(L)$ such that $|(-\log r_L) - f| \le \varepsilon(L)$ for all $L \ge 2$.

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 $\mathcal{L}_L^*$ approximates the marginal $\rho_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):

**Theorem 3.2.** There exist BHN-constants $\mathcal{G}, \mathcal{G}'$ and a superexponentially decaying function $\varepsilon(L)$ such that the following is true. Let $v_L \in \mathcal{A}_{[1,L-1]}$ be the eigenvector of the map $\mathcal{L}_L^*$ corresponding to its spectral radius. For $L \ge e^\mathcal{G}$, $v_L$ is unique, strictly positive, and

$$ \|v_L - \rho_L\|_1 \le e^{\mathcal{G}'\varepsilon(L)}. $$

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* and the corresponding eigenvector. By making this choice and its dependence on the problem input explicit, we obtain the following corollary.

**Corollary 3.3.** For some numerical constant C > 0, there exists an algorithm that takes as input the local dimension of a quantum system d, its 2-local translation-invariant Hamiltonian given by h, an additive error ε < 1/e, runs in time at most

$$ \exp \left( \log(d) \exp \left( C \left( \|h\| + 1 \right) \right) \frac{\log(1/\varepsilon)}{\log(\log(1/\varepsilon))} \right), \quad (9) $$

and outputs an approximation $\tilde{f}$ such that $|f(h) - \tilde{f}| \le \varepsilon$, where $f(h)$ is defined in (3).

Furthermore, for some numerical constants $C_1, \dots, 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

$$ \exp \left( \log(d) \max \left\{ \exp \left( C_1 e^{C_2 \|h\|} \right) \frac{\log(1/\varepsilon)}{\log(\log(1/\varepsilon))}, e^{C_3 e^{C_4 \|h\|}}, C_5 k \right\} \right), \quad (10) $$

and outputs an approximation $v \in A_{[1,k]}$ such that $\|v - \rho_{k+1}\|_1 \le \varepsilon$, where $\rho_{k+1}$ is defined in (8).

While we have made both the dependence on $\varepsilon$ 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 $\beta = 1$, for general $\beta$, 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/\varepsilon$ and $\|h\|$, it could also approximate the ground state energy efficiently by taking $\|h\| \to \infty$. 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/\varepsilon))$.

We note that the temperature dependence in the second part of the corollary is worse than in the first part<fnref n="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

<footnote n="5">⁵To give more insight into the complexity estimate (10), the second term in the maximum in (10) comes from the lower bound on L in Theorem 3.2. The last term comes from the fact that we have to compute the eigenvector of at least the map $L_{k+1}^*$ for it to be defined on at least k sites (for larger maps we take a partial trace in the end).</footnote>
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 $\mathcal{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 $\mathcal{L}_L^*$ defined in (5). The key to our proof is to show that $\rho_L$, the marginal of the Gibbs state on $[1, L-1]$ (see (7)), is an approximate eigenvector of $\mathcal{L}_L^*$, in the sense that

$$ (1 - \varepsilon(L))e^{-f} \rho_L \le \mathcal{L}_L^*(\rho_L) \le (1 + \varepsilon(L))e^{-f} \rho_L \quad (11) $$

where $\varepsilon(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:

$$ f = - \lim_{N \to \infty} \log \frac{Z_{N+1}}{Z_N}. \quad (12) $$

Moreover $e^{-f}, e^f \le e^{\|h\|}$.

*Proof.* Using the continuity of the logarithm, we have to prove the convergence of $Z_{N+1}/Z_N$. For $N > M$, we have

$$ \begin{aligned} \frac{Z_{N+1}}{Z_N} &= \frac{\operatorname{tr}[e^{-H_{[1,N+1]}}]}{\operatorname{tr}[e^{-H_{[1,N]}}]} = \frac{\operatorname{tr}[\mathbf{1} \otimes e^{-H_{[1,N]}} E_{N+1}^{\dagger} E_{N+1}]}{\operatorname{tr}[e^{-H_{[1,N]}}]} \\ &= \frac{\operatorname{tr}[e^{-H_{[2,N+1]}} E_{M+1}^{\dagger} E_{M+1}]}{\operatorname{tr}[e^{-H_{[2,N+1]}}]} - \frac{\operatorname{tr}[\mathbf{1} \otimes e^{-H_{[1,N]}} (E_{M+1}^{\dagger} E_{M+1} - E_{N+1}^{\dagger} E_{N+1})]}{\operatorname{tr}[e^{-H_{[1,N]}}]} \\ &= \operatorname{tr}[\rho_{M+1} \operatorname{tr}_1 [E_{M+1}^{\dagger} E_{M+1}]] + \operatorname{tr}[(\rho_{M+1,N} - \rho_{M+1}) E_{M+1}^{\dagger} E_{M+1}] \\ &\quad - \operatorname{tr}[\mathbf{1} \otimes \rho_{N+1,N} (E_{M+1}^{\dagger} E_{M+1} - E_{N+1}^{\dagger} E_{N+1})]. \end{aligned} $$

Thus for $N, N' > M$ we have

$$ \begin{aligned} \frac{Z_{N+1}}{Z_N} - \frac{Z_{N'+1}}{Z_{N'}} &= \operatorname{tr}[(\rho_{M+1,N} - \rho_{M+1,N'}) E_{M+1}^{\dagger} E_{M+1}] \\ &\quad - \operatorname{tr}[\mathbf{1} \otimes \rho_{N+1,N'} (E_{M+1}^{\dagger} E_{M+1} - E_{N+1}^{\dagger} E_{N+1})] \\ &\quad + \operatorname{tr}[\mathbf{1} \otimes \rho_{N'+1,N'} (E_{M+1}^{\dagger} E_{M+1} - E_{N'+1}^{\dagger} E_{N'+1})]. \end{aligned} $$

Let $\varepsilon > 0$. By Lemma 2.2, there is a large enough M such that for any $N, N' > M$ the last two terms have magnitude $\le \varepsilon$. Furthermore, since the sequence $\rho_{M+1,N}$ converges
as $N \to \infty$ (see (8)), there is a large enough $N_0$ such that for any $N, N' \ge N_0 \|\rho_{M+1,N} - \rho_{M+1,N'}\|_1 \le \varepsilon$. Together with the boundedness of $E_{M+1}^\dagger 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

$$ \left| \log \left( \frac{Z_{m+1}}{Z_m} \right) \right| \le \|h\| $$

by [30, Lemma 3.6].

<figure><img src="image_1.png" /></figure>

The second technical lemma we need gives a bound on the operator norm of $\rho_{L,m}$, the marginal of the finite Gibbs state defined in (7), and its inverse.

**Lemma 4.2.** There is a BHN-constant $\mathcal{G}$ such that for any $L, m$, $\|\rho_{L,m}\|, \|\rho_{L,m}^{-1}\| \le \mathcal{G}e^{2\|h\|L}$.

*Proof.* We consider the operator

$$ \tilde{E}_{L,m} = e^{-H_{[1,m]}/2} e^{H_{[L,m]}/2} e^{H_{[1,L-1]}/2}. $$

By [29, Proposition 4.2], we know that $\mathcal{G}^{-1}\mathbf{1} \le \tilde{E}_{L,m}^\dagger \tilde{E}_{L,m} \le \mathcal{G}\mathbf{1}$ where $\mathcal{G} > 0$ is some BHN-constant. This implies that

$$ \mathcal{G}^{-1} e^{-H_{[1,L-1]}} \le e^{H_{[L,m]}/2} e^{-H_{[1,m]}} e^{H_{[L,m]}/2} \le \mathcal{G} e^{-H_{[1,L-1]}} $$

and, since $e^{-\|h\|L}\mathbf{1} \le e^{-H_{[1,L-1]}} \le e^{\|h\|L}\mathbf{1}$, we then get

$$ \mathcal{G}^{-1} e^{-\|h\|L} e^{-H_{[L,m]}} \le e^{-H_{[1,m]}} \le \mathcal{G} e^{\|h\|L} e^{-H_{[L,m]}}. $$

Taking the partial trace (which is a positive map), we get, with $\gamma = \text{tr} e^{-H_{[L,m]}} > 0$, and recalling that $\bar{\rho}_{L,m} = \text{tr}_{[L,m]} e^{-H_{[1,m]}}$,

$$ \mathcal{G}^{-1} e^{-\|h\|L} \gamma \mathbf{1} \le \bar{\rho}_{L,m} \le \mathcal{G} e^{\|h\|L} \gamma \mathbf{1}. $$

We then divide by $Z_m$ which normalizes $\bar{\rho}_{L,m}$

$$ \mathcal{G}^{-1} e^{-\|h\|L} \frac{\gamma}{Z_m} \mathbf{1} \le \frac{\bar{\rho}_{L,m}}{Z_m} \le \mathcal{G} e^{\|h\|L} \frac{\gamma}{Z_m} \mathbf{1}. $$

Using translation invariance, the factor $\gamma/Z_m$ can be bounded as follows

$$ \left| \log \left( \frac{\gamma}{Z_m} \right) \right| = \left| \log(Z_{m-L+1}) - \log(Z_m) \right| \le \|H_{[1,L-1]}\| \le \|h\|(L-1) \quad (13) $$

by [30, Lemma 3.6], so since $\rho_{L,m} = \bar{\rho}_{L,m}/Z_m$

$$ \mathcal{G}^{-1} e^{-2\|h\|L} \mathbf{1} \le \rho_{L,m} \le \mathcal{G} e^{2\|h\|L} \mathbf{1}. $$

or equivalently $\|\rho_{L,m}\|, \|\rho_{L,m}^{-1}\| \le \mathcal{G}e^{2\|h\|L}$ as desired.

<figure><img src="image_2.png" /></figure>

The next lemma is crucial, and shows that $\rho_L$ is an approximate eigenvector of the linear map $\mathcal{L}_L^*$ with eigenvalue $e^{-f}$, for the trace distance.

**Lemma 4.3.** There is a superexponentially decaying function $\varepsilon(L)$ such that

$$ \left\| \mathcal{L}_L^*(\rho_L) - e^{-f} \rho_L \right\|_1 \le \varepsilon(L). \quad (14) $$
Proof. For any $m > L$, we have, with $X = E_L E_{m+1}^{-1} \in \mathcal{A}_{[1,m+1]}$

$$
\begin{align*}
\mathcal{L}_L^*(\bar{\rho}_{L,m}) &= \mathrm{tr}_L(E_L \mathrm{tr}_{[L+1,m+1]} e^{-H_{[2,m+1]}} E_L^\dagger) \\
&= \mathrm{tr}_{[L,m+1]}(E_L e^{-H_{[2,m+1]}} E_L^\dagger) \\
&= \mathrm{tr}_{[L,m+1]}(X E_{m+1} e^{-H_{[2,m+1]}} E_{m+1}^\dagger X^\dagger) \\
&= \mathrm{tr}_{[L,m+1]}(X e^{-H_{[1,m+1]}} X^\dagger).
\end{align*}
$$

By Lemma 2.2 we have $\|X - \mathbf{1}\| \le \varepsilon(L)$ for some superexponentially decaying function $\varepsilon(L)$ and $\|X\| \le G^2$ for some BHN-constant $G$. Let $Z_m = \mathrm{tr} e^{-H_{[1,m]}}$ so that $\rho_{L,m} = \bar{\rho}_{L,m}/Z_m$. Let $S$ be an arbitrary Hermitian operator in $\mathcal{A}_{[1,L-1]}$. Then, interpreting as necessary $S \in \mathcal{A}_{[1,m+1]}$ by tensoring with the identity and denoting the Hilbert-Schmidt inner product as $\langle A, B \rangle = \mathrm{tr}[A^*B]$, we have

$$
\begin{align}
\langle \mathcal{L}_L^*(\rho_{L,m}), S \rangle &= \frac{1}{Z_m} \langle e^{-H_{[1,m+1]}}, X^\dagger S X \rangle \nonumber \\
&= \frac{Z_{m+1}}{Z_m} \left\langle \frac{e^{-H_{[1,m+1]}}}{Z_{m+1}}, X^\dagger S X \right\rangle \tag{15} \\
&= \frac{Z_{m+1}}{Z_m} \left\langle \rho_{L,m+1}, S \right\rangle + \frac{Z_{m+1}}{Z_m} \left\langle \frac{e^{-H_{[1,m+1]}}}{Z_{m+1}}, X^\dagger S X - S \right\rangle. \nonumber
\end{align}
$$

Note that

$$
\begin{align*}
\| X^\dagger S X - S \| &= \|(X^\dagger - \mathbf{1})S X + S(X - \mathbf{1})\| \\
&\le \|X - \mathbf{1}\| \|S\| (\|X\| + 1) \le \varepsilon'(L) \|S\|
\end{align*}
$$

where $\varepsilon'(L) = \varepsilon(L)(1 + G^2)$ is a superexponentially decaying function. Since $Z_{m+1}/Z_m$ is bounded by $\exp(\|h\|)$ (see Equation (13)) and $\mathrm{tr}\,\rho_{L,m+1} = 1$, the second term in the last line of (15) has a magnitude which is at most $\varepsilon''(L)\|S\|$, where $\varepsilon''(L)$ is a superexponentially decaying function. Thus this shows that for any $S$, we have

$$
\left| \left\langle \mathcal{L}_L^*(\rho_{L,m}) - \frac{Z_{m+1}}{Z_m} \rho_{L,m+1}, S \right\rangle \right| \le \varepsilon''(L) \|S\|
$$

which is the same as

$$
\| \mathcal{L}_L^*(\rho_{L,m}) - \frac{Z_{m+1}}{Z_m} \rho_{L,m+1} \|_1 \le \varepsilon''(L).
$$

Letting $m \to \infty$ gives us (14) as desired. $\square$

<figure><img src="image_3.png" /></figure>

The following corollary shows that $\rho_L$ is an approximate eigenvector of $\mathcal{L}_L^*$ in the positive semidefinite order.

**Corollary 4.4.** There is a superexponentially decaying function $\varepsilon(L)$ such that

$$
(1 - \varepsilon(L))e^{-f} \rho_L \le \mathcal{L}_L^*(\rho_L) \le (1 + \varepsilon(L))e^{-f} \rho_L.
$$

Proof. From (14) we get,

$$
\| \mathcal{L}_L^*(\rho_L) - e^{-f} \rho_L \| \le \| \mathcal{L}_L^*(\rho_L) - e^{-f} \rho_L \|_1 \le \varepsilon_1(L)
$$

The above can be rewritten as

$$
-\varepsilon_1(L)\mathbf{1} \le \mathcal{L}_L^*(\rho_L) - e^{-f}\rho_L \le \varepsilon_1(L)\mathbf{1}.
$$
Writing $1 \le \|\rho_L^{-1}\|\rho_L$ the above implies

$$-\varepsilon_1(L)\|\rho_L^{-1}\|\rho_L \le \mathcal{L}_L^*(\rho_L) - e^{-f}\rho_L \le \varepsilon_1(L)\|\rho_L^{-1}\|\rho_L.$$

Using Lemma 4.2 to bound $\|\rho_L^{-1}\|$ and $e^f \le e^{\|h\|}$, we obtain that $\varepsilon(L) := e^f\varepsilon_1(L)\|\rho_L^{-1}\|$ is still superexponentially decaying, and we have

$$e^{-f}(1 - \varepsilon(L))\rho_L \le \mathcal{L}_L^*(\rho_L) \le e^{-f}(1 + \varepsilon(L))\rho_L. \quad (16)$$

□

We can now finish the proof of Theorem 3.1.

*Proof of Theorem 3.1.* Since $\mathcal{L}_L$ is a positive map with the same spectral radius as $\mathcal{L}_L^*$, there exists $w_L \ge 0$ such that $\mathcal{L}_L(w_L) = r_L w_L$ [31, Theorem 2.5]. Then, we have

$$r_L \langle w_L, \rho_L \rangle = \langle \mathcal{L}_L(w_L), \rho_L \rangle = \langle w_L, \mathcal{L}_L^*(\rho_L) \rangle \le e^{-f}(1 + \varepsilon(L)) \langle w_L, \rho_L \rangle. \quad (17)$$

Since $\langle w_L, \rho_L \rangle > 0$ (because $w_L \ge 0$ and $\rho_L > 0$), we get from the above that $r_L \le e^{-f}(1 + \varepsilon(L))$. In a similar way, we get $r_L \ge e^{-f}(1 - \varepsilon(L))$, and so

$$e^{-f}(1 - \varepsilon(L)) \le r_L \le e^{-f}(1 + \varepsilon(L)),$$

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 $\mathcal{T}$, known as the Collatz-Wielandt formula:

$$\begin{aligned} r(\mathcal{T}) &= \sup_{w>0} \sup\{\lambda \ge 0 : \mathcal{T}(w) \ge \lambda w\} \\ &= \inf_{w>0} \inf\{\lambda \ge 0 : \mathcal{T}(w) \le \lambda w\}. \end{aligned}$$

It holds for any irreducible positive operator $\mathcal{T}$ and we have denoted $r(\mathcal{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 $\rho_L$ is positive definite. However to get that $\rho_L$ is close to the eigenvector of $\mathcal{L}_L^*$ associated to its spectral radius, we do establish irreducibility (and even the stronger property of primitivity) of $\mathcal{L}_L^*$ for sufficiently large $L$ in the next section, see Lemma 5.5.

# 5 Convergence speed for the eigenvector of $\mathcal{L}_L^*$ (proof of Theorem 3.2)

In this section, we prove that the eigenvector of $\mathcal{L}_L^*$ associated to its spectral radius is superexponentially close in trace distance to $\rho_L$, the marginal of the Gibbs state on $[1, L-1]$. The key to obtain this result is to prove that the linear maps $\mathcal{L}_L^*$ are primitive, i.e., that for a sufficiently large integer $k$, $(\mathcal{L}_L^*)^k = \mathcal{L}_L^* \circ \cdots \circ \mathcal{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.
**Theorem 5.1.** Let $\mathcal{E} \cong \mathbb{C}^n$ be a finite-dimensional Hilbert space, and let $\mathcal{T}: \mathcal{B}(\mathcal{E}) \to \mathcal{B}(\mathcal{E})$ be a Hermitian-preserving positive linear map, i.e., $\mathcal{T}(\mathcal{B}_+(\mathcal{E})) \subset \mathcal{B}_+(\mathcal{E})$, where $\mathcal{B}_+(\mathcal{E})$ is the cone of Hermitian positive semidefinite operators on $\mathcal{E}$.

Assume there exists an integer $k$ such that the following is true: for any $\psi \in \mathcal{E}$, $\mathcal{T}^k(\psi\psi^*) = \mathcal{T} \circ \dots \circ \mathcal{T}(\psi\psi^*) > 0$, and more precisely

$$ \frac{\lambda_{\max}(\mathcal{T}^k(\psi\psi^*))}{\lambda_{\min}(\mathcal{T}^k(\psi\psi^*))} \le C \quad \forall \psi \in \mathcal{E}. \qquad (18) $$

Then the eigenvector $v \in \mathcal{B}_+(\mathcal{E})$ of $\mathcal{T}$ with $\mathrm{tr}(v) = 1$ associated to its spectral radius is positive definite and unique. Furthermore, if $x > 0$ satisfies

$$ (1 + \varepsilon)^{-1}rx \le \mathcal{T}x \le (1 + \varepsilon)rx \qquad (19) $$

for some constant $r > 0$ and $\mathrm{tr} x = 1$, then $\|x - v\|_1 \le 2kC\varepsilon$.

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 $\mathcal{L}_L^*$ satisfy the assumption (18) and use this to complete the proof of Theorem 3.2.

## 5.1 Proof of Theorem 5.1

We start by recalling the definition of the Hilbert projective metric on $\mathcal{B}_+(\mathcal{E}) \setminus \{0\}$. We adopt notations from [33].

Given $x, y \in \mathcal{B}_+(\mathcal{E})$, $x, y \ne 0$ we let

$$ \sup(x/y) = \inf\{\lambda \ge 0 : x \le \lambda y\} $$

$$ \inf(x/y) = \sup\{\lambda \ge 0 : \lambda y \le x\} $$

Note that $\inf(x/y) = 1/\sup(y/x)$. The Hilbert metric is defined by

$$ d_H(x, y) = \log \left( \frac{\sup(x/y)}{\inf(x/y)} \right). $$

It is easy to check that the following properties are satisfied:

* $d_H(x, y) = d_H(y, x)$

* $d_H(\lambda x, \mu y) = d_H(x, y)$ for any $\lambda, \mu > 0$

* Triangle inequality: $d_H(x, z) \le d_H(x, y) + d_H(y, z)$

* Quasi-convexity (see e.g., [34, Lemma 6.2]): $d_H(x + y, z) \le \max(d_H(x, z), d_H(y, z))$

If $\mathcal{T}$ is a positive map, it is immediate to check that $\sup(\mathcal{T}x/\mathcal{T}y) \le \sup(x/y)$ and $\inf(\mathcal{T}x/\mathcal{T}y) \ge \inf(x/y)$, which implies $d_H(\mathcal{T}x, \mathcal{T}y) \le d_H(x, y)$ for any $x, y > 0$. For primitive maps, one can show that $\mathcal{T}$ is in fact a contraction for the Hilbert metric.

**Theorem 5.2** (Birkhoff, see e.g., [33, Theorem 4]). Let $\mathcal{T}: \mathcal{B}(\mathcal{E}) \to \mathcal{B}(\mathcal{E})$ be a positive map and let

$$ \Delta(\mathcal{T}) = \sup_{\substack{x,y \ge 0 \\ x, y \ne 0}} d_H(\mathcal{T}x, \mathcal{T}y). \qquad (20) $$

Then for any $x, y \ge 0$,

$$ d_H(\mathcal{T}x, \mathcal{T}y) \le \tanh(\Delta(\mathcal{T})/4)d_H(x, y). \qquad (21) $$
It is easy to check that if $\mathcal{T}$ is primitive and satisfies assumption (18), then $\Delta(\mathcal{T}) < \infty$. Indeed, we can easily prove

**Proposition 5.3.** If $\mathcal{T} : \mathcal{B}(\mathcal{E}) \to \mathcal{B}(\mathcal{E})$ is a positive map satisfying (18), then $\Delta(\mathcal{T}) \le 2\log(C)$.

*Proof.* Since $d_H$ is quasi-convex in each of its arguments, one can restrict $x$ and $y$ in the definition of $\Delta(\mathcal{T})$ to be rank-one. Furthermore, we have by the triangle inequality

$$d_H(\mathcal{T}x, \mathcal{T}y) \le d_H(\mathcal{T}x, \mathbf{1}) + d_H(\mathbf{1}, \mathcal{T}y).$$

Note that $d_H(a, \mathbf{1}) = \log(\lambda_{\max}(a)/\lambda_{\min}(a))$, and so equation (18) says that $d_H(\mathcal{T}x, \mathbf{1}) \le \log(C)$ for any rank-one positive operator $x$. This proves the claim. $\square$

Since the Hilbert metric is projective, we know that $v$ is a positive eigenvector of $\mathcal{T}$ if, and only if, $d_H(\mathcal{T}v, v) = 0$. If $\Delta(\mathcal{T}) < \infty$, then it is easy to see that $\mathcal{T}$ can only have one positive eigenvector, for if $v, w$ were two positive eigenvectors of $\mathcal{T}$, then $d_H(\mathcal{T}v, \mathcal{T}w) = d_H(v, w)$ which would contradict (21). One can make this statement more quantitative: if $d_H(\mathcal{T}x, x) \le \delta$, then $x$ is necessarily close to the unique positive eigenvector of $\mathcal{T}$ in the Hilbert metric.

**Proposition 5.4.** Let $\mathcal{T}$ be a positive map with $\Delta = \Delta(\mathcal{T}) < \infty$, and let $v$ be its (unique) positive eigenvector. Assume that $x > 0$ satisfies $d_H(\mathcal{T}x, x) \le \delta$. Then $d_H(x, v) \le \delta/(1 - \tanh(\Delta/4))$.

*Proof.* This is a simple triangle inequality. We have

$$d_H(x, v) \le d_H(x, \mathcal{T}x) + d_H(\mathcal{T}x, v) \\ = d_H(x, \mathcal{T}x) + d_H(\mathcal{T}x, \mathcal{T}v) \le d_H(x, \mathcal{T}x) + \tanh(\Delta/4)d_H(x, v)$$

and so $d_H(x, v) \le \delta/(1 - \tanh(\Delta/4))$ as desired. $\square$

We are now ready to prove Theorem 5.1.

*Proof of Theorem 5.1.* By iterating (19) $k$ times we get

$$r^k(1 + \varepsilon)^{-k}x \le \mathcal{T}^k x \le r^k(1 + \varepsilon)^k x.$$

This tells us that $d_H(\mathcal{T}^k x, x) \le 2k \log(1 + \varepsilon)$. Applying Prop. 5.4 with $\mathcal{T}^k$ we get that

$$d_H(x, v) \le \frac{2k \log(1 + \varepsilon)}{1 - \tanh(\Delta(\mathcal{T}^k)/4)} \le 2kC \log(1 + \varepsilon). \quad (22)$$

where in the last inequality we used Prop. 5.3, which asserts that $\Delta(\mathcal{T}^k) \le 2\log C$, which implies $1 - \tanh(\Delta(\mathcal{T}^k)/4) \ge 2/(C+1) \ge 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 $\mathrm{tr} x = \mathrm{tr} v = 1$.

Let $\eta = 2kC \log(1 + \varepsilon)$. Equation (22) tells us that

$$ce^{-\eta/2}v \le x \le ce^{\eta/2}v$$

for some constant $c > 0$. By taking the trace, and using the fact that $\mathrm{tr} x = \mathrm{tr} v = 1$, we see that $e^{-\eta/2} \le c \le e^{\eta/2}$ and so the inequality above implies

$$e^{-\eta}v \le x \le e^{\eta}v.$$
This implies that $(e^{-\eta} - 1)v \le x - v \le (e^{\eta} - 1)v$. Since $\eta = 2kC \log(1 + \epsilon)$ and $2kC \ge 1$, it is easy to check that $e^{\eta} - 1 \le 2kC\epsilon$ and $e^{-\eta} - 1 \ge -2kC\epsilon$. This gives us $-2kC\epsilon v \le x - v \le 2kC\epsilon v$ and this implies that $\|x - v\|_1 \le 2kC\epsilon$.<fnref n="6" />

## 5.2 Proof of Theorem 3.2

To apply Theorem 5.1 to our map $\mathcal{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$.

**Lemma 5.5.** There exist BHN-constants $\mathcal{G}, \mathcal{G}', \mathcal{G}'', \mathcal{G}'''$, such that for $L \ge \exp(\mathcal{G}'')$ and for all $\psi \in \mathcal{H}_{[1,L-1]}$ with $\|\psi\| = 1$ we have

$$ \mathcal{L}_L^{*\mathcal{G}L}(\psi\psi^*) \ge \frac{1}{\mathcal{G}'}e^{-\mathcal{G}'''L}1 $$

**Remark 5.6.** We believe that Lemma 5.5 can be improved. In particular, it is natural to think that $\mathcal{L}_L^*$ is primitive for all $L \ge 2$. In fact, we conjecture that for any $L \ge 2$, and any $\psi \in \mathcal{H}_{[1,L-1]}$ such that $\|\psi\| = 1$, $(\mathcal{L}_L^*)^{L-1}(\psi\psi^*) \ge 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 $\mathcal{L}_L^*$ and their adjoint, which we recall now.

**Limit of $\mathcal{L}_L$** Note that the adjoint of the map $\mathcal{L}_L^*$ is given by

$$ \mathcal{L}_L(Q) = \text{tr}_1(E_L^\dagger(Q \otimes 1)E_L). \qquad (23) $$

The key additional fact that we will require from [24] is that if $\mathcal{L}$ is the limit of the maps $(\mathcal{L}_L)$, then the iterations $\mathcal{L}^n$ converge, as $n \to \infty$ to a “rank-one” operator (property (iii) below). More precisely, we have<fnref n="7" />:

**Lemma 5.7 ([29] Theorem 4.4, 4.7 and 4.12).** There exists a map $\mathcal{L} : \mathcal{A}_N \to \mathcal{A}_N$, a positive operator $g \in \mathcal{A}_N$, a BHN-constant $\mathcal{G}$, and a superexponentially decaying function $\epsilon(L)$ such that for all $L$ and $X \in \mathcal{A}_{[1,L-1]}$

$$ (i) \quad \|\mathcal{L}(X)\|, \|\mathcal{L}_L(X)\| \le \mathcal{G}\|X\| $$

$$ (ii) \quad \|\mathcal{L}(X) - \mathcal{L}_L(X)\| \le \epsilon(L)\|X\| $$

$$ (iii) \quad \|e^{nf}\mathcal{L}^n(X) - \text{tr}[\rho_L X]g\| \le \mathcal{G}e^{-n/\mathcal{G}}2^{L+1}\|X\| $$

$$ (iv) \quad \|g\|, \|g^{-1}\| \le \mathcal{G} $$

We are now ready to prove Lemma 5.5.

<footnote n="6">⁶Indeed, if $-B \le A \le B$ then necessarily $\|A\|_1 \le \text{tr} B$. This follows by noting that if $A = \sum_\lambda \lambda P_\lambda$ is a spectral decomposition of $A$, then $-\text{tr}(BP_\lambda) \le \text{tr}(AP_\lambda) \le \text{tr}(BP_\lambda)$ and so $\|A\|_1 = \sum_\lambda |\text{tr}(AP_\lambda)| \le \sum_\lambda \text{tr}(BP_\lambda) = \text{tr}(B)$.</footnote>

<footnote n="7">⁷We use the notation $\mathcal{A}_N$ for the closure, with respect to the operator norm, of the union of all local algebras $\mathcal{A}_{[1,n]}$, where operators in $\mathcal{A}_{[1,n]}$ are embedded into $\mathcal{A}_{[1,m]}$ for $m > n$ by tensoring with the identity.</footnote>
Proof of Lemma 5.5. Throughout the proof $\mathcal{G}_i$ are BHN-constants. As a first step we note that we can equivalently prove the statement for the adjoint channel:

$$ 
\begin{align*}
\forall \psi \in \mathcal{H}_{[1,L-1]} \text{ with } \|\psi\| = 1, \quad \mathcal{L}_L^{*\mathcal{G}_i}(\psi\psi^*) &\ge a \mathbf{1} \\
\Leftrightarrow \forall \psi, \phi \in \mathcal{H}_{[1,L-1]}, \quad \|\psi\| = \|\phi\| = 1, \quad \langle \mathcal{L}_L^{*\mathcal{G}_i}(\psi\psi^*), \phi\phi^* \rangle &\ge a \\
\Leftrightarrow \forall \psi, \phi \in \mathcal{H}_{[1,L-1]}, \quad \|\psi\| = \|\phi\| = 1, \quad \langle \psi\psi^*, \mathcal{L}_L^{\mathcal{G}_i}(\phi\phi^*) \rangle &\ge a \\
\Leftrightarrow \forall \phi \in \mathcal{H}_{[1,L-1]} \text{ with } \|\phi\| = 1, \quad \mathcal{L}_L^{\mathcal{G}_i}(\phi\phi^*) &\ge a \mathbf{1}.
\end{align*}
$$

Now using $\mathcal{L}$ and $g$ from Lemma 5.7 and introducing an integer $M$ that we choose later,

$$ 
\begin{align*}
\mathcal{L}_L^M(\phi\phi^*) &= \mathcal{L}^M(\phi\phi^*) + (\mathcal{L}_L^M(\phi\phi^*) - \mathcal{L}^M(\phi\phi^*)) \\
&= e^{-fM}(\mathrm{tr}[\rho_L\phi\phi^*]g + (e^{fM}\mathcal{L}^M(\phi\phi^*) - \mathrm{tr}[\rho_L\phi\phi^*]g) + e^{fM}(\mathcal{L}_L^M(\phi\phi^*) - \mathcal{L}^M(\phi\phi^*)) \\
&\ge e^{-\|h\|^M} \mathbf{1}(\mathcal{G}_1^{-1}e^{-2\|h\|^L} - \mathcal{G}_2e^{-M/\mathcal{G}_2}2^{L+1} - \mathcal{G}_3^M\varepsilon(L)),
\end{align*}
$$

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 $\varepsilon(L)$, which follows from the following telescope sum:

$$ 
\begin{aligned}
\|\mathcal{L}_L^M(\phi\phi^*) - \mathcal{L}^M(\phi\phi^*)\| &= \left\| \sum_{k=0}^{M-1} \mathcal{L}^k(\mathcal{L}_L(\mathcal{L}_L^{M-k-1}(\phi\phi^*)) - \mathcal{L}_L(\mathcal{L}_L^{M-k-1}(\phi\phi^*)) ) \right\| \\
&\le M \mathcal{G}_4^M \varepsilon'(L).
\end{aligned}
$$

We used Lemma 5.7, which proves the existence of a BHN-constant $\mathcal{G}_4$ and a superexponentially decaying function $\varepsilon'$ such that $\|\mathcal{L}_L(Q)\|, \|\mathcal{L}(Q)\| \le \mathcal{G}_4\|Q\|$ and $\|\mathcal{L}_L(Q) - \mathcal{L}(Q)\| \le \varepsilon'(L)$. The factor $M$ can also be upper bounded by $\mathcal{G}_4^M$, which we can combine into a new BHN-constant $\mathcal{G}_3 = \mathcal{G}_4^2 e^f$.

We choose $M = \lceil \mathcal{G}_2(2\|h\| + \log(4)) \rceil L := \mathcal{G}L$ and express the function $\varepsilon$ explicitly to obtain

$$ 
\mathcal{L}_L^M(\phi\phi^*) \ge e^{-\|h\|^M} \mathcal{G}_1^{-1} \mathbf{1} \\
\times \left( e^{-2\|h\|^L} - 2\mathcal{G}_1\mathcal{G}_2e^{-2\|h\|^L} \left(\frac{1}{2}\right)^L - \mathcal{G}_1e^{(\mathcal{G}\log(\mathcal{G}_5)+C_3(\|h\|+\log(d))-\frac{1}{2}\log(L/2))L} \right).
$$

where $C_3$ is a numerical constant. We now impose that

$$ L \ge \log_2(8\mathcal{G}_1\mathcal{G}_2) $$

to bound the second term in brackets by $\frac{1}{4}e^{-2\|h\|^L}$ and similarly

$$ 
\begin{align*}
L &\ge 2 \exp(2(\mathcal{G} \log(\mathcal{G}_5) + C_3(\|h\| + \log(d)) + 2\|h\| + \log(2))) \\
L &\ge \log_2(4\mathcal{G}_1)
\end{align*}
$$

to ensure the same bound for the last term. We combine the above bounds into a new BHN-constant $\mathcal{G}''$ such that $L \ge \exp(\mathcal{G}'')$ implies all the above inequalities. Under this condition we conclude

$$ 
\mathcal{L}_L^{\mathcal{G}''}(\phi\phi^*) \ge e^{-\|h\|^M} \frac{1}{2\mathcal{G}_1} e^{-2\|h\|^L} \ge \frac{1}{\mathcal{G}'} e^{-\mathcal{G}'''L}
$$

with $\mathcal{G}' = 2\mathcal{G}_1$ and $\mathcal{G}''' = 3\|h\|\mathcal{G}$.

□
Proof of Theorem 3.2. We consider again BHN-constants $\mathcal{G}_i$. We know from Lemma 5.7 that $\mathcal{L}_L^*(Q) \le \mathcal{G}_1\|Q\|\mathbf{1}$, and so this implies that for the iterated map $(\mathcal{L}_L^*)^{\mathcal{G}_2L}(Q) \le \mathcal{G}_1^{\mathcal{G}_2L}\|Q\|\mathbf{1}$. Together with Lemma 5.5 this implies that for any pure state $Q = \psi\psi^*$,

$$ \frac{\lambda_{\max}((\mathcal{L}_L^*)^{\mathcal{G}_2L}(Q))}{\lambda_{\min}((\mathcal{L}_L^*)^{\mathcal{G}_2L}(Q))} \le \mathcal{G}_3 e^{\mathcal{G}_4L}. $$

Furthermore, from Corollary 4.4 we know that there is a superexponentially decaying function $\varepsilon(L)$ such that

$$ (1 + \varepsilon(L))^{-1} e^{-f} \rho_L \le \mathcal{L}_L^*(\rho_L) \le (1 + \varepsilon(L)) e^{-f} \rho_L. $$

Applying Theorem 5.1, we get that for $L \ge \exp(\mathcal{G}'')$, $\|v_L - \rho_L\|_1 \le 2\mathcal{G}_2L\mathcal{G}_3e^{\mathcal{G}_4L}\varepsilon(L) \le e^{\mathcal{G}L}\varepsilon(L)$ for some BHN-constant $\mathcal{G}$ and the superexponentially decaying function $\varepsilon(L)$. $\square$

# 6 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 $\mathcal{L}_L^*$ is at most

$$ \mathcal{G} \frac{\exp(C' L)}{(L/2)!}, \qquad (24) $$

where $\mathcal{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 $\varepsilon$.

As the expression in (24) cannot be inverted analytically to get $L$ in terms of $\varepsilon$, we use Stirling’s approximation to upper bound (24). Thus a sufficient condition on our $L$ is that it satisfies:

$$ \varepsilon \ge \mathcal{G} \exp \left( \left( C - \frac{1}{2} \log \frac{L}{2} \right) L \right). $$

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

$$ \log \left( \frac{\varepsilon}{\mathcal{G}} \right) \ge \exp \left( 2 \left( C + \gamma - \frac{1}{2} \right) \right) - \gamma L $$

is fulfilled for some $\gamma$. Using $\gamma = \log(\log(1/\varepsilon))/2$, one can check that the value

$$ L = \left[ \frac{\log(1/\varepsilon)(2 + 2 \exp(2C - 1)) + 2 \log(\mathcal{G})}{\log(\log(1/\varepsilon))} \right] $$

will work.

The algorithm is now given by the computation of the spectral radius of a square matrix of size $d^{2(L-1)} \times d^{2(L-1)}$, which can be done in time polynomial in its size. Making the constants explicit again, using the assumption $\varepsilon < 1/e$ to include the $\log(\mathcal{G})$ in the prefactor of $\log(1/\varepsilon)$, and absorbing all numbers into the new constant $C_1 > 0$ the resulting runtime then reads

$$ \exp \left( \log(d) \left( \frac{\log(1/\varepsilon) \exp(C_1(\|h\| + 1))}{\log(\log(1/\varepsilon))} \right) \right). $$
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, C_2 > 0$

$$ \exp\left(\log(d) \left( \frac{\log(1/\varepsilon) \exp(C_1 e^{C_2 \|h\|})}{\log(\log(1/\varepsilon))} \right) \right). $$

In addition, we have to choose $L \ge \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. □

# 7 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

$$ H = -\beta \sum_{i=1}^{N/2} \left[ (S_{2i-1}^x S_{2i}^x + S_{2i-1}^y S_{2i}^y) + \gamma (S_{2i}^x S_{2i+1}^x + S_{2i}^y S_{2i+1}^y) \right], \quad (25) $$

where the spin operators are given by the Pauli matrices $S^{x,y} = \sigma^{x,y}/2$. In the case $\gamma = 1$, this model is translation-invariant and is given by the interaction term $h = -\frac{\beta}{4}(\sigma^x \otimes \sigma^x + \sigma^y \otimes \sigma^y)$. For $\gamma \ne 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

$$ h = -\frac{\beta}{4} ((\sigma^x \otimes \sigma^x + \sigma^y \otimes \sigma^y) \otimes \mathbf{1} \otimes \mathbf{1} + \gamma \mathbf{1} \otimes (\sigma^x \otimes \sigma^x + \sigma^y \otimes \sigma^y) \otimes \mathbf{1}). $$

Note also that the system is gapless for $\gamma = 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 $\beta f_\beta(h)$ and its dependence on temperature and the parameter $L$ in Figure 1 for $\gamma = 1$. The value of the log-partition function itself as well as the expectation value of the energy is depicted in Figure 2 for $\gamma = 2$. We compute the energy

$$ e_{\beta} = \lim_{N \to \infty} \frac{\text{tr}[h_{1,2} e^{-H_{[-N,N]}}]}{\text{tr}[e^{-H_{[-N,N]}}]} $$

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 $\mathcal{L}_L^*$. In Figure 4, we show the decay of the mutual information

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

between two particles depending on their distance and the same for the conditional mutual information

$$ I(A:C|B) = S(\rho_{AB}) + S(\rho_{BC}) - S(\rho_{ABC}) - S(\rho_B) $$

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 $\mathcal{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><img src="image_4.png" /><figcaption>Figure 1: (a) Errors on logarithmic scale plotted against the inverse temperature for the model (25) with $\gamma = 1$. The errors decay with $L$ but grow with $\beta$. While our estimates are doubly exponential in $\beta$ the shown curves are sublinear suggesting that in practice errors only grow subexponentially with $\beta$. The local minimum for $L = 5$ is due to a change of sign in the error. (b) The decay of the error with $L$ for $\gamma = 1$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 2: (a) Log-partition function plotted against the inverse temperature for the model (25) with $\gamma = 2$. The curves for all $L$ are not distinguishable from the exact solution for high temperatures but deviate for lower temperatures, where higher $L$ still give better approximations. Due to the choice of $\gamma$ we use a blocking procedure resulting in a local dimension $d = 4$, which restricts us to smaller values for the parameter $L$ compared to the case $\gamma = 1$ of Figure 1. (b) Energy plotted against temperature for $\gamma = 1$. Again, we obtain accurate results at high temperatures.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: The error of the energy at $\gamma = 1$ when comparing to the exact solution. The dips again arise from a change of sign of the error. The growth with $\beta$ does not seem to be worse than for the free energy.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Decay of the conditional mutual information (a) and the mutual information (b) for the two-sided chain and $\gamma = 1$ between two particles with distance $d(\cdot, \cdot)$. The conditioning is on the remaining particles in the 8-site marginal. We observe a growth of the decay rate with temperature.</figcaption></figure>

## Acknowledgements

We thank Álvaro Alhambra for useful discussions, Arsalan Motamedi for pointing out an error in Lemma 4.1, and the anonymous referees for their careful assessment and helpful feedback. HF acknowledges funding from UK Research and Innovation (UKRI) under the UK government's Horizon Europe funding guarantee EP/X032051/1. OF acknowledges funding by the European Research Council (ERC Grant AlgoQIP, Agreement No. 851716). SOS acknowledges support from the UK Engineering and Physical Sciences Research Council (EPSRC) under grant number EP/W524141/1.

## References

[1] Julia Kempe, Alexei Kitaev, and Oded Regev. “The complexity of the local hamiltonian problem”. SIAM Journal on Computing 35, 1070–1097 (2006).

[2] Sevag Gharibian, Yichen Huang, Zeph Landau, Seung Woo Shin, et al. “Quantum hamiltonian complexity”. Foundations and Trends® in Theoretical Computer Science 10, 159–282 (2015).

[3] Daniel Gottesman and Sandy Irani. “The quantum and classical complexity of translationally invariant tiling and hamiltonian problems”. Theory of Computing 9, 31–116 (2013).

[4] Johannes Bausch, Toby Cubitt, and Maris Ozols. “The complexity of translationally invariant spin chains with low local dimension”. Annales Henri Poincaré 18, 3449–3513 (2017).

[5] Álvaro M. Alhambra. “Quantum many-body systems in thermal equilibrium” (2022). arXiv:2110.15466.

[6] Anders W. Sandvik. “Computational studies of quantum spin systems”. In AIP Conference Proceedings. Volume 1297, pages 135–338. American Institute of Physics (2010).
[7] James D. Watson and Toby S. Cubitt. "Computational complexity of the ground state energy density problem". In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing. Page 764-775. STOC 2022New York, NY, USA (2022). ACM.

[8] Dorit Aharonov and Sandy Irani. "Hamiltonian complexity in the thermodynamic limit". In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing. Page 750-763. STOC 2022New York, NY, USA (2022). ACM.

[9] Tomotaka Kuwahara, Kohtaro Kato, and Fernando G. S. L. Brandão. "Clustering of conditional mutual information for quantum gibbs states above a threshold temperature". Physical Review Letters 124, 220601 (2020).

[10] Aram W. Harrow, Saeed Mehraban, and Mehdi Soleimanifar. "Classical algorithms, correlation decay, and complex zeros of partition functions of quantum many-body systems". In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing. Page 378-386. ACM (2020).

[11] Ryan L. Mann and Tyler Helmuth. "Efficient algorithms for approximating quantum partition functions". Journal of Mathematical Physics 62, 022201 (2021).

[12] Tomotaka Kuwahara and Keiji Saito. "Polynomial-time classical simulation for one-dimensional quantum gibbs states" (2018). arXiv:1807.08424.

[13] M. B. Hastings. "Quantum belief propagation: An algorithm for thermal quantum systems". Physical Review B 76, 201102 (2007).

[14] Andreas Bluhm, Ángela Capel, and Antonio Pérez-Hernández. "Exponential decay of mutual information for gibbs states of local hamiltonians". Quantum 6, 650 (2022).

[15] Tomotaka Kuwahara, Álvaro M. Alhambra, and Anurag Anshu. "Improved thermal area law and quasilinear time algorithm for quantum gibbs states". Physical Review X 11, 011047 (2021).

[16] Álvaro M. Alhambra and J. Ignacio Cirac. "Locally accurate tensor networks for thermal states and time evolution". PRX Quantum 2, 040331 (2021).

[17] Sergey Bravyi, Anirban Chowdhury, David Gosset, and Pawel Wocjan. "On the complexity of quantum partition functions" (2021). arXiv:2110.15466.

[18] Matthias Troyer, Stefan Wessel, and Fabien Alet. "Flat histogram methods for quantum systems: Algorithms to overcome tunneling problems and calculate the free energy". Physical Review Letters 90, 120201 (2003).

[19] M Suzuki. "Quantum monte carlo methods in condensed matter physics". World Scientific. (1993).

[20] Steven R. White. "Density matrix formulation for quantum renormalization groups". Physical Review Letters 69, 2863-2866 (1992).

[21] Ulrich Schollwöck. "The density-matrix renormalization group in the age of matrix product states". Annals of Physics 326, 96-192 (2011).

[22] Zeph Landau, Umesh Vazirani, and Thomas Vidick. "A polynomial time algorithm for the ground state of one-dimensional gapped local hamiltonians". Nature Physics 11, 566-569 (2015).

[23] Sacha Friedli and Yvan Velenik. "Statistical mechanics of lattice systems: A concrete mathematical introduction". Cambridge University Press. (2017).
[24] Huzihiro Araki. “Gibbs states of a one dimensional quantum lattice”. Communications in Mathematical Physics 14, 120–157 (1969).

[25] Rajendra Bhatia. “Matrix analysis”. Volume 169. Springer Science & Business Media. (2013).

[26] Allan Sly. “Computational transition at the uniqueness threshold”. In 2010 IEEE 51st Annual Symposium on Foundations of Computer Science. Pages 287–296. (2010).

[27] Allan Sly and Nike Sun. “The computational hardness of counting in two-spin models on d-regular graphs”. In 2012 IEEE 53rd Annual Symposium on Foundations of Computer Science. Pages 361–369. (2012).

[28] Masanori Ohya and Denes Petz. “Quantum entropy and its use”. Theoretical and Mathematical Physics. Springer. Berlin, Germany (2004).

[29] David Pérez-García and Antonio Pérez-Hernández. “Locality estimates for complex time evolution in 1d”. Communications in Mathematical Physics 399, 929–970 (2023).

[30] Marco Lenci and Luc Rey-Bellet. “Large deviations in quantum lattice systems: One-phase region”. Journal of Statistical Physics 119, 715–746 (2005).

[31] David E. Evans and Raphael Høegh-Krohn. “Spectral Properties of Positive Maps on $C^*$-Algebras”. Journal of the London Mathematical Society s2-17, 345–355 (1978).

[32] Michael M Wolf. “Quantum channels & operations: Guided tour” (2012).

[33] David Reeb, Michael J Kastoryano, and Michael M Wolf. “Hilbert’s projective metric in quantum information theory”. Journal of mathematical physics 52, 082201 (2011).

[34] Elon Kohlberg and John W Pratt. “The contraction mapping approach to the Perron-Frobenius theory: Why Hilbert’s metric?”. Mathematics of Operations Research 7, 198–210 (1982).

[35] L. N. Bulaevskii. “Theory of non-uniform antiferromagnetic spin chains”. Journal of Experimental and Theoretical Physics 17, 684 (1963).

[36] R J Bursill, T Xiang, and G A Gehring. “The density matrix renormalization group for a quantum spin chain at non-zero temperature”. Journal of Physics: Condensed Matter 8, L583–L590 (1996).

[37] Walter Rudin. “Principles of mathematical analysis”. International series in pure and applied mathematics. McGraw-Hill Professional. New York, NY (1976). 3 edition.

## 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 $\Phi$.

$$ \begin{aligned} \Omega_0 &= \Omega_1 = \Omega_2 = 2\|h\| \\ 0 &= \Omega_3 = \Omega_4 = \dots \\ \|\Phi\|_{\lambda} &= 2\|h\|(1 + e^{\lambda} + e^{2\lambda}) \end{aligned} $$
$$ \Omega_k^*(2) = \sum_{n=1}^{\infty} \left( \sum_{\substack{\alpha \in \mathbb{N}^n \\ |\alpha|=k}} \prod_{j=1}^{k} \Omega_{\alpha_j} \right) \frac{2^n}{n!} $$

The following constant appears in several results cited in our main text and is a BHN-constant due to the bound in [29, Section 4.2]:

$$ \mathcal{G} = \exp \left( \sum_{k \ge 1} e^{2\Omega_0 k} \Omega_k^*(2) \right) \le \exp(e^{2\|\Phi\|_\lambda}) $$

To obtain the superexponential decay of $\mathcal{G}_l$ we start by upper bounding $\Omega_k^*(2)$. To that end a generally useful formula is the Lagrange remainder bound for the tail of the exponential series:

$$ \sum_{k \ge l} \frac{x^k}{k!} = \exp(x) - \sum_{k=1}^{l-1} \frac{x^k}{k!} \le \frac{e^x x^k}{k!} $$

Note that the terms in the outer summation in the definition of $\Omega_k^*(2)$ are zero whenever $n < \lceil k/2 \rceil$ as for those one of the factors has to be zero. Furthermore, the number of vectors in $\mathbb{N}^n$ with entries no larger than 2 is bounded by $3^n$ and we obtain

$$ \Omega_k^*(x) \le \sum_{n=\lceil k/2 \rceil}^{\infty} 3^n (2\|h\|)^n \frac{x^n}{n!} \le e^{6x\|h\|} \frac{(6x\|h\|)^{\lceil k/2 \rceil}}{\lceil k/2 \rceil!} $$

Plugging this into the expression for $\mathcal{G}_l$ yields

$$ \begin{aligned} \mathcal{G}_l &\le e^{2\|\Phi\|} \sum_{k \ge l} e^{2\Omega_0 k} e^{12\|h\|} \frac{(12\|h\|)^{\lceil k/2 \rceil}}{\lceil k/2 \rceil!} \\ &\le e^{2\|\Phi\|_\lambda + 12\|h\|} 2 \sum_{k \ge \lceil l/2 \rceil} e^{4\Omega_0 k} \frac{(12\|h\|)^k}{k!} \\ &\le e^{2\|\Phi\|_\lambda + 12\|h\|} 2 \frac{(e^{8\|h\|} 12\|h\|)^{\lceil l/2 \rceil}}{\lceil l/2 \rceil!} \end{aligned} $$

which gives us the superexponential-decay of this quantity.

We finally need a bound on the constants $K_2$ and $\delta_2$ from [29, Theorem 4.12] which relies on [29, Theorem 4.11]. In the later theorem we have

$$ C_2 = 2K_2 $$

if $r = \frac{\log(a)}{\log(2)}$ and thereby $N = 3r = \frac{3\log(a)}{\log(2)}$. $K_2$ is bounded by

$$ K_2 \le 2\mathcal{G}^4 \left( 1 + \sum_{l \ge 1} \mathcal{G}_l 2^l \right) \mathcal{G}^4 (2\mathcal{G}^3 \mathcal{G}_r 2^r + 2\mathcal{G}^4) + 2\mathcal{G}^8 \left( 1 + \sum_{l \ge 1} \mathcal{G}_l 2^l \right), $$

which is a BHN-constant if $1 + \sum_{l \ge 1} \mathcal{G}_l 2^l$ is as well. This is the case because, looking at the $l$-dependent term in $\mathcal{G}_l$

$$ 1 + \sum_{l \ge 1} \frac{e^{8\|h\|} 12\|h\|^{\lceil l/2 \rceil}}{\lceil l/2 \rceil!} \le 2 \exp(e^{8\|h\|} 12\|h\|) $$

which is BHN and the coefficient in $\mathcal{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

$$\|L^r\|_{1,2} = \sup_{\|Q\|_{1,2} \le 1} \|L^r(Q)\|_{1,2} = \sup_{\|Q\|_{1,2} \le 1} \|L^r(Q)\| + \sum_{l \ge 1} \|L(Q)\|_l,$$

we apply [29, Corollary 4.9] to bound

$$\|L^r\| \le \sup_{\|Q\|_{1,2}} \mathcal{G}^4 \|Q\| + \sum_{l \ge 1} 2 \mathcal{G}^3 \mathcal{G}_l \|Q\| + \sum_{l \ge 1} \mathcal{G}^4 \|Q\|_{n+l}$$

where the first sum is BHN due to the decay of the $\mathcal{G}_l$ as before and the second sum is upper bounded by $\mathcal{G}^4 \|Q\|_{1,2} \le \mathcal{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 \le 2\mathcal{G}^3\mathcal{G}_l$, which we can sum to obtain

$$\|h\|_{1,2} \le \|h\| + \sum_{l \ge 1} \|h\|_l 2^l$$

which is now BHN as the sum converges as before. This closes the argument for $K$. The constant $\delta$ is given as

$$\delta = -\log\left(1 - \frac{1}{2C_2}\right)/N$$

with $C_2$ and $N = \frac{3}{\log(2)} \log(a) = \frac{3}{\log(2)} \log(C_2)$ We estimate using the concavity of the logarithm and its derivative

$$-\log\left(1 - \frac{1}{2C_2}\right) = \log(2C_2) - \log(2C_2 - 1) \ge \frac{1}{2C_2},$$

so finally

$$\delta \ge \frac{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} \in \mathcal{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

$$\mu = \lim_{N \to \infty} \frac{\mathrm{tr} [P_{1,2} \exp(-H_{[-N,N]})]}{\mathrm{tr}[\exp(-H_{[-N,N]})]}, \quad (26)$$

where we assume $\|P\| \le 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 $\epsilon$, runs in time at most

$$ \exp \left( \log(d) \exp(C_1 e^{C_2 \|h\|}) \frac{\log(1/\epsilon)}{\log(\log(1/\epsilon))} \right), $$

and outputs an approximation $\tilde{\mu}$ such that $|\mu - \tilde{\mu}| \le \epsilon$, where $\mu$ 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.

Given a 2-local Hamiltonian, we define the time-evolution operator

$$ \Gamma^s_{[a,b]}(A) = e^{isH_{[a,b]}} A e^{-isH_{[a,b]}}, $$

Also, we define an $l$-neighbourhood of a particle $j \in [1, N]$ as $\Lambda_j^l = [\max\{1, j-l\}, \min\{N, j+1+l\}]$.

We also need the following notation for the thermal state of the finite system on $[a, b]$:

$$ \tilde{\rho}_{[a,b]} = \frac{e^{-H_{[a,b]}}}{\mathrm{tr}[e^{-H_{[a,b]}}]}. $$

The results needed for our proof are summarized in the following Lemma.

**Lemma B.2** ([29, Theorems 4.16 and 2.4] and [14, Theorem 6.2]). Consider a 2-local translation invariant Hamiltonian of the form $H_{[a,b]} = \sum_{i=a}^{b-1} h_{i,i+1}$. Then the following holds:

(i) There exist constants $C, \delta > 0$, such that for any normalized operator $P_{[a,b]} \in \mathcal{A}_{[a,b]}$, and any integer $n \ge 1$, we have

$$ \left| \mathrm{tr} [P_{[a,b]} \tilde{\rho}_{[a-n,b-n]}] - \lim_{N \to \infty} \mathrm{tr} [P_{[a,b]} \tilde{\rho}_{[-N,N]}] \right| \le C e^{-\delta n}. $$
(ii) There is a superexponentially decaying function $\eta(l)$, such that for any normalized observable $P_{j,j+1} \in \mathcal{A}_{[j,j+1]}$

$$ \left\| \Gamma_{[1,N]}^{is}(P_{j,j+1}) - \Gamma_{\Lambda_j^l}^{is}(P_{j,j+1}) \right\| \le \eta(l). $$

(iii) There are BHN-constants $\mathcal{G}$, $K$, $1/\alpha$ such that for any observables $O_A, O_B$ with supports separated by $l$ sites we have

$$ \left| \mathrm{tr}[O_A O_B \tilde{\rho}_{[1,N]}] - \mathrm{tr}[O_A \tilde{\rho}_{[1,N]}] \mathrm{tr}[O_B \tilde{\rho}_{[1,N]}] \right| \le \|O_A\| \|O_B\| \left( K e^{-\alpha l} + \mathcal{G} \frac{\mathcal{G}^l}{(l/2)!} \right). $$

## B.2 Proof of Lemma B.1

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

$$ H_{[1,N]}(\varepsilon) = \sum_{i=1}^{N-1} (h_{i,i+1} + \varepsilon P_{i,i+1}). $$

We denote the corresponding partition function $Z_N(\varepsilon) = \mathrm{tr} e^{-H_{[1,N]}(\varepsilon)}$ and free energy per site $f_N(\varepsilon) = -\frac{1}{N} \log Z_N(\varepsilon)$, as well as the limit

$$ f(\varepsilon) = \lim_{N \to \infty} f_N(\varepsilon). \qquad (27) $$

We also make use of the perturbed thermal state associated to the finite system on the interval $[a, b]$ as

$$ \tilde{\rho}_{\varepsilon, [a, b]} = \frac{\exp(-H_{[a, b]}(\varepsilon))}{\mathrm{tr}[\exp(-H_{[a, b]}(\varepsilon))]} $$

For convenience, we let $\tilde{\rho}_{\varepsilon, N} = \tilde{\rho}_{\varepsilon, [1, N]}$.

The first lemma shows that the expectation value $\mu$ defined in (26) is equal to the derivative of the free energy for the perturbed Hamiltonian $f(\varepsilon)$ defined in (27).

**Lemma B.3.** For any 2-local observable $P \in \mathcal{A}_{[1,2]}$ we have

$$ \frac{\mathrm{d}}{\mathrm{d}\varepsilon} f(0) = \mu := \lim_{N \to \infty} \frac{\mathrm{tr} [P_{1,2} \exp(-H_{[-N,N]})]}{\mathrm{tr}[\exp(-H_{[-N,N]})]}, $$

where $f(\varepsilon)$ 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$,

$$ \frac{\mathrm{d}}{\mathrm{d}\varepsilon} f_N(\varepsilon) = \mathrm{tr} \left[ \frac{1}{N-1} \sum_{i=1}^{N-1} P_{i,i+1} \tilde{\rho}_{\varepsilon, N} \right]. $$

Define $\mu(\varepsilon) = \lim_{N \to \infty} \mathrm{tr}[P_{1,2} \tilde{\rho}_{\varepsilon, [-N,N]}]$. Lemma B.2 (i) tells us that for any $i \in [1, N]$,

$$ |\mathrm{tr}[P_{i,i+1} \tilde{\rho}_{\varepsilon, N}] - \mu(\varepsilon)| \le C e^{-\delta \min(i, N-i)} $$
for some constants $C, \delta > 0$. It thus follows that

$$ \left| \frac{1}{N-1} \sum_{i=1}^{N-1} \mathrm{tr}[P_{i,i+1} \tilde{\rho}_{\varepsilon,N}] - \mu(\varepsilon) \right| \le \frac{C}{N-1} \sum_{i=1}^{N-1} e^{-\delta \min(i,N-i)} \\ \le \frac{2C}{N-1} \sum_{i=1}^{N/2} e^{-\delta i} \le \frac{2C}{N-1} (M + (N/2 - M)e^{-\delta M}) $$

where in the last step $M$ is any integer in $[1, N/2]$. By taking $M = \sqrt{N/2}$ we see that $\frac{d}{d\varepsilon}f_N(\varepsilon) \to \mu(\varepsilon)$ as $N \to \infty$. To prove that $\frac{d}{d\varepsilon}f_N \to \frac{d}{d\varepsilon}f$, we need the convergence to be uniform on all compact sets. While the constants $C$ and $\delta$ do depend on the interaction strength $\|h\| + \varepsilon$, an upper bound on these constants grows monotonically with $\|h\| + \varepsilon$. Therefore, we just choose the constants for interaction strength $\|h\| + \varepsilon'$ to prove uniform convergence for $\varepsilon \in [0, \varepsilon']$. 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) \to \mathbb{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) \le C$, then the same bound holds for the limit $g'(x) \le C$.

*Proof.* Let us assume for contradiction that there exists $x_0 \in (a, b)$, such that $g'(x_0) > C$. Due to the continuity of $g'$ there exist $\delta, \varepsilon$ such that $g'(x) > C + \delta$ for all $x \in [x_0, x_0 + \varepsilon]$. We now choose $N$ sufficiently large such that $|g_N(x) - g(x)| \le \varepsilon\delta/3$ for $x \in [x_0, x_0 + \varepsilon]$. Then

$$ \frac{g_N(x_0 + \varepsilon) - g_N(x_0)}{\varepsilon} \ge \frac{\varepsilon(C + \delta) - 2\varepsilon\delta/3}{2\varepsilon} \ge C + \frac{1}{3}\delta $$

By the mean value theorem this implies that there exists $x_1 \in [x_0, x_0 + \varepsilon]$ such that $g'_N(x_1) \ge C + \frac{1}{3}\delta$, which contradicts our assumption and closes the proof. □

We can now prove a bound on the second derivative of the free energy using the quasi-locality of complex time-evolved operators and the decay of correlations in the thermal state.

**Lemma B.5.** For each $\varepsilon_0$ there exists a BHN-constant $G$ for Hamiltonian strength $\|h\| + \varepsilon_0$ such that

$$ \frac{d^2}{d\varepsilon^2} f(\varepsilon) \le e^G $$

for all $\varepsilon \in [0, \varepsilon_0]$.

*Proof.* Let us consider the parameter derivatives of the free energy density of a finite system
$$f_N(\varepsilon) = -\frac{1}{N} \log \operatorname{tr} \exp(-H_{[1,N]}(\varepsilon))$$

$$\frac{d}{d\varepsilon} f_N(\varepsilon) = \operatorname{tr} \left[ \frac{1}{N-1} \sum_{i=1}^{N-1} P_{i,i+1} \tilde{\rho}_{\varepsilon,N} \right]$$

$$\frac{d^2}{d\varepsilon^2} f_N(\varepsilon) = \frac{1}{N-1} \left[ \frac{\operatorname{tr} \left[ \sum_{i=1}^{N-1} P_{i,i+1} \int_0^1 e^{-sH_{[1,N]}(\varepsilon)} \sum_{i=1}^{N-1} P_{i,i+1} e^{-(1-s)H_{[1,N]}(\varepsilon)} ds \right]}{\operatorname{tr}[\exp(-H_{[1,N]}(\varepsilon))]} - \left( \operatorname{tr} \left[ \sum_{i=1}^{N-1} P_{i,i+1} \tilde{\rho}_{\varepsilon,N} \right] \right)^2 \right],$$

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

$$\frac{d^2}{d\varepsilon^2} f_N(\varepsilon) = \frac{1}{N-1} \sum_{j,k=1}^{N-1} \operatorname{tr} \left[ P_{j,j+1} \int_0^1 \Gamma_{[1,N]}^{is,\varepsilon}(P_{k,k+1}) ds \tilde{\rho}_{\varepsilon,N} \right] - \operatorname{tr}[P_{j,j+1} \tilde{\rho}_{\varepsilon,N}] \operatorname{tr}[P_{k,k+1} \tilde{\rho}_{\varepsilon,N}],$$

where we used the time evolution operator with respect to the perturbed Hamiltonian

$$\Gamma_{[a,b]}^{s,\varepsilon}(A) = e^{isH_{[a,b]}(\varepsilon)} A e^{-isH_{[a,b]}(\varepsilon)}.$$

We fix $j \in [1, N-1]$ and $s \in [0, 1]$ and focus on the sum

$$\sum_{k=1}^{N-1} \operatorname{tr} \left[ P_{j,j+1} \Gamma_{[1,N]}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] - \operatorname{tr}[P_{j,j+1} \tilde{\rho}_{\varepsilon,N}] \operatorname{tr}[P_{k,k+1} \tilde{\rho}_{\varepsilon,N}].$$

The idea is that the imaginary-time evolved operator $\Gamma_{[1,N]}^{is,\varepsilon}(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 apply Lemma B.2 (ii) and (iii) and choose $l(k) = \max\{(\lfloor j-k \rfloor - 2)/2, 0\}$ to estimate

$$\begin{aligned} & \left| \sum_{k=1}^{N-1} \left( \operatorname{tr} \left[ P_{j,j+1} \Gamma_{[1,N]}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] - \operatorname{tr}[P_{j,j+1} \tilde{\rho}_{\varepsilon,N}] \operatorname{tr}[P_{k,k+1} \tilde{\rho}_{\varepsilon,N}] \right) \right| \\ & \le \left| \sum_{k=1}^{N-1} \left( \operatorname{tr} \left[ P_{j,j+1} \Gamma_{[1,N]}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] - \operatorname{tr} \left[ P_{j,j+1} \Gamma_{\Lambda_k^{l(k)}}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] \right) \right| \\ & \qquad + \left| \sum_{k=1}^{N-1} \left( \operatorname{tr} \left[ P_{j,j+1} \Gamma_{\Lambda_k^{l(k)}}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] - \operatorname{tr} [P_{j,j+1} \tilde{\rho}_{\varepsilon,N}] \operatorname{tr} \left[ \Gamma_{\Lambda_k^{l(k)}}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] \right) \right| \\ & \qquad + \left| \sum_{k=1}^{N-1} \left( \operatorname{tr} [P_{j,j+1} \tilde{\rho}_{\varepsilon,N}] \operatorname{tr} \left[ \Gamma_{\Lambda_k^{l(k)}}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N} \right] - \operatorname{tr}[P_{j,j+1} \tilde{\rho}_{\varepsilon,N}] \operatorname{tr}[P_{k,k+1} \tilde{\rho}_{\varepsilon,N}] \right) \right| \\ & \le \sum_{k=1}^{N-1} \left( \eta(l(k)) + G \frac{G^{l(k)}}{|l(k)|!} + Ke^{-\alpha(l(k))} + \eta(l(k)) \right) \\ & \le e^{G'}. \end{aligned}$$

We used the quasi-locality to bound the first and third term after the first inequality (note that $\operatorname{tr}[P_{k,k+1} \tilde{\rho}_{\varepsilon,N}] = \operatorname{tr}[\Gamma_{[1,N]}^{is,\varepsilon}(P_{k,k+1}) \tilde{\rho}_{\varepsilon,N}]$ as the time evolution operator commutes with
the thermal state). The second term is due to the decay of correlations and the bound on $||\Gamma_{\Lambda_j^{l(k)}}^{is,\varepsilon}(P_{k,k+1})||$ that is also implied by Lemma B.2 (ii). Note that the sum over a superexponentially decaying function is BHN by definition of the exponential series, as is the sum of an exponential when the inverse of the decay rate is BHN. The second term in the sum is responsible for the exponential bound in a BHN-constant again by definition of the exponential series.

By integrating over $s$, summing over $j$, and dividing by $N-1$ this implies that $\frac{d^2}{d\varepsilon^2}f_N(\varepsilon)$ is bounded uniformly in $N$ by the exponential of a BHN constant (for the Hamiltonian $H(\varepsilon)$ whose interaction term has norm up to $\|h\| + \varepsilon$). Using the continuous second-order differentiability of the limit function $f(\varepsilon)$ from [24, Lemma 9.3], the claim follows from Lemma B.4 applied to $\frac{d}{d\varepsilon}f(\varepsilon)$. $\square$

We now have all the tools to estimate thermal expectation values using an approximation of the free energy.

*Proof of Lemma B.1.* We are looking for an approximation $\tilde{\mu}$ such that

$$|\mu - \tilde{\mu}| \le \varepsilon.$$

We have for the second derivative of the free energy from Lemma B.5

$$\left| \frac{d^2}{d\varepsilon^2} f(\varepsilon) \right| \le e^{\mathcal{G}}$$

for all $\varepsilon \in [0, 1]$ and for some constant $\mathcal{G}$ that is BHN for the interaction strength $\|h\| + 1$. To achieve the bound we pick $\varepsilon' = \min\{1, \frac{\varepsilon}{2\exp(\mathcal{G})}\}$ and compute

$$\tilde{\mu} := \frac{\tilde{f}(\varepsilon') - \tilde{f}(0)}{\varepsilon'}$$

where we use an error of $|f(\varepsilon') - \tilde{f}(\varepsilon')|, |f(0) - \tilde{f}(0)| \le \delta := \varepsilon\varepsilon'/4$. Using a standard calculus formula

$$\left| \frac{d}{d\varepsilon} f(0) - \frac{f(\varepsilon') - f(0)}{\varepsilon'} \right| \le \frac{\varepsilon'}{2} \max_{0 \le s \le \varepsilon'} \frac{d^2}{ds^2} f(s),$$

we conclude using our bound on the second derivative and the chosen errors

$$|\mu - \tilde{\mu}| \le \left| \frac{d}{d\varepsilon} f(0) - \frac{f(\varepsilon') - f(0)}{\varepsilon'} \right| + \left| \frac{f(\varepsilon') - \tilde{f}(\varepsilon')}{\varepsilon'} \right| + \left| \frac{f(0) - \tilde{f}(0)}{\varepsilon'} \right| \le \varepsilon.$$

Plugging the desired error into the runtime bounds for the free energy of Corollary 3.3, we conclude that the runtime for the thermal expectation value is bounded by

$$\exp\left(\log(d) \exp(C_1 e^{C_2 \|h\|}) \frac{\log(1/\varepsilon)}{\log(\log(1/\varepsilon))}\right)$$

for some numerical constants $C_1, C_2$. $\square$
C Reducing the Ground Energy Problem to the Free Energy

We reduce the ground-state problem to the free energy problem to establish QMA$_{EXP}$-hardness of the following problem.

**Definition C.1.** The *Free Energy for Infinite Translation-Invariant Hamiltonian FE-ITIH* is defined as follows

**Problem parameter:** Three polynomials $\bar{p}, \bar{q}, \bar{r}$ and $d$ the dimension of a particle

**Problem input:** $N$ specified in binary, $\beta \ge 0$ and the matrix $h$, each specified in binary with at most $\log \bar{r}(N) + 1$ bits

**Promise:** The free energy density of the infinite system defined by $h$, i.e., $f_{\beta}(h) = \lim_{N \to \infty} -\frac{1}{\beta N} \log \text{tr} \exp(-\beta H_{[1,N]})$ and lies outside the interval $[1/\bar{p}(N), 1/\bar{p}(N) + 1/\bar{q}(N)]$

**Output:** Determine if free energy density is at most $1/\bar{p}(N)$ or at least $1/\bar{p}(N) + 1/\bar{q}(N)$

**Lemma C.2.** There exist parameters $\bar{p}, \bar{q}, \bar{r}, \bar{d}$ such that the above problem is QMA$_{EXP}$-hard.

*Proof.* We start by giving a quantitative convergence bound for the free energy to the ground-state energy for $\beta \to \infty$.

Using the well-known variational formula for the free energy, we can write the free energy density for a finite system as

$$f_{\beta,N}(h) = \frac{1}{N} \min_{\rho \ge 0, \text{tr}[\rho]=1} \left( \text{tr}[H_{[1,N]}\rho] - \frac{1}{\beta} S(\rho) \right),$$

where $S(\rho) = -\text{tr}[\rho \log_2 \rho]$ denotes the von Neumann entropy. As $0 \le S(\rho) \le N \log(d)$ we can bound the free energy as

$$\min_{\rho \ge 0, \text{tr}[\rho]=1} \text{tr} \left[ \frac{H_{[1,N]}}{N} \rho \right] - \frac{1}{\beta} \log(d) \le f_{\beta,N} \le \min_{\rho \ge 0, \text{tr}[\rho]=1} \text{tr} \left[ \frac{H_{[1,N]}}{N} \rho \right].$$

Let $e_{0,N} = \frac{1}{N} \min_{\rho \ge 0, \text{tr}[\rho]=1} \text{tr}[H_{[1,N]}\rho]$ so that we have

$$|f_{\beta,N} - e_{0,N}| \le \frac{\log(d)}{\beta}.$$

The existence of the limit $\lim_{N \to \infty} f_{\beta,N}$ shows the existence of the limit ground-state energy $e_0 = \lim_{N \to \infty} e_{0,N}$. As a result, we have

$$|f_{\beta} - e_{0}| \le \frac{\log(d)}{\beta}.$$

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: $\bar{p} = p, \bar{q} = 2 \cdot q, \bar{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, \beta, h)$ where we let $\beta = 2q(N) \log d$. Note that $\beta \le \bar{r}(N)$ and thus can be specified using at most $\log \bar{r}(N) + 1$ bits. We now check that the promise is satisfied: we know that $e_0$ is either $\le \frac{1}{p(N)}$ or at least $\frac{1}{p(N)} + \frac{1}{q(N)}$. Now for the free energy density, as $f_{\beta} \le e_0$, we have either $f_{\beta} \le \frac{1}{p(N)}$ or $f_{\beta} \ge \frac{1}{p(N)} + \frac{1}{q(N)} - \frac{\log d}{\beta} = \frac{1}{p(N)} + \frac{1}{2q(N)}$. It then follows that for the instance we constructed, $e_0 \le \frac{1}{p(N)}$ if and only if $f_{\beta} \le \frac{1}{p(N)}$ which proves the desired result. $\square$

We can interpret this hardness result by saying that we cannot expect to have an algorithm that has a dependence of the form $\exp(\text{polylog}(\beta, 1/\varepsilon))$, unless $QMA_{EXP} = \text{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 Hamil-
tonian.

We consider a local dimension $d$ and 2-particle Hamiltonian $h \in \mathcal{B}(\mathbb{C}^d \otimes \mathbb{C}^d)$ as given. We construct a system where the local Hilbert space corresponds to two copies of the original one, i.e., $\mathcal{H}' = \mathbb{C}^d \otimes \mathbb{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' \in \mathcal{B}(\mathbb{C}^{d'} \otimes \mathbb{C}^{d'})$ is then chosen as $h'_{1,2} = h_{2u,1u} + h_{1d,2d} + h_{1u,1d} - h_{2u,2d}$, 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 $\rho_{i_u, \dots, i_{2u}, 1_u, 1_d, 2_d, 3_d, \dots, i_d}$.

<figure><img src="image_8.png" /><figcaption>Figure 5: The nodes in each column correspond to the bipartite new local Hilbert space. Each black line corresponds to a Hamiltonian term, while the red lines represent a zero due to the cancellation from two neighbouring terms.</figcaption></figure>