Observables

AbstractObservable

Abstract base type for observable specifications.

Observables are CALLABLE STRUCTS: an instance holds the observable's parameters (e.g. the entanglement cut) and is invoked on a state as obs(state) -> Float64. Register one for recording with track!(state, :name => obs); record! then appends obs(state) (via the record_value hook) to state.observables[:name].

Backend support is per-method multiple dispatch: the generic call (obs::MyObs)(state) implements the MPS path, and backend-specific methods (obs::MyObs)(state::SimulationState{StateVectorBackend}) / (obs::MyObs)(state::SimulationState{CliffordBackend}) override it. An observable that cannot be computed on a backend defines a method throwing an informative ArgumentError instead (e.g. StringOrder on the Clifford backend) — see docs/src/devdocs/backend_interface.md for the full support matrix.

User-defined observables subtype AbstractObservable, implement (obs::MyObs)(state) -> Float64, and are then accepted by track! and recorded like the built-ins. Observables needing extra recording-time context override the record_value hook (see DomainWall).

Subtyping is OPTIONAL for custom observables: track! accepts ANY callable f(state) -> Number | AbstractVector (a closure, a function, or a callable struct) — see track! for the full callable contract. Subtype AbstractObservable when you want the built-in scalar Float64 storage and the ability to override record_value.

EntanglementEntropy(; cut::Int, renyi_index::Int=1, threshold::Float64=1e-16, base::Real=2)

Entanglement entropy observable.

Computes the entanglement entropy across a specified cut in the MPS.

Arguments

• cut::Int: Physical site where the cut is made (must satisfy 1 <= cut < L)
• renyi_index::Int=1: Rényi index for entropy (must be >= 1)
  • renyi_index=1: von Neumann entropy S₁ = -Σᵢ λᵢ log(λᵢ) (default)
  • renyi_index=2: Rényi-2 entropy S₂ = log(Σᵢ λᵢ²) / (1-2)
  • renyi_index=n: Rényi-n entropy Sₙ = log(Σᵢ λᵢⁿ) / (1-n)
• threshold::Float64=1e-16: Minimum threshold for singular values (default: 1e-16)
• base::Real=2: Base of logarithm for entropy computation (default: 2 for bits)

Implementation

The entropy is computed by:

1. Converting the physical cut position to RAM ordering
2. Orthogonalizing the MPS at the cut site
3. Performing SVD to obtain Schmidt values
4. Computing the entropy from the Schmidt spectrum

Example

ee = EntanglementEntropy(; cut=2, renyi_index=1)
entropy = ee(state)

Magnetization(axis::Symbol)

Average single-site expectation value along axis.

Computes Mₐ = (1/L) Σᵢ ⟨axis_i⟩ where axis ∈ {:X, :Y, :Z}.

Example

track!(state, :Mz => Magnetization(:Z))

BornProbability(site::Int, outcome::Int)

Observable for Born rule probability P(measurement outcome | state) at a physical site.

born_probability(state::SimulationState{GaussianBackend}, site::Int, outcome::Int) -> Float64

Compute the Born probability of measuring occupation outcome (0 = unoccupied, 1 = occupied) at physical site for a fermionic Gaussian state.

For a Gaussian state the on-site occupation probability is an affine function of a SINGLE covariance-matrix element g = Γ[2i−1, 2i] (with i the RAM-mapped mode index):

P(0) = (1 + g) / 2        P(1) = (1 − g) / 2

(equivalently P(1) = ⟨cᵢ†cᵢ⟩, verified convention: vacuum has g = +1 ⇒ P(0) = 1.0 exactly). This is a NON-DESTRUCTIVE, read-only query — a single matrix-element read, no copies, state.backend.corr is untouched.

The raw affine value is returned without clamping: for a purified state |g| ≤ 1 up to machine precision, so any excursion outside [0, 1] is at the 1e-15 level and harmless to the rand() < p₀ threshold convention.

born_probability(state::SimulationState{MPSBackend}, physical_site::Int, outcome::Int) -> Float64

Compute Born probability P(outcome | state) at a physical site. Converts physical site to RAM index for MPS access.

This is the MPS-backend implementation. SimulationState{StateVectorBackend} and SimulationState{CliffordBackend} have their own, more specific overrides (see src/StateVector/measurement.jl, src/Clifford/measurement.jl); narrowing this signature to MPSBackend ensures any future/unknown backend gets a clear MethodError here instead of silently crashing on a backend-specific field (state.backend.mps) that doesn't exist.

