Quantum many-body dynamics · classical simulation

Pauli
propagationfor real-time correlators

To predict what a neutron-scattering experiment will see, you need how a material responds over time. That calculation is notoriously expensive. This paper evolves the question instead of the state — and keeps the cost in check with two ideas.

The strategy in four stages Read the paper

F1 The pipeline (after Fig. 1)

The full strategy. 1 Evolve the operator in the Heisenberg picture as a growing sum of Pauli strings, truncated to stay tractable. 2 A Clifford-augmented MPS supplies a cheap, accurate ground state to measure against. 3 A positivity-based fit super-resolves the short-time signal out to late times, dodging the Pauli explosion. 4 The result is the dynamical structure factor S(k, ω) — exactly what experiments measure.

Evolve the operator

Work in the Heisenberg picture: the observable carries the time dependence, so a single-site probe can be tracked without ever writing down the full wavefunction.

Two brakes on the growth

Truncation on coefficient and weight controls the short-time cost; a positivity super-resolution reaches the long-time spectrum from a handful of points.

Where it pays off

1D as a proving ground, but the real target is 2D and beyond, where the usual tensor-network methods run out of room — here, an 8×8 lattice.

What we actually want

Almost everything an experiment measures about a material — how it scatters neutrons, absorbs light, conducts a current — is some version of the same quantity: a correlation function. Poke the system at one place and time, ask how it answers somewhere else, a little later. Written down, it is

$$C(t) \;=\; \operatorname{Tr}\big[\rho\, A(t)\, B\big]$$

(1)

where $A$ and $B$ are the pokes — say, a spin flip at two sites — and $\rho$ is the state the system starts in. Its Fourier transform $C(\omega)$ is what a spectrometer actually plots. Correlation functions are, in the authors' phrase, the connective tissue between the theory on the page and the data in the lab.

The catch is the $A(t)$. Getting there means running the system's dynamics in real time, and real-time quantum dynamics is one of the hardest things we know how to ask a computer to do. The three standard tools each hit a wall:

Why real-time correlators are hard. Each classical workhorse fails for a different reason, and all three fail worst in more than one dimension.
Method	Where it breaks
Exact diagonalization	Memory: the Hilbert space doubles with every spin.
Real-time quantum Monte Carlo	The sign (phase) problem: contributions cancel catastrophically.
Tensor networks (DMRG/TDVP)	Entanglement growth: the bond dimension explodes past 1D.

So the usual move — evolve the state $|\psi(t)\rangle$ forward in time — is expensive precisely because the state gets complicated fast. This paper takes the other door.

Evolve the question, not the state

Quantum mechanics offers two equivalent bookkeeping schemes. In the Schrödinger picture, the state moves and the observable sits still. In the Heisenberg picture, you freeze the state and let the observable carry the time dependence instead:

$$A(t) \;=\; e^{iHt}\, A\, e^{-iHt}$$

(2)

Same physics, different labor. And for a local probe — a single spin flip — the Heisenberg side has a real advantage. An operator does not instantly smear across the whole system; it spreads outward at a finite speed, its reach bounded by a light cone (the Lieb–Robinson bound). For a while, $A(t)$ stays small. The trick is to work inside that window.

The reframing

Instead of tracking an ever-more-entangled wavefunction, track the operator $A(t)$ and write it in a convenient basis. If it stays modestly sized, the whole correlation function is cheap. The basis that makes this work is the one built from Pauli strings.

A Pauli string is just a choice, at every site, of one of four one-qubit operators — identity $I$, or $X$, $Y$, $Z$ — multiplied together, like $X_0\,Z_1\,I_2\,Y_3\cdots$. Any operator can be written as a weighted sum of them:

$$A(t) \;=\; \sum_i a_i(t)\, P_i$$

(3)

Pauli strings are a beautiful basis for this job. Any two of them either commute or anticommute; multiply two and you get a third (up to a phase); and all the algebra reduces to flipping bits. The number of non-identity sites in a string is its weight — a direct measure of how far the operator has spread.

