arXiv:2608.19448

Efficient Classical Simulation of Weakly Interacting Fermion Dynamics

Chu Zhao, Iman Marvian, Yu Tong

quant-phmath-phmath.NASubmitted Aug 19, 2026• Announced Aug 21, 2026• 0 citations

Abstract

We consider the task of simulating the real-time dynamics of weakly interacting fermionic systems. In particular, we focus on computing the expectation value of a local observable $A$ at time $t$. By analyzing the convergence of the perturbative expansion in the interaction strength $λ$ for the Heisenberg-picture observable, we propose a polynomial-time algorithm for estimating this expectation value in the weakly interacting regime $λ|t|^{2D+1}=\mathcal{O}(1)$, when the Hamiltonian is geometrically local on a $D$-dimensional lattice. Importantly, this condition is independent of the system size. If the goal is instead to approximate the time-evolved observable in normalized Frobenius norm, we extend the convergence regime to $λ|t|=\mathcal{O}(1)$ with quasi-polynomial runtime. When the non-interacting part exhibits Anderson localization, our polynomial-time algorithm can be extended up to $λ|t|=\mathcal{O}(1)$, modulo polylogarithmic factors. Our algorithm brings together ideas from continuous-time QMC, diagrammatic QMC, and Majorana Propagation, but with a new Heisenberg-picture operator-growth analysis that makes the sampling complexity rigorously controllable. This leads to provably efficient classical algorithms in regimes where the interaction is weak enough that the sampling variance remains bounded independently of system size. Together, these results identify broad regimes in which weak interactions, locality, and localization can be leveraged to make real-time fermionic dynamics classically tractable.

Related papers

Fast convergence of Majorana Propagation for weakly interacting fermions

Giorgio Facelli, Hamza Fawzi, Omar Fawzi

quant-pharXiv:2601.05226Jan 9, 20260 citations

Fast convergence of propagation algorithms for open quantum systems

Giorgio Facelli, Hamza Fawzi, Omar Fawzi +3 more

quant-pharXiv:2610.06769Oct 6, 20260 citations

Convergence of the Cumulant Expansion and Polynomial-Time Algorithm for Weakly Interacting Fermions

Hongrui Chen, Cambyse Rouzé, Jielun Chen +6 more

quant-pharXiv:2512.12010Dec 16, 20250 citations

Fast classical simulation algorithms for free-fermion dynamics with magic input

Jiwon Heo, Oliver Reardon-Smith, Michał Oszmaniec +2 more

quant-pharXiv:2609.40187Oct 1, 20260 citations

Majorana string simulation of nonequilibrium dynamics in two-dimensional lattice fermion systems

Matteo D'Anna, Jannes Nys, Juan Carrasquilla

cond-mat.quant-gasCross-listedarXiv:2511.02809Physical review. B./Physical review. BNov 5, 20250 citations

Polynomial Time Quantum Gibbs Sampling for Fermi-Hubbard Model at any Temperature

Štěpán Šmíd, Richard J. Meister, Mario Berta +1 more

arXiv:2501.01412Jan 2, 20250 citations

arXiv:2608.19448

Efficient Classical Simulation of Weakly Interacting Fermion Dynamics

Chu Zhao, Iman Marvian, Yu Tong

quant-phmath-phmath.NASubmitted Aug 19, 2026• Announced Aug 21, 2026• 0 citations

arXiv HTML PDF

Abstract

We consider the task of simulating the real-time dynamics of weakly interacting fermionic systems. In particular, we focus on computing the expectation value of a local observable

A

at time

t

. By analyzing the convergence of the perturbative expansion in the interaction strength

λ

for the Heisenberg-picture observable, we propose a polynomial-time algorithm for estimating this expectation value in the weakly interacting regime

λ|t|^{2D+1}=\mathcal{O}(1)

, when the Hamiltonian is geometrically local on a

D

-dimensional lattice. Importantly, this condition is independent of the system size. If the goal is instead to approximate the time-evolved observable in normalized Frobenius norm, we extend the convergence regime to