born_probability(state::SimulationState{StateVectorBackend}, physical_site::Int, outcome::Int) -> Float64

Compute Born probability P(outcome | ψ) at physical_site for the state-vector backend.

For a state vector ψ of length d^L, sums |ψ_n|² over all basis integers n (0-indexed, 0 to d^L - 1) whose base-d digit at position physical_site equals outcome. Site 1 is MSB (slowest index).

Digit extraction: for 0-indexed basis integer n, physical site s (1-indexed, site 1 = MSB), and local dimension d: digit = (n ÷ d^(L-s)) % d (see _sv_digit; the loop-invariant stride d^(L-s) is hoisted out of the loop).

born_probability(state::SimulationState{CliffordBackend}, site::Int, outcome::Int) -> Float64

Compute the Born (Z-measurement) probability of outcome (0 or 1) at physical site for a stabilizer state.

For a stabilizer state, a single-qubit Z measurement is either DETERMINISTIC (probability exactly 0.0 or 1.0) or perfectly UNDETERMINED (probability exactly 0.5 for either outcome) — no other value is possible, a mathematical fact of the stabilizer formalism. This function is a NON-DESTRUCTIVE, read-only query: it operates on a copy of the tableau via QuantumClifford.projectZ!, so state.backend.tableau is left unmodified.

StringOrder(i::Int, j::Int; order::Int=1)

String order parameter observable for spin-1 chains.

Formulas

order=1 (nearest-neighbor AKLT):

O¹(i,j) = ⟨Sz[i] * exp(iπ Σ_{k=i+1}^{j-1} Sz[k]) * Sz[j]⟩

Expected: |O¹| ≈ 4/9 ≈ 0.444 for NN AKLT ground state

order=2 (next-nearest-neighbor AKLT):

O²(n,m) = ⟨Sz[n]·Sz[n+1] * exp(iπ Σ_{k=n+2}^{m-2} Sz[k]) * Sz[m-1]·Sz[m]⟩

Expected: |O²| ≈ (4/9)² ≈ 0.198 for NNN AKLT ground state

Arguments

• i: First site index (physical indexing)
• j: Second site index (physical indexing, must be j > i)
• order: 1 for NN formula (default), 2 for NNN formula with paired endpoints

Notes

• order=2 requires j >= i+4 (for non-overlapping endpoint pairs)
• NNN AKLT creates two decoupled chains; paired endpoints project onto both

Example

s = SimulationState(L=8, bc=:periodic, site_type="S=1")
initialize!(s, ProductState(binary_int=0))
# ... apply AKLT protocol ...
so1 = compute(StringOrder(1, 5), s)           # order=1 (default)
so2 = compute(StringOrder(1, 7, order=2), s)  # order=2

References

• AKLT (1987): Rigorous results on valence-bond ground states
• String order distinguishes Haldane phase from trivial phases

DomainWall(; order::Int, i1_fn::Union{Function, Nothing}=nothing)

Domain wall observable for CT model. Only supports xj=Set([0]) case (first "1" in bit string).

The domain wall computes: DW = Σ_j (L-j+1)^order * P(first "1" at position j starting from i1)

where position j is measured cyclically starting from i1.

Parameters:

• order: The order of the domain wall (≥ 1)
• i1_fn: Optional function that returns the sampling site i1 when called If provided, record! can be called without i1 parameter

PauliString(ops::Pair{Int,Symbol}...)

Expectation value ⟨∏ᵢ Pᵢ⟩ of a product of single-qubit Pauli operators.

Each argument is a site => pauli pair with pauli ∈ (:X, :Y, :Z); identity is implied on every site not listed. Sites must be distinct and positive; range validation against the system size happens at evaluation time. Supported on all three backends (MPS, state vector, Clifford).

Qubit-only in v0.4.0: evaluating on a non-qubit state (e.g. site_type="S=1") throws an ArgumentError (spin-S operator strings are on the roadmap).

Sign convention

Matches Magnetization: ⟨Zᵢ⟩ = +1 on |0⟩ and -1 on |1⟩ at site i.

Examples

state = SimulationState(L=4, bc=:open,
    rng=RNGRegistry(gates_spacetime=1, gates_realization=2, born_measurement=3))
initialize!(state, ProductState(binary_int=0))

PauliString(1 => :Z)(state)                    # ⟨Z₁⟩ = +1.0 on |0000⟩
PauliString(1 => :Z, 4 => :Z)(state)           # ⟨Z₁Z₄⟩ = +1.0
track!(state, :zz => PauliString(1 => :Z, 2 => :Z))  # record during simulate!

