Polynomial-time quantum Gibbs sampling for the weak and strong coupling regime of the Fermi-Hubbard model at any temperature

Abstract

Quantum computers hold the potential to revolutionise the simulation of quantum many-body systems, with profound implications for fundamental physics and applications like molecular and material design. However, demonstrating quantum advantage in simulating quantum systems of practical relevance remains a significant challenge. In this work, we introduce a quantum algorithm for preparing Gibbs states of interacting fermions on a lattice with provable polynomial resource requirements. Our approach builds on recent progress in theoretical computer science that extends classical Markov chain Monte Carlo methods to the quantum domain. We derive a bound on the mixing time for quantum Gibbs state preparation by showing that the generator of the quantum Markovian evolution is gapped at any temperature up to a maximal interaction strength. This enables the efficient preparation of low-temperature states of weakly interacting fermions and the calculation of their free energy. We present exact numerical simulations for small system sizes that support our results and identify well-suited algorithmic choices for simulating the Fermi-Hubbard model beyond our rigorous guarantees.

While quantum computers have strong potential for quantum-many-body simulations, demonstrating an advantage for systems of practical relevance is still a challenge. Here, the authors show that quantum computers can efficiently sample thermal states of weakly and strongly interacting fermions – which is notoriously hard for classical Monte Carlo methods due to the fermionic sign problem.

Introduction

Quantum computers promise to have a transformative impact on computing, as quantum algorithms are believed to solve certain computational problems significantly faster—potentially offering super-polynomial speed-ups over their classical counterparts. However, it is important to recognise that such substantial quantum advantages are far from generic. Among the most promising areas of application is the simulation of quantum many-body systems¹. While quantum approaches to understanding the zero-temperature ground state physics of quantum many-body systems date back to the early days of quantum computing², recent progress in theoretical computer science has introduced new algorithmic techniques for preparing quantum Gibbs states at non-zero temperature. In particular, simulated Lindbladian evolution within the quantum circuit model has emerged as a powerful tool for this task³. These so-called quantum Gibbs samplers represent the non-commutative analogues of the classically successful Markov chain Monte Carlo (MCMC) methods⁴, and are expected to offer efficient access to non-zero temperature properties of quantum many-body systems, especially in regimes that are classically intractable.

Starting from the algorithmic breakthrough result³, there is the physical intuition that quantum Gibbs samplers can perform well for some specific case instances relevant to computational physics and computational chemistry. This is contrasted to known worst-case hardness results^{2,5, 6–7}, as well as to previously proposed intricate non-zero temperature quantum methods that are (partially) missing rigorous guarantees^{8, 9, 10, 11, 12, 13–14} or are believed to be computationally expensive on relevant finite-size instances for finite temperature^{15, 16, 17, 18–19}. The guiding idea behind the latest algorithmic Lindbladian constructions is to efficiently simulate thermalisation processes in nature, as, e.g., modelled by the Davies generator^20,21. In particular, it is possible to bring together algorithmic efficiency with an exact notion of quantum detailed balance, and we refer to the recent quantum Gibbs sampler frameworks^22,23 as well as references therein for an extended discussion.

Classical MCMC algorithms are termed efficient when they converge to the Gibbs state in time polynomial or even logarithmic in system size⁴. In contrast, the recently proposed quantum Gibbs samplers based on algorithmic Lindbladian evolution are solely efficient in the sense that all algorithmic steps are implemented efficiently, but the complexity or overall runtime crucially relies on the so-called mixing time, which specifies for how long we need to conduct the evolution to get ϵ-close to the Gibbs state. For an efficient algorithm, we would then require the Lindbladian to undergo fast mixing in time scaling polynomially with system size. This then has to be resolved for each physical system of interest on a case-by-case basis. The situation is akin to quantum adiabatic algorithms, where the runtime is governed by the Hamiltonian gap along the adiabatic path²⁴ and careful analytical studies are typically challenging.

With classical MCMC methods it is well-understood that rigorous bounds on the convergence time as studied in mathematical physics^{25, 26, 27, 28–29} often vastly overestimate the observed convergence times when running the algorithms in practice and reading out physical information, such as the partition function, relevant observables, or correlation functions. However, it is by design challenging to numerically run quantum Gibbs samplers to estimate practical mixing times—after all, we do not (yet) have reliable large-scale quantum computers and the classical simulation thereof is believed to be hard. As such, analytical bounds on mixing times are the first main tool to start our efficiency analysis. They would also serve as valuable benchmarks for future empirical studies of these models on fault-tolerant quantum computers, extending beyond the regimes with analytical guarantees.

In order to witness decisive quantum advantages, our goal is to provide quantum many-body systems with instances of quantum Gibbs samplers where

1. Classical methods are not sufficient to conclusively determine the physics at non-zero temperature.

2. We can derive rigorous and efficient bounds on the mixing time, at least for certain non-trivial parameter regimes.

3. We can subsequently give informed heuristics on how to fine-tune quantum Gibbs samplers to see efficient mixing times for relevant specific case parameter regimes beyond analytical worst-case guarantees.

As an example, previous general results on efficient quantum mixing times have been recently obtained in the high-temperature limit^30,31, where, however, classical methods are proven to also be efficient³². Other problem-specific abstract results on efficient mixing time bounds from the theoretical computer science literature include random sparse Hamiltonians³³, random local Hamiltonians³⁴, the toric code³⁵, and parent Hamiltonians of shallow quantum circuits^36,37.

In contrast, here we focus on the quantum simulation of fermionic systems of practical relevance—specifically the Fermi-Hubbard model on D-dimensional lattices—and establish rigorous bounds on its corresponding mixing times in both the weakly and strongly interacting regimes. This model enjoys widespread applications in science, for example, to the Mott metal-insulator transition and to high-temperature superconductivity³⁸, while at the same time remaining challenging for classical methods^39,40. As such, it further serves as a standard benchmark for computational methods⁴¹, including quantum simulations using ultra-cold atoms⁴². There is also evidence that, in contrast to quantumly hard, glassy spin systems, interacting fermionic systems are computationally challenging due to their entanglement structure⁴³—which is believed to be amenable to quantum methods (in regimes where classical ones might not be). Lastly, we emphasise again that the sought-after efficiency for specific cases and typical finite system sizes would not be in contradiction with the known worst-case hardness results for fermionic systems^6,7.

Results

Fermionic systems

We start with the recent algorithmic work³, which derives that the mixing time of Lindbladian quantum Gibbs samplers can be upper bounded by giving lower bounds on the spectral gap of the Lindbladian. This spectral gap, in turn, is equal to the spectral gap of a corresponding parent Hamiltonian, which we analyse. Our approach is to first study free fermions, governed by the quadratic Hamiltonian $$H=H_0={\sum }_{i,j=1}^{2n}{\omega }_i{h}_{ij}{\omega }_j,\text{with}\left\{{\omega }_i,{\omega }_j\right\}=2{\delta }_{ij},{\omega }_i^{†}={\omega }_i,$$ written in terms of Majorana fermions ω_i. Note that the single-particle Hamiltonian h includes the chemical potential terms, which effectively control the number of fermions in the system. Since our ensuing results allow for arbitrary chemical potentials, these methods can be used to study systems with desired levels of hole doping. We follow the general algorithmic quantum Gibbs sampler framework of ref. ²², which we review in detail in Supplementary Information, but the key point of the quantum algorithm is to simulate dynamics generated by an algorithmic Lindbladian of the form $${\mathcal{L}}^{†}\left[\rho \right]=-i\left[G,\rho \right]+{\sum }_{a\in \mathcal{A}}\left(L_a\rho L_a^{†}-\frac{1}{2}\left\{L_a^{†}{L}_a,\rho \right\}\right),$$ which is fully specified after choosing a suitable set of jump operators A^a and corresponding filter functions f^a. The Lindblad operators L_a are then filtered operator Fourier transforms of the jump operators, while the coherent term G is uniquely determined from the detailed-balance condition; this is the condition that guarantees the Gibbs state to be a fixed point of the Lindbladian evolution. Intuitively, these choices are akin to establishing the set of proposals and their acceptance probabilities in the Metropolis-Hastings algorithm. We make the specific design choice of Majorana jump operators, ${\left\{A^a={\omega }_a\right\}}_{a=1}^{2n}$ , and equal Gaussian filter functions f^a = f, obeying the required conditions. In particular, this allows us to exactly compute the finite spectral gap of the free fermionic parent Hamiltonian via Prosen’s third quantisation formalism⁴⁴ as (see Methods)