λ|t|=\mathcal{O}(1)

with quasi-polynomial runtime. When the non-interacting part exhibits Anderson localization, our polynomial-time algorithm can be extended up to

λ|t|=\mathcal{O}(1)

, modulo polylogarithmic factors. Our algorithm brings together ideas from continuous-time QMC, diagrammatic QMC, and Majorana Propagation, but with a new Heisenberg-picture operator-growth analysis that makes the sampling complexity rigorously controllable. This leads to provably efficient classical algorithms in regimes where the interaction is weak enough that the sampling variance remains bounded independently of system size. Together, these results identify broad regimes in which weak interactions, locality, and localization can be leveraged to make real-time fermionic dynamics classically tractable.

Related papers

Fast convergence of Majorana Propagation for weakly interacting fermions

Giorgio Facelli, Hamza Fawzi, Omar Fawzi

Simulating the time dynamics of an observable under Hamiltonian evolution is one of the most promising candidates for quantum advantage as we do not expect efficient classical algorithms for this problem except in restricted settings. Here, we introduce such a setting by showing that Majorana Propagation, a simple algorithm combining Trotter steps and truncations, efficiently finds a low-degree approximation of the time-evolved observable as soon as such an approximation exists. This provides the first provable guarantee about Majorana Propagation for Hamiltonian evolution. As an application of this result, we prove that Majorana Propagation can efficiently simulate the time dynamics of any sparse quartic Hamiltonian up to time $t_{\text{max}}(u)$ depending on the interaction strength $u$ . For a time horizon $t \leq t_{\text{max}}(u)$ , the runtime of the algorithm is $N^{O(\log(t/\varepsilon))}$ where $N$ is the number of Majorana modes and $\varepsilon$ is the error measured in the normalized Frobenius norm. Importantly, in the limit of small $u$ , $t_{\text{max}}(u)$ goes to $+\infty$ , formalizing the intuition that the algorithm is accurate at all times when the Hamiltonian is quadratic.Simulating the time dynamics of an observable under Hamiltonian evolution is one of the most promising candidates for quantum advantage as we do not expect efficient classical algorithms for this problem except in restricted settings. Here, we introduce such a setting by showing that Majorana Propagation, a simple algorithm combining Trotter steps and truncations, efficiently finds a low-degree approximation of the time-evolved observable as soon as such an approximation exists. This provides the first provable guarantee about Majorana Propagation for Hamiltonian evolution. As an application of this result, we prove that Majorana Propagation can efficiently simulate the time dynamics of any sparse quartic Hamiltonian up to time $t_{\text{max}}(u)$ depending on the interaction strength $u$. For a time horizon $t \leq t_{\text{max}}(u)$, the runtime of the algorithm is $N^{O(\log(t/\varepsilon))}$ where $N$ is the number of Majorana modes and $\varepsilon$ is the error measured in the normalized Frobenius norm. Importantly, in the limit of small $u$, $t_{\text{max}}(u)$ goes to $+\infty$, formalizing the intuition that the algorithm is accurate at all times when the Hamiltonian is quadratic.

quant-pharXiv:2601.05226Jan 9, 20260 citations

Fast convergence of propagation algorithms for open quantum systems

Giorgio Facelli, Hamza Fawzi, Omar Fawzi +3 more