Pauli propagation, one rotation at a time

Now evolve. Chop the time evolution $e^{-iHt}$ into a long sequence of tiny rotations, one per term in the Hamiltonian, each of the form $e^{-i\theta P_k}$ with a small angle $\theta$. Feeding a Pauli string $P$ through one such rotation does one of exactly two things:

If $P$ commutes with $P_k$, nothing happens. The string passes through untouched. No new terms, no cost.
If $P$ anticommutes with $P_k$, it splits in two. One copy stays as $P$; a second, brand-new string $P_kP$ is born alongside it, and the weight $\cos 2\theta$ / $\sin 2\theta$ is shared between them.

$$e^{i\theta P_k}\, P\, e^{-i\theta P_k} \;=\; \cos(2\theta)\,P \;+\; i\sin(2\theta)\,P_k P$$

(4)

A concrete case makes it tangible. Rotate the single-site operator $X$ by an angle about $Z$. Since $X$ and $Z$ anticommute, $X$ splits — and because $ZX = iY$, the newborn term is a $Y$:

$$e^{i\theta Z}\, X\, e^{-i\theta Z} \;=\; \cos(2\theta)\,X \;-\; \sin(2\theta)\,Y$$

One operator became two. Do this again and again, across every rotation in the sequence, and each anticommuting step branches a string into two. The result is a binary tree whose leaves are all the Pauli strings making up $A(t)$. That is Pauli propagation: no wavefunction anywhere, just an operator growing through its own space.

F2 Playground — watch a string branch

angle θ 8°

rotations 4

Leaves (Pauli strings)—grows like 2ᵈ

Kept above ε—the sin branches fade

Biggest branch weight—cos vs sin split

Each anticommuting rotation splits a string into a cos 2θ branch (it stays the same) and a smaller sin 2θ branch (a new, heavier string). For small angles the sin branches are faint, so the tree is dominated by a handful of thick paths. Turn θ up and the split evens out — the tree fills in and the cost explodes. The green count keeps only branches above a coefficient cutoff ε ≈ 10⁻³, which is exactly what makes the method affordable.

The catch, and two brakes

The tree is the whole story, good and bad. Left alone, the number of Pauli strings grows exponentially — after enough time it would fill the entire operator space, which is itself exponentially large. Nobody can store that. The method only works because two things can be thrown away without much harm:

Coefficient cutoff ε. After each step, discard any string whose coefficient has shrunk below ε. These are the dynamically negligible terms — the faint sin branches. Exact as ε → 0.
Weight cutoff w. Discard any string spread over more than w sites. This caps how far the operator is allowed to reach in operator space. Exact as w → N.

Both are honest knobs: tighten them and you converge to the exact answer; loosen them and you trade accuracy for speed. In the paper's 1D tests, a coefficient cutoff of $\varepsilon = 10^{-4}$ already reproduces the exact correlator; loosen it to $10^{-2}$ and the early oscillations are still right, but the signal soon peels away from the truth.

F3 The cost of honesty — string count vs cutoff

coefficient cutoff ε 1e−3

weight cutoff w 6

Peak Pauli strings—at the chosen cutoffs

Accuracy window—how long it stays exact

A schematic of the central trade-off. The curve is the number of retained Pauli strings against time; tightening either cutoff pushes the peak up (more strings, more cost) but buys a longer window in which the propagation stays faithful. Loose cutoffs are cheap but go wrong early; tight cutoffs are accurate but expensive. Schematic, illustrating the behavior reported for the N = 64 chain — not the raw counts.

This buys a good short-time signal. But a spectrum needs long times — so how do you get a sharp $C(\omega)$ out of a correlator you can only trust for a little while? That is the second idea, and it is the clever one.

Super-resolution from a handful of points

Two facts about real correlators rescue the calculation. First, they decay — the system leaks information to its surroundings, so $C(t)$ dies away, and a spectrum only ever needs the times before it does. The paper enforces this honestly by multiplying in a damping window $e^{-t/\tau}$, which keeps the Fourier transform anchored on the early, accurate data.