1

$${\mathrm{\Delta }}_0=2\cdot e^{-4{\beta }^2\parallel h{\parallel }^2}\cosh \left(2\beta \parallel h\parallel \right),$$ where β = T⁻¹ is the inverse temperature. This is lower bounded due to the locality of the system ensuring $\parallel h\parallel \le \mathcal{O}\left(1\right)$ . On top of this, recognising that the dynamics is restricted to that of Gaussian states if we start from a Gaussian state, we calculate the covariance matrix of the evolved state, which, together with optimal trace norm bounds⁴⁵ allows us to prove a logarithmic upper bound on the mixing time for free fermions $$t_{\mathrm{mix}}\le \frac{1}{2{\mathrm{\Delta }}_0}\log \left(\frac{2n}{ϵ}\right).$$ For the first time, this thus proves rapid mixing in logarithmic time of non-commuting Hamiltonians at any temperature.

Then, after observing that the free fermionic parent Hamiltonian in third quantisation itself simplifies to a (different) free fermionic Hamiltonian, we make use of Hastings’ stability result^{46, 47–48} on the spectral gap of free fermions under perturbation in order to quantitatively extend the finite spectral gap in the thermodynamic limit to the interacting parent Hamiltonian (Theorem 1). To lift the locality of interactions from the fermionic Hamiltonian to the Lindbladian’s parent Hamiltonian, we use Lieb-Robinson bounds and employ matrix analysis methods and identities such as Duhamel’s formula.

Our main finding is phrased most generally for lattice fermionic Hamiltonians with exponentially decaying interactions, henceforth termed quasi-local. Namely, we show that for such systems and at any constant temperature T > 0, there exists a constant maximal interaction strength below which the corresponding purified Gibbs states can be prepared efficiently. This result can be summarised in the following theorem and corollary:

Theorem 1

For any interacting quasi-local fermionic Hamiltonian H = H₀ + λV at any temperature T > 0, there exists a positive system-size-independent value ${\mathrm{\lambda }}_{\max }$ such that the Lindbladian ${\mathcal{L}}^{†}$ for Gibbs state preparation has a spectral gap Δ lower bounded by a constant for any $\mid \mathrm{\lambda }\mid \le {\mathrm{\lambda }}_{\max }$ , closing at most linearly in ∣λ∣ from that of the non-interacting case H₀.

Corollary 1.1

The mixing time of the Gibbs state preparation is upper bounded by $t_{\mathrm{mix}}=\mathcal{O}\left(n+\log \left(1/ϵ\right)\right)$ , and the overall quantum complexity to create the quantum Gibbs state is $\tilde{\mathcal{O}}\left(n^3\text{polylog}\left(1/ϵ\right)\right)$ and requires $\mathcal{O}\left(n\right)$ qubits; where n denotes the system size and ϵ > 0 the desired accuracy in trace norm.

Crucially, the constant on the interaction strength is independent of the system size, and thus we conclude that (weakly) interacting fermionic systems in fixed dimension and at any constant temperature can be efficiently simulated on quantum computers. Although free fermions are efficiently solvable, to the best of our knowledge, there is no provably efficient classical algorithm for the weakly-interacting regime, leading to a potential exponential quantum advantage. This is in contrast to the previously studied case of high temperatures^30,31, where they have shown stability around the infinite temperature limit, hence offering only a polynomial quantum advantage^32,49.

Our results apply, in particular, to the Fermi-Hubbard model. Its spinful version on a D-dimensional lattice is governed by the Hamiltonian $$H_{\mathrm{FH}}=-t{\sum }_{⟨i,j⟩,\sigma }\left(a_{i,\sigma }^{†}{a}_{j,\sigma }+a_{j,\sigma }^{†}{a}_{i,\sigma }\right)+U{\sum }_i{a}_{i,\uparrow }^{†}{a}_{i,\uparrow }{a}_{i,\downarrow }^{†}{a}_{i,\downarrow },$$ where ⟨ ⋅ , ⋅ ⟩ denotes neighbouring sites on the lattice, σ ∈ {↑, ↓} the spins, and $a_{i,\sigma }^{\left(†\right)}$ the fermionic annihilation and creation operators on site i with spin σ. The weakly-interacting limit corresponds to U/t ≲ 1^{39, 40–41}, which can then serve as an analytical starting point for further numerical investigations. Furthermore, while the Fermi-Hubbard model is exactly solvable for the D = 1 case⁵⁰, we emphasise that our results are equally valid in any dimension. We refer to ref. ³⁹ for a recent review of exact and heuristic results for the Hubbard model. In particular, the state-of-the-art in powerful heuristic methods, such as tensor networks, allows one to study the ground state of the Fermi-Hubbard model up to lattices of size 16 × 16⁵¹. Notably, in two-dimensional systems with next-nearest-neighbour hopping, even the weak-coupling regime exhibits rich behaviour as the chemical potential is varied, with multiple competing instabilities emerging near the van Hove singularity⁴⁰. Our algorithm enables the study of such phenomena with, in principle, arbitrarily high theoretical precision.

Applications, extensions, and numerical simulations

As an application of Gibbs sampling, we adapt ref. ³¹, [Theorem 8] to calculate partition functions for the considered systems. Hence, this provides the means to resolve the physics in thermal equilibrium of the underlying quantum many-body system in an end-to-end fashion:

Proposition 2

For any interacting quasi-local fermionic Hamiltonian H = H₀ + λV at any temperature T > 0, there exists a positive system-size-independent value ${\mathrm{\lambda }}_{\max }$ such that for any $\mid \mathrm{\lambda }\mid \le {\mathrm{\lambda }}_{\max }$ the corresponding partition function Z is determined in quantum complexity $\tilde{\mathcal{O}}\left(n^5{ϵ}^{-2}\right)$ , with success probability at least 3/4 and up to a relative error ϵ > 0.

Further, our methods equally apply to the opposite regime U/t ≫ 1 of the Fermi-Hubbard model. Here, we start by exactly solving the atomic limit t = 0 case of no hopping with the same choices of jump operators and filter functions (see Supplementary Information B3), after which we use adapted eigenvalue perturbation techniques to control finite t > 0. This result can be summarised in the following theorem:

Theorem 3

For any quasi-local perturbed atomic Hamiltonian H = H_atomic + λV at any temperature T > 0, there exists a positive system-size-independent value ${\mathrm{\lambda }}_{\max }$ such that the Lindbladian ${\mathcal{L}}^{†}$ for Gibbs state preparation has a spectral gap Δ lower bounded by a constant for any $\mid \mathrm{\lambda }\mid \le {\mathrm{\lambda }}_{\max }$ ; closing at most as $\mathcal{O}\left(\mid \mathrm{\lambda }{\mid }^{\alpha }\right)$ with arbitrary α < 1 from that of the atomic case H_atomic.

