Backends

PRISM-Q ships nine CPU backends, an optional CUDA path attached to the statevector and stabilizer backends, and a feature-gated distributed statevector backend that shards the dense state across MPI ranks. The simulation engine picks a backend automatically (the density matrix, tensor network, and distributed backends are explicit-dispatch only), or you can select explicitly. For a task-oriented version of this material, see the Backends Deep Dive guide.

The diagrams below are rendered directly from PRISM-Q's own SVG circuit renderer.

Reset semantics

reset is the channel rho -> |0⟩⟨0| ⊗ tr_q rho on every backend: the qubit is traced out and replaced by |0⟩, leaving the rest of the register in the mixture the trace produces. Projecting onto |0⟩ and renormalizing is not equivalent. The two agree only when the reset qubit is unentangled; when it is entangled, projection also collapses its partners into the branch correlated with the |0⟩ outcome. Resetting qubit 1 of a Bell pair leaves ⟨Z0⟩ = 0 under the channel and ⟨Z0⟩ = 1 under projection.

A backend holding a single pure state cannot represent the resulting mixture, so it runs one trajectory of the channel: sample the measurement outcome, collapse onto it, and apply X when the outcome is 1. Averaged over shots that reproduces the channel, and a reset consumes one draw from the backend's RNG stream. The density-matrix backend holds the mixture and applies the channel directly, with no draw. tests/reset_channel.rs pins the contract across backends against the density-matrix oracle.

Memory budget

A circuit that does not fit in memory is an error, not a fallback. No backend silently hands the work to a different one when its state would not fit: it returns PrismError::IncompatibleBackend naming itself, the qubit count, the cap, and the environment variable that overrides it. Choosing a different backend is the caller's decision, and BackendKind::Auto makes it from circuit structure before any backend is constructed.

The check lives in Backend::init, which is the one point every execution path passes through before reserving its state. Putting it there means a caller that drives a backend directly, through run_on rather than simulate, gets the same guard as one that goes through dispatch.

The five growth caps are deliberately independent of PRISM_MAX_SV_QUBITS: the sparse, factored, MPS, tensor network, and factored stabilizer backends exist to run above the statevector cap, so lowering that cap to steer routing must not shrink what they may hold. Their defaults come from the same detected-memory budget.

The factored stabilizer cap is the one that is not a 2^n amplitude count. A stabilizer cluster costs O(n^2 / 64) words, so a dense cap is the wrong scale here: it would hold a cluster to the dense backends' qubit ceiling, far below the widths reached by the Clifford circuits at 128 qubits and above that dispatch selects that backend for.

The density-matrix cap is the tighter of its own override and half the statevector cap, computed in one place so dispatch-time validation and the backend's init guard cannot disagree about where the ceiling is. Raising PRISM_MAX_DM_QUBITS past that bound needs PRISM_MAX_SV_QUBITS raised with it, which is what the rejection says: the backend reports it itself rather than surfacing an error naming the statevector it allocates internally. When physical memory cannot be detected the caps are disabled and a warning is printed, because guessing a budget is worse than saying the budget is unknown.

Parallel noisy trajectories are the one path holding more than one state at a time: each Rayon thread runs its own backend, so peak memory is threads * state(n). That path is restricted to circuits below 14 qubits, where a statevector replica is 256 KiB and a full thread pool stays in the tens of megabytes. Above it trajectories run serially with one live backend, bounded by the ordinary state cap.

The three growth paths that once sat outside this contract are bounded at their growth events, each rejecting with an error naming its own backend before the allocation: a factored sub-state merge checks the merged block width against the statevector cap, the sparse map checks a branching gate's worst-case fan-out against the entry cap (so a rejection can fire one gate early on a state that would have deduplicated below it), and MPS gate application checks its live contraction-buffer total against the statevector budget. The MPS check bounds the workspace, not max_bond_dim itself: a large cap on a circuit whose bonds stay small is fine, and the rejection fires only when the bonds actually grow past what memory holds.

Statevector

Full-state simulation in a flat Vec<Complex64> of 2^n amplitudes. The primary backend for circuits up to ~28 qubits.

