paulistrings-rs

Library for simulating quantum circuits in the Pauli basis, built to scale across threads, processes and nodes. Pauli strings are stored as symplectic bitvectors, and terms are assigned to workers by a hash function that is linear over the bitfield. As a result, the library can efficiently parallelize the work and compute the communication patterns for a given circuit. Additionally, for standard gate and noise sets, the amount of communication remains bounded and no global synchronization or reduction is required.

Inspired by PauliStrings.jl.

Quickstart

In many quantum simulations, the object of interest is typically the expectation value $⟨O⟩=\text{tr}(\rho U^{†}OU)$. Here we start from some initial density matrix $\rho$ which is evolved through a circuit $U$ and then measured with an operator $O$. For Pauli propagation methods, instead we work in the Heisenberg picture and work backwards: we start from the operator $O$, which is then evolved backwards through the circuit $U$, and finally measured against a density matrix $\rho$.

Simple circuit

As a simple example, we may look at a four-qubit initial state $|++++⟩$, and measure the total magnetization after propagating through the following simple circuit:

import math
from paulistrings import Circuit, PauliSum

observable = PauliSum.from_strings(
    {"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4
)

circuit = Circuit(4)
circuit.rz(math.pi / 8, 0)
circuit.cnot(0, 1)
circuit.h(2)

evolved = observable.propagate(circuit, direction="heisenberg")
print(len(evolved), evolved.expectation("x+").real)

5 0.7309698831278217

Quenched time-evolution

For a slightly more involved example, we can consider the 2D transverse-field Ising model, and measure the magnetization after a quench from the $|+{⟩}^{\otimes N}$.

$$H=-J\sum _{⟨i,j⟩}{Z}_i{Z}_j-h\sum _i{X}_i$$

This is achieved by propagating the total magnetization $M={\sum }_i{X}_i$ through a Trotterized circuit and measuring against the initial state $|+⟩⟨+{|}^{\otimes N}$. However, in order to keep this computation tractable at longer times, we truncate the intermediate sums of strings, for example by keeping only the largest-magnitude terms up to a fixed count. Finally, we can extrapolate to the limit of no truncation to validate our results.

Two field strengths h (color) and five TopN term caps 10², 10³, 10⁴, 10⁵, 3×10⁵ (opacity, faintest at the smallest cap) on a 6×6 periodic lattice, J = 1. At h = 0.3 the three largest caps already sit on top of each other through t = 0.8: that overlap is the trusted regime, where raising the cap further would not move the curve. At h = 3.1 even the two largest caps are still visibly tightening rather than fully flat by t = 0.8 — the same term budget buys less trustworthy time here, because a bigger single-qubit rotation per Trotter step spreads the operator across more of the Pauli sum faster. The smallest cap, TopN = 10², peels away from the trusted band by t ≈ 0.2 at both field strengths, well before the run ends. Figure script: docs/figures/ising-quench-convergence/quench_convergence.py.

from paulistrings import truncation


def x_magnetization(lx, ly):
    n = lx * ly
    terms = {}
    for site in range(n):
        key = ["I"] * n
        key[site] = "X"
        terms["".join(key)] = 1.0 / n
    return PauliSum.from_strings(terms, num_qubits=n)


def trotter_step(lx, ly, dt, J=1.0, h=1.0):
    n = lx * ly
    circuit = Circuit(n)

    def idx(x, y):
        return (y % ly) * lx + (x % lx)

    for y in range(ly):
        for x in range(lx):
            i = idx(x, y)
            for nx, ny in ((x + 1, y), (x, y + 1)):
                j = idx(nx, ny)
                circuit.pauli_rotation("ZZ", [i, j], 2 * J * dt)
    for site in range(n):
        circuit.rx(2 * h * dt, site)
    return circuit


lx, ly = 3, 3
dt = 0.1
step_circuit = trotter_step(lx, ly, dt)
observable = x_magnetization(lx, ly)
for k in range(1, 6):
    observable = observable.propagate(step_circuit, truncation.coeff(1e-6), direction="heisenberg")
    print(f"t={k * dt:.2f}  <X>={observable.expectation('x+').real:+.6f}  terms={len(observable)}")

t=0.10  <X>=+0.922619  terms=144
t=0.20  <X>=+0.720483  terms=7323
t=0.30  <X>=+0.472221  terms=23067
t=0.40  <X>=+0.271814  terms=57493
t=0.50  <X>=+0.183984  terms=98698

The magnetization falls and the tracked term count climbs by more than two orders of magnitude over five steps at fixed qubit count.

Sections

Source, issues and the full research record: github.com/lkdvos/paulistrings-rs. Dual-licensed MIT OR Apache-2.0.

Installation

pip install paulistrings-rs        # or: uv add paulistrings-rs

The distribution is paulistrings-rs and the import is import paulistrings. Wheels cover Linux x86_64 and aarch64 (manylinux_2_28, glibc 2.28 or newer, which includes every Rusty node), macOS x86_64 and arm64, and Windows x86_64, and one abi3 wheel per platform serves every Python 3.9 or newer. Anywhere else pip builds the source distribution, which needs a Rust toolchain from rustup.

Core propagation needs nothing but numpy. The examples extra adds stim and qiskit, which interop.circuit_from_stim, interop.circuit_from_qiskit and interop.stabilizers_from_stim import on first use, plus qiskit-aer and matplotlib for the example suite:

pip install "paulistrings-rs[examples]"

With the package installed, the Manual is where to go next; it opens with a four-line run and explains each part in turn. First propagation carries that same run further, into term inspection and a validation sweep.

GPUs

The Linux and Windows wheels include the CUDA backend, and need no CUDA toolkit to install. At runtime a device needs the NVIDIA driver’s libcuda and a CUDA 12 libnvrtc on the library search path: module load cuda/12.8.0 on a Flatiron host, pip install nvidia-cuda-nvrtc-cu12 elsewhere (CUDA devices). Without them paulistrings.cuda_available() is False and everything else works.

MPI and multi-node runs

The mpi feature is never in a wheel, since no MPI build is portable across clusters, so it is a source build against the cluster’s MPI. The commands below are for Rusty; on another cluster, substitute its MPI, LLVM and Python modules.

Where to build

• On a Rocky 9 workstation or login node. The build fetches crates from crates.io, and the modules/2.4-20250724 stack is Rocky 9 only.
• With a Rust toolchain from rustup. It is needed at build time only; the nodes run the compiled extension.
• Into a venv on the shared filesystem (/mnt/home or /mnt/ceph), which every node mounts. Local /home and /tmp are invisible to the nodes.

Build

mpi4py must be built against the same MPI as the extension (its pip wheel bundles a different one), so both are built from source:

module load modules/2.4-20250724 openmpi/5.0.6 llvm/19.1.7 python-mpi/3.12.9
export LIBCLANG_PATH=$(llvm-config --libdir)    # rsmpi's bindgen needs libclang

python3 -m venv ~/venvs/paulistrings-mpi
source ~/venvs/paulistrings-mpi/bin/activate
pip install "mpi4py>=4,<5" --no-binary mpi4py
pip install paulistrings-rs --no-binary paulistrings-rs --config-settings=build-args="--features mpi"
python -c "import paulistrings; print(paulistrings.mpi_available())"    # True

A source build replaces the wheel, so --features names every feature the build keeps: --features cuda,mpi for one GPU per rank. The extension records the MPI library directory in its rpath, so import paulistrings works without the modules loaded; load them in the job anyway, since srun and mpi4py need the same MPI.

Check it

Save this as check_mpi.py; it evolves Σ Z under a layer of Hadamards and reads ⟨+|Σ X|+⟩ = n back from the ranks’ shares:

import mpi4py

mpi4py.rc.thread_level = "serialized"  # before mpi4py.MPI is imported
from mpi4py import MPI

import paulistrings
from paulistrings import Circuit, PauliSum, truncation

comm = MPI.COMM_WORLD
n = 64
observable = PauliSum.from_strings({"I" * q + "Z" + "I" * (n - q - 1): 1.0 for q in range(n)}, num_qubits=n)
circuit = Circuit(n)
for q in range(n):
    circuit.h(q)
local = observable.propagate(circuit, truncation.approx_topn(1000), direction="heisenberg", comm=comm, result="local")
value = comm.allreduce(local.expectation("x+"))
if comm.Get_rank() == 0:
    print(f"paulistrings {paulistrings.__version__}: {comm.Get_size()} ranks, <+|sum X|+> = {value.real:g} (expect {n})")

On the build workstation, mpirun -n 2 python check_mpi.py runs it on two ranks, and the job below runs it across nodes.

Run across nodes

One rank per NUMA domain, bound to it, and a power-of-two rank count; the Rocky 9 constraint keeps the job on nodes that can load a Rocky 9 build. Rusty’s CPU nodes have two sockets, each one NUMA domain unless the node reports more (lscpu | grep 'NUMA node(s)' inside a job), so two icelake nodes give four ranks of 32 cores:

#!/bin/bash
#SBATCH --nodes=2
#SBATCH --exclusive
#SBATCH --constraint=icelake&rocky9
#SBATCH --time=01:00:00

module load modules/2.4-20250724 openmpi/5.0.6 python-mpi/3.12.9
source ~/venvs/paulistrings-mpi/bin/activate
srun --ntasks-per-node=2 --cpus-per-task=32 --cpu-bind=ldoms --mpi=pmix python check_mpi.py

For rome nodes use --cpus-per-task=64, for genoa 48. Without --cpu-bind=ldoms every rank’s thread pool spans the whole node and the ranks compete for cores. Every rank must build the same observable and circuit, since each keeps its own share of a sum it holds whole; MPI ranks lists the requirements and the launch options.

GPU and MPI features

Built with --features cuda,mpi, each rank runs its share on one GPU and the ranks exchange over NCCL (MPI ranks). The job loads the CUDA and NCCL libraries, which the extension opens at runtime, and makes every GPU of a node visible to every rank on it, which --gpus-per-node does and --gpus-per-task does not; without that, NCCL falls back from NVLink to copies through host memory:

#!/bin/bash
#SBATCH --partition=gpu
#SBATCH --nodes=1
#SBATCH --gpus-per-node=4
#SBATCH --ntasks-per-node=4
#SBATCH --cpus-per-task=8
#SBATCH --constraint=a100-80gb&rocky9
#SBATCH --time=01:00:00

module load modules/2.4-20250724 openmpi/5.0.6 python-mpi/3.12.9 cuda/12.8.0 nccl/2.23.4-1
source ~/venvs/paulistrings-mpi/bin/activate
srun --cpu-bind=cores --mpi=pmix python script.py    # propagate(..., comm=comm, device="auto")

--constraint=h100-sxm5 selects the H100 nodes instead. Above one node, the NCCL exchange is not yet tested (#40). Every rank holds the whole input in host memory before the scatter (#35), so size --mem-per-cpu for that.

Rust

cargo add paulistrings

The crate is the engine the Python package wraps, with its features off by default: cuda, mpi (which needs an MPI installation and libclang to build) and phase-timing. Its API documentation is a developer guide that states each item’s contracts. Performance numbers in these docs are measured with lto = "fat" and codegen-units = 1 in the release profile.

Developing

A checkout builds with uv, which creates ./.venv from the locked dependencies and builds the extension into it:

git clone https://github.com/lkdvos/paulistrings-rs
cd paulistrings-rs
./scripts/setup.sh                        # uv sync, plus the examples extra
uv run pytest python/paulistrings/tests

The Python version comes from .python-version and the Rust toolchain is pinned in rust-toolchain.toml, so neither needs choosing. uv run rebuilds the extension whenever a Rust source changed, so there is no separate rebuild step; from an activated .venv, uv sync does the same. uv sync removes whatever an extra it is not given had installed, so name every extra each time (uv sync --extra examples --extra bench); uv run only adds. Feature builds of a checkout are in CUDA devices and MPI ranks.

Manual

The Manual is the reading path through the library: what each object is, how to use it, and where it bites. Every computation here has the same shape: write an observable, build a circuit, propagate the observable through the circuit, measure the result against a state, then validate that the truncation did not decide the answer.

import math
from paulistrings import Circuit, PauliSum

observable = PauliSum.from_strings(
    {"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4
)

circuit = Circuit(4)
circuit.rz(math.pi / 8, 0)
circuit.cnot(0, 1)
circuit.h(2)

evolved = observable.propagate(circuit, direction="heisenberg")
print(len(evolved), evolved.expectation("x+").real)

5 0.7309698831278217

from_strings accumulates the observable, the three gates conjugate it in reverse under direction="heisenberg", and expectation reads the result against |++++⟩. Four lines, four chapters: PauliSum.from_strings is Operators, Circuit/rz/cnot/h is Circuits, propagate(..., direction=...) is Engine and propagation, and expectation(...) is Measurements. The four chapters follow that same order.

• Operators — the Pauli sum: how an observable or Hamiltonian is written, stored, inspected and saved.
• Circuits — gates and noise channels, composition, importing from stim, qiskit and task JSON, and why one gate is one truncation point.
• Engine and propagation — the propagate call: direction, truncation, validation, incremental runs, stats and memory, the engine internals, and scaling out over NUMA partitions and MPI ranks.
• Measurements — reading a propagated sum against a product state, a stabilizer state, or another Pauli sum, and which state the label denotes in each direction.

Read the chapters in order once; each assumes the ones before it and does not repeat them. Afterwards, enter by topic from the sidebar: every section is self-contained enough to answer one question, and ends with a pointer to the Library page holding the exact signatures and to an Example that uses the same call at scale.

• Against other tools says whether Pauli propagation is the right method for your problem at all, and names the two places a state-vector or stabilizer simulator is strictly better.
• Library is the signature-and-table reference the Manual cites throughout; it holds no rationale of its own.
• First propagation carries the run above further, into the surviving terms’ identities and a cutoff-sweep validation.

Operators

The primary objects of this library are operators acting on $L$ qubits. For a single qubit, the identity along with the Pauli matrices forms an orthonormal basis. Therefore, we can generate any multi-qubit operator as a linear combination of so-called Pauli strings. These are strings of length $L$ containing the letters $\{I,X,Y,Z\}$, such that all possible strings carry all possible basis elements.

Pauli strings

A single Pauli string $P$ is stored in the so-called symplectic encoding. This means that for each qubit, we store two separate bits that dictate what operator is used. Concretely, for bits $(x,z)\in {\mathbb{F}}_2^2$, we encode the operator in a Hermitian convention as:

$$P(x,z)=i^{xz}{X}^x{Z}^z.$$

In particular, this gives $X=(1,0)$, $Z=(0,1)$, and $Y=iXZ=i(1,1)$, equivalently $XZ=-iY$.

A string on $L$ qubits is stored in $2L$ bits, collected in the symplectic vector $\mathbf{p}=({\mathbf{p}}_x|{\mathbf{p}}_z)\in {\mathbb{F}}_2^{2L}$ with $x$- and $z$-parts ${\mathbf{p}}_x,{\mathbf{p}}_z\in {\mathbb{F}}_2^L$, and

$$P(\mathbf{p})=i^{N_Y(\mathbf{p})}\underset{j=0}{\overset{L-1}{⨂}}{X}^{p_{x,j}}{Z}^{p_{z,j}},N_Y(\mathbf{p})=|{\mathbf{p}}_x\wedge {\mathbf{p}}_z|.$$

Here $\wedge$ is the bitwise AND and $|\mathbf{v}|$ is the Hamming weight of a bitvector $\mathbf{v}$, the number of its nonzero bits, so $N_Y(\mathbf{p})$ counts the $Y$ factors of the string. `All arithmetic on bits and bitvectors is over $\mathbb{F}_2$, with $\oplus$ denoting addition (XOR). The weight of a string is the number of sites on which it acts nontrivially,

$$ w(\boldsymbol{p}) = |\boldsymbol{p}_x \vee \boldsymbol{p}_z|, $$

with $\vee$ the bitwise OR.

The convention used here is that sites are numbered from $0$: site $j$ is character $j$ of a label and qubit index $j$ of every method.

The PauliString type

A single Pauli string is its own type, PauliString, that is, the object the symplectic encoding above describes, before any coefficient is attached to it. It is immutable and carries num_qubits, a weight (the number of non-identity factors, the same quantity truncation.weight caps), and a label (the IXYZ string). str(p) and repr(p) both print the label directly.

from paulistrings import PauliString, p

identity = PauliString.identity(3)
x0 = PauliString.x(0, 3)
y1 = PauliString.y(1, 3)
z2 = PauliString.z(2, 3)

print(identity, identity.weight)
print(x0, y1, z2)
print(p("XYZ") == x0.mul(y1)[1].mul(z2)[1])

III 0
XII IYI IIZ
True

identity, x, y, z and from_label are the five constructors on the class itself; num_qubits is a plain positional argument, not a required keyword, and each single-site constructor takes the qubit index first. p(label) is a shorthand for PauliString.from_label(label), for the common case of writing one down by hand. Two strings compare equal when their labels and num_qubits agree. The constructor, accessor and operation tables are in PauliString.

Single string operations

Four operations on bare Pauli strings are closed-form from the symplectic encoding: whether two strings (anti-)commute, their product, and their (anti-)commutator.

Multiplication

Multiplying two strings is an XOR of their symplectic vectors together with a phase:

$$ P(\boldsymbol{p}), P(\boldsymbol{q}) = i^{k(\boldsymbol{p}, \boldsymbol{q})}, P(\boldsymbol{p} \oplus \boldsymbol{q}), $$

with

$$ k(\boldsymbol{p}, \boldsymbol{q}) = N_Y(\boldsymbol{p}) + N_Y(\boldsymbol{q}) - N_Y(\boldsymbol{p} \oplus \boldsymbol{q}) + 2, |\boldsymbol{p}_z \wedge \boldsymbol{q}_x|, $$

computed modulo 4. mul returns the phase and the product string:

phase, product = p("XZY").mul(p("ZXY"))
print(phase, product)

(1+0j) YYI

It is a good exercise for anyone new to Pauli strings to manually check this output. Regarding package notation for multiplication: operator * between two strings raises a TypeError; * scales by a number, and mul is the product.

Commutation

Two strings either commute or anticommute,

$$ P(\boldsymbol{p}), P(\boldsymbol{q}) = (-1)^{s(\boldsymbol{p}, \boldsymbol{q})}, P(\boldsymbol{q}), P(\boldsymbol{p}), $$

with the bilinear form

$$ s(\boldsymbol{p}, \boldsymbol{q}) = \boldsymbol{p}_x \cdot \boldsymbol{q}_z \oplus \boldsymbol{q}_x \cdot \boldsymbol{p}_z, $$

where the dot products are taken over $\mathbb{F}_2$. They commute when $s = 0$ and anticommute when $s = 1$. In matrix form, $s(\boldsymbol{p}, \boldsymbol{q}) = \boldsymbol{p}^T J \boldsymbol{q}$ with $J = \begin{pmatrix} 0 & \mathbb{I}_L \ \mathbb{I}_L & 0 \end{pmatrix}$, the symplectic form that gives the encoding its name. The commutator and anticommutator are each either $2PQ$ or zero:

$$ [P, Q] = \big(1 - (-1)^{s(\boldsymbol{p}, \boldsymbol{q})}\big), PQ, \qquad {P, Q} = \big(1 + (-1)^{s(\boldsymbol{p}, \boldsymbol{q})}\big), PQ. $$

commutes_with returns True if and only if $s(\boldsymbol{p}, \boldsymbol{q}) = 0$; commutator and anticommutator return the coefficient together with the string $PQ$:

x, z = p("X"), p("Z")
print(x.commutes_with(z), x.commutator(z), x.anticommutator(z))

False (-2j, Y) (0j, Y)

Here $[X, Z] = 2XZ = -2iY$ and ${X, Z} = 0$. The string is reported in both cases, so a zero coefficient is an exact algebraic zero.

Orthonormality

Pauli strings are orthonormal in the normalized Hilbert–Schmidt inner product,

$$ \langle P(\boldsymbol{p}), P(\boldsymbol{q}) \rangle = 2^{-L} \operatorname{Tr}\big(P(\boldsymbol{p})^\dagger P(\boldsymbol{q})\big) = \delta_{\boldsymbol{p}, \boldsymbol{q}}. $$

The expansion $O = \sum_i c_i P_i$ is therefore unique, with coefficients

$$ c_i = \langle P_i, O \rangle = 2^{-L} \operatorname{Tr}(P_i, O), $$

and $O$ is Hermitian if and only if every $c_i$ is real. For $O_1 = \sum_i a_i P(\boldsymbol{p}_i)$ and $O_2 = \sum_j b_j P(\boldsymbol{q}_j)$, only strings present in both operators contribute:

$$ \langle O_1, O_2 \rangle = \sum_{i, j \colon \boldsymbol{p}_i = \boldsymbol{q}_j} \bar{a}_i b_j, \qquad \lVert O \rVert_2^2 = \langle O, O \rangle = \sum_i |c_i|^2. $$

overlap computes $\langle O_1, O_2 \rangle$:

a, b = 1.0 * p("XY"), 1.0 * p("YX")
print(a.overlap(a), a.overlap(b))

(1+0j) 0j

Unitary conjugation preserves $\lVert O \rVert_2$; one can use this to validate results, i.e. measure what truncation removed.

Pauli sums

Storage and construction

The PauliSum type stores the expansion as three columns, one row per term: the $x$-parts $\boldsymbol{p}_x$, the $z$-parts $\boldsymbol{p}_z$, and the coefficients $c_i$. Combining two sums therefore requires the union of their strings with the coefficients of equal strings added, a reduce-by-key operation called the deduplication step, after which terms whose coefficient is exactly zero are dropped. Every constructor and every addition of sums ends with this step.

How does one create a Pauli sum?

from paulistrings import PauliSum, p

observable = 0.25 * p("XIII") + 0.25 * p("IXII") + 0.25 * p("IIXI") + 0.25 * p("IIIX")
print(observable)
print(PauliSum.from_strings({"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}))

0.25*XIII + 0.25*IXII + 0.25*IIXI + 0.25*IIIX
0.25*XIII + 0.25*IXII + 0.25*IIXI + 0.25*IIIX

A string times a number is a one-term sum, and + adds sums. Each + builds a new sum, so assembling $N$ terms one at a time costs $O(N^2)$. PauliSum.from_strings takes all terms at once, as a dict from label to coefficient or as two equal-length sequences of labels and coefficients, in which a repeated label accumulates. num_qubits is inferred from the first label when it is not given.

Arithmetic

+ and - combine two sums with the deduplication step, += and -= do the same in place, and * scales every coefficient by a number:

a = PauliSum.from_strings({"ZZII": -1.0})
b = PauliSum.from_strings({"IZZI": -1.0})
print(a + a)
print(2.0 * (a + b))

-2*ZZII
-2*ZZII + -2*IZZI

Of course, both operands must have the same num_qubits. Note, the product of two sums is not implemented, and * between two PauliSum objects raises TypeError.

A Hamiltonian

The transverse-field Ising chain with open boundaries,

$$ H = -J \sum_{j=0}^{L-2} Z_j Z_{j+1} - h \sum_{j=0}^{L-1} X_j, $$

is a sum of $2L - 1$ strings:

L, J, h = 4, 1.0, 0.5
labels, coefficients = [], []
for j in range(L - 1):
    labels.append("I" * j + "ZZ" + "I" * (L - j - 2))
    coefficients.append(-J)
for j in range(L):
    labels.append("I" * j + "X" + "I" * (L - j - 1))
    coefficients.append(-h)

hamiltonian = PauliSum.from_strings(labels, coefficients)
print(hamiltonian)

-1*ZZII + -1*IZZI + -1*IIZZ + -0.5*XIII + ... (3 more terms)

The same $H$ generates time evolution, where it enters as Pauli rotations rather than as a PauliSum. For $H = \sum_k h_k Q_k$, each Trotter factor is a rotation,

$$ \exp(-i h_k, \delta t, Q_k) = R_{Q_k}(2 h_k, \delta t), \qquad R_Q(\theta) = \exp(-i \theta Q / 2), $$

which is the convention of every rotation in Circuits. For the Ising chain, $h_k = -J$ on the bonds and $h_k = -h$ on the fields, so the angles are $-2J,\delta t$ and $-2h,\delta t$.

Symplectic arrays

x_array() and z_array() return the $x$- and $z$-parts of all $N$ terms as uint64 arrays of shape (N, width), with site $j$ in bit j % 64 of column j // 64. coefficients_array() returns the $c_i$ as a complex128 array of length $N$. All three are copies, and all three list the terms in the same canonical order: partition-bucket index, then lexicographic order of $(\boldsymbol{p}_x, \boldsymbol{p}_z)$. This order is neither insertion order nor label order, and propagate can change it, so rows are paired only within one snapshot. to_arrays() returns all three at once, and PauliSum.from_arrays is the inverse: duplicate rows are summed, arrays narrower than the sum’s width are zero-padded, and a bit set at or beyond site $L$ raises ValueError.

import numpy as np

x, z, c = hamiltonian.to_arrays()
print(x.dtype, x.shape, c.dtype, c.shape)

rebuilt = PauliSum.from_arrays(x, z, c, num_qubits=L)
print(len(rebuilt), np.allclose(rebuilt.coefficients_array(), c))

uint64 (7, 1) complex128 (7,)
7 True

A label is recovered from the two bits of each site, using the encoding table above:

LETTER = np.array([["I", "Z"], ["X", "Y"]])  # LETTER[x bit, z bit]

def bit(words, j):
    return (int(words[j // 64]) >> (j % 64)) & 1

for row in range(len(c)):
    letters = "".join(LETTER[bit(x[row], j), bit(z[row], j)] for j in range(L))
    print(letters, c[row].real)

ZZII -1.0
IZZI -1.0
IIZZ -1.0
XIII -0.5
IXII -0.5
IIXI -0.5
IIIX -0.5

The three $ZZ$ strings come first because their $x$-parts are zero.

Inspecting a sum

len(...) is the number of terms $N$, num_qubits is $L$, and width is the number of words per part. overlap is the inner product $\langle O_1, O_2 \rangle$, so O.overlap(O) is $\lVert O \rVert_2^2$:

print(len(hamiltonian), hamiltonian.num_qubits, hamiltonian.width)
print(hamiltonian.overlap(hamiltonian).real, np.sum(np.abs(c) ** 2))

7 4 1
4.0 4.0

PauliSum(L) is the empty sum on $L$ qubits, the zero operator, and prints as 0.

Saving and loading

A sum is saved as a .npz archive holding the three columns, num_qubits, and a format tag:

from paulistrings import io as psio

psio.save("hamiltonian.npz", hamiltonian)
reloaded = psio.load("hamiltonian.npz")
print(len(reloaded), reloaded.num_qubits, psio.FORMAT)

7 4 paulistrings-npz-v1

save appends .npz when the path lacks it, and load raises ValueError on an archive without the paulistrings-npz-v1 tag.

Circuits

Building a circuit

A Circuit is an ordered list of channels on a fixed register, and every method call on it pushes exactly one channel. Circuit(n) is the empty circuit; the named gates are h, s, sdg, x, y, z, cnot, cz, swap, the rotations rx, ry, rz and pauli_rotation, and the checked matrices unitary_1q and unitary_2q.

Every rotation implements exp(-i·θ·P/2) for its generator P, the same convention as qiskit’s rz. A Trotter step written as exp(-i·s·P), with no /2, therefore needs theta = 2*s; getting this wrong halves every angle and the run still completes, just for a different Hamiltonian. pauli_rotation(pauli, qubits, theta) rotates about a Pauli string of any weight: pauli[k] acts on qubits[k], and identity positions are left out rather than written as I, so pauli_rotation("ZZ", [0, 1], theta) is the two-qubit ZZ rotation.

import math
from paulistrings import Circuit

n, J, h, dt = 4, 1.0, 0.5, 0.05
step = Circuit(n)
for i in range(n - 1):
    step.pauli_rotation("ZZ", [i, i + 1], theta=2 * -J * dt)   # exp(-i * (-J dt) * ZZ)
for i in range(n):
    step.rx(2 * -h * dt, i)                                     # exp(-i * (-h dt) * X)
print(len(step), step.gates[0])

7 {'name': 'pauli_rotation', 'qubits': [0, 1], 'pauli': 'ZZ', 'theta': -0.1}

circuit.gates lists the channels as plain dicts, the task-JSON gate vocabulary, and is the quickest way to check what was actually pushed. unitary_1q(qubit, matrix) and unitary_2q(q0, q1, matrix) take a 2x2 or 4x4 complex matrix, reject a non-unitary one with ValueError, and order the two-qubit tensor factors as |q0 q1⟩, q0 more significant — so a CNOT matrix in the textbook basis is unitary_2q(control, target, CNOT). The same gates exist as factories in the gates module, returning a channel for circuit.append, for building a channel once and pushing it onto several circuits; circuit.append(gates.h(0)) is circuit.h(0).

The full table, with per-gate caveats, is Circuit; C — Deep Trotter circuits is the same shape of step — ZZ rotations then X kicks — on a 127-qubit heavy-hex lattice, 5 to 20 steps deep.

Composition

Circuits compose by channel list. a + b is a new circuit with a’s channels then b’s, a.extend(b) appends in place, both require equal num_qubits, and len, indexing and slicing work on the channel list. circuit.adjoint() is a new circuit with the channel order reversed and every gate replaced by its dagger: rotations negate theta, s becomes sdg and back, matrices are conjugate-transposed, and the self-adjoint Cliffords are unchanged.

front = Circuit(2)
front.h(0)
front.s(0)
back = Circuit(2)
back.rz(0.3, 1)

whole = front + back
print(len(front), len(back), len(whole))
print([g["name"] for g in whole.adjoint().gates], whole.adjoint().gates[0]["theta"])

2 1 3
['rz', 'sdg', 'h'] -0.3

The adjoint is what turns one direction into the other: propagating through c in the Heisenberg picture agrees, to floating-point tolerance, with propagating through c.adjoint() in the forward picture. It exists only for unitary circuits — a noise channel has no adjoint in the time-reversal sense, so adjoint() on a circuit containing one raises ValueError naming the channel; take the adjoint of the unitary part and add the noise to the reversed circuit explicitly. How push order and direction interact is Push order; the composition table is under Circuit.

One gate, one truncation point

Throughout this Manual, a layer is one applied channel — one gate, or one noise channel on one qubit or pair — never a brickwork layer of parallel gates. Truncation runs after every layer in that sense, so len(circuit) is the number of truncation points a propagation will pass, and the per-layer records in PropagationStats have one entry per channel.

Two consequences follow. Fusing two gates into one unitary_2q changes the answer, because it removes a truncation point between them; every circuit in the example suite is built one gate per push for that reason, and it is what makes a per-layer comparison against another engine meaningful at all. And a broadcast call is still one channel per target: depolarize(p, [0, 1]) pushes two channels, not one bundled one.

noisy = Circuit(3)
noisy.h(0)
noisy.depolarize(0.01, [0, 1, 2])
print(len(noisy))

4

Why a truncation point between gates matters, and why noise makes a fixed threshold bite harder rather than softer, is A gate is a truncation point.

Noise channels

Noise is pushed like a gate, one channel per qubit in the list, with the noise module offering the same channels as factories for append.

from paulistrings import noise

circuit = Circuit(4)
circuit.h(0)
circuit.depolarize(0.01, [0])                    # p spread evenly over X, Y, Z
circuit.cnot(0, 1)
circuit.depolarize2(0.01, [(0, 1)])              # p over the 15 non-identity two-qubit Paulis
circuit.dephase(0.01, [1])                       # Z error with probability p
circuit.amplitude_damping(0.01, [2])             # relaxation toward |0> with probability gamma
circuit.pauli_channel(0.002, 0.002, 0.008, [3])  # X, Y, Z with px, py, pz; px + py + pz <= 1
circuit.append(noise.depolarize(0.01, qubit=0))
print(len(circuit))

8

In the Pauli basis a Pauli channel is a coefficient rescale, not a fan-out: depolarize(p) multiplies every non-identity Pauli on its qubit by 1 - 4p/3, and dephase(p) leaves Z alone and scales X and Y by 1 - 2p. pauli_channel(p/3, p/3, p/3, q) is depolarize(p, q) and pauli_channel(0, 0, p, q) is dephase(p, q). amplitude_damping is the one channel here that is not a Pauli channel: it moves Z partly onto the identity, so it changes the key set, and it is not self-adjoint, which is why its Heisenberg and forward results differ in more than a sign.

Every noise channel is non-unitary, so a circuit that carries one has no adjoint(). The rescale is why noise makes a simulation cheaper rather than dearer under a coefficient threshold, the effect B2 measures; the semantics table is Circuit.

Importing circuits

Three importers in the interop module build a Circuit from an outside description. Each maps what it can onto the gate vocabulary above and raises ValueError naming the instruction on anything it cannot; nothing is ever skipped silently. circuit_from_stim and circuit_from_qiskit import stim/qiskit lazily on first call — install the examples or bench extra from Installation to use them.

circuit_from_stim takes a stim circuit object, stim program text, or a path to a .stim file, expands REPEAT blocks, and returns (circuit, observable), where the observable comes from OBSERVABLE_INCLUDE instructions with Pauli targets or is None. Measurements, resets, detectors and sweep or record targets are errors.

import stim
from paulistrings import interop

circuit, observable = interop.circuit_from_stim(
    stim.Circuit("H 0\nCNOT 0 1\nDEPOLARIZE1(0.01) 0 1\nOBSERVABLE_INCLUDE(0) Z0 Z1")
)
print(len(circuit), [g["name"] for g in circuit.gates], observable.coefficients())

4 ['h', 'cnot', 'depolarize', 'depolarize'] [(1+0j)]

circuit_from_qiskit maps the named gates h s sdg x y z cx cz swap rz rx ry directly, rzz/rxx/ryy to pauli_rotation, and anything else with a one- or two-qubit unitary through the checked unitary_1q/unitary_2q path, reordering the two-qubit basis from qiskit’s little-endian convention for you. Measurements, resets, conditioned instructions and gates on more than two qubits are errors; qubit indices come from find_bit, so multi-register circuits import correctly.

from qiskit import QuantumCircuit

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.rzz(0.3, 0, 1)
qc.t(1)
print([g["name"] for g in interop.circuit_from_qiskit(qc).gates])

['h', 'cnot', 'pauli_rotation', 'unitary_1q']

circuit_from_json builds from the frozen task-JSON schema’s "circuit" object, the interchange format shared with the cross-engine comparison, and load_task reads a whole task file — circuit, observable, truncation policy and direction — into one object.

task = interop.load_task({
    "version": 1,
    "n_qubits": 2,
    "circuit": {"gates": [{"name": "h", "qubits": [0]}, {"name": "rz", "qubits": [1], "theta": 0.3}]},
    "observable": {"ZI": 1.0},
    "truncation": {"min_abs_coeff": 1e-8},
    "run": {"direction": "heisenberg"},
})
print(len(task.circuit), len(task.observable), task.direction)

2 1 heisenberg

Unknown keys, unknown gate names and a missing run.direction are errors. Importer signatures are under Circuit; the stim importer’s sibling, stabilizers_from_stim, reads a Clifford circuit’s output state rather than its gates and belongs to Stabilizer states.

Hard edge: support on more than two qubits

The engine applies a channel through a lookup on its support, and that support is capped at two qubits: a channel acting on three or more makes propagate panic, with no fallback path. pauli_rotation is the one exemption and accepts a generator of any weight, so a ZZZ or weight-ten rotation is fine.

From Python the constructors keep you inside the edge: every named gate and noise method is one- or two-qubit, unitary_2q is the widest matrix accepted, and the importers reject a three-qubit instruction with ValueError before it reaches a circuit. The edge is therefore what a three-qubit gate becomes when you need one: a decomposition into one- and two-qubit gates, each its own truncation point, or a Pauli rotation if it is one.

See it in use

• B2 — Noisy circuit verification — a depolarizing channel on every qubit a gate just touched, at 127 qubits, then the same run with each of the other channels swapped in, and the collapse in term count that follows.
• B5 — Hybrid depth reduction — a circuit split into a front and a tail, the tail’s propagated observable written to a task JSON and read back with load_task.
• A — Clifford point — the same recorded gate list handed to the engine and, as schema-v1 task JSON, to PauliPropagation.jl, for a term-for-term parity check.
• C — Deep Trotter circuits — a kicked-Ising Trotter step on 127 qubits, one rotation per push, 5 420 channels at 20 steps.

Engine and propagation

Propagation is the one thing this library does: take a Pauli sum, push it through a circuit one channel at a time, and truncate after every channel so the sum stays tractable. Everything else on this site — observables, circuits, measurements — exists to feed that loop or read its output. This chapter is the loop itself: what propagate computes, which knobs change what it computes and which only change how fast, and how to tell whether the truncated result can be trusted.

The propagate call

from paulistrings import Circuit, PauliSum, truncation

observable = PauliSum.from_strings({"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4)

circuit = Circuit(4)
circuit.rz(0.3, 0)
circuit.cnot(0, 1)

evolved = observable.propagate(circuit, truncation.coeff(1e-10), direction="heisenberg")
print(len(observable), "->", len(evolved), "terms")

4 -> 5 terms

Three arguments carry the physics and one carries the cost. circuit is the channel list to apply, and it must share num_qubits with the sum. direction chooses between U†OU and UOU†, which are different operators and answer different questions — Direction. policy is the truncation applied after every channel, and it is the only argument that trades accuracy for time — Truncation. propagate returns a new sum and leaves the input untouched, so the same observable can be propagated several ways and compared. The GIL is released for the duration of the call.

propagate_with_stats takes the same arguments and returns (evolved, stats), where stats is a PropagationStats with per-layer term counts and timings; enabling it does not change the evolved sum — Stats, memory and logging. The remaining keyword arguments (partitions, comm, device, bucket sizing) change how the sum is stored and where the work runs, never the operator being computed, and are covered on the settings, NUMA partitions, MPI ranks and CUDA devices pages. The full signature is in the Library.

The loop

A layer is one applied channel, not a brickwork layer of parallel gates — One gate, one truncation point is the definition, and it is the sense the word carries throughout this book and in PropagationStats.

For each layer the engine does the same three things. It fans every term out over the channel’s outputs: a Clifford gate moves each term to exactly one new Pauli string, a Pauli rotation splits each term into two (the identity part and the generator part), a depolarizing or dephasing channel rescales coefficients in place, and a dense two-qubit unitary can split a term into up to 16 outputs. It merges terms that landed on the same Pauli string, adding their coefficients. Then it applies the truncation policy to the merged sum, and only what survives is fed to the next layer.

Truncation runs after every layer, not once per propagate call. That single fact has consequences the rest of this chapter returns to: fusing two gates into one channel changes the answer, noise makes a fixed threshold cheaper to satisfy at every depth, and a circuit can be split into consecutive propagate calls with no change in the result.

The order in which the layers are applied depends on direction: forward walks the circuit as written, Heisenberg walks it in reverse and applies each channel’s adjoint. Direction states that rule once and precisely; the Incremental propagation page is where it bites. How the engine makes each layer fast — the bucketed layout and the write-disjoint parallel decomposition — is on Inside the engine.

Cost is terms, not qubits

The cost of a propagation is not the qubit count. It is the number of Pauli strings the operator spreads over, which grows with circuit depth and entangling strength until truncation holds it. A local observable starts as a handful of terms and only spreads as far as its causal cone, so 127 qubits and a shallow circuit can be cheap while 12 qubits and a deep one are not. The qubit count enters only through the number of channels a circuit contains and through the per-term storage width — Estimate the memory a run needs.

Term growth is the cost signal to watch. propagate_with_stats reports it per layer in terms_out, and a short prefix of the circuit is usually enough to see whether the growth is saturating or still exponential. What every truncated result then has to report is how much of the operator the truncation deleted, and whether the answer still moves when the cutoff is tightened — Validating a result.

Always pass direction

The Python binding defaults direction=None to "forward", but the two pictures compute different operators, and most examples on this site use "heisenberg". A result computed in the wrong picture is a plausible number with no error raised, so treat the default as a one-off exploration convenience and pass direction explicitly in anything you keep — Direction has the two numbers side by side.

In this chapter

• Direction — U†OU versus UOU†, which state the expectation label then denotes, and push order.
• Truncation — why the policy is the accuracy knob, which policy to reach for, and why a truncated sum has no variational bound.
• Validating a result — the cutoff sweep, the retained norm, the plateau criterion, and exact references.
• Incremental propagation — splitting a circuit is free, which makes a time series one pass, with the Heisenberg-ordering trap.
• Stats, memory and logging — PropagationStats, the bytes-per-term arithmetic, logging, thread count and engine selection.
• Inside the engine — the bucketed layout, the GF(2) hash, and write-disjoint cosets.
• NUMA partitions — one pinned pool and one share of the sum per NUMA domain.
• MPI ranks — the same split across processes, for capacity.
• CUDA devices — the whole sum on one GPU, and a resident GpuPauliSum for loops of short calls.

First propagation runs this whole chapter once on four qubits; Benchmark D is the cost model measured, with term growth against depth and chain length.

Direction

direction selects which conjugation propagate performs, and the two choices compute different operators. Under "heisenberg" the evolved sum is U†OU; under "forward" it is UOU†. Both are valid, both are cheap or expensive for their own reasons, and neither raises when it is the wrong one for the question being asked — which is why this is the first page of the chapter.

Heisenberg answers “what does this circuit measure?” — you have an observable O and a circuit U, and want ⟨ψ|U†OU|ψ⟩ for a product state |ψ⟩ the circuit acts on. A local observable starts as a handful of terms and only spreads as far as its causal cone, which is why this is the direction almost every example on this site uses.

Forward answers “what does this circuit do to this operator?” — evolving a density matrix or a Hamiltonian in the Schrödinger picture. It starts from a wide operator, so the cost profile is different from the outset.

from paulistrings import Circuit, PauliSum, truncation

observable = PauliSum.from_strings({"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4)
circuit = Circuit(4)
circuit.rz(0.3, 0)
circuit.cnot(0, 1)
policy = truncation.coeff(1e-10)

heisenberg = observable.propagate(circuit, policy, direction="heisenberg")
forward = observable.propagate(circuit, policy, direction="forward")

The default is "forward", and most examples on this site do not use it. The Python binding maps direction=None to "forward", so a call that omits the argument silently computes UOU†. Treat the default as a one-off exploration convenience, not an idiom, and pass direction explicitly in anything you keep.

Which state the label is read against

expectation(state) evaluates ⟨s|A|s⟩ for whatever sum A it is handed, so direction decides what the label s refers to:

• under "heisenberg" the evolved sum is U†OU, so expectation(s) is O’s expectation in U|s⟩ — s is the input state, the one the circuit acts on;
• under "forward" the evolved sum is UOU†, so expectation(s) is O’s expectation in U†|s⟩ — s is a state at the circuit’s output end.

Getting it backwards is silent. The same observable, the same circuit and the same state label give two perfectly plausible numbers:

from paulistrings import Circuit, PauliSum, truncation

observable = PauliSum.from_strings({"ZII": 1.0}, num_qubits=3)

both = Circuit(3)
both.h(0)
both.pauli_rotation("ZZ", [0, 1], 0.6)
both.rx(0.9, 0)
both.rx(0.9, 2)
both.cnot(1, 2)

exact = truncation.coeff(1e-12)
heisenberg = observable.propagate(both, exact, direction="heisenberg")
forward = observable.propagate(both, exact, direction="forward")

print("heisenberg", heisenberg.expectation("x+").real)
print("forward   ", forward.expectation("x+").real)

heisenberg 0.6216099682706644
forward    0.8253356149096783

0.6216… is ⟨Z_0⟩ measured after running this circuit on |+++⟩ — the Heisenberg answer, and the one nearly every question of the form “what does this circuit measure?” wants. 0.8253… is the same observable pushed forward through the circuit and then read against |+++⟩ at the far end, which answers a different question. Both were checked against a dense 2³×2³ construction; neither is a bug.

The rule of thumb: decide which end of the circuit your state lives at, and that decides the direction. A state prepared before the circuit is a Heisenberg run; an operator you want to see after the circuit has acted on it is a forward run. Measurements restates this from the measurement side, for readers who arrive there first.

Push order

Under "heisenberg" the engine walks the channel list in reverse and applies each channel’s adjoint; under "forward" it walks the list as written and applies each channel. So pushing ZZ rotations before X rotations builds U = U_X · U_ZZ, and Heisenberg evolution computes U_ZZ† U_X† O U_X U_ZZ — the X rotations act on O first. Build the circuit in the order the gates act on the state; the direction flag does the rest.

Two places this rule surfaces later in the chapter. PropagationStats lists its layers in application order, so under "heisenberg" the first entry is the circuit’s last gate — circuit_index recovers the position as written (Stats). And an incremental propagation extends the evolved operator at the application end, which under "heisenberg" means a second propagate call prepends its circuit rather than appending it — The Heisenberg-ordering trap.

The full two-row table is in the Library. Every cross-engine comparison on this site is run in the Heisenberg direction, because PauliPropagation.jl defines no forward map for several of the gates in the shared vocabulary — Against other tools has the asymmetry. Showcase B5 is a Heisenberg run whose evolved observable is then handed to a shorter front circuit, which only works because the label denotes the input state.

Truncation

Truncation is what makes Pauli propagation work, and it is the one argument to propagate that trades accuracy for time. It runs after every layer — after every applied channel, not once per propagate call (The loop) — and what it deletes is gone for every later layer. This page is about what that does to the result and which policy to reach for; Validating a result is how to find out whether the result can be quoted.

Why truncate

Each layer fans terms out and merges them, and without a filter the number of distinct Pauli strings grows with depth until it saturates at 4ⁿ or exhausts memory, whichever comes first. For a generic circuit at a generic angle that growth is exponential in depth, and the sum is only tractable because most of the coefficients it produces are tiny. A truncation policy names which terms count as tiny and deletes them at every layer.

The cost of a run is then set by the policy rather than by the circuit alone: the same circuit under coeff(1e-2) and under coeff(1e-5) can differ by orders of magnitude in term count, wall time and memory. The cost model is terms, and the policy is the knob on it.

Choosing a policy

from paulistrings import truncation

policy = truncation.coeff(1e-10)

coeff(eps) is the default choice, and the only knob whose sweep gives a convergence statement that is comparable across engines. It keeps a term when |c| > eps — strictly greater, so a coefficient exactly equal to the threshold is dropped. Reach for it unless you have a specific reason not to.

weight(k) is a blunt instrument. It keeps terms of Pauli weight at most k, which is useful when the causal cone is the real constraint and measurably useless when it is not: on a 2D lattice at high entangling strength it can delete most of the operator for no time saved, and at Clifford angles the operator passes through high weight mid-circuit even when it lands on a low-weight string, so a tight cap can truncate the whole sum to zero terms.

topn(k) bounds memory, but changes what “converged” means. It keeps at most k terms, largest |c| first, so the error is set by the discarded tail rather than by a threshold; a topn run cannot carry a coeff-sweep convergence panel, and it has no equivalent in PauliPropagation.jl for cross-engine comparison. A sum of more than k terms keeps exactly k: when a group of exactly-equal magnitudes straddles the cut, as symmetric problems produce, part of the group survives, taken in storage order, so the kept part depends on bucketing and partitioning. approx_topn(n) is the cheap sibling: it bins terms by octave of |c|², keeps whole octaves from the top down while they fit in n, and fills the rest from the next octave in storage order, so exactly n terms survive but not necessarily the n largest. Exact topn has no partitioned form, so approx_topn is the policy for NUMA partitions and MPI ranks.

collapse_sample(cache, seed) samples instead of filtering. Whenever a gate leaves more than cache strings it draws one string with probability |c|² / Σ|c|² and replaces the sum with it at coefficient 1, so a run is one trajectory of a Monte Carlo estimator and its result is a random variable of the seed. It is unbiased only for weight-diagonal read-outs such as the diagonal echo; average many seeds and report the standard error. PropagationStats.collapses counts the draws.

sampler = truncation.collapse_sample(1_000, 7)

Policies compose with & (keep when both agree, and run every layer pass) and | (keep when either does); | combines per-term filters only, and raises ValueError on a topn, approx_topn or collapse_sample:

policy = truncation.weight(6) & truncation.coeff(1e-10)
budgeted = truncation.coeff(1e-10) & truncation.topn(200_000)

A rotation by a multiple of π/2 is an exact Clifford, one term with no residual branch, so a Clifford circuit needs no threshold. policy=None disables truncation entirely, dropping only exact zeros; it is the exact reference where the problem is small enough to afford it (Untruncated runs).

Whichever policy you pick, the cutoff is a guess until it has been swept. The full policy table is in the Library.

A gate is a truncation point

Because truncation runs after every layer, the granularity of the circuit is part of the computation. Fusing two gates into one channel — replacing rx followed by rz by a single unitary_1q, say — removes a truncation point, so the fused circuit and the two-gate circuit produce different sums under the same policy. Neither is wrong; they are different approximations, and only the untruncated results agree.

import numpy as np
from paulistrings import Circuit, PauliSum, truncation

def rz(theta):
    return np.diag([np.exp(-1j * theta / 2), np.exp(1j * theta / 2)])

def rx(theta):
    c, s = np.cos(theta / 2), np.sin(theta / 2)
    return np.array([[c, -1j * s], [-1j * s, c]])

observable = PauliSum.from_strings({"X": 0.6, "Z": 0.6}, num_qubits=1)

two_gates = Circuit(1)
two_gates.rx(1.0, 0)
two_gates.rz(0.15, 0)

fused = Circuit(1)
fused.unitary_1q(0, rz(0.15) @ rx(1.0))

for policy in (None, truncation.coeff(0.1)):
    a = observable.propagate(two_gates, policy, direction="heisenberg")
    b = observable.propagate(fused, policy, direction="heisenberg")
    print(f"{str(policy):22s} two gates: {len(a)} terms  <Z>={a.expectation('z+').real:+.6f}   fused: {len(b)} terms  <Z>={b.expectation('z+').real:+.6f}")

None                   two gates: 3 terms  <Z>=+0.399630   fused: 3 terms  <Z>=+0.399630
Truncation(Coeff(0.1)) two gates: 3 terms  <Z>=+0.324181   fused: 3 terms  <Z>=+0.399630

The small rz angle produces a Y branch of magnitude 0.6·sin(0.15) ≈ 0.09, which the two-gate circuit truncates at its intermediate point and the fused circuit never sees in isolation. The term counts agree and the answers do not — a 20% difference in ⟨Z⟩ at a cutoff that looks harmless — and no amount of tightening on one circuit says anything about the other.

That is why every circuit in this repository’s example suite is built one gate per push (One gate, one truncation point), and why the cross-engine comparison can be compared per layer at all: both engines truncate once per gate object, so they see the same schedule. It is also why the next page’s cutoff sweep must be run on the circuit as it will be used — a sweep on a fused circuit says nothing about the unfused one.

A noise channel is a truncation point too

A noise channel is a layer like any other, and in the Pauli basis the common ones are coefficient rescales: a single-qubit depolarizing channel multiplies a term’s coefficient by 1 − 4p/3 for each non-identity Pauli on its support, dephasing rescales X and Y, amplitude damping rescales and adds an identity contribution. So adding noise makes a fixed coeff threshold bite harder at every depth — a weight-w string crossing d noise layers is down by (1 − 4p/3)^{wd} before any gate has touched it.

The consequence is the reverse of what a density-matrix simulator sees. There, noise costs extra: a 4ⁿ object instead of 2ⁿ, and each channel more expensive than a gate. Here, noise makes the run cheaper: the filter is exponential in weight, so the threshold acts as an effective weight cap that tightens with depth, and the high-weight tail that dominates a scrambling circuit’s term count is exactly what gets deleted. Showcase B2 measures this on a 127-qubit circuit: at p = 3e-2 the tracked set peaks at 651× fewer terms and finishes 1078× faster than the noiseless run, at the same cutoff.

Two cautions follow. The retained-norm diagnostic on the next page no longer isolates truncation under noise, because the norm falls whether or not anything was deleted. And the granularity rule above applies to noise as well: a broadcast depolarize(p, [0, 1]) is two layers with a truncation between them, not one (Noise channels).

No variational bound

A truncated Pauli sum has no variational bound. Discarded terms carry signs, so a partial sum can sit on either side of the truth, and the error need not fall monotonically as the cutoff tightens.

Benchmark B records a weight-10 observable at θ_h = π/4 whose error against an exact reference is 1.10e-3 at coeff(1e-2), 2.31e-3 at 1e-3 and 2.00e-3 at 1e-4 — the loosest cutoff is the closest. Benchmark D records the same thing on an XXZ chain: 3.3e-6 at 1e-4 then 9.0e-6 at 1e-5 in the free regime, 2.9e-7 at 1e-7 then 4.3e-7 at 1e-8 in the interacting one, while the trend across the whole grid still converges.

Read a convergence sweep as a trend across the grid, never as a point-to-point improvement, and never read agreement between two adjacent cutoffs as an error bound. That is why a single number at a single cutoff is not a result, and why the next page exists. The same absence of a bound also means the sweep’s self-estimated uncertainty can be wrong in either direction; Benchmark C measures a case where a budget-stopped sweep reports an uncertainty 15.7× smaller than the true error.

See it in use

• Benchmark B — Kick-angle sweep and Benchmark D — XXZ chain both record non-monotone truncation error against an exact reference.
• Benchmark C — Deep Trotter circuits shows the self-estimated uncertainty flipping from conservative to falsely confident at depth.
• Showcase B2 — Noisy circuit verification is noise making truncation cheaper, with a dense noisy reference to score it against.
• Benchmark E — Random SU(4) brickwork is the generic worst case for a coeff policy: rise, plateau, then collapse to zero terms.
• First propagation runs a first cutoff sweep on four qubits.

Validating a result

One expectation value at one cutoff is not a result. Truncation has no variational bound, so nothing in the number itself says how far it is from the truth. Three checks turn it into one: sweep the cutoff, read the retained norm, and score against an exact reference wherever that is affordable.

Sweep the cutoff

Run the same propagation at a few decreasing truncation.coeff thresholds and watch the answer, the term count and the retained norm together.

import numpy as np
from paulistrings import Circuit, PauliSum, truncation

n, steps = 12, 6

step = Circuit(n)
for q in range(n - 1):
    step.pauli_rotation("ZZ", [q, q + 1], 1.0)
for q in range(n):
    step.rx(1.1, q)
for q in range(n):
    step.rz(0.6, q)

circuit = Circuit(n)
for _ in range(steps):
    circuit = circuit + step

observable = PauliSum.from_strings({"Z" + "I" * (n - 1): 1.0}, num_qubits=n)
norm_in = float(np.sum(np.abs(observable.coefficients_array()) ** 2))

for eps in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5]:
    evolved = observable.propagate(circuit, truncation.coeff(eps), direction="heisenberg")
    retained = float(np.sum(np.abs(evolved.coefficients_array()) ** 2)) / norm_in
    print(f"{eps:7.0e}  terms={len(evolved):6d}  <Z0>={evolved.expectation('z+').real:+.9f}  retained={retained:.6f}")

  1e-01  terms=     0  <Z0>=+0.000000000  retained=0.000000
  1e-02  terms=   708  <Z0>=+0.301649751  retained=0.856704
  1e-03  terms=  3466  <Z0>=+0.323737940  retained=0.998545
  1e-04  terms=  4392  <Z0>=+0.324826803  retained=0.999998
  1e-05  terms=  4478  <Z0>=+0.324877392  retained=1.000000

Three things to read off, in order. The 1e-1 row truncated the sum to zero terms and reports 0.0; a zero-term sum is not a converged answer however stable it looks, and the run has to be thrown away rather than quoted. The value then moves by 2.2·10⁻² from 1e-2 to 1e-3, by 1.1·10⁻³ to 1e-4, and by 5.1·10⁻⁵ to 1e-5 — a falling sequence of differences, which is what a converging sweep looks like. The retained norm rounds to 1.000000 at 1e-5, but it is not exactly 1 — that row is one term short of the 4 479 an untruncated run returns (see below) — so read it as “the sweep has bottomed out”, not as “this row is exact”.

The differences are the uncertainty estimate, not the cutoff. Quote the tightest row’s value with the last difference as its error bar, and never quote a value from a cutoff you did not sweep past.

The sweep must be run on the circuit exactly as it will be used — same gate granularity, same noise, same direction — because a gate is a truncation point and changing any of them changes the schedule the sweep is measuring.

Read the retained norm

The retained Hilbert–Schmidt norm is N = Σ_P |c_P|², which coefficients_array() gives directly:

def retained_norm(pauli_sum):
    return float(np.sum(np.abs(pauli_sum.coefficients_array()) ** 2))

N is conserved under exact unitary evolution, so the fraction of the operator a truncated run deleted is 1 − N_out / N_in, measured against the same quantity on the input sum. For a single Pauli seed N_in is 1 and the deleted fraction is just 1 − N_out, which is the form the showcases and benchmarks quote.

Two caveats. A noise channel rescales coefficients by design, so under a noisy circuit N falls whether or not anything was truncated and the diagnostic no longer isolates truncation — see A noise channel is a truncation point too. And a large retained norm does not by itself mean a correct expectation value: the deleted terms carry signs and may be exactly the ones the state projects onto, which is why the norm is read alongside the sweep rather than instead of it.

Apply the plateau criterion

Where no exact reference exists, the sweep is the only reference, and the criterion for accepting one is stricter than “two successive values agree”.

A self-converged value may be quoted only when both of the following hold:

• the two most recent successive differences are small against the accuracy you need, and
• either the term count has stopped growing — the sum has saturated, and the plateau is the exact answer — or both of those differences are strictly nonzero.

A flat value on a still-growing sum is rejected: at a small rotation angle, loosening the cutoff by a decade can admit thousands of terms none of which contributes to the state being measured, so the value does not move while the sum keeps growing, and an exactly-zero difference there means “no relevant term has arrived yet”. A sum truncated to zero terms is rejected outright, however flat it looks.

That distinction is not a nicety. The obvious test — tighten the cutoff until two successive values agree to tol — was run against an exact reference at a small kick angle in Benchmark B and declared convergence with an estimated uncertainty of exactly zero while the value was still 5.6·10⁻⁷ from the truth, because at that angle the only terms contributing to ⟨0|O|0⟩ are those rotated all the way to pure Z, and a decade of extra terms contained none. Replacing it with the criterion above is worth a measured 190× in the accuracy of the reported uncertainty, and the benchmark suite imports that criterion as one function object rather than re-implementing it per page.

Truncation has no variational bound, so read a sweep as a trend across the whole grid, never as a point-to-point improvement: a partial sum can sit on either side of the truth and the error need not fall monotonically — see No variational bound. The criterion does not survive every regime either: Benchmark C measures a depth at which a budget-stopped sweep reports a falsely confident uncertainty, and the honest verdict there is “no usable estimate”, not a weaker one.

Score against an exact reference where you can

The sweep estimates its own uncertainty; an exact reference measures the error. Reach for one whenever the problem admits it:

• Drop the policy entirely. A small system or a shallow circuit may fit untruncated, in which case policy=None is the exact answer and the sweep can be scored against it (Untruncated runs). The sweep above bottoms out at 4 478 terms, one short of the 4 479 that observable.propagate(circuit, None, direction="heisenberg") returns.
• Restrict to the causal cone. A local observable only spreads as far as its light cone, so the exact problem is often far smaller than the full register.
• Use a dense simulator. A state-vector reference (qiskit Aer, or a hand-rolled Kronecker construction) is exact to ~26–30 qubits and does not care about depth — see Against other tools.
• Use stim at a Clifford point. Where the circuit is Clifford, a tableau simulator gives the exact ±1 integer at any qubit count — Benchmark A exists to be scored that way.

An exact reference at a reachable corner of the parameter space also calibrates the sweep at the corners you cannot reach.

A time series is a grid of truncated results and carries the same obligation at every point — Incremental propagation. The term counts and memory behind a sweep are on Stats, memory and logging, and the sweep as a cross-engine convergence panel is what every page under Benchmarks ends with.

Incremental propagation

propagate returns a new sum, and that sum is a perfectly good input to another propagate call. Because truncation runs after every layer rather than once per call, chaining calls costs nothing in accuracy, and that identity is what turns a time series from one run per time point into a single pass read out along the way. It also has one trap, which follows directly from push order.

Splitting a circuit is free

Propagating through a and then b applies the same sequence of (apply, merge, truncate) steps as propagating through a + b in one call, so the two agree exactly, in application order:

import numpy as np
from paulistrings import Circuit, PauliSum, truncation

observable = PauliSum.from_strings({"ZIII": 1.0}, num_qubits=4)
policy = truncation.coeff(1e-10)

front = Circuit(4)
front.h(0)
front.cnot(0, 1)

back = Circuit(4)
back.rz(0.3, 1)

split = observable.propagate(front, policy, direction="forward").propagate(
    back, policy, direction="forward"
)
whole = observable.propagate(front + back, policy, direction="forward")

assert len(split) == len(whole)
assert np.allclose(split.coefficients_array(), whole.coefficients_array())

Application order is circuit order under direction="forward", as above, but the reverse of circuit order under direction="heisenberg" (Push order). A Heisenberg split must therefore propagate back before front, not front before back, to reproduce observable.propagate(front + back, ..., direction="heisenberg").

The identity is what Showcase B5 rests on: the tail of a circuit is back-propagated classically, and the shorter front circuit plus the evolved observable are handed to a device, at no cost in accuracy over propagating the whole thing. It is also what makes the next section one pass rather than many.

A time series in one pass

A Trotterized time series costs one propagation, not one per time point: propagate a single step, read the expectation value, propagate the next step from where the last one stopped.

from paulistrings import Circuit, PauliSum, truncation

n, steps = 12, 6

step = Circuit(n)
for q in range(n - 1):
    step.pauli_rotation("ZZ", [q, q + 1], 1.0)
for q in range(n):
    step.rx(1.1, q)
for q in range(n):
    step.rz(0.6, q)

observable = PauliSum.from_strings({"Z" + "I" * (n - 1): 1.0}, num_qubits=n)
policy = truncation.coeff(1e-4)

evolved = observable
series = [(0, len(evolved), evolved.expectation("z+").real)]
for k in range(1, steps + 1):
    evolved = evolved.propagate(step, policy, direction="heisenberg")
    series.append((k, len(evolved), evolved.expectation("z+").real))

for k, terms, value in series:
    print(f"step {k}  terms {terms:5d}  <Z0> {value:+.6f}")

step 0  terms     1  <Z0> +1.000000
step 1  terms     3  <Z0> +0.453596
step 2  terms    14  <Z0> +0.022743
step 3  terms    69  <Z0> +0.241941
step 4  terms   279  <Z0> +0.531583
step 5  terms  1116  <Z0> +0.498573
step 6  terms  4392  <Z0> +0.324827

The last point is identical to propagating the whole six-step circuit in one call — same 4 392 terms, same value, and the same row as the 1e-4 line of the cutoff sweep on the previous page — because truncation runs after every layer either way. The series is therefore free: it is the same single pass, read out at every step boundary.

Each expectation call is one masked pass over the surviving terms and does not consume the sum, so recording several observables per time point costs only those passes. The quenched time-evolution on the landing page is the same pattern at 36 qubits.

The Heisenberg-ordering trap

Under direction="heisenberg" each propagate call conjugates the outside of what is already there, so an incremental call prepends its circuit rather than appending it. Propagating a and then b gives what the single circuit b + a gives, not a + b.

This is invisible above because every Trotter step is the same channel list: U conjugated k times is (U^k)† O U^k whichever way the steps are ordered. It is not invisible for a time-dependent schedule — a ramp, a varying angle, a noise rate that changes with depth. There the loop above computes the steps in the wrong time order, and each time point has to be propagated from the seed with the steps walked in reverse: O(t_k) needs step k applied first and step 1 last. A time-dependent Heisenberg series is therefore not one pass; it is one pass per time point, each from the seed, and the cost difference is the price of asking “what does the circuit up to t_k measure?” for several k.

Under direction="forward" the ordering is the ordinary one — an incremental call appends — so a forward time series over a time-dependent schedule does loop in time order and stays one pass. The two-row table on Direction is the rule this follows from; nothing here is specific to time series.

Keeping the series honest

A time series is a grid of truncated results, so every point carries the same obligation a single number does: sweep the cutoff and check the retained norm, and expect the later points to converge last, since they sit at the largest term count. The sweep on Validating a result is the series above at step 6, swept over five cutoffs, and shows that 1e-4 is one row short of bottoming out.

Term growth across the series is the cost signal — propagate_with_stats reports it per layer rather than per step, and running the loop on a short prefix is the cheapest way to see whether the growth is saturating (Stats). Showcase B1 reads support growth and OTOCs off exactly this kind of series, and Benchmark D measures how the term count at each step scales with the chain length.

Stats, memory and logging

The pages before this one were about what propagate computes. This one is about watching it run: the per-layer counters, how to turn a term count into a memory bill before committing to a run, the log stream, and the handful of arguments and environment variables that change how fast the answer arrives without changing the answer.

Stats

propagate_with_stats returns the same evolved sum propagate would, alongside a PropagationStats with one entry per layer in application order. Enabling it does not change the propagated sum.

from paulistrings import Circuit, PauliSum, truncation

observable = PauliSum.from_strings({"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4)
circuit = Circuit(4)
circuit.rz(0.3, 0)
circuit.cnot(0, 1)
circuit.depolarize(0.01, [0, 1])
policy = truncation.coeff(1e-10)

evolved, stats = observable.propagate_with_stats(circuit, policy, direction="heisenberg")

print(stats.layers, stats.peak_terms, stats.final_terms)
print(stats.terms_in, stats.terms_out)   # one entry per layer, post-truncation
print(stats.gate_name)
print(stats.circuit_index)

4 5 5
[4, 4, 4, 4] [4, 4, 4, 5]
['depolarizing', 'depolarizing', 'cnot', 'PauliRotation']
[3, 2, 1, 0]

A layer is one applied channel (The loop), so the broadcast depolarize(p, [0, 1]) shows up as two entries and len(circuit) == stats.layers. The lists run in application order, which under direction="heisenberg" is the reverse of the circuit as written: the first entry above is the depolarizing channel pushed last. circuit_index recovers the position as written regardless of direction, and application_index is the loop counter.

terms_out[k] is the count after layer k’s truncation, terms_in[k + 1] == terms_out[k], and peak_terms is the maximum of those resident counts. nanos is wall-clock time per layer; for a partitioned or distributed run it is the maximum over partitions, and stats.partition holds the per-partition detail, one partition for an unpartitioned run (NUMA partitions, MPI ranks). The full field list is in the Library.

Estimate the memory a run needs

A term costs 16 × width + 16 bytes: two uint64 key words per width word (the x and z columns) plus a complex128 coefficient. PauliSum.width is the compile-time word tier the qubit count picked — 1 word up to 64 qubits, 2 up to 128, then 4, 8, 16 — so a 127-qubit sum is 48 B/term and a 1024-qubit sum 272 B/term.

stats.peak_terms is the peak resident term count between layers, so peak_terms × (16 × width + 16) is the sum’s own high-water mark:

print(stats.peak_terms, evolved.width, stats.peak_terms * (16 * evolved.width + 16))

5 1 160

Two things that estimate does not include. Each layer also allocates a transient in-layer expansion — a term fans out over the channel’s delta set before the merge collapses it — bounded by the fanout of the widest channel in the circuit (2 for a Pauli rotation, up to 16 for a dense two-qubit unitary) times one coset’s working set per worker, not times the whole sum. And buckets keep their capacity across layers, so a sum that shrinks after its peak does not give the memory back.

Run a short prefix of the circuit first and extrapolate peak_terms from the per-layer growth in terms_out: the growth is what sets the bill, and the truncation policy is the knob on it. A topn or approx_topn budget is the direct way to cap peak_terms when the bill has to be known in advance. The layout this arithmetic comes from is on Inside the engine; Showcase B7 reports this engine’s own measured peak RSS (10.2 GiB single-threaded) against a run this arithmetic would estimate for.

Logging

For running logs instead of a post-hoc summary, the library logs through the log facade under the target paulistrings.propagate — INFO on entry and exit of each propagate call, DEBUG once per layer:

import logging
import paulistrings

logging.basicConfig(level=logging.DEBUG)
logging.getLogger("paulistrings.propagate").setLevel(logging.DEBUG)
paulistrings.reset_log_cache()   # pyo3-log caches each logger's effective level

The INFO lines read propagate: N terms through n channels (direction), scope and propagate: n layers applied, before -> after terms, t s, scope, where scope is 1 partition, P partitions [placement], or partition r/P [placement] for one rank of a comm= run, whose counts are that rank’s own. A layer’s DEBUG line reads layer k/n [name]: before -> after terms, t ms (partition r/P, d remote deltas, m rows in) at every partition count, one line per partition. Placement, pinning and transport setup log under paulistrings.partitioned, and device setup (kernel compiles and the kernel cache, device choice, the NCCL and peer wires) under paulistrings.gpu, so each filters on its own.

Call reset_log_cache() again after changing the level mid-process; the cache is otherwise stale for the rest of the process. Leave logging off when timing — an enabled DEBUG filter adds a clock read per layer.

Every PAULISTRINGS_* environment variable is read once per process, the first time it is needed, and its effective value is logged at DEBUG under the target of the subsystem it tunes; an unparsable value logs a warning and the default applies. Set them before the first call, and alike on every rank of a comm= run.

Threads

In a single unpartitioned process, set the thread count via RAYON_NUM_THREADS before the interpreter starts — Rayon builds its global pool at the first propagate call and never resizes it, so setting the variable from inside an already-running script does not reliably reach it:

RAYON_NUM_THREADS=32 python my_script.py

RAYON_NUM_THREADS does not reach a partitioned run; there the thread count comes from the placement (NUMA partitions).

More threads are not always faster, and the layer class decides. On the reference two-socket host, dense-PTM-heavy circuits (general two-qubit unitaries) saturate the machine’s write bandwidth at about 16 threads and extra threads only add contention, while sparse-rotation and Clifford circuits stay latency-bound and take the full thread count profitably — 11–13× at 32 threads. Hyperthreads add no bandwidth and help only the latency-bound phases. The measured roofline behind those numbers is on Engine performance.

Every benchmark on this site is run single-threaded, RAYON_NUM_THREADS=1 exported before the interpreter starts, so that term counts and wall times are comparable across pages (Comparability rules).

Bucket sizing

target_bucket_len and min_buckets size the bucketed engine’s storage per layer; PauliSum.num_buckets reads back what was realized, which only ever grows. They are tuning knobs for the engine’s own benchmarks rather than for users, and the defaults are what every number on this site was measured with. target_bucket_len must be at least 1 and min_buckets at least 16; anything smaller is a ValueError.

Determinism

The correctness bar is agreement to floating-point tolerance, not bit-identical output. Equal-key summation order is unspecified: two runs of the same propagation can differ in the last bits of a coefficient, and a partitioned run can differ from an unpartitioned one, all within tolerance and all equally correct. Compare evolved sums with numpy.allclose on coefficients_array(), as Splitting a circuit is free does, never with equality.

Term counts, on the other hand, are exact integers and are the load-independent quantity every benchmark on this site quotes and gates on.

Untruncated runs

policy=None applies no truncation at all; only coefficients that are exactly zero after a merge are dropped. The result is the exact evolved operator, and it is the reference the cutoff sweep is scored against wherever the problem is small enough to afford it.

exact = observable.propagate(circuit, None, direction="heisenberg")
print(len(exact), "terms untruncated")

5 terms untruncated

Afford it deliberately: term growth without a filter is exponential in depth for a generic circuit. Estimate peak_terms from a prefix and the bytes-per-term arithmetic above before running a deep circuit untruncated, and expect the untruncated run to be the most expensive point of any sweep by a wide margin.

Inside the engine

paulistrings evolves a weighted sum of Pauli strings through a circuit layer by layer, forward or in the Heisenberg picture, with truncation keeping the sum tractable. This page explains the engine design that makes that loop fast: the sum never maintains a global sorted order, every layer’s data movement is known before any term is touched, and the parallel decomposition needs no locks, no atomics, and no synchronization inside a layer. The companion Engine performance page puts measured numbers to each claim.

The propagation loop

A PauliSum stores each term as a symplectic key: qubit i’s Pauli is a bit pair, packed into fixed-width x/z bit arrays, one machine word per 64 qubits, with a complex coefficient. Under this encoding Pauli multiplication is bitwise XOR of keys plus a phase, and the phase is folded into the coefficient at the boundary, so a key is a plain element of the vector space GF(2)^{2n}.

Each circuit layer applies one channel to every term. A channel with support on k qubits maps an input key v to outputs v ⊕ d, where d ranges over a small delta set D determined by the channel alone: a Pauli rotation of any generator weight has |D| = 2 (identity and the generator), a two-qubit Clifford at most 4, a dense two-qubit unitary at most 16. Applying a layer means fanning every term out over D, then recombining terms that landed on the same key.

The recombination is the interesting part. Done naively it is a global sort or a concurrent hash map per layer, at millions to hundreds of millions of terms.

No global sort

The sum is partitioned into buckets, and its canonical order is bucket index ascending, then lexicographic key within a bucket, with no duplicate keys anywhere. Deduplication therefore only ever needs a canonical order within one bucket, and every bucket is small and cache-resident. Nothing in the propagation loop ever sorts the whole sum.

Storage is structure-of-arrays per bucket: parallel columns for the x key bits, the z key bits, and the complex coefficients. Coefficient-only scans (truncation, expectation values) and key-only scans (weight, commutation) each stream exactly the bytes they use, and every column maps directly to a GPU device buffer. Buckets own their columns, so each retains its capacity across layers: the steady state of a propagation loop allocates nothing.

GF(2)-linear bucket hash

The bucket function is h(v) = H·v for a fixed dense random matrix H over GF(2), giving B = 2^b buckets. Linearity is the property everything else follows from:

h(v ⊕ d) = h(v) ⊕ h(d)

A term in bucket i fans out only into buckets i ⊕ h(d) for d ∈ D — so the layer’s entire data-movement pattern is a function of the channel and the hash, known before touching a single term. Equal keys always share a bucket (h is a function), so duplicates can never straddle buckets. Because H is dense and random, bucket loads stay uniform whatever structure truncation imposes on the keys, and the map from delta to bucket offset is almost always one-to-one.

Bucket count follows the sum’s size, targeting about a thousand terms per bucket so a bucket and its scratch fit in L2. Refining the partition when the sum grows is a single parity pass per term, not a re-sort: adding one row to H splits each bucket in two with the within-bucket order inherited.

Write-disjoint cosets

Take the span of h(D) in bucket-index space. Its cosets partition the buckets, and every output bucket in a coset reads only input buckets in that same coset: a coset is a closed task. That is the engine’s whole parallel decomposition — one coset per work-stealing task, and by construction no two tasks touch the same bucket, so a layer runs with no atomics, no locks, and no cross-thread reconciliation of any kind. Load balance comes from the random hash plus work-stealing.

Output agrees to floating-point tolerance across runs — equal-key summation order is deliberately unspecified.

In-place layer execution

Each coset task, independently:

1. Swap its bucket columns into worker-persistent scratch, leaving emptied, capacity-retaining columns behind as the write destinations. The layer is in-place: peak memory is the sum itself plus one coset’s working set per worker, never a second full-size copy.
2. Size every output run exactly from the swapped-out bucket lengths, before any work.
3. Gather input-major: each term is loaded exactly once and its whole fanout scattered to output runs via an O(1) index identity, with no read amplification. The inner loop per (term, delta) is a few bit extractions, one table lookup, one XOR mask, one complex multiply.
4. Per run, sort by key and merge straight into the bucket’s live slot, restoring the strictly-ascending no-duplicates invariant.

Two structural tricks keep the moved bytes down. The identity delta’s stream keeps its source keys, so it inherits the source bucket’s sorted order and is never sorted; and when the identity amplitude is dense (every rotation and general unitary), the gather materializes only the 16-byte coefficients and the merge borrows the key columns from the source bucket in place — identity keys are neither written nor re-read. The rest-stream sort itself exploits that a gather run is a concatenation of already-sorted per-delta streams; a radix variant takes over for dense two-qubit PTMs, where the sort dominates the layer. The merge is fused with the segmented reduction: one two-pointer walk sums equal-key coefficients and applies the truncation filter to each fully summed coefficient as it passes.

Channels that never change keys at all (depolarizing, dephasing, Pauli gates) bypass the whole pipeline as a parallel coefficient rescale.

Channel preparation

Applying a channel through its general single-term interface per term would pay an indirect call and re-derived tables millions of times. Instead the engine prepares each channel once per layer into a local Pauli-transfer-matrix form: per delta, a bucket offset, full-width XOR masks, and an amplitude table with all phases pre-folded. The preparation is derived automatically for any channel with support on at most two qubits by probing the channel’s own single-term action on the local Pauli basis — a custom channel defines only that action and gets the bucketed engine for free. Pauli rotations above two qubits prepare specially, since their delta set is just {0, generator} at any weight.

Truncation

A truncation policy has two hooks split by cost. The per-term test runs on every merged output — potentially billions of times — sees the fully summed coefficient inside the merge, and is compiled directly into it, so a threshold test inlines to a compare with no call. The per-layer step runs once per layer and may be non-local; TopN and its histogram-based approximation live there. Policies compose with and/or combinators; or combines per-term filters only. Magnitude comparisons use |c|² rather than |c| (the same ordering, no hypot on the hot path), and TopN/ApproxTopN keep exactly n, filling from the boundary in bucket order.

Further reading

The rustdoc under /api/ documents the public types; the module docs carry the engine’s contracts next to the code that keeps them.

The same loop runs with the sum split across NUMA domains or MPI ranks — NUMA partitions, MPI ranks and GPU.

Sources: module docs engine/coset.rs, pauli_sum/hash.rs, engine/bucketed.rs, engine/merge.rs, pauli_sum/mod.rs.

NUMA partitions

A partitioned run splits the sum across NUMA domains instead of sharing one thread pool across the whole box. Each domain gets its own pinned pool and its own share of the terms, selected by designated rows of the GF(2) hash, and a layer exchanges only the rows that cross a domain boundary. The engine it runs on is Inside the engine; this page is when to reach for it and how. To split across processes instead, for capacity rather than bandwidth, see MPI ranks.

When partitioning pays

Locality decides everything, in the sense Engine performance measures in full: a dense two-qubit layer with no row crossing runs 4–18% faster at P = 2 than P = 1, while any layer that exports rows runs 2–8.4× slower.

Random partition rows put roughly half of a dense two-qubit gate’s deltas across a boundary, so the default draw lands a mixed circuit in the slow case. In-process partitioning is worth it when the exchange is rare, and there are two ways to make it rare. When the circuit has a known geometry, partition_row_blocks= replaces the random draw with an explicit cut — one contiguous block of qubits per partition — so only gates that straddle the cut exchange; the kwarg is described under MPI ranks and works identically under partitions=. On a heavy-hex kicked-Ising step, cut rows made P = 2 15–21% faster per step than the single-process engine (research/FINDINGS.md §Cut partition rows). Choosing rows automatically for a circuit with no obvious geometry is open research (research/FINDINGS.md §Partition rows without a known lattice).

Python

The snippets below are skipped by the doc checker because they need a multi-node NUMA host to do anything.

evolved = observable.propagate(
    circuit,
    truncation.approx_topn(1_000_000),
    direction="heisenberg",
    partitions="auto",            # one partition per NUMA node, rounded down to a power of two
)

partitions takes "auto", an integer (a power of two, at most 64), or an explicit list of CPU lists — one per partition, in the "0-7,16-23" spelling — when the automatic split is not the one you want:

partitions=["0-7,16-23", "8-15,24-31"]

partitions=None or partitions=1 is the ordinary unpartitioned run. On a single-node box "auto" also falls back to the unpartitioned path, and an integer larger than the node count this process may run on is a ValueError that points at explicit CPU lists.

paulistrings.numa_nodes() reports what the machine offers, as one CPU list per node intersected with this process’s affinity mask, so a script can decide for itself before choosing. pin_memory=False pins the threads but not their allocations; the default binds both, which is the point of the placement.

PropagationStats.partition holds, per layer, the terms each partition held, the rows and bytes it sent, and the imbalance across partitions. An unpartitioned run reports one partition. local[k] says whether layer k exchanged at all, which is the figure the cost/benefit split above turns on.

_, stats = observable.propagate_with_stats(circuit, policy, partitions="auto")
print(stats.partition.rows_exported)   # rows sent across a boundary, per layer
print(stats.partition.local)           # True where a layer exchanged nothing

The other fields are under PartitionStats.

What changes under partitioning

RAYON_NUM_THREADS is ignored. The partitioned engine builds one pool per partition from the placement, not from Rayon’s global pool. An explicit CPU-list placement takes each partition’s thread count from its list length; "auto" and an integer count take it from the CPUs of the partition’s NUMA node inside this process’s affinity mask. The variable still governs an unpartitioned run (Stats, memory and logging).

Exact topn is unavailable. Choosing the n-th largest magnitude across partitions is a distributed selection, and the engine has no collective form for it, so a partitioned run raises NotImplementedError on truncation.topn. Use truncation.approx_topn(n), which is partition-exact: its histogram is all-reduced, so the retained set is exactly the set the single-partition run would have kept. coeff and weight are per-term filters and need nothing.

Limits

• Pinning is Linux-only. Elsewhere the topology module reports one node and pins nothing, so the run is correct but unplaced.
• P is a power of two, at most 64. A partition is named by log2(P) GF(2) hash rows, and the engine caps that at 6 (P_MAX_BITS).
• Peak memory is above a single-pool run’s. A partitioned run holds P shares of the sum plus one layer’s exchange buffers, which are pooled and reused rather than freed between layers.
• Results agree to floating-point tolerance, not bit for bit — as they do across bucket counts.

See it in use

Engine performance has the P = 2 versus P = 1 tables per host and thread count, the per-layer times of an exchange-free layer, the bandwidth counters behind the gain, and the benchmarking protocol for reproducing them.

MPI ranks

A distributed run is a partitioned run with one partition per MPI process: each rank holds a disjoint share of the terms, selected by designated rows of the GF(2) hash, and a layer exchanges only the rows that cross a rank boundary. Everything on the NUMA page still applies — the same layer loop, the same truncation rules, the same things that change. What is new is the transport (point-to-point MPI instead of in-process channels) and the launch.

Why: capacity

The reason to reach for it is capacity. A sum that does not fit one node’s memory fits P of them, and the per-rank overhead is bounded and does not grow with the rank count. Measured at one rank per NUMA domain on Icelake ccq nodes, a layer with no row crossing costs about the same at 2, 4 and 8 ranks, and a Pauli rotation whose generator crosses costs 3.5× that intra-node and 4.4–4.7× inter-node — see Engine performance for the full table. A remote layer is transfer-bound, so the lever is fewer bytes (locality rows, a lower truncation), not more threads.

It is an off-by-default build option the released wheel omits, so the default wheel cannot do it. paulistrings.mpi_available() says whether this build can, and comm= in a build without it raises RuntimeError.

Building from source

The mpi feature needs an MPI installation (for mpicc, which rsmpi’s build script probes) and a libclang for its bindgen; Installation has the pip build against a cluster’s MPI, a check script and Slurm job scripts.

mpi4py must be built against the same MPI the extension links; the boundary checks the width of MPI_Comm and refuses a mismatch rather than corrupting a handle. In a checkout, scripts/sync-mpi-venv.sh builds a venv ./.venv-mpi with uv: the extension with --features mpi, and mpi4py compiled from source against the loaded MPI.

module load modules/2.4-20250724 openmpi/5.0.6 llvm/19.1.7 python-mpi/3.12.9 uv
export LIBCLANG_PATH=$(llvm-config --libdir)

scripts/sync-mpi-venv.sh          # re-run after any change; --cuda adds the cuda feature
source .venv-mpi/bin/activate

The install is non-editable, so the venv keeps the build it was synced from. The built extension carries an rpath to that MPI’s library directory, so import paulistrings works from a shell with no modules loaded. scripts/mpi-test.sh --ranks 2,4 --python syncs that venv and runs its installed test_mpi.py under mpirun at each rank count.

Python

The snippets on this page are skipped by the doc checker because they need an mpi build and an MPI launcher.

import mpi4py
mpi4py.rc.thread_level = "serialized"      # before mpi4py.MPI is imported
from mpi4py import MPI

import paulistrings
from paulistrings import truncation

comm = MPI.COMM_WORLD
observable = build_observable(seed=20260908)   # the SAME terms on every rank
circuit = build_circuit()

evolved = observable.propagate(
    circuit,
    truncation.approx_topn(1_000_000),
    direction="heisenberg",
    comm=comm,
)
if comm.Get_rank() == 0:
    print(len(evolved), evolved.expectation("z+"))

comm and partitions are alternatives — passing both is a ValueError. A distributed run places one partition per process, so the placement is the launcher’s job (Launching).

Partition row selection

By default the rows are a GF(2)-random draw from the sum’s own hash seed, which spreads a gate’s deltas over every rank whatever the qubits’ geometry. partition_row_seed= picks the draw; partition_row_blocks= replaces it with an explicit locality cut, one contiguous list of qubit indices per rank, so a gate whose qubits share a block never exchanges at all:

half = observable.num_qubits // 2
blocks = [list(range(half)), list(range(half, observable.num_qubits))]
evolved = observable.propagate(circuit, policy, comm=comm, partition_row_blocks=blocks)

Both kwargs work the same way under partitions=, and the two are alternatives to each other. The block count must equal the rank count, the blocks must be disjoint, and they must be identical on every rank — the split is a local filter each rank computes for itself. A term’s rank is then the XOR of the blocks it has odd Z-weight in, so a term supported inside one block belongs to that block’s rank.

partition_row_exclude= keeps the random draw but forbids chosen coordinates, as a dict of "x" and "z" qubit lists:

everything = list(range(observable.num_qubits))
evolved = observable.propagate(circuit, policy, comm=comm, result="local", partition_row_exclude={"z": everything})

A gate that never changes the excluded bits never moves a string between ranks: CZ and rz leave x-bits alone, so rows drawn on x-bits only keep a CZ/rz circuit’s exchange to its rx gates. Excluding coordinates also keeps strings that differ only there on one rank, which a per-rank echo read-out needs. It cannot be combined with partition_row_blocks=.

result=“gather” versus result=“local”

The shares are disjoint and cover the result, so reductions over comm are exact:

local = observable.propagate(circuit, policy, comm=comm, result="local")
terms = comm.allreduce(len(local))
value = comm.allreduce(local.expectation("z+"))

Prefer "local" when the answer is a scalar: gathering materializes the whole sum on rank 0, and that rank transiently holds one extra rank’s copy on top of it.

Launching

One rank per NUMA domain, bound to it — that is the placement the in-process engine would have made with explicit CPU lists, done by the launcher instead:

mpirun -n 4 --map-by ppr:1:numa --bind-to numa python script.py
srun --ntasks-per-node=2 -c 32 --cpu-bind=ldoms --mpi=pmix python script.py

--ntasks-per-node is the node’s NUMA-domain count and -c its cores per domain. Each rank’s Rayon pool sizes itself from the CPUs the launcher left in its affinity mask, so an unbound launch gives every rank a pool over the whole node and they fight. RAYON_NUM_THREADS is ignored, as in any partitioned run.

On Rusty, scripts/slurm/mpi-ranks.sbatch does the arithmetic: it reads the node’s domain count, rounds nodes × domains down to a power of two, and runs the differential net and then the probe at that rank count. See scripts/slurm/README.md.

One GPU per rank

Built with both mpi and cuda, the Rust driver gpu::MpiGpuSum holds each rank’s share on one CUDA device, with the same contract as MpiSum: replicated input, collective calls, rank 0 gathers. gpu::local_device_for_comm picks the device collectively over the ranks sharing a node: each rank gets a device on a NUMA node its CPUs are on where it can, ranks on one node get distinct devices while there are enough, and without readable NUMA facts it falls back to the rank’s node-local index modulo the visible devices:

#![allow(unused)]
fn main() {
use paulistrings::gpu::{local_device_for_comm, MpiGpuSum};
use paulistrings::mpi::{rsmpi, MpiTransport};

let (universe, _) = rsmpi::initialize_with_threading(rsmpi::Threading::Serialized).unwrap();
let world = universe.world();
let device = local_device_for_comm(&world)?;
let transport = MpiTransport::from_communicator(&world);
let mut split = MpiGpuSum::<2>::scatter_to_device(&observable, transport, device)?;
split.propagate(&circuit, ApproxTopN(10_000_000), Direction::Heisenberg)?;
if let Some(evolved) = split.gather()? {
    println!("{} terms", evolved.len());
}
}

scatter_to_device_with(&observable, transport, GpuScatterOptions { rows, ..GpuScatterOptions::new(&[device])? }) splits by other rows or layer options; gpu::propagate_mpi_gpu is the one-shot form. A device failure on any rank fails the call on every rank, with GpuError::Poisoned naming the failing rank on its peers, so the group never falls out of step. Launch with every node GPU visible to every task (srun --gpus-per-node=<k> --mpi=pmix, as scripts/slurm/mpi-gpu.sbatch does), where the locality pick matters.

A remote layer’s columns move device to device over NCCL (CUDA devices), which needs the NCCL library at run time (Installation). The communicator starts once, collectively, at the scatter that builds the split; a rank that cannot load NCCL, or two ranks sharing a device, fail the scatter on every rank. --gpus-per-node rather than --gpus-per-task=1 is what lets NCCL find every rank’s device and use NVLink between them.

From Python, propagate(..., comm=comm, device="auto") is the same run, in an extension built with --features cuda,mpi: device= takes this rank’s ordinal or "auto", which is local_device_for_comm, and result=, partition_row_seed= and partition_row_blocks= keep their host meanings (CUDA devices). scripts/mpi-test.sh --ranks 2,4 --python --cuda builds that extension and runs test_mpi.py’s device cases, which skip on every rank unless every rank sees a device and, above one rank, the group can start NCCL.

Requirements

Thread level at least MPI_THREAD_SERIALIZED. The layer loop runs inside a pinned Rayon pool, so MPI is called from a pool worker rather than the process’s main thread — and not necessarily the same worker on every layer. Only ever one at a time, which is exactly SERIALIZED; FUNNELED would be a false claim. Set mpi4py.rc.thread_level before from mpi4py import MPI (mpi4py’s default, "multiple", is also fine). A weaker level is a RuntimeError at the boundary.

A power-of-two rank count. A partition is named by log2(P) GF(2) hash rows, so mpirun -n 3 is a ValueError.

A replicated input. Every rank must call propagate with the same observable, circuit, policy, direction and options. The scatter is a local filter of a sum every rank already holds, not a distribution of rank 0’s copy — build the terms from the same seed, or load the same file, on every rank. Nothing checks the terms; the run’s shape is checked, and a group that disagrees about the circuit aborts with a fingerprint mismatch rather than hanging.

propagate is collective: every rank of comm must reach it, in the same order. A rank that skips one hangs the rest.

Stats

propagate_with_stats(..., comm=comm) fills PropagationStats.partition with a PartitionStats whose rank and size name the group (both are None for a run that is not distributed). Its per-layer lists hold this rank’s entry only — gathering the group’s counters would mean a collective per layer for a diagnostic — so terms_in[k] is a one-element list, rows_exported[k] is what this rank sent, and imbalance[k] is always 1.0. Reduce over comm for the group’s picture:

_, stats = observable.propagate_with_stats(circuit, policy, comm=comm)
exported = comm.allreduce(sum(stats.partition.rows_exported))

The field list is under PartitionStats.

Limits

• Exact topn is unavailable, as in any partitioned run: choosing the n-th largest magnitude across ranks is a distributed selection. truncation.approx_topn(n) all-reduces its octave histogram and retains exactly the set a single-process run would have.
• A replicated input caps the sum at what one rank can build. Reaching the capacity the ranks together have needs a driver that ingests already distributed; the scatter itself is a local filter and costs nothing extra.
• Homogeneous groups only. The wire carries raw host bytes, so a mixed-architecture job, or one mixing compile-time width tiers across ranks, is silently wrong.
• One partition per rank. The hybrid — several NUMA domains inside one rank — is not implemented. Intra-node, that is what would take the 2-rank remote layer below the 3.5× in Engine performance.
• The library never initializes MPI. The application owns MPI_Init and MPI_Finalize.

CUDA devices

A device run holds the sum on CUDA GPUs and applies every layer there, with the same layer loop, truncation rules and trace as a partitioned run: one partition on one device, one partition per listed device, or one device per MPI rank. Each output bucket is one thread block that builds, sorts and merges its incoming rows in shared memory, so a layer never materializes a gather row or runs a global sort. The engine it runs is Inside the engine; this page is when to reach for it and how.

When a device pays

Dense two-qubit layers at a million terms and up are where the device wins. On an RTX A6000 a saturated su4 layer runs 11× faster on wall than the 16-thread host engine at 1.41e7 terms, and 15× at 5.65e7 (research/HARDWARE.md § ccqlin038 — GPU). Sparse layers (cnot, a Pauli rotation) gain 2–5×, because each output bucket receives few rows and a block’s fixed sort and synchronization cost weighs more. Every layer carries a fixed cost of a few launches and synchronizations that the host engine does not pay, so the advantage shrinks below about 1e4 terms.

The upload and download are paid once per propagate(device=...) call. A loop of many short calls — a Trotter time series, a cutoff sweep over the same prefix — should hold the sum on the device with PauliSum.to_device instead.

The Linux and Windows wheels include the backend; the macOS wheels, with no CUDA driver to load, omit it. paulistrings.cuda_available() says whether this build can run on a device and one is visible, and device= in a build without the feature raises RuntimeError.

Runtime libraries

The backend needs no CUDA toolkit at build time: the kernels are CUDA C++ compiled by NVRTC when a width is first used, and the driver and NVRTC libraries are loaded at runtime. At runtime the process needs the NVIDIA driver’s libcuda and a CUDA 12 libnvrtc on the library search path. On a Flatiron host module load cuda/12.8.0 provides the latter; elsewhere pip install nvidia-cuda-nvrtc-cu12 does, with its lib directory added to LD_LIBRARY_PATH:

pip install nvidia-cuda-nvrtc-cu12
export LD_LIBRARY_PATH="$(python -c 'import nvidia.cuda_nvrtc as m; print(list(m.__path__)[0])')/lib:$LD_LIBRARY_PATH"

On Windows the libraries are found through PATH instead.

A missing library makes cuda_available() return False rather than fail at import. The first propagation at each width compiles its kernels, which takes several seconds; later calls in the process reuse them, and an on-disk cache ($PAULISTRINGS_KERNEL_CACHE, default ~/.cache/paulistrings/kernels, off to disable) makes the first call in a later process reuse them too.

Building from source

A checkout builds the feature into its uv venv:

export MATURIN_PEP517_ARGS="--features cuda"   # uv rebuilds whenever this changes, so keep it set for this venv
uv sync

Python

The snippets on this page are skipped by the doc checker because they need a CUDA device.

import paulistrings
from paulistrings import truncation

if paulistrings.cuda_available():
    evolved = observable.propagate(
        circuit,
        truncation.approx_topn(10_000_000),
        direction="heisenberg",
        device=0,
    )

device takes a device ordinal, a list of ordinals, or "auto". device= is an alternative to partitions=, and passing both, or result="local" without comm=, is a ValueError. target_bucket_len and min_buckets still set the host-side bucket schedule the device refines on top of; engine is ignored.

propagate_with_stats(..., device=0) fills PropagationStats.partition as a one-partition run, and PartitionStats.devices names the device:

_, stats = observable.propagate_with_stats(circuit, policy, device=0)
print(stats.partition.devices)   # [0]
print(stats.nanos)               # wall time per layer, on the device

Resident sums

PauliSum.to_device(device) uploads a sum and returns a GpuPauliSum, which propagate steps in place and to_host copies back:

resident = observable.to_device(0)
values = []
for _ in range(steps):
    resident.propagate(trotter_step, policy, direction="heisenberg")
    values.append(resident.to_host().expectation("z+"))

GpuPauliSum.propagate and propagate_with_stats take policy, direction, target_bucket_len and min_buckets with the meanings above; len(), num_qubits, device and num_buckets read the resident sum without a download. A device error mid-run leaves the resident sum holding the last completed layer’s output, and a later call resumes from it.

Several devices

A list places one partition on each entry, split by GF(2) partition rows and exchanged between layers as in a partitioned run:

evolved, stats = observable.propagate_with_stats(
    circuit, truncation.approx_topn(10_000_000), direction="heisenberg", device=[0, 1, 2, 3]
)
print(stats.partition.devices)   # [0, 1, 2, 3]

The list length must be a power of two, since a partition is named by log2(P) rows; device=[0, 1, 2] is a ValueError. An ordinal may repeat, device=[0, 0] putting two partitions on device 0, which runs the exchange on one GPU. "auto" takes devices 0..k for the largest power of two k visible, and is the one-device run when only one is. partition_row_seed= and partition_row_blocks= choose the rows as they do for partitions=, the block count equal to the list length. stats.partition.partitions is the list length, devices the list, and the per-layer lists carry one entry per partition. The call scatters, propagates and gathers every time; a resident multi-device sum is Rust-only (below).

One device per MPI rank

With comm=, device= names this rank’s one device: an ordinal, or "auto" for a device near the rank’s CPUs, distinct per rank on a node while there are enough (see MPI ranks):

evolved = observable.propagate(circuit, policy, direction="heisenberg", comm=comm, device="auto")

Everything else is the host comm= contract: replicated input, collective calls, result="gather" or "local", and stats.partition holding this rank’s entry, with devices this rank’s device. It needs the extension built with both features (maturin develop --release --features cuda,mpi); a build lacking either raises RuntimeError naming it. A list of several ordinals under comm= is a ValueError. A rank that cannot use its device fails the call on every rank, the peers raising the same exception type naming it.

Rust: several devices, and one device per MPI rank

GpuPartitionedSum splits a sum across device partitions, one per listed device, and holds it across calls; GpuPartitionedSum::from_host(&sum, ordinal) is the one-device form:

#![allow(unused)]
fn main() {
use paulistrings::gpu::GpuPartitionedSum;

let mut split = GpuPartitionedSum::scatter_to_devices(&observable, &[0, 1, 2, 3])?;
split.propagate(&circuit, ApproxTopN(10_000_000), Direction::Heisenberg)?;
let evolved = split.gather()?;
}

scatter_to_devices_with takes a GpuScatterOptions { devices, rows, layer }: a runtime resolved from Placement::Devices, the partition rows, and the device layer’s knobs (GpuLayerOptions); GpuScatterOptions::new(&devices) fills in the defaults.

The exchanged columns move device to device. With enable_trace, every layer’s PartitionLayerRecord carries the device partitions’ counters as device (a DeviceLayerRecord each: records, block size, arena batches, fallbacks, path taken) beside the shared refine_passes and receive_chunks. A device driver takes its policy as a BuiltinTruncation, which every builtin and combinator converts into (ApproxTopN(n) or &ApproxTopN(n) alike). A runtime of Placement::Devices { per_device > 1 } puts several partitions on one device, which is how the exchange is tested on a single GPU — give each partition its own GPU in practice; sharing one is a testing configuration, not a performance one. With features mpi and cuda, gpu::MpiGpuSum is the MPI ranks driver with each rank’s share on its own device.

Exchange mode, merge and chunking

Every device group runs one protocol: the block headers and offsets cross the transport, the group agrees to go ahead, and the columns move straight into the receiver’s device memory, by device-to-device copies in process and over NCCL in an MPI group of more than one rank (gpu::MpiGpuSum). The NCCL communicator starts once, collectively, at scatter: a rank that cannot load NCCL, or two ranks sharing a device (which NCCL refuses), fail the scatter on every rank. PAULISTRINGS_WIRE_TIMEOUT_S (default 300) bounds every wait on the exchange, and every collective of an in-process device group (a host group waits 10 s); a timeout or a failed transfer poisons the split (and aborts the NCCL communicator), surfacing as GpuError::Timeout, GpuError::Wire or GpuError::Nccl (RuntimeError from Python). Running it needs the nccl/2.23.4-1 module alongside cuda,mpi — see Installation and MPI ranks.

A device sender merges one partner’s rows by key before the exchange (GpuLayerOptions::premerge), so two remote deltas that land on the same receiver key ship as one row. It runs for a dense two-qubit unitary, whose remote entries share an output support pattern, and never for a Pauli rotation, which has one remote entry, or a Clifford, whose remote entries never share an output pattern. It costs a second fused pass on the sender that a same-device exchange (two virtual partitions on one card) does not repay at low merge ratios; a real deployment, one partition per GPU, is where it pays.

A device receive moves in chunks of destination positions rather than all at once (GpuLayerOptions::exchange_bytes, PAULISTRINGS_GPU_EXCHANGE_BYTES, K/M/G suffixes), capping the receive volume resident on the device at once; unbounded by default, one chunk. The send side is never chunked: its export volume stays resident until the layer’s last chunk moved.

A Clifford layer whose deltas stay inside the partition skips the sort. A Clifford maps every key to exactly one key, so nothing can collide; the device counts and scatters such a layer directly (GpuLayerOptions::clifford off returns it to the fused layer). A Clifford layer that crosses a partition boundary keeps the fused layer, which merges the received rows.

Peer access and NVLink

Before its first copy from device src into device dst, an in-process group (GpuPartitionedSum) lets dst reach src’s memory directly, so the copy goes over NVLink or PCIe peer-to-peer, granting both the driver’s peer-context access and the source’s memory-pool access list a pooled allocation needs to be reachable from a peer at all. On a pair the driver reports as unsupported the copy stages through the host instead of faulting, and the outcome is logged.

What changes on a device

Results agree to floating-point tolerance, not bit for bit. The device sums equal keys in a different order from the host, as a different bucket count does.

Storage order is the host’s canonical order, but for the device’s own bucket count. The download re-sorts each bucket to the host’s lexicographic order, and the device chooses its bucket count by rows per block rather than terms per bucket, so x_array() of a device result can list the terms in a different order from a host run’s. Compare two results by key, never by position.

Device failures are exceptions, not aborts: an exhausted device raises MemoryError, a request the backend does not implement NotImplementedError, a device placement that does not resolve ValueError, and every other device failure RuntimeError. Under comm=, a failure on one rank’s device fails the call on every rank, its peers raising a RuntimeError that names the failing partition and layer.

Limits

• The resident GpuPauliSum is one device. A multi-device or per-rank run from Python scatters and gathers on every call; the Rust API above keeps the split resident.
• Exact topn runs only on one device. truncation.topn matches the host term for term on a lone device=<int>, to_device, or a device= that resolves to one ordinal; a device list of more than one entry or comm= with device= raises NotImplementedError, since the n-th largest of a split sum has no collective form — use truncation.approx_topn(n) there instead.
• Only the built-in policies run on a device. Every truncation factory and its &/| compositions run on the device; a custom Rust TruncationPolicy has no conversion to BuiltinTruncation, so passing one to a device driver does not compile.
• Memory caps the sum at about 5e7 terms per 48 GB card at 128 qubits, since a layer holds its input, its output and a staging arena at once.
• Widths W ≥ 8 (more than 256 qubits) are correct but untuned.
• A multi-device or MPI group holds its exchange in device memory. One export volume and one receive volume stay resident on a partition’s device during a remote layer, on top of its sum, so two virtual partitions on one card can run out of memory at a term count a single device still fits.
• A device partition in a group cannot refine mid-run. It runs every remote layer at the group’s agreed bucket count and raises rather than growing the count when a block or a received segment exceeds the fused kernel’s cap.
• The sender-side merge (premerge) does not pay on a same-device exchange at low merge ratios. Two virtual partitions sharing one card are a testing configuration for that reason; give each partition its own GPU.
• Multi-node runs are unmeasured. Peer copies between the GPUs of one node run at about 94 GB/s over NVLink, and the NCCL exchange has run within one node only.
• One chunk’s transfer does not overlap the fused layer of the chunk before it, and the chunked receive has not run over a real NCCL communicator.

See it in use

python/paulistrings/tests/test_cuda.py runs every policy against the host engine at three widths, and crates/paulistrings/tests/propagate_gpu.rs is the Rust differential net under the same channels.

Measurements

Product states

A propagated sum is read out with expectation(state), which evaluates ⟨s|A|s⟩ for a single-qubit product state |s⟩ in one masked pass over the terms. The state is either a uniform shorthand — "x+", "y+", "z+", the +1 eigenstate of that Pauli on every qubit — or a per-qubit label string of exactly num_qubits characters from 0 1 + - r l, qiskit’s Statevector.from_label alphabet, with 0/1 the Z eigenstates, +/- the X ones and r/l the Y ones. The default is "x+", which is the state most pages on this site read against.

from paulistrings import Circuit, PauliSum, truncation

observable = PauliSum.from_strings(
    {"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4
)
circuit = Circuit(4)
circuit.rz(0.3, 0)
evolved = observable.propagate(circuit, truncation.coeff(1e-10), direction="heisenberg")

print(evolved.expectation("x+"))        # |++++>
print(evolved.expectation("z+").real)   # |0000>
print(evolved.expectation("0+1r").real) # |0>|+>|1>|+i> per qubit

(0.9888341222814014+0j)
0.0
0.25

The return type is Python complex; take .real for a Hermitian observable, and treat a non-negligible imaginary part as a sign the sum is not Hermitian rather than as noise to discard. A label of the wrong length, or a character outside the alphabet, is a ValueError that restates the rule.

The pass is cheap next to the propagation that produced the sum, and the state is the thing that is cheap to vary: one propagation, any number of product states. The observable is not: the propagated object is the observable, so a second observable is a second propagation, not a second pass. The label table is Measurement; First propagation is this readout on a full run.

Which state the label denotes

expectation does not know which direction produced its sum, so the meaning of the label is fixed by the direction you passed to propagate.

Swapping the two is silent: both return a number of the right size and type, and nothing raises. One rx rotation on Z, read against the Y eigenstate "r", shows the two answers differing by a sign:

z0 = PauliSum.from_strings({"Z": 1.0}, num_qubits=1)
rotate = Circuit(1)
rotate.rx(0.7, 0)
print(z0.propagate(rotate, direction="heisenberg").expectation("r").real)
print(z0.propagate(rotate, direction="forward").expectation("r").real)

0.644217687237691
-0.644217687237691

Most physics questions — a state prepared, a circuit run, an observable measured — are the Heisenberg row, and the label is the prepared state. The worked example with both directions on the same circuit is Direction, and the table with the engine’s view of each direction is Direction semantics.

Stabilizer states

A product-state label cannot express an entangled state; expectation_stabilizer(generators) reads the sum against any stabilizer state — Bell, GHZ, cluster, the output of any Clifford circuit — given as a list of exactly num_qubits signed generators. Each is a Pauli string in the from_strings alphabet with an optional sign prefix, "+XX", "-ZZ", or bare "ZZ" for +; the generators must pairwise commute and be independent, and anything else is a ValueError.

bell = PauliSum.from_strings({"YY": 1.0}, num_qubits=2)
print(bell.expectation_stabilizer(["XX", "ZZ"]).real)     # |Φ+>, YY = -XX·ZZ
print(bell.expectation_stabilizer(["-XX", "+ZZ"]).real)   # flip one sign, flip the state

-1.0
1.0

Rather than write generators by hand, interop.stabilizers_from_stim reads them from a Clifford circuit’s output state — a stim circuit, program text, tableau or tableau simulator — in the sign convention this library shares with stim. It imports stim lazily on first call — install the examples or bench extra from Installation to use it. num_qubits= pads the register on the right with +Z generators, leaving those qubits in |0⟩, so a stim circuit on a few qubits can be read against a wider sum.

import stim
from paulistrings import interop

generators = interop.stabilizers_from_stim(stim.Circuit("H 0\nCNOT 0 1"))
print(generators, bell.expectation_stabilizer(generators).real)
print(interop.stabilizers_from_stim(stim.Circuit("H 0\nCNOT 0 1"), num_qubits=3))

['+XX', '+ZZ'] -1.0
['+XXI', '+ZZI', '+IIZ']

The cost is O(terms · num_qubits² / 64) after a one-off O(num_qubits³ / 64) reduction of the generators, so a few times a product-state pass per term and still negligible next to the propagation. This is the split B7 is built on: stim prepares a 36-qubit cluster state in polynomial time, and Pauli propagation handles the non-Clifford tail it cannot touch. Signatures and the generator rules are under Measurement.

Overlaps and the identity coefficient

Two sums on the same register have a Hilbert–Schmidt overlap a.overlap(b) = tr(a† b) / 2^n, which is Σ conj(a_P) b_P over the strings they share, since distinct Pauli strings are orthogonal. Its two everyday uses are the norm Σ|c|² as a.overlap(a), introduced under Inspecting a sum, and comparing two propagated results term for term without decoding either. identity_coefficient() is the coefficient of I…I, or 0 when the sum has no identity term, which equals tr(O) / 2^n and hence the expectation of O in the maximally mixed state.

truncated = observable.propagate(circuit, truncation.coeff(1e-1), direction="heisenberg")
print(evolved.overlap(evolved).real, truncated.overlap(evolved).real)
print(evolved.identity_coefficient())

0.25 0.24454173796592743
0j

A num_qubits mismatch between the two sums is a ValueError; there is no implicit padding. Both signatures are under Measurement.

Echo read-outs

An echo experiment perturbs the evolved observable A by V = ⊗_q exp(-iδ G_q) and reads 2^-n tr(A V† A V). The library has no echo read-out of its own: to_arrays() hands the propagated sum to numpy, and Showcase B8 computes the exact echo and its diagonal approximation there, per rank under MPI.

From a number to a result

A single expectation value at a single cutoff is a number, not a result. The truncation that made the sum tractable deleted terms carrying signs, so the value can sit on either side of the truth and need not move monotonically as the cutoff tightens — No variational bound is the statement and the measured cases. What turns the number into a result is Validating a result: sweep the cutoff, read the retained norm alongside, apply the plateau criterion, and score against an exact reference wherever one is affordable.

See it in use

• First propagation — expectation("x+") after a Heisenberg run, then the same value across three cutoffs with the retained norm alongside.
• B7 — Stabilizer-state preparation — stabilizers_from_stim on a 36-qubit cluster state, expectation_stabilizer on the propagated observable.
• B1 — Operator scrambling — the OTOC and the conserved norm computed from the coefficient arrays rather than an expectation call, the array route of Symplectic arrays.
• B2 — Noisy circuit verification — one Z expectation under four noise channels, and why amplitude damping moves it where the Pauli channels only shrink it.

Examples

Contents

• First propagation — the site’s guided walkthrough: build an observable, build a circuit, propagate, validate.
• Showcases — five measured applications, each checked against an independent oracle.
• Benchmarks — five benchmarks against PauliPropagation.jl and exact references, two with negative results.
• Engine performance — roofline analysis, layer-time breakdown, and partitioned/distributed scaling, all measured.
• Against other tools — where this method fits against state-vector, stabilizer and MPO simulation.

First propagation

Manual: Operators, Engine and propagation

This walks through one Heisenberg-picture propagation end to end: build an observable, build a circuit, propagate it, read out an expectation value.

import math
from paulistrings import Circuit, PauliSum, truncation

# Observable: average X magnetization on 4 qubits. Coefficients multiply the
# literal Hermitian Pauli string — `Y` carries no phase of its own.
observable = PauliSum.from_strings(
    {"XIII": 0.25, "IXII": 0.25, "IIXI": 0.25, "IIIX": 0.25}, num_qubits=4
)

circuit = Circuit(4)
circuit.rz(math.pi / 8, 0)
circuit.cnot(0, 1)
circuit.h(2)

evolved = observable.propagate(
    circuit,
    truncation.coeff(1e-10),      # drop |c| <= 1e-10 after every channel
    direction="heisenberg",       # U† O U — always pass this explicitly
)
print(len(evolved), evolved.expectation("x+").real)

5 0.7309698831278217

Five Pauli strings survived, and the average X magnetization in |++++⟩ after the circuit is 0.731.

Step by step

PauliSum.from_strings builds the observable directly from a dict of Pauli strings to coefficients. Here it is the average X magnetization on 4 qubits: a quarter weight on each single-qubit X string.

Circuit(4) starts an empty 4-qubit circuit. circuit.rz(math.pi / 8, 0) pushes a Z rotation on qubit 0. circuit.cnot(0, 1) pushes a CNOT entangling qubits 0 and 1. circuit.h(2) pushes a Hadamard on qubit 2. Each call appends one gate, so the circuit is exactly the three-gate sequence written above, in that order.

observable.propagate(...) evolves the observable through the circuit. truncation.coeff(1e-10) is the truncation policy: after every channel, drop any term with |c| <= 1e-10. direction="heisenberg" says to compute U† O U, walking the circuit in reverse and applying each channel’s adjoint — pass this explicitly, since the two directions answer different questions.

evolved is the propagated PauliSum. len(evolved) is the number of surviving terms. evolved.expectation("x+") evaluates the expectation value in the uniform product state |+…+⟩; .real drops the (here zero) imaginary part.

Look at the terms that survived

PauliSum has no term-listing accessor: it hands back the symplectic columns the engine stores, and a Pauli character is the pair of bits for that qubit. Qubit q’s Pauli is bit q of the x word and bit q of the z word — (0,0) is I, (1,0) is X, (0,1) is Z, (1,1) is Y. Decode the five terms and sort them by magnitude:

import numpy as np

PAULI = np.array([["I", "Z"], ["X", "Y"]])

x = evolved.x_array()
z = evolved.z_array()
coefficients = evolved.coefficients_array()

for row in np.argsort(-np.abs(coefficients)):
    label = "".join(PAULI[(x[row, 0] >> q) & 1, (z[row, 0] >> q) & 1] for q in range(4))
    print(f"{label}  {coefficients[row].real:+.9f}")

IIZI  +0.250000000
IXII  +0.250000000
IIIX  +0.250000000
XXII  +0.230969883
YXII  -0.095670858

Two of the four seed terms, IXII and IIIX, came through untouched, on the qubits the cnot and the rz never reach. IIXI picked up circuit.h(2) and became IIZI, since a Hadamard swaps X and Z on the qubit it acts on. XIII is the one the cnot and the rz both acted on, since qubit 0 is the CNOT’s control: cnot(0, 1) turns it into XXII, and rz(math.pi / 8, 0) then splits that into XXII and YXII. The [row, 0] index picks word 0 of the term, which is the whole key here because 4 qubits fit in one 64-bit word — see Symplectic arrays for the general layout.

Check the cutoff

truncation.coeff(1e-10) was a guess. Run the same propagation at a few looser cutoffs and watch the term count, the answer and the retained norm Σ|c|² together:

norm_in = float(np.sum(np.abs(observable.coefficients_array()) ** 2))

for eps in [1e-1, 1e-2, 1e-10]:
    trial = observable.propagate(circuit, truncation.coeff(eps), direction="heisenberg")
    retained = float(np.sum(np.abs(trial.coefficients_array()) ** 2)) / norm_in
    print(f"{eps:7.0e}  terms={len(trial):3d}  <O>={trial.expectation('x+').real:+.9f}  retained={retained:.6f}")

  1e-01  terms=  4  <O>=+0.730969883  retained=0.963388
  1e-02  terms=  5  <O>=+0.730969883  retained=1.000000
  1e-10  terms=  5  <O>=+0.730969883  retained=1.000000

The answer is the same to every printed digit at all three cutoffs, and 1e-2 already retains the whole operator — this circuit is small enough that the truncation never bit.

Read the top row carefully, though. At 1e-1 the sum lost 3.7% of its norm and the answer did not move at all, because the deleted YXII term has a Y on qubit 0 and contributes nothing to ⟨++++|O|++++⟩. A flat answer over a changing sum is not evidence of convergence, and on a real problem that is exactly the trap to check for — Validating a result is the full recipe.

Where to go next

• Manual — operators, circuits, engine and propagation, measurements, read start to finish or by topic.
• Validating a result — the cutoff sweep, the retained norm, and when a number may be quoted.
• Direction and Direction semantics — what "heisenberg" and "forward" compute, and which state expectation is then read against.
• Measurements — reading a sum against a stabilizer state or another sum, not just a product state.
• Incremental propagation — the same walkthrough, stepped out into a Trotterized time series.
• Against other tools — whether Pauli propagation is the right method for your problem at all.
• examples/common/observables.py and examples/common/circuits.py — more observable and circuit constructions.

Showcases

Five measured applications.

B3 (variational pre-training) and B4 (QML/QCNN) are not part of this suite.

Contents common to every showcase page

A convergence panel on every truncated result: a single number from a single cutoff is not a result here. The retained Hilbert–Schmidt norm N = Σ|c_P|² alongside it, conserved under exact unitary evolution and equal to 1 for a single Pauli seed, so 1 − N is exactly the deleted fraction of the operator under truncation. Named dependencies rather than silent approximations: where a reference was not reachable, the page says so and what it would cost. Validating a result is that panel as a recipe, for your own runs.

Reproducing a showcase

./scripts/setup.sh                # .venv with the examples extra
source .venv/bin/activate
uv sync --extra examples          # after any Rust change: rebuilds the extension

RAYON_NUM_THREADS=1 python examples/b1_operator_scrambling/run_b1_1d.py
RAYON_NUM_THREADS=1 python examples/b2_noisy_verification/run_b2.py
RAYON_NUM_THREADS=1 python examples/b5_operator_backpropagation/run_b5.py
RAYON_NUM_THREADS=1 python examples/b6_resource_probes/run_b6.py
RAYON_NUM_THREADS=1 python examples/b7_stabilizer_prep/run_b7.py
RAYON_NUM_THREADS=8 python examples/b8_ole/run_b8.py --alpha 0.15 --L 6 --cache 5e6 --seeds 0:4 --out results/b8

Each script rewrites every figure and JSON file next to itself. Each showcase also has a CI-visible correctness gate under examples/tests/ that runs in about a second on numpy alone — the physics is checked on every commit even though the full runs are manual.

Source: examples/README.md.

B1 — Operator scrambling

Manual: Truncation, Direction

Heisenberg-evolve a single-site Z through a kicked Ising circuit and watch the operator spread. The light cone, OTOC, butterfly velocity, and two-point function are read off the evolved Pauli sum in one pass, checked against an independent dense 2ⁿ×2ⁿ construction.

Operator scrambling and the question

A local operator evolved backward through a chaotic circuit, O(t) = U(t)† O U(t), spreads across a growing set of sites: a light cone whose edge advances at a butterfly velocity v_B. The out-of-time-order correlator (OTOC) C(r,t) tracks how much of the operator anticommutes with a probe at site r, approaching C → 1 as the signature of full scrambling. In the Pauli basis this is all one object, O(t) = Σ_P c_P(t) P, so support, OTOC, and two-point function are three sums over the same coefficients — the question is whether truncation converges to a single, contour-independent v_B, and whether a weight cap, effective elsewhere, survives a 2D quench.

Running it

source .venv/bin/activate
python examples/b1_operator_scrambling/run_b1_1d.py     # 1D chain, a few minutes
python examples/b1_operator_scrambling/run_b1_2d.py     # 2D quench, ~15 minutes

Evolution is one propagate call per Trotter step, truncating after each — an excerpt from the actual example script, verified by pytest examples/tests, not standalone:

from paulistrings import truncation
from common import observables

evolved = observables.single_z(center, n).propagate(
    step_circuit, truncation.coeff(min_abs_coeff), direction="heisenberg"
)

Every quantity below is read from three arrays exported once — evolved.x_array(), evolved.z_array(), evolved.coefficients_array(): support w_q(t) sums |c|² per site from the bit columns (plotted against (site, t) this is the light cone); the OTOC C(r,t) = 2·Σ_{P anti W_r}|c_P|² tests one symplectic bit per site (X_r→z_r(P), Z_r→x_r(P), Y_r→x_r(P)≠z_r(P)); the two-point function G(r,t) = ⟨Z_r,O(t)⟩ is the coefficient of the weight-one string Z_r, found by popcount. N(t) = Σ_P|c_P(t)|² = ⟨O(t),O(t)⟩ equals 1 for the seed and is conserved exactly. Truncation only lowers N(t), and 1 − N(t) is exactly the norm discarded — the retained-norm diagnostic behind every curve below.

The 1D chain

Open n = 61 chain, kicked-Ising entangler θ_zz = −π/2, kick angle θ_h = 0.9, seed Z_30, 12 Trotter steps, Heisenberg direction. Because this engine truncates after every channel, t successive calls on the one-step circuit reproduce exactly the (apply-adjoint, truncate) sequence of one call on the t-step circuit, so the whole time series costs one 12-step propagation.

At θ_h = π/2 every gate is Clifford, so one Pauli string evolves into one Pauli string with unit coefficient: its support is the exact causal cone, measured rather than assumed. Over steps 1–8 the sum stays at exactly one string, weight running 1, 3, 5, …, 15 and cone radius 0, 1, 2, …, 7 — 1.000 sites per Trotter step, the structural bound of a single-site X layer followed by one commuting ZZ layer moving the boundary by at most one bond.

Cutoff breakdown

Four cutoffs, twelve steps. At the final step:

At 10⁻³ the sum shrinks after step 9 (59 140 → 11 900 terms) while N collapses to 2.7% — meaningless past step ~8, and N(t) is the only signal that says so. At 3·10⁻⁶ the retained norm stays above 0.999 through step 9.

The OTOC

C(r, t) at r = centre + 5, t = 12, as the cutoff tightens: 0.031137 (10⁻³), 1.003293 (10⁻⁴), 1.346556 (10⁻⁵), 1.368813 (3·10⁻⁶) — the loosest cutoff wrong by 40×, the two tightest agreeing to 2.2·10⁻² (1.6% relative). The cone interior saturates at C ≈ 1.3–1.4, heading for C → 1.

Butterfly velocity

Fitting front distance against step (steps 3–12), for every cutoff and contour level:

Verdict at θ_h = 0.9: v_B = 1.000 sites per Trotter step, converged — the structural bound measured above. The w > 10⁻² column is the cautionary tale: still climbing (0.776 → 0.848 → 0.903), it reads v_B ≈ 0.9 to anyone who picks one contour and stops.

Scanning the kick angle at two cutoffs (n = 61, 12 steps, 27 s total, w > 10⁻⁴ contour, min_abs_coeff = 10⁻⁵) gives a velocity rising with kick strength and saturating the causal bound, 0.327 → 0.503 → 0.776 → 1.000 for θ_h = 0.2, 0.4, 0.6, 0.9. Only the first two are converged (N(12) = 1.0000); at θ_h = 0.6 and 0.9 velocity moves by 0.05 between cutoffs (N(12) = 0.9995 and 0.7158).

The 2D quench

Open 7×7 lattice (49 sites), seed Z at the centre, and a physical Trotter step of H = J Σ Z_iZ_j + h Σ X_i with J = h = 1, dt = 0.15, up to t = 1.5 — physical time, not a Floquet period, so the step count resolves the dynamics rather than merely advancing it.

2D as the hard case

In 1D the causal cone holds O(t) sites; in 2D it is an area, O(t²), and Pauli strings inside it grow as 4^O(t²). A weight cap does not fix this: on a degree-4 lattice one exp(iπ/4 Z_iZ_j) layer takes weight-1 to weight 5, so max_weight = 4 gives N = 0.386 at step 2, 0.058 at step 4 (sum never above two strings); max_weight = 12 holds N = 1.0 through step 3, then 0.523 — the cap deletes the dynamics, not a tail.

So the 2D sweep uses min_abs_coeff only, and each cutoff runs until a term ceiling (1.2·10⁸ stored terms) stops it, so the converged window is measured rather than assumed:

A cutoff loose enough to reach t = 1.5 has thrown away 19% of the operator; one tight enough to keep 99.9997% of it hits the term ceiling at t = 0.9.

Magnetization

Verdict: the 2D quench magnetization is converged to ~10⁻³ out to t = 0.90 (6 Trotter steps), and not converged beyond it. The two tightest cutoffs differ by 0, 0, 6.5·10⁻⁵, 1.3·10⁻⁴, 3.1·10⁻⁴ and 1.2·10⁻³ at steps 1–6; past step 6 only the two loosest cutoffs reach at all, disagreeing by up to 3.8·10⁻² — shown in the figure’s right-hand panel.

The maps are diamonds, not squares — a Manhattan ball, since one Trotter step moves the boundary by one bond. Mean per-site weight at step 6, min_abs_coeff = 10⁻⁷, by graph distance 0…6: 0.985, 0.504, 0.0402, 4.3·10⁻⁴, 8.4·10⁻⁷, 1.2·10⁻¹¹, 0 — one to two orders of magnitude per site.

The two-point function is a sensitivity floor, not an error: a coefficient G(r,t) below min_abs_coeff reads back as exactly 0.0. At step 6 the only sites with |G| > 0 are the centre (G = 0.3477) and the four diagonal distance-2 sites (G = 5.235·10⁻⁵ each, as symmetry requires); every other site reads exactly zero. A tighter 10⁻¹⁴ probe on a 5×5 lattice shows the centre alone through step 3, so at least part of this is a real selection rule, unresolved here.

Performance

Wall times above are 32-thread on the default Rayon pool, taken on a shared, loaded workstation (~2× an idle host for the 2D sweep); nothing here is a benchmark claim. Peak resident memory: 19 GB for the 2D script; the 1D script’s tightest step holds 2.3·10⁸ strings, ~7 GB in the engine plus about the same again in the numpy export.

3D scalability. On a 3×3×3 cubic lattice (27 sites, 54 bonds), same quench and term ceiling as the 2D sweep, the converged cutoff (10⁻⁶) reaches t = 0.90 with 4.9·10⁸ strings — 88 s for that step, growing ~6× per step; step 7 projects to ~3·10⁹ strings (≈100 GB) and ~10 minutes. Truncation error is directly visible: ⟨Z_c⟩ must lie in [−1, 1], and the looser 10⁻⁵ cutoff reports 1.010 and 1.088 at steps 7–8, exactly where N has fallen to 0.97 and 0.91. This 27-site pilot is a lower bound — its centre is 3 bonds from the farthest corner, so by step 4 the cone already covers the whole lattice.

Validation

scrambling.dense_* builds U(t) as an explicit 2ⁿ×2ⁿ matrix by Kronecker products and evaluates the same quantities with dense linear algebra — no PauliSum, no engine, no qiskit. At n = 9, 4 Trotter steps, untruncated (1430 terms), against a 10⁻¹⁰ bar, the worst gap between engine and dense is 6.7·10⁻¹⁶ on every coefficient, 3.1·10⁻¹⁵ on ⟨O,O⟩, 4.0·10⁻¹⁵ on the support profile, 5.8/5.3/3.8·10⁻¹⁵ on the three OTOC probes, 4.9·10⁻¹⁷ on the two-point function, and 2.2·10⁻¹⁶ on the probe-average identity. The 2D script repeats this on a 3×3 lattice with the quench magnetization added: worst gap 2.1·10⁻¹⁴.

At min_abs_coeff = 0.05 the truncated leg keeps 54 of 1430 terms, discards 0.318 of the norm, and its worst per-site weight error is 0.280 — a ratio of 0.88, licensing 1 − N(t) as the error proxy used throughout: a calibration, not a theorem, so the script asserts a factor-of-10 band.

The CI-visible gate (pytest examples/tests/test_showcase_b1.py) is numpy-only and runs in under a second. Reference host: Intel Xeon Gold 6244 @ 3.60 GHz (ccqlin038), rustc 1.94.0, Python 3.11.11.

Numbers: computed by run_b1_1d.py / run_b1_2d.py, with full per-(cutoff, step) records in results_1d.json / results_2d.json and the headline table in the README. No number here is from literature.

B2 — Noisy circuit verification

Manual: Noise channels, Truncation

The 127-qubit heavy-hex kicked-Ising circuit from Benchmark C gets per-gate depolarizing noise added, and the noise makes the simulation cheaper. At p = 3e-2 the tracked set peaks at 651× fewer terms than at p = 0, and the run finishes 1078× faster, at the same cutoff. Every density-matrix method scales the opposite way.

Tracked-set shrinkage from noise

A Kraus/density-matrix simulator pays for noise twice: it carries a 4ⁿ object instead of a 2ⁿ one, and each channel costs more than a gate. Pauli propagation pays negatively: a single-qubit depolarizing channel rescales a coefficient by 1 − 4p/3 per non-identity qubit on its support, so a weight-w string crossing one noise layer loses (1 − 4p/3)^w, and after d layers a string that held weight w throughout is down by (1 − 4p/3)^{wd}. The filter is exponential in weight, so a fixed min_abs_coeff becomes an effective weight cap that tightens with depth, and the tracked set of a scrambling circuit is dominated by its high-weight tail. One channel runs on every qubit in the support of the gate that just ran, 686 channels per step at n = 127, 13 720 across 20 steps. p = 0 runs as depolarize(0.0) rather than a separate circuit, so the rescale is exactly 1.0 and the p = 0 leg is the same computation as the noiseless one, on the noisy legs’ channel schedule — pinned term for term by a test.

Running it

source .venv/bin/activate
RAYON_NUM_THREADS=1 python examples/b2_noisy_verification/run_b2.py           # ~26 min
RAYON_NUM_THREADS=1 python examples/b2_noisy_verification/run_b2.py --quick   # ~1 s, 20 qubits, writes nothing
python examples/b2_noisy_verification/run_b2.py --figures-only                # re-render the SVGs
pytest examples/tests/test_showcase_b2.py                          # the CI gate

RAYON_NUM_THREADS=1 must be exported before the interpreter starts; the driver refuses to run otherwise. The driver adds one noise call per gate (circuit.depolarize/dephase/amplitude_damping/pauli_channel) and propagates with harness.run_propagation(circuit, observable, TruncationSpec(min_abs_coeff=eps), direction="heisenberg").

Terms and time across the noise sweep

θ_h = 5π/16, 20 Trotter steps, Z_62, Heisenberg, contracted against |0…0⟩, single-threaded, fixed min_abs_coeff = 2⁻¹⁴:

Maximum surviving weight falls 19, 19, 14, 13, 0 and mean weight 10.3, 9.8, 8.6, 7.4, 0, the exponential-in-weight filter acting on the tail directly. Wall time falls faster than the term count (1078× against 651×) because the merge phase sorts a resident sum far out of cache at the larger sizes, so a smaller tracked set pays twice. A density-matrix method at n = 127 is not on the table (4^127 amplitudes), so the comparison is within this engine: p = 0 against the noisy legs in the same sweep.

Convergence across cutoffs

At p = 3e-2 the fixed-cutoff answer above reads exactly 0; the converged answer is 7.8·10⁻⁵.

A fixed cutoff is not a fixed accuracy: as p grows the signal shrinks, so the cutoff must tighten with it, and at p = 3e-2 a 2⁻¹⁴ cutoff sits above every contributing coefficient. Noise buys the cutoff reach that fixes this cheaply: the whole p = 3e-2 sweep to 2⁻²⁰ costs 76 s, an eighth of the single p = 0 point at 2⁻¹⁴, resolving the answer to 4.2·10⁻⁷ against 6.5·10⁻⁴ at p = 0; at p = 0 the next dyadic pair is out of reach, projecting ~1.9·10⁸ terms and ~2 h single-threaded. The plateau verdict reuses Benchmark B’s criterion unchanged, and rejects p = 1e-3: its last difference is 2.61·10⁻⁴ but its previous one is 1.03·10⁻³, just over tolerance, so that leg is unresolved rather than quoted.

The same collapse, three other channels

Same circuit, same 2⁻¹⁴ cutoff, same 20 steps, with the depolarizing channel swapped for each of the others the bindings expose:

All four collapse the tracked set by one to two orders of magnitude against the p = 0 baseline’s 1.44·10⁷ peak, a property of Pauli-basis noise generally, not of depolarizing noise specifically. amplitude_damping moves ⟨Z₆₂⟩ upward, to +0.125, since it relaxes toward |0⟩, the +1 eigenstate of Z, and is the only channel here neither self-adjoint nor key-preserving, its Heisenberg dual Z → (1−γ)Z + γI fanning out — checked against the comparisons page’s cross-engine noise-channel parity results. dephase(1e-2) damps ⟨Z₆₂⟩ less since dephasing leaves Z untouched.

Claimable results

At p = 0 the noise model must vanish and this showcase must return Benchmark C’s claimable numbers, read from C’s committed summary.json at run time rather than transcribed: the two 5-step configurations reproduce C’s exact causal-cone references (a 19-qubit reduction cross-checked between an Aer statevector and an untruncated Pauli propagation) to gaps of 3.9·10⁻¹⁵ and 2.2·10⁻¹⁶ against a 0.01 bar, and the 20-step configuration reproduces C’s self-converged 20-step run to a gap of exactly 0.0. Against the exact 5-step reference, ⟨Z₆₂⟩ departs from it linearly in p over two decades (2.2e-16 at p = 0, 1.1e-5 at 1e-6, 1.1e-3 at 1e-4) while the term count holds at 2 146 424 throughout, since at 5 steps the sum has already saturated.

The circuit is ingested, not reconstructed: the 144-edge Eagle coupling map comes from a checked-in file generated from qiskit-ibm-runtime’s FakeSherbrooke, and Z_62 is loaded through a provenance-tagged JSON of published operator supports. The noiseless configuration is claimable at θ_h = 5π/16, 20 steps (⟨Z₆₂⟩ = +0.016131374386, converged from 2⁻⁸ to 2⁻¹⁴ with a last difference of 6.50·10⁻⁴, agreeing with Benchmark C’s independent reference to 0.0), and so are the noisy configurations whose sweeps plateaued, the ones a hardware comparison would actually want, since a device is noisy:

p = 1e-3 at 20 steps is not claimable, since its cutoff sweep never plateaued, and neither is θ_h = 7π/32 at 20 steps noiseless: it is Benchmark C’s measured reachability boundary.

Reachability boundary, with and without noise

Running that 7π/32, 20-step circuit noisily at p = 1e-2 shows what noise does to the boundary. At the same 2⁻¹⁶ cutoff, the noiseless run (Benchmark C) peaks at 45.4 million terms and ends at 38.8 million, last difference 1.44e-1; the noisy run peaks at 0.94 million, ends at 74 717, last difference 9.32e-4 — 48× fewer peak terms, 519× fewer resident at the end, 154× smaller last difference. It still does not pass, since the second-to-last difference is 1.08·10⁻³, over the 1e-3 bar, so converged = false. ⟨Z₆₂⟩ under the noisy channel is also a different quantity from the unitary-circuit value, agreeing only as p → 0, but it is the quantity a hardware verification claim actually needs.

Performance

Peak resident terms range from 14 396 463 at p = 0 down to 22 105 at p = 3e-2, and wall time from 531.0 s down to 0.49 s. Process RSS peaks at 1.25 GiB across the full driver run, which covers 31 records (the noise grid, the channel variants, the noiseless-limit and reachability legs, and the convergence sweeps) in 25.7 min single-threaded end to end, on ccqlin038, rustc 1.94.0, Python 3.11.11.

Validation: an independent dense noisy reference

Every oracle in the suite refuses a noise channel by construction (the statevector oracle because it is unitary-only, the stim oracle because a tableau simulation samples one Pauli error rather than averaging over them), so B2’s correctness gate hand-rolls a dense reference instead: gates applied to an explicit 2ⁿ × 2ⁿ density matrix, forward from |0…0⟩⟨0…0|, unitaries by conjugation and channels as Σ_k K_k ρ K_k† with Kraus operators from each channel’s definition, not the engine’s dual. Two-qubit Cliffords are sums of embedded Pauli products (CNOT = (I + Z_c + X_t − Z_c X_t)/2, SWAP = (II + XX + YY + ZZ)/2), and with no noise gate the reference must reproduce the engine’s exact untruncated propagation. At 1e-10, for all five channels mixed with h/s/cnot/cz/swap, non-Clifford rx/ry/rz, and a weight-two XY rotation:

29 tests, ~1 s, 97 MiB peak RSS, numpy-only.

Caveats

A truncated Pauli sum has no variational bound, since discarded terms carry signs and the error need not shrink monotonically with the cutoff; converged = false means no usable uncertainty estimate, not a slightly weaker one. The noise model is a model, not a device: a real Eagle device has different error rates, coherent errors, crosstalk, and readout error, none of which is here. Peak, not final, term count is the cost, since the sum grows mid-circuit and collapses by the end.

Numbers: raw records in examples/b2_noisy_verification/results.json, verdicts in summary.json next to it, source and provenance in examples/b2_noisy_verification/README.md.

B5 — Hybrid depth reduction

Manual: Incremental propagation, Truncation

Split a circuit near the end, back-propagate the observable through the tail classically, and hand a QPU the shorter front circuit plus an evolved observable. The composed expectation is exactly the full-circuit one — so the trade is QPU depth against classical term count, and this page measures both sides of it.

As the tail depth k grows the residual front circuit shrinks linearly in Trotter steps (186 → 0 gates) while the weight-capped evolved observable grows roughly exponentially (1 → 12 413 terms).

The idea

Split a circuit at layer k from the end:

|0...0> --[ front circuit, depth L-k ]--[ tail circuit, depth k ]--> measure O

Back-propagate O through the tail classically — direction="heisenberg" on the last k layers, with the same engine this repository ships for everything else — to get O' = U_tail† O U_tail. Composing the halves, ⟨psi_front| O' |psi_front⟩ (with |psi_front⟩ = U_front |0...0⟩) equals ⟨ψ|O|ψ⟩ for the full circuit exactly, since Heisenberg conjugation composes. The QPU only runs the shorter front circuit; O' is computed once, off the QPU.

The artifact is a schema-v1 task file — the residual front circuit plus the evolved observable, ready for a QPU-side runner. task_exact.json and task_truncated.json in the showcase directory are exactly that, and they are the same schema the cross-engine comparison drives both engines from.

Running it

source .venv/bin/activate
python examples/b5_operator_backpropagation/run_b5.py

Key calls — an excerpt from the actual example script, verified by pytest examples/tests/test_showcase_b5.py, not standalone:

tail_evolved = observable.propagate(tail_circuit, policy, direction="heisenberg")
psio.save(npz_path, tail_evolved)                 # round-trip through disk
loaded = interop.load_task(task_path)
composed = loaded.observable.propagate(loaded.circuit, loaded.truncation,
                                        direction=loaded.direction)

This regenerates every artifact in the directory — both task JSONs, the .npz, the CSV, both SVG figures — and the whole sweep runs in about 3 seconds on a laptop-class machine; this page is about depth reduction, not speed. The CI-visible correctness gate is numpy-only and runs well under a second: pytest examples/tests/test_showcase_b5.py.

Validation

8-qubit, 4-layer seeded hardware_efficient_ansatz, observable Z_4, tail depth k = 1 (front: 3 layers / 69 gates; tail: 1 layer / 23 gates):

The round-trip gap is exactly 0.0, well inside the 1e-12 bound, and it is a genuine round trip through disk: the evolved observable is saved to .npz via paulistrings.io and read back before being embedded in the task file.

Sweeping min_abs_coeff against the exact reference above:

Truncation is applied after every channel, not after every Python call (see Truncation), so splitting a circuit and truncating separately on each half agrees exactly, for any split point, with truncating the full circuit in one shot — a test pins it. The table above is therefore independent of where k is chosen. The gap is not monotone in the cutoff (2.9e-2 at 3e-2, then 6.3e-2 at 1e-2) — dropped terms carry signs, and a truncated Pauli sum has no variational bound.

The actual trade-off: depth against term count

16-qubit heavy-hex kicked-Ising sublattice, 6 Trotter steps, θ_h = 0.6 (a generic, non-Clifford kick), observable Z_8, truncated at weight <= 6 throughout — the weight cap is what makes the classical half tractable in practice. Sweeping the tail depth k from 0 (nothing back-propagated) to 6 (front empty):

The residual circuit shrinks linearly in Trotter steps while the weight-capped evolved-observable term count grows roughly exponentially before the cap bites — the classic operator-spreading-vs-light-cone shape, and exactly the cost paid classically for every layer moved off the QPU.

Numbers: every value on this page comes from examples/b5_operator_backpropagation/run_b5.py, committed alongside its outputs — task_exact.json, task_truncated.json, evolved_observable_exact.npz, depth_vs_terms.csv, depth_vs_terms.svg, and convergence_panel.svg — in the showcase directory.

B6 — Resource probes

Manual: Cost is terms, not qubits, Inside the engine

Two diagnostics read off the evolved sum's numpy export, each answering "how hard is this operator?" under a different cost model: Pauli-spectrum entropy (a magic-adjacent diagnostic for truncation-based Pauli propagation) and operator entanglement, the cost model for matrix-product-operator methods. Both are zero on a single Pauli string.

Exact (untruncated) sweep of the kick angle on a 16-qubit kicked-Ising chain. Both diagnostics vanish at both Clifford endpoints and are strictly positive at all 15 interior angles. Dotted verticals mark the Clifford angles.

The two quantities

All entropies are in nats, computed in pure Python over the symplectic arrays PauliSum.x_array() / z_array() / coefficients_array() — no core additions.

Pauli-spectrum entropy

This is the operator quantity, not the state stabilizer Rényi entropy (SRE) of Leone, Oliviero and Hamma, built on a pure state and normalized to vanish on stabilizer states — the quantity for which Leone and Bittel proved α ≥ 2 a magic monotone, for pure states only. Neither result transfers here, and nothing in this showcase is presented as a magic monotone.

For O = Σ_P c_P P in this repository’s Hermitian convention, p_P = |c_P|²/Σ|c|² is a probability distribution over the Pauli basis, S_α its Rényi-α entropy, L the linear (purity) form of α=2, and S_2 = −ln(1 − L) identically. Shao, Cheng and Liu name the operator-side quantity the Operator Stabilizer Rényi entropy (OSE) and prove it governs Pauli-propagation truncation error and the Top-K budget for a target accuracy; their c_i is this repository’s c_P, and when Σc_i² = 1 (true for a single Pauli string, preserved by exact unitary Heisenberg evolution), their definition is algebraically identical to the renormalized S_α here. Truncation breaks Σ|c|² = 1 by (α/(1−α)) ln(Σ|c|²); this showcase always renormalizes and reports the surviving weight Σ|c|² in every table, so the difference is never hidden.

Zero on a single Pauli string, invariant under Clifford conjugation, raised by any non-Clifford rotation, additive over tensor factors, and basis-dependent by construction — a cost model for this representation, not a basis-independent resource measure.

Operator entanglement

Every Pauli string factorizes across a cut, P = P_A ⊗ P_B, so the coefficient vector reshapes into M[a,b] = c_P, and M’s SVD is the operator Schmidt decomposition: Zanardi’s operator entanglement, whose entropy is the operator space entanglement entropy of Prosen and Pižorn, introduced there as the statement that simulating observables, unlike typical states, is efficient for initially local operators. Unlike the Pauli-spectrum entropy it is not Clifford-invariant (a CNOT across the cut raises it), but a Clifford circuit still maps a single Pauli string to a single Pauli string, whose operator entanglement is zero for any cut.

M is dense-allocated and SVD’d: O(n_left·n_right·min(n_left,n_right)) time and 16·n_left·n_right bytes, with n_left·n_right able to approach T² in the worst case. The routine refuses past max_entries (4·10⁶ entries = 64 MiB), and every table below prints the Schmidt shape; the largest across the depth sweep was 1027 × 3328 (3.4 M entries).

Dense cross-check

At each size the dense 2ⁿ × 2ⁿ matrix is rebuilt with numpy.kron, and both diagnostics are recomputed by routes sharing no code with the array-based probes — brute force over all 4ⁿ traces for the Pauli spectrum, SVD of the reshaped dense matrix for the Schmidt spectrum. Bound: 1e-10.

Every gap is at or below 8.9·10⁻¹⁶, twelve orders inside the bound; hs_weight is a fourth identity, Σ_P|c_P|² = tr(O†O)/2ⁿ, confirming that untruncated unitary evolution preserves spectral weight exactly. Depth and cut vary on purpose — at fixed depth the light cone makes the operator identical at n = 6, 8, 10, so a fixed (depth, cut) would give three oracles on one number. The exhaustive 4ⁿ spectrum runs at n = 6 only (0.6 s; 99 s at n = 8, refused above it). At the two Clifford points (θ_h = 0, θ_h = π/2) a single-Pauli seed must stay a single Pauli string, and it does to floating-point dust: both diagnostics land at or below 5.4·10⁻³⁰ against bounds of 0 and 1e-25.

Kick-angle sweep

1D kicked-Ising chain, n = 16, 5 Trotter steps, seed Z_8, bipartition [0,8) | [8,16) — crossing exactly one lattice bond (a chain, not a heavy-hex sublattice, keeps the cut count from confounding the reading). policy=None throughout: exact, nothing truncated, no convergence panel owed.

(Nine of the seventeen measured angles; full table in theta_sweep.csv.)

• Both diagnostics vanish at both Clifford endpoints and are strictly positive at all 15 interior angles; term count does not follow the same pattern — 16 796 terms at θ_h → 0⁺, 4181 at π/2, 1 exactly at 0. Stored-term count tracks how coefficients round; entropy tracks where the weight actually is.
• The interior is non-monotone: a first local maximum at θ_h ≈ 0.687, a dip near 0.884, and the global maximum at θ_h = 3π/8 (S_2 = 4.68 nats ≈ 108 effective Pauli strings out of 16 796 stored). S_op never exceeds 1.40 nats while S_2 reaches 4.68 — different cost models, same operator.
• S_op^(2) ≤ S_op and S_2 ≤ S_1 at every point, as Rényi entropies must be.

Depth sweep: exact against truncated

Same chain at n = 20, generic θ_h = 0.6, seed Z_10, one cut bond. Exact through depth 6 (208 012 terms, a 462×1715 Schmidt matrix); depth 7 exact is 2.67 M terms and past the guard, so it is truncated-only.

• S_2 grows steadily with depth (0.57 → 4.61 nats, i.e. 1.8 → 100 effective Pauli strings), while stored terms grow far faster (2 → 2.67 M exact): propagation cost tracks term count, achievable accuracy at a budget tracks S_2.
• S_op saturates around 1.3 nats from depth 5 on, even dipping at depth 7 — operator entanglement across a fixed 1D cut is generated only by the gates crossing that cut, so it need not grow the way the Pauli spectrum does. The same operator gets steadily harder for one method and not for the other.
• Truncation at min_abs_coeff = 1e-6 reproduces every exact value to five decimals at depth 6, keeping 6.6× fewer terms with 0.999999996 of the spectral weight surviving.

Truncation convergence

Depth 6 converges monotonically against the exact value over five orders of min_abs_coeff, ending at 4.6·10⁻¹¹ (S_2) and 3.7·10⁻¹¹ (S_op); at 1e-7 the truncated sum already has 4.5× fewer terms than exact while agreeing to eleven digits. Depth 7 has no exact reference, so it is shown self-converging instead: successive drift falls by roughly 1.5 orders per decade of cutoff, converging the quoted depth-7 values to about 10⁻⁷. Both statements are asserted by the script, not eyeballed. Full tables are in the source README.

Performance

Reproducing the whole showcase (both CSVs, the cross-check JSON, and all three SVGs) takes well under a minute. That is not a performance claim about the engine; it does depend on one control, though: the script pins OMP/OPENBLAS/MKL_NUM_THREADS=1 before importing numpy, because with LAPACK left to spawn its own pool on a busy shared host the same 462×1715 SVD was observed at 56 s instead of 0.14 s.

source .venv/bin/activate
python examples/b6_resource_probes/run_b6.py

The CI gate is 18 tests, numpy-only, under a second:

pytest examples/tests/test_showcase_b6.py

References

1. L. Leone, S. F. E. Oliviero, A. Hamma, “Stabilizer Rényi Entropy”, Phys. Rev. Lett. 128, 050402 (2022); arXiv:2106.12587. State SRE — the quantity this showcase’s Pauli-spectrum entropy resembles but is not.
2. L. Leone, L. Bittel, “Stabilizer entropies are monotones for magic-state resource theory”, Phys. Rev. A 110, L040403 (2024); arXiv:2404.11652. Monotonicity for α ≥ 2, restricted to pure states.
3. Y. Shao, S. Cheng, Z. Liu, “Characterizing Pauli Propagation via Operator Complexity”, arXiv:2510.22311 (2025). Operator Stabilizer Rényi entropy (OSE), Definition 1; the truncation-error and Top-K budget bounds.
4. P. Zanardi, “Entanglement of quantum evolutions”, Phys. Rev. A 63, 040304(R) (2001). Operator Schmidt decomposition / operator entanglement.
5. T. Prosen, I. Pižorn, “Operator space entanglement entropy in transverse Ising chain”, Phys. Rev. A 76, 032316 (2007). OSEE as the cost model for simulating observables.

Numbers: raw sweeps in examples/b6_resource_probes/theta_sweep.csv / depth_sweep.csv, cross-check in exact_cross_check.json, full convergence tables and reproduction detail in examples/b6_resource_probes/README.md.

B7 — Stabilizer-state preparation

Manual: Stabilizer states

stim prepares a 36-qubit 2D cluster state as a stabilizer tableau; a non-Clifford tail is propagated in the Heisenberg picture; the expectation is contracted against the stabilizer state at O(m·n²/64) cost, never building the 1.0 TiB state vector the same check would need as a dense array.

The idea

stim runs a Clifford circuit of any depth in polynomial time and cannot touch a non-Clifford gate; Pauli propagation eats non-Clifford rotations for a living and gains nothing from re-deriving a Clifford circuit’s output. Chaining them across the observable gets both:

|0...0> --[ deep Clifford prep U_C ]--> |psi>   --[ non-Clifford tail U_T ]--> measure O

In the Heisenberg picture the tail acts on the observable and the state is only needed at the very end, ⟨ψ|U_T† O U_T|ψ⟩ with |ψ⟩ = U_C|0...0⟩. A stabilizer state is fully described by its n signed Pauli generators s_i·G_i, and for a Hermitian Pauli string P

⟨ψ|P|ψ⟩ = σ    if σ·P is in the group ⟨s_i G_i⟩ for some σ = ±1
        = 0    otherwise

— a group-membership test, O(n²/64) word operations per term after an O(n³/64) GF(2) echelon reduction of the generators. Contracting the evolved m-term sum is O(m·n²/64), with no 2^n object ever built. At n = 36 a state vector would hold 6.9·10¹⁰ amplitudes, 1.0 TiB.

Three API calls carry the whole pipeline: interop.stabilizers_from_stim(...) turns a stim tableau into signed generator strings, PauliSum.propagate(tail, policy, direction="heisenberg") evolves the observable, and PauliSum.expectation_stabilizer(generators) contracts it.

Results

The setup is an open 6×6 square lattice, 36 qubits. Preparation is the 2D cluster state — H on every qubit then one CZ per edge, 96 Clifford gates — a universal resource state that is not a product state in any local basis. stabilizers_from_stim reads the 36 signed generators out of the tableau in under a millisecond. The observables are two sites’ own cluster stabilizers, K_q = X_q ∏_{q'∈N(q)} Z_{q'}, at the lattice centre (degree 4) and a corner (degree 2); both read exactly +1 before the tail runs. The tail is two kicked-Ising Trotter steps (rx(θ_h) then ZZ(-π/2) per step), with θ_h swept from 0 to π/2:

(all 17 sweep points: theta_sweep.csv)

Preparation depth is free on the Heisenberg side: it never sees the state, so one evolved observable is contracted against the same cluster state prepared through four circuits of growing depth (96 to 50 236 Clifford gates, via provably-identity padding rounds). The generator strings come back byte-identical at every depth and the estimate is unchanged to all twelve printed digits — a 523× longer preparation costs 2.1 ms instead of 0.5 ms on the stim side and nothing on the propagation side, because it is literally the same evolved sum contracted four times (prep_depth.csv).

Scope of the check

Preparation-depth invariance holds for structured preparations. A generic stabilizer state’s group holds 2^n of the 4^n Pauli strings, so a random low-weight operator term lands in the group with probability 2^-n (1.5·10⁻¹¹ at n = 36) — an unstructured deep Clifford preparation drives the estimate to exactly zero, by membership, not by truncation. The pipeline is not restricted to graph states, but the observable and the state need group elements in common; for a structured resource state (cluster, surface-code, GHZ) that is the normal case, and such a state’s own generators remain estimable regardless of preparation depth.

The closed-form check

Every point in the sweep equals cos^deg(q)(θ_h) to 1.1e−16, derived rather than fit: pushing K_q through the two-step circuit U = ZZ₂X₂ZZ₁X₁, the first ZZ layer is Clifford at θ_zz = -π/2 and turns K_q into ±X_q (weight one) whenever deg(q) is even; X_q commutes with the X rotations; the second ZZ layer reverses the first; the second X layer splits each Z_{q'} factor into a cos/sin pair, and only the all-cos branch survives contraction because every other branch carries a Y and contracts to zero. Odd-degree sites pick up one extra cos factor and a sign and are excluded from the closed form, not from the pipeline.

Performance

Contraction cost was measured directly against the O(m·n²/64) bound, on a 1D-chain cluster state at each n, 20 000 random full-support terms, setup cost subtracted:

Local slope d log(per-term time) / d log n climbs from 1.22 (64→128) to 1.85 (512→1024) against the bound’s 2.0 — below n ≈ 256 a per-term fixed cost (pivot scan, term decode) dilutes the n²/64 row work, and the last decade sits within 8% of the bound. Per-word-operation cost falls from 6.2 ns to 2.1 ns and flattens near the cost of a load, an xor pair and a parity. Linearity in m at n = 256 holds flat to 7% over two decades, from 1 000 to 100 000 terms (scaling.csv).

The full run takes 116 s and 10.2 GiB peak RSS, single-threaded (RAYON_NUM_THREADS=1), the minimum of five repeats after warm-up. Every timing here is one campaign’s numbers on a shared workstation; term counts and expectation values reproduce to the digit, milliseconds do not — read the exponents and the flatness, not the last digit.

Validation

Three independent routes agree: a dense 2^n state-vector reference built straight off the Circuit object the engine propagates, a projector Π = ∏_i(I + s_i G_i)/2 built from the generator strings alone at n = 8, and qiskit Aer where installed. Worst gap over every case and route: 2.220e−15, against a 1e-10 bar; qiskit Aer’s worst gap is 1.55e-15.

A truncation-convergence panel (tail depths 4 and 5, sweeping min_abs_coeff) shows depth 4 self-converged inside the grid — the last two points agree to all twelve digits — while depth 5 is still moving by 5.2e-05 between its last two points, reported as the residual it is. The estimate converges faster than the sum: between 1e-4 and 1e-6 at depth 4 the term count grows 1.4× while the answer moves 2.3e-5, because only terms landing in the stabilizer group contribute at all.

Reproducing

source .venv/bin/activate
RAYON_NUM_THREADS=1 python examples/b7_stabilizer_prep/run_b7.py          # 116 s, 10.2 GiB peak RSS
RAYON_NUM_THREADS=1 python examples/b7_stabilizer_prep/run_b7.py --quick  #  40 s,  1.5 GiB peak RSS
pytest examples/tests/test_showcase_b7.py                        # 0.5 s

RAYON_NUM_THREADS=1 must be exported before the interpreter starts; the script asserts the pin and refuses to run unpinned. A full run regenerates every artifact in the showcase directory. --quick writes into the gitignored benchmarks/results/b7-quick/ instead, and --out-dir DIR picks another directory.

Numbers: every value on this page is sourced from examples/b7_stabilizer_prep/README.md and its committed artifacts — theta_sweep.csv, scaling.csv, prep_depth.csv, validation_b7.json, results_b7.json.

B8 — Operator Loschmidt echo

Source: examples/b8_ole/

The 56-qubit operator Loschmidt echo of arXiv:2607.25998 on ibm_boston, estimated by the paper's hybrid Pauli-propagation Monte Carlo (PP-MC), with caches that extend past one node through the MPI engine.

The circuit

The circuit is the Quantum Advantage Tracker’s operator_loschmidt_echo_56x1488: a CZ kicked-Ising Floquet echo U V_δ U† on a heavy-hex patch. Its 11 scattering qubits carry an rx angle shifted by η = 3πα/2. echo_half(L, η) rebuilds the measured half C for any depth and reproduces the six vendored tracker files gate for gate. The echo is S_δ = 2⁻ⁿ Tr(A V_δ† A V_δ), where A = C†OC, O = Z^{⊗12}, and V_δ = ⊗ e^{-iδX} on 35 qubits with δ = 0.3.

PP-MC

import ole
from paulistrings import truncation
policy = truncation.collapse_sample(cache=500_000_000, seed=s)
A = O.propagate(C, policy, direction="heisenberg")
S = ole.diagonal_echo(ole.anticommute_histogram(A, sites, "x"), 0.3)

• collapse_sample propagates exactly until a gate leaves more than cache strings.
• It then draws one string with probability b_P²/Σb² and restarts from it.
• The last block is kept whole.
• The diagonal estimator Σ b_P² cos(2δ)^{n_P} / Σ b² is the resummed moment series of the paper’s Eq. (C34).
• ole.rotated_overlap gives the exact echo of the same A, off-diagonal terms included.
• Both estimators are numpy over A.to_arrays(); under comm= each rank reads its result="local" share and comm.allreduce sums the values.
• collapse_sample draws the same string whatever the thread count.

Cost

One L = 6 trajectory at η = 9π/40 on 8 threads of ccqlin038 (Xeon Gold 6244):

For comparison, the paper’s Table VI reports about 8 minutes on 10 cores and 8–10 GB at cache 5e7. scripts/slurm/ole-ppmc-disbatch.sbatch runs the 800-trajectory sweep, and scripts/slurm/ole-mpi.sbatch runs caches beyond one node’s memory.

Benchmarks

Five benchmarks, each with a setup, an oracle, and a result. Two of the five results are negative.

Comparability rules

Four of them are worth stating up front, because they are what makes the tables comparable:

1. One gate per channel, everywhere. Truncation is applied after every channel, so fusing two gates into one channel changes the answer. Every circuit in the suite is built one gate per Circuit method call, which is also what makes a per-layer comparison against another engine meaningful.
2. Term-count parity blocks timing. No cross-engine wall time is reported for a configuration whose evolved Pauli sums diverge term-for-term at matched truncation. The parity gate runs first, untimed, and compares every per-layer count in application order, not just the final one — a divergence that cancels by the end is exactly the bug the check exists to catch.
3. Single-threaded, warm, with input generation outside the timed region. RAYON_NUM_THREADS=1 exported before the interpreter starts. References are exempt: an oracle is not a timing measurement, so reference sweeps are allowed threads (and are run in a spawned child, which also confines qiskit-aer’s persistent OpenMP pool).
4. Every truncated result ships with a convergence panel, and a self-converged reference may only be quoted if its plateau test passed and its reported uncertainty is inside half the accuracy bar. Rows failing either test are printed as not claimable and no value is quoted from them.

The plateau criterion

The obvious self-convergence test is “tighten the cutoff until two successive values agree to tol”. That test is wrong here, and Benchmark B caught it. Run against an exact reference at a small kick angle, it declared convergence with an estimated uncertainty of exactly zero while the value was still 5.6·10⁻⁷ from the truth, because at a small kick angle the only terms contributing to ⟨0|O|0⟩ are those rotated all the way to pure Z, so loosening the cutoff by a decade admits thousands of new terms none of which is pure Z, and the expectation does not move at all while the sum keeps growing. An exactly-zero difference there means “no relevant term has arrived yet”.

So the criterion in force requires the two small successive differences and one of: the term count has stopped growing (the sum has saturated, and the plateau is the exact answer), or both differences are strictly nonzero (an ordinary slowly-converging series). A flat value with a still-growing sum is rejected. A sum truncated to zero terms is rejected outright, however flat it looks. The fix is worth a measured 190× in accuracy.

Benchmarks C and B2 import that criterion as a function object rather than re-implementing it, and a test asserts it is the same object. Validating a result states it as user guidance, with the sweep that feeds it.

Reproducing a benchmark

./scripts/setup.sh && source .venv/bin/activate
uv sync --extra examples --extra bench   # also rebuilds the extension after a Rust change

RAYON_NUM_THREADS=1 pytest benchmarks/python/bench_a_clifford.py --benchmark-only
RAYON_NUM_THREADS=1 python benchmarks/python/bench_b_theta_sweep.py --validate-convergence
RAYON_NUM_THREADS=1 python benchmarks/python/bench_c_deep_trotter.py --validate-convergence
RAYON_NUM_THREADS=1 python examples/xxz_chain/run_benchmark_d.py all
RAYON_NUM_THREADS=1 python benchmarks/python/bench_e_su4.py

None of these is in CI. Each has a CI-safe correctness gate at smaller scale under benchmarks/python/tests/, so the physics is checked on every commit.

The cross-library construction/conjugation comparison against qiskit.SparsePauliOp and openfermion.QubitOperator lives in benchmarks/python/bench_baseline.py.

Reading a wall time

Term counts, expectation values, parity outcomes and convergence verdicts are load-independent — quote those. A wall-time delta under ~10% needs scripts/ab-compare.sh (two prebuilt binaries alternated adjacent in time, paired per-run deltas, acceptance by direction consistency across every pair), not these tables.

Sources: benchmarks/README.md and the per-benchmark READMEs linked from each page.

Benchmark A — Clifford point

Manual: Validating a result

The one point in the suite where the exact answer is an integer. At a Clifford kick angle a single Pauli string evolves into a single Pauli string with coefficient ±1, so stim can be asked for the truth and the engine has nowhere to hide.

Setup

127-qubit heavy-hex kicked Ising, 5 Trotter steps, θ_h = π/2 — the Clifford point of the utility experiment — timing the Heisenberg propagation of the published weight-10 and weight-17 observables against |0…0⟩. One pytest-benchmark group per observable. min_abs_coeff = 1e-8, single-threaded, seeded fixtures built outside the timed region.

The cutoff is 1e-8, not a dyadic value. This engine drops |c| <= eps while PauliPropagation.jl keeps |c| == eps — a genuine, measured cross-engine divergence. Clifford angles produce exact dyadic coefficients, so a dyadic cutoff is exactly where that boundary is likely to be hit bit-for-bit. 1e-8 is far from any dyadic value, so it changes nothing about which branch survives.

Oracle

stim, at the Clifford point, exact: the tableau simulation gives the ±1 integer directly. Every paulistrings entry also asserts the Clifford invariant on its own result — a single term with coefficient exactly ±1.0 — so a correctness regression fails the benchmark run too, not only the dedicated test file.

Result

The timing numbers for this benchmark live in benchmarks/results/bench_a.json, which is gitignored — the repository commits no timing table for A, so this page quotes none. What is committed, and what this benchmark actually establishes, is correctness and parity:

• The exact integers are reproduced bit-exactly. Benchmark B scores the same Clifford endpoints across a full eight-point coefficient sweep and finds a worst deviation of 0 for every observable at both endpoints: +1 for weight-10 and −1 for weight-17 at θ_h = π/2, matching A. See Benchmark B.
• Cross-engine parity is required before either engine’s timed entry runs: both are run once, untimed, and every one of the 1355 per-layer term counts must be identical. The schema-v1 task JSON handed to PauliPropagation.jl is built from the same recorded gate list the paulistrings side runs, so neither engine gets a transcription of the other’s circuit.

Regenerate the timings with:

RAYON_NUM_THREADS=1 pytest benchmarks/python/bench_a_clifford.py \
    --benchmark-only --benchmark-json=benchmarks/results/bench_a.json

The correctness gate

The CI-safe half lives in benchmarks/python/tests/test_benchmark_a_clifford.py and importorskips stim, so the numpy-only CI job stays green without it:

pytest benchmarks/python/tests/test_benchmark_a_clifford.py

Sources: benchmarks/README.md, the driver’s own module docstring in benchmarks/python/bench_a_clifford.py, and, for the reproduced integers, benchmarks/python/theta_sweep/README.md.

B — Kick-angle sweep

Manual: Validating a result, Truncation

Six kick angles, three published observables, 5 Trotter steps: how far a truncated Pauli sum drifts from the exact answer, and what an exact reference costs at each cone size.

The plot the benchmark exists for: |error| against warm wall time, over the truncation grid, for every observable and angle.

Setup

Heavy-hex kicked Ising, n = 127, 5 Trotter steps, θ_zz = −π/2, six kick angles θ_h ∈ {0, 0.2, π/8, π/4, 3π/8, π/2}, three observables from the utility experiment (Z_62, a weight-10 operator, a weight-17 operator). Heisenberg picture, contracted against |0…0⟩, single-threaded, warm timings, 1355 channels per circuit at one gate per channel. Recorded run: ccqlin038 (Xeon Gold 6244 @ 3.60 GHz), rustc 1.94.0, Python 3.11.11, numpy 2.4.6, qiskit 2.5.2, stim 1.16.0, julia 1.12.6 + PauliPropagation 0.8.2, 39.6 min end to end, 236 records.

Oracle

Reachable exactness is set by cone size, which the driver recomputes from the gate list on every run rather than assuming it:

Both endpoints (θ_h ∈ {0, π/2}) are Clifford for every observable and get exact integers from stim in under 0.1 s.

Weight-10 uses a statevector instead of Pauli propagation over the same cone on cost alone: Pauli terms over the 30-qubit cone hit 4.3e8 at 160 applied gates (87 s) and abort past a 26 GiB cap before the 305-gate cone completes, while the statevector is a flat ~150 s and 16.1 GiB regardless of depth. Z_62’s 19-qubit cone is small enough for the Pauli path to be affordable instead (12.8 s against 2.8 s), which is why it gets both methods behind one reference. Weight-17 has neither option: 59 qubits rules out any dense method, and untruncated propagation over its 59-qubit, 403-gate cone is far past the weight-10 wall above, hence the self-converged reference.

The self-convergence criterion

A criterion that tightens the cutoff until two successive expectation values agree to tol reports convergence too early here. Run against the exact Z_62 reference at θ_h = 0.2, it settles on an estimated uncertainty of exactly zero while the value is still 5.6·10⁻⁷ from the truth:

eps=1e-03  <O>=+0.981076800  Δ=  —        terms=     59
eps=1e-04  <O>=+0.981076800  Δ= 0.00e+00  terms=    225
...   bit-identical through eps=1e-07, then
eps=1e-08  <O>=+0.981077278  Δ= 4.78e-07  terms=  21595
exact                +0.981077357894

At a small kick angle only terms rotated all the way to pure Z contribute to ⟨0|O|0⟩, and each rotation costs a factor sin θ_h. Loosening the cutoff by a decade admits thousands of new terms, but for several decades none of them is pure Z, so the expectation does not move: a zero difference here means no relevant term has arrived yet, not convergence.

The criterion used in this benchmark also requires the term count to have stopped growing (the sum has saturated and the plateau is exact), or both successive differences to be strictly nonzero (ordinary slow convergence). A flat value on a still-growing sum is rejected, and a sum truncated to zero terms is rejected outright however flat it looks. With this criterion, Z_62 at θ_h = 0.2 self-converges to a true error of 2.98·10⁻⁹ with an estimate of 4.78·10⁻⁷ — a ~190× tighter estimate than the naive rule, which reports the same run’s 5.58·10⁻⁷ true error as zero.

Reference values

Z_62 and weight-10 carry exact references at every interior angle (12 decimal digits, in the committed summary.json). Weight-17’s four interior references (7.67e-08, 0, 3.62e-05, 1.57e-03 uncertainty at θ_h = 0.2, π/8, π/4, 3π/8) did not converge: the plateau test refused every one, and the budget guard stopped each sweep first. Reaching ⟨0|O|0⟩ from a weight-17 seed means rotating all seventeen X/Y sites to Z, so coefficients sit far below any affordable cutoff except near θ_h = π/2 — at π/8 the expectation is identically zero out to 1e-5 and 5.2·10⁷ terms, consistent with zero at the resolution reached rather than a measurement. An error scored against a self-converged reference is not an error against the truth: for weight-17 the reference is this engine’s own tightest run, so the tightest point of each sweep has an error of ~0 by construction.

Clifford endpoints

Scored over the coefficient sweep (8 cutoffs from 1e-2 to 1e-9) at each endpoint:

Every endpoint reproduces bit-exactly at every cutoff, matching Benchmark A. The weight-cap column is a genuine finding: at a Clifford angle the back-evolved operator passes through weight ~30–40 mid-circuit even though it lands on a single low-weight string, so max_weight <= 8 (<= 12 for weight-17) empties the sum; weight_10 recovers the exact +1 as soon as the cap reaches 10.

Time to |error| < 1e-3

z62 and weight-10 clear the bar at every angle in under 3 s, except weight-10 at π/4 and 3π/8, which never clear it inside this run’s grid (cut off at 1e-4: errors 2.00·10⁻³ and 2.29·10⁻³; the next decade would cost ~212 s and ~9.3·10⁷ terms). z62’s priciest passing point is π/4: 2.44 s at 2 146 372 terms. Three weight-17 rows pass vacuously — references of ~10⁻⁷–10⁻⁴ already sit inside 1e-3, so a 0.0 returned with zero terms “passes” without measuring anything; the fourth, at 3π/8, needs 101 s and 4 291 840 terms.

Non-monotone truncation error

weight-10 at the two largest interior angles, against exact references:

The loosest cutoff is closest at π/4: discarded terms carry signs, so a partial sum can sit nearer the truth than a larger one, and a truncated Pauli sum has no variational bound against it. The CI test asserts convergence across its grid rather than point to point.

Self-convergence validation

The weight-17 procedure, run at every z62/weight-10 angle against exact references and scored by whether its estimated uncertainty stayed honest (true_error ≤ max(10 × estimate, 1e-12)), is 8/8 conservative, the estimate over-stating the true error by 1.6×–1000× (full row breakdown in the committed summary.json). The two zero estimates are plateaus reached by saturation, the remaining error being summation rounding (1.3·10⁻¹⁴, 6.9·10⁻¹⁸), which sets the floating-point floor in the comparison. The two budget-stopped rows are still honest: a larger uncertainty, not a falsely confident one. That property does not survive to Benchmark C’s depth — see C.

Cross-engine parity

Matched truncation, one gate per channel on both sides, Heisenberg, |0…0⟩, single-threaded, θ_h = 0.2, min_abs_coeff ∈ {1e-3, 1e-4, 1e-5} — three cutoffs rather than one, compared per applied layer, all 1355 of them.

Per-layer counts are 1355/1355 identical for all nine (observable, cutoff) pairs. Final terms and |Δ⟨O⟩|:

9/9 pass: all 12 195 compared per-layer term counts are identical, and eight of the nine expectation values agree to the last bit. Cutoffs are powers of ten, deliberately not dyadic, and no eps perturbation was needed. See Comparisons for why that choice matters.

Two side findings. On the heaviest case (2.85·10⁶ resident terms) PauliPropagation.jl’s dict backend peaked at 67.6 GiB RSS against ~1.2 GiB for this engine’s bucketed columns, not a clean like-for-like (part of the jl figure is a separate untimed counting pass), but worth knowing before anyone reruns it on a shared workstation. It also does not reproduce: a per-process sampler in Benchmark C measured the same quantity 30–50× lower. Wall time is not compared here — the jl timings are a single warm repeat on a loaded machine.

Reproducing

source .venv/bin/activate
RAYON_NUM_THREADS=1 python benchmarks/python/bench_b_theta_sweep.py --validate-convergence
pytest benchmarks/python/tests/test_benchmark_b_sweep.py    # the CI gate, on a 20-qubit sublattice

RAYON_NUM_THREADS=1 must be exported before the interpreter starts; the driver refuses to run otherwise.

Caveats

• Reference computations run in a spawned child process, since qiskit-aer’s statevector simulator spawns an OpenMP pool that would trip the single-thread assertion later; only plain data crosses back. The weight-17 self-convergence child additionally gets 16 Rayon workers, since a reference is an oracle rather than a timing measurement; the workers buy the cutoff reach that reached 1.0·10⁸ terms at θ_h = 0.2.
• Timings were taken on a shared workstation with other work running: the same weight-10 configuration measured 13.8 s and 28.6 s minutes apart. Term counts, expectation values and parity outcomes are load-independent.
• min_abs_coeff ≥ 1e-12 everywhere.
• Two derived fields in the committed summary.json were recomputed from the recorded measurements; that file’s notes field says so.

Numbers: every figure on this page traces to benchmarks/python/theta_sweep/README.md, with the 236 raw records in results.json and the verdicts in summary.json next to it.

Benchmark C — Deep Trotter circuits

Manual: Validating a result

At 5 Trotter steps the 0.01 accuracy target is met in a tenth of a second against an exact reference. At 15–20 steps in the hard interior neither the target nor a reference to score it against is reachable, and the published record agrees: the paper that set the bar publishes no exact 20-step value either. That reachability boundary, not a single accuracy number, is this benchmark's headline.

|error| against warm wall time, one curve per depth in the ladder, with the 0.01 accuracy target drawn. Read the claimable column below, not the curve’s crossing point.

Setup

Heavy-hex kicked Ising, n = 127, θ_zz = −π/2, observable Z_62 (the weight-1 operator whose 20-step point is the marquee number of the utility experiment), a depth ladder of 5 / 9 / 15 / 20 Trotter steps, two kick angles in the hard interior θ_h ∈ {7π/32, 5π/16}, and a dyadic truncation grid min_abs_coeff ∈ {2⁻¹⁴, 2⁻¹⁶, 2⁻¹⁸} extended with three looser dyadics. Heisenberg, contracted against |0…0⟩, single-threaded, warm timings, one gate per channel: 5 420 channels at 20 steps.

Recorded run: ccqlin038 (Intel Xeon Gold 6244 @ 3.60 GHz), rustc 1.94.0, Python 3.11.11, julia 1.12.6 + PauliPropagation 0.8.2. 85.7 min end to end, 42 records.

Angle, depth, grid and reference

7π/32 = 0.687223… and 5π/16 = 0.981748… sit on the k·π/32 grid of the published exact benchmarks (Begušić, Gray & Chan, arXiv:2308.05077), which also fixes the rung at 9 steps rather than 10: that data’s exact.csv has a 5a column (⟨Z₆₂⟩ after 9 steps) and no 5b column (after 20), so the paper that set the 0.01-accuracy bar publishes no exact 20-step value either — independent corroboration of the finding below. No reference file is checked in, so every number here is computed by an oracle in this repository. The full six-dyadic grid runs at 5 steps; at 9, 15 and 20 steps it stops at 2⁻¹⁴, since past that point at 20 steps the term count projects to ~6.2·10⁸, over the 4·10⁸ guard. The two angles bracket the behaviour: 7π/32 keeps an O(0.5) signal at 20 steps, 5π/16 decays to O(0.01).

The reference is the commutation-aware backward cone of Z_62: exact at 5 steps (19 q, 83 gates), where light_cone_exact(method="both") requires an Aer statevector and an untruncated Pauli propagation to agree, and self-converged beyond (65 q at 9 steps, 127 q at 15 and 20). Past 5 steps no dense method or published exact value exists, so the tolerance is 1e-3 (10× headroom under the 0.01 bar) and the grid extends two dyadic powers past the tightest timed point, measuring every timed run against something strictly tighter than itself — the machinery is Benchmark B’s.

A reference is claimable if it is exact, or self-converged with the plateau test converged and the uncertainty inside 0.01 / 2 — otherwise an “achieved” row is circular, since the tightest timed point agrees with its own tightest run by construction. Only claimable rows are quoted below.

Results

Reference values, with convergence evidence

7π/32 at 20 steps is the clearest negative result here: the reference sweep reached 2⁻¹⁶ (3.9·10⁷ terms, 276 s at 16 threads) and never settled.

The partial sums swing between 0.25 and 0.52 as the cutoff tightens by factors of four — the signature of discarded terms that carry signs, so a truncated Pauli sum has no variational bound. The grid is nowhere near resolving this point: uncertainty 1.44·10⁻¹, 14× the target.

5π/16 at 20 steps is the one deep point that converges (0.011941, 0.014730, 0.015481, 0.016131 at 2⁻⁸ … 2⁻¹⁴, two successive differences below 1e-3), because the observable has decayed to ⟨Z₆₂⟩ ≈ 0.016 and the sum collapses from a 1.4·10⁷-term transient to 34 698 resident terms.

Non-bound uncertainty estimate

Its bias flips sign. Measured on the 20-qubit sublattice at 20 steps, where a dense statevector gives the truth, over the full grid it is conservative; stopped after two points, what a budget guard does, it is not:

The last row is past the 10× slack Benchmark B’s validation allows: a budget-truncated sweep no longer just reports a weaker uncertainty here, it can report a falsely confident one. converged still stays false, since two successive small differences were never seen, and at this depth that must be read as no usable estimate, not a weaker one.

Time to |error| < 0.01

Read the claimable column, not reached?: every row reaches the bar, and most of those passes are meaningless.

• The two exact-reference rows are the real measurements: the 0.01 target is met in 0.11 s / 5.9·10⁴ terms and 0.14 s / 7.2·10⁴ terms against a reference with no truncation anywhere.
• 7π/32 at 15 and 20 steps “pass” at 2⁻¹⁰ by luck. The reference is itself unresolved and the error against it is not monotone in the cutoff (at 20 steps: 1.24e-1, 1.28e-3, 1.51e-1, 1.44e-1 through 2⁻⁸ … 2⁻¹⁴); the 2⁻¹⁰ pass is one partial sum crossing the wrong reference, exactly what the claimability test exists to catch.
• 5π/16 at 15 and 20 steps pass vacuously: the signal is 0.041 and 0.016, so a 0.01 absolute bar is 25% and 62% of the whole answer, and 11 resident terms clear it at 20 steps.

The 0.01 target is met, against an exact reference, in a tenth of a second at 5 steps; at 9 steps it needs 2⁻¹⁴, 2.2·10⁶ terms and 37 s against a reference resolved only to 1e-2; at 15–20 steps neither the target nor a reference to score it against is reachable inside this benchmark’s box.

Reachability requirements at 20 steps, θ_h = 7π/32

Extrapolating the measured growth (15.9× in terms and ~15× in wall time per factor of four in the cutoff, against a measured error law of only ~2.2× error reduction per the same step):

The 0.01 bar needs roughly 2⁻²⁰–2⁻²² here, out of reach of a workstation. Showcase B2 revisits this point with noise: 48× fewer peak terms, a 154× smaller last difference, and still not a pass.

Sanity envelope and methodology check

12 records fall on the three cutoffs tracked here, scored against a 1.2·10⁶–9.3·10⁶ tracked-set envelope: nine of twelve are inside or expectedly below, and the three “above ceiling” readings are the same 1.44·10⁷-term transient at 5π/16 and 2⁻¹⁴. The envelope lands squarely on the 20-step, 2⁻¹⁴ point it was quoted for, 3 108 582 peak terms, dead centre: the accuracy shortfall above is what that tracked-set size actually buys at this depth, not a mis-specified circuit.

Running the self-convergence procedure at the real system size, at the two depths where an exact reference exists, both converge by saturation to the exact answer at floating-point precision: true error 3.44e-15 against an estimated 3.71e-05 at 7π/32, and 2.22e-16 against 0 at 5π/16 — both conservative.

The dyadic cutoffs and the one-ulp mitigation

Exact dyadic cutoffs are the one case where this engine and PauliPropagation.jl provably disagree: this repository drops |c| <= eps, jl keeps |c| == eps, and at a Clifford θ_zz the coefficients are exact dyadics too, so an exact straddle is not measure-zero, unlike Benchmark B’s powers of ten. The mitigation: paulistrings runs use the dyadic verbatim, and jl runs get math.nextafter(eps, inf) — since no float sits strictly between eps and eps′, |c| < eps′ is exactly |c| <= eps, jl’s rule becomes this engine’s rule bit for bit, with no coefficient touched. The perturbation is one ulp: 1.1·10⁻¹⁹ absolute at 2⁻¹⁴.

Cross-engine parity at the deepest point

Matched truncation at θ_h = 7π/32, 20 Trotter steps, the deepest point in the benchmark, one gate per channel on both sides, Heisenberg, |0…0⟩, single-threaded, at the three dyadic cutoffs the memory gate allowed, per applied layer, all 5 420 of them:

3/3 pass: every one of the 16 260 compared per-layer term counts is identical, final and peak counts agree exactly, and the expectations agree to ≤ 5.6·10⁻¹⁷ against a 1e-9 bar. This is the first place in the suite where the one-ulp mitigation is load-bearing.

Memory

A figure from Benchmark B does not reproduce here: extrapolating this section’s fit to B’s case gives ~2–3 GiB, not the 67.6 GiB B reported. Every jl leg ran, none skipped for memory; the gate’s affine model, refitted from each leg’s directly-sampled RSS, converges on:

So jl’s dict backend costs ~0.44–0.74 KiB per resident term on this host plus a ~0.7 GiB fixed footprint, ~30–50× lower than the 24 KiB/term implied by B’s figure. The two measurements are not directly comparable, but the gap is worth naming: this figure samples /proc/<pid>/status VmRSS on the runner.jl process directly, twice a second, while getrusage(RUSAGE_CHILDREN).ru_maxrss is a process-lifetime maximum over all reaped children, dominated here by multi-gigabyte reference children — the same 1 925-term jl task read 3.68 GiB by getrusage and 0.66 GiB by direct sampling. Re-measure with a per-process sampler before quoting B’s figure again. For scale, this engine’s 2⁻¹⁴ sum is ~0.15 GiB by construction (3.1·10⁶ terms × 48 B at W=2), and the whole 42-run campaign’s process high-water was 1.11 GiB.

Wall time — reported, not claimed

jl is 1.4×, 2.3×, 2.3× slower here, a single warm repeat per point on a shared workstation, not a benchmark claim — but the gap is well outside the ±5–8% noise band and consistent in direction across two orders of magnitude in problem size. The numbers to quote from this benchmark are the term counts and accuracy rows, which are load-independent; Benchmark D characterises the crossover properly.

Reproducing

source .venv/bin/activate
RAYON_NUM_THREADS=1 python benchmarks/python/bench_c_deep_trotter.py --validate-convergence
pytest benchmarks/python/tests/test_benchmark_c_deep.py    # the CI gate: 25 tests, ~50 s

Caveats

• References run in a spawned child with 16 Rayon workers (measured 11.2× cutoff reach at 32 threads on the 20-step, 2⁻¹⁴ point), since a reference is an oracle, not a timing measurement; the child also confines qiskit-aer’s persistent OpenMP pool.
• Timings were taken on a shared workstation. The recorded run started on a quiet box (load 12, 238 GiB free) and stayed there.
• min_abs_coeff ≥ 1e-12 everywhere; the tightest cutoff used here is 2⁻²² ≈ 2.4·10⁻⁷, far above the floor.
• Peak vs final term count. At 5π/16 and 20 steps the sum peaks at ~1.5·10⁷ resident terms and lands on ~2·10⁴; every term-count figure here uses the peak, since that is what a run has to hold.

Numbers: every value on this page is computed by benchmarks/python/deep_trotter/bench_c_deep_trotter.py, with the raw records in results.json and the full convergence evidence in summary.json next to it.

D — XXZ chain scaling

Manual: Truncation, Stats, memory and logging

The only benchmark in the suite with an analytic prediction to check against, and the one that finds the cross-engine ranking changes sign as the tracked set grows. Below roughly 10⁴ tracked terms PauliPropagation.jl is 3–4× faster; above it this engine is ~1.5× faster and pulling away. A single-point comparison would have "shown" either engine winning by 3–4×.

Setup

The open XXZ chain

H = Σ_{i=0}^{n-2} ( X_i X_{i+1} + Y_i Y_{i+1} + Jz Z_i Z_{i+1} )

is first-order-Trotterized at dt = 0.1: three pauli_rotation channels per bond, even bonds then odd bonds, one gate per channel. Two regimes: Jz = 0 (free) and Jz = 0.5 (interacting). Observables: the central Z_c (weight 1) and Z_c Z_{c+1} (weight 2). Direction is always Heisenberg, initial state always a domain wall |0…01…1⟩.

|0…0⟩ is an eigenstate of H at every Jz, so every expectation would be a constant; |+…+⟩ gives ⟨Z_c⟩ = 0 by symmetry. The domain wall is a computational basis state (so PauliPropagation.jl can contract it too), and ⟨Z_c⟩ starts at exactly −1 and moves. D’s deliverable is a scaling sweep (36 scaling points plus two truncation grids), which is why it lives in examples/xxz_chain/ rather than with the other benchmarks (no pytest-benchmark entry; see Limitations).

Quadratic growth law at Jz = 0

The measured log-log slope of untruncated non-zero term count against Trotter steps is exactly 2.0000 at n = 40, 60, 80, 100, weight-1 seed. The counts are 16 s² exactly for every unsaturated point, independent of dt (0.05, 0.1, 0.37 all identical) and of the seed site:

The n = 40 row is deliberate: past s = 10 its light cone reaches both chain ends and the count saturates (the law’s boundary of validity), with a fit window of s ≤ n/4 − 1; a test asserts the excluded tail breaks the slope. At Jz = 0 the XX+YY chain is free-fermion: Jordan–Wigner maps it to hopping Majoranas, and every Trotter gate is Gaussian, its conjugation acting on the Majorana operators linearly with O orthogonal. A single-site Z_c = −i γ_{2c} γ_{2c+1} is a Majorana bilinear, so it stays a sum of bilinears, and each γ_a γ_b is exactly one Pauli string. The reachable (a, b) pairs number the square of the reachable Majorana indices, hence O(s²), and with this Trotter decomposition’s cone exactly (4s)² = 16 s². The CI test pins the slope, not the counts: the prefactor 16 is a property of the bond-sweep ordering, not of the engine.

Z_c Z_{c+1} is a Majorana quartic, so the same argument gives O(s⁴); measured counts are exactly (8s² + 6s − 1)² (169, 1849, 7921, 22801, 52441, 104329, 187489, 312481), whose fitted log-log slope over 2 ≤ s ≤ 8 is 3.70, the asymptotic 4 approached from below. Jz = 0.5 for contrast (same seed, n = 40, untruncated): 40 → 9512 → 2 453 872 terms at s = 1, 2, 3, then exponential growth as the interaction breaks the bilinear closure — the regime truncation exists for, and why D’s large-n results are truncated and convergence-checked rather than exact.

Oracle: statevector at n ≤ 26

qiskit Aer dense statevector against the engine at min_abs_coeff = 1e-12, for both regimes, both observables, n ∈ {20, 24, 26}, up to three Trotter steps — 30 cases. All 30 cases agree, worst |Δ| = 2.2·10⁻¹² against a 10⁻⁹ bar (median 2.0·10⁻¹⁵). The two worst cases are the deepest interacting ones (~818 000 terms after the 10⁻¹² cutoff), so the error is the cutoff’s, not a convention’s; the free-regime shallow cases are at 0 or 2·10⁻¹⁶. The two omitted case families (n = 26 with s = 3, and Jz = 0.5 with the weight-2 seed at s = 3) are cost, not doubt: both keep ~4·10⁷ untruncated terms.

A measured gotcha, worth knowing before adding any oracle-using timed run: one statevector_expectation call takes the process from 32 to 97 threads (Aer’s own OpenMP pool), and the harness’s thread pin counts threads gained since import, so any threads=1 run after an Aer call fails the assertion — both oracle-using modes therefore run every propagation first.

Time and peak memory against n

6 Trotter steps, matched min_abs_coeff = 1e-6 in both regimes, n = 20…100. One Python subprocess per point, because VmHWM is a process-lifetime high-water mark: a single process would report every later point’s memory as the largest earlier one’s. Warm propagation time (s), one discarded warmup per point:

Term count does not grow with n at all: ~257 (Jz=0, Z_c), ~7841 (Jz=0, weight 2), ~206 000 (Jz=0.5, Z_c), 517 000–822 000 (Jz=0.5, weight 2), essentially flat from n = 20 to n = 100. At fixed depth the operator’s light cone, not the chain, sets the size of the tracked set: six Trotter steps reach at most ~25 sites, so every n ≥ ~30 run is the same physics padded with identity. Time is nevertheless linear in n (Jz=0.5, Z_c: 0.40 s → 1.47 s for a 5.2× channel count): the circuit has 3(n−1) channels per step and the engine pays a pass over the sum for each one regardless of support, so cost is channels × terms. A support-aware channel skip would flatten this curve; nothing in the current engine has one. The counts alternate along the sweep (7841 vs 6631; 821 750 vs 517 398): the seed sits on bond c = n//2, an even bond when n ≡ 0 (mod 4) and odd otherwise, so the even-then-odd sweep hits it in a different half-step each time, changing the truncation schedule (a property of the bond ordering, not noise). Peak memory tracks the term count and the compile-time width tier, and nothing else: for the Jz=0.5, Z_c series (~206 000 terms at every n), growth is 11.2 MiB for n ≤ 60 and 17.5 MiB for n ≥ 70, 55 B/term against 87 B/term — the W = 1 → W = 2 boundary at 64 qubits (64-bit vs 128-bit symplectic keys: 32 B/term of key becomes 48 B/term).

Convergence panels

n = 24, 6 Trotter steps, Z_c, against the exact statevector value:

Exact reference: ⟨Z_c⟩ = +0.129334949228 (Jz=0), +0.024779796796 (Jz=0.5). The free regime reaches 10⁻⁸ with 365 terms and 4 ms; the interacting regime needs 2.7·10⁶ terms and 5 s to reach 3·10⁻⁷. The error is not monotone in the cutoff (Jz=0: 3.3e−6 at 1e−4, then 9.0e−6 at 1e−5; Jz=0.5: 2.9e−7 at 1e−7, then 4.3e−7 at 1e−8), since dropped terms carry signs and can cancel, so a tighter cutoff can land slightly further from the exact value while the trend still converges. The last decade buys nothing measurable: below ~3·10⁻⁷ the comparison stops resolving, since 2.7·10⁶ floating-point coefficients summed in an unspecified order, contracted against a dense reference with an error budget of its own, is the floor of the comparison, not of the truncation.

The self-converged panels at n = 60 and n = 100 produce bit-identical values at every cutoff (+0.0247793658 at 10⁻⁸ for both) and identical term counts, differing only in wall time (9.6 s vs 16.5 s at 10⁻⁸, the channel-count effect above): a 60-site chain and a 100-site chain are the same calculation at six Trotter steps, so the n = 100 self-converged answer is corroborated by an n = 24 run with an exact reference, where the same sweep converges to 4·10⁻⁷ of the statevector value.

Cross-engine timing

Blocking parity gate first, then warm times: per-layer term counts are compared index by index in gate-application order, and the task JSON is built from the same recorded gate list the engine runs, so neither side gets a transcription of the other’s circuit. Cutoffs are non-dyadic and strictly positive (1e-5, 1e-6), which avoids both measured divergences; no eps perturbation was needed. Parity: 4/4 cases pass, on 171–702 layers each — every per-layer term count identical index by index, final counts identical, expectations agreeing to ≤ 1.7·10⁻¹⁶ against a 10⁻¹² bar. Only then were the times below recorded (PauliPropagation.jl 0.8.2, julia 1.12.6, dict backend, -t1, warm minimum of 3 repeats; this engine warm after one discarded run).

(“ratio” = jl warm / paulistrings; > 1 means this engine is faster.) The ranking changes sign with the size of the tracked set, somewhere between 3·10³ and 3·10⁴ terms: below the crossover PauliPropagation.jl’s hash-map backend is 3–4× faster; above it this engine is ~1.5× faster and pulling away. These circuits have 3(n−1) channels per Trotter step (702 channels for the last row), and this engine pays a bucketed per-layer pass per channel: a fixed cost a 257-term sum cannot amortize. Both directions are far outside the ±5–8% single-threaded noise floor, so no A/B protocol is needed here: the first three rows were measured twice and reproduced their ratios as 0.24/0.28, 0.31/0.32 and 1.45/1.48, same ranking and order of magnitude, well under the sign changes being reported.

Reproducing

source .venv/bin/activate
RAYON_NUM_THREADS=1 python examples/xxz_chain/run_benchmark_d.py all
# or one mode at a time: growth | statevector | scaling | convergence | julia | figures
pytest benchmarks/python/tests/test_benchmark_d_xxz.py    # the CI gate: 11 tests, ~4 s

Results are committed as results/*.json, one file per mode, overwritten on rerun.

Limitations

No TDVP / tensor-network baseline exists at large n; no such package is a dependency of this repository, so the interacting large-n numbers are self-converged and labelled as a follow-up rather than silently approximated. Nor is there a pytest-benchmark entry: its per-point recalibration would rerun each point for no gain, and the memory curve needs one process per point.

Numbers: every value on this page traces to examples/xxz_chain/README.md and the raw sweeps in examples/xxz_chain/results/. Host: ccqlin038, Intel Xeon Gold 6244 @ 3.60 GHz, RAYON_NUM_THREADS=1, RUST_LOG unset.

E — Random SU(4) brickwork

Manual: Truncation

The deliberate worst case. Every other benchmark here is built from a physically structured circuit with commuting layers, Clifford points, or both, so truncation has real structure to exploit. This one draws an independent Haar-random SU(4) block for every brickwork site: no stabilizer structure, no commuting sublattice, no light-cone shortcut past nearest-neighbour causality. Every number below shows how the engine behaves with no help from problem structure.

Setup

n = 36 qubits, brickwork layers alternating on (0,1),(2,3),… and (1,2),(3,4),…, an independent Haar-random SU(4) block per site, observable Z_18 (the central qubit), Heisenberg direction. Seed 20260831, fixed throughout. Single-threaded.

Results

There is no Clifford savings anywhere in this circuit. Two regimes are visible in the numbers:

• Pre-saturation growth. At min_abs_coeff = 1e-4, peak term count goes 15 → 253 → 4038 → 58 079 → 667 088 → 2 296 294 → 4 110 454 for depth 1…7 — roughly an order of magnitude every one to two layers, faster than 4^depth: a generic 2-qubit Haar block’s adjoint action is a dense ~15×15 real orthogonal map on the traceless Pauli algebra, not a 4-way branch.
• Post-saturation collapse. At looser cutoffs the peak plateaus quickly (~1300–1700 terms at 1e-2, ~80 000–110 000 at 1e-3) and the final term count then falls to zero: by depth 10 (1e-2) or depth 20 (1e-3), no coefficient in the operator’s support survives the cutoff. This is the expected fate of a local operator diffusing through a non-integrable, non-Clifford brickwork — ⟨Z_18(t)⟩ reaches the maximally mixed value at truncation resolution.

The 1e-4 grid stops at depth 7 (25.0 s; depth 6 was 7.9 s) — per-step cost grows ~2.3–3× at this cutoff, so depth 8 alone would cost roughly a minute. The two looser truncations already show the full rise–plateau–collapse shape cheaply.

Error against runtime

n = 16, depth 6 (statevector-checkable, deep enough for fanout to be well underway). Oracle: −0.0497243601.

Error falls roughly a decade per decade of cutoff until ~1e-6, where term count and runtime saturate: two more decades of tightness add under 4% more terms and no measurable runtime, but still buy three more decades of accuracy, since the residual error there is dominated by which low-weight terms survive rather than raw term count.

The eps = 1e-2 row is non-monotone: its error (8.1e-2) is larger than the empty-sum row above it (4.97e-2, just |oracle|). Truncation error is not guaranteed monotone in the cutoff for a generic, non-sign-coherent operator.

Time and memory against n

Depth 6, min_abs_coeff = 1e-4, n from 8 to the 36-qubit headline. Oracle-checked up to n = 24; n ∈ {28, 32, 36} is self-converged only (oracle_checked: false in results.json). Wall time tracks term count, not n directly: this circuit’s causal structure means the light cone from the central qubit sets how many terms accumulate by a fixed depth, not n itself.

Every point’s peak_memory_kb_delta reads 0.0 in results.json. RSS is cumulative across a process; the figure is an upper bound when other suites ran first — a real per-point number needs one process per point, as in Benchmark D.

Cross-engine comparison

Schema v1’s unitary_2q gate is exactly the SU(4) block this circuit is built from, and PauliPropagation.jl 0.8.2 defines no Schrödinger transfer map for a TransferMapGate, so only direction="heisenberg" is comparable — the same restriction applies to any unitary_1q/unitary_2q/noise gate against this version.

min_abs_coeff = 1e-4 is deliberately non-dyadic: the one known coefficient-boundary divergence needs exact dyadic coefficients landing on the cutoff, which requires a Clifford point. Every coefficient in a Haar-random SU(4) circuit is an irrational float, so the boundary case cannot arise here.

Per-layer term-count parity holds exactly — matched at both sizes, including the final expectation to 1e-12:

Timed comparison (n = 10, depth 5, single-threaded both sides, warm/JIT-excluded on the jl side):

The two engines are within measurement noise of each other (±5–8% single-threaded is expected between any two runs on this host) — a parity-gated same-ballpark result, not a speed claim in either direction.

Validation

n ∈ {8, 12, 16, 20, 24}, depth 3, no truncation — the light cone at depth 3 has at most 15 non-identity Pauli letters, giving 4095 terms at every size. Agreement is to double-precision tolerance against qiskit Aer’s dense statevector, two engines sharing no simulation code:

Same seed twice gives identical term counts and expectation to 1e-12 (test_benchmark_e_su4.py::test_same_seed_twice_gives_identical_term_counts_and_expectation, always in CI), with a negative control confirming distinct seeds generically disagree.

Numbers: every figure on this page traces to benchmarks/python/su4_staircase/README.md, including the run command, with the raw records in results.json next to it.

Engine performance

Where the engine’s measured limits are: a roofline analysis of the current engine against the reference host’s measured memory bandwidth, then the partitioned and distributed measurements. Every number is copied from a committed source, named inline and in the footer.

Manual: Stats, memory and logging, NUMA partitions, MPI ranks

The design these numbers measure is Inside the engine; the thread-count advice they support is in Stats, memory and logging.

Layer time breakdown

Measured on the reference host (2× Xeon Gold 6244, 16c/32t) with the repo’s phase-timing probe, per layer class, as a share of summed worker busy time (gather / sort / merge):

Gather + merge dominate the sparse classes; the sort dominates only dense two-qubit PTMs. Parallel efficiency (busy time over coset-loop wall × threads) is 0.99 for su4 at 16 threads: load balance is a solved problem, and the scaling limits below are memory-system effects, not imbalance.

Per-term cost is flat in the sum size over the measured range: 30.5–30.6 ns/term for the ZZ rotation across m = 1.50e6 → 4.50e6, 138.5–141.3 for the sparse 2q unitary, 322–327 for su4 (single thread).

Memory wall

The ceiling comes from a STREAM-style probe (crates/membench via scripts/bandwidth.sh), nominal bytes with plain write-allocating stores — the same store pattern the engine uses. Best GB/s on the reference host:

Two structural facts from the ceiling measurement: hyperthreads add no bandwidth (39.0 → 39.2 GB/s), and the second socket adds only 15–25% under first-touch page placement with work-stealing, not 2×.

The modeled traffic side prices each layer from its row counts at T = 48 B/term (benchmarks/PROFILING.md §Roofline model):

bytes/layer = terms_in×T + 2×(rows_gathered−rows_id)×T + 2×rows_id×16
              + 2×rows_sorted×T + terms_out×T

Reading the ratio of measured DRAM traffic to ceiling as a classification: at or above ~70% of the measured ceiling a phase is bandwidth-bound; modeled traffic far above measured means the working set is cache-served; far below ceiling with a high LLC miss rate points at latency, not bandwidth.

Single-thread roofline

Measured on the current engine (fact sheet: research/HARDWARE.md §Engine roofline, single thread; ceilings per the 1-core row above):

The byte model over-counts DRAM traffic by 2.5–12.8× at one thread because most modeled traffic is served from cache. The sparse classes are limited by load latency (IPC ≈ 2.2–2.6 with 34–41% LLC load-miss); su4 is compute-bound in its sort (IPC 2.98, 2.1% LLC load-miss).

Multi-thread roofline

su4 at m = 1.41e7 (both-socket ceilings: read 45.0 / write 25.3 at 8–16t, 48.8 / 23.1 at 32t):

By 16 threads the dense-PTM class sits at the machine’s write ceiling — the figures above 100% are against the STREAM nominal-store ceiling, which the engine’s write-allocating stores share, so “at the ceiling” is the correct reading. 32 threads buy zero additional bandwidth (67.0 vs 67.0 GB/s attributable) and cost 12% of wall time.

The sparse classes at 32 threads stay latency-bound and keep scaling:

| cell | ns/term | speedup vs 1t | read GB/s | write GB/s | % write ceil | verdict | |—|—:|—:|—:|—:|—:|—:|—| | ZZ rotation, m = 4.50e6 | 2.7 | 11.3× | 8.9 | 11.5 | 50% | latency-bound | | sparse 2q, m = 3.0e6 | 10.8 | 13.1× | 11.7 | 11.3 | 49% | latency-bound |

They stop at half the write ceiling with 70–90% of their modeled traffic cache-served: bandwidth is not what limits them, and they take the full thread count profitably.

Partitioned: P = 2 vs P = 1

Two facts above — the second socket adding 15–25% rather than 2×, and the dense-PTM class pinned to the write ceiling — are the same effect: pages are placed by first touch and the thread pool then steals work across sockets, so roughly half of the second socket’s reads are remote. The engine can instead be run partitioned, one pinned pool and one share of the sum per NUMA domain, with only the rows a layer moves across a domain boundary exchanged (NUMA partitions).

Locality decides everything. A layer whose gates leave the partition rows fixed runs entirely inside its own domain and gains; a layer that moves any row across a boundary pays for the export pass, the copy and the receiver’s larger merge stream. Measured at P = 2 against P = 1 on the same binary, on two-socket hosts from 8 to 48 cores per socket:

The gain grows with cores per socket, as a memory-system effect should: −4% on an 8-core-per-socket workstation, −18% on a 48-core Genoa at 96 threads. Load imbalance across partitions stays within 1.000–1.004 with random rows. Random partition rows put roughly half of a dense two-qubit gate’s deltas across a boundary, so the default draw lands a mixed circuit in the bottom row of that table.

The per-host cells behind the summary, copied from research/HARDWARE.md §Partitioned engine, P=1 vs P=2 in-process: scripts/ab-compare.sh runtime-knob A/B, 5 pairs abba, --n 1000000 --qubits 128 (W = 2), random partition rows, each partition pinned to its socket with MPOL_BIND; steady-state m ≈ 1.5e6 for rotations, 1.0e6 for gu2q, 1.4e7 for su4/su4_local. Wall time per layer, P=2 vs P=1 (median Δ%, pairs agreeing):

Exchange-free rotation layer, per layer (from the probe sidecars):

Counters, ccqlin038, dense layers at 16 threads (8 per socket at P=2); scripts/perf-stat.sh, shared box:

partition_imbalance is 1.000–1.004 everywhere; exchange volume is bytes ≈ rows × 48.

Benchmarking

P = 1 versus P = N is a runtime-knob A/B of one binary, not a code A/B:

scripts/ab-compare.sh partitions-1v2 --a . --b . \
  --probe   '--n 1000000 --threads 32 --layers su4_local --partitions 1' \
  --probe-b '--n 1000000 --threads 32 --layers su4_local --partitions 2 \
             --partition-cpus "0-15,32-47;16-31,48-63"'

--threads is the total at every P. A partitioned cell runs under no placement prefix: numactl and taskset both defeat the split the engine is making. On an exclusive cluster node, scripts/slurm/ab-campaign.sbatch builds the CPU lists from the node’s own topology. The protocol and the sidecar fields are in benchmarks/PROFILING.md.

Distributed: weak scaling

The same split across processes is one partition per MPI rank (MPI ranks), and there the goal is capacity rather than bandwidth. Measured on Icelake ccq nodes (2 × 32 cores), one rank per NUMA domain, 32 threads per rank, UCX shared memory within a node and InfiniBand between nodes, at 6·10⁶ terms per rank:

Weak scaling is flat once the exchange leaves the node: the remote layer costs the same at 4 and 8 ranks, and local layers cost ~10.5 ms whatever the rank count. The step from 2 to 4 is shared memory giving way to InfiniBand. A remote layer is transfer-bound — the pipeline runs the coset loop under the transfer, so what is left is the export pass plus the bytes on the wire — which means the lever is fewer bytes (locality rows, a lower truncation), not more threads.

The full per-layer breakdown, copied from research/HARDWARE.md §Partitioned engine, MPI weak scaling: phase_breakdown --mpi, replicated input --n 4e6 × ranks → 6.0e6 terms per rank (4.0e6 for cnot), --reps 8; medians over ranks, ms per layer.

The remote layer is transfer-bound: ~6–8 ms compute, 7–8 ms export, ~30 ms transfer for 192 MB per rank (two ranks share a NIC → ~13 GB/s per node). vmhwm growth with rank count is the probe’s replicated input, not the engine; engine-side peak per rank is flat.

Noise floor

Single-shot campaign noise on the reference host is ±5–8% single-threaded and ±10–26% at 8–32 threads; untouched code moves that much between campaigns. Effects below that need the interleaved A/B protocol (scripts/ab-compare.sh), whose acceptance criterion is direction consistency across every pair.

Sources: research/HARDWARE.md (roofline tables, phase shares, bandwidth ceilings, partitioned and distributed tables, all achieved numbers and verdicts); benchmarks/PROFILING.md (byte model, interpretation rules, noise floor); the Rust API (layout and dispatch).

Against other tools

Manual: Truncation, Engine and propagation

Against PauliPropagation.jl, the other mature Pauli-propagation engine, the comparison is term for term: same circuits, same truncation, parity-gated timing. Against state-vector, stabilizer and MPO methods the tools do not overlap, and this page records where each applies.

Which method fits which problem

The sections below give the measured basis for each row. Everything under Methodology is the provenance behind the cross-engine numbers rather than part of choosing a method.

vs state-vector simulation

Not a competitor but a complementary oracle, and this suite uses it as one everywhere it reaches. Every exact reference on this site is a dense statevector (usually qiskit Aer); where two exact routes were affordable the reference is both of them, required to agree.

Two measured illustrations of where the boundary sits:

• Benchmark B needed an exact reference for a weight-10 observable whose causal cone is 30 qubits. Untruncated Pauli propagation over that cone exceeded a 26 GiB address-space cap at 4.3 × 10⁸ terms; the statevector over the same cone took ~150 s and does not care about depth.
• The same benchmark’s weight-17 observable has a 59-qubit cone: 2^59 amplitudes rules out any dense method, and untruncated propagation is far past the wall above. Neither method reaches it, which is why those references are self-converged and reported as not converged.

A state-vector simulator gives the answer; Pauli propagation gives the answer plus a truncation error to bound. That asymmetry is why every page on this site carries a convergence panel.

vs stabilizer (stim) simulation

Also an oracle, not a competitor. At a Clifford point stim gives the exact ±1 integer in under 0.1 s at any qubit count, and Benchmark A exists to be scored against it. Benchmark B reproduces those integers bit-exactly at every one of eight cutoffs, for three observables at both Clifford endpoints.

• stim: Clifford circuits, exactly, at enormous scale. Where a circuit is Clifford it is strictly the better tool.
• Pauli propagation: non-Clifford circuits, where the tableau method has nothing to say. The kicked-Ising kick angle separates the two: at θ_h ∈ {0, π/2} the circuit is Clifford and stim answers; at the hard interior angles the operator spreads over millions of Pauli strings.

A stabilizer simulator cannot serve as a noisy oracle: a tableau simulation samples one Pauli error rather than averaging over them, which is why Showcase B2 carries a hand-rolled Kraus density-matrix reference instead.

vs tensor-network / MPO methods

No measured head-to-head, and this site does not claim one. What it has is a cost-model comparison on the same operators: Showcase B6 computes the Pauli-spectrum entropy (the quantity governing truncation error for this engine) alongside the operator entanglement across a bipartition (the quantity governing MPO bond dimension), and finds them saying different things about the same operator: S_2 grows steadily with depth while S_op saturates around 1.3 nats from depth 5 on. A TDVP baseline at large n is a named limitation of Benchmark D, not silently approximated.

vs qiskit.SparsePauliOp / openfermion.QubitOperator

Different scope: these are Pauli-operator containers with algebraic manipulation, not propagation engines with truncation — no crossover concept applies. The committed comparison (baseline_comparison/README.md) benchmarks construction from string terms and one-layer Clifford conjugation, seeded inputs, n_terms ∈ {100, 1000, 10 000}; medians in µs, ratio = library / paulistrings:

openfermion has no equivalent conjugation operation and is not in the second group. PauliStrings.jl, the library that inspired this one, is excluded for the same reason PauliPropagation.jl is driven by subprocess: no PyJulia wiring anywhere.

vs PauliPropagation.jl

The comparison baseline is subprocess-driven and out of CI, pinned to PauliPropagation.jl 0.8.2 on julia 1.12.6 in a committed Project.toml/Manifest.toml. There is no PyJulia or juliacall anywhere: the entry points are a Julia script that reads a task JSON and emits a result JSON, and a subprocess wrapper that skips cleanly when no julia is on PATH.

Performance vs tracked-set size

There is no single ratio: the ranking changes sign, and where it changes sign depends on the workload by an order of magnitude.

The source is a dedicated head-to-head study (jl_performance/README.md): single-threaded core versus core, five interleaved abba pairs per configuration, accepted on direction consistency, never on a difference of two independently-noisy means. Every configuration passes a per-layer term-count parity gate before any timing is reported. ratio > 1 means this engine is faster.

Ranking crossover

The crossover spans 7× across these three workloads and the SU(4) matrix-gate path has none at all, so no single global crossover is quoted anywhere on this site.

Above the crossover

Memory, from the same study: process floors are 37.8 MB against Julia’s 0.601 GiB, a factor of 16; at the largest SU(4) configuration peak RSS is 0.239 GiB against 1.625 GiB, or 95 floor-subtracted bytes per peak term against 479. Both engines sample their own /proc/self/status, never a driver-side getrusage(RUSAGE_CHILDREN).

Below the crossover

Below the crossover jl’s hash-map backend is faster, by up to 3.6× at 68 terms: a hash-map insert per term costs little at 10² terms, while the bucketed per-layer pipeline costs nearly the same whatever the term count.

Methodology

How the PauliPropagation.jl comparison is made, and what had to be pinned before any number above could be quoted.

Parity discipline

Term-count parity blocks timing: no cross-engine wall time is reported for a configuration whose evolved Pauli sums diverge term-for-term at matched truncation, so every timed comparison runs the parity check first, untimed.

Four properties make the check strong:

1. One description, two engines. Both sides are driven from the same schema-v1 task JSON, built from the same recorded gate list; unknown keys, gate names and fields are hard errors on both sides.
2. Per applied layer, not just the final count: a divergence that cancels by the end is exactly the truncation-schedule bug the check exists to catch.
3. Term by term, not just the contracted expectation. Every gate name in the schema vocabulary (plus reversed-qubit variants of cnot and unitary_2q) gets a single-gate task compared coefficient by coefficient to 1e-12.
4. Along a sweep, not at one point: the parity legs use three cutoffs.

The results, all from committed benchmark READMEs:

The truncated rows of the last entry are non-vacuous: the same circuit with no policy keeps 3 881 terms, so the 1e-4 row exercises coefficient truncation and the max_weight=4 row exercises weight truncation.

Convention agreement

Both engines are Hermitian-Y: a real coefficient multiplies the literal Pauli string, Y carries no phase of its own, and the coefficient type stays real under every gate in the vocabulary. This was verified by hand-derived probe: S X S† = +Y and S† X S = −Y come out identically on both sides, and a cross-engine test encodes the sign.

Index conventions map cleanly: jl is 1-based with the leftmost Pauli character qubit 1, this repository 0-based with the leftmost character qubit 0. Observable keys map verbatim and gate indices map with a +1. Direction maps exactly: "heisenberg" ↔ heisenberg=true (jl’s default), "forward" ↔ heisenberg=false.

That mapping is not fully symmetric in practice: jl 0.8.2 defines no Schrödinger-picture transfer map at all for unitary_1q, unitary_2q, amplitude_damping, pauli_channel or depolarize2, so direction="forward" has no jl counterpart for any circuit using those gates. Every cross-engine comparison on this site is therefore run in the Heisenberg direction, where both sides agree.

Semantic divergences

Established by probes whose expected values are hand-derived in comments, never read back from the library.

The one real divergence: the coefficient boundary

jl truncates on abs(coeff) < min_abs_coeff, so it keeps a coefficient exactly equal to the threshold; this repository’s CoefficientThreshold keeps |c| > eps, so it drops it. Measured on both sides:

For generic angles the divergence is a measure-zero event and every parity row above passes untouched. It is not measure-zero for dyadic cutoffs at Clifford angles, where coefficients are exact dyadics and can land on the cutoff bit-exactly (Benchmark C hits this). The mitigation, when it bites: perturb the threshold on one side by one ulp and report it, never a coefficient. jl gets nextafter(eps, ∞); no float lies between, so jl’s |c| < eps′ becomes exactly this engine’s |c| <= eps. A test pins the divergence so a version bump cannot change it silently.

The second, narrower divergence: exact zeros

With min_abs_coeff = 0.0, abs(c) < 0 is never true, so jl keeps an exactly-zero coefficient; this engine’s merge kernels drop exact zeros unconditionally. Pinned with amplitude_damping(γ=1), whose X → √(1−γ)·X = 0 is bit-exact: this engine keeps 0 terms, jl keeps 1. Not measure-zero in practice, since Clifford-point angles produce exact cancellations. Mitigation for comparative runs: a strictly positive min_abs_coeff (any eps > 0 kills jl’s zeros too), stated in the results file.

Verified agreements

• The weight boundary: both engines keep weight == max_weight.
• When truncation is applied: jl truncates once per gate; this engine once per channel. The two coincide iff one gate object is one channel, which is the suite’s construction rule and structural in the schema. Measured: rz(0.05) on X with min_abs_coeff = 0.1 kills the sin branch immediately, so the second gate sees 1 term, not 2.
• Noise-channel parameter scales: jl damps by 1 − λ where this repository takes a probability p, so depolarize(p) maps to λ = 4p/3 and dephase(p) to λ = 2p; amplitude_damping(γ) is 1:1. jl has no native general Pauli channel or two-qubit depolarizing gate, so the runner builds each as a single diagonal-PTM gate with the exact dual — one gate, one truncation point. Both verified term by term.
• Two-qubit matrix ordering for unitary_2q is undocumented upstream; it was pinned against a known CNOT in both qubit orderings.

Noise-channel parity

Noise-channel semantics agree between the engines term by term, including the orientation of amplitude_damping under direction="heisenberg": the Heisenberg map is the unital dual Φ† (Φ†(I) = I), so ⟨Z⟩ for a qubit in |0⟩ — the channel’s fixed point — stays at 1. The shared fixture gives 9 terms on both engines with identical labels and bit-exact coefficients, a test pins the orientation from both sides, and Showcase B2 carries the physics.

Known gaps

• direction="forward" with unitary_1q, unitary_2q, amplitude_damping, pauli_channel, depolarize2: PauliPropagation.jl 0.8.2 defines no Schrödinger transfer map for those, and the runner rejects such a task up front. Every benchmark in the suite is Heisenberg.
• Non-computational, non-uniform product states: jl evaluates against |0…0⟩-style states only; this repository’s per-qubit label alphabet is strictly larger, so such a state cannot be compared against jl.
• Stim-sourced circuits must be expanded into an inline gate list on the Python side; jl has no Stim parser, and the runner makes that a hard error.
• topn truncation is absent from the interchange schema: jl has no equivalent. Likewise jl’s max_freq / max_sins truncations are excluded.
• jl’s experimental fused rotation kernel has no parity established, because it truncates during gate application.

Sources for this page: benchmarks/julia/README.md (the probes, the divergences, the gaps) and benchmarks/README.md, plus the per-benchmark READMEs linked inline. Cross-engine performance numbers: jl_performance/README.md (protocol and headline tables, engine 81c568a) and post-optimization/README.md. qiskit/openfermion numbers: baseline_comparison/README.md and its committed results.json.

Library

Signatures and tables for the Python API. Usage and rationale are in the Manual.

PauliString

One IXYZ string in the symplectic encoding, with no coefficient attached. Immutable and hashable; str and repr both print the label.

Constructors

Character i addresses qubit i, as in PauliSum.from_strings, and Y is the (x=1, z=1) key with no phase of its own. A qubit at or beyond num_qubits is a ValueError, as is a character outside IXYZ.

Accessors

Operations

All five require both operands to have the same num_qubits; anything else is a ValueError. The product string is mul’s either way, so a 0 coefficient is an exact algebraic zero. See Pauli strings for the encoding and worked examples.

Arithmetic

c must be a complex or real number — p * q between two PauliStrings is a TypeError pointing at .mul(other), the Pauli product with its phase, since * on a PauliString is scalar-only, matching PauliSum. This is what lets a sum be built directly out of strings, e.g. p("XYZ") * 2 + p("YZI") * 3; see The PauliSum type.

PauliSum

Constructors

num_qubits is inferred from the first label’s length when omitted from either from_strings form; an empty terms/labels with no explicit num_qubits is a ValueError, since there is nothing to infer it from. from_strings and from_arrays both use the crate’s Hermitian convention: a coefficient multiplies the literal Pauli string, and Y carries no phase of its own. Duplicate keys/rows/labels sum their coefficients; exact-zero coefficients are dropped.

from_arrays parameters:

Accessors

All five array/list accessors return the sum’s canonical order (partition-bucket index ascending, then lexicographic (x, z)), consistently across calls. str/repr show that same order’s prefix, never sorted by coefficient magnitude — a magnitude sort would cost O(len log len) just to print a preview, on a type whose whole point is staying cheap at a huge term count.

Arithmetic

Both operands of +/- must have the same num_qubits; anything else is a ValueError. Terms whose coefficients cancel exactly are dropped, as is every term when c is exactly zero. * is scalar-only: multiplying two sums is a full operator product, which this library does not implement.

Saving, loading, importing

See Saving and loading for the recipe.

Circuit

Constructor

Circuit(num_qubits) — an empty circuit.

Conventions

One gate per channel. Every Circuit method pushes exactly one channel per qubit/pair — a broadcast call like depolarize(p, [0, 1]) pushes two channels, not one bundled channel. Truncation runs after every channel, never after a bundle, so fusing gates changes the answer; build circuits one gate per call — see Truncation.

Rotation angle convention. Every rotation method implements U = exp(-i*theta*P/2) for its generator P. A generator coefficient s in exp(-i*s*P) (no /2) needs theta = 2*s.

Gates

Each of these is a Circuit method, a gates.<name>(...) factory returning a Channel for circuit.append(...), or both:

A channel with support on more than two qubits, other than pauli_rotation, makes propagate panic; there is no fallback path.

circuit.append(gates.h(0)) is equivalent to circuit.h(0).

Noise channels

Every noise channel is non-unitary: circuit.adjoint() raises ValueError naming it if the circuit contains one.

Coefficient rescale per channel: depolarize(p) scales every non-identity Pauli on its qubit by 1 - 4p/3; dephase(p) scales X and Y by 1 - 2p and leaves Z untouched; pauli_channel(px, py, pz) scales X by 1 - 2(py + pz), Y by 1 - 2(px + pz), Z by 1 - 2(px + py); depolarize2(p) scales every non-identity two-qubit Pauli on its pair by 1 - 16p/15.

pauli_channel(p/3, p/3, p/3, q) is depolarize(p, q); pauli_channel(0, 0, p, q) is dephase(p, q). See Noise channels for a worked example.

Composition

circuit.adjoint() satisfies obs.propagate(c, direction="heisenberg") == obs.propagate(c.adjoint(), direction="forward") to floating-point tolerance. Per gate: rz/rx/ry/pauli_rotation negate theta, s becomes sdg (and back), unitary_1q/unitary_2q conjugate-transpose, and h/x/y/z/cnot/cz/swap are self-adjoint. A non-unitary channel (a noise channel) raises ValueError naming it.

Importing a circuit