We again find that at any finite temperature, there exists a constant (system-size-independent) maximal t below which the corresponding Gibbs states can be algorithmically created with the complexities as stated in Corollary 1.1. We emphasise that the flexibility of our proof techniques naturally lends itself to future explorations of other quantum many-body systems in various regimes, such as spin systems in external fields, where the spin-spin interaction is treated as the perturbation. Here, one has to take into account that provably efficient classical algorithms exist for quantum perturbations of atomic Hamiltonians also at low temperatures⁵².

We perform small-scale exact classical simulations for the weakly-interacting spinless and spinful Fermi-Hubbard model in order to trial the hidden asymptotic constants in our analytical result. We find reasonable finite-gap behaviour and confirm, in particular, the predicted system size-independent scaling. Next, we aim to heuristically improve the algorithmic choices in order to shorten the mixing times beyond our theoretical guarantees. In general, the temperature dependence scales unfavourably in our analytical result, and this becomes at first equally visible in the numerical analysis. However, varying the choice of the jump operators from Majoranas to Paulis and the filter functions from Gaussian to Metropolis (see Supplementary Information Equation (A8)), we observe a spectral gap of the Lindbladian whose magnitude is no longer decreasing when lowering the temperature for the small system sizes we simulate. This would make the overall algorithmic dependence on the inverse temperature β only quasi-quadratic, significantly broadening the practical applicability of the algorithm.

To probe beyond the weakly-interacting regime covered by our analytical bounds, we investigate intermediate coupling strengths with 2 ≲  U/t ≲ 6 for spinless D = 1, 2 systems at different temperatures. Our numerical results suggest that, for one-dimensional and nearly one-dimensional cases, the algorithm remains efficient even at arbitrary coupling strength, incurring only a small additional cost that scales poly-logarithmically with U/t. For higher dimensions, however, the diversity of physical phenomena exceeds the scope of our simulations, currently preventing any definitive conclusions beyond our theoretical analysis. Details are provided in Figs. 1 and 2, and the corresponding simulation code is available at ref. ⁵³.

Fig. 1 Numerical results beyond our analytical guarantees. [Images not available. See PDF.]

Plotting the gap Δ of the full Lindbladian ${\mathcal{L}}^{†}$ associated with the spinless D = 2 Fermi-Hubbard thermal state with design choices beyond our analytical results—when using the Metropolis filter function and single site Pauli jump operators instead—as a function of the coupling strength U. Here we plot different system sizes separately, for the case β = 1 and t = 1, demonstrating a large spectral gap in the regime of intermediate coupling 2 ≲ U/t ≲ 6, which also does not seem to close with growing inverse temperature β (see Supplementary Figs. 7 and 9). At what coupling strength the gap closes is controlled by the support of the filter function, incurring only polylogarithmic additional algorithmic cost in U/t.

Fig. 2 Numerical results confirming our analytical findings. [Images not available. See PDF.]