We establish general conditions for the rapid convergence of propagation algorithms for computing expectation values of local observables in quantum many-body systems, covering both spin and fermionic systems. Our results apply to evolutions that are sufficiently contractive, providing a general mechanism by which contractivity controls the growth generated by local interactions and enables efficient classical simulation. As a first application, we consider noisy time evolution generated by local Hamiltonians on arbitrary interaction graphs. We show that the dynamics can be efficiently simulated when the noise strength $λ$ is sufficiently large compared to the degree of the interaction graph. More generally, for systems with an interaction strength $u$ , our bounds determine a time horizon $t_{\text{max}}$ as a function of the interaction strength $u$ , the noise strength $λ$ , and the interaction structure below which the algorithm is efficient. In the noiseless limit, $λ=0$ , our analysis extends the efficiently simulable time scale from $\log(1/u)$ proved in Facelli, Fawzi, and Fawzi (2026) to $t_{\text{max}} \sim 1/u$ matching the recent improvement by Zhao, Marvian and Tong (2026). As a second application, we apply the same framework to the computation of local observables in Gibbs states of local Hamiltonians. We consider a pseudo-Lindbladian whose steady state is the Gibbs state and show that, at sufficiently high temperature, its propagation algorithm converges rapidly. In particular, for weakly-interacting systems with interaction strength $u$ , our result applies up to inverse temperatures $β\sim\log(1/u)$ .We establish general conditions for the rapid convergence of propagation algorithms for computing expectation values of local observables in quantum many-body systems, covering both spin and fermionic systems. Our results apply to evolutions that are sufficiently contractive, providing a general mechanism by which contractivity controls the growth generated by local interactions and enables efficient classical simulation. As a first application, we consider noisy time evolution generated by local Hamiltonians on arbitrary interaction graphs. We show that the dynamics can be efficiently simulated when the noise strength $λ$ is sufficiently large compared to the degree of the interaction graph. More generally, for systems with an interaction strength $u$, our bounds determine a time horizon $t_{\text{max}}$ as a function of the interaction strength $u$, the noise strength $λ$, and the interaction structure below which the algorithm is efficient. In the noiseless limit, $λ=0$, our analysis extends the efficiently simulable time scale from $\log(1/u)$ proved in Facelli, Fawzi, and Fawzi (2026) to $t_{\text{max}} \sim 1/u$ matching the recent improvement by Zhao, Marvian and Tong (2026). As a second application, we apply the same framework to the computation of local observables in Gibbs states of local Hamiltonians. We consider a pseudo-Lindbladian whose steady state is the Gibbs state and show that, at sufficiently high temperature, its propagation algorithm converges rapidly. In particular, for weakly-interacting systems with interaction strength $u$, our result applies up to inverse temperatures $β\sim\log(1/u)$.

quant-pharXiv:2610.06769Oct 6, 20260 citations

Convergence of the Cumulant Expansion and Polynomial-Time Algorithm for Weakly Interacting Fermions

Hongrui Chen, Cambyse Rouzé, Jielun Chen +6 more

We propose a randomized algorithm to compute the log-partition function of weakly interacting fermions with polynomial runtime in both the system size and precision. Although weakly interacting fermionic systems are considered tractable for many computational methods such as the diagrammatic quantum Monte Carlo, a mathematically rigorous proof of polynomial runtime has been lacking. In this work we first extend the proof techniques developed in previous works for proving the convergence of the cumulant expansion in periodic systems to the non-periodic case. A key equation used to analyze the sum of connected Feynman diagrams, which we call the tree-determinant expansion, reveals an underlying tree structure in the summation. This enables us to design a new randomized algorithm to compute the log-partition function through importance sampling augmented by belief propagation. This approach differs from the traditional method based on Markov chain Monte Carlo, whose efficiency is hard to guarantee, and enables us to obtain a algorithm with provable polynomial runtime.We propose a randomized algorithm to compute the log-partition function of weakly interacting fermions with polynomial runtime in both the system size and precision. Although weakly interacting fermionic systems are considered tractable for many computational methods such as the diagrammatic quantum Monte Carlo, a mathematically rigorous proof of polynomial runtime has been lacking. In this work we first extend the proof techniques developed in previous works for proving the convergence of the cumulant expansion in periodic systems to the non-periodic case. A key equation used to analyze the sum of connected Feynman diagrams, which we call the tree-determinant expansion, reveals an underlying tree structure in the summation. This enables us to design a new randomized algorithm to compute the log-partition function through importance sampling augmented by belief propagation. This approach differs from the traditional method based on Markov chain Monte Carlo, whose efficiency is hard to guarantee, and enables us to obtain a algorithm with provable polynomial runtime.

quant-pharXiv:2512.12010Dec 16, 20250 citations

Fast classical simulation algorithms for free-fermion dynamics with magic input

