Pauli Path Simulators Overview

• Pauli Path Simulators are quantum simulation methods that evolve observables in the Heisenberg picture by propagating Pauli strings.
• They employ both stochastic sampling and deterministic coefficient tracking to handle branching and noise in quantum circuits.
• PPS are applied to variational state preparation and Hamiltonian dynamics, supported by robust software frameworks and theoretical resource bounds.

Searching arXiv for relevant papers on Pauli Path Simulators and closely related Pauli propagation work. Pauli Path Simulators (PPS) are a family of quantum-simulation methods in which the primary computational object is not the evolving statevector or density matrix but the Heisenberg-evolved observable, represented in the Pauli basis and manipulated as a branching set of Pauli strings. In this framework, gates, channels, or Trotter steps either permute Pauli strings, damp them, or split them into short linear combinations, so that the target expectation value becomes a sum, truncation, or sampling over “Pauli paths” rather than a direct traversal of the full Hilbert space. Across the recent literature, PPS has been developed as a qubit-circuit simulation method based on Pauli propagation and stabilizer norm, as a truncation-based simulator for noisy Clifford+$T$ circuits, as a practical software framework for large-scale utility experiments, as a classical engine for variational state preparation, and as an observable-centric simulator for real-time spin dynamics whose cost is governed by operator complexity rather than entanglement (Rall et al., 2019, González-García et al., 2024, Rudolph et al., 27 May 2025, Lin et al., 2 Oct 2025, Shao et al., 25 Oct 2025).

1. Definition and formal viewpoint

The common formal viewpoint is Heisenberg-picture propagation. Instead of simulating $U\rho U^\dagger$, PPS evaluates expectation values by moving the observable backward through the circuit or channel and then overlapping the result with a simple initial state. A representative expression is

$$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$

with the observable expanded as a Pauli sum,

$O=\sum_{i=1}^{N_O} c_i P_i.$

At the level of elementary updates, the central mechanism is that a Pauli rotation is either non-branching or 2-branching: $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$ This branching rule is the operational origin of a “Pauli path”: each successive commutation or anticommutation choice defines a branch through a propagation tree (Rudolph et al., 27 May 2025).

In the noisy-circuit setting, the path language becomes explicit. For a circuit with layers $U_1,\dots,U_d$, a Pauli path is a sequence

$$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$

and the noisy expectation value can be grouped by total path weight,

$\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$

Here depolarizing noise suppresses high-weight paths exponentially, which makes truncation by path weight natural in regimes where only low-weight contributions survive (González-García et al., 2024).

The older Pauli-propagation literature also distinguishes two unbiased Monte Carlo viewpoints: Schrödinger propagation, which advances a Pauli proxy for the state forward, and Heisenberg propagation, which propagates a Pauli proxy for the observable backward. Both replace direct state simulation with randomized Pauli estimators, and both frame circuit evaluation as repeated local updates of Pauli labels and scalar coefficients (Rall et al., 2019).

2. Algorithmic structure and propagation mechanics

PPS implementations divide into two broad algorithmic styles. One style is stochastic: sample Pauli representatives from a distribution derived from Pauli overlaps and use Monte Carlo to estimate the observable. In the qubit formulation, the sampling gadget is

$\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$

with

$\mathbb E[\hat c(A)\hat \sigma(A)] = A.$

The quantity $U\rho U^\dagger$0, called the stabilizer norm, controls estimator magnitude and therefore sample complexity. Clifford gates are efficient because they map Pauli strings to Pauli strings and satisfy $U\rho U^\dagger$1, so propagation does not branch (Rall et al., 2019).

A second style is deterministic or semi-deterministic propagation with explicit coefficient tracking. In variational state-preparation work, the observable is written as

$U\rho U^\dagger$2

and for a gate $U\rho U^\dagger$3,

$U\rho U^\dagger$4

After each gate, Pauli terms with coefficients below a threshold $U\rho U^\dagger$5 are discarded, yielding a sparse approximation $U\rho U^\dagger$6. This thresholded branching rule is the operational core of “Pauli path” simulation in utility-scale variational training (Lin et al., 2 Oct 2025).

