# Certified algorithms for equilibrium states of local quantum Hamiltonians

Hamza Fawzi, Omar Fawzi, Samuel O Scalet

► To cite this version:

Hamza Fawzi, Omar Fawzi, Samuel O Scalet. Certified algorithms for equilibrium states of local quantum Hamiltonians. Nature Communications, 2024, 15 (1), pp.7394. 10.1038/s41467-024-51592-3. hal-04682901

HAL Id: hal-04682901

https://hal.science/hal-04682901v1

Submitted on 8 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Certified algorithms for equilibrium states of local quantum Hamiltonians

Received: 4 March 2024

Accepted: 12 August 2024

Published online: 27 August 2024

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

Check for updates

Hamza Fawzi<fnref n="1" />, Omar Fawzi ② & Samuel O. Scalet ①✉

Predicting observables in equilibrium states is a central yet notoriously hard question in quantum many-body systems. In the physically relevant thermodynamic limit, certain mathematical formulations of this task have even been shown to result in undecidable problems. Using a finite-size scaling of algorithms devised for finite systems often fails due to the lack of certified convergence bounds for this limit. In this work, we design certified algorithms for computing expectation values of observables in the equilibrium states of local quantum Hamiltonians, both at zero and positive temperature. Importantly, our algorithms output rigorous lower and upper bounds on these values. This allows us to show that expectation values of local observables can be approximated in finite time, contrasting related undecidability results. When the Hamiltonian is commuting on a 2-dimensional lattice, we prove fast convergence of the hierarchy at high temperature and as a result for a desired precision $\epsilon$, local observables can be approximated by a convex optimization program of quasi-polynomial size in $1/\epsilon$.

A central question in physics is to determine the properties of a many-body quantum system as a function of the interaction between its constituents. The topic of Hamiltonian complexity<fnref n="1" />,² studies this question from a complexity-theoretic point of view. The results in Hamiltonian complexity suggest that efficient algorithms answering this question are unlikely to exist. In fact, determining the ground energy is hard for the complexity class QMA³, and this even holds for translation-invariant systems⁴. This means that we do not expect polynomial-time (classical or quantum) algorithms for this problem. Computing expectation values of observables in the ground state is even harder than the ground energy⁵. In the thermodynamic limit, when the number of systems is taken to infinity, finding good approximations to the energy can be hard for fixed Hamiltonians⁶,⁷, and the spectral gap is even uncomputable⁸.

These complex results put severe limitations on provably efficient and correct classical and quantum algorithms for the fundamental questions in many-body physics. However, these limitations only apply to algorithms that have the required properties for all valid instances of the problem. Moreover, typically the instances showing hardness are highly contrived. In order to avoid these limitations, we consider here certified algorithms where we require the correctness condition

for all instances, but relax the condition of provable efficiency for all instances. We say that an algorithm for computing $p^*(H)$ is certified if on input $H$, it outputs a pair of numbers $(p_l^{\min}, p_l^{\max})$ such that we have $p^*(H) \in [p_l^{\min}, p_l^{\max}]$ for all $H$ and all $l$, i.e., the algorithm provides upper and lower bounds on the quantity of interest. Even without having any performance guarantees a priori, the error is bounded by the size of the interval a posteriori.

We further require that as $l \to \infty$, the interval $[p_l^{\min}, p_l^{\max}]$ converges to the desired value $p^*(H)$. Here, $l$ is a parameter that governs the runtime of the algorithm such that fast convergence in $l$ leads to an efficient algorithm. The above cited complexity results show, e.g., that if $p^*(H)$ is the ground energy of $H$, then for small $l$, the interval $[p_l^{\min}, p_l^{\max}]$ has to be large for some $H$. However, for any input $H$ of interest, we can always run the algorithm and the returned information will be correct and certifiably so. As such, no a priori analysis of the convergence speed is needed in order to obtain rigorous approximations of the quantity of interest.

A well-studied general way of determining properties of thermal states is by preparing such states on a quantum computer. This includes, for example, works on quantum Metropolis sampling⁹,¹⁰ and, more generally, quantum Gibbs samplers¹¹,¹². Here, we focus on