Jiwon Heo, Oliver Reardon-Smith, Michał Oszmaniec +2 more

The classical cost of simulating free-fermion (matchgate) circuits depends on both the input state and the computational task. We develop classical algorithms for circuits with non-Gaussian inputs, addressing exact sampling and both additive-error and exact estimation of expectation values. For inputs consisting of $n$ copies of the four-mode magic state, we give exact sampling algorithms, generalizing the Clifford and Clifford algorithm for Boson sampling, for passive and active free-fermion dynamics with worst-case arithmetic cost $O(2^n)$ per sample, without any multiplicative polynomial prefactor. The passive algorithm enables an exact simulation of the non-interacting regime of a recent trapped-ion experiment. We also extend recently developed additive-error estimation of number correlators and related observables from passive to active dynamics with runtime polynomial in $n$ and the inverse additive-error, for input states formed by products of $n$ four-mode even parity states. Our estimator uses a Gaussian-state decomposition adapted to each sampled state to control its second moment. This includes estimation of individual output probabilities. Finally, we show that for even-parity product inputs composed of constant-size blocks, expectation values of Majorana monomials of logarithmic weight can be computed exactly in polynomial time, improving on prior Majorana propagation-style algorithms with quasi-polynomial runtime. The exact algorithm processes the input blocks by dynamic programming, sharing calculations across subsets of Majorana factors. Our suite of algorithms provide classical benchmarks for fermionic quantum simulations across a broad spectrum of computational tasks.The classical cost of simulating free-fermion (matchgate) circuits depends on both the input state and the computational task. We develop classical algorithms for circuits with non-Gaussian inputs, addressing exact sampling and both additive-error and exact estimation of expectation values. For inputs consisting of $n$ copies of the four-mode magic state, we give exact sampling algorithms, generalizing the Clifford and Clifford algorithm for Boson sampling, for passive and active free-fermion dynamics with worst-case arithmetic cost $O(2^n)$ per sample, without any multiplicative polynomial prefactor. The passive algorithm enables an exact simulation of the non-interacting regime of a recent trapped-ion experiment. We also extend recently developed additive-error estimation of number correlators and related observables from passive to active dynamics with runtime polynomial in $n$ and the inverse additive-error, for input states formed by products of $n$ four-mode even parity states. Our estimator uses a Gaussian-state decomposition adapted to each sampled state to control its second moment. This includes estimation of individual output probabilities. Finally, we show that for even-parity product inputs composed of constant-size blocks, expectation values of Majorana monomials of logarithmic weight can be computed exactly in polynomial time, improving on prior Majorana propagation-style algorithms with quasi-polynomial runtime. The exact algorithm processes the input blocks by dynamic programming, sharing calculations across subsets of Majorana factors. Our suite of algorithms provide classical benchmarks for fermionic quantum simulations across a broad spectrum of computational tasks.

quant-pharXiv:2609.40187Oct 1, 20260 citations

Majorana string simulation of nonequilibrium dynamics in two-dimensional lattice fermion systems

Matteo D'Anna, Jannes Nys, Juan Carrasquilla

The study of real-time dynamics of fermions remains one of the last frontiers beyond the reach of classical simulations and is key to our understanding of quantum behavior in chemistry and materials, with implications for quantum technology. Here we introduce a Heisenberg-picture algorithm that propagates observables expressed in a Majorana-string basis using a truncation scheme that preserves Trotter accuracy and aims at maintaining computational efficiency. The framework is exact for quadratic Hamiltonians -- remaining restricted to a fixed low-weight sector determined by the physical observable -- admits variational initial states, and can be extended to interacting regimes via systematically controlled truncations. We benchmark our approach on one- and two-dimensional Fermi-Hubbard quenches, comparing against tensor network methods (MPS and fPEPS) and recent experimental data. The method achieves high accuracy on timescales comparable to state-of-the-art variational techniques and experiments, demonstrating that controlled Majorana-string truncation is a practical tool for simulating two-dimensional fermionic dynamics.The study of real-time dynamics of fermions remains one of the last frontiers beyond the reach of classical simulations and is key to our understanding of quantum behavior in chemistry and materials, with implications for quantum technology. Here we introduce a Heisenberg-picture algorithm that propagates observables expressed in a Majorana-string basis using a truncation scheme that preserves Trotter accuracy and aims at maintaining computational efficiency. The framework is exact for quadratic Hamiltonians -- remaining restricted to a fixed low-weight sector determined by the physical observable -- admits variational initial states, and can be extended to interacting regimes via systematically controlled truncations. We benchmark our approach on one- and two-dimensional Fermi-Hubbard quenches, comparing against tensor network methods (MPS and fPEPS) and recent experimental data. The method achieves high accuracy on timescales comparable to state-of-the-art variational techniques and experiments, demonstrating that controlled Majorana-string truncation is a practical tool for simulating two-dimensional fermionic dynamics.

