Private / Internal API

Largest spin for which "S=k/2" site types are pre-defined (S = 3/2 to 10 in half-integer steps; "S=1/2" and "S=1" are native to ITensors).

POSTSELECTION_PROB_TOL

Smallest Born probability $\lVert P\psi\rVert^2 / \lVert\psi\rVert^2$ of the branch a gate with needs_normalization(gate) == true may select (1e-14). The probability is the squared norm of the projected state relative to the input's, so the test is independent of how the input is normalized — the MPS backend does not renormalize after truncated unitary layers, so $\lVert\psi\rVert^2 < 1$ is common there and carries no physics. Below this threshold there is no normalized conditional state (renormalizing would only amplify rounding noise into a fake state), so the MPS and state-vector backends reject the operation with an ArgumentError and leave the state unchanged — see _check_postselection. SpinSectorMeasurement applies the same threshold to the total Born probability of its allowed sectors.

AbstractBackend

Abstract base type for all simulation backends. A backend is the mutable struct that owns the numerical representation of the quantum state; the user-facing SimulationState is parametrically typed on it (SimulationState{B<:AbstractBackend}), so backend selection is pure multiple dispatch.

Concrete backends: MPSBackend (ITensor MPS), StateVectorBackend (dense Vector{ComplexF64}), CliffordBackend (QuantumClifford.jl stabilizer tableau).

A new backend must implement, for SimulationState{MyBackend}:

• initialize!(state, init) for each supported AbstractInitialState
• _apply_single!(state, gate, phy_sites) — the gate-application primitive (either one generic method that resolves gates to matrices/operators, or per-gate methods plus an informative-ArgumentError fallback)
• born_probability(state, site, outcome) — required by the default measurement path and the BornProbability observable
• _measure_single_site!(state, site) — only if the default (Born-sample + Projection) implementation in src/Core/apply.jl cannot work for the representation
• a method (or an explicit ArgumentError rejection) for every exported observable callable — unhandled observables fall through to MPS-assumed generics and crash with a FieldError

The full contract — required fields, method tables per backend, RNG-stream expectations, and indexing conventions — is documented in docs/src/devdocs/backend_interface.md.

AbstractFeedback

Supertype for measurement-feedback specifications attached to a Measure gate via its feedback= keyword.

Concrete subtypes:

• OnOutcome: declarative outcome → gate map
• CallbackFeedback: wraps a raw function (state, sites::Vector{Int}, outcome) -> Any (raw functions passed to Measure(...; feedback=f) are auto-wrapped)

Feedback runs INSIDE with_guarded_stream(registry, :gates_spacetime): any attempt to draw from :gates_spacetime during feedback throws (fixed-draw contract — see RNGRegistry). Feedback randomness must come from :gates_realization (e.g. HaarRandom) or :born_measurement (nested measurements). Feedback gates do NOT advance engine counters (gate_idx) and emit no GateApplied events; a nested Measure still emits its MeasurementOutcome.

Abstract type for initial state specifications.

AbstractStaircase <: AbstractGeometry

Base type for staircase geometries with internal pointer.

CallbackFeedback(f)

Closure escape hatch for measurement feedback. f is called as f(state, sites::Vector{Int}, outcome) after the Born sample; its return value is ignored. sites is always a Vector{Int} (single-site Measure passes [site]).

Raw functions passed to Measure(...; feedback=f) are wrapped in this type automatically — you rarely construct it yourself.

Recursion (feedback applying another Measure) is allowed and is the user's responsibility to terminate.

CircuitBuilder

Internal mutable structure for recording circuit operations during do-block construction.

Users never see this type directly - they interact via:

circuit = Circuit(L=4, bc=:periodic) do c
    apply!(c, Reset(), SingleSite(1))
    apply_with_prob!(c; outcomes=[
        (probability=0.5, gate=PauliX(), geometry=SingleSite(1))
    ])
end

The builder records operations as NamedTuples which are then passed to the Circuit constructor.

Fields

• L::Int: Number of physical sites
• bc::Symbol: Boundary conditions (:periodic or :open)
• operations::Vector{NamedTuple}: Accumulated operation records
• params::Dict{Symbol,Any}: User-defined parameters passed from outer constructor

CircuitEvent

Abstract supertype for all typed circuit events recorded in a SimulationState's opt-in event log (see SimulationState(...; log_events=true)).

Concrete subtypes:

• GateApplied: a gate was executed on specific sites
• MeasurementOutcome: a Born-sampled measurement produced an outcome

Access recorded events via events and measurements.

Clifford (stabilizer) backend: holds a QuantumClifford.jl stabilizer tableau.

Uses MixedDestabilizer (tracks both stabilizer and destabilizer generators), which enables efficient $O(n^2)$ measurement via project!.

GateApplied(step, op_idx, element_idx, gate_label, sites) <: CircuitEvent

Recorded when a gate is applied during circuit execution.

Fields:

• step::Int: circuit step index (1-based; 0 if unavailable at emission site)
• op_idx::Int: index of the operation within circuit.operations
• element_idx::Int: index of the geometry element within the operation (1-based)
• gate_label::String: human-readable gate label (see gate_label)
• sites::Vector{Int}: physical sites the gate acted on

MPS-based backend: holds the underlying matrix product state, site indices, and truncation parameters (SVD cutoff, maximum bond dimension).

retruncate selects what happens to the bonds a projective gate (one with needs_normalization(gate) == true) did NOT touch:

• :never (default): only the bonds inside the gate's range are SVD-truncated, exactly as for a unitary gate; the rest of the chain is left as is until the next gate that touches it. The MPS is still renormalized after every projective gate.
• :gate: additionally run a full-chain truncate!(mps; cutoff) after every projective gate (two sweeps over all L sites per gate — the behavior before this option existed).

MeasurementOutcome(step, op_idx, sites, outcome) <: CircuitEvent

Recorded when a Born-sampled projective measurement produces an outcome.

Fields:

• step::Int: circuit step index (0 when emitted from the low-level measurement primitive, which has no step context — Task 9 engine rewire threads real indices)
• op_idx::Int: operation index (0 when unavailable, see above)
• sites::Vector{Int}: measured physical sites
• outcome::Int: measurement outcome (0 or 1 for qubits)

SentinelRNG <: Random.AbstractRNG

Poison-pill RNG installed by with_guarded_stream in place of a guarded stream. Any attempt to draw from it throws an ErrorException explaining that the stream is forbidden in the current context (e.g. :gates_spacetime during measurement feedback, whose draw count must stay data-independent — see the fixed-draw contract in RNGRegistry).

State-vector backend: holds the full statevector as a dense complex vector.

engine selects the gate-application algorithm used by _apply_single!:

• :optimized (default): Tier 2, hand-written in-place stride-loop kernels for 1- and 2-site gates (see apply_gate_sv_optimized! in src/StateVector/optimized.jl). No per-gate copy of ψ.
• :builtin: Tier 1, reshape + permutedims + matmul (see apply_gate_sv! in src/StateVector/StateVector.jl). Ground truth the kernels are verified against; the two engines agree to rounding (~1e-15), not bitwise.

Base.copy(circuit::Circuit) -> Circuit

Deep copy of a circuit, including its geometries. This is THE per-trajectory pattern for thread safety: staircase/Pointer geometries are mutable (their positions advance during simulate!), so concurrent trajectories must each own a private copy of the circuit.

copy is implemented as deepcopy, which preserves INTRA-circuit geometry aliasing: if the same staircase object appears in several operations of the original, the copies also share one (new) object — required for CIPT-style shared random-walk positions.

Threads.@threads for seed in seeds
    c = copy(circuit)                       # private geometry state
    st = SimulationState(...; rng=RNGRegistry(gates_spacetime=seed, ...))
    initialize!(st, ProductState(binary_int=0))
    simulate!(c, st; n_steps=n)
end

Base.getproperty(state::SimulationState{MPSBackend}, name::Symbol)

Property forwarding for the MPS backend — SUPPORTED API (not a deprecation shim). For state::SimulationState{MPSBackend}, the four property names

• state.mps → state.backend.mps
• state.sites → state.backend.sites
• state.cutoff → state.backend.cutoff
• state.maxdim → state.backend.maxdim

forward transparently to the MPSBackend fields, for both reads and writes (Base.setproperty! forwards the same four names). All other property names resolve to SimulationState's own fields.

This convenience is MPS-only by design: SimulationState{StateVectorBackend} and SimulationState{CliffordBackend} have no such forwarding (their payloads are reached explicitly via state.backend.ψ / state.backend.tableau), so accessing state.mps on them raises a FieldError — a loud signal that MPS-specific code received a non-MPS state.

rand(registry::RNGRegistry, stream::Symbol) -> Float64

Draw a random Float64 from the specified stream. Mid-level API.

randn(registry::RNGRegistry, stream::Symbol, dims...) -> Array{Float64}

Draw random normal array from the specified stream. For Haar random matrices.

randn(registry::RNGRegistry, stream::Symbol) -> Float64

Draw a random normal Float64 from the specified stream.

Define $\exp(i\pi S_z)$ operator for $S=1$ sites.

\[S_z = \operatorname{diag}(+1, 0, -1) \text{ in basis } |{+1}\rangle, |0\rangle, |{-1}\rangle\]

\[\exp(i\pi S_z) = \operatorname{diag}(e^{i\pi}, e^{0}, e^{-i\pi}) = \operatorname{diag}(-1, 1, -1)\]

This is a diagonal operator in the Sz basis.

_apply_single!(state::SimulationState, gate::AbstractGate, phy_sites::Vector{Int})

Apply gate to specific physical sites. Internal workhorse.

Steps:

1. Validate support matches site count
2. Convert physical sites to RAM indices
3. Build operator with physical site indices
4. Apply operator to MPS
5. Normalize iff needs_normalization(gate) (trait, Contract 3.5), by rescaling the gauge-center tensor left by step 4 (no full-chain inner product). Such a gate is rejected with an ArgumentError — before any tensor is written back, so the state is unchanged — when it annihilates the state (postselection onto a zero-probability outcome; see _check_postselection). The probability is judged relative to the input norm, ‖Pψ‖²/‖ψ‖²: the MPS is not renormalized after truncated unitary layers, so ‖ψ‖² alone says nothing about the physics.
6. With state.backend.retruncate == :gate, follow the normalization with a full-chain truncate!(mps; cutoff), recompressing the bonds the gate did not touch (a projection can lower the Schmidt rank at every bond). With :never (default) those bonds are left alone until the next gate that touches them, like after a unitary; this saves two O(L) sweeps per projective gate.

_apply_single!(state::SimulationState, gate::BondParity, phy_sites::Vector{Int})

Generic (MPS/state-vector-backend) rejection fallback: BondParity has no gate_matrix/build_operator representation and is only implemented on backend=:gaussian. The Gaussian backend provides the Gaussian-specific _apply_single!/execute! override (more specific on state), which Julia's dispatch prefers whenever the state actually uses GaussianBackend. The CliffordBackend's own AbstractGate fallback independently rejects BondParity on the Clifford backend with its own message.

_apply_single!(state::SimulationState, gate::GaussianHaar, phy_sites::Vector{Int})

Generic (MPS/state-vector-backend) rejection fallback: GaussianHaar has no gate_matrix/build_operator representation and is only implemented on backend=:gaussian. The Gaussian backend provides the Gaussian-specific _apply_single!(state::SimulationState{GaussianBackend}, ::GaussianHaar, ...) override; Julia's dispatch prefers that method (more specific on state) over this fallback whenever the state actually uses GaussianBackend. The CliffordBackend's own AbstractGate fallback (src/Clifford/Clifford.jl) independently rejects GaussianHaar on the Clifford backend with its own message (more specific on state there too).

_apply_single!(state::SimulationState{GaussianBackend}, gate::AbstractGate, phy_sites::Vector{Int})

Fallback for any gate NOT handled by one of the specific _apply_single! methods above. The Gaussian (free-fermion covariance-matrix) backend can only represent fermionic Gaussian operations; generic qubit gates (e.g. Hadamard, CNOT, Haar-random qubit unitaries, and PauliX — see the file header) are not Gaussian and have no covariance-matrix representation. Throws an informative ArgumentError naming the offending gate type and suggesting the dense-backend alternatives.

_apply_single!(state::SimulationState{GaussianBackend}, gate::BondParity, phy_sites::Vector{Int})