<footnote n="1">¹Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, United Kingdom. ²Univ Lyon, Inria, ENS Lyon, UCBL, LIP, Lyon, France. ✉ e-mail: sos25@cam.ac.uk</footnote>
classical algorithms, which clearly cannot prepare quantum thermal states, but can still compute specified properties of general quantum thermal states. We note that our algorithms have a similar feature to (quantum) Metropolis sampling and, more generally, (quantum) Markov Chain Monte Carlo methods, in that they can be applied for any local Hamiltonian, but this comes at the (unavoidable) cost of not always having fast convergence.

In this work, we consider a hierarchy of convex optimization problems that provide certified algorithms for computing observables in-ground and thermal states in finite and infinite systems. As a main result, we prove the decidability of a formulation of the equilibrium observable problem, which contrasts with previous closely related undecidability results. The main technical ingredient is a formulation of constraints relaxing the Gibbs condition in a convex and algorithmically feasible manner. Furthermore, we are able to prove efficiency in a more restricted setting. Preliminary numerical results also suggest that the method can be of use in practice.

## Results

We now present the setup in an informal way. We refer to the Supplementary information for more formal statements. We assume that our system is described by a local Hamiltonian on a discrete set of sites Γ, which can be formally written as

$$H = \sum_{X \in \Gamma} h_X \quad (1)$$

where $h_X$ is the interaction term acting only on the finite set of sites $X \subset \Gamma$. We assume that the Hilbert space describing each site has finite dimension $d$ and the Hamiltonian $H$ is local, meaning that $h_X$ is nonzero only when $X$ has a size smaller than a constant. For finite systems, the set of equilibrium states at temperature $T \ge 0$ is the set of states $\rho$ that minimize $\text{tr}(H\rho) - TS(\rho)$, where $\text{tr}(H\rho)$ is the energy, and $S(\rho)$ is the entropy. At $T > 0$, this set reduces to a single equilibrium state, also called the Gibbs state given by $\rho = \frac{e^{-H/T}}{\text{tr}e^{-H/T}}$. However, for infinite systems, the thermal equilibrium state is not necessarily unique, and this property lies at the heart of the existence of phase transitions.

Given the description of such a system, our objective is to determine the physical properties of the corresponding equilibrium states. Consider a local observable $O$, e.g., the magnetization in the $z$-direction for a spin-1/2 system for some site $x \in \Gamma$. We are interested in the set of values $\text{tr}(\rho O)$ that an equilibrium state at temperature $T = 1/\beta$ can take:

$$[\langle O \rangle_{\beta}^{\text{min}}, \langle O \rangle_{\beta}^{\text{max}}] := \{\text{tr}(\rho O) : \rho \text{ is an equilibrium state at temperature } T=1/\beta\}. \quad (2)$$

Note that for a finite system and $\beta < \infty$, we have $\langle O \rangle_{\beta}^{\text{min}} = \langle O \rangle_{\beta}^{\text{max}}$ as the thermal state is unique. However, we can have $\langle O \rangle_{\beta}^{\text{min}} < \langle O \rangle_{\beta}^{\text{max}}$ for $\beta = \infty$ if the ground space of $H$ is degenerate and for infinite systems when there are multiple thermal equilibrium states.

## Main result

We formulate a hierarchy of convex optimization programs which produce converging lower (resp. upper) bounds on $\langle O \rangle_{\beta}^{\text{min}}$ (resp. $\langle O \rangle_{\beta}^{\text{max}}$) for any $\beta \ge 0$. Our optimization problems are formulated in terms of the matrix-valued relative entropy function

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

defined for arbitrary positive definite matrices $A, B$, and which is jointly convex in $(A, B)$. Consider a finite subset $\Lambda$ of the lattice sites $\Gamma$ containing the support of the local observable $O$. Then we can formulate the following convex optimization program over density operators supported on $\bar{\Lambda} = \Lambda \cup \partial^{\text{ex}}\Lambda$ where $\partial^{\text{ex}}\Lambda$ is the external

boundary of $\Lambda$ (see the Supplementary information Equation (4) for the precise definition):

$$\min_{\rho} \max_{\rho} \text{tr}(\rho O) \quad (4)$$

$$\text{s.t.} \quad \rho \text{ density operator on } (\mathbb{C}^d)^{\otimes |\bar{\Lambda}|} \quad (5)$$

$$\text{tr}_{\bar{\Lambda} \setminus \Lambda}(H_{\bar{\Lambda}} \rho - \rho H_{\bar{\Lambda}}) = 0 \quad (6)$$