Gate kernels use enum dispatch with specialized routines for CX, CZ, SWAP, Cu, MCU, Rzz, BatchRzz, BatchPhase, DiagonalBatch, MultiFused, and PauliRot (one pass over (j, j ^ xmask) pairs for exp(-i θ P / 2), a parity-phase sweep when the string is Z-only; backends without the kernel receive the CNOT-ladder lowering from expand_pauli_rotations). Single-qubit gates go through PreparedGate1q with FMA-vectorized SIMD. MultiFused gates use a three-tier tiled kernel (L2 16K / L3 131K / individual passes) for cache locality. MultiFused batches where all gates are diagonal dispatch to a dedicated fast path (1 complex multiply/element vs 4+2 for full 2×2).

Rayon parallelism at ≥14 qubits with par_chunks_mut and MIN_PAR_ELEMS = 4096 per task. BMI2 _pext_u64 accelerates BatchPhase, BatchRzz, and DiagonalBatch LUT indexing.

Deferred measurement normalization: pending_norm accumulates normalization factors without full-state scaling passes. Zero-cost for circuits without measurements.

The Quantum Fourier Transform is a representative statevector workload, dense with controlled-phase gates that the fusion pipeline batches:

Stabilizer

Aaronson-Gottesman bit-packed tableau for Clifford circuits. O(n²) time and space. Scales to thousands of qubits. Gate kernels use wordwise bitwise ops and popcount for phase computation. Supports H, S, Sdg, SX, SXdg, X, Y, Z, Id, CX, CZ, SWAP, plus measurement, reset, and classical conditionals.

Word-group batching fuses multiple 1q gate flushes into single tableau passes. Type-grouped masks apply all gates of the same Pauli type with one wordwise op instead of per-gate dispatch. Sparse Generator Indexing (SGI) tracks per-qubit active generator lists, enabling targeted row operations instead of full-tableau scans. Lazy destabilizer materialization defers destabilizer rows until probabilities are requested.

Probability extraction uses coset-based enumeration with GF(2) Gaussian elimination. O(2^k) where k is the number of non-diagonal generators, rather than O(2^n).

Factored Stabilizer (FactoredStabilizerBackend): Per-cluster tableaux with dynamic merging. Starts with one qubit per cluster. Cross-cluster 2q gates merge tableaux. Measurement and reset can split independent sub-tableaux again. Independent subsystems avoid full-tableau work when product structure is preserved.

Sparse

HashMap<usize, Complex64> for states with few non-zero amplitudes. O(k) memory. Entries at or below a pruning threshold on |a|² (1e-16) are removed after gates that can shrink or cancel amplitudes, and the kept entries are rescaled so the state keeps its norm; a run whose threshold was raised above the default reports Approximate with a fidelity bound derived from the dropped weight. Best for circuits whose support stays concentrated in computational-basis states at large qubit counts.

The map's per-entry gate cost is about 16x the statevector's per-amplitude cost on a mixed diagonal and permutation workload (the sparse/densify bench rows), so a state that densifies past roughly 1/16 load factor runs slower than a dense vector at the same width would. There is deliberately no mid-run handoff to the statevector: automatic dispatch selects this backend only above the statevector memory cap, where the dense state exceeds the memory budget, and a run that branches past the entry cap rejects the gate rather than degrading silently. The map is keyed by a usize basis index, so a circuit wider than usize::BITS qubits is rejected at init and automatic dispatch sends it to MPS instead. An explicitly selected sparse run on a densifying circuit degrades in place, measured at up to 24x the dense cost when fully dense at 20 qubits.

MPS (Matrix Product State)

Chain of rank-3 tensors with adaptive bond dimension (default max 256). O(n·χ²) memory. Single-qubit gates absorb via FMA-vectorized SIMD over bond-dimension slices. Two-qubit gates contract adjacent sites, apply the gate, then SVD-truncate back. Non-adjacent gates route through SWAP chains.

Hybrid SVD dispatch: faer (bidiag+D&C) for matrices with m×n ≥ 256, hand-rolled Jacobi for small matrices.

Product State

Per-qubit [Complex64; 2] storage. O(n) memory, O(1) per single-qubit gate. Rejects entangling gates. Selected automatically for circuits with no 2q gates.

Shots and Pauli expectations answer from the per-qubit states rather than the 2^n probability vector, so both stay O(n) and the backend runs queries at widths no dense route reaches. See Sampling Architecture.

Tensor Network

Deferred contraction planned on metadata: a greedy min-size pass picks the pair order, seeded noisy restarts rerun it when the greedy tree's peak intermediate grows large, and the kernel replays the winner. Gates append tensors; contraction happens lazily at probability extraction, where the PRISM_MAX_PROB_QUBITS cap guards the dense readout and an explicit run past it errors naming the cap rather than reporting probabilities: None. Every contraction, dense or doubled, checks its planned peak intermediate against PRISM_MAX_TN_PEAK_QUBITS before allocating.

