When do propagation algorithms for open quantum systems converge fast?

Contractivity, not locality alone, decides runtime

A propagation algorithm represents an observable as a sum of Pauli or Majorana strings and evolves it in short steps, truncating all monomials above a degree cutoff ℓ after each step. Truncation is what keeps the representation storable: an operator of degree at most ℓ has at most N^O(ℓ) coefficients. The difficulty is that coherent interactions can keep raising the degree of the strings, so the discarded tail may grow rather than shrink.

The paper's central mechanism is contraction in norms that attach exponentially larger weights to higher-degree monomials. Writing N for the degree operator, the basic estimate is that the weighted norm of the evolved observable does not exceed the weighted norm of the initial one. Because high-degree components carry a large weight, this forces their unweighted contribution to be small: for an initial observable of degree at most d and cutoff ℓ ≥ d, the discarded tail is bounded by e^(−a(ℓ−d)) times the initial norm. The truncation error therefore falls exponentially with the cutoff, and combining this tail bound with stability estimates and a short-time approximation error yields guarantees for the full algorithm.

This reframes the question in the title. Fast convergence is not a property of locality by itself but of whether the dynamics contracts strongly enough to beat the degree growth generated by local interactions. The paper establishes the required contraction using the Frobenius norm for dissipative Hamiltonian dynamics and weighted coefficient ℓ1 norms for thermal expectation values; when dissipation is weak or absent, a time-dependent exponential weight controls degree growth only over a finite interval.

Noise strong enough removes the time limit

The first application studies the evolution generated by L(A) = i[H₀ + uV, A] − λN(A), where H₀ is quadratic, V is a quartic interaction with bounded coefficients, and u and λ ≥ 0 set interaction and dissipation strengths. The Hamiltonian H = H₀ + uV is assumed Δ-sparse: each Majorana mode appears in at most Δ interaction terms. No geometric locality assumption is imposed on the interaction graph.

The three terms act differently on the Majorana expansion. Evolution under the quadratic H₀ preserves degree, the interaction uV can increase it, and the noise suppresses a degree-k monomial at rate λk. The contraction estimate quantifies this competition. When λ is sufficiently large compared with Δu, a fixed exponential weight gives a truncation-tail bound that is uniform in time. The threshold is λ_c(u) = 2Δu sinh(1); at or above it, t_max is infinite.

Below that threshold the weight must be allowed to decrease with time, and the guarantee holds only up to a finite horizon: t_max(u, λ) = [1/(8Δu cosh(1))] log(1 + Δu/(λ_c(u) − λ)). For any finite 0 ≤ t ≤ t_max(u, λ) and 0 < ε < 1, the algorithm produces a degree-ℓ approximation with Frobenius-norm error at most ε‖A‖_F using ℓ = d + O(log(1 + t/ε)) and time N^O(log((1+t)/ε)). The practical reading is that dissipation buys simulation time, and enough dissipation buys unlimited time.

Noiseless limit extends the horizon to 1/u

Setting λ = 0 isolates the effect of the interaction alone. In that limit the simulation horizon scales as 1/(Δu), which the paper describes as extending the logarithmic dependence on 1/u established in Facelli, Fawzi, and Fawzi (2026) and matching the scaling recently obtained by Zhao, Marvian and Tong (2026).

This is a statement about the efficiently simulable time scale, not about the accuracy of the algorithm at a fixed time. The degree cutoff still grows only logarithmically in (1 + t)/ε, and the runtime remains N^O(log((1+t)/ε)); what changes is how far in time the guarantee reaches. The comparison matters because it locates the paper's contribution relative to prior bounds on the same quantity: the horizon improves from log(1/u) to about 1/u as the interaction weakens.

The horizon is also where the method's boundary sits. For λ below λ_c(u), the guarantee is finite-time by construction, and the paper does not claim efficient simulation beyond t_max. The bound depends on Δ, u, and λ, so a denser interaction graph or a stronger interaction shortens the window.

Thermal observables through a pseudo-Lindbladian

The second application computes Tr(Aρ_{β,H}) for ρ_{β,H} = e^(−βH)/Tr(e^(−βH)) with H = H₀ + uV. Rather than simulating relaxation, the paper introduces an auxiliary evolution generated by a pseudo-Lindbladian K_β, inspired by earlier work. Its semigroup P_t = e^(tK_β) is unital and its adjoint fixes ρ_{β,H}, so if P_t A converges to a scalar multiple of the identity, that scalar must be Tr(Aρ_{β,H}). The algorithm propagates A and extracts the coefficient of the identity.

The auxiliary evolution need not be completely positive; its role is to provide a computational route to the expectation value. Contraction does double duty here: it controls the high-degree tail and it establishes exponential decay of the nonidentity component, proving that P_t A converges to Tr(Aρ_{β,H})I. The convergence speed is explicit: for an even degree-d observable X, |tr(P_t X) − Tr(ρ_{β,H} X)| ≤ e^(−t/4) e^(a₀d)‖X‖₁.