$$D_{op}(A_{\rho} \parallel B_{\rho}) \le \beta C_{\rho}. \quad (7)$$

Here, $H_{\bar{\Lambda}}$ is the truncated Hamiltonian acting on $\bar{\Lambda}$, and $A_{\rho}, B_{\rho}, C_{\rho}$ are Hermitian matrices that depend linearly on $\rho$, defined by:

$$(A_{\rho})_{ij} = \text{tr}(\rho a_i^* a_j), \quad (B_{\rho})_{ij} = \text{tr}(\rho a_j a_i^*), \quad (C_{\rho})_{ij} = \text{tr}(\rho a_i^* [H_{\bar{\Lambda}}, a_j]) \quad (8)$$

where $\{a_i\}$ is a basis of the space of operators acting on $(\mathbb{C}^d)^{\otimes |\bar{\Lambda}|}$. Note that the program (4-7) involves a matrix variable of dimension $d^{|\bar{\Lambda}|} \times d^{|\bar{\Lambda}|}$ and convex constraints involving matrices of dimension at most $d^{2|\bar{\Lambda}|} \times d^{2|\bar{\Lambda}|}$.

By taking larger and larger subsets $\Lambda \uparrow \Gamma$, one can prove that the solutions of the convex optimization programs will converge to the expectation values $\langle O \rangle_{\beta}^{\text{min}}$ and $\langle O \rangle_{\beta}^{\text{max}}$. This is the content of the following theorem.

**Theorem 2.1.** (Certified algorithms for expectation values of equilibrium states). Let $\Lambda_0 \subset \Lambda_1 \subset \dots \subset \Gamma$ be an increasing sequence such that the support of the local observable $O$ is contained in $\Lambda_0$. For any $\ell \in \mathbb{N}$, let $p_{\ell}^{\text{min}}$ and $p_{\ell}^{\text{max}}$ be respectively the minimum and maximum values of the convex optimization problems (4-7) with $\Lambda = \Lambda_{\ell}$. Then we have

$$p_{\ell}^{\text{min}} \le \langle O \rangle_{\beta}^{\text{min}} \le \langle O \rangle_{\beta}^{\text{max}} \le p_{\ell}^{\text{max}}. \quad (9)$$

Furthermore, $p_{\ell}^{\text{min}} \uparrow \langle O \rangle_{\beta}^{\text{min}}$ and $p_{\ell}^{\text{max}} \downarrow \langle O \rangle_{\beta}^{\text{max}}$ as $\ell \to \infty$.

We note that for finite systems, convergence happens after a finite number of steps, namely for $\ell$ such that $\Lambda_{\ell} = \Gamma$. While the convergence result becomes trivial in the finite case, the formulation of a convex program containing constraints for the Gibbs state still adds a novel tool for the numerical treatment of finite-sized quantum many-body systems. The more interesting case for our results is however that of infinite systems. Given a fixed computational budget, one can choose $\ell$ and run the corresponding program (4-7) and obtain a superset $[p_{\ell}^{\text{min}}, p_{\ell}^{\text{max}}]$ of the target interval $[\langle O \rangle_{\beta}^{\text{min}}, \langle O \rangle_{\beta}^{\text{max}}]$. In the cases where $\langle O \rangle_{\beta}^{\text{min}} = \langle O \rangle_{\beta}^{\text{max}}$, e.g., for finite systems at $T > 0$, one can also run the program (4-7) with increasing values of $\ell$ until $p_{\ell}^{\text{max}} \le p_{\ell}^{\text{min}} + \varepsilon$, where $\varepsilon$ is some desired accuracy $\varepsilon$. We note that, as previously mentioned, the problem of computing expectation values for equilibrium states is at least QMA-hard so we cannot hope to have fast convergence for all choices of $H$. But we stress that no a priori analysis of convergence speed is required to obtain some desired accuracy: as soon as we get $p_{\ell}^{\text{max}} \le p_{\ell}^{\text{min}} + \varepsilon$, an additive $\varepsilon$ approximation is guaranteed for this instance. Numerical results illustrating this algorithm can be found in the Supplementary information Section 3.

## Translation-invariant infinite systems