Measurement and reset do not contract to the dense state: the outcome draws from the single-qubit reduced density matrix and the renormalizing projector is absorbed into the tensor holding the measured qubit's output leg, so the network keeps its deferred form, mid-circuit measurement carries no width ceiling, and the tensor count does not grow across measurements.

Two further queries stay off the dense route by contracting the network against its conjugate. The bra copy's legs are shifted clear of the ket index space, and each qubit's boundary is either closed against its twin, which is a trace, or joined through an operator. A one-qubit reduced density matrix leaves that qubit's ket and bra indices open and returns a 2x2; a Pauli expectation joins every non-identity factor through its operator and contracts to a scalar. Both follow the doubled network's cost rather than the qubit count. An identity factor is a closed leg rather than an appended tensor, so a weight-k observable adds k tensors and not n.

Nothing about the planner changed: an index is open when exactly one tensor holds it, which the greedy ordering already carries through to its result.

The reduced density matrix is the half of general-noise support the backend was missing, so the trajectory engine now runs amplitude damping, phase damping, thermal relaxation, and custom Kraus channels here under explicit dispatch.

Factored

Dynamic split-state simulation. Starts with n independent 1-qubit states, merges via tensor product only when 2q gates bridge groups. Parallel kernels match statevector patterns for sub-states ≥14 qubits. Selected when subsystem decomposition detects partial independence.

Density Matrix

Exact mixed-state evolution. Stores the full density operator rho for n qubits as a 4^n Complex64 buffer laid out row-major: index (r << n) | c holds ⟨r|rho|c⟩. That layout is isomorphic to a 2n-qubit statevector whose high n qubits index the ket (row) and low n qubits index the bra (column), so gate application reuses the statevector kernels. A unitary U on the ket register gives the left product U rho; the right product rho U^dagger takes the gate's conjugate form on the bra register where one exists, and otherwise conjugates the buffer around the pass. So U rho U^dagger costs two statevector passes, plus two conjugations only for the variants with no conjugate form. Rzz is the exception that carries gate math of its own: both factors are diagonal, so the ket and bra phases cancel wherever the two registers agree on the target pair's parity, and the sandwich collapses to a single pass over a combined table.

Memory is 16 * 4^n bytes, so the ceiling is about 14 qubits on a 16 GiB host and 15 on 32 GiB (PRISM_MAX_DM_QUBITS moves it within the statevector budget). With a device attached (BackendKind::DensityMatrixGpu, or with_gpu on the backend) the buffer lives in VRAM, budgeted against free memory at init: an 11 GiB card holds 13 qubits, so the device lifts the ceiling by about one qubit and is a throughput arm rather than a width arm. On the device the unitary half still runs the dense statevector kernels over the embedded buffer, while the channels, projection, reset, and the diagonal and Pauli readouts have kernels of their own. Both kinds are explicit-dispatch only; Auto and AutoGpu never select them, and the device kind has no host fallback.

Selecting it with a noise model attached is the exact route for every Simulate terminal except the adjoint gradient: the mixture is evolved once and observables, marginals, probabilities, and shots all read that one evolution, and the parameter-shift gradient evaluates the mixture once per shifted angle. Readout error is the one part of a model that no evolution holds, so run and marginals reject a model carrying it and point at sample_counts. The adjoint stays excluded because it backpropagates against a pure state and a channel has no reverse evolution to walk. See Noise across the terminals for what that route accepts and what stays on trajectory averaging.

Pauli Path

Heisenberg propagation of an observable through a noisy circuit, as a weighted sum of Pauli strings. It holds no state of any kind, so it is not a Backend: like SPP and SPD it is an engine the dispatcher reaches directly, and it serves expectation_values and observable_expectation (and the parameter-shift gradient built on them) and nothing else. run, shots, marginals, and probabilities are rejected naming the two terminals it does serve.

The observable starts as one term and propagates backward. Clifford gates conjugate it one term at a time. Rz and Rzz split every anticommuting term into a cos branch and a -i sin branch, which is what grows the sum. A noise channel scales each term by the channel's action on the Pauli letters it touches, which is what shrinks it. The engine is polynomial exactly where the shrinking wins, so the useful regime is a circuit whose noise rate outpaces its density of non-Clifford rotations, and the term count at a given width is the thing to watch rather than the width itself.