cond-mat.quant-gasCross-listedarXiv:2511.02809Physical review. B./Physical review. BNov 5, 20250 citations

Polynomial Time Quantum Gibbs Sampling for Fermi-Hubbard Model at any Temperature

Štěpán Šmíd, Richard J. Meister, Mario Berta +1 more

Recently, there have been several advancements in quantum algorithms for Gibbs sampling. These algorithms simulate the dynamics generated by an artificial Lindbladian, which is meticulously constructed to obey a detailed-balance condition with the Gibbs state of interest, ensuring it is a stationary point of the evolution, while simultaneously having efficiently implementable time steps. The overall complexity then depends primarily on the mixing time of the Lindbladian, which can vary drastically, but which has been previously bounded in the regime of high enough temperatures [Rouzé et al. arXiv:2403.12691 and arXiv:2411.04885]. In this work, we calculate the spectral gap of the Lindbladian for free fermions using third quantisation, and also prove a logarithmic bound on its mixing time by analysing corresponding covariance matrices. Then we prove a constant gap of the perturbed Lindbladian corresponding to interacting fermions up to some maximal coupling strength. This is achieved by using theorems about stability of the gap for lattice fermions. Our methods apply at any constant temperature and independently of the system size. The gap then provides an upper bound on the mixing time, and hence on the overall complexity of the quantum algorithm, proving that the purified Gibbs state of weakly interacting (quasi-)local fermionic systems of any dimension can be prepared in $\widetilde{\mathcal{O}} (n^3 \operatorname{polylog}(1/ε))$ time on $\mathcal{O}(n)$ qubits, where $n$ denotes the size of the system and $ε$ the desired accuracy. As an application, we explain how to calculate partition functions for the considered systems. We provide exact numerical simulations for small system sizes supporting the theory and also identify different suitable jump operators and filter functions for the sought-after regime of intermediate coupling in the Fermi-Hubbard model.Recently, there have been several advancements in quantum algorithms for Gibbs sampling. These algorithms simulate the dynamics generated by an artificial Lindbladian, which is meticulously constructed to obey a detailed-balance condition with the Gibbs state of interest, ensuring it is a stationary point of the evolution, while simultaneously having efficiently implementable time steps. The overall complexity then depends primarily on the mixing time of the Lindbladian, which can vary drastically, but which has been previously bounded in the regime of high enough temperatures [Rouzé et al. arXiv:2403.12691 and arXiv:2411.04885]. In this work, we calculate the spectral gap of the Lindbladian for free fermions using third quantisation, and also prove a logarithmic bound on its mixing time by analysing corresponding covariance matrices. Then we prove a constant gap of the perturbed Lindbladian corresponding to interacting fermions up to some maximal coupling strength. This is achieved by using theorems about stability of the gap for lattice fermions. Our methods apply at any constant temperature and independently of the system size. The gap then provides an upper bound on the mixing time, and hence on the overall complexity of the quantum algorithm, proving that the purified Gibbs state of weakly interacting (quasi-)local fermionic systems of any dimension can be prepared in $\widetilde{\mathcal{O}} (n^3 \operatorname{polylog}(1/ε))$ time on $\mathcal{O}(n)$ qubits, where $n$ denotes the size of the system and $ε$ the desired accuracy. As an application, we explain how to calculate partition functions for the considered systems. We provide exact numerical simulations for small system sizes supporting the theory and also identify different suitable jump operators and filter functions for the sought-after regime of intermediate coupling in the Fermi-Hubbard model.