At the software and systems level, deterministic PPS is naturally formulated as a tree-search problem. Depth-first search has low memory but cannot merge equivalent strings that arise from different histories; breadth-first search is more memory intensive but supports coefficient merging. The practical recommendation is “merging-BFS”: after each gate, identical Pauli strings are merged by adding coefficients. This choice is reinforced by a bit-level encoding in which $U\rho U^\dagger$7 are mapped to $U\rho U^\dagger$8, so Pauli strings become unsigned integers and can be manipulated by shifts, masks, XOR, and population counts. On top of this encoding, hash tables are presented as a strong general-purpose container because they support amortized $U\rho U^\dagger$9 insertion, deletion, and lookup, while still enabling efficient merging of repeated Pauli strings (Rudolph et al., 27 May 2025).

3. Resource measures, truncation rules, and analytical guarantees

The central resource quantity is not unique across PPS variants. In Monte Carlo Pauli propagation, runtime is governed by products of local stabilizer norms $$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$0, rather than by explicit entanglement measures. This produces linear-time per-sample simulation for local circuits, but sample complexity can still grow exponentially if many components have $$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$1. The same formalism identifies a broad efficiently simulatable family of non-stabilizer input states, namely the hyper-octahedral states defined by $$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$2 (Rall et al., 2019).

For noisy Clifford+$$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$3 circuits, the principal control parameter is the competition between non-Clifford branching and depolarizing suppression. If a 2D circuit is $$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$4-sparse in $$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$5 gates, then the path-count bound

$$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$6

implies efficient truncation whenever

$$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$7

The informal theorem states that if a 2D Clifford+$$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$8 circuit is $$\langle O\rangle=\langle 0|U^\dagger O\,U|0\rangle=\langle 0|O(t)|0\rangle,$$9-sparse for some constant $O=\sum_{i=1}^{N_O} c_i P_i.$0, then for any

$O=\sum_{i=1}^{N_O} c_i P_i.$1

the noisy expectation value of any observable $O=\sum_{i=1}^{N_O} c_i P_i.$2 can be classically computed to precision $O=\sum_{i=1}^{N_O} c_i P_i.$3 in polynomial time. This extends earlier average-case Pauli-path guarantees to fixed structured circuits (González-García et al., 2024).

For Hamiltonian real-time dynamics, the dominant quantity is the Operator Stabilizer Rényi entropy (OSE),

$O=\sum_{i=1}^{N_O} c_i P_i.$4

where $O=\sum_{i=1}^{N_O} c_i P_i.$5 are the Pauli coefficients. In this formulation, one keeps only the $O=\sum_{i=1}^{N_O} c_i P_i.$6 largest-magnitude coefficients,

$O=\sum_{i=1}^{N_O} c_i P_i.$7

and the truncation error is controlled by the Top-$O=\sum_{i=1}^{N_O} c_i P_i.$8 tail weight through

$O=\sum_{i=1}^{N_O} c_i P_i.$9

The paper then bounds $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$0 in terms of OSE and gives the sufficient condition

$R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$1

for target accuracy $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$2. In the 1D XY Heisenberg case with $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$3, it further proves that the number of nonzero Pauli coefficients in $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$4 scales as $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$5, establishing compressibility of the Heisenberg-evolved observable. The stated total cost is $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$6 for $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$7 Trotter steps, $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$8 Hamiltonian terms, system size $R_G(\theta)[P] = \begin{cases} P, & [P,G]=0,\[4pt] \cos(\theta)\,P+\sin(\theta)\,P', & \text{otherwise}, \end{cases} \qquad P' = i[G,P]/2.$9, and truncation budget $U_1,\dots,U_d$0 (Shao et al., 25 Oct 2025).

4. Software, engineering, and practical workflow

The maturation of PPS has included full software stacks. A general computational framework is provided by PauliPropagation.jl, which exposes PauliString and PauliSum types, a propagate() routine for evolving observables through circuits, overlap routines such as overlapwithzero(), overlapwithplus(), overlapwithcomputational(), and overlapwithpaulisum(), and circuit/topology helpers such as bricklayertopology(), rectangletopology(), and tfitrottercircuit(). The package also exposes multiple truncation controls, including min_abs_coeff, max_weight, and customtruncfunc, and is designed to be differentiable through Julia’s automatic differentiation ecosystem (Rudolph et al., 27 May 2025).

