This document teaches the stabilizer formalism the way one would teach it at a chalkboard with a kernel open: every claim is immediately computed. We start from the observation that a single qubit can be specified not by its two amplitudes but by the one operator that fixes it, build that idea into the Aaronson-Gottesman tableau, and ride it all the way to thousand-qubit Clifford circuits, error-correcting codes, graph states, magic states, and Clifford channels. We close with a measured, honest comparison against the specialized C++ and Julia stabilizer simulators. The arc is: a state is what stabilizes it (single qubit),
n
generators replace
n
2
amplitudes (the tableau), entanglement and measurement become linear algebra over
𝔽
2
(the field of two elements
{0,1}
with all arithmetic done mod 2, so that
1+1=0
and vector addition is the bitwise XOR), and then scale, codes, and the frontier. The audience is graduate students who know quantum mechanics and the Pauli group but have not yet seen why a stabilizer state of a thousand qubits is something you can hold in your hand. Every construction in this document is built explicitly so you can change it, break it, and extend it. The Wolfram QuantumFramework symbols we exercise live in the Wolfram QuantumFramework paclet; the central one is PauliStabilizer.

Setup

We load the framework from the working tree so the fast paths and recent fixes are present. On a current release, the plain
Needs
line alone suffices.
PacletInstall["https://www.wolfr.am/DevWQCF",ForceVersionInstall->True];​​Needs["Wolfram`QuantumFramework`"]

The core functions at a glance