Observable weight dominates that term count, ahead of both width and depth. On a two-layer hardware-efficient ansatz under 1% depolarizing, Z on one qubit holds the sum at 11 terms from 20 qubits to 100, while Z on two adjacent qubits reaches a 16384-term budget by 30 qubits and a full-width Z chain reaches it at every width. term_count_is_width_independent_at_unit_weight pins both halves.

Channels enter through their adjoint on the Pauli basis, not through a twirl, so nothing is approximated at the channel. A unital Pauli channel (Pauli, Depolarizing, PhaseDamping, TwoQubitDepolarizing) scales each letter by one eigenvalue. AmplitudeDamping and ThermalRelaxation are not unital: their adjoint sends Z to (1 - gamma) Z + gamma I, and the sum carries that identity branch as a second term. Custom Kraus, Kraus2q, and readout error have no Pauli-basis form and are rejected naming the density matrix.

With max_terms = 0 the run is exact and errors at the shared term ceiling rather than truncating silently. With a budget set, terms below epsilon are dropped once the sum exceeds it, and the discarded coefficient mass bounds the error: every channel and every Clifford conjugation is a contraction in the Pauli 1-norm, so a dropped term contributes at most its own magnitude to the terminal value. A run that truncated nothing reports itself exact whatever budget it was given.

State diagnostics

Simulate::reduced_density_matrix, Simulate::entanglement_entropy and Simulate::overlap read the output state once the circuit has been applied, so all three require a unitary circuit: a measurement, reset or conditional leaves one seeded branch of several, not the state the diagnostic is defined on. Each resolves to a single backend, as the native expectation path does, and asks that backend for the answer in its own representation; an explicitly selected backend with no kernel for one of them reports BackendUnsupported naming itself and the diagnostic rather than falling back to a dense export.

Under BackendKind::Auto the route is the dispatcher's choice, not the caller's, so a resolved backend that cannot answer is replaced by the statevector while the circuit fits its cap. A partially independent circuit routed to the factored backend and a sparse-friendly one both read their entropy that way, and the same circuits decline when the backend is named explicitly. A Clifford circuit keeps its tableau, which answers the entropy and the marginal without expanding anything.

The entropy is the von Neumann entropy in nats, so a Bell pair reads ln 2, and the Schmidt values come back descending with their squares summing to one whatever norm the representation carried. A stabilizer cut of rank r has 2^r equal weights, so its tableau reads the rank off one elimination and builds the list from it; past the dense export cap those values no longer fit while the rank still does, and the entropy comes back alone with EntropyResult::schmidt_values at None. That is the width where the old fallback to the statevector could not answer at all. The reduced density matrix is row major with side 2^k and trace one, and its 4^k entries are priced as a 2k-qubit statevector against the dense export cap. A noise model sends the marginal and the entropy to the density matrix, which answers the marginal of the exact mixture and declines the entropy.

Simulate::overlap takes a second seeded builder, so each side carries its own backend, seed and start state, and the two circuits must declare the same width. The result is the modulus squared of the inner product over the two normalized states; the amplitude is not reported, since a tableau keeps no global phase and every MPS truncation moves one. Normalization divides by both norms on every route, so an unnormalized chain answers the same as a unit-norm state. Two states in the same representation take the native route at any width, and every other pair is served by a dense export of both, which reaches exactly as far as the export cap does. A noise model on either side is rejected, since the fidelity of two mixtures is a different computation.

What a backend reports about its own result

Three Backend methods carry provenance onto every result: resolved names the engine, exactness says whether its representation can discard state weight and how much this run discarded, and placement says whether the state lived on the device. All three have defaults, so an out-of-tree backend compiles unchanged and is named by Backend::name.

These are reports, not predictions. exactness is read after the circuit has been applied, so the MPS bound reflects the singular values this run actually discarded, and placement reflects where the amplitudes ended up after any device fallback. The MPS accumulates discarded weight per SVD and returns 1 - total as a fidelity lower bound; the sum is over relative discarded weights, so the bound is conservative.

The decomposed route runs one backend per independent block and merges: its exactness is the weakest of the parts, its fidelity bound is the product, and its placement is Device only when every block was. Per-shot routes evolve one state per shot and keep the weakest claim across them, with the bound a minimum rather than a product.

The GPU backend is documented as a user guide. The distributed statevector backend is covered in the Capability and Support Matrix.