A separate line of work addresses the empirical question of whether a given PPS run is feasible before attempting an expensive fine-threshold simulation. The proposed protocol runs a circuit at a few coarse thresholds $U_1,\dots,U_d$1 and exploits an empirical power-law model for coefficient magnitudes,

$U_1,\dots,U_d$2

From this, the number of tracked Paulis and the runtime are extrapolated as functions of $U_1,\dots,U_d$3. The reported practical guidance is that starting around $U_1,\dots,U_d$4 and stepping by $U_1,\dots,U_d$5 already gives good extrapolations; one example achieved a memory prediction in under a minute on a 6-core CPU with less than $U_1,\dots,U_d$6 error; and for a 2D Ising benchmark the method estimated roughly $U_1,\dots,U_d$7 billion tracked terms and $U_1,\dots,U_d$8 hours runtime, while the estimate itself took under 7 minutes on one CPU core. The same work presents the BlueQubit SDK as an implementation vehicle, using device="pauli-path" and the option $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$6 for threshold sweeps (Gharibyan et al., 14 Jul 2025).

5. Applications to variational circuits, many-body dynamics, and noisy optimization

One major PPS application is classical training of parametrized circuits for quantum state preparation. In this setting, the variational objective

$U_1,\dots,U_d$9

is evaluated classically via PPS, with $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$0 for ground-state preparation. The method is applied to the 1D quantum Ising model, the 2D quantum Ising model on square and heavy-hex lattices, and the Kitaev honeycomb model at system sizes of one hundred qubits or more, specifically $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$1 for the 1D Ising, 2D square-lattice Ising, and Kitaev honeycomb cases, and $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$2 for the heavy-hex Ising case. Validation against exact solutions or DMRG gives relative energy error below $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$3 in the 1D transverse-field Ising model, below $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$4 in the gapped tilted-field case against DMRG, and below $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$5 in the gapped A phase of the $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$6 Kitaev honeycomb model for the best depths. On Quantinuum’s System Model H2-2, a PPS-trained $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$7 Kitaev circuit achieved hardware relative energy errors of $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$8 and $$s=(s_0,s_1,\dots,s_d), \qquad s_i \in \{I,X,Y,Z\}^{\otimes n},$$9, summarized as roughly $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$0, without error mitigation, and supported interferometric braiding of Abelian anyons beyond fixed-point models (Lin et al., 2 Oct 2025).

A second major application is Hamiltonian real-time dynamics. For the XXZ Heisenberg chain benchmark on a length-$\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$1 chain with open boundaries, initialized in the Néel state and measured through staggered magnetization $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$2, PPS shows high accuracy with very small $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$3 in the free regime $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$4, with good results already around $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$5. In the interacting regime $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$6, operator complexity grows faster and larger $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$7 is needed, but Pauli propagation remains competitive with TDVP; the comparison is described as reasonably balanced because an MPS with bond dimension $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$8 and PPS with $\langle O\rangle_{\mathrm{err}}=\sum_w F_w(1-p)^w.$9 have comparable parameter counts. The same study reports that OSE increases with time, OSE increases with $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$0, and the number of significant Pauli words grows roughly linearly in time for $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$1 but more rapidly, with a quadratic trend, for $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$2 (Shao et al., 25 Oct 2025).

PPS has also been used to analyze noisy variational circuits beyond the average case. In particular, for 2D QAOA-like circuits implementing Ising energies

$\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$3

the structured-circuit analysis shows that when geometrically non-local Ising instances are embedded into 2D using SWAP networks, constant depolarizing noise makes the noisy output energy efficiently classically approximable. The upshot is that such noisy 2D QAOA implementations on hard non-local instances do not provide asymptotic quantum advantage over classical algorithms under the stated noise model (González-García et al., 2024).

6. Hybrid and spatiotemporal extensions