Before the narrative, the whole working surface, grouped by task. Every call is exercised in context where it is introduced; this is the map.
◼
  • Construct a register, from signed Pauli strings, a named code, or a random Clifford:
    PauliStabilizer[n]
    ,
    PauliStabilizer[{"XX","ZZ"}]
    ,
    PauliStabilizer["SteaneCode"]
    ,
    PauliStabilizer["Random"[n]]
    .
    • ◼
  • Apply a Clifford gate one at a time:
    ps["H",j]
    ,
    ps["S",j]
    ,
    ps["CNOT",j,k]
    ,
    ps["CZ",j,k]
    ,
    ps["SWAP",j,k]
    .
    • ◼
  • Apply a whole circuit as a gate list, routed through the compiled engine (bare names take default targets,
    ->
    sets explicit ones):
    ps[{"H","CNOT","CNOT"->{2,3}}]
    .
    • ◼
  • Inspect:
    ps["Stabilizers"]
    ,
    ps["Tableau"]
    ,
    ps["Matrix"]
    ,
    ps["State"]
    .
    • ◼
  • Measure:
    ps["M",q]
    ,
    ps["M","ZZ"]
    ,
    ps["Expectation","ZZ"]
    ,
    ps["SymbolicMeasure",q]
    .
    • ◼
  • Entanglement entropy across a cut:
    ps["Entropy",{1,2}]
    .
    • ◼
  • Sibling heads:
    StabilizerFrame
    ,
    CliffordChannel
    ,
    GraphState
    with
    LocalComplement
    , and the predicate
    StabilizerStateQ
    .

    • Why stabilizers: the exponential wall

    A pure state of
    n
    qubits is a unit vector in
    n
    2
    ℂ
    , so writing it down means storing
    n
    2
    complex amplitudes. At
    n=20
    that is already a million numbers; at
    n=50
    it is about
    15
    10
    , a petabyte of amplitudes even at one byte each; by
    n=100
    it passes the number of atoms in a building (
    ∼
    30
    10
    ). Yet many of the states we care most about, including every state any error-correcting code protects, are not generic: they sit at the fixed points of a large group of Pauli operators, and that structure can be stored in
    O(
    2
    n
    )
    bits.
    To feel the wall, ask for the dense state vector of a 20-qubit GHZ state and look only at its length:
    In[]:=
    Length@QuantumState["GHZ"[20]]["StateVector"]
    Out[]=
    1048576
    That is
    20
    2
    amplitudes for a state that is pinned down completely by 20 short Pauli strings. The stabilizer formalism is the bookkeeping that exploits exactly this gap.

    A state is what fixes it: the single qubit

    The central move is to stop describing a state by its amplitudes and start describing it by the operators that leave it invariant. The state
    |0〉
    is the unique
    +1
    eigenvector of the Pauli operator
    Z
    . We say
    Z
    stabilizes
    |0〉
    , and since one operator already determines a one-qubit state up to normalization, that single fact is the state.
    A caution that defines the whole subject: this works only for special states. Of all one-qubit pure states, exactly six are
    +1
    eigenstates of a Pauli:
    |0〉
    and
    |1〉
    for
    ±Z
    ,
    |±〉
    for
    ±X
    , and
    |i〉
    for
    ±Y
    . These six are the stabilizer states, and the formalism is exactly their theory and its multi-qubit generalization; a generic
    cosθ|0〉+sinθ|1〉
    is fixed by no Pauli and falls outside it. What makes the restriction worth a whole subject is that the states it does cover are exactly the ones error correction and large parts of quantum algorithms are built from.
    This deserves a closer look, because "fixed by no Pauli" is not the same as "fixed by nothing." Every pure state is the
    +1
    eigenstate of some operator. The reflection
    2|ψ〉〈ψ|-I
    does the job for any
    |ψ〉
    at all, leaving
    |ψ〉
    alone and flipping the sign of everything orthogonal to it:
    In[]:=
    psi={Cos[1/3],Sin[1/3]};​​reflection=2Outer[Times,psi,psi]-IdentityMatrix[2];​​Simplify[reflection.psi]==psi
    Out[]=
    True
    That stabilizing operator is a dense, structureless thing that turns continuously as the state turns:
    In[]:=
    Simplify[reflection]//MatrixForm
    Out[]=
    Cos
    2
    3
    
    Sin
    2
    3
    
    Sin
    2
    3
    
    -Cos
    2
    3
    
    So "a state is what fixes it" is, on its own, vacuous: this reflection takes the full
    n
    2
    amplitudes to write down, exactly the cost we are fleeing, and a gate turns it into some other arbitrary operator. The power comes entirely from restricting the fixing operators to the Pauli group, which is special on four counts that only pay off together:
    ◼
  • Cheap to name: a Pauli string is
    O(n)
    bits, two per qubit and a sign, not
    n
    2
    amplitudes.
    • ◼
  • Clean spectrum: it squares to the identity, so its eigenvalues are
    ±1
    , and "fixed by
    P
    " sharpens to the binary fact "
    +1
    eigenstate."
    • ◼
  • Commute or anticommute: any two Paulis do one or the other, with nothing in between, so
    n
    independent commuting Paulis halve the
    n
    2
    -dimensional space
    n
    times and leave a unique joint
    +1
    eigenstate, which is why a list of
    n
    generators pins down one state exactly.
    • ◼
  • Closed under Clifford gates: conjugation by
    H
    ,
    S
    , and
    CNOT
    sends Paulis to Paulis, so the description never leaves the Pauli language as the state evolves.
    • The reflection
    2|ψ〉〈ψ|-I
    has none of these properties. Drop any one of them and naming a state by what fixes it stops paying its way; keep all four and it becomes the
    O(
    2
    n
    )
    tableau the rest of this document runs on.
    Verify that
    Z
    fixes
    |0〉
    , recalling that the built-in
    PauliMatrix[1]
    ,
    PauliMatrix[2]
    ,
    PauliMatrix[3]
    are
    X
    ,
    Y
    ,
    Z
    :
    In[]:=
    PauliMatrix[3].{1,0}=={1,0}
    Out[]=
    True
    In the framework, the one-qubit register
    |0〉
    is
    PauliStabilizer[1]
    . The object prints as a summary box that shows its generators and, for small systems, its tableau, so the cheap representation is the thing you see:
    In[]:=
    PauliStabilizer[1]
    Out[]=
    PauliStabilizer
    
    {Z}
    (
    0
    1
    )
    
    We read off the operator that stabilizes it with the
    "Stabilizers"
    property:
    In[]:=
    PauliStabilizer[1]["Stabilizers"]
    Out[]=
    {Z}
    The answer is the single Pauli string
    "Z"
    , exactly as the chalkboard argument predicted. The
    "State"
    property is the bridge back to the rest of the framework: it materializes the stabilizer object into a full
    QuantumState
    , at the cost of the
    n
    2
    blow-up the formalism exists to avoid.
    In[]:=
    PauliStabilizer[1]["State"]
    Out[]=
    QuantumState
    Pure state
    Qudits: 1
    Type: Vector
    Dimension: 2
    
    What comes back is an ordinary
    QuantumState
    object, so a stabilizer state drops straight into the rest of the framework. Chaining
    "StateVector"
    on that object reads off the raw
    n
    2
    amplitudes, here the single column
    |0〉
    :
    In[]:=
    PauliStabilizer[1]["State"]["StateVector"]//Normal
    Out[]=
    {1,0}
    Different states are fixed by different operators. The state
    |+〉=(|0〉+|1〉)
    2
    is the
    +1
    eigenstate of
    X
    , and a Hadamard turns one into the other. Apply
    H
    to
    |0〉
    and read the new stabilizer:
    In[]:=
    PauliStabilizer[1]["H",1]["Stabilizers"]
    Out[]=
    {X}
    The stabilizer flips from
    Z
    to
    X
    , the signature of
    |+〉
    . The state vector confirms it:
    In[]:=
    PauliStabilizer[1]["H",1]["State"]["StateVector"]//Normal
    Out[]=
    
    1
    2
    ,
    1
    2
    
    The same move with
    X
    in place of
    H
    carries
    |0〉
    to
    |1〉
    , and here a sign appears for the first time:
    |1〉
    is the
    +1
    eigenstate of
    -Z
    , not of
    Z
    .
    In[]:=
    PauliStabilizer[1]["X",1]
    Out[]=
    PauliStabilizer
    
    {-Z}
    -
    (
    0
    1
    )
    
    Reading off the stabilizer makes the sign explicit:
    In[]:=
    PauliStabilizer[1]["X",1]["Stabilizers"]
    Out[]=
    {-Z}
    The answer is the signed string
    "-Z"
    , the operator that fixes
    |1〉
    . The sign is part of the stabilizer data, which is exactly why the formalism can tell
    |0〉
    from
    |1〉
    even though the same unsigned
    Z
    pins both.
    A short word of honesty about what the tableau leaves out. The tableau tracks the stabilizer group, including the signs, but it is blind to the overall phase of the state. Reconstructing a state vector from a gate-updated tableau is therefore exact only up to a global phase. You can see this by acting with
    Z
    on
    |+〉
    : physically
    Z|+〉=|-〉
    , and the framework returns
    |-〉
    up to an overall sign.
    In[]:=
    PauliStabilizer[1][{"H","Z"}]["State"]["StateVector"]//Normal
    Out[]=
    -
    1
    2
    ,
    1
    2
    
    For everything observable, probabilities, expectation values, entanglement, the global phase is irrelevant, so this costs us nothing in practice. We flag it now because it is the one place where the cheap representation is genuinely less informative than the dense one.
    The essentials of the one-qubit picture:
    ◼
  • A one-qubit stabilizer state is named by the single Pauli operator that fixes it:
    Z
    for
    |0〉
    ,
    X
    for
    |+〉
    ,
    -Z
    for
    |1〉
    .
    • ◼
  • PauliStabilizer
    carries that operator, not the amplitudes; the
    "State"
    property hands back a full
    QuantumState
    on demand, and
    "StateVector"
    on that object gives the raw
    n
    2
    amplitudes.
    • ◼
  • The reconstructed vector is exact up to a global phase, which never affects an observable.

    • The tableau:
    n
    generators instead of
    n
    2
    amplitudes

    For
    n
    qubits the idea scales because the stabilizer group is generated by just
    n
    commuting, independent Pauli strings, and each Pauli string of length
    n
    is two bits per qubit: an
    X
    -exponent and a
    Z
    -exponent, since
    Y=iXZ
    . So a stabilizer state is a binary matrix of shape
    n×2n
    together with
    n
    signs, which is
    O(
    2
    n
    )
    bits in place of
    n
    2
    amplitudes. This binary
    (XZ)
    encoding is the symplectic representation of the Pauli group.
    Look at the two-qubit register
    |00〉
    . The object itself previews the generators and a tableau grid:
    In[]:=
    PauliStabilizer[2]
    Out[]=
    PauliStabilizer
    
    {ZI,IZ}
    
    0
    0
    1
    0
    0
    0
    0
    1
    
    
    The raw symplectic data is
    "Tableau"
    , a rank-3 array stacking the
    X
    block above the
    Z
    block (shape
    2×n×2n
    , here
    2×2×4
    : two blocks,
    n=2
    qubits,
    2n=4
    rows):
    In[]:=
    PauliStabilizer[2]["Tableau"]
    Out[]=
    {{{1,0,0,0},{0,1,0,0}},{{0,0,1,0},{0,0,0,1}}}
    Its dimensions confirm the
    {2,n,2n}={2,2,4}
    layout directly:
    In[]:=
    PauliStabilizer[2]["Tableau"]//Dimensions
    Out[]=
    {2,2,4}
    The same content reads more easily flattened to a binary matrix. Each row is one generator written across as
    (
    x
    1
    …
    x
    n
    
    z
    1
    …
    z
    n
    )
    : the left half of the columns are the
    X
    -exponents, the right half the
    Z
    -exponents. For
    |00〉
    it comes out as the
    2n×2n
    identity:
    In[]:=
    PauliStabilizer[2]["Matrix"]//MatrixForm
    Out[]=
    1
    0
    0
    0
    0
    1
    0
    0
    0
    0
    1
    0
    0
    0
    0
    1
    There is a subtlety worth pausing on: the matrix has
    2n
    rows, not
    n
    , because the framework keeps two interleaved bases. The bottom
    n
    rows are the stabilizer generators proper, the operators that fix the state; the top
    n
    are destabilizer generators, a symplectic-dual partner basis introduced by Aaronson and Gottesman, each anticommuting with exactly one stabilizer and commuting with all the rest. The destabilizers cost nothing asymptotically and turn measurement from an exponential search into an
    O(
    2
    n
    )
    tableau update.
    The rows are bits, but each one decodes back to a signed Pauli string, and pairing them makes the correspondence concrete. Read in storage order, destabilizers first and then stabilizers:
    In[]:=
    MapThread[Rule,{​​PauliStabilizer[2]["Matrix"],​​Join[PauliStabilizer[2]["Destabilizers"],PauliStabilizer[2]["Stabilizers"]]​​}]//Column
    and read off its amplitudes:

    Two pictures of one circuit: Schrödinger and Heisenberg

    Its diagram makes the gate order visible:
    Its amplitudes are the answer the formalism has to reproduce:

    Entanglement for free: Bell and GHZ

    Build the Bell state:
    and read its generators:
    Its single-qubit entanglement entropy is still one bit:
    Time the entropy of this thousand-qubit GHZ state to see that it is essentially instantaneous:
    What we have established:

    Measurement: randomness, back-action, and correlation

    A measurement on a stabilizer state returns an association from outcome to post-state. For the Bell state, measuring qubit 1 is random, so both branches appear:
    Measure all three qubits and keep just the outcome keys:
    Its generators now carry the formal outcome in a sign:
    Two signs now carry the free bits; the third generator is still the bare parity check. A small helper reads the bit-string a resolved sample collapses to, by locating its one nonzero amplitude:
    Substituting the two free bits over all four assignments traces out the determined third outcome, which is in every case the XOR of the first two:
    Only the four even-parity strings appear, each about a quarter of the shots: the randomness is genuine, but the parity the circuit built is enforced on every sample, carried from the one traversal by substitution alone.

    Scaling up: random Cliffords and compiled circuits

    then count its gates, which stay polynomial in the number of qubits:
    This inverts gate-folding, which goes from circuit to tableau, and it is the workhorse behind code encoders: a code state's preparation circuit is exactly its encoding circuit.
    then fold it and watch the wall-clock:
    run it through the stabilizer backend:
    and read its generators, the familiar Bell pair:
    Run it densely:
    run it in the tableau:

    Error-correcting codes as stabilizer states

    A bit-flip on qubit 2 anticommutes with the first check alone:
    Every error produces a different, nonzero syndrome. In fact the fifteen errors realize exactly the fifteen nonzero four-bit strings, none missing and none repeated:

    The frontier: graph states, magic, and channels

    Three directions take the formalism to the edge of current research, and the framework has a head for each. This section gives each one enough room to compute with, not just name.

    Graph states

    and read its generators:
    The graph is updated correctly, but the accompanying single-qubit Clifford corrections, the vertex operators that would make the equivalence exact at the state level, are not yet tracked. Use this for graph surgery rather than for exact state identity.

    Magic and stabilizer rank

    At twenty qubits the frame carries eight stabilizer terms where the amplitudes would need over a million: it is exponentially smaller than the state it represents, sized by the magic content alone. That compactness is the entire reason stabilizer-rank methods exist, and the relating Pauli on each component keeps the superposition phase-coherent, so the frame is also exact.

    Clifford channels

    A pure stabilizer state is itself a state-preparation Clifford channel, whose tableau is the state's stabilizer half. Converting the Bell state to its channel recovers its two parity-check rows:

    An honest comparison with other packages

    The stabilizer simulators built for raw throughput are Stim, a C++ engine, and QuantumClifford.jl, a Julia package; both are excellent and both are narrower in scope than a general quantum framework. The fair question is two-sided: how close is the framework on the task they specialize in, and what does it offer that they do not.
    On gate-folding throughput the gap is now a small constant factor. The figures here are the measured wall-clock for an identical Clifford stream on an Apple M2 Pro, measured June 2026 from the framework's own bottleneck audit, with Stim for reference; treat them as a dated snapshot, not a guarantee:
    At a thousand qubits the framework also runs ahead of QuantumClifford.jl on the same stream (58.0 ms against 108.4 ms). The remaining gap to Stim is SIMD lane width, Stim packs 256-bit vector lanes against the framework's 62-bit machine words, plus a small per-gate circuit-encoding pass; it is a constant factor, not an algorithmic deficit. Getting there took compiling the gate fold to a single C kernel: the gate-by-gate interpreted path is 15 to 30 times slower, around 74, 402, and 955 ms at 100, 500, and 1000 qubits.
    What the framework offers that the dedicated engines do not:
    ◼
  • Symbolic signs. Outcomes can be carried as formal symbols through a circuit (the SymPhase method), so a single traversal serves arbitrarily many samples and the correlations between outcomes come out as explicit polynomials; in the dedicated engines the same device propagates Pauli faults to syndromes, which is the sampling bottleneck the technique was built to break.
    • ◼
  • Exact arithmetic. Coefficients, signs, and the random-Clifford sampler use exact rationals, so seeded output is bit-reproducible and there is no floating-point drift.
    • ◼
  • Breadth of heads. Graph states with local complementation, stabilizer frames for magic accounting, and Clifford channels for open-system Clifford dynamics are first-class objects, not features one bolts on.
    • Where the framework genuinely trails, stated plainly:
    ◼
  • It has no constant-cost bulk Pauli-frame sampler, the trick that makes Stim's shot-by-shot sampling so cheap.
    • ◼
  • Reconstructing an exact state vector from a gate-updated tableau, or from a gate-built stabilizer frame, is correct up to a single global phase, which is never observable; only the relative phase between two separately-built objects is gauge-limited.
    • The summary is the one a physicist would want: for raw Clifford throughput the specialized engines are a few times faster, and if that is the only thing one needs, they are the right tool. For stabilizer simulation that has to interoperate with the rest of a quantum-mechanics toolkit, carry symbols, stay exact, and reach into codes, graph states, magic, and channels, the framework buys all of that for a small constant in speed.

    Where this leaves us

    The key references underpinning the implementation are Gottesman's stabilizer-codes thesis (1997), Aaronson and Gottesman on the tableau and its measurement update (2004), Fattal and collaborators on stabilizer entanglement entropy (2004), Anders and Briegel on graph states and local complementation (2005), Koenig and Smolin on the random-Clifford sampler (2014), Garcia and Markov on stabilizer frames (2015), Fang and Ying on symbolic phases (2023), and Yashin on Clifford channels (2025).

    CITE THIS NOTEBOOK

    Stabilizer formalism via computation: Aaronson-Gottesman tableaux, Clifford circuits and QEC codes​
    by Mohammad Bahrami​
    Wolfram Community, STAFF PICKS, June 24, 2026
    https://community.wolfram.com/groups/-/m/t/3738980