Dispatch disambiguator (always throws). BondParity is a projective measurement: on the Gaussian backend it is executed through the execute! measurement protocol (Gaussian execute! override), never through the _apply_single! gate path. This method exists so the AbstractGate catch-all below (specializing on state's type parameter) is not ambiguous against the un-parameterized _apply_single!(state::SimulationState, gate::BondParity, ...) rejection fallback in src/Gates/bond_parity.jl (specializing on gate's type) — the same ambiguity bug class previously fixed for the StateVector/Clifford backends (see the disambiguating overrides in that file).

_apply_single!(state::SimulationState{GaussianBackend}, gate::GaussianHaar, phy_sites::Vector{Int})

Apply a Haar-random SO(n) Majorana rotation to the two sites in phy_sites, where n is the total number of Majorana indices carried by the two sites (resolved granularity-aware via site_majoranas):

• fermionic-mode granularity (default, majoranas_per_site == 2): the two sites carry 4 Majoranas ix = [2a-1, 2a, 2b-1, 2b] → Haar-$SO(4)$.
• Majorana-chain granularity (site_type="Majorana", majoranas_per_site == 1): the two sites ARE two Majoranas ix = [a, b] → Haar-$SO(2)$. Haar on SO(2) is EXACTLY the uniform-angle rotation $\exp(\theta\,\gamma_a\gamma_b)$ with $\theta \sim U[0,2\pi)$ ($SO(2)\cong U(1)$, Haar = uniform angle) — the class-DIII unitary K_U of Pan, Shapourian, Jian, arXiv:2411.04191 (Eq. S-III.1; Python reference GTN.measure_all_tri_op's Υ = kraus((0, cos φ, sin φ)) branch). The φ ↔ rotation convention, matching the reference Python implementation (cross-checked in test/gaussian/test_majorana_chain.jl): contracting kraus((0, cos φ, sin φ)) on the Majorana pair (a, b) equals direct conjugation $\Gamma \leftarrow R\,\Gamma\,R^T$ with $R = \begin{pmatrix}\cos\varphi & -\sin\varphi\\ \sin\varphi & \cos\varphi\end{pmatrix}$ on rows/columns (a, b) (exact to machine precision). Since $\varphi \sim U[0,2\pi)$ makes R uniform over SO(2), the two parameterizations define the SAME ensemble.

The orthogonal matrix O = haar_orthogonal(rng, n) is drawn from the :gates_realization RNG stream (one draw per application, mirroring RandomClifford on the Clifford backend) and conjugated DIRECTLY onto the covariance matrix at the Majorana rows/columns ix:

\[\Gamma[\mathrm{ix},:] = O\,\Gamma[\mathrm{ix},:] \\ \Gamma[:,\mathrm{ix}] = \Gamma[:,\mathrm{ix}]\,O^T\]

i.e. $\Gamma \leftarrow R\,\Gamma\,R^T$ with $R = O \oplus I$ — exact for Gaussian unitaries (no contraction kernel involved). The result is re-antisymmetrized and, if the purity diagnostic $\max\lvert\mathrm{diag}(\Gamma^2)+1\rvert$ exceeds state.backend.purify_tol, re-purified via purify!.

_apply_single!(state::SimulationState{CliffordBackend}, gate::AbstractGate, phy_sites::Vector{Int})

Fallback for any gate NOT handled by one of the specific _apply_single! methods above. The Clifford (stabilizer-tableau) backend can only represent Clifford-group operations; non-Clifford gates (e.g. arbitrary rotations, Haar-random unitaries, or non-Clifford projections) have no native tableau representation. Throws an informative ArgumentError naming the offending gate type and suggesting the dense-backend alternatives.

_apply_single!(state::SimulationState{CliffordBackend}, gate::CNOT, phy_sites::Vector{Int})

Apply a CNOT gate to the stabilizer tableau via QuantumClifford's sCNOT.

_apply_single!(state::SimulationState{CliffordBackend}, gate::CZ, phy_sites::Vector{Int})

Apply a controlled-Z gate to the stabilizer tableau. NOTE: QuantumClifford.jl does NOT have a function named sCZ — the correct symbolic gate is sCPHASE, verified (via stabilizer conjugation of a +X_ state) to match the standard CZ = diag(1,1,1,-1) convention.

_apply_single!(state::SimulationState{CliffordBackend}, gate::Hadamard, phy_sites::Vector{Int})

Apply a Hadamard gate to the stabilizer tableau via QuantumClifford's sHadamard.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PauliX, phy_sites::Vector{Int})

Apply a Pauli-X gate to the stabilizer tableau via QuantumClifford's sX.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PauliY, phy_sites::Vector{Int})

Apply a Pauli-Y gate to the stabilizer tableau via QuantumClifford's sY.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PauliZ, phy_sites::Vector{Int})

Apply a Pauli-Z gate to the stabilizer tableau via QuantumClifford's sZ.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PhaseGate, phy_sites::Vector{Int})

Apply an S (phase) gate to the stabilizer tableau via QuantumClifford's sPhase.

_apply_single!(state::SimulationState{CliffordBackend}, gate::RandomClifford, phy_sites::Vector{Int})

Sample a random n-qubit Clifford operator from the :gates_realization RNG stream via QuantumClifford.random_clifford(rng, gate.n), then apply it NATIVELY to the stabilizer tableau (via QuantumClifford.apply!(tableau, qubit_indices, op)) — no dense-matrix conversion (no QuantumOpticsBase involvement), unlike the MPS/state-vector backends' gate_matrix path.

_apply_single!(state::SimulationState{CliffordBackend}, gate::SWAP, phy_sites::Vector{Int})

Apply a SWAP gate to the stabilizer tableau via QuantumClifford's sSWAP.

_apply_single!(state::SimulationState{StateVectorBackend}, gate::AbstractGate, phy_sites::Vector{Int})

Apply gate to specific physical sites on a dense state-vector backend. Internal workhorse for both engines (see apply_gate_sv! and apply_gate_sv_optimized!).

Steps:

1. Validate support matches site count
2. Convert physical sites to RAM indices (identity for SV, kept for code-path consistency with the MPS backend)
3. Resolve the gate's dense matrix (HaarRandom needs rng + local_dim)
4. Apply via apply_gate_sv_optimized! (Tier 2, the default, in place) or apply_gate_sv! (Tier 1, state.backend.engine == :builtin, returns a fresh vector), then store the result in state.backend.ψ
5. Normalize iff needs_normalization(gate) (trait, Contract 3.5) — there is NO truncate! equivalent for state vectors (no bond dimension). Such a gate is rejected with an ArgumentError that leaves the state unchanged when it annihilates the state — postselection onto a zero-probability outcome, see _check_postselection. The guard gets the Born probability ‖Pψ‖²/‖ψ‖² (same contract as the MPS backend); a state vector stays normalized, so the ratio only differs from ‖Pψ‖² by rounding here. Because the Tier 2 kernel mutates ψ in place, the guard runs BEFORE it when ‖Pψ‖² has a closed form (_sv_projected_norm2; Projection, i.e. Measure/Reset), and the kernel works on a copy of ψ otherwise (SpinSectorProjection, user-defined projective gates). Tier 1 always returns a fresh vector, so it checks after applying.

_bit_pattern_string(init::ProductState, L::Int, local_dim::Int) -> String

Derive the physical-site bit-pattern string (site 1 = MSB, site L = LSB) from a ProductState's binary_int/binary_decimal/bitstring specification, padded or truncated to L characters (a binary_int that needs more than L binary digits is rejected with an ArgumentError instead of truncated, so every backend sees the same L-site state). Shared verbatim by all three backends' initialize!(::SimulationState, ::ProductState) methods (MPS: src/State/initialization.jl, state-vector: src/StateVector/initialization.jl, Clifford: src/Clifford/initialization.jl) — previously copy-pasted identically in all three.

Does NOT handle init.spin_state — callers must branch on spin_state before calling this function, since each backend treats it differently (MPS/state-vector build a uniform product state from it directly; Clifford rejects it since it is qubit-only). Assumes exactly one of binary_int/binary_decimal/bitstring is set (enforced by the ProductState inner constructor); the final else branch is an unreachable-in-practice defensive fallback preserved verbatim from the pre-refactor code.

local_dim is currently UNUSED by the derivation itself (only L affects padding/truncation) but is accepted as part of the signature for forward compatibility with an arbitrary spin-S generalization of ProductState (planned future work) — kept so that extension does not require touching call sites again.

_build_gate_operator(state, gate, phy_sites, ram_sites) -> ITensor

Build the operator tensor for the gate.

_check_postselection(norm2_out::Real, norm2_in::Real, gate::AbstractGate, sites::Vector{Int})

Guard shared by the MPS and state-vector _apply_single! methods for gates with needs_normalization(gate) == true: norm2_out is the squared norm $\lVert P\psi\rVert^2$ of the state after applying gate on the physical sites, norm2_in the squared norm $\lVert\psi\rVert^2$ of the state before. Their ratio is the Born probability of the requested branch whatever the input's normalization. Throws an ArgumentError naming the gate, the sites, and both norms when the ratio is below POSTSELECTION_PROB_TOL — or when norm2_in is zero, where no probability is defined at all; returns nothing otherwise. Callers run it before committing the new tensors/vector, so a rejected gate leaves the state exactly as it was.

_check_string_order_site_type(state::SimulationState)

StringOrder is defined for spin-1 chains only: its $S_z$ and $\exp(i\pi S_z)$ factors are the spin-1 operators, so on any other site type the formula is not the advertised observable. Throw an ArgumentError unless state.site_type == "S=1". Shared by the MPS and state-vector methods.

_contract_op_block!(mps::MPS, op::ITensor, sites::Vector{Index}) -> (i_list, block)

First phase of apply_op_internal!: orthogonalize mps to the first RAM site touched by op, contract the tensors spanning the operator's RAM range i_list (sorted) into one block, and apply op to it. Mutates only the gauge of mps — no tensor is replaced — so the represented state is unchanged until _write_op_block! commits the block. Because the MPS is gauged at i_list[1], norm(block) is the norm of the full post-operator state.

_ee_check_common(renyi_index, threshold, base) -> Float64

Shared construction-time validation of the non-cut parameters, identical for the bipartition and region forms. Returns the NORMALIZED renyi_index (a Float64), which both callers store in the struct's Float64 field.

renyi_index is NORMALIZED THEN VALIDATED, in that order:

1. Bool is rejected outright — true/false are Real in Julia but are not meaningful Rényi indices, and Float64(true) == 1.0 would otherwise be silently accepted as von Neumann.
2. The value is converted to Float64 (conversion failure is re-thrown as an ArgumentError, never a bare InexactError/MethodError).
3. The CONVERTED value must be finite and > 0.

Validating after normalization is what keeps the struct's Float64 field invariant honest for wide inputs: BigFloat("1e400") converts to Inf and BigFloat("1e-400") to 0.0, so both are rejected here instead of being stored as a nonsense index. n = 0 (Hartley) stays excluded.

_ee_validate_region(region::Vector{Int}, L::Int) -> nothing

Call-time validation of a stored (sorted, duplicate-free, positive) region against the system size L, shared by every region-capable backend method.

Throws an ArgumentError when

• any site exceeds L (maximum(region) > L), or
• the region covers the whole system (length(region) >= L) — a trivial bipartition with zero entropy, mirroring the bipartition path's 1 <= cut < L call-time check.

_evaluate_recording(record_when, ctx, step_idx, n_steps) -> (set_flag::Bool, record_now::Bool)

Evaluate recording criteria and return whether to set the recording flag and/or record immediately.

For compound geometries (called inside element loops), :every_gate triggers immediate recording. For simple geometries, the caller handles the immediate recording after setting the flag.

Returns a tuple:

• set_flag: Whether to set should_record_this_step = true
• record_now: Whether to call record!(state) immediately

_gaussian_region_majoranas(state, region) -> Vector{Int}

Map a (sorted, duplicate-free) collection of PHYSICAL sites to the sorted vector of its Majorana indices, granularity-aware via site_majoranas: fermionic-mode granularity — each site s contributes $(2r-1, 2r)$ with r = state.phy_ram[s] (2·|region| indices); Majorana-chain granularity (site_type="Majorana") — each site contributes its own single Majorana index (|region| indices).

_haar_unitary(N::Int, rng::AbstractRNG) -> Matrix{ComplexF64}

Haar random N×N unitary via QR of a complex Ginibre matrix. EXACT reproduction of the CT.jl U() algorithm (lines 585-592): two separate randn(rng, N, N) calls (real then imaginary parts), QR, then phase correction by diag(R)/|diag(R)|. Do NOT change the call shape/order — the :gates_realization stream consumption is an API contract (bit-exact regression tests).

_infer_matrix_gate_dims(N::Int) -> (d, n)

Infer (local dimension, site count) from a MatrixGate size N per the documented rule: 2^n (n ≥ 1) → qubits; 3^n (n ≥ 2) → spin-1; otherwise ArgumentError.

_is_static_geometry(geo::AbstractGeometry) -> Bool

Trait: true for geometries whose elements(geo, L, bc) enumeration is provably STEP-INVARIANT (immutable structs with no advancing/movable position): Bricklayer, AllSites, EachSite, Sites, SingleSite, AdjacentPair.

false for mutable geometries — StaircaseLeft/StaircaseRight (positions advance during simulate!) and Pointer (moved explicitly via move!) — and, CONSERVATIVELY, for any geometry type not explicitly listed here (the fallback is false: correctness over performance). Any future compound geometry containing a member must only return true if EVERY member is static.

Used by simulate! to decide whether an op's elements() result may be cached across steps within one simulate! call. Never used to cache mutable geometries — no invalidation scheme exists by design.

_kraus(n::NTuple{3,<:Real}) -> Matrix{Float64}

$4\times4$ Gaussian gate tensor $\Upsilon$ parameterized by $n = (n_1,n_2,n_3)$ — port of GTN.kraus (GTN.py:105-118) with c = [1,1,1]. Unitary 2-Majorana rotation for $\|n\| = 1$; parity projection for $n = (\pm1,0,0)$.

_local_basis_vector(site_type::String, local_dim::Int, name::String) -> Vector{Float64}

Resolve a single-site ITensor state name (e.g. "0", "1", "Up", "Z0", "Dn") into a dense local basis vector, WITHOUT constructing a full MPS. Used to build a state-vector-backend product state site-by-site.

Verified outputs: Qubit "0"->[1,0], "1"->[0,1]; S=1 "Up"->[1,0,0], "Z0"->[0,1,0], "Dn"->[0,0,1].

_matrix_gate_site_count(N::Int, d::Int) -> n

Exact site count n with d^n == N (n ≥ 1), or ArgumentError if N is not a positive power of d. Integer arithmetic only (no float log rounding).

_measure_single_site!(state::SimulationState, site::Int) -> Int

Perform Born-sampled projective measurement on a single site. Returns the measurement outcome (a level index 0 .. local_dim-1; 0 or 1 for qubits).

This is the FUNDAMENTAL measurement operation:

1. Compute Born probabilities $P(k \mid \psi)$ for levels k = 0 .. local_dim-2
2. Sample ONE categorical outcome using a single scalar draw from the :born_measurement RNG stream (at local_dim=2 this reduces EXACTLY to the historical binary draw rand < P(0) ? 0 : 1 — same draw count, same float comparison, bitwise-identical qubit trajectories)
3. Apply the per-level Projection operator
4. Return outcome (for conditional logic in Reset / feedback)

_measure_single_site!(state::SimulationState{GaussianBackend}, site::Int) -> Int

Override the default (Projection-based) _measure_single_site! for the Gaussian backend. Born-samples the occupation of site and collapses the covariance matrix onto the observed parity sector via parity_projection! (the closed-form rank-2 evaluation of gaussian_contraction! with the projector parity_projection_upsilon), mutating state.backend.corr in place.

Draw contract (REDUNDANT-DRAW, cross-backend lockstep): exactly ONE scalar :born_measurement draw is consumed per measured site — always, unconditionally, BEFORE any probability computation — matching the generic MPS/SV _measure_single_site! (Core/apply.jl) and the Clifford override (src/Clifford/measurement.jl) draw-for-draw. The outcome uses the shared threshold convention outcome = r < p₀ ? 0 : 1; probabilities are continuous here, so no deterministic branch is needed — the draw is always genuinely consumed by the comparison (for a deterministic site p₀ ∈ {0,1} and the comparison is vacuous, i.e. the draw is redundant, exactly as on the Clifford backend).

Outcome → parity sign mapping: parity_projection!(Γ, a, b, s) projects onto the parity sector $i\gamma_{2i-1}\gamma_{2i} = -s$ and leaves the post-measurement state with $\Gamma[2i-1,2i]=-s$. To end with Γ[2i−1,2i] = +1 (outcome 0, unoccupied) we project with s = −1; for Γ[2i−1,2i] = −1 (outcome 1, occupied), s = +1. I.e. $s = 2\cdot\mathrm{outcome}-1$.

Throws ArgumentError if the sampled outcome has probability ≤ 1e-15 (probability-zero branch — cannot be projected onto); parity_projection! additionally raises its own vanishing-state ArgumentError as a second line of defense.

Logs a MeasurementOutcome event exactly like the default implementation (Core/apply.jl) does, and returns the outcome (0 or 1).

_measure_single_site!(state::SimulationState{CliffordBackend}, site::Int) -> Int

Override the default (Projection-based) _measure_single_site! for the Clifford backend. Uses QuantumClifford.projectZ! to measure site site in the Z basis, mutating state.backend.tableau in place.

Draw contract (REDUNDANT-DRAW, cross-backend lockstep): exactly ONE scalar :born_measurement draw is consumed per measured site — always, unconditionally — matching the generic MPS/SV _measure_single_site! (Core/apply.jl) draw-for-draw:

• Deterministic case (anticom_index == 0): the outcome is read directly off the returned phase byte. The draw is REDUNDANT here (the stabilizer formalism already fixes the outcome) and its value is DISCARDED — it is consumed purely to keep the :born_measurement stream position identical to the MPS/SV backends, where the (equally ignored) draw is structurally unavoidable.
• Undetermined case (anticom_index > 0): the outcome is genuinely random (50/50) and uses the same draw via the shared rand(rng) < p₀ ? 0 : 1 threshold convention (p₀ = 0.5 exactly for a stabilizer state); the phase of the anticommuting stabilizer row is set accordingly.

This deliberate redundant draw is what guarantees ABSOLUTE cross-backend reproducibility: same seeds ⇒ same measurement record on MPS, state-vector, and Clifford alike (test/audit/cross_backend.jl (b), test/audit/born_measurement.jl (e)).

Logs a MeasurementOutcome event exactly like the default implementation (Core/apply.jl) does, and returns the outcome (0 or 1).

_mi_contiguous_region(r, which::String) -> UnitRange{Int}

Normalize a region spec (Int, range, or vector of integers) into a UnitRange{Int}, throwing an informative ArgumentError unless it is a non-empty, positive, strictly contiguous ascending set of sites. Used where contiguity is a CONSTRUCTION-time requirement (e.g. TripartiteMutualInformation's region layout).

_mi_entropy_from_probs(p, n, base, threshold) -> Float64

Entropy from a Schmidt/RDM probability spectrum p (need not be normalized; tiny negative eigenvalues from numerical RDMs are clamped at threshold^2).

$n=1$ gives the von Neumann entropy $-\sum q\log_b(q)$; a general real $n \neq 1$ (a normalized Float64) gives the Rényi entropy $\log_b(\sum q^n)/(1-n)$ via a LOG-DOMAIN evaluation. This is a deliberately separate copy of the same three-branch structure used by _von_neumann_entropy (src/Observables/entanglement.jl) and the state-vector bipartition kernel (src/StateVector/entanglement.jl) — the bodies are independent, only the branch thresholds _RENYI_SHUNT / _RENYI_NEAR1 are shared, so every kernel switches branches at the same place.

_mi_is_contiguous(v::Vector{Int}) -> Bool

true iff the (sorted) site vector is a plain ascending unit-step range.

_mi_region_vector(r, which::String) -> Vector{Int}

Normalize a region spec (Int, range, or vector of integers) into a sorted, duplicate-free Vector{Int} of positive sites, throwing an informative ArgumentError if empty, repeated, or non-positive. Contiguity is NOT required here — it is enforced per backend at evaluation time (see _validate_mutual_information); the Gaussian backend accepts arbitrary site subsets.

_mi_validate_bounds(mi::MutualInformation, state) -> nothing

Evaluation-time bounds check shared by ALL backends (including Gaussian): every region site must lie within 1:L. Regions are stored sorted, so checking last suffices.

_mps_norm2(state::SimulationState{MPSBackend}) -> Float64

Squared norm $\langle\psi|\psi\rangle$ of the backend MPS, used by the MPS observables that contract $\langle\psi|O|\psi\rangle$ with inner (PauliString, StringOrder, DomainWall) to return the physical expectation value $\langle\psi|O|\psi\rangle / \langle\psi|\psi\rangle$.

The MPS backend does NOT renormalize after unitary gates (only after gates with needs_normalization(gate) == true), so a truncating unitary layer leaves $\langle\psi|\psi\rangle < 1$. The retained state is still a valid (unnormalized) representation, and every reported expectation must refer to its normalized version — born_probability, Magnetization (via ITensorMPS.expect), EntanglementEntropy, and MutualInformation already divide out the norm; this helper gives the inner-based observables the same contract.

Non-mutating. Cheap when the MPS has a well-defined orthogonality center (the norm of a single tensor), a full $\langle\psi|\psi\rangle$ contraction otherwise. Throws an ArgumentError on a zero-norm MPS, for which no expectation value is defined.

_mps_subset_rdm_probs(state, phys_sites) -> Vector{Float64}

Eigenvalues of the reduced density matrix of the given PHYSICAL sites of an MPS state. Sites are mapped to RAM positions via phy_ram (handling the folded-PBC layout transparently); the RDM is built by sweeping the RAM span lo:hi, keeping open (primed) site indices on the requested sites and tracing all others. Left/right environments collapse to identity because the MPS is orthogonalized to lo (sites < lo left-orthonormal, sites > lo right-orthonormal).

_parse_spin_site_type(site_type::AbstractString) -> Union{Nothing, Rational{Int}}

Parse a spin site-type string of the form "S=<n>" or "S=<k>/2" and return the spin S as a Rational{Int} (e.g. "S=1" → 1//1, "S=3/2" → 3//2). Returns nothing for any string that is not a spin site type (e.g. "Qubit", "Qudit"), so callers can use it as a pattern test.

_product_region(pg::ProductGate, L::Int, bc::Symbol) -> Vector{Int}

Canonical region of a ProductGate: the sorted union of elements(pg.region_geometry, L, bc). This is the ONLY region a ProductGate may be applied to (API contract; see the ProductGate docstring).

_product_state_vector(site_type, local_dim, state_names_physical) -> Vector{ComplexF64}

Build the full dense product-state vector via Kronecker product of the L per-site local basis vectors, with site 1 FIRST in the kron chain (i.e. kron(v1, v2, ..., vL)). This makes physical site 1 the slowest-varying/MSB tensor index, matching the MPS backend's documented MSB convention (bit_pattern_str[1] is site 1's bit) and the SAME Kronecker convention documented in src/Gates/matrix_gate.jl (used by gate_matrix).

_projection_matrix(outcome::Int, d::Int) -> Matrix{ComplexF64}

Dense d×d projector $\lvert \text{outcome} \rangle\langle \text{outcome} \rvert$ (0-based level index), validated against the local dimension d.

_purity_drift(Γ::AbstractMatrix{Float64}) -> Float64

Purity diagnostic $\max_i\lvert(\Gamma^2)_{ii} + 1\rvert$ (zero for a pure state, $\Gamma^2 = -I$) of an antisymmetric covariance matrix, computed from the diagonal only: $(\Gamma^2)_{ii} = \sum_j \Gamma_{ij}\Gamma_{ji} = -\sum_j \Gamma_{ji}^2$, i.e. $\max_i\lvert 1 - \lVert\Gamma_{:,i}\rVert^2\rvert$ — $O(N^2)$, allocation-free, and a contiguous (vectorizable) column sweep, instead of forming the full $O(N^3)$ product $\Gamma\Gamma$. Same quantity as maximum(abs.(diag(Γ * Γ) .+ 1)) up to summation order and up to any roundoff-level asymmetry of Γ (which enters only at first order, far below any sensible purify_tol).

_random_clifford_unitary(n::Int, rng::AbstractRNG; local_dim::Int=2) -> Matrix{ComplexF64}

Sample a Haar-random n-qubit Clifford operator via QuantumClifford.random_clifford(rng, n) and convert it to a dense local_dim^n × local_dim^n unitary matrix.

Conversion path: QuantumClifford.CliffordOperator → QuantumOpticsBase.Operator (via the QuantumCliffordQOpticsExt package extension, triggered by having both QuantumClifford and QuantumOpticsBase loaded) → dense .data matrix. This is the officially documented conversion route in QuantumClifford.jl (no direct Matrix(::CliffordOperator) method exists in QuantumClifford.jl itself).

_renyi_scaled_tails(logs, n::Real) -> (t_n, s)

Rescaled tail sums of a log-probability vector, shared by the GENERAL-n branch of all three entropy kernels (and by nothing else — no n = 1 path calls it).

With lmax = maximum(logs) attained at index imax, returns

\[t_n = \sum_{k \neq i_{\max}} \exp\big(n(\ell_k - \ell_{\max})\big) = \sum_{k \neq i_{\max}} (p_k/p_{\max})^n\]

\[s = \sum_{k \neq i_{\max}} \exp(\ell_k - \ell_{\max}) = (\textstyle\sum p - p_{\max})/p_{\max}\]

Every exponent is ≤ 0, so neither sum can overflow at ANY accepted n, and both lie in [0, length(logs) − 1].

The max index is located, and the sums accumulated, with PLAIN INTEGER indexing on purpose. logs is not always a Vector: on the MPS path it is an NDTensors.DenseTensor (from diag(S) of the bond SVD), for which eachindex yields CartesianIndex{1} while argmax returns an Int64 — so the natural for k in eachindex(logs) if k != argmax(logs) filter silently excludes NOTHING and leaves the maximum in the tail, inflating a flat rank-2 spectrum's Rényi entropy from log 2 to log 3. Integer indexing is consistent across both conventions.

_reset_circuit_geometries!(circuit::Circuit)

Reset all staircase geometry positions to their initial values. Called at the start of simulate! to ensure deterministic behavior.

_s_dot_s(s::Rational) -> Matrix{Float64}

The $S_1 \cdot S_2$ operator for two spin-s particles as a (2s+1)² × (2s+1)² matrix in the descending-m tensor product basis (generic version of s1_dot_s2, built from spin_operators(s)).

_spin_m_label(m::Rational) -> String

Format a magnetic quantum number m for use in "Z<m>" state names: integers render without denominator ("1", "-2"), half-integers as "<2m>/2" fractions ("3/2", "-1/2").

_staircase_max_position(geo::AbstractStaircase, L::Int) -> Int

Largest position an open-boundary staircase may occupy: L - range, so the pair (pos, pos+range) always fits inside 1:L. Throws ArgumentError when no position fits (range >= L).

_sv_digit(n0, stride, d) -> Int

Extract the base-d digit of the 0-indexed basis integer n0 at the physical site whose precomputed stride is stride = d^(L - s) (site 1 = MSB convention shared by all state-vector observables). Callers hoist the loop-invariant d^(L - s) power out of their inner loops and pass it in.

_sv_matvec!(out, U, buf, m) -> out

out = U * buf for the m-dimensional gate block (m = d^n), as explicit loops. mul! would dispatch to a BLAS gemv for every block, and the per-call overhead of that dominates a 2×2 or 4×4 product; the stride kernels call this once per block, d^(L-n) times per gate.

_sv_projected_norm2(gate, state::SimulationState{StateVectorBackend}, phy_sites) -> Union{Float64, Nothing}

Squared norm $\lVert P\psi\rVert^2$ of the state after gate acts on phy_sites, computed WITHOUT applying the gate, or nothing when gate has no closed form (the caller then applies it to a copy of $\psi$). Projection(k) keeps exactly the amplitudes whose digit at the site is k, so its value is the Born probability born_probability(state, site, k): one read-only pass over $\psi$, no allocation. This is the Measure/Reset hot path.

_sv_spin_magnetization(state, s::Rational) -> Float64

Spin-site state-vector Magnetization: $M_z = \frac{1}{L}\sum_i \langle S_{z,i}\rangle$ with $\langle S_z\rangle = \sum_k (S-k)\,P(\mathrm{digit}=k)$, summed over all 2S+1 levels (level k, 0-based, has eigenvalue m = S-k; site 1 = MSB digit). One O(L·dᴸ) pass — spin SV states are small, so the simple double loop is fine.

_sv_subset_probs(state, phys_sites) -> Vector{Float64}

Schmidt probabilities (RDM eigenvalues) of an arbitrary subset of physical sites of a dense state vector: permute the kept sites' tensor dimensions to the front, reshape to (d^m × d^(L-m)), and square the singular values. The ordering of kept dimensions does not affect the spectrum.

_validate_mutual_information(mi::MutualInformation, state) -> nothing

Shared evaluation-time validation for the MPS, state-vector, and Clifford backends: region sites within 1:L AND each region CONTIGUOUS (these backends keep the historical contiguous-region contract). The Gaussian override calls _mi_validate_bounds only and accepts arbitrary site subsets (see src/Gaussian/mutual_information.jl).

_validate_pauli_string(obs::PauliString, state) -> nothing

Shared evaluation-time validation for all backends: qubit-only (v0.4.0) and site indices within 1:L.

_von_neumann_entropy(mps::MPS, i::Int; n::Real=1, threshold::Float64=1e-16, base::Float64=2.0) -> Float64

Compute entanglement entropy at bond i of an MPS.

Arguments:

• mps: The MPS state
• i: The bond index (site index) where entropy is computed
• n: Rényi index (1=von Neumann, n=Rényi-n); any finite real > 0
• threshold: Minimum threshold for singular values to avoid log(0)
• base: Base of logarithm for entropy computation (default: 2.0 for bits)

Returns:

• Entanglement entropy value

The function:

1. Orthogonalizes the MPS to site i
2. Performs SVD on the tensor to extract Schmidt values
3. Computes probabilities from Schmidt values (squared)
4. Returns entropy based on Rényi index, in three branches of δ = n − 1:
   • |δ| <= _RENYI_SHUNT: von Neumann entropy $S_1 = -\sum p \log_b(p)$ — the exact n = 1 expression, used as the continuous limit (error O(1e-8))
   • _RENYI_SHUNT < |δ| <= _RENYI_NEAR1: cancellation-safe near-1 form
   • |δ| > _RENYI_NEAR1: Rényi entropy $S_n = \log_b\!\big(\sum p^n\big)/(1-n)$, evaluated in a scale-before-overflow log domain

Hartley entropy (n = 0) is NOT a case here: EntanglementEntropy rejects renyi_index <= 0 at construction (see _ee_check_common).

_warn_bricklayer_odd_pbc(geo::Bricklayer, L::Int, bc::Symbol)

Emit a one-time warning (once per (parity, L) combination per session, via maxlog=1 with a per-combination log _id) when a single brickwork layer under periodic boundary conditions sits on an odd-length ring and therefore is not a disjoint pair layer. Two cases:

• :odd / :even with odd L. An odd ring has no valid brickwork tiling: the layer leaves one site unpaired (:odd → site L, :even → site 1 — no wrap pair is added, see elements(::Bricklayer, ...)), and the wrap bond (L, 1) is gated by neither single layer, so an alternating :odd/:even circuit is effectively open across that bond.
• :nnn_odd_2 / :nnn_even_2 with L ≥ 6 and L % 4 == 2. The NNN bonds of an even ring form two rings of length L/2 (odd sites, even sites); when L/2 is odd neither ring splits into two disjoint pair layers, and the sublayer's wrap pair ((L-1, 1) resp. (L, 2)) shares its first site with the last bulk pair ((L-3, L-1) resp. (L-2, L)). The two gates on that site are applied sequentially in enumeration order, so the sublayer is not a depth-1 layer. Bond coverage is unaffected: :nnn still gates every NNN bond exactly once.

Enumeration behavior is NOT changed by this helper — it only warns.

Called from the circuit-builder recording path (apply!(builder, ...), apply_with_prob!(builder; ...) — fires at circuit-definition time) and the immediate-mode dispatch path (_apply_dispatch!(state, gate, ::Bricklayer)). Deliberately NOT called inside elements() itself, which sits in performance-critical loops (benchmarks, per-element expansion).

:nn and :nnn (ALL NN / NNN bonds — bond enumerations, not single layers) and the wrap-free sublayers :nnn_odd_1 / :nnn_even_1 are exempt. At odd L the NNN sublayers are left as they are (odd-L NNN coverage policy is a separate open question).

_write_op_block!(mps::MPS, i_list::Vector{Int}, mps_ij::ITensor, cutoff::Float64, maxdim::Int)

Second phase of apply_op_internal!: write the contracted block mps_ij back into mps over the sorted RAM range i_list — direct assignment for a single site, an SVD chain (respecting cutoff/maxdim) for a multi-site range.

On return the gauge center of mps is at i_list[end], and the MPS ortho limits say so. For a single site this is where _contract_op_block! left it. For a multi-site range the chain writes left-orthonormal U factors on i_list[1]:i_list[end]-1 and the remaining S * V on i_list[end]; ITensorMPS's setindex! only ever widens the ortho limits (to i_list[1]:i_list[end] here), so they are reset explicitly. Without this the next orthogonalize! would re-factorize the already-orthonormal Us, and normalize!/norm would fall back to a full-chain inner product.

_xlogx(x::Real) -> Float64

x * log(x) with the exact limit $0 \cdot \log(0) = 0$ (returns 0.0 for any x <= 0). Natural log. Used by subsystem_entropy so that exactly (un)occupied modes ($\lambda \in \{0, 1\}$) contribute zero entropy with no NaN/Inf and no imaginary-epsilon regularization.

advance!(geo::StaircaseLeft, L::Int, bc::Symbol)

Advance staircase left by one position. Internal use by apply!.

• StaircaseLeft: pos -= 1, wraps 1 → L (PBC) or 1 → L-range (OBC)

advance!(geo::StaircaseRight, L::Int, bc::Symbol)

Advance staircase right by one position. Internal use by apply!.

• StaircaseRight: pos += 1, wraps L → 1 (PBC) or L-range → 1 (OBC)

antisymmetrize!(Γ::AbstractMatrix{Float64}, rows) -> Γ

antisymmetrize! restricted to the pairs $(i, j)$ with $i$ or $j$ in rows. For a Γ that was exactly antisymmetric before an update touching only those rows and columns (a Majorana rotation $\Gamma \leftarrow R\Gamma R^T$ with $R = O \oplus I$), this is bit-for-bit the full pass: on the untouched block the full pass computes $(a - (-a))/2 = a$, an identity, so skipping it changes nothing. $O(N\,\lvert\mathrm{rows}\rvert)$ instead of $O(N^2)$.

antisymmetrize!(Γ::AbstractMatrix{Float64}) -> Γ

In-place $\Gamma \leftarrow (\Gamma - \Gamma^T)/2$: every off-diagonal pair becomes $(\Gamma_{ij} - \Gamma_{ji})/2$ and $(\Gamma_{ji} - \Gamma_{ij})/2$, the diagonal becomes zero. Bit-for-bit the same result as the broadcast Γ .= (Γ .- transpose(Γ)) ./ 2, without the temporary copy of transpose(Γ) that broadcast allocates to avoid aliasing. Both halves of a pair are computed from their own subtraction rather than as v and -v: the two agree for every nonzero value, but a pair of exact zeros gives +0.0/+0.0 under subtraction and +0.0/-0.0 under negation, and the sign of a zero changes LAPACK's Householder sign choices downstream (eigvals of $i\Gamma_A$ moved by 1e-14 before this was pinned). Returns Γ.

apply_feedback!(fb::AbstractFeedback, state, sites::Vector{Int}, outcome)

Internal: dispatch a feedback specification after a measurement. OnOutcome looks up the gate registered for outcome (no-op when absent) and executes it on sites; CallbackFeedback calls the wrapped function.

apply_gate_sv!(ψ, U, target_sites, L, d) -> Vector{ComplexF64}

Apply a d^n × d^n gate matrix U to target_sites (physical site numbers, 1-indexed, site 1 = MSB) of a length-d^L state vector ψ. Returns the new state vector (caller reassigns state.backend.ψ = result, does NOT mutate ψ's contents in-place — matches the existing state.backend.mps = ... reassignment pattern used elsewhere in the codebase).

Algorithm (verified numerically across multiple test cases, including an asymmetric CNOT gate applied in both site-argument orders):

1. Reshape ψ into an L-dimensional tensor of shape (d,d,...,d). Due to Julia's column-major array layout, tensor dimension k corresponds to PHYSICAL SITE L - k + 1 (site 1/MSB is the LAST dimension, site L/LSB is the FIRST dimension).
2. Convert each target physical site s to its Julia array dimension via L - s + 1, then REVERSE the resulting list (required so the merged multi-site index matches gate_matrix's "first site in the gate's argument list = slowest/most-significant digit" convention).
3. Permute dimensions so the (reversed) target dims come first, followed by all other dims (in ascending order).
4. Reshape the permuted tensor to 2D: (d^n, d^(L-n)) where n = length(target_sites).
5. Left-multiply by U: new_mat = U * A_mat.
6. Reshape back to L dims (d,d,...,d), then permutedims with the INVERSE of the original permutation.
7. vec(...) the result back to a flat Vector{ComplexF64}.

Handles 1-site AND n-site gates (adjacent or non-adjacent) UNIFORMLY — no special-casing needed for adjacent vs non-adjacent target sites.

apply_gate_sv_optimized!(ψ, U, target_sites, L, d) -> ψ

Tier 2 entry point: dispatches to the zero-allocation stride-loop kernels for 1-site and 2-site gates, or the reshape/permutedims fallback for n-site (n>=3) gates. Mutates ψ in place and returns it (caller reassigns state.backend.ψ = apply_gate_sv_optimized!(...), matching Tier 1's reassignment call-site pattern even though Tier 2 itself mutates in place).

apply_gate_sv_optimized_1site!(ψ, U, site, L, d) -> ψ

In-place, zero-inner-loop-allocation 1-site gate application via the bit/digit-stride pattern (Yao.jl's u1apply!, generalized from d=2 to arbitrary local dimension d).

step = d^(L - site) is the stride between the d amplitudes that get mixed together by U at each combined "outer" index — this matches Tier 1's convention that physical site s corresponds to Julia tensor dimension L - s + 1 (site 1 = MSB = slowest-varying digit).

Mutates ψ in place and also returns it (for chaining / consistency with the rest of the Tier 2 API).

apply_gate_sv_optimized_2site!(ψ, U, s1, s2, L, d) -> ψ

In-place 2-site gate application (adjacent or not) by the same digit-stride pattern as the 1-site kernel. The d² amplitudes mixed by U at a base index n0 (both target digits zero) sit at n0 + i*st1 + j*st2, where st = d^(L - s) is the stride of each target site; the merged gate index is i*d + j with i the digit of s1, the FIRST site in the gate's argument list — the "first site = most-significant digit" convention of gate_matrix and Tier 1's apply_gate_sv!. The three nested loops enumerate exactly the base indices, so there is no per-amplitude digit test and no copy of ψ.

apply_gate_sv_optimized_nsite!(ψ, U, target_sites, L, d) -> ψ

n-site (n >= 3) fallback for the optimized engine: adjacent AND non-adjacent target sites, using the SAME reshape + permutedims + matmul algorithm as Tier 1's apply_gate_sv! (guaranteeing bitwise-identical results by construction), but mutating ψ IN PLACE via ψ .= vec(new_A) rather than allocating and returning a fresh vector. No built-in gate has support ≥ 3 (only a user MatrixGate can), so this path is off the hot loop; the 1- and 2-site stride kernels above cover every gate the circuits ship with.

apply_op_internal!(mps::MPS, op::ITensor, sites::Vector{Index}, cutoff::Float64, maxdim::Int)

Apply operator to MPS following CT.jl algorithm (lines 147-172): _contract_op_block! (gauge + contract + apply) followed by _write_op_block! (SVD chain write-back). _apply_single! runs the two phases itself for projective gates, to reject an impossible postselection between them.

Contract 3.6: Index matching via Index comparison, NOT tag parsing.

build_operator(gate, site_or_sites, local_dim; kwargs...) -> ITensor

Build the operator tensor for a gate acting on given site(s).

build_operator(gate::CNOT, sites::Vector{Index}, local_dim::Int) -> ITensor

Build CNOT gate operator. Control = sites[1], target = sites[2]. Qubit-only (local_dim == 2).

Vectorized construction: reshape the fixed gate_matrix(::CNOT) into a 4-index tensor, same convention as CZ/MatrixGate above. Bitwise identical to the loop-based construction it replaced (element-by-element ==).

build_operator(gate::CZ, sites::Vector{Index}, local_dim::Int) -> ITensor

Build CZ gate operator. Qubit-only (local_dim == 2).

Vectorized construction: reshape the fixed gate_matrix(::CZ) into a 4-index tensor and hand it directly to ITensor(...), following the same output-primed-first, input-unprimed-second, reverse-site-order convention as MatrixGate (see matrix_gate.jl:98-110) — eliminates the 16 scalar setindex! calls of the previous per-application loop. Bitwise identical to the loop-based construction it replaced (element-by-element ==).

build_operator(gate::HaarRandom, site::Index, local_dim::Int; rng) -> ITensor

Single-site (n = 1) Haar random unitary. Same Ginibre-QR algorithm on a d × d matrix, consuming from :gates_realization.

build_operator(gate::HaarRandom, sites::Vector{Index}, local_dim::Int; rng::RNGRegistry) -> ITensor

Build an n-site Haar random unitary operator from the :gates_realization stream. For gate.n == 2 this is bit-identical to the historical implementation (same randn call shape/order, same index ordering).

build_operator(gate::MatrixGate, site::Index, local_dim::Int) -> ITensor

Single-site MatrixGate: ITensor(U, site', site) with primed = output.

build_operator(gate::MatrixGate, sites::Vector{Index}, local_dim::Int) -> ITensor

Reshape the stored d^n × d^n matrix into a 2n-index ITensor. Column-major reshape puts the LAST site of the region on the fastest matrix digit, so the tensor dims map to $(s_n', \ldots, s_1', s_n, \ldots, s_1)$ with primed = output (row) and unprimed = input (column).

build_operator(gate::PauliX, site::Index, local_dim::Int) -> ITensor

Build Pauli X operator tensor.

build_operator(gate::PauliY, site::Index, local_dim::Int) -> ITensor

Build Pauli Y operator tensor.

build_operator(gate::PauliZ, site::Index, local_dim::Int) -> ITensor

Build Pauli Z operator tensor.

build_operator(gate::PhaseGate, site::Index, local_dim::Int) -> ITensor

Build Phase (S) operator tensor.

build_operator(gate::Projection, site::Index, local_dim::Int) -> ITensor

Build projection operator $\lvert \text{outcome} \rangle\langle \text{outcome} \rvert$ via the site type's per-level "Proj<k>" op string (defined for Qubit/"S=1/2" by ITensors and for all spin site types by src/Core/spin_sites.jl).

build_operator(gate::RandomClifford, site::Index, local_dim::Int; rng) -> ITensor

Single-site (n = 1) random Clifford unitary. Same conversion path on a d × d matrix, consuming from :gates_realization.

build_operator(gate::RandomClifford, sites::Vector{Index}, local_dim::Int; rng) -> ITensor

Build an n-site random Clifford unitary operator from the :gates_realization stream. Follows the same index-ordering convention as HaarRandom's build_operator (see two_qubit.jl): output-primed-first, input-unprimed-second, reverse-site-order.

build_operator(gate::SWAP, sites::Vector{Index}, local_dim::Int) -> ITensor

Build SWAP gate operator. Qubit-only (local_dim == 2).

Vectorized construction: reshape the fixed gate_matrix(::SWAP) into a 4-index tensor, same convention as CZ/CNOT/MatrixGate above. Bitwise identical to the loop-based construction it replaced (element-by-element ==).

build_operator(gate::SpinSectorMeasurement, sites::Vector{Index}, local_dim::Int; rng, mps, ram_sites) -> ITensor

Build measurement operator for two spin-1 sites with proper Born sampling.

Implements the Born rule measurement:

1. Compute Born probabilities $p(S) = \langle \psi \rvert P_S \lvert \psi \rangle$ for each allowed sector S
2. Normalize probabilities over allowed sectors
3. Sample outcome S with probability p(S) using the provided RNG
4. Return the projector P_S onto the sampled sector

Arguments

• gate: SpinSectorMeasurement specifying allowed sectors
• sites: ITensor site indices for the two sites
• local_dim: Local Hilbert space dimension (d = 2s+1; 3 for spin-1)
• rng: RNG registry for reproducible sampling (uses :born_measurement stream)
• mps: Current MPS state for computing Born probabilities
• ram_sites: RAM indices of the two sites

Returns

• ITensor projector onto the randomly sampled spin sector

build_operator(gate::SpinSectorProjection, sites::Vector{Index}, local_dim::Int; kwargs...) -> ITensor

Build projector operator for two spin-1 sites. Returns ITensor representation of the projector matrix.

build_operator(gate::Union{Rx,Ry,Rz,Hadamard}, site::Index, local_dim::Int) -> ITensor

Build the exact 2×2 matrix operator: ITensor(M, site', site) with primed = output. Qubit-only (local_dim == 2).

build_template_groups_ascii(circuit; n_steps::Int=1) -> Vector{Vector{Vector{ExpandedOp}}}

Build operation groups for ASCII rendering by iterating circuit.operations directly, showing ALL outcomes of stochastic operations unconditionally (circuit template).

For deterministic ops: one group per apply! call (same as expand_circuit_grouped). For stochastic ops: one group per outcome (ALL outcomes shown, no random selection).

Returns the same steps → groups → ops structure as expand_circuit_grouped.

Compute physical-to-RAM and RAM-to-physical site mappings.

For OBC: identity mapping (1:L -> 1:L) For PBC: folded mapping - interleaves sites zig-zagging outward from pbc_fold_start

Parameters:

• L: system size
• bc: boundary condition (:open or :periodic)
• pbc_fold_start: physical site the PBC zig-zag fold starts from (default L÷4+1, the middle-aligned choice giving a contiguous half-cut). Ignored for bc == :open. Must satisfy 1 <= pbc_fold_start <= L for bc == :periodic.

Returns: (phy_ram, ram_phy) where:

• phy_ram[physical_site] = ram_index
• ram_phy[ram_index] = physical_site

compute_pair_staircase(pos::Int, L::Int, bc::Symbol) -> Vector{Int}

Convert a single staircase position to an adjacent pair [pos, pos+1].

Arguments

• pos: Position (1-indexed)
• L: System size (number of sites)
• bc: Boundary condition (:periodic or :open)

Returns

• Two-element vector [pos, second] where second wraps at L→1 for PBC

Validation

• For OBC, pos must be in range 1:(L-1) (cannot form pair at position L)
• For PBC, pos can be 1:L (position L pairs with 1)

Example

compute_pair_staircase(3, 5, :periodic)  # Returns [3, 4]
compute_pair_staircase(5, 5, :periodic)  # Returns [5, 1] (wrap)
compute_pair_staircase(4, 5, :open)      # Returns [4, 5]
compute_pair_staircase(5, 5, :open)      # Throws error (no site 6)

compute_projector_product_expectation(state, sites_zero::Vector{Int}, site_one::Int) -> Float64

Compute $\langle\psi| \big(\textstyle\prod_k P_{0,k}\big) P_1 |\psi\rangle$ where $P_0$ projects to $|0\rangle$ and $P_1$ projects to $|1\rangle$. sites_zero: physical sites that should be "0" site_one: physical site that should be "1"

This uses MPO construction for the projector product. The returned matrix element is NOT divided by $\langle\psi|\psi\rangle$; domain_wall computes the norm once and divides every projector expectation by it.

compute_site_staircase_left(start::Int, step::Int, L::Int, bc::Symbol) -> Int

Compute the position of a StaircaseLeft geometry after step iterations.

Arguments

• start: Starting position (1-indexed)
• step: Step number (1 = initial position, 2 = after first advance, etc.)
• L: System size (number of sites)
• bc: Boundary condition (:periodic or :open)

Returns

• Position after (step-1) advances

Wrapping behavior

• PBC: Position cycles L → ... → 2 → 1 → L
• OBC: Position cycles (L-1) → ... → 2 → 1 → (L-1)

Example

# Start at position 3, system size 5, periodic BC
compute_site_staircase_left(3, 1, 5, :periodic)  # Returns 3 (initial)
compute_site_staircase_left(3, 2, 5, :periodic)  # Returns 2 (after 1 advance)
compute_site_staircase_left(3, 3, 5, :periodic)  # Returns 1 (after 2 advances)
compute_site_staircase_left(3, 4, 5, :periodic)  # Returns 5 (after 3 advances, wrapped)

compute_site_staircase_right(start::Int, step::Int, L::Int, bc::Symbol) -> Int

Compute the position of a StaircaseRight geometry after step iterations.

Arguments

• start: Starting position (1-indexed)
• step: Step number (1 = initial position, 2 = after first advance, etc.)
• L: System size (number of sites)
• bc: Boundary condition (:periodic or :open)

Returns

• Position after (step-1) advances

Wrapping behavior

• PBC: Position cycles 1 → 2 → ... → L → 1
• OBC: Position cycles 1 → 2 → ... → (L-1) → 1

Example

# Start at position 3, system size 5, periodic BC
compute_site_staircase_right(3, 1, 5, :periodic)  # Returns 3 (initial)
compute_site_staircase_right(3, 2, 5, :periodic)  # Returns 4 (after 1 advance)
compute_site_staircase_right(3, 3, 5, :periodic)  # Returns 5 (after 2 advances)
compute_site_staircase_right(3, 4, 5, :periodic)  # Returns 1 (after 3 advances, wrapped)

compute_sites(geo::AdjacentPair, step::Int, L::Int, bc::Symbol) -> Vector{Int}

Compute sites for AdjacentPair geometry.

Arguments

• geo: AdjacentPair geometry
• step: Step number (unused for static geometry)
• L: System size
• bc: Boundary condition (:periodic or :open)

Returns

• Two-element vector [geo.first, second] where second wraps at L→1 for PBC

compute_sites(geo::AbstractStaircase, step::Int, L::Int, bc::Symbol, gate::AbstractGate) -> Vector{Int}

Compute sites for a StaircaseRight/StaircaseLeft geometry based on gate support, at the staircase's CURRENT position (step is ignored; advance! manages progression).

Arguments

• geo: StaircaseRight or StaircaseLeft geometry
• step: Step number (unused)
• L: System size
• bc: Boundary condition (:periodic or :open)
• gate: Gate to determine support (1-site or 2-site)

Returns

• For single-site gates (support == 1): [pos]
• For two-site gates (support == 2): [pos, pos+range] — the same region elements reports (PBC wrap via mod1; OBC out-of-bounds throws ArgumentError)

compute_sites(geo::SingleSite, step::Int, L::Int, bc::Symbol) -> Vector{Int}

Compute sites for SingleSite geometry (always returns [geo.site]).

Arguments

• geo: SingleSite geometry
• step: Step number (unused for static geometry)
• L: System size (unused for static geometry)
• bc: Boundary condition (unused for static geometry)

Returns

• Single-element vector [geo.site]

compute_sites_dispatch(geo::AbstractGeometry, gate::AbstractGate, step::Int, L::Int, bc::Symbol) -> Vector{Int}

Dispatch compute_sites with appropriate arguments based on geometry type.

For StaircaseRight/StaircaseLeft: requires gate parameter to determine support. For SingleSite/AdjacentPair/Sites: gate parameter not needed.

compute_two_site_born_probability(mps::MPS, projector::Matrix{Float64}, ram_sites::Vector{Int}, local_dim::Int) -> Float64

Compute the Born probability $\langle \psi \rvert P \lvert \psi \rangle / \langle \psi \vert \psi \rangle$ for a two-site projector.

Computes it by:

1. Orthogonalizing MPS to the target region
2. Contracting the two-site block into a local tensor
3. Computing $\langle T \rvert P \lvert T \rangle$ using proper index contraction
4. Dividing by $\langle T \vert T \rangle = \langle \psi \vert \psi \rangle$ (the norm of the gauge-center tensor), so the result is the physical probability even when truncation has left the MPS un-normalized — the MPS backend does not renormalize after unitary layers.

Arguments

• mps: Current MPS state
• projector: 9×9 projector matrix (d²×d² where d=local_dim)
• ram_sites: RAM indices of the two sites [i, j]
• local_dim: Local Hilbert space dimension (3 for spin-1)

Returns

• Probability $p = \langle \psi \rvert P \lvert \psi \rangle / \langle \psi \vert \psi \rangle$ (real, non-negative; NaN on a zero-norm MPS, for which no probability is defined)

conjugate_rows_cols!(Γ, O, ix, scratch) -> Γ

In-place $\Gamma \leftarrow R\,\Gamma\,R^T$ with $R = O \oplus I$ acting on the Majorana indices ix (O is m×m, m = length(ix)): Γ[ix, :] = O * Γ[ix, :] followed by Γ[:, ix] = Γ[:, ix] * O', which is the direct (exact) action of a Gaussian unitary on the covariance matrix. The row and column blocks are staged in scratch (the $N\times N$ state.backend.scratch) so that no temporaries are allocated; the two products are the same BLAS gemm calls as the slice form and give the same floating-point result. When the buffer cannot hold the $2m \times N$ / $N \times 2m$ staging blocks ($N < 2m$: fewer than four fermionic modes, or fewer than four Majoranas on a Majorana chain), the allocating slice form is used instead. Only rows and columns ix are modified; the result is not re-antisymmetrized here (see antisymmetrize!(Γ, ix)).

current_position(p::Pointer) -> Int

Get the current position (read-only).

current_position(geo::AbstractStaircase) -> Int

Get the current position of the staircase (READ-ONLY accessor).

domain_wall(state::SimulationState, i1::Int, order::Int) -> Float64

Compute domain wall observable at sampling site i1 with given order. Ports CT.jl's dw_FM algorithm for xj=Set([0]) case.

The domain wall measures where the first "1" appears in the bit string when scanning cyclically from position i1.

Each probability is $\langle\psi| \big(\textstyle\prod_{k<j} P_{0,k}\big) P_{1,j} |\psi\rangle / \langle\psi|\psi\rangle$: the MPS is not renormalized after truncating unitary layers, so the norm (computed once per call, see _mps_norm2) is divided out to keep the weights true probabilities of the retained state.

domain_wall(state::SimulationState{StateVectorBackend}, i1::Int, order::Int) -> Float64

State-vector implementation of the domain_wall function (mirrors the formula of the existing MPS domain_wall in src/Observables/domain_wall.jl exactly: cyclic scan from i1, weight (L-j+1)^order, term j is $\langle\psi\rvert\left(\prod_{k<j}P0_k\right)P1_j\lvert\psi\rangle$ over the cyclic site list).

Implementation: a single O(d^L) pass instead of L independent O(d^L) scans (one per term j). The L projector-product conditions are mutually exclusive in the computational basis: basis state n satisfies the term-j condition iff j is the FIRST position in the cyclic site list whose digit is nonzero AND that digit equals 1. So one sweep over the amplitudes finds, for each basis state, its unique contributing term (if any) by scanning for the first nonzero digit (amortized O(1) per basis state — a depth-k scan occurs for a d^-k fraction of basis indices). Each per-term accumulator receives exactly the same addends in the same ascending basis-index order as the former per-term scans, and the final weighted sum runs in the same j order, so the result is bitwise identical to the O(L·d^L) implementation.

draw(registry::RNGRegistry, stream::Symbol) -> Float64
draw(state::SimulationState, stream::Symbol) -> Float64

Draw ONE scalar uniform Float64 from the named RNG stream. This is the canonical way to consume a stream (replaces the removed, type-pirating rand(state, stream) extension; behavior is identical).

All engine coin draws are scalar — see the SCALAR-DRAW CONTRACT note in the RNGRegistry docstring.

Example:

if draw(state, :gates_spacetime) < p_ctrl
    # control branch
end

draw(state::SimulationState, stream::Symbol) -> Float64

Draw ONE scalar uniform Float64 from the named stream of the state's RNG registry (see draw(::RNGRegistry, ::Symbol)). Errors if no registry is attached.

execute!(state::SimulationState, gate::AbstractGate, region::Vector{Int})

Execute gate on the physical sites region (already resolved from a geometry). This is the v0.1 gate-execution protocol:

• Default implementation: build_operator → apply_op_internal!, then normalize (and, with retruncate == :gate, a full-chain truncate!) iff needs_normalization(gate) (see _apply_single!).
• Gate-specific overrides: gates with non-operator semantics (Born sampling + classical logic) override this method — see Measure and Reset below. User-defined gates may override it the same way:

function QuantumCircuitsMPS.execute!(state::SimulationState, g::MyGate, region::Vector{Int})
    # custom behavior; region is a vector of physical sites
end

Related traits: needs_normalization(gate) (post-apply renormalization) and is_measurement(gate) (gate Born-samples via :born_measurement).

region must not repeat a site: every backend's _apply_single! assumes distinct sites (a repeated site gives a non-unitary tableau update on the Clifford backend, for instance), so a repeated site is an ArgumentError here, before any backend is touched. This covers apply!(state, gate, sites::Vector{Int}), every geometry, and the circuit engine.

execute!(state::SimulationState, gate::Measure, region::Vector{Int})

Measure (v0.1, feedback-capable): Born-sample the single site in region (via _measure_single_site!, which emits the MeasurementOutcome event with the engine's real step/op indices), then dispatch the gate's feedback — if any — with the observed outcome.

Feedback runs INSIDE with_guarded_stream(registry, :gates_spacetime): it can never consume spacetime coins (fixed-draw contract), while :gates_realization and :born_measurement remain freely usable. Feedback gates execute through this same execute! protocol but are NOT engine op slots: they advance no counters and emit no GateApplied events. Recursion (feedback applying another Measure) is allowed — user responsibility.

execute!(state::SimulationState, pg::ProductGate, region::Vector{Int})

Execute a ProductGate: validate that region covers exactly the union of elements(pg.region_geometry, state.L, state.bc), then apply pg.inner to each element in canonical order via the uniform execute! protocol (one fresh :gates_realization draw per element for random inner gates). Logs one GateApplied per element (inner gate's label) when the event log is enabled.

execute!(state::SimulationState, gate::Reset, region::Vector{Int})

Reset (DERIVED - measurement + conditional X): Born-sample the single site in region; if the outcome is 1, flip it back to $|0\rangle$ with PauliX.

Defined for two-level sites only (local_dim == 2: "Qubit", "S=1/2", a local_dim=2 "Qudit"). On any other site type the conditional X cannot return the site to $|0\rangle$ (a spin-1 site measured in level 2 would be left there), so the gate is rejected with an ArgumentError before the measurement — the state and the :born_measurement stream are untouched.

execute!(state::SimulationState{GaussianBackend}, gate::BondParity, phy_sites::Vector{Int})

Projective measurement of the BOND fermion parity i γ_a γ_b on the two INNER Majoranas of an adjacent-site bond, for the fermionic Gaussian backend.

Inner-Majorana convention: site i carries Majoranas (2i−1, 2i) (RAM-mapped). For the bond between sites (i, i+1) the measured pair is ix = [2i, 2i+1] — site i's SECOND Majorana and site i+1's FIRST Majorana. For the periodic wrap bond (L, 1) (only valid when state.bc == :periodic) the pair is ix = [2L, 1] — site L's second Majorana and site 1's first. Sites may be passed in either order ([i, i+1] or [i+1, i]); non-adjacent pairs throw ArgumentError, as does the wrap bond under bc = :open.

Majorana-chain granularity (site_type="Majorana", majoranas_per_site == 1): each site IS one Majorana, so the bond (i, i+1) measures the pair ix = [i, i+1] directly — this is the i γ_i γ_{i+1} parity measurement of the class-DIII monitored Majorana chain (wrap bond (L, 1) → ix = [L, 1] under PBC). Same Born rule, same kernel, same draw contract — only the site→Majorana index mapping (via site_majoranas) differs.

Born rule (cross-validated against the exact-diagonalization/Pfaffian reference in test/gaussian/oracle.jl): the covariance element $g = \Gamma[ix_1,ix_2]$ is the parity expectation value, $\langle i\hat\gamma_{ix_1}\hat\gamma_{ix_2}\rangle = g$, and the outcome encoding is $\mathrm{outcome} = 0$ ↔ parity eigenvalue $+1$, $\mathrm{outcome} = 1$ ↔ parity eigenvalue $-1$ (parity eigenvalue $= 1 - 2\cdot\mathrm{outcome}$), so:

\[P(\mathrm{outcome}=0)=\frac{1+g}{2},\qquad P(\mathrm{outcome}=1)=\frac{1-g}{2}\]

— the SAME affine structure as the on-site occupation measurement (_measure_single_site!), because on-site occupation is itself the bond parity of the intra-mode pair ($n = (1 - i\gamma_{2i-1}\gamma_{2i})/2$). The kernel sign $s = 2\cdot\mathrm{outcome}-1$ is MINUS the measured parity eigenvalue: parity_projection!(Γ, ix₁, ix₂, s) leaves $\Gamma[ix_1,ix_2] = -s$, i.e. $\langle i\hat\gamma\hat\gamma\rangle = -s = 1 - 2\cdot\mathrm{outcome}$ post-measurement.

Draw contract (REDUNDANT-DRAW, cross-backend lockstep): exactly ONE scalar :born_measurement draw is consumed per BondParity application — always, unconditionally, BEFORE any probability computation — matching _measure_single_site! draw-for-draw.

Throws ArgumentError if the sampled outcome has probability ≤ 1e-15. Logs a MeasurementOutcome event exactly like Measure's execute! (via _measure_single_site!) does. Returns nothing.

execute!(state::SimulationState{GaussianBackend}, gate::Reset, region::Vector{Int})

Always throws an ArgumentError: Reset is not available on the Gaussian backend, on either site granularity.

Reset is "measure, then PauliX if the outcome is 1", and the flip is the problem. A covariance matrix can only flip one occupation by reflecting a single Majorana (conjugation by $\gamma_{2i}$), which is a parity-odd operator: it changes the fermion parity of the state, something no closed fermionic system can do. Rather than ship a flip that is physically meaningless for fermions, both PauliX and Reset are rejected.

The generic execute!(::SimulationState, ::Reset, ...) in src/Core/apply.jl would Born-sample the site FIRST and only then fail on PauliX (for outcome 1 — and silently succeed for outcome 0). This override exists so the rejection happens before anything changes: the covariance matrix and the :born_measurement stream are untouched when it throws.

To project a mode onto a definite occupation use Measure(:Z) (random outcome); to prepare an occupation pattern use initialize! with ProductState(bitstring=...). For the qubit Reset, use backend=:mps or backend=:statevector.

gate_label(gate::AbstractGate) -> String

Return a short visualization label for a gate type.

Labels

• Reset → "Rst"
• HaarRandom → "Haar"
• RandomClifford → "Cl"
• Projection → "Prj"
• PauliX → "X"
• PauliY → "Y"
• PauliZ → "Z"
• CZ → "CZ"
• Other → Type name as string

Examples

gate_label(Reset())       # Returns "Rst"
gate_label(HaarRandom())  # Returns "Haar"
gate_label(CZ())          # Returns "CZ"

gate_matrix(gate::AbstractGate) -> Matrix{ComplexF64}

Return the explicit matrix of a gate whose action is defined by a fixed matrix (e.g. MatrixGate, Rx, Ry, Rz, Hadamard). Matrix convention: $U_{\text{out},\text{in}} = \langle \text{out} \rvert U \lvert \text{in} \rangle$ with Kronecker (row-major site) ordering — the FIRST site of the region is the slowest basis digit, so kron(A, B) acts with A on the first site.

gate_matrix(g::HaarRandom, rng::AbstractRNG; local_dim::Int=2) -> Matrix{ComplexF64}

State-vector-path equivalent of build_operator(gate::HaarRandom, ...): draws a fresh d^n × d^n Haar random unitary (d = local_dim, n = g.n) by reusing the same _haar_unitary core used by the MPS build_operator path. Consumes from whichever RNG stream rng is (caller is responsible for passing the appropriate stream, e.g. :gates_realization).

gate_matrix(g::RandomClifford, rng::AbstractRNG; local_dim::Int=2) -> Matrix{ComplexF64}

State-vector-path equivalent of build_operator(gate::RandomClifford, ...): draws a fresh d^n × d^n Haar-random Clifford unitary (d = local_dim, n = g.n) by reusing the same _random_clifford_unitary core used by the MPS build_operator path. Consumes from whichever RNG stream rng is (caller is responsible for passing the appropriate stream, e.g. :gates_realization).

gate_matrix(gate::SpinSectorProjection) -> Matrix{ComplexF64}

State-vector-path equivalent of build_operator(gate::SpinSectorProjection, ...): the 9×9 projector as a dense complex matrix (this method was previously missing, so apply! on backend=:statevector threw MethodError and the README AKLT protocol could not run on the SV backend).

Basis ordering matches the state-vector engine convention ("first site in the gate's argument list = slowest/most-significant digit", see apply_gate_sv!): projector row/col = (m1-1)*3 + m2 with m1 = first site (slow), m2 = second site (fast). Renormalization after application is provided by the needs_normalization trait above (honored by the SV _apply_single! path).

gaussian_contraction!(Γ, Υ, ix; scratch=nothing, purify_tol=1e-10) -> Γ

Contract the Gaussian gate tensor Υ (2m×2m for m = length(ix) Majorana legs) into the covariance matrix Γ at Majorana indices ix (1-based), in place — line-by-line port of P_contraction_2 (GTN.py:1063-1102):

\[C = (\Gamma_{RR}\Upsilon_{LL}+I)^{-1}\\ \Gamma[\bar{ix},\bar{ix}] \mathrel{+}= \Gamma_{LR}(\Upsilon_{LL}C)\Gamma_{LR}^T\\ \Gamma[ix,\bar{ix}] = \Upsilon_{RL}C\Gamma_{LR}^T\\ \Gamma[ix,ix] = \Upsilon_{RR}+\Upsilon_{RL}(\Gamma_{RR}C^T)\Upsilon_{RL}^T\\ \Gamma[\bar{ix},ix] = -\Gamma[ix,\bar{ix}]^T\]

Throws ArgumentError when $\Gamma_{RR}\Upsilon_{LL}+I$ is singular (a probability-zero measurement outcome) — no least-squares fallback. If the result drifts off the pure-state manifold beyond purify_tol ($\max\lvert\mathrm{diag}(\Gamma^2)+1\rvert$, see _purity_drift), calls purify! and re-antisymmetrizes. scratch (same size as Γ) is an optional preallocated buffer for the ix̄-block update.

get_all_sites(geo::AllSites, state) -> Vector{Int}

Get all physical sites (1:L).

get_compound_elements(geo::AbstractGeometry, L::Int, bc::Symbol) -> Vector{Vector{Int}}

Get elements for compound geometry iteration.

Legacy name — delegates to the canonical elements; enumeration order is identical (bit-for-bit API contract).

Examples

# Bricklayer with odd parity on L=4 system
geo = Bricklayer(:odd)
get_compound_elements(geo, 4, :open)  # [[1, 2], [3, 4]]

# AllSites on L=3 system
geo = AllSites()
get_compound_elements(geo, 3, :open)  # [[1], [2], [3]]

get_op_ram_sites(op::ITensor, sites::Vector{Index}) -> Vector{Int}

Get RAM site indices from operator indices using Index comparison (Contract 3.6). Does NOT parse tags.

get_pairs(geo::Bricklayer, state) -> Vector{Tuple{Int,Int}}

Get all pairs for bricklayer pattern. Returns pairs of physical sites.

Legacy name — delegates to the canonical elements(geo, L, bc) (see Geometry/elements.jl); enumeration order is identical (bit-for-bit API contract).

get_sites(geo::AbstractGeometry, state) -> Vector{Int}

Get the physical sites for this geometry. Returns physical site indices. For iterating geometries (Bricklayer, AllSites), returns the first/next set. For staircases, returns current position pair.

get_sites(p::Pointer, state) -> Vector{Int}

Get the current pair of sites [pos, pos+1] for two-qubit gates. Handles PBC wrap when pos == L.

get_sites(geo::AbstractStaircase, state) -> Vector{Int}

Get current pair of physical sites for the staircase. Returns [pos, pos+range] with proper boundary condition handling.

haar_orthogonal(rng::AbstractRNG, n::Int) -> Matrix{Float64}

Exact Haar-random special-orthogonal matrix $Q \in SO(n)$: QR-decompose a Ginibre matrix, fix the R-diagonal sign ambiguity ($Q \leftarrow Q\,\mathrm{diag}(\mathrm{sign}(\mathrm{diag}(R)))$, required for Haar on O(n)), then flip the first column if det(Q) < 0 to land in SO(n).

is_aliased(registry::RNGRegistry) -> Bool

true when the :gates_spacetime and :gates_realization streams are the SAME RNG object (===), as in RNGRegistry(Val(:ct_compat); ...). Aliased registries are exempt from the fixed-draw invariant and from the with_guarded_stream sentinel guard (guarding one alias would also block the other; interleaved consumption is CT.jl-faithful by design).

is_compound_geometry(geo::AbstractGeometry) -> Bool

Check if geometry requires element-by-element iteration.

Compound geometries (Bricklayer, AllSites, EachSite) need to be expanded into multiple individual gate applications, one for each element/site pair.

Equivalent to the v0.1 trait is_broadcast; kept as a legacy alias until engine call sites are rewired.

is_measurement(gate::AbstractGate) -> Bool

Trait: does this gate perform a Born-rule measurement (consuming the :born_measurement RNG stream and collapsing the state)? Default false.

true for Measurement and SpinSectorMeasurement. Note Reset is a derived operation (measure + conditional flip) implemented via its own execute! method; it reports false here because the gate as a whole is a deterministic-outcome channel, not an outcome-reporting measurement.

is_record_mark(op::ExpandedOp) -> Bool

True when op is a record! marker pseudo-op (gate-less annotation produced from a (type=:record_mark, ...) circuit operation). Marker ops carry no gate (gate === nothing) and no sites; visualization renders them as a marker row (▽ in Unicode mode, [R] in ASCII mode).

log_event!(state::SimulationState, ev::CircuitEvent)

Append a CircuitEvent to the state's event log. Silent no-op when the log is disabled (the default), so emission sites can call this unconditionally. Returns nothing.

majorana_indices(site::Int) -> (Int, Int)

Majorana indices (2site-1, 2site) for 1-based fermionic mode site.

needs_normalization(gate::AbstractGate) -> Bool

Trait: does applying this gate require renormalizing (and truncating) the MPS afterwards? Default false (unitaries preserve the norm).

true for projective/collapsing gates (Projection, SpinSectorProjection, SpinSectorMeasurement). User-defined non-unitary gates opt in with:

QuantumCircuitsMPS.needs_normalization(::MyProjectiveGate) = true

Renormalization is only possible when the gate leaves something to renormalize: if the applied operator annihilates the state (a postselection onto an outcome of zero Born probability), the backends throw an ArgumentError and leave the state unchanged — see POSTSELECTION_PROB_TOL.

occupation_covariance(bits::AbstractVector{Bool}) -> Matrix{Float64}

Product-state covariance matrix for occupations bits: mode i gets block [[0,1],[−1,0]] when bits[i] == false ($\langle c^\dagger c\rangle = 0$) and the sign-flipped block $\begin{pmatrix}0&-1\\1&0\end{pmatrix}$ when bits[i] == true ($\langle c^\dagger c\rangle = 1$). Sign ↔ occupation mapping matches the Python reference get_C_f (GTN.py:1518-1523).

pack_ops_into_layers(ops::Vector{ExpandedOp}) -> Vector{Vector{ExpandedOp}}

Greedily pack a flat list of expanded operations into non-overlapping layers.

Each layer contains operations that act on disjoint sets of qubits. Operations are placed using a first-fit strategy: each op is assigned to the first existing layer where it does not conflict with any already-placed op, or a new layer is opened if no such layer exists.

Conflict detection uses set intersection on op.sites, so wrapped pairs like [8, 1] are handled correctly.

Arguments

• ops: Flat list of ExpandedOp values to pack (e.g., from expand_circuit).

Returns

A Vector{Vector{ExpandedOp}} where each inner vector is one non-overlapping layer. Returns an empty vector when ops is empty.

Examples

# Bricklayer(:nn) on L=8 periodic → 2 layers of 4 ops each
ops = [ExpandedOp(1, HaarRandom(), [1,2], "Haar"),
       ExpandedOp(1, HaarRandom(), [3,4], "Haar"),
       ExpandedOp(1, HaarRandom(), [5,6], "Haar"),
       ExpandedOp(1, HaarRandom(), [7,8], "Haar"),
       ExpandedOp(1, HaarRandom(), [2,3], "Haar"),
       ExpandedOp(1, HaarRandom(), [4,5], "Haar"),
       ExpandedOp(1, HaarRandom(), [6,7], "Haar"),
       ExpandedOp(1, HaarRandom(), [8,1], "Haar")]
layers = pack_ops_into_layers(ops)
length(layers)       # 2
length(layers[1])    # 4
length(layers[2])    # 4

parity_projection!(Γ, a, b, s; purify_tol=1e-10) -> Γ

Collapse the covariance matrix Γ onto the parity sector $i\gamma_a\gamma_b = -s$ ($s \in \{+1,-1\}$) of the Majorana pair (a, b), in place. This is gaussian_contraction! with Υ = parity_projection_upsilon(s) evaluated in closed form: for that projector $\Upsilon_{RL} = 0$ and $\Gamma_{RR}\Upsilon_{LL} + I = (1 - s\,g)\,I$ with $g = \Gamma_{ab}$, so the general Schur-complement update collapses to a rank-2 update of the complement block,

\[\Gamma_{ij} \leftarrow \Gamma_{ij} + \frac{s}{1 - s\,g}\left(\Gamma_{ia}\Gamma_{jb} - \Gamma_{ib}\Gamma_{ja}\right) \qquad (i, j \notin \{a, b\}),\]

after which rows and columns a, b are zeroed except $\Gamma_{ab} = -s$, $\Gamma_{ba} = s$. The denominator $1 - s\,g = 2\,P(\text{outcome})$ is twice the Born probability of the selected sector. $O(N^2)$ and allocation-free: no index-set complement, no slices, no lu/inv/cond on 2×2 blocks, and no BLAS calls — one contiguous sweep over the complement block plus the purity check. The closed form reads only $g = \Gamma_{ab}$ and assumes the antisymmetric structure $\Gamma_{ba} = -g$, $\Gamma_{aa} = \Gamma_{bb} = 0$ (a roundoff-level asymmetry perturbs the result at roundoff level); every operation in this backend leaves Γ exactly antisymmetric, and so does this one: the $(i,j)$ and $(j,i)$ updates use the same two products in opposite order, which round identically.

Throws ArgumentError for $1 - s\,g \le 0$ (a probability-zero outcome — the "vanishing state" of the general kernel); like gaussian_contraction!, calls purify! when the purity drift (_purity_drift) exceeds purify_tol. Agrees with the general kernel to roundoff ($\sim 10^{-15}$, not bitwise: different floating-point operation order).

parity_projection_upsilon(s::Int) -> Matrix{Float64}

4×4 Υ projecting a Majorana pair (i,j) onto the parity sector $i\gamma_i\gamma_j = -s$, $s \in \{+1,-1\}$, i.e. _kraus((s,0,0)). The post-measurement state has $\Gamma[i,j] = \langle i\gamma_i\gamma_j\rangle = -s$ (matching the Python reference's sign convention for s: the vacuum Γ[2i−1,2i] = +1 is the s = −1 outcome; contracting the vacuum pair with s = +1 is the probability-zero outcome and throws).

purify!(Γ::Matrix{Float64}) -> Matrix{Float64}

Project a (numerically drifted) covariance matrix back onto the manifold of pure Gaussian states ($\Gamma^2 = -I$) — port of purify (GTN.py:1131-1136, see App. B2 of PhysRevB.106.134206). Eigendecomposes the Hermitian matrix Γ/i, clamps eigenvalues to ±1, reconstructs $-\mathrm{Im}(V\,\mathrm{diag}(\mathrm{vals})\,V^\dagger)$, then antisymmetrizes. In-place; returns Γ.

record_value(obs, state; i1=nothing) -> Number | AbstractVector

Observable-level recording hook: compute the value record! should append for obs. The default is simply obs(state) — both for AbstractObservable instances and for generic callables (the untyped fallback), so ANY callable registered via track! is recorded without further ceremony.

Observables that need extra recording-time context override this method instead of being special-cased inside record! — e.g. DomainWall resolves its sampling site from i1_fn (registration-time closure) or the explicit i1 keyword. User-defined observables may override it the same way:

QuantumCircuitsMPS.record_value(obs::MyObs, state; i1=nothing) = ...

This is a MECHANISM hook only — it must return the value the observable's call protocol defines: a Float64 for the built-in scalar observables, or more generally any Number or AbstractVector (see track! for the storage contract).

reset!(::AbstractGeometry)

No-op base method for reset! on general geometries.

reset!(geo::AbstractStaircase)

Reset staircase position to its initial value.

resolve_probability_schedule(outcomes, K::Int) -> Matrix{Float64}

Resolve every outcome's probability field into a dense N_outcomes × K Matrix{Float64} (row i = outcome i, column k = element k). probability may also be an AbstractVector{<:Real} with one entry per geometry element (in elements(geo, L, bc) order), so each gate location gets its own probability; a scalar broadcasts across all K elements; mixed outcomes and K=1 vectors are fine.

Validates, in order: nonempty outcomes; vector length == K; type/shape (Real or AbstractVector{<:Real}); finiteness and [0,1] range; and per-column totals <= 1 + 1e-10 — the staircase/Pointer unit-sum guard stays with callers.

s1_dot_s2()

Compute the $S_1 \cdot S_2$ operator for two spin-1 particles. Returns a 9×9 matrix in the tensor product basis.

\[S_1 \cdot S_2 = S_z^{(1)}\otimes S_z^{(2)} + \frac{1}{2}\left(S_+^{(1)}\otimes S_-^{(2)} + S_-^{(1)}\otimes S_+^{(2)}\right)\]

select_outcome_index(rng::AbstractRNG, probs::AbstractVector{<:Real}) -> Int

SINGLE SOURCE OF TRUTH for the v0.1 unified stochastic rule's categorical selection. Draws exactly ONE scalar coin from rng and returns the 1-based index of the selected outcome, or 0 for the identity remainder ($r \ge \sum p$).

Semantics (bit-identical to test/reference_rule.jl's reference_select for a single element):

• one scalar rand(rng) per call — never vectorized draws
• cumulative walk over probs with strict <
• cumsum snapping: when abs(sum(probs) - 1) <= 1e-10, the LAST cumulative boundary is snapped to exactly 1.0, so float dust in $\sum p$ cannot leak spurious identity selections

probs accepts any AbstractVector{<:Real} — including a SubArray/ column view into a per-element probability-schedule Matrix{Float64} (see simulate!'s per-call schedule cache below) — in addition to a plain Vector{Float64}.

expand.jl's visualization path also delegates to this function — there is exactly one categorical-selection implementation in the package.

set_event_context!(state::SimulationState, step::Int, op_idx::Int, element_idx::Int)

Record the engine's current execution position on the state. simulate! calls this before every execute! invocation so that events emitted from deep call sites (e.g. MeasurementOutcome from the Born-measurement primitive) carry real (step, op_idx) indices. Feedback gates (Task 10) run inside their measuring gate's execute! and inherit its context — they never advance counters or context themselves.

Costs three integer stores; never allocates.

site_majoranas(state, site::Int) -> NTuple{N,Int}

Majorana indices carried by physical site of a Gaussian-backend SimulationState, RAM-mapped through state.phy_ram (identity on the Gaussian backend, kept for protocol uniformity). THE single source of truth for the site → Majorana index mapping — every Gaussian gate / measurement / observable resolves sites through this helper:

• fermionic-mode granularity (state.backend.majoranas_per_site == 2, default): returns (2r−1, 2r) — the pair of one fermionic mode.
• Majorana-chain granularity (majoranas_per_site == 1, site_type="Majorana"): returns (r,) — the site IS one Majorana mode.

Duck-typed on state (needs only phy_ram and backend.majoranas_per_site) so this pure-kernel file keeps no SimulationState dependency.

spin1_operators()

Return the spin-1 operators Sz, S+, S- as 3×3 matrices. Basis ordering: $\lvert +1 \rangle, \lvert 0 \rangle, \lvert -1 \rangle$ (descending m).

spin_operators(s) -> (Sz, Sp, Sm)

Return the spin-s operators Sz, S+, S- as dense $(2s+1)\times(2s+1)$ matrices in the descending-m basis $|s,s\rangle, |s,s-1\rangle, \ldots, |s,-s\rangle$ (matching ITensors' spin site conventions and spin1_operators).

Standard ladder formulas:

\[\begin{aligned} S_z|s,m\rangle &= m|s,m\rangle \\ S_\pm|s,m\rangle &= \sqrt{s(s+1) - m(m\pm1)}\,|s,m\pm1\rangle \end{aligned}\]

s may be any positive integer or half-integer (3//2, 1.5, 2, ...).

subsystem_entropy(Γ::AbstractMatrix{<:Real}, majorana_idx::AbstractVector{Int};
                  renyi_index::Real = 1) -> Float64

Rényi-n entropy IN NATS (natural log) of the reduced state of a fermionic Gaussian state on the subsystem spanned by the Majorana indices majorana_idx, computed from the full 2L×2L covariance matrix Γ. The default renyi_index = 1 is the von Neumann entropy.

Port of von_Neumann_entropy_m (GTN.py:753-759):

1. Γ_A = Γ[majorana_idx, majorana_idx] — covariance submatrix (plain fancy indexing; majorana_idx need not be contiguous, which is what makes this helper reusable for MutualInformation on disjoint regions).
2. ξ = eigvals(Hermitian(im .* Γ_A)) — real spectrum in [−1, 1] (i·Γ_A is exactly Hermitian since Γ is exactly antisymmetric).
3. λ = clamp.((1 .- ξ) ./ 2, 0.0, 1.0) — occupation eigenvalues in [0, 1] (clamped: float noise can push ξ marginally outside [−1, 1]).
4. $S = -\sum \left[\lambda \log \lambda + (1-\lambda)\log(1-\lambda)\right] / 2$ at n = 1, and generally $S_n = \sum \log(\lambda^n + (1-\lambda)^n) / (1-n) / 2$ — the division by 2 compensates the double-counting of the Majorana eigenvalues (they come in ±ξ pairs; Python divides the full sum by 2 the same way).

Equivalently, over the PAIRED part of $\mathrm{spec}(i\Gamma_A) = \{\pm\xi_k\}$ (each pair gives the occupation pair $\{\lambda, 1-\lambda\} = \{(1-\xi_k)/2, (1+\xi_k)/2\}$ from step 3),

\[S_n = \frac{1}{1-n} \sum_k \ln\!\left[\left(\frac{1+\xi_k}{2}\right)^n + \left(\frac{1-\xi_k}{2}\right)^n\right]\]

and an ODD-dimensional Γ_A (a Majorana-chain region) carries one extra UNPAIRED ξ = 0 mode contributing exactly ln(2)/2 to Sₙ for EVERY n — the per-eigenvalue sum above produces that with no special-casing.

renyi_index is evaluated in three branches of $\delta = n - 1$, at the SAME thresholds as the MPS / state-vector / mutual-information kernels (_RENYI_SHUNT, _RENYI_NEAR1): the exact n = 1 body, a cancellation-safe expm1/log1p near-1 form, and a scale-before-overflow log-domain form that stays finite at n = 2048 or n = floatmax(Float64).

Callers wanting a different log base divide the result by log(base).

support(gate::AbstractGate) -> Int

Return the number of sites this gate acts on (1 for single-qubit, 2 for two-qubit).

sync_staircase_positions!(outcomes, selected_geo::AbstractStaircase)

After advancing the selected staircase in a stochastic group, sync all other AbstractStaircase geometries in the group to the selected geometry's new position. This implements shared random-walk position for CIPT circuits.

vacuum_covariance(L::Int) -> Matrix{Float64}

2L×2L Majorana covariance matrix $\Gamma_0 = \bigoplus_i \begin{pmatrix}0&1\\-1&0\end{pmatrix}$ of the all-modes- unoccupied vacuum (port of correlation_matrix, GTN.py:29-55, non-op branch). Matches the Python reference get_C_f: every mode has $\langle c^\dagger c\rangle = 0$.

validate_geometry(geo::AbstractGeometry)

Validate that a geometry type is supported for circuit expansion.

Supported:

• StaircaseRight, StaircaseLeft
• SingleSite, AdjacentPair, Sites
• Bricklayer, AllSites, EachSite

Throws

ArgumentError if geometry is not supported (e.g. Pointer, whose position depends on runtime measurement outcomes and cannot be expanded statically).

validate_support(gate::AbstractGate, geo::Sites, L::Int, bc::Symbol)

Validate that support(gate) equals the region size of a Sites geometry; throws ArgumentError naming both numbers on mismatch. For all other geometries this is a no-op — their regions are resolved from the gate support itself (staircases) or are fixed-size by construction.

with_guarded_stream(f, registry::RNGRegistry, stream::Symbol)

Run f() with registry's stream temporarily replaced by a SentinelRNG: any draw from the guarded stream inside f throws an ErrorException ("... forbidden ..."), while all other streams remain fully usable. The original stream object is ALWAYS restored (try/finally), even if f throws. Returns f()'s value.

Used by the feedback system to guarantee that measurement feedback cannot consume :gates_spacetime coins (fixed-draw contract).