Second, and more powerfully: many correlators contain only a few frequencies. A recent result (by one of the authors) says that a correlator of the form $\operatorname{Tr}[\rho\,A^\dagger(t)A]$ is a positive function. Line its time samples up into a matrix and that matrix is positive semi-definite with Toeplitz structure — and a classical theorem, running back to Carathéodory, Toeplitz and Pisarenko, says such a matrix is fixed uniquely by just a few numbers. Concretely: if the signal is a sum of $r$ oscillations, then $2r$ clean samples pin down all of them.

The payoff

You do not need the long-time signal — you need to identify the frequencies, and a short accurate window is enough. Fit a few oscillations to the early data (a Prony-style / low-rank decomposition), then extend the correlator to any later time for free. In the 2D calculation, seven time points are stretched out to a full spectrum.

F4 Playground — reconstruct from short data

points measured 7

frequencies (rank r) 1

Samples needed (2r)—to fix the frequencies

Reconstruction—enough data?

The blue markers are the short, accurate window Pauli propagation actually computes; the solid curve is the full correlator recovered by fitting r oscillations and extending. Once you have 2r good samples the reconstruction locks on and runs out to any time you like. Ask for more frequencies than the data can support and the fit becomes underdetermined — the readout turns red. This is why momentum space is so friendly: there, each mode is usually a single sharp frequency, so r = 1 often suffices.

Where it fails, honestly: when the spectrum is a continuum rather than a few sharp lines — like the spinon excitations of the 1D antiferromagnet — there is no small $r$ to find, and the low-rank trick loses its grip. The method is at its best when the physics really is a handful of modes.

One more obstacle: the ground state

There is a subtlety unique to equilibrium correlators. Equation (1) asks for a trace against $\rho$ — usually the true ground state. If that ground state is a simple product state (as for the 1D ferromagnet), measuring a giant sum of Pauli strings against it is trivial. But for an antiferromagnet, or anything in 2D, the ground state is itself richly entangled, and evaluating those expectation values becomes the new bottleneck — sometimes a worse one than the propagation.

The fix is a recent tool called a Clifford-augmented matrix product state (CAMPS). Clifford operations are the "free" moves of quantum computing — they shuffle Pauli strings into other Pauli strings without ever increasing their number. CAMPS peels off a Clifford circuit that absorbs much of the entanglement, leaving a far cheaper (low-bond-dimension) matrix product state to represent the rest. Because the Clifford part acts directly on the Pauli operators for free, it is an unusually good match for Pauli propagation. That pairing is what lets the authors reach an $8\times 8$ antiferromagnet — the edge of what other methods can do.

What comes out

The end product is the dynamical structure factor $S(k,\omega)$ — the intensity of the system's response, sorted by momentum $k$ and frequency $\omega$. It is precisely the picture a neutron-scattering experiment builds up. Across the test cases, the recipe recovers the physics:

1D Heisenberg ferromagnet. Even loose truncations reproduce the smoothly dispersing single-magnon band, $\omega(k) = 4(1-\cos k)$. Momenta appear at the right frequencies; short-time data plus the positivity extension nails the spectrum.
1D Heisenberg antiferromagnet. A harder, entangled ground state. Pauli propagation captures the strongly dispersing mode and its spinon continuum near the band top — the honest place where the few-frequency assumption strains, and where tensor-network methods still hold the lead in 1D.
2D antiferromagnet, 8×8. The showcase. Seven time points, extended by a rank-1 fit, give a full $S(k,\omega)$ matching linear spin-wave theory: a Goldstone mode at $\Gamma$, a rise out to $X$, a softening at $M$. This is the regime where tensor networks struggle and the approach earns its keep.

F5 The dynamical structure factor

damping τ 1.5

Dispersion—peak frequency vs k

Linewidth ≈ 1/τ—set by damping

