Gaussian Backend

QuantumCircuitsMPS.jl also ships a fermionic Gaussian (free-fermion / FLO — "fermionic linear optics") backend for circuits built entirely out of Gaussian-preserving gates and parity measurements. Instead of an MPS, a dense state vector, or a stabilizer tableau, the state is stored as a $2L \times 2L$ real antisymmetric Majorana covariance matrix $\Gamma$, and gate/measurement application is a Schur-complement tensor contraction on $\Gamma$ rather than any operation on a Hilbert-space vector. apply!, track!, record!, and simulate! all work exactly as on the other three backends; only the SimulationState(...) constructor call changes.

When to use it: circuits built entirely out of fermionic Gaussian unitaries and parity measurements — measurement-induced phase transitions in free-fermion (class-D/DIII) circuits, entanglement transitions with free fermions, and the staggered class-DIII Majorana-chain protocol of Pan et al. Polynomial-time simulation ($O(L^3)$ per gate), no truncation error, exact for any circuit depth.

When not to use it: any circuit that needs a non-Gaussian gate — Hadamard, CNOT, HaarRandom, RandomClifford, generic Projection, or anything from the qubit-gate vocabulary. Use backend=:mps (see MPS Backend) or backend=:statevector (see State Vector Backend) for those. The Gaussian backend is fermionic-mode-only (local_dim=2); arbitrary spin-S sites (see Arbitrary Spin-S Support) are MPS/state-vector-only.

Quick Example

using QuantumCircuitsMPS

state = SimulationState(L=32, bc=:periodic, backend=:gaussian,
    rng=RNGRegistry(gates_spacetime=1, gates_realization=2, born_measurement=3, state_init=4))
initialize!(state, ProductState(binary_int=0))          # vacuum
apply!(state, GaussianHaar(), Bricklayer(:odd))          # random Gaussian circuit
apply!(state, Measure(:Z), SingleSite(5))                # on-site occupation parity measurement
apply!(state, BondParity(), AdjacentPair(3))             # bond parity measurement
track!(state, :entropy => EntanglementEntropy(cut=16))
record!(state)

entropies = state.observables[:entropy]
println("Final entropy: $(entropies[end])")

Physics: GaussianHaar() draws an independent Haar-random $O(4)$/$O(2)$ rotation for each bond (from :gates_realization) and conjugates it directly onto the Majorana covariance matrix — no dense Hilbert-space unitary is ever built. Measure(:Z) and BondParity() are projective parity measurements of, respectively, the on-site occupation parity $i\gamma_{2i-1}\gamma_{2i} = 1 - 2c_i^\dagger c_i$ and the bond parity spanning two adjacent sites; both collapse $\Gamma$ with the closed-form rank-2 update that the fermionic-linear-optics contraction formula reduces to for a parity projector ($O(N^2)$ per measurement, no linear solve; parity_projection! in the internals), so a projective measurement costs about the same as a Haar gate.

What a Gaussian State Is

A pure fermionic Gaussian state on L modes is completely characterized (no exponentially large amplitude vector needed) by its $2L \times 2L$ real antisymmetric Majorana covariance matrix

\[\Gamma_{ab} = \frac{i}{2}\langle [\gamma_a, \gamma_b] \rangle \, ,\]

where $\gamma_1, \ldots, \gamma_{2L}$ are the Majorana operators of the fermionic modes,

\[\gamma_{2i-1} = c_i + c_i^\dagger \, , \qquad \gamma_{2i} = -i\,(c_i^\dagger - c_i) \, ,\]

for mode i (with a Jordan–Wigner string on lower-indexed modes) — the convention of [1]. For $a \neq b$ this is simply $\Gamma_{ab} = \langle i\gamma_a\gamma_b\rangle$, and the occupation of mode i is $c_i^\dagger c_i = (1 - i\gamma_{2i-1}\gamma_{2i})/2$, so an unoccupied mode has $\Gamma_{2i-1,2i} = +1$ (see Conventions below). A pure Gaussian state satisfies the invariant

\[\Gamma^2 = -I\]

