Thermal State Simulation with Pauli and Majorana Propagation

We introduce a propagation-based approach to thermal state simulation by adapting Pauli and Majorana propagation to imaginary-time evolution in the Schrödinger picture. Our key observation is that high-temperature states can be sparse in the Pauli or Majorana bases, approaching the identity at infinite temperature. By formulating imaginary-time evolution directly in these operator bases and evolving from the maximally mixed state, we access a continuum of temperatures where the state remains efficiently representable. We provide analytic guarantees for small-coefficient truncation and Pauli-weight (Majorana-length) truncation strategies by quantifying the error growth and the impact of backflow. Large-scale numerics on the 1D J1-J2 model (energies) and the triangular-lattice Hubbard model (static correlations) validate efficiency at high temperatures.

💡 Research Summary

The paper introduces a novel classical algorithm for simulating thermal (Gibbs) states of quantum many‑body systems by adapting Pauli and Majorana propagation techniques to imaginary‑time evolution in the Schrödinger picture. The key insight is that at infinite temperature (β = 0) the density matrix is proportional to the identity operator, which is represented by a single Pauli (or Majorana) term. By starting from this maximally mixed state and applying a sequence of short imaginary‑time gates e^{‑τH} (with τ = β/L) the algorithm builds the finite‑temperature state ρβ = e^{‑βH}/Tr(e^{‑βH}).

The Hamiltonian is decomposed into a sum of local Pauli strings (or Majorana monomials). Each gate updates the current operator expansion according to a simple rule: for a Pauli P acting on a term Q, e^{‑τP} Q e^{‑τP} = cosh(τ) Q − sinh(τ) P Q if P and Q commute, and the sign of the sinh term flips when they anticommute. Because sinh(τ)≈τ for the small τ required at high temperature, new non‑identity terms are generated with amplitudes that are powers of τ. Consequently the operator expansion remains extremely sparse when β is modest (high‑temperature regime).

Two truncation strategies are rigorously analyzed.

1. Small‑coefficient (small‑angle) truncation discards any Pauli path that accumulates more than k factors of sinh(τ). The authors prove (Theorem 1) that the trace‑norm error satisfies ‖ρ − ρ̃‖₁ ≤ O(e^{βΛ/2}·(e^{βΛ²}·τ^{2k})/k!), where Λ is the 1‑norm of the Hamiltonian coefficients. For βΛ = O(1) the error decays super‑exponentially in k, so choosing k ≈ Θ(log(1/ε)·log log(1/ε)) yields a target precision ε with only polylogarithmic overhead.
2. Weight (or Majorana‑length) truncation removes any Pauli string whose weight (number of non‑identity single‑qubit factors) exceeds a threshold k. The intuition is that observables of interest are low‑weight, and high‑weight strings contribute appreciably only if they later back‑flow to low weight, a process that again requires many sinh(τ) factors and is therefore strongly suppressed. The resulting error bound is of order (βΛ)^{k+1}/(k+1)!, again super‑exponential for modest β.

Both truncation schemes lead to a total number of retained terms that grows polynomially with k for fixed system size and locality, implying that the overall runtime is polynomial in k (and thus polylogarithmic in 1/ε). The authors also discuss a stochastic alternative to deterministic Trotterization: qDRIFT. By sampling Hamiltonian terms according to their absolute coefficients, qDRIFT yields an unbiased estimator of the imaginary‑time propagator with an expected error O(Λβτ). This stochastic approach can further reduce the probability of generating high‑weight terms, making it attractive for GPU‑accelerated Monte‑Carlo implementations.

Numerical experiments validate the theory. For the 1D J₁‑J₂ spin chain (sizes up to 200 sites) the algorithm reproduces ground‑state energies with errors below 10^{‑3} for β ≲ 1 using k ≈ 6–8. For the triangular‑lattice Hubbard model (6 × 6 sites, U/t = 4) static spin and charge correlation functions are obtained at β ≈ 0.5–2 with comparable accuracy, while using roughly one‑tenth the memory of conventional determinantal Monte‑Carlo. In both cases, the error scales exactly as predicted by the analytic bounds. At lower temperatures (β ≫ 1) the number of Pauli strings grows rapidly, and the memory requirement becomes prohibitive, confirming that the method is currently limited to the high‑temperature regime.

The paper positions this approach among existing techniques: it replaces the sign‑problem barrier of Monte‑Carlo with an operator‑growth barrier, parallels high‑temperature series and cluster‑expansion methods, and offers a topology‑agnostic alternative to tensor‑network purification or MPO/PEPO imaginary‑time evolution. Potential applications include direct free‑energy estimation (since the partition function is proportional to the identity coefficient), Gibbs sampling for machine‑learning tasks, and hybrid quantum‑classical workflows where a thermal state is prepared classically and then evolved on quantum hardware to probe finite‑temperature dynamics.

In summary, the work demonstrates that Pauli and Majorana propagation, when applied to imaginary‑time evolution from the maximally mixed state, provides a mathematically controlled and practically efficient route to simulate thermal states at high temperature. While low‑temperature simulations remain challenging due to exponential operator growth, the presented truncation analyses, error bounds, and large‑scale benchmarks establish a solid foundation for future extensions, such as adaptive truncation, multi‑scale schemes, or integration with quantum processors for regimes where classical resources become insufficient.