MutualInformation(regionA, regionB; renyi_index=1, threshold=1e-16, base=ℯ)

Mutual information I(A:B) = S(A) + S(B) - S(A∪B) between two site regions.

Each region is a non-empty, duplicate-free collection of physical sites (a UnitRange like 2:4, a single Int, or a Vector{Int}; stored sorted), and the two regions must be DISJOINT. Regions refer to physical sites on all backends and under all boundary conditions.

On the MPS, state-vector, and Clifford backends each region must additionally be CONTIGUOUS (a plain ascending range): non-contiguous input is accepted at construction (so one observable object can serve every backend) but rejected with an ArgumentError when EVALUATED on those backends. The Gaussian (covariance-matrix) backend supports arbitrary site subsets — including non-contiguous and PBC-wrapped regions such as [7, 8, 1, 2] at L=8.

Arguments

• regionA, regionB: the two site regions (contiguous, disjoint)
• renyi_index::Int=1: entropy index used for all three terms
  • renyi_index=1: von Neumann entropies (default)
  • renyi_index=n≥2: Rényi-n entropies. NOTE: the Rényi "mutual information" Iₙ = Sₙ(A)+Sₙ(B)−Sₙ(A∪B) is NOT guaranteed non-negative for n≠1 — it is a commonly used diagnostic, not a proper mutual information. Documented, not forbidden.
• threshold::Float64=1e-16: singular-value floor (probabilities are clamped at threshold^2 before taking logs), mirroring EntanglementEntropy
• base::Real=ℯ: logarithm base (default natural log, so a Bell pair gives I = 2·log(2) ≈ 1.386; use base=2 for bits)

Backend cost

• MPS: S(A∪B) for disjoint A, B requires a two-block reduced density matrix, contracted with cost/memory ~ χ²·d^(|A|+|B|) (χ = bond dimension). A size guard throws an informative ArgumentError when d^(|A|+|B|) > 256 (e.g. more than 8 qubits combined).
• StateVector: exact dense partial trace — practical for L ≲ 20 only (memory/time scale as d^L).
• Clifford: poly-time GF(2)-rank stabilizer entropies; every Rényi index gives the same value (flat entanglement spectrum).
• Gaussian: three covariance-submatrix eigendecompositions, O((|A|+|B|)³) — arbitrary site subsets, von Neumann only.

Properties (analytic anchors)

• Product state: I = 0
• Pure global state with B = complement(A): I = 2·S(A)
• Bell pair, A = {1}, B = {2}: I = 2·log(2)
• GHZ(4), A = {1}, B = {4}: I = log(2)

Examples

mi = MutualInformation(1:2, 5:6)            # blocks {1,2} and {5,6}
mi = MutualInformation(1, 4; base=2)        # single sites, result in bits
value = mi(state)
track!(state, :I => MutualInformation(1, 4))

Correlator(pi::Pair{Int,Symbol}, pj::Pair{Int,Symbol})

Connected two-point correlator C(i,j) = ⟨PᵢPⱼ⟩ − ⟨Pᵢ⟩⟨Pⱼ⟩.

Each argument is a site => pauli pair with pauli ∈ (:X, :Y, :Z), matching PauliString's argument style. The two sites must be DISTINCT (the self-correlation C(i,i) with Pᵢ² = I is identically 1 − ⟨Pᵢ⟩² and is not computed here — an ArgumentError is thrown for i == j).

Supported on all three backends (composition of three PauliString evaluations, which dispatch per backend). Qubit-only in v0.4.0 (inherited from PauliString).

Properties (analytic anchors)

• Product state |0…0⟩: C(i,j) = 1 − 1·1 = 0 for Z operators
• Bell pair, Correlator(1 => :Z, 2 => :Z) = 1 − 0·0 = 1

Examples

c = Correlator(1 => :Z, 4 => :Z)
value = c(state)
track!(state, :czz => Correlator(1 => :Z, 2 => :Z))

EntropyProfile(; renyi_index=1, threshold=1e-16, base=ℯ)

Entanglement-entropy profile: the vector [S(cut=x) for x in 1:L-1] of bipartite entropies at every cut, computed by the existing per-cut EntanglementEntropy on each backend.

Arguments