exactly (the eigenvalues of $i\Gamma$ are all $\pm 1$). Every gate and measurement in this backend updates $\Gamma$ while preserving this invariant to machine precision; a re-purification step (purify! — eigen-clamp $i\Gamma$'s spectrum back to $\pm 1$) fires automatically whenever floating-point drift exceeds state.backend.purify_tol (default 1e-10).

Supported Gates

Any gate outside this set — Hadamard, CNOT, HaarRandom, RandomClifford, SWAP, PhaseGate, PauliX, PauliY, PauliZ, CZ, Projection, SpinSectorProjection, and Reset — has no covariance-matrix representation and raises an informative error rather than being silently approximated:

using QuantumCircuitsMPS

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

ArgumentError: Gaussian backend only supports fermionic Gaussian operations (GaussianHaar, Measure(:Z), BondParity). Received: CNOT. Please switch to backend=:mps or backend=:statevector for non-Gaussian gates.

PauliX and Reset deserve a word, since they are available on the Clifford backend. The only occupation flip a covariance matrix can represent is the reflection of a single Majorana, i.e. conjugation by $\gamma_{2i}$. That operator is parity-odd: it changes the fermion parity of the state, which no closed fermionic system can do (fermion-parity superselection). The Jordan–Wigner image of the qubit $X_i$ fares no better — it is a product of $2i-1$ Majoranas, parity-odd as well. So there is no physically meaningful "flip" for a fermionic mode, and Reset (measure, then flip if occupied) goes with it. Reset is rejected by a dedicated execute! override before its Born draw, so the covariance matrix and the :born_measurement stream are untouched. To put a mode into a definite occupation, use Measure(:Z) (random outcome) or prepare the pattern with initialize!(state, ProductState(bitstring=...)).

Measure(:Z) does not get its own _apply_single!/execute! override on the fermionic-mode granularity: it flows through the generic engine, which calls this backend's born_probability + _measure_single_site!. BondParity gets a dedicated execute! override (see Backend Interface Contract) because it measures two sites at once — outside the single-site born_probability contract.

Majorana Chain (site_type="Majorana")

A second site granularity turns each site into one Majorana mode rather than one fermionic mode (2 Majoranas). This models the class-DIII monitored Majorana-chain circuit of Pan, Shapourian, Jian et al. directly, with no new gate types: GaussianHaar on 2 Majorana sites is $\exp(\theta\gamma_a\gamma_b)$, $\theta \sim U[0,2\pi)$ (Haar on $SO(2) \cong U(1)$ is exactly the uniform-angle rotation), and BondParity on 2 Majorana sites is the $i\gamma_i\gamma_{i+1}$ parity measurement — the staggered odd/even-link protocol maps directly onto Bricklayer(:odd)/Bricklayer(:even).

using QuantumCircuitsMPS

mstate = SimulationState(L=8, bc=:periodic, backend=:gaussian, site_type="Majorana",
    rng=RNGRegistry(gates_spacetime=1, gates_realization=2, born_measurement=3, state_init=4))
initialize!(mstate, ProductState(binary_int=0))
apply!(mstate, GaussianHaar(), Bricklayer(:odd))    # odd links: exp(θγγ), θ~U[0,2π)
apply!(mstate, BondParity(), Bricklayer(:even))     # even links: iγγ parity measurement
size(mstate.backend.corr)   # (8, 8) — Γ is L×L, not 2L×2L, on the Majorana chain

Key facts:

• site_type="Majorana" requires even L (a pure Gaussian state needs an even number of Majoranas): odd L throws ArgumentError at construction.
• $\Gamma$ is $L \times L$ (one Majorana per site) instead of $2L \times 2L$.
• ProductState bit patterns have length L÷2: bit k sets the parity sign of the consecutive Majorana pair $(\gamma_{2k-1}, \gamma_{2k})$ (dimerized vacuum when all bits are 0).
• Rejected on the Majorana chain (informative ArgumentError, each naming "Majorana"): Measure(:Z), Magnetization, BornProbability. There is no single-Majorana occupation or $\langle Z \rangle$ — parity lives on a pair. (PauliX and Reset are rejected on both granularities, see above.)
• EntanglementEntropy, MutualInformation, and TripartiteMutualInformation work unchanged (the site→Majorana index mapping is the identity on this granularity, and arbitrary/wrapped site subsets are still supported).

Physics sanity check on the dimerized vacuum (ProductState(binary_int=0), all pairs unoccupied): even cuts see zero entanglement (the cut falls between dimerized pairs), odd cuts split a pair and see log(2)/2 nats (half a fermion's worth), and two dimerized-paired Majorana sites have MI = log(2):

using QuantumCircuitsMPS

mstate = SimulationState(L=8, bc=:periodic, backend=:gaussian, site_type="Majorana",
    rng=RNGRegistry(gates_spacetime=1, gates_realization=2, born_measurement=3, state_init=4))
initialize!(mstate, ProductState(binary_int=0))
EntanglementEntropy(cut=2, base=ℯ)(mstate)   # ≈ 0.0
EntanglementEntropy(cut=1, base=ℯ)(mstate)   # ≈ log(2)/2 ≈ 0.3466
MutualInformation([1], [2])(mstate)          # ≈ log(2) ≈ 0.6931

Supported and Rejected Observables

Every rejection is an informative ArgumentError naming the observable and suggesting backend=:mps/:statevector — never a raw field-access crash:

using QuantumCircuitsMPS

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

ArgumentError: StringOrder is not supported on the Gaussian backend: its formula requires spin-1 Sz-expectation-string MPO/MPS contractions, which have no native fermionic-Gaussian (covariance-matrix) representation. Please use backend=:mps or backend=:statevector for StringOrder.

Region Entanglement Entropy

Like the state-vector and Clifford backends, EntanglementEntropy accepts an arbitrary physical-site region on the Gaussian backend — cut=c1:c2 or cut=[s1, s2, ...], including non-contiguous and PBC-wrapped selections — by mapping the region to its Majorana indices and reusing the same covariance-matrix subsystem_entropy formula the bipartition path uses:

using QuantumCircuitsMPS

state = SimulationState(L=6, bc=:periodic, backend=:gaussian,
    rng=RNGRegistry(gates_spacetime=1, gates_realization=2, born_measurement=3, state_init=4))
initialize!(state, ProductState(binary_int=0))
apply!(state, GaussianHaar(), Bricklayer(:odd))
apply!(state, GaussianHaar(), Bricklayer(:even))

EntanglementEntropy(cut=[6, 1])(state)   # PBC-wrapped region {6, 1}

As with the bipartition form, the region path supports real renyi_index (normalized to Float64, finite and > 0): both paths compute $S_n = \frac{1}{1-n}\sum_k \ln\!\big[((1+\lambda_k)/2)^n + ((1-\lambda_k)/2)^n\big]$ over the paired eigenvalues $\lambda_k$ of the region's covariance matrix, with an extra $\ln(2)/2$ contribution per unpaired zero mode on odd-sized Majorana-chain regions (n = 1 recovers the existing von Neumann formula bit-for-bit). Region sites are always physical sites, well-defined under both open and periodic boundary conditions (including wrapped regions like [L, 1]), and the region path is more expensive than the bipartition path: it builds an arbitrary Majorana-index submatrix rather than reusing the contiguous prefix structure of cut::Int.

Fermionic vs. spin subsystems

Every entropy on this backend is the entropy of a set of fermionic modes: the reduced state of a region $A$ is the restriction $\Gamma_A$ of the covariance matrix to the Majoranas of $A$, and its Rényi entropies come from the spectrum of $i\Gamma_A$. This is the quantity used in the free-fermion literature (the antipodal mutual information of [2], for instance), and it is the only one a covariance matrix can produce.

It is not always the number a qubit backend would return for the same sites after a Jordan–Wigner transformation. A spin region's reduced density matrix is a partial trace over the complementary spins; the fermionic one is defined by the operator algebra generated by the Majoranas of $A$, whose Jordan–Wigner strings reach through everything to the left of each mode. The two agree whenever $A$ or its complement is a single contiguous block of sites: every cut::Int bipartition, every range c1:c2, and PBC-wrapped regions such as [L, 1] (whose complement is contiguous). For a region that is non-contiguous and has a non-contiguous complement, the strings of the second block run through the gap — which belongs to the complement — and the two entropies differ in general. MutualInformation(A, B) inherits this through $S(A \cup B)$: for two non-adjacent blocks it is the fermionic mutual information, not the spin one.

A minimal example, $L = 4$: the state $\tfrac{1}{2}(1 + c_1^\dagger c_3^\dagger)(1 + c_2^\dagger c_4^\dagger)\lvert 0\rangle$ is a pure Gaussian state in which modes $\{1, 3\}$ form one pair and $\{2, 4\}$ another, so EntanglementEntropy(cut=[1, 3]) is exactly 0. Its Jordan–Wigner image is $\tfrac{1}{2}(\lvert 0000\rangle + \lvert 0101\rangle + \lvert 1010\rangle - \lvert 1111\rangle)$, whose spin sites $\{1, 3\}$ share one bit of entanglement with $\{2, 4\}$ (the Jordan–Wigner sign on the last term is what entangles them): the state-vector backend would report $\ln 2$ nats. Neither number is wrong; they answer different questions, and only the fermionic one is the entropy of a free-fermion subsystem. On the Majorana chain the point is sharper still: a region with an odd number of Majorana sites is not a tensor factor of any Hilbert space, and its entropy carries the unpaired-mode $\ln(2)/2$ described above.

Example: Reproducing the Class-DIII Phase Diagram

examples/gaussian_example.ipynb reproduces Fig. 1b of Pan et al., "Topological Modes in Monitored Quantum Dynamics": the staggered class-DIII monitored Majorana chain, mutual information vs. measurement probability p, at demo system sizes finishing in a few minutes. It uses the Majorana-chain granularity documented above with Bricklayer(:odd)/Bricklayer(:even) staggering.

Initialization

Three ways to prepare a Gaussian state, all via initialize!:

using QuantumCircuitsMPS, LinearAlgebra

mkstate() = SimulationState(L=4, bc=:open, backend=:gaussian,
    rng=RNGRegistry(gates_spacetime=1, gates_realization=2, born_measurement=3, state_init=4))

# 1. Vacuum: all modes unoccupied
s1 = mkstate(); initialize!(s1, ProductState(binary_int=0))
s1.backend.corr[1, 2]   # +1.0 (unoccupied convention)

# 2. Arbitrary occupation pattern (site 1 = MSB, same convention as every other backend)
s2 = mkstate(); initialize!(s2, ProductState(bitstring="0101"))   # sites 2, 4 occupied
s2.backend.corr[3, 4]   # -1.0 (occupied)

# 3. Haar-random pure Gaussian state
s3 = mkstate(); initialize!(s3, RandomGaussianState())
norm(s3.backend.corr * s3.backend.corr + I) < 1e-10   # true — purity holds

On the Majorana chain (site_type="Majorana"), ProductState bit patterns have length L÷2 (one bit per dimerized Majorana pair — see the Majorana Chain section above) and RandomGaussianState draws O ∈ SO(L) instead of SO(2L).

Conventions

• Covariance and Majorana definitions: $\Gamma_{ab} = \frac{i}{2}\langle[\gamma_a,\gamma_b]\rangle$ with $\gamma_{2i-1} = c_i + c_i^\dagger$ and $\gamma_{2i} = -i(c_i^\dagger - c_i)$ — the convention of [1]. Every $\Gamma$ in this package (the stored state.backend.corr, the docstring formulas, and the exact-diagonalization test oracle) uses it. Off-diagonal elements are $\Gamma_{ab} = \langle i\gamma_a\gamma_b\rangle$, so a projective measurement of the parity $i\gamma_a\gamma_b$ with eigenvalue $\pm 1$ leaves $\Gamma_{ab} = \pm 1$.
• Mode ↔ Majorana mapping: fermionic mode i (1-indexed) ↔ Majorana indices $(2i-1, 2i)$. On the Majorana chain, site i IS Majorana i directly.
• Occupation sign: $\Gamma_{2i-1,2i} = \langle i\gamma_{2i-1}\gamma_{2i}\rangle = 1 - 2\langle c_i^\dagger c_i\rangle$, so $\Gamma_{2i-1,2i} = +1$ ↔ mode i unoccupied; $\Gamma_{2i-1,2i} = -1$ ↔ occupied — i.e. $\langle c_i^\dagger c_i\rangle = (1 - \Gamma_{2i-1,2i})/2$. Verified empirically against the Python GTN reference implementation's get_C_f.
• Measurement outcome: outcome = 0 ↔ unoccupied / parity $i\gamma_a\gamma_b = +1$ ($\Gamma_{ab} = +1$ after the collapse); outcome = 1 ↔ occupied / parity $-1$ ($\Gamma_{ab} = -1$). ProductState bit 1 ↔ occupied, matching this convention.
• Exact-Haar sampler note: haar_orthogonal draws Haar-random SO(n) matrices exactly via QR decomposition of a Ginibre matrix (sign-fixed, det-corrected) — this is a deliberate departure from the Python GTN reference implementation, whose get_O uses expm of a random skew-symmetric matrix (an approximation the GTN codebase's own notes flag as not proven exactly Haar). All Gaussian-backend randomness (GaussianHaar, RandomGaussianState) uses the exact sampler.

Cross-Backend RNG Reproducibility

Like the other three backends, every measurement (Measure(:Z) and BondParity()) consumes exactly one :born_measurement draw per measured site/bond — even for a deterministic outcome, where the draw is discarded (the REDUNDANT-DRAW CONTRACT; see the Backend Interface Contract). Random gate content (GaussianHaar) draws from :gates_realization; random initial states (RandomGaussianState) draw from :state_init. Same seeds ⇒ bitwise-identical $\Gamma$ and identical measurement records, within the Gaussian backend:

using QuantumCircuitsMPS

mkrng(k) = RNGRegistry(gates_spacetime=k, gates_realization=k+10, born_measurement=k+20, state_init=k+30)

sA = SimulationState(L=6, bc=:open, backend=:gaussian, rng=mkrng(7))
initialize!(sA, ProductState(binary_int=0))
apply!(sA, GaussianHaar(), Bricklayer(:odd))

sB = SimulationState(L=6, bc=:open, backend=:gaussian, rng=mkrng(7))
initialize!(sB, ProductState(binary_int=0))
apply!(sB, GaussianHaar(), Bricklayer(:odd))

sA.backend.corr == sB.backend.corr   # true — bitwise identical

References

• The fermionic-linear-optics / Majorana-covariance-matrix formalism this backend implements: [1].
• The Haar-random SO(2n) Gaussian-unitary ensemble and the general contraction formula: [3]; conventions follow [2].
• The class-DIII staggered monitored Majorana-chain protocol reproduced by examples/gaussian_example.ipynb: [2].