Translation-invariant Hamiltonians on the infinite lattice $\Gamma = \mathbb{Z}^D$, i.e., satisfying $h_X = h_{X+x}$ for all $x \in \mathbb{Z}^D$, play a specifically important role in statistical physics, in particular for understanding phase transitions. For such systems, one often considers the expectation value per site of an observable, for example, the energy per particle also called the energy density. One way to define expectation values per site is to compute the observable $O$ on any fixed site on a translation-invariant
## BOX 1

# Equilibrium observable problem

**Input:** Local dimension *d*, local Hamiltonian term *h* with rational coefficients, local observable *O* with rational coefficients, rational number *a*, temperature *β* ∈ [0, +∞]

**Promise:** Either $\langle O \rangle_{\beta}^{\min, \mathrm{TI}} > a$ or $\langle O \rangle_{\beta}^{\max, \mathrm{TI}} < a$

**Question:** Output YES if $\langle O \rangle_{\beta}^{\min, \mathrm{TI}} > a$, NO otherwise

equilibrium state of *H*. More precisely, we define the average expectation value of observable *O* per site to be the interval

$$ \begin{aligned} & [\langle O \rangle_{\beta}^{\min, \mathrm{TI}}, \langle O \rangle_{\beta}^{\max, \mathrm{TI}}] \\ &= \{\mathrm{tr}(\rho O) : \rho \text{ is a translation-invariant equilibrium state at temperature } T = 1/\beta\}. \end{aligned} \quad (10) $$

By adding to the program (4-7) a translation-invariance constraint in $\bar{\Lambda}$ given by

$$ \mathrm{tr}_{\Sigma^c}(\rho) = \mathrm{tr}_{(\Sigma+x)^c}(\rho) \quad \text{for all } \Sigma \subset \bar{\Lambda}, x \in \mathbb{Z}^d \text{ such that } \Sigma+x \subset \bar{\Lambda}, \quad (11) $$

Theorem 2.1 can be adapted for translation-invariant states, as follows:

**Theorem 2.2.** Let $\Lambda_0 \subset \Lambda_1 \subset \dots \subset \Gamma$ be an increasing sequence such that the support of the local observable *O* is contained in $\Lambda_0$. For any $\ell \in \mathbb{N}$, let $p_{\ell}^{\min, \mathrm{TI}}$ and $p_{\ell}^{\max, \mathrm{TI}}$ be respectively the minimum and maximum values of the convex optimization problems given by (4-7) together with the constraint (11). Then we have

$$ p_{\ell}^{\min, \mathrm{TI}} \le \langle O \rangle_{\beta}^{\min, \mathrm{TI}} \le \langle O \rangle_{\beta}^{\max, \mathrm{TI}} \le p_{\ell}^{\max, \mathrm{TI}}. \quad (12) $$

Furthermore, $p_{\ell}^{\min, \mathrm{TI}} \uparrow \langle O \rangle_{\beta}^{\min, \mathrm{TI}}$ and $p_{\ell}^{\max, \mathrm{TI}} \downarrow \langle O \rangle_{\beta}^{\max, \mathrm{TI}}$ as $\ell \to \infty$.

One difficulty we would like to highlight regarding the thermodynamic limit concerns the definition of equilibrium states at any given temperature $T \ge 0$. It was actually shown in ref. 13 that computing local observables on the ground state of infinite systems is undecidable. This might seem to contradict our result. This is not the case, however, due to a subtle point in the definition of ground-state observables in the infinite limit. The authors of ref. 13 define those as limits of ground-state observables in finite systems, i.e., the observable is computed for the ground-state of a truncated Hamiltonian with open boundary condition $H_A = \sum_{X \subset A} H_X$, where $\Lambda \subset \mathbb{Z}^D$ is finite, and afterward the limit $\Lambda \uparrow \mathbb{Z}^D$ is taken. The problem definition we use in our work, which is standard in the operator algebraic framework, is more general and, in particular, contains limits of ground-state observables of finite system Hamiltonians with any choice of boundary conditions. To avoid confusion, it should be added that this does not mean that our algorithm can compute the value of an observable for a specific boundary condition. Instead, it gives an outer relaxation of the interval given by all boundary conditions. A boundary condition is not input to the algorithm. It is simple to construct a Hamiltonian where fixing the boundary condition while taking the limit excludes natural ground states (see Remark 1.1 and Remark 2.3 in the supplemental information concerning thermal states of the 2D Ising model). The comparison, however, illustrates that even the existence of any convergent algorithm for observables in ground states is far from obvious, which is what we achieve in Theorem 2.2. Considering the operator algebraic definition of equilibrium states leads to the natural computational problem defined in Box 1.