• renyi_index::Int=1: Rényi index for all cuts (1 = von Neumann)
• threshold::Float64=1e-16: singular-value floor (see EntanglementEntropy)
• base::Real=ℯ: logarithm base (default natural log — NB: EntanglementEntropy itself defaults to base=2; pass base=2 for bits)

Backend cost

• MPS: O(L) orthogonalized MPS copies — O(L²·χ³) total
• StateVector: O(L) dense reshapes + SVDs
• Clifford: O(L) tableau copies + GF(2)-rank computations

PBC caveat (cross-backend semantics)

On the MPS backend under bc=:periodic, each cut is the RAM bond index of the FOLDED-PBC MPS (src/Observables/entanglement.jl), NOT the physical prefix bipartition {1..cut} used by the state-vector and Clifford backends (only cut = L÷2 is fold-aligned). A periodic-BC MPS profile is therefore in RAM-bond coordinates and is NOT directly comparable to the other backends' physical-cut profiles. Cross-backend profile comparisons must use bc=:open, where all three backends agree on the physical bipartition.

Recording

Returns a Vector{Float64} (one entry per cut). When tracked, each record point appends the whole vector as ONE entry; the observable storage is transparently widened from Vector{Float64} to Vector{Any} at the first record (see track!).

Examples

ep = EntropyProfile(; base=2)
profile = ep(state)                 # Vector{Float64} of length L-1
track!(state, :Sx => EntropyProfile())

TripartiteMutualInformation(A, B, C; renyi_index=1, threshold=1e-16, base=ℯ)

Tripartite mutual information I₃ = I(A:B) + I(A:C) − I(A:BC) between three site regions (Gullans–Huse MIPT convention).

Each region must be a CONTIGUOUS, non-empty range of physical sites, the three regions must be pairwise DISJOINT, and B ∪ C must itself be contiguous (i.e. B and C adjacent) — inherited from MutualInformation's contiguous-region constraint on the I(A:BC) term. Note that A and C need NOT be adjacent: I(A:C) uses MutualInformation's two disjoint blocks.

Sign convention (MIPT usage)

The standard diagnostic partitions the chain into four QUARTERS A, B, C, D and computes I₃(A:B:C) with D traced out — NOT equal thirds of the whole system: for ANY pure state, a full tripartition (A∪B∪C = everything) gives I₃ ≡ 0 identically, so it carries no information (useful only as a consistency check). With a fourth region traced out, scrambled volume-law states give I₃ ≤ 0 (the MIPT scrambling diagnostic), while e.g. GHZ-type global correlations give I₃ > 0: GHZ(L=8) with A=1:2, B=3:4, C=5:6 (D=7:8 traced) gives I₃ = +log 2.

Arguments

• A, B, C: the three site regions (contiguous, pairwise disjoint, B and C adjacent)
• renyi_index::Int=1, threshold::Float64=1e-16, base::Real=ℯ: forwarded to all three MutualInformation terms (see its docstring; Rényi-n I₃ for n≥2 is a diagnostic, not a proper mutual information)

Backend cost

Inherits MutualInformation: on the MPS backend the largest term I(A:BC) requires d^(|A|+|B|+|C|) ≤ 256 (e.g. ≤ 8 qubits combined); state vector is exact dense (L ≲ 20); Clifford is poly-time.

Examples

tmi = TripartiteMutualInformation(1:2, 3:4, 5:6)   # quarters of L=8
value = tmi(state)
track!(state, :I3 => TripartiteMutualInformation(1:2, 3:4, 5:6; base=2))

MagnetizationFluctuations(region; axis=:Z)

Variance of the total magnetization M = Σ_{i∈R} Pᵢ over the site region R, with P the Pauli operator selected by axis (:X, :Y, or :Z).

Computed as

Var(M) = |R| + Σ_{i≠j} ⟨PᵢPⱼ⟩ − (Σᵢ ⟨Pᵢ⟩)²

using Pᵢ² = I for the diagonal of ⟨M²⟩ (the analytic constant |R|), so the formula never evaluates a same-site Pauli-string product. This is O(|R|²) PauliString evaluations. Supported on all three backends (composition of PauliString, which dispatches per backend); qubit-only in v0.4.0 (inherited from PauliString).

region may be any collection of distinct positive sites (a range like 1:6, an Int, or a vector — contiguity is NOT required, since M is a sum of single-site operators).

Properties (analytic anchors)

• Product Z-basis state, axis=:Z: Var = 0 (M is sharp)
• GHZ(L=6), R=1:6, axis=:Z: ⟨ZᵢZⱼ⟩=1 for all i≠j and ⟨Zᵢ⟩=0, so Var = 6 + 30 − 0 = 36 exactly