To obtain the contraction, the generator is split into a reference part associated with the quadratic H₀ and a perturbation from uV. The reference evolution contracts a suitable exponentially weighted coefficient ℓ1 norm, and a Dyson expansion bounds the perturbation, showing that contraction survives when the interaction is weak relative to the temperature. The result is a deterministic propagation algorithm for thermal expectation values.

Weak interactions push thermal computation to log(1/u)

For a Δ-sparse Hamiltonian H = H₀ + uV with V quartic and H₀ quadratic and geometrically local on a D-dimensional lattice, and for an even degree-d Majorana monomial A, the algorithm returns an estimate μ̃ with |μ̃ − Tr(Aρ_{β,H})| ≤ ε‖A‖_∞ in time (N + 1)^O(d + log(1/ε)) poly(d, 1/ε). The guarantee applies whenever uΔ e^(Cβ)(e^(2β) − 1) ≤ c, where c > 0 is an absolute constant and C > 0 depends only on D.

For fixed Δ, that condition permits inverse temperatures of order log(1/u) as u → 0. This is the thermal analogue of the dynamical result: weak interactions let the propagation algorithm reach beyond a fixed high-temperature regime. The restriction to even observables is described in the paper as a convenience only.

The comparison with prior thermal algorithms is instructive. A cluster-expansion approach gives a deterministic scheme for partition functions with runtime poly(N, 1/ε) at fixed β and bounded interaction degree and rank, also permitting inverse temperatures of order log(1/u), but it requires the noninteracting part H₀ to be on-site and does not allow hopping terms of the form γ_iγ_j for unpaired i, j. Allowing hopping is nontrivial because diagonalizing H₀ generally requires a nonlocal change of basis that can destroy the locality of V. A cumulant-expansion approach handles hopping and gives randomized algorithms with runtime poly(N, 1/ε) at fixed β, with a sufficient condition on a D-dimensional lattice of β ≲ u^(−1/(D+2)).

The paper is explicit that its runtime guarantees are weaker than those of the cumulant-expansion algorithm. What it offers instead is a simple, deterministic propagation algorithm that accommodates hopping in H₀ and requires no geometric locality of V. The lattice assumption in the thermal theorem can also be replaced by summability conditions similar to those used for the noninteracting Green's function on bounded-degree graphs.

Where the guarantees stop and what could improve

The bounds are sufficient conditions, not characterizations. In the dissipative setting, efficiency is guaranteed only when λ is large enough relative to Δu, or only up to t_max when it is not. In the thermal setting, the condition uΔ e^(Cβ)(e^(2β) − 1) ≤ c ties the admissible temperature to the interaction strength and the sparsity, and the runtime is polynomial in d and 1/ε but exponential in the degree-dependent exponent.

The truncation rule itself is a place where the paper expects improvement. The authors state that they believe the runtime guarantees can be improved by analyzing more selective truncation rules that exploit coefficient magnitudes or spatial structure, and that establishing guarantees for such refinements is a natural direction for future work. That is a concrete statement about the method rather than a generic caveat: the current analysis discards all monomials above a fixed degree, and a rule that keeps only the coefficients that matter could tighten the bounds.

The paper also notes that the general framework covers both spin and fermionic systems, while the main results are presented for fermionic systems with N Majorana operators for concreteness. The two applications therefore share one mechanism — weighted contraction controlling degree growth — but they differ in the norm used, the object computed, and the conditions under which the guarantee holds.

The threshold to remember is λ_c

The single most decision-relevant quantity in this paper is the noise threshold λ_c(u) = 2Δu sinh(1). Above it, the dissipative propagation algorithm has a truncation-tail bound uniform in time and t_max is infinite. Below it, the guarantee is confined to t_max(u, λ) = [1/(8Δu cosh(1))] log(1 + Δu/(λ_c(u) − λ)), which grows as λ approaches λ_c from below and collapses toward the noiseless horizon 1/(Δu) as λ → 0.

Read together with the thermal result, the paper's priority is clear: the interesting regime is weak interaction, where the noiseless horizon scales as 1/(Δu) rather than log(1/u), and where the thermal algorithm reaches inverse temperatures of order log(1/u) while still allowing hopping in H₀ and no geometric locality in V. The trade-off the reader should weigh is that this deterministic, hopping-friendly guarantee comes with runtime bounds weaker than those of the cumulant-expansion algorithm, and with a temperature condition that couples u, Δ, and β.

For anyone deciding whether a propagation algorithm is the right tool, the paper's own variables give the checklist: compare the noise strength λ to 2Δu sinh(1); if it is below, compute t_max before trusting the simulation; and for thermal observables, check whether uΔ e^(Cβ)(e^(2β) − 1) stays below the constant c. Those three comparisons, not the asymptotic statements alone, determine whether the fast-convergence guarantee applies to a given system.

Sources used in this answer