Schematic reconstruction of the extracted spectra. The bright ridge is where the response concentrates — the magnon dispersion ω(k) = 4(1 − cos k) for the 1D chain, and the linear-spin-wave curve ω(k) = A√(1 − γₖ²) across the 2D Brillouin zone, with its Goldstone zero at Γ and soft mode at M. The damping τ sets the linewidth: shorter windows blur the ridge (a Lorentzian of width 1/τ), longer ones sharpen it. Curves follow the functional forms reported in the paper; intensities are illustrative.

The honest summary the authors give: in 1D, mature tensor-network methods are hard to beat and this will not displace them. The interest is higher dimensions, where those methods lose efficiency and the combination here — Pauli-propagated dynamics, a CAMPS ground state, and positivity-based signal processing — opens a path to correlators on large 2D lattices that were previously out of reach. The operators' limited spreading even suggests far bigger systems are within range; a $100\times 100$ ferromagnet in the appendix recovers the full magnon spectrum.

Read it yourself

Efficient computation of real-time correlators using Pauli propagation — Alexander F. Kemper, Raghav G. Jha, Arnab Bachhar, Goksu C. Toga, Mariano Guerrero Perez, and Nicholas J. Mayhall (NC State University and Indiana University). arXiv:2607.24924 [quant-ph].

Open the PDF

Selected references

The load-bearing citations, from the paper's bibliography.

LR72
E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems. Commun. Math. Phys. 28 (1972) 251. The light cone that keeps operators local.
RLCK19
P. Rall, D. Liang, J. Cook and W. Kretschmer, Simulation of qubit quantum circuits via Pauli propagation. Phys. Rev. A 99 (2019) 062337. Where the method originates.
Kim23
Y. Kim et al., Evidence for the utility of quantum computing before fault tolerance. Nature 618 (2023) 500. The IBM Eagle experiment that sparked the classical race.
BC23
T. Begušić and G. K.-L. Chan, Fast classical simulation of evidence for the utility of quantum computing before fault tolerance (2023). arXiv:2306.16372.
BC25
T. Begušić and G. Kin-Lic Chan, Real-time operator evolution in two and three dimensions via sparse Pauli dynamics. PRX Quantum 6 (2025) 020302. arXiv:2409.03097.
RJTAH25
M. S. Rudolph, T. Jones, Y. Teng, A. Angrisani and Z. Holmes, Pauli propagation: a computational framework for simulating quantum systems (2025). arXiv:2505.21606.
SBJVM26
C. Shrikhande, A. Bachhar, A. R. Jimenez, E. F. Valeev and N. J. Mayhall, Pauli propagation enables fast classical simulation of strongly correlated quantum systems (2026). arXiv:2511.21651.
KYG24
A. F. Kemper, C. Yang and E. Gull, Denoising and extension of response functions in the time domain. Phys. Rev. Lett. 132 (2024) 160403. The positivity result behind the super-resolution.
Car11
C. Carathéodory, Über den Variabilitätsbereich der Fourier'schen Konstanten von positiven harmonischen Funktionen. Rend. Circ. Mat. Palermo 32 (1911) 193.
Pis73
V. F. Pisarenko, The retrieval of harmonics from a covariance function. Geophys. J. Int. 33 (1973) 347. Few-frequency retrieval, a.k.a. Prony.
QHQ24
X. Qian, J. Huang and M. Qin, Augmenting density matrix renormalization group with Clifford circuits. Phys. Rev. Lett. 133 (2024) 190402. CAMPS, the ground-state engine.
KKFK24
E. Kökcü, H. A. Labib, J. Freericks and A. F. Kemper, A linear response framework for quantum simulation of bosonic and fermionic correlation functions. Nat. Commun. 15 (2024) 3881.
KJBTGM26
A. F. Kemper, R. G. Jha, A. Bachhar, G. C. Toga, M. Guerrero Perez and N. J. Mayhall, Efficient computation of real-time correlators using Pauli propagation (2026). arXiv:2607.24924. The paper this page explains.