Examples

vm = MagnetizationFluctuations(1:6)          # Var of Σ Zᵢ over sites 1..6
value = vm(state)
track!(state, :varM => MagnetizationFluctuations(2:5; axis=:X))

track!(state::SimulationState, spec::Pair{Symbol, <:Any})

Register an observable to be tracked. Values are stored in state.observables[name], one entry per record point.

The pair's value may be a built-in observable instance (any AbstractObservable, e.g. EntanglementEntropy, Magnetization, PauliString) — or any callable f(state): a closure, a plain function, or a user-defined callable struct.

The callable contract

At every record point (an eager record!(state) call, or the points selected by simulate!'s record_when policy), each tracked observable is invoked as f(state) (through the record_value hook) and the returned value is appended to state.observables[name]:

• f MUST be side-effect-free with respect to the quantum state: it may freely READ state (e.g. via born_probability, PauliString, EntanglementEntropy, measurements(state)), but must not apply gates, measure, or consume RNG draws.
• f may return a scalar (Number) or a vector (AbstractVector), e.g. a per-site profile. Each record point appends exactly ONE entry (a returned vector is stored as a single element, not splatted).
• Errors thrown by f propagate out of record!. For non-AbstractObservable callables the error is wrapped in an ErrorException naming the observable key, with the underlying exception attached as caused by:.

Storage

• AbstractObservable specs record into a Vector{Float64} (the built-in scalar contract). Should such an observable return a non-scalar value, the storage is transparently widened to Vector{Any} at that record.
• Generic callables record into a Vector{Any} (their returns are not constrained to Float64).

Examples

track!(state, :dw1 => DomainWall(order=1))            # built-in observable
track!(state, :p1 => s -> born_probability(s, 1, 0))  # custom closure
track!(state, :zprofile =>
    s -> [PauliString(i => :Z)(s) for i in 1:s.L])    # vector-valued closure

See the "Custom Observables" documentation page for worked examples and the public building blocks to compose.

record!(state::SimulationState; i1::Union{Int,Nothing}=nothing, only=nothing)

Compute tracked observables and append values to state.observables.

Each observable's value is obtained through the record_value hook (default obs(state); DomainWall uses i1_fn/i1 — see record_value). Errors thrown by a tracked observable propagate; for generic callables (non-AbstractObservable specs, see track!) they are wrapped in an ErrorException naming the observable key.

Keyword Arguments

• i1: Explicit DomainWall sampling site (required when the tracked DomainWall was registered without an i1_fn)
• only: Optional collection of tracked-observable names (Symbols). When given, ONLY those observables are recorded — the others' vectors do not grow. Naming an untracked observable throws an ArgumentError. Used by selective record!(c, names...) circuit markers.

record!(builder::CircuitBuilder)
record!(builder::CircuitBuilder, names::Symbol...)

Insert a recording MARKER into the circuit at this position (v0.1).

Stores a pseudo-operation: (type=:record_mark, names=Symbol[names...]).

Markers are pure annotations: they never draw from any RNG stream, never advance the gate_idx element-slot counter, and never touch staircase positions. They fire only under record_when=:marks (which records the tracked observables exactly at each marker position) or inside custom predicates (which receive RecordingContexts with at_mark=true).

With no names, ALL tracked observables are recorded at the marker. With explicit names (e.g. record!(c, :entropy)), only the named observables are recorded there — each observable's vector grows at its own cadence. Naming an observable that is not tracked on the state raises an ArgumentError at simulate! time.

simulate! refuses (ArgumentError) to run a marker-containing circuit under record_when ∈ (:every_step, :every_gate, :final_only) — those policies would silently ignore the markers.

Example

circuit = Circuit(L=L, bc=:periodic) do c
    apply!(c, HaarRandom(), Bricklayer(:even))
    apply_with_prob!(c; outcomes=[(probability=p, gate=Measure(:Z), geometry=EachSite(2:L-1))])
    record!(c)              # record all tracked observables here
    apply!(c, HaarRandom(), Bricklayer(:odd))
    record!(c, :entropy)    # record only :entropy here
end
simulate!(circuit, state; n_steps=25, record_when=:marks)

list_observables() -> Vector{String}

Return a list of available observable type names.

Returns the names of all observable types that can be used with the tracking API.

Example:

obs_types = list_observables()
# Returns: ["DomainWall", "BornProbability"]