Plotting the slope $\tilde{d}=\mp {\left(\frac{\partial \mathrm{\Delta }}{\partial U}\right.∣}_{U=0^{\pm }}$ under which the spectral gap Δ of the full Lindbladian ${\mathcal{L}}^{†}$ closes from that of the unperturbed Lindbladian ${\mathcal{L}}_0^{†}$ in the analytically bounded regime, as the system size n grows, at β = 1. As per our main result (Theorem 1), leading to the complexities as stated in Corollary 1.1, this quantity has to be upper-bounded uniformly in n. We refer to Supplementary Figs. for data for more sets of parameters exhibiting different types of behaviours.

Discussion

We have shown for the first time that the Gibbs states of weakly-interacting fermions in any dimension and at any constant temperature can be efficiently prepared in polynomial time on quantum computers. We have further demonstrated that this allows us to directly determine physical properties of the underlying systems—such as the partition function—in an end-to-end fashion while staying in polynomial quantum complexity. This is in contrast to, e.g., phase estimation-based quantum methods, which suffer from the so-called state preparation complexity bottleneck that generically hides large complexity overheads^1,41. As such, we believe that the presented methods have great potential for resolving the relevant non-zero-temperature physics of the Fermi-Hubbard and other fermionic models in classically challenging regimes, hence shedding light on the unknown parts of their phase diagrams. Other applications of efficient quantum Gibbs samplers to be explored are in quantum approaches to optimisation^54,55 and quantum machine learning⁵⁶.

Going forward, it would be interesting to improve the exact dependencies on all the relevant parameters in our main result, as well as to work out all the hidden constants therein. One might then also fine-tune the choice of jump operators and filter functions for specific systems and parameter regimes, hopefully even improving on the current cubic complexity in system size, potentially all the way down to quasi-linear. The rapid mixing of free fermions serves as a good indication for the weakly-interacting case to also mix rapidly, which would already bring down one factor of n. As we have seen, the cubic dependency on the system size is then mainly due to the quadratic dependence on the number of jump operators taken, which needs to be linear so that the Lindbladian dynamics can be irreducible and ergodic. However, we could potentially use different sets of jump operators at different times, consisting of as few as a single jump operator at any given time. The generated dynamics then could not be described by a single quantum Markov semigroup, and any one of the corresponding Lindbladians would not be able to get us from an arbitrary starting position to the Gibbs state; however, the combination of all of them could create a complicated path in the state space eventually getting us to the Gibbs state, with potentially significantly better dependency on the system size.

In order to determine the classical-quantum efficiency boundary for the Fermi-Hubbard model in the absence of reasonable quantum computers, it would be important to perform larger-scale classical simulations of the quantum Gibbs samplers, with varying design parameters to estimate the relevant spectral gap (e.g., based on tensor network methods^51,57,58). Especially for translationally invariant systems in D = 1, for which the parent Hamiltonian would also inherit the invariance, one could use the imaginary-time iTEBD algorithm to simulate the evolution or iDMRG to calculate the spectral gap for infinite-sized systems. Any such findings should then be compared to the state-of-the-art classical results^39,40 to make statements about large quantum advantages for practically relevant problems.

Finally, our presented proof methods based on eigenvalue perturbation techniques also seem promising to explore other quantum many-body systems in different regimes, including, for example, bosonic systems. As a first step, the extensions presented in Supplementary Information B3 are easily shown to hold for any Hamiltonians that are separable in the lattice sites.

Methods

Free fermions

Let’s consider a free fermionic Hamiltonian $H_0={\mathit{\omega }}^T\cdot h\cdot \mathit{\omega }={\sum }_{i,j=1}^{2n}{\omega }_i{h}_{ij}{\omega }_j$ given in terms of Majorana fermions, with h being a Hermitian and anti-symmetric matrix. We denote these operators using bold matrix-vector notation for convenience. We start by choosing the jump operators in the algorithmic Lindbladian to be Majorana fermions, ${\left\{A^a={\omega }_a\right\}}_{a=1}^{2n}$ , and the filter functions f^a to be equal and real in the Fourier space. Note that the conditions on the filter function then say that ${\widehat{f}}^a\left(\nu \right)=\widehat{f}\left(\nu \right)=q\left(\nu \right)\cdot e^{-\beta \nu /4}$ with q(ν) real and even.

We readily find the Heisenberg time-evolved jump operators as $$\mathbf{A}\left(t\right)=e^{i{H}_{0}t}\mathbf{A}{e}^{-i{H}_{0}t}=e^{i{H}_{0}t}\mathit{\omega }{e}^{-i{H}_{0}t}=e^{-4iht}\cdot \mathit{\omega },$$ and the Lindblad operators are then just $$\mathbf{L}={\int }_{-\infty }^{\infty }f\left(t\right)e^{-4iht}\mathrm{d}t\cdot \mathit{\omega }=\widehat{f}\left(-4h\right)\cdot \mathit{\omega }.$$ Hence we find that $${\sum }_{a\in \mathcal{A}}{L}_a^{†}{L}_a={\mathit{\omega }}^T\cdot q{\left(4h\right)}^2\cdot \sinh \left(2\beta h\right)\cdot \mathit{\omega }+\text{Tr}\left(\widehat{f}{\left(-4h\right)}^2\right),$$ but since [ω^T ⋅ B ⋅ ω, ω^T ⋅ C ⋅ ω] = ω^T ⋅ [B, C] ⋅ ω, this sum will commute with the Hamiltonian H₀. This, in turn, means it is invariant under time evolution, and so the coherent term G will be proportional to $\widehat{g}\left(0\right)$ , as $$G={\int }_{-\infty }^{\infty }g\left(t\right)\cdot {\sum }_{a\in \mathcal{A}}{L}_a^{†}{L}_a\mathrm{d}t\propto \widehat{g}\left(0\right).$$ But since $\widehat{g}\left(\nu \right)\propto \tanh \left(\beta \nu /4\right)$ , this means that the coherent term vanishes, G = 0.

For future convenience, before analysing the spectrum of the Lindbladian, let’s consider the similarity transformation ${\mathcal{H}}_0\left[\rho \right]={\sigma }_{\beta }^{-1/4}\cdot {\mathcal{L}}_0^{†}\left[{\sigma }_{\beta }^{1/4}\cdot \rho \cdot {\sigma }_{\beta }^{1/4}\right]\cdot {\sigma }_{\beta }^{-1/4}$ into the parent Hamiltonian. This superoperator is Hermitian due to the quantum detailed-balance condition. As this is a similarity transformation, the spectrum of ${\mathcal{H}}_0$ will be the same as of ${\mathcal{L}}_0^{†}$ . We can calculate that $${\sigma }_{\beta }^{-1/4}\mathbf{L}{\sigma }_{\beta }^{1/4}=\widehat{f}\left(-4h\right)\cdot e^{-\beta h}\cdot \mathit{\omega },$$ and the QDB condition also ensures ${\sigma }_{\beta }^{-1/4}{L}_a{\sigma }_{\beta }^{1/4}={\sigma }_{\beta }^{1/4}{L}_a^{†}{\sigma }_{\beta }^{-1/4}$ . Bringing these calculations together and following Prosen’s third quantisation⁴⁴, which we review in Supplementary Section A3, the parent Hamiltonian simplifies to $$\begin{array}{ccc}\hfill {\mathcal{H}}_0& \hfill \cong \hfill & -{\mathbf{c}}^{†}\cdot S\cdot \mathbf{c}+\mathbf{c}\cdot S\cdot {\mathbf{c}}^{†}+{\mathbf{c}}^{†}\cdot A\cdot {\mathbf{c}}^{†}+\mathbf{c}\cdot A\cdot \mathbf{c}-\text{Tr}\left(\sqrt{S^2+A^2}\right).\hfill \end{array}$$ Here we have restricted the Hilbert space to that of physical states with even numbers of Majorana fermions; and S = q(4h)², $A=q{\left(4h\right)}^2\sinh \left(2\beta h\right)$ , and ${\left\{c_i^{†},c_i\right\}}_{i=1}^{2n}$ is a set of 2n canonical fermionic creation and annihilation operators. This is just a quadratic fermionic operator, and hence its complete spectrum, which is the same as that of ${\mathcal{L}}_0^{†}$ , is straightforwardly calculable as $$\text{spec}\left({\mathcal{L}}_0^{†}\right)={\left\{{\sum }_{i=1}^{2n}\left(-1+{\left(-1\right)}^{x_i}\right)\cdot q{\left(4{ϵ}_i\right)}^2\cosh \left(2\beta {ϵ}_i\right)\right\}}_{x\in {\left\{0,1\right\}}^{2n}},$$ where ϵ_i ∈ spec(h). In particular, the corresponding spectral gap between the highest and second-highest eigenvalue is $${\mathrm{\Delta }}_0=2\cdot {\min }_{i}q{\left(4{ϵ}_i\right)}^2\cosh \left(2\beta {ϵ}_i\right).$$ Taking the Gaussian filter function, $q\left(\nu \right)=e^{-{\beta }^2{\nu }^2/8}$ , the spectral gap would be monotonically decaying with ∥h∥, and hence lower bounded by a constant for local systems obeying $\parallel h\parallel \le \mathcal{O}\left(1\right)$ . This argument also assures that the Gibbs state is the unique fixed point of the dynamics generated by ${\mathcal{L}}_0^{†}$ .

On top of this, by recognising that when we start with a Gaussian state ρ₀ and evolve it with a quadratic Lindbladian, we will stay within the subspace of Gaussian states, we can restrict our view to the evolution of the covariance matrix ${\mathrm{\Gamma }}_{ij}\left(t\right)=\frac{i}{2}\text{Tr}\left(\left[{\omega }_i,{\omega }_j\right]\rho \left(t\right)\right)$ . Denoting the initial covariance matrix by Γ₀, we can straightforwardly solve its equation of motion⁵⁹ and get that $\mathrm{\Gamma }\left(t\right)=\frac{i}{2}\tanh \left(2\beta h\right)+e^{-2q{\left(4h\right)}^2\cosh \left(2\beta h\right)\cdot t}\cdot \left({\mathrm{\Gamma }}_0-\frac{i}{2}\tanh \left(2\beta h\right)\right)\cdot e^{-2q{\left(4h\right)}^2\cosh \left(2\beta h\right)\cdot t}$ . We can also check that the covariance matrix of the Gibbs state σ_β is ${\mathrm{\Gamma }}_{{\sigma }_{\beta }}=\frac{i}{2}\tanh \left(2\beta h\right)=\mathrm{\Gamma }\left(\infty \right)$ , and so the evolution indeed converges to the Gibbs state.

Finally, we can use optimal trace norm bounds obtained in ref. ⁴⁵, which tell us that ${∥\rho \left(t\right)-{\sigma }_{\beta }∥}_{\mathrm{Tr}}\le \frac{1}{2}{\parallel \mathrm{\Gamma }\left(t\right)-{\mathrm{\Gamma }}_{{\sigma }_{\beta }}\parallel }_{\mathrm{Tr}}$ , which we can set smaller to ϵ and solve for t and hence deduce that $$t_{\mathrm{mix}}\le \frac{1}{4{\min }_{i}q{\left(4{ϵ}_i\right)}^2\cosh \left(2\beta {ϵ}_i\right)}\log \left(\frac{2n}{ϵ}\right).$$

Locality of parent Hamiltonians

Ref. ²² shows that if H is a geometrically local Hamiltonian, A_a are local, and the filter function is Gaussian, then the Lindblad operators L_a are quasi-local and G is a sum of quasi-local terms. Here, we extend this result and discuss the locality properties of the parent Hamiltonian for fermionic systems and systems with exponentially decaying interactions. The quasi-locality of the parent Hamiltonian will be an important ingredient in the proofs of gap stability we present below.

We shall again consider the (general) transformation into the parent Hamiltonian

4.1

$$\mathcal{H}\left[\rho \right]={\sigma }_{\beta }^{-1/4}\cdot {\mathcal{L}}^{†}\left[{\sigma }_{\beta }^{1/4}\cdot \rho \cdot {\sigma }_{\beta }^{1/4}\right]\cdot {\sigma }_{\beta }^{-1/4},$$ and also define the transformed operators appearing therein by ${\tilde{L}}_a={\sigma }_{\beta }^{-1/4}{L}_a{\sigma }_{\beta }^{1/4}$ and $\tilde{G}={\sigma }_{\beta }^{-1/4}G{\sigma }_{\beta }^{1/4}$ .

Proposition 4

Consider a Hamiltonian H with interactions that decay at least exponentially, and local jump operators A^a with Gaussian filter functions. Then the parent Hamiltonian (4.1) is a sum of quasi-local terms.

Proof

We define local approximations of ${\tilde{L}}_a$ and $\tilde{G}$ using $${\tilde{L}}_a^{\left(r\right)}={\int }_{-\infty }^{\infty }{f}^a\left(t+i\beta /4\right)e^{i{H}_{B_r\left(a\right)}t}{A}^a{e}^{-i{H}_{B_r\left(a\right)}t}\mathrm{d}t,{\tilde{G}}_a^{\left(r\right)}={\int }_{-\infty }^{\infty }g\left(t+i\beta /4\right)e^{i{H}_{B_r\left(a\right)}t}\left(L_a^{\left(r\right)†}{L}_a^{\left(r\right)}\right)e^{-i{H}_{B_r\left(a\right)}t}\mathrm{d}t,$$ where B_r(a) is a ball of radius r around the support of A^a and $H_{\mathrm{\Omega }}={\sum }_{I\mid I\cap \mathrm{\Omega }\ne \mathrm{\varnothing }}{h}_I$ is the truncated Hamiltonian to the region Ω. The conjugation by the Gibbs state translates into the shift of the functional arguments in the integrands, and is explained in Supplementary Information A2, D.

Here we shall use a weaker version of the Lieb-Robinson bound than the one for local systems from ref. ⁶⁰ [Lemma 5] used in ref. ²² [Prop. 20], which also holds for exponentially decaying Hamiltonian interactions, and tells us that

4.2

$$\begin{array}{c}∥e^{iHt}{A}^a{e}^{-iHt}-e^{i{H}_{B_r\left(a\right)}t}{A}^a{e}^{-i{H}_{B_r\left(a\right)}t}∥\hfill \\ \le \parallel A^a\parallel \min \left\{2,J{e}^{-\mu r}\left(e^{\mu v\mid t\mid }-1\right)\right\}\hfill \end{array}$$ for some constants J, v, and μ. From now on, we shall assume ∥A^a∥ ≤ 1. Using the Gaussian filter, it will follow that $$\begin{array}{c}\hfill \parallel {\tilde{L}}_a-{\tilde{L}}_a^{\left(r\right)}\parallel \le C{e}^{-\mu r},\\ \hfill \parallel {\tilde{G}}_a-{\tilde{G}}_a^{\left(r\right)}\parallel \le \tilde{C}{e}^{-\tilde{\mu }r};\end{array}$$ calculations of which are detailed in the Supplementary Information A4.

Note that the Lieb-Robinson bound, which follows from the bound on the commutator with the Hamiltonian terms, is true in our fermionic setting independently of whether A_a is even or odd in the number of fermions, since the constituent Hamiltonian terms are always even, so that (4.2) still holds. We refer to ref. ⁶¹ for more on Lieb-Robinson bounds and locality for fermions.

Bounds on the strength of the perturbation

In this section, we shall think of the parent Hamiltonian $\mathcal{H}$ of the interacting system as a perturbation of the free fermionic parent Hamiltonian ${\mathcal{H}}_0$ . The expectation is that for a small enough interaction strength λ of the system’s Hamiltonian H = H₀ + λV, the spectral gap of the parent Hamiltonian $\mathcal{H}$ should remain constant. To prove this idea, we will make use of results about stability of spectral gap under perturbation for free fermions^{46, 47–48}, but before that we would need to understand how the strength of the parent Hamiltonian perturbation $\mathcal{V}=\mathcal{H}-{\mathcal{H}}_0$ depends on λ.

We start by deriving the following Lemma from Duhamel’s formula (see Supplementary Information D):

Lemma 5

For any operator O, Hermitian operators H₀, V, and $\mathrm{\lambda },\alpha \in \mathbb{C}$ , we have $$\begin{array}{c}∥e^{\alpha \left(H_0+\mathrm{\lambda }V\right)}O{e}^{-\alpha \left(H_0+\mathrm{\lambda }V\right)}-e^{\alpha H_0}O{e}^{-\alpha H_0}∥\hfill \\ \le \mid \mathrm{\lambda }\mid \mid \alpha \mid {\max }_{s\in \left[0,1\right]}∥\left[V,e^{s\alpha H_0}O{e}^{-s\alpha H_0}\right]∥.\hfill \end{array}$$

As previously, we shall denote the operators appearing in the full parent Hamiltonian by ${\tilde{L}}_a$ and $\tilde{G}$ , and we shall also denote their corresponding counterparts appearing in the free fermionic parent Hamiltonian by ${\tilde{L}}_a^0$ and ${\tilde{G}}^0$ . Then to understand the strength of $\mathcal{V}=\mathcal{H}-{\mathcal{H}}_0$ , we will need to bound $\parallel {\tilde{L}}_a-{\tilde{L}}_a^0\parallel$ and $\parallel {\tilde{G}}_a-{\tilde{G}}_a^0\parallel$ .

Using our Lemma together with the exact solution for time evolution in the free fermionic case, we find that $$\begin{array}{c}\hfill ∥e^{H\left(\beta /4+it\right)}{A}^a{e}^{-H\left(\beta /4+it\right)}-e^{H_0\left(\beta /4+it\right)}{A}^a{e}^{H_0\left(\beta /4+it\right)}∥\\ \hfill \le c_1\mid \mathrm{\lambda }\mid \cdot \mid \beta /4+it\mid \cdot e^{c_2\mid \beta /4+it\mid },\end{array}$$ where we have assumed that $\parallel h{\parallel }_{\infty }\le \mathcal{O}\left(1\right)$ and that the Hamiltonian contains only terms with even numbers of fermions. Importantly, the constants appearing in this expression are independent of the system size.

Hence, it immediately follows that $$\parallel {\tilde{L}}_a-{\tilde{L}}_a^0\parallel \le \mid \mathrm{\lambda }\mid {\int }_{-\infty }^{\infty }\mid f^a\left(t\right)\mid c_1\mid \beta /4+it\mid e^{c_2\mid \beta /4+it\mid }\mathrm{d}t=c_3\mid \mathrm{\lambda }\mid ,$$ for some constant c₃ independent of the system size. As the dissipative part of the perturbation $\mathcal{V}$ is a sum over different products or tensor products of this or equivalent expressions, we can conclude that the strength of this part grows linearly in ∣λ∣, independently of the system size (see Supplementary Information B2 for details).

Now looking at the coherent term, we shall split it up into quasi-local contributions G = ∑_aG_a with $G_a={\int }_{-\infty }^{\infty }g\left(t\right)e^{iHt}{L}_a^{†}{L}_a{e}^{-iHt}\mathrm{d}t$ . Then we similarly need to bound $\parallel {\tilde{G}}_a-{\tilde{G}}_a^0\parallel$ , which we will split into a part depending on $∥L_a^{†}{L}_a-L_a^{0†}{L}_a^0∥$ and a part depending on $$∥e^{H\left(\beta /4+it\right)}{L}_a^{0†}{L}_a^0{e}^{-H\left(\beta /4+it\right)}-e^{H_0\left(\beta /4+it\right)}{L}_a^{0†}{L}_a^0{e}^{-H_0\left(\beta /4+it\right)}∥.$$ Here we can again use our Lemma to bound this second contribution by $$\mid \mathrm{\lambda }\mid \mid \beta /4+it\mid {\max }_{s\in \left[0,1\right]}∥\left[V,e^{s\left(\beta /4+it\right)H_0}{L}_a^{0†}{L}_a^0{e}^{-s\left(\beta /4+it\right)H_0}\right]∥.$$ Using the exact solution $L_a^0={\sum }_i\widehat{f}{\left(-4h\right)}_{ai}{\omega }_i$ , we can upper-bound this further like $$\begin{array}{c}\hfill {\max }_{s\in \left[0,1\right]}∥\left[V,e^{s\left(\beta /4+it\right)H_0}{L}_a^{0†}{L}_a^0{e}^{-s\left(\beta /4+it\right)H_0}\right]∥\\ \hfill \le 2c_2{e}^{c_3\beta /4}\cdot w_h\left(t\right)\cdot \parallel \widehat{f}\left(-4h\right){\parallel }_{\infty }^2,\end{array}$$ where w_h(t) is system-size-independent function growing subexponentially in t (as discussed in Supplementary Information E). Observe that the previous argument for bounding the conjugated expression $∥{\sigma }_{\beta }^{-1/4}{L}_a^{†}{L}_a{\sigma }_{\beta }^{1/4}-{\sigma }_{\beta ,0}^{-1/4}{L}_a^{0†}{L}_a^0{\sigma }_{\beta ,0}^{1/4}∥$ also shows that $∥L_a^{†}{L}_a-L_a^{0†}{L}_a^0∥\le c_4\mid \mathrm{\lambda }\mid$ . Finally, this means that $$\parallel {\tilde{G}}_a-{\tilde{G}}_a^0\parallel \le {\int }_{-\infty }^{\infty }\mid g\left(t\right)\mid \cdot \left(c_4\mid \mathrm{\lambda }\mid +\mid \mathrm{\lambda }\mid \mid \beta /4+it\mid c_5{w}_h\left(t\right)\right)\mathrm{d}t=c_6\mid \mathrm{\lambda }\mid ,$$ where the convergence is ensured by the decay bounds of g(t) obtained in ref. ²² [Lemma 30].

This proves that the strength of the perturbation of the parent Hamiltonian is upper bounded by a constant multiple of the strength of the perturbation of the system’s Hamiltonian, uniformly in system size, i.e., that ${\mathcal{V}}_a$ , where $\mathcal{V}={\sum }_{a\in \mathcal{A}}{\mathcal{V}}_a$ , is upper bounded like $\parallel {\mathcal{V}}_a\parallel \le c\mid \mathrm{\lambda }\mid$ . To match the locality definition we will need to use exactly, we would further make a standard argument by writing this perturbation as a telescoping sum over different radii (see Supplementary Information B2).

In Supplementary Lemma E.1, we also present a slightly weaker notion of this result, with the strength bounded by ∣λ∣^α for an arbitrary constant α < 1 for small enough ∣λ∣, which works for general Hamiltonians, and hence can be used to obtain our secondary result for perturbations around the atomic limit.

Constant gap and fast mixing

In our previous sections, we have proven that for a quasi-local interacting fermionic Hamiltonian H = H₀ + λV and with our algorithmic choices, the parent Hamiltonian ${\mathcal{H}}_0$ has [J, ν]-decay and the perturbation $\mathcal{V}$ has (c∣λ∣, μ)-decay as per Definitions 1 and 2 of ref. ⁴⁶ (see Supplementary Information B2 for review). Hence, we can finish our discussion of the stability of the spectral gap by using [ref. ⁴⁶, Corollary 1] to show that the gap of $\mathcal{H}$ closes at most linearly in ∣λ∣ from that of ${\mathcal{H}}_0$ – uniformly in system size—and hence is lower bounded by a constant as per our main Theorem 1.

Finally, this result allows us to bound the mixing time t_mix of the Lindbladian ${\mathcal{L}}^{†}$ , defined as the minimal time necessary to get us ϵ-close to the Gibbs state from an arbitrary initial position, and show fast mixing in polynomial time bounded by $$t_{\mathrm{mix}}\le \frac{\log \left(\frac{2}{ϵ}\parallel {\sigma }_{\beta }^{-1/2}\parallel \right)}{\mathrm{\Delta }}=\mathcal{O}\left(n+\log \left(1/ϵ\right)\right).$$ Together with the complexity analysis from ref. ²² [Theorem 34], this gives us the overall time complexity of the algorithm as stated in Corollary 1.1.

Calculating the partition function

As a possible application of the efficient Gibbs state preparation, we adapt the strategy from ref. ³¹ for calculating partition functions $Z_{\beta }\left({\mathrm{\lambda }}_i\right)=\text{Tr}\left(e^{-\beta H\left({\mathrm{\lambda }}_i\right)}\right)$ to the case of interacting fermionic systems, where we have denoted H(λ_i) = H₀ + λ_iV. Note that we can calculate the non-interacting partition function explicitly as $$Z_{\beta }\left(0\right)={\prod }_{i=1}^{n}2\cosh \left(2\beta {ϵ}_i\right),$$ where the product is taken only over one ϵ_i ∈ spec(h) from each symplectic pair  ± ϵ_i. By measuring the observable $e^{\beta H\left({\mathrm{\lambda }}_i\right)}{e}^{-\beta H\left({\mathrm{\lambda }}_{i+1}\right)}$ in the state σ_β(λ_i), we would obtain the ratio $\frac{Z_{\beta }\left({\mathrm{\lambda }}_{i+1}\right)}{Z_{\beta }\left({\mathrm{\lambda }}_i\right)}$ . Preparing this observable and the Gibbs state will require access to block encodings of H₀ and V, from which we get a block encoding for H(λ_i) via LCU, and hence a block encoding for the observable and the Hamiltonian simulation via QSVT. By choosing a schedule 0 = t₁ ≤ t₂ ≤  ⋯  ≤ t_l−1 ≤ t_l = ∣λ∣ and denoting ${\mathrm{\lambda }}_i=t_i\frac{\mathrm{\lambda }}{\mid \mathrm{\lambda }\mid }$ , we can calculate Z_β(λ) as a telescoping product $$Z_{\beta }\left(\mathrm{\lambda }\right)=Z_{\beta }\left(0\right){\prod }_{i=1}^{l-1}\frac{Z_{\beta }\left({\mathrm{\lambda }}_{i+1}\right)}{Z_{\beta }\left({\mathrm{\lambda }}_i\right)}.$$ We refer to ref. ³¹ [Appendix C] for the details of these calculations, the gist of which lies in choosing the schedule such that t_i+1 − t_i = Θ(n⁻¹), and so l = Θ(n) as λ = Θ(1). Then we would prepare the Gibbs states σ_β(λ_i) and measure the expectation values of the observables $e^{\beta H\left({\mathrm{\lambda }}_i\right)}{e}^{-\beta H\left({\mathrm{\lambda }}_{i+1}\right)}$ for each i ∈ [l−1] at least Θ(nϵ⁻²) times. Averaging these measurements and calculating the partition function using the telescoping product would yield the result in Proposition 2.

Acknowledgements

During the write-up of our manuscript, we have become aware that Yu Tong and Yongtao Zhan were independently working on the same topic⁶². We thank them for agreeing on a date for uploading to the arXiv together. All authors acknowledge support from the EPSRC Grant number EP/W032643/1. M.B. acknowledges funding by the European Research Council (ERC grant agreement no. 948139) and the Excellence Cluster Matter and Light for Quantum Computing (ML4Q). The authors acknowledge the use of the Imperial College Research Computing Service⁶³ to obtain some of the results in this work.

Author contributions

The initial directions were conceptualised by R.B., M.B., and Š.Š. The main calculations and proofs were developed by Š.Š. and R.B., with Š.Š. focusing on weakly-interacting fermions and R.B. on the atomic limit. R.M. and Š.Š. developed the numerical simulations. All authors contributed to conceptual discussions and the write-up of the manuscript.

Peer review

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.

Data availability

We have not analysed any datasets as our work proceeds within a theoretical and mathematical approach. The figures shown contain all the available data, which was generated using the provided code.

Code availability

Source code for simulating the spectral gap is available at https://github.com/Quantum-AI-Lab-ICL/Quantum-Gibbs-Sampling.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-65765-1.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Dalzell, A. M. et al. Quantum algorithms: a survey of applications and end-to-end complexities (Cambridge University Press, 2025).

2. Kitaev, A. Y., Shen, A. & Vyalyi, M. N. Classical and Quantum Computation (American Mathematical Society, 2002).

3. Chen, C-F; Kastoryano, M; Brandão, FGSL; Gilyén, A. Efficient quantum thermal simulation. Nature; 2025; 646, pp. 561-566.2025Natur.646.561C [DOI: https://dx.doi.org/10.1038/s41586-025-09583-x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/41094249][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC12527926]

4. Levin, D. A. & Peres, Y. Markov Chains and Mixing Times (American Mathematical Society, 2017).

5. Aharonov, D; Gottesman, D; Irani, S; Kempe, J. The power of quantum systems on a line. Commun. Math. Phys.; 2009; 287, pp. 41-65.2009CMaPh.287..41A2480741 [DOI: https://dx.doi.org/10.1007/s00220-008-0710-3]

6. Schuch, N; Verstraete, F. Computational complexity of interacting electrons and fundamental limitations of density functional theory. Nat. Phys.; 2009; 5, pp. 732-735. [DOI: https://dx.doi.org/10.1038/nphys1370]

7. O’Gorman, B; Irani, S; Whitfield, J; Fefferman, B. Intractability of electronic structure in a fixed basis. PRX Quantum; 2022; 3, 020322.2022PRXQ..3b0322O [DOI: https://dx.doi.org/10.1103/PRXQuantum.3.020322]

8. Temme, K; Osborne, TJ; Vollbrecht, KG; Poulin, D; Verstraete, F. Quantum Metropolis sampling. Nature; 2011; 471, pp. 87-90.2011Natur.471..87T [DOI: https://dx.doi.org/10.1038/nature09770] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21368829]

9. Yung, M-H; Aspuru-Guzik, A. A quantum–quantum Metropolis algorithm. Proc. Natl Acad. Sci.; 2012; 109, pp. 754-759.2012PNAS.109.754Y [DOI: https://dx.doi.org/10.1073/pnas.1111758109] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22215584][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3271878]

10. Shtanko, O. & Movassagh, R. Preparing thermal states on noiseless and noisy programmable quantum processors. arxivhttps://arxiv.org/abs/2112.14688 (2021).

11. Moussa, J. E. Low-depth quantum Metropolis algorithm. arxivhttps://arxiv.org/abs/1903.01451 (2022).

12. Rall, P; Wang, C; Wocjan, P. Thermal state preparation via rounding promises. Quantum; 2023; 7, 1132. [DOI: https://dx.doi.org/10.22331/q-2023-10-10-1132]

13. Wocjan, P; Temme, K. Szegedy walk unitaries for quantum maps. Commun. Math. Phys.; 2023; 402, pp. 3201-3231.2023CMaPh.402.3201W4630494 [DOI: https://dx.doi.org/10.1007/s00220-023-04797-4]

14. Jiang, J. & Irani, S. Quantum Metropolis sampling via weak measurement. arxivhttps://arxiv.org/abs/2406.16023 (2024).

15. Chowdhury, AN; Somma, RD. Quantum algorithms for Gibbs sampling and hitting-time estimation. Quantum Inf. Comput.; 2017; 17, pp. 0041-0064.3676655

16. van Apeldoorn, J; Gilyén, A; Gribling, S; de Wolf, R. Quantum SDP-Solvers: better upper and lower bounds. Quantum; 2020; 4, 230. [DOI: https://dx.doi.org/10.22331/q-2020-02-14-230]

17. van Apeldoorn, J. & Gilyén, A. Improvements in quantum SDP-Solving with applications. In: 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), vol. 132, 99:1–99:15 (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2019).

18. Gilyén, A., Su, Y., Low, G. H. & Wiebe, N. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019, 193–204 (Association for Computing Machinery, New York, NY, USA, 2019).

19. An, D., Childs, A. M. & Lin, L. Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters. arxivhttps://arxiv.org/abs/2312.03916 (2023).

20. Mozgunov, E; Lidar, D. Completely positive master equation for arbitrary driving and small level spacing. Quantum; 2020; 4, 227. [DOI: https://dx.doi.org/10.22331/q-2020-02-06-227]

21. Nathan, F; Rudner, MS. Universal Lindblad equation for open quantum systems. Phys. Rev. B; 2020; 102, 115109.2020PhRvB.102k5109N [DOI: https://dx.doi.org/10.1103/PhysRevB.102.115109]

22. Ding, Z; Li, B; Lin, L. Efficient quantum Gibbs samplers with Kubo–Martin–Schwinger detailed balance condition. Commun. Math. Phys.; 2025; 406, 67.2025CMaPh.406..67D4871911 [DOI: https://dx.doi.org/10.1007/s00220-025-05235-3]

23. Gilyén, A., Chen, C.-F., Doriguello, J. F. & Kastoryano, M. J. Quantum generalizations of Glauber and Metropolis dynamics. arxivhttps://arxiv.org/abs/2405.20322 (2024).

24. Albash, T; Lidar, DA. Adiabatic quantum computation. Rev. Mod. Phys.; 2018; 90, 015002.2018RvMP..90a5002A3788424 [DOI: https://dx.doi.org/10.1103/RevModPhys.90.015002]

25. Temme, K. & Kastoryano, M. J. How fast do stabilizer Hamiltonians thermalize? 1505.07811 (2015).

26. Kastoryano, MJ; Brandão, FGSL. Quantum Gibbs samplers: the commuting case. Commun. Math. Phys.; 2016; 344, pp. 915-957.2016CMaPh.344.915K3508165 [DOI: https://dx.doi.org/10.1007/s00220-016-2641-8]

27. Brandão, FGSL; Kastoryano, MJ. Finite correlation length implies efficient preparation of quantum thermal states. Commun. Math. Phys.; 2019; 365, pp. 1-16.2019CMaPh.365..1B3900825 [DOI: https://dx.doi.org/10.1007/s00220-018-3150-8]

28. Bardet, I et al. Entropy decay for Davies semigroups of a one dimensional quantum lattice. Commun. Math. Phys.; 2024; 405, 42.2024CMaPh.405..42B4703456 [DOI: https://dx.doi.org/10.1007/s00220-023-04869-5]

29. Kochanowski, J., Alhambra, A. M., Capel, A. & Rouzé, C. Rapid thermalization of dissipative many-body dynamics of commuting Hamiltonians. arxivhttps://arxiv.org/abs/2404.16780 (2024).

30. Rouzé, C., França, D. S. & Alhambra, A. M. Efficient thermalization and universal quantum computing with quantum Gibbs samplers. In Proceedings of the 57th Annual ACM Symposium on Theory of Computing, STOC ’25, 1488–1495 (Association for Computing Machinery, New York, NY, USA, 2025).

31. Rouzé, C., França, D. S. & Alhambra, Á. M. Optimal quantum algorithm for Gibbs state preparation. arxivhttps://arxiv.org/abs/2411.04885 (2024).

32. Bakshi, A., Liu, A., Moitra, A. & Tang, E. High-temperature Gibbs states are unentangled and efficiently preparable. arxivhttps://arxiv.org/abs/2403.16850 (2024).

33. Ramkumar, A. & Soleimanifar, M. Mixing time of quantum Gibbs sampling for random sparse Hamiltonians. arxivhttps://arxiv.org/abs/2411.04454 (2024).

34. Basso, J., Chen, C.-F. & Dalzell, A. M. Optimizing random local Hamiltonians by dissipation. arxivhttps://arxiv.org/abs/2411.02578 (2024).

35. Ding, Z., Li, B., Lin, L. & Zhang, R. Polynomial-time preparation of low-temperature Gibbs states for 2D toric code. arxivhttps://arxiv.org/abs/2410.01206 (2024).

36. Bergamaschi, T., Chen, C.-F. & Liu, Y. Quantum computational advantage with constant-temperature Gibbs sampling. arxivhttps://arxiv.org/abs/2404.14639 (2024).

37. Rajakumar, J. & Watson, J. D. Gibbs sampling gives quantum advantage at constant temperatures with O(1)-local Hamiltonians. arxivhttps://arxiv.org/abs/2408.01516 (2024).

38. Bärtschi, A. et al. Potential applications of quantum computing at Los Alamos National Laboratory. arxivhttps://arxiv.org/abs/2406.06625 (2024).

39. Arovas, DP; Berg, E; Kivelson, SA; Raghu, S. The Hubbard model. Annu. Rev. Condens. Matter Phys.; 2022; 13, pp. 239-274.2022ARCMP.13.239A [DOI: https://dx.doi.org/10.1146/annurev-conmatphys-031620-102024]

40. Qin, M; Schäfer, T; Andergassen, S; Corboz, P; Gull, E. The Hubbard model: a computational perspective. Annu. Rev. Condens. Matter Phys.; 2022; 13, pp. 275-302.2022ARCMP.13.275Q [DOI: https://dx.doi.org/10.1146/annurev-conmatphys-090921-033948]

41. Chan, G. K.-L. Quantum chemistry, classical heuristics, and quantum advantage. arxivhttps://arxiv.org/abs/2407.11235 (2024).

42. Gross, C; Bloch, I. Quantum simulations with ultracold atoms in optical lattices. Science; 2017; 357, pp. 995-1001.2017Sci..357.995G [DOI: https://dx.doi.org/10.1126/science.aal3837] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28883070]

43. Anschuetz, ER; Chen, C-F; Kiani, BT; King, R. Strongly interacting fermions are nontrivial yet nonglassy. Phys. Rev. Lett.; 2025; 135, 030602.2025PhRvL.135c0602A [DOI: https://dx.doi.org/10.1103/cbqf-d24r] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/40758023]

44. Prosen, T. Third quantization: a general method to solve master equations for quadratic open Fermi systems. N. J. Phys.; 2008; 10, 043026.2410814 [DOI: https://dx.doi.org/10.1088/1367-2630/10/4/043026]

45. Bittel, L., Mele, A. A., Eisert, J. & Leone, L. Optimal trace-distance bounds for free-fermionic states: testing and improved tomography. arxivhttps://arxiv.org/abs/2409.17953 (2025).

46. Hastings, M. B. The stability of free Fermi Hamiltonians. arxivhttps://arxiv.org/abs/1706.02270 (2017).

47. De Roeck, W; Salmhofer, M. Persistence of exponential decay and spectral gaps for interacting fermions. Commun. Math. Phys.; 2018; 365, pp. 773-796.3907957 [DOI: https://dx.doi.org/10.1007/s00220-018-3211-z]

48. Koma, T. Stability of the spectral gap for lattice fermions. arxivhttps://arxiv.org/abs/2005.04548 (2020).

49. Mann, R. L. & Helmuth, T. Efficient algorithms for approximating quantum partition functions. J. Math. Phys.62, https://arxiv.org/abs/2004.11568 (2021).

50. Essler, F. H., Frahm, H., Göhmann, F., Klümper, A. & Korepin, V. E. The one-dimensional Hubbard model (Cambridge University Press, 2005).

51. Liu, W.-Y., Zhai, H., Peng, R., Gu, Z.-C. & Chan, G.K.-L. Accurate simulation of the Hubbard model with finite fermionic projected entangled pair states. Phys. Rev. Lett.134, 256502 (2025).

52. Helmuth, T; Mann, RL. Efficient algorithms for approximating quantum partition functions at low temperature. Quantum; 2023; 7, 1155. [DOI: https://dx.doi.org/10.22331/q-2023-10-25-1155]

53. Šmíd, Š., Meister, R., Berta, M. & Bondesan, R. Polynomial-time quantum Gibbs sampling for the weak and strong coupling regime of the Fermi-Hubbard model at any temperature. Quantum-Gibbs-Samplinghttps://doi.org/10.5281/zenodo.17390410 (2025).

54. Krzakała, F; Montanari, A; Ricci-Tersenghi, F; Semerjian, G; Zdeborová, L. Gibbs states and the set of solutions of random constraint satisfaction problems. Proc. Natl. Acad. Sci.; 2007; 104, pp. 10318-10323.2007PNAS.10410318K2317690 [DOI: https://dx.doi.org/10.1073/pnas.0703685104] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17567754]

55. Somma, RD; Boixo, S; Barnum, H; Knill, E. Quantum simulations of classical annealing processes. Phys. Rev. Lett.; 2008; 101, 130504.2008PhRvL.101m0504S [DOI: https://dx.doi.org/10.1103/PhysRevLett.101.130504] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18851429]

56. Amin, M. H., Andriyash, E., Rolfe, J., Kulchytskyy, B. & Melko, R. Quantum Boltzmann machine. Phys. Rev. X8, 021050 (2018).

57. Biamonte, J. & Bergholm, V. Tensor networks in a nutshell. https://arxiv.org/abs/1708.00006 (2017).

58. Mortier, Q. et al. Fermionic tensor network methods. https://arxiv.org/abs/2404.14611 (2024).

59. Barthel, T; Zhang, Y. Solving quasi-free and quadratic Lindblad master equations for open fermionic and bosonic systems. J. Stat. Mech.; 2022; 2022, 113101.4535557 [DOI: https://dx.doi.org/10.1088/1742-5468/ac8e5c]

60. Haah, J; Hastings, MB; Kothari, R; Low, GH. Quantum algorithm for simulating real time evolution of lattice Hamiltonians. SIAM J. Comput.; 2021; 52, FOCS18–250–FOCS18–284.4679966

61. Nachtergaele, B; Sims, R; Young, A. Lieb-Robinson bounds, the spectral flow, and stability of the spectral gap for lattice fermion systems. Math. Probl. Quantum Phys.; 2018; 717, pp. 93-115.3869135 [DOI: https://dx.doi.org/10.1090/conm/717/14443]

62. Tong, Y; Zhan, Y. Fast mixing of weakly interacting fermionic systems at any temperature. PRX Quantum; 2025; 6, 030301.2025PRXQ..6c0301T [DOI: https://dx.doi.org/10.1103/h1dx-ps5p]

63. Harvey, M. Imperial College research computing service. https://www.re3data.org/repository/r3d100011965 (2017).