PPS has also become the basis for hybrid classical-quantum methods. Quantum Enhanced Pauli Propagation (QuEPP) is explicitly built on approximate classical Pauli-path simulation and Clifford perturbation theory (CPT). Low-order CPT terms give a truncated classical approximation

$\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$4

where each $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$5 is a classically simulable Clifford circuit arising from a Pauli path. The hybrid step is to run the target circuit and the CPT ensemble on hardware, infer a global rescaling factor from the noisy ensemble, and use it to correct the omitted higher-order remainder. On IBM Heron, QuEPP is demonstrated on 2D random mirror circuits of up to 49 qubits, on a 32-qubit harder random mirror instance with depth 80, and on 10-qubit Trotterized Hamiltonian evolution, with consistent improvement over both raw noisy execution and truncated classical CPT (Majumder et al., 15 Mar 2026).

The notion of a Pauli path has also been generalized from circuit branches to full spatiotemporal fault trajectories. Spatiotemporal Pauli Processes (SPPs) are defined by applying a multi-time Pauli twirl to a process tensor, yielding a convex mixture over Pauli trajectories in spacetime,

$\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$6

This construction turns non-Markovian open-system dynamics into a classical joint distribution over Pauli trajectories while retaining temporal memory and spatial correlations. The resulting objects inherit tensor-network structure, with bond dimensions bounded by the environment’s Liouville-space dimension and, for finite environments, $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$7. The framework is demonstrated in surface-code simulations up to distance $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$8, including a temporally correlated “storm” model and a 2D quantum cellular automaton bath that exhibits a pseudo-critical regime with macroscopic error avalanches and complete breakdown of the usual surface-code distance scaling (Kam et al., 5 Mar 2026).

7. Limitations, failure modes, and terminological issues

PPS is not a uniform complexity class but a family of approximations whose efficiency depends on structural features of the circuit, the observable, and the noise model. In Hamiltonian dynamics, efficiency depends on low OSE and sparse Pauli support; strongly interacting dynamics or highly nonlocal observables can make the Pauli expansion dense; and Trotterization introduces a second source of error in addition to truncation. The OSE-based bound is rigorous but may be loose in practice, and numerical $\hat \sigma(A) = \sigma \quad \text{with prob. } \frac{|\Tr(\sigma A)|}{2^n \mathcal D(A)}, \qquad \mathcal D(A)=\frac{1}{2^n}\sum_{\sigma\in\mathcal P_n}|\Tr(\sigma A)|,$9 can be much smaller than the worst-case requirement (Shao et al., 25 Oct 2025).

The literature also contains explicit negative results. For noisy circuits, a simple circuit with one layer of 3-qubit unitaries can cause naive path-weight truncation to fail badly: for certain observables and $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$0, the truncated error can scale as $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$1 when $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$2, so $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$3 is not a universal guarantee. At the practical level, threshold sweeps show that smaller $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$4 does not always improve the answer, that convergence can be non-monotonic, and that deeper circuits may in some cases be easier to simulate than shallower ones. This has led to an operational distinction between problems with “apparent convergence,” where PPS acts as a reliable verifier, and problems without apparent convergence, where it is better interpreted as a Monte Carlo-like estimate rather than a standalone predictor (González-García et al., 2024, Gharibyan et al., 14 Jul 2025).

Some additional weaknesses are method-specific. In the original Monte Carlo Pauli-propagation framework, many components with $\mathbb E[\hat c(A)\hat \sigma(A)] = A.$5 lead to exponentially increasing sample complexity; Schrödinger propagation is weak for local observables on many output qubits because the observable cost can blow up; Heisenberg propagation cannot exploit noisy input states to reduce runtime; and adaptive channels are generally not efficiently closed under Pauli propagation unless the Pauli transfer matrix is diagonal (Rall et al., 2019).

Finally, the acronym “PPS” is historically ambiguous. In an earlier and unrelated optics literature, PPS denoted “pseudorandom phase sequence,” used in a single-photon-field construction for effective simulation of Bell and GHZ states. That usage concerns phase-modulated optical simulation and not the Heisenberg-picture Pauli-path framework of modern quantum-circuit and quantum-dynamics simulation (Fu et al., 2010).