arXiv:2501.01412Jan 2, 20250 citations

We consider the task of simulating the real-time dynamics of weakly interacting fermionic systems. In particular, we focus on computing the expectation value of a local observable

A

at time

t

. By analyzing the convergence of the perturbative expansion in the interaction strength

λ

for the Heisenberg-picture observable, we propose a polynomial-time algorithm for estimating this expectation value in the weakly interacting regime

λ|t|^{2D+1}=\mathcal{O}(1)

, when the Hamiltonian is geometrically local on a

D

λ|t|=\mathcal{O}(1)

with quasi-polynomial runtime. When the non-interacting part exhibits Anderson localization, our polynomial-time algorithm can be extended up to

λ|t|=\mathcal{O}(1)

t_{\text{max}}(u)

depending on the interaction strength

u

. For a time horizon

t \leq t_{\text{max}}(u)

, the runtime of the algorithm is

N^{O(\log(t/\varepsilon))}

where

N

is the number of Majorana modes and

\varepsilon

is the error measured in the normalized Frobenius norm. Importantly, in the limit of small

u

t_{\text{max}}(u)

goes to

+\infty

, formalizing the intuition that the algorithm is accurate at all times when the Hamiltonian is quadratic.

λ

is sufficiently large compared to the degree of the interaction graph. More generally, for systems with an interaction strength

u

, our bounds determine a time horizon

t_{\text{max}}

as a function of the interaction strength

u

, the noise strength

λ

, and the interaction structure below which the algorithm is efficient. In the noiseless limit,

λ=0

, our analysis extends the efficiently simulable time scale from

\log(1/u)

proved in Facelli, Fawzi, and Fawzi (2026) to

t_{\text{max}} \sim 1/u

matching the recent improvement by Zhao, Marvian and Tong (2026). As a second application, we apply the same framework to the computation of local observables in Gibbs states of local Hamiltonians. We consider a pseudo-Lindbladian whose steady state is the Gibbs state and show that, at sufficiently high temperature, its propagation algorithm converges rapidly. In particular, for weakly-interacting systems with interaction strength

u

, our result applies up to inverse temperatures

β\sim\log(1/u)

n

copies of the four-mode magic state, we give exact sampling algorithms, generalizing the Clifford and Clifford algorithm for Boson sampling, for passive and active free-fermion dynamics with worst-case arithmetic cost

O(2^n)

per sample, without any multiplicative polynomial prefactor. The passive algorithm enables an exact simulation of the non-interacting regime of a recent trapped-ion experiment. We also extend recently developed additive-error estimation of number correlators and related observables from passive to active dynamics with runtime polynomial in

n

and the inverse additive-error, for input states formed by products of

n

four-mode even parity states. Our estimator uses a Gaussian-state decomposition adapted to each sampled state to control its second moment. This includes estimation of individual output probabilities. Finally, we show that for even-parity product inputs composed of constant-size blocks, expectation values of Majorana monomials of logarithmic weight can be computed exactly in polynomial time, improving on prior Majorana propagation-style algorithms with quasi-polynomial runtime. The exact algorithm processes the input blocks by dynamic programming, sharing calculations across subsets of Majorana factors. Our suite of algorithms provide classical benchmarks for fermionic quantum simulations across a broad spectrum of computational tasks.

\widetilde{\mathcal{O}} (n^3 \operatorname{polylog}(1/ε))

time on

\mathcal{O}(n)

qubits, where

n

denotes the size of the system and

ε

the desired accuracy. As an application, we explain how to calculate partition functions for the considered systems. We provide exact numerical simulations for small system sizes supporting the theory and also identify different suitable jump operators and filter functions for the sought-after regime of intermediate coupling in the Fermi-Hubbard model.