Note that this is a promise problem where the objective is to decide if all the equilibrium states have an expectation value per site $> a$ or all the equilibrium states have an expectation value per site $< a$.

**Theorem 2.3.** (Decidability of translation-invariant ground state and thermal observables). The equilibrium observable problem is decidable.

## Main ingredients

For clarity of the discussion, we focus here on finite systems. A key ingredient to obtain our algorithms is to use an operator algebraic characterization of equilibrium states. These are expressed solely in terms of tr($\rho b$) for operators $b$. For $T=0$, the common definition of an equilibrium state is a state supported on the eigenspace of $H$ with the minimum eigenvalue. It turns out that an equivalent operator algebraic formulation of this condition is:

$$ \mathrm{tr}(\rho a^*[H, a]) \ge 0 \quad \forall a \quad (13) $$

where $[H, a]$ is the commutator and $a$ is an arbitrary observable. Intuitively, the condition above expresses the fact that the energy of $\rho$ has to increase under any perturbation. There are two crucial facts about (13): First, if $a$ is supported on a small set $\Lambda$ of sites then the condition (13) only depends on tr($\rho b$) for operators $b \in \bar{\Lambda}$, where we recall that $\bar{\Lambda} = \Lambda \cup \partial^{\mathrm{ex}}\Lambda$. Second, the inequalities in (13) can be concisely captured by a convex constraint involving the positivity of a Hermitian matrix that depends linearly on $\rho$. By restricting the operators $a$ in (13) to be supported on $\Lambda$, this leads to (4-7) for $T=0$.

At positive temperature $T=1/\beta > 0$, the situation is more complicated. As previously mentioned, the common definition of a thermal equilibrium state is given by $\rho = \frac{e^{-\beta H}}{\mathrm{tr}(e^{-\beta H})}$. It turns out that an equivalent operator algebraic formulation is via the Kubo-Martin-Schwinger (KMS) condition:

$$ \mathrm{tr}(\rho ba) = \mathrm{tr}(\rho a e^{-\beta H} b e^{\beta H}) \quad \forall a, b. \quad (14) $$

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

$$ \mathrm{tr}(\rho|\sigma\rangle\langle\sigma|) = \mathrm{tr}(\rho|\sigma'\rangle\langle\sigma'|) \exp(-\beta(\langle\sigma|H|\sigma\rangle - \langle\sigma'|H|\sigma'\rangle)), \quad (15) $$

for all $\sigma, \sigma'$. When $\sigma'$ is obtained from $\sigma$ by flipping a single spin, these are known as the spin-flip equations and are exploited in the bootstrap approach for the classical Ising model¹⁴.

To obtain Theorem 2.1, a natural idea (as was done for $T=0$) is to relax the set of $\beta$-KMS states and impose the condition (14) for a subset of observables $a, b$ supported on some small $\Lambda$. However, the main obstruction one is faced with is that even if $a, b$ are local observables, the expression $ae^{-\beta H}be^{\beta H}$ is, in general, not local, except for commuting Hamiltonians. As such, even though (14) forms a set of linear equations on $\rho$, they involve the expectation of $\rho$ on nonlocal observables. We circumvent this issue by using another characterization of thermal equilibrium states via so-called energy-entropy balance (EEB) inequalities¹⁵: this is an infinite set of scalar convex inequalities, each indexed by an operator $a$, which carve out the set of $\beta$-KMS states. On the one hand, these inequalities are better suited than the $\beta$-KMS condition because they preserve locality, i.e., they only require the expectation value of the state $\rho$ on a finite region around the support
of a. On the other hand, a drawback of these inequalities is that there are infinitely many of them, and unlike the inequalities (13), it is not possible to express them as a linear positive semidefinite constraint. A key ingredient to prove our theorem is to formulate a matrix generalization of such inequalities using the matrix-valued relative entropy function (3). We show that an infinite set of scalar energy-entropy balance inequalities can be compactly formulated by a single nonlinear convex matrix inequality of the form $D_{op}(A_{\rho}||B_{\rho}) \le \beta C_{\rho}$ for some suitable matrices $A_{\rho}, B_{\rho}, C_{\rho}$ that depend linearly on $\rho$ (see (7)). To the best of our knowledge, this is the first application of the matrix relative entropy function in the design of convex relaxations in quantum information. Optimization problems involving this function can be solved using interior-point methods¹⁶ or via semidefinite programming approximations¹⁷.

## Quantifying convergence speed

Theorems 2.1 and 2.2 do not provide quantitative guarantees on the convergence speed. As previously mentioned, one cannot hope to prove general fast convergence guarantees as even the special case where the observable $O$ corresponds to energy is unlikely to have a polynomial-time quantum algorithm, even when $D=1^4$.

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

**Theorem 2.4.** (Quantitative convergence rate). Consider a translation-invariant local Hamiltonian $H$ on $\Gamma = \mathbb{Z}^D$ with $D \le 2$. Assume furthermore that $H$ is commuting, i.e., $[h_X, h_Y] = 0$ for all $X, Y \subset \mathbb{Z}^D$. For $D=1$ and $\beta_1 = \infty$ or for $D=2$ and some $\beta_1 > 0$, we have for all $0 \le \beta < \beta_1$, and for any local observable $O$, $\langle O \rangle_{\beta}^{\min, \mathrm{TI}} = \langle O \rangle_{\beta}^{\mathrm{max, TI}} = \langle O \rangle_{\beta}^{\mathrm{TI}}$ and, using the same notation as in Theorem 2.1

$$ \max\{p_{\ell}^{\max} - \langle O \rangle_{\beta}^{\mathrm{TI}}, \langle O \rangle_{\beta}^{\mathrm{TI}} - p_{\ell}^{\min}\} \le c_1 \|O\| e^{-c_2 \ell} \quad (16) $$

for some constants $c_1, c_2 > 0$ depending on the dimension and the interaction.

This result addresses an open problem raised in¹⁸ about the speed of convex optimization hierarchies for the (classical) Ising model, in particular, whether exponential convergence holds away from criticality. Theorem 2.4 establishes such a statement for the more general class of commuting local Hamiltonians.

## Discussion

A special case of the problem considered in this work is when the observable $O$ is the energy. Certified algorithms do exist for the ground energy of local Hamiltonians: one can combine semidefinite programming relaxations¹⁹⁻²² for lower bounds and variational methods such as tensor networks²³⁻²⁵ for upper bounds. Using such two-sided bounds on the ground energy, the recent works²⁶,²⁷ obtain bounds for expectation values of local observables in the ground state, although no convergence guarantees are given. The approach based on imposing the additional constraint (13), which leads to convergence guarantees, was proposed independently and concurrently in²⁸, where it was derived as a special case of a method to strengthen semidefinite relaxations for noncommutative polynomial optimization problems by incorporating optimality conditions as additional constraints. These papers, however, do not discuss the case of positive temperature since the equilibrium states are not characterized by a noncommutative polynomial optimization problem which is linear in the state. One could try to use lower

bounds for the free energy coming from convex relaxations²⁹, but it is not clear how to use such bounds for observables. We note, however, that for classical systems, thermal observables can be obtained via the bootstrap approach¹⁴,³⁰, by directly imposing the KMS conditions (15). The resulting hierarchies were shown to be asymptotically convergent¹⁸, even if only a weaker set of constraints are imposed. For quantum systems, to the best of our knowledge, Theorem 2.1 provides the first certified algorithms for general observables in thermal states.

Note that, in some restricted settings, such as 1D gapped systems³¹ at zero temperature or 1D systems at positive temperature³²⁻³⁵, there are provably efficient algorithms computing representations of equilibrium states and thus expectation values, but such algorithms are tailored to these settings. In addition, for arbitrary dimensions and high temperatures, provably efficient algorithms for computing the free energy exist³⁶⁻³⁸, which can be used to compute observables³⁹.

The preliminary numerical experiments in the supplemental information demonstrate that the approach proposed in this paper is not only theoretical, and with additional efforts on the computational side, can play an important role alongside other classical algorithms. For example, because our algorithm produces certified bounds, it can be used to rigorously benchmark variational algorithms for quantum many-body systems⁴⁰. In addition, one can easily identify several directions for improving the accuracy of the algorithm presented here: First, one can exploit additional convex constraints that are known to hold for marginals of equilibrium states. For example, one can add entropy-based inequalities as proposed in ref. 41. Furthermore, other valid inequalities tailored to specific models can be added to the convex relaxation such as reflection positivity which was used in ref. 14 for the Ising model. Second, it would be very interesting to combine the methods from this paper with variational methods (e.g., from tensor networks) to obtain more accurate bounds, such as the recent work on the ground energy problem⁴². Another natural question is whether one can use our algorithms within hybrid classical-quantum algorithms for quantum simulation problems⁴³. On the analytical side, it remains open whether convergence guarantees for other classes of Hamiltonian can be proven. Promising candidates would be regimes in which other classical algorithms are efficient such as thermal states in 1D and at high-temperature (for general noncommuting Hamiltonians) as well as gapped ground states in 1D.

## Data availability

No datasets were generated or analysed during the current study.

## References

1. Osborne, T. J. Hamiltonian complexity. Rep. Progr. Phys. 75, 022001 (2012).

2. Gharibian, S. et al. Quantum Hamiltonian Complexity (Now Publishers, 2015).

3. Kempe, J., Kitaev, A. & Regev, O. The complexity of the local Hamiltonian problem. SIAM J. Comput. 35, 1070-1097 (2006).

4. Gottesman, D. & Irani, S. The quantum and classical complexity of translationally invariant tiling and Hamiltonian problems. In 2009 50th Annual IEEE Symposium on Foundations of Computer Science 95-104 (IEEE, 2009).

5. Ambainis, A. On physical problems that are slightly more difficult than QMA. In 2014 IEEE 29th Conference on Computational Complexity (CCC) 32-43 (IEEE, 2014).

6. Aharonov, D. & Irani, S. Hamiltonian complexity in the thermodynamic limit. In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing 750-763 (ACM, 2022).

7. Watson, J. D. & Cubitt, T. S. Computational complexity of the ground state energy density problem. In Proc. 54th Annual ACM SIGACT Symposium on Theory of Computing 764-775 (ACM, 2022).

8. Cubitt, T. S., Pérez-García, D. & Wolf, M. M. Undecidability of the spectral gap. Nature 528, 207-211 (2015).
9. Terhal, B. M. & DiVincenzo, D. P. Problem of equilibration and the computation of correlation functions on a quantum computer. Phys. Rev. A **61**, 022301 (2000).

10. Temme, K., Osborne, T. J., Vollbrecht, K. G., Poulin, D. & Verstraete, F. Quantum metropolis sampling. *Nature* **471**, 87–90 (2011).

11. Rall, P., Wang, C. & Wocjan, P. Thermal state preparation via rounding promises. *Quantum* **7**, 1132 (2023).

12. Chen, C.-F., Kastoryano, M. J., Brandão, F. G. & Gilyén, A. Quantum thermal state preparation. Preprint at arXiv:2303.18224 (2023).

13. Bausch, J., Cubitt, T. S. & Watson, J. D. Uncomputability of phase diagrams. *Nat. Commun.* **12**, 452 (2021).

14. Cho, M. et al. Bootstrapping the ising model on the lattice. Preprint at arXiv:2206.12538 (2022).

15. Araki, H. & Sewell, G. L. KMS conditions and local thermodynamical stability of quantum lattice systems. *Commun. Math. Phys.* **52**, 103–109 (1977).

16. Fawzi, H. & Saunderson, J. Optimal self-concordant barriers for quantum relative entropies. *SIAM J. Optim.* **33**, 2858–2884 (2023).

17. Fawzi, H., Saunderson, J. & Parrilo, P. A. Semidefinite approximations of the matrix logarithm. *Found. Comput. Math.* **19**, 259–296 (2019).

18. Cho, M. & Sun, X. Bootstrap, Markov Chain Monte Carlo, and LP/SDP hierarchy for the lattice Ising model. *J. High Energ. Phys.* **2023**, 47 (2023).

19. Mazziotti, D. A. Quantum chemistry without wave functions: two-electron reduced density matrices. *Account. Chem. Res.* **39**, 207–215 (2006).

20. Pironio, S., Navascués, M. & Acín, A. Convergent relaxations of polynomial optimization problems with noncommuting variables. *SIAM J. Optim.* **20**, 2157–2180 (2010).

21. Baumgratz, T. & Plenio, M. B. Lower bounds for ground states of condensed matter systems. *New J. Phys.* **14**, 023027 (2012).

22. Barthel. T. & Hübener, R. Solving condensed-matter ground-state problems by semidefinite relaxations. *Phys. Rev. Lett.* **108**, 200404 (2012).

23. Bridgeman, J. C. & Chubb, C. T. Hand-waving and interpretive dance: an introductory course on tensor networks. *J. Phys. A Math. Theor.* **50**, 223001 (2017).

24. Cirac, J. I., Pérez-García, D., Schuch, N. & Verstraete, F. Matrix product states and projected entangled pair states: concepts, symmetries, theorems. *Rev. Mod. Phys.* **93**, 045003 (2021).

25. White, S. R. Density matrix formulation for quantum renormalization groups. *Phys. Rev. Lett.* **69**, 2863–2866 (1992).

26. Wang, J. et al. Certifying Ground-State Properties of Many-Body Systems. *Phys. Rev. X* **14**, 031006 (2024).

27. Han, X. Quantum many-body bootstrap. Preprint at arXiv:2006.06002 (2020).

28. Araújo, M., Klep, I., Vértesi, T., Garner, A. J. P. & Navascués, M. Karush-Kuhn-Tucker conditions for non-commutative optimization problems. (2023).

29. Poulin, D. & Hastings, M. B. Markov entropy decomposition: a variational dual for quantum belief propagation. *Phys. Rev. Lett.* https://doi.org/10.1103/physrevlett.106.080403 (2011).

30. Poland, D., Rychkov, S. & Vichi, A. The conformal bootstrap: theory, numerical techniques, and applications. *Rev. Modern Phys.* **91**, 015002 (2019).

31. Landau, Z., Vazirani, U. & Vidick, T. A polynomial time algorithm for the ground state of one-dimensional gapped local Hamiltonians. *Nat. Phys.* **11**, 566–569 (2015).

32. Hastings, M. B. Solving gapped Hamiltonians locally. *Phys. Rev. B* https://doi.org/10.1103/physrevb.73.085115 (2006).

33. Molnar, A., Schuch, N., Verstraete, F. & Cirac, J. I. Approximating Gibbs states of local Hamiltonians efficiently with projected entangled pair states. *Phys. Rev. B* **91**, 045138 (2015).

34. Kuwahara, T., Alhambra, A. M. & Anshu, A. Improved thermal area law and quasilinear time algorithm for quantum Gibbs states. *Phys. Rev. X* https://doi.org/10.1103/physrevx.11.011047 (2021).

35. Fawzi, H., Fawzi, O. & Scalet, S. O. A subpolynomial-time algorithm for the free energy of one-dimensional quantum systems in the thermodynamic limit. *Quantum* **7**, 1011 (2023).

36. Kuwahara, T., Kato, K. & Brandão, F. G. Clustering of conditional mutual information for quantum gibbs states above a threshold temperature. *Phys. Rev. Lett.* **124**, 220601 (2020).

37. Harrow, A. W., Mehraban, S. & Soleimanifar, M. Classical algorithms, correlation decay, and complex zeros of partition functions of quantum many-body systems. In *Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing* (ACM, 2020).

38. Mann, R. L. & Helmuth, T. Efficient algorithms for approximating quantum partition functions. *J. Math. Phys.* https://doi.org/10.1063/5.0013689 (2021).

39. Bravyi, S., Chowdhury, A., Gosset, D. & Wocjan, P. Quantum Hamiltonian complexity in thermal equilibrium. *Nat. Phys.* **18**, 1367–1370 (2022).

40. Wu, D. et al. Variational benchmarks for quantum many-body problems. Preprint at arXiv:2302.04919 (2023).

41. Fawzi, H., Fawzi, O. & Scalet, S. O. Entropy constraints for ground energy optimization. *J. Math. Phys.* **65**, 032201 (2024).

42. Kull, I., Schuch, N., Dive, B. & Navascués, M. Lower bounding ground-state energies of local Hamiltonians through the renormalization group. *Phys. Rev. X* **14**, 021008 (2024).

43. Zhang, R., Wang, G. & Johnson, P. Computing ground state properties with early fault-tolerant quantum computers. *Quantum* **6**, 761 (2022).

## Acknowledgements

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

## Author contributions

H.F., O.F., and S.O.S. contributed equally to this work. H.F., O.F., and S.O.S. contributed to discussions, proofs, and preparation of the manuscript.

## Competing interests

The authors declare no competing interest.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-51592-3.

**Correspondence** and requests for materials should be addressed to Samuel O. Scalet.

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024