Fermionic Gaussian States

• Fermionic Gaussian states are quantum states fully described by their two-point correlation (covariance) matrices, making them central to free-fermion models.
• They employ covariance-matrix techniques and Wick’s theorem to reconstruct higher-order correlators, ensuring analytical and numerical tractability even in volume-law entangled regimes.
• Their framework underpins scalable numerical methods, tensor-network representations, and rigorous entanglement and channel analyses in complex many-body systems.

Fermionic Gaussian states are quantum states of fermionic many-body systems whose full correlation structure is fixed by two-point functions. Equivalently, they are quasifree states, thermal or pure states of quadratic fermionic Hamiltonians, and—within a geometric formulation—states specified by a linear complex structure on the classical fermionic phase space (Surace et al., 2021, Hackl et al., 2020, Carollo et al., 2019). They constitute the natural state class of free-fermion models, matchgate circuits, and fermionic linear optics, and they remain analytically and numerically tractable through covariance-matrix techniques and Wick’s theorem even in regimes with volume-law entanglement and nontrivial multipartite structure (Collura et al., 2024).

1. Algebraic definition and covariance-matrix formalism

An $n$-mode fermionic system is generated by creation and annihilation operators obeying the canonical anticommutation relations, or equivalently by $2n$ Majorana operators such as

$$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$

with

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$

Across the literature summarized here, a fermionic Gaussian state is written in equivalent exponential forms,

$$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$

where $h$ and $\Omega$ are real antisymmetric matrices and $K$ is quadratic in the fermionic fields (Negari et al., 25 Aug 2025, Carollo et al., 2019, Ares et al., 17 Mar 2026).

The central object is the covariance matrix. In one common convention it is

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$

while in another it is

$$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$

These matrices are real or purely imaginary antisymmetric matrices depending on convention, and physicality is encoded by conditions such as

$2n$0

Purity is characterized by saturation,

$2n$1

or, equivalently, by covariance eigenvalues of unit modulus (Negari et al., 25 Aug 2025, Surace et al., 2021, Spee et al., 2017).

Gaussianity is equivalent to Wick/Pfaffian closure of higher moments. For even correlators one has Pfaffian factorization, and all odd moments vanish for parity-even physical states. This is the operative reason why covariance matrices suffice: once the second moments are known, all higher-order correlators are algorithmically reconstructible (Carollo et al., 2019, Spee et al., 2017). A crucial structural constraint is parity superselection: physical fermionic states commute with the parity operator,

$2n$2

which makes separability, local operations, and channel theory more subtle than in qubit or bosonic settings (Negari et al., 25 Aug 2025).

This covariance-matrix closure explains the dual status of fermionic Gaussian states. They are the free sector of fermionic many-body theory and of matchgate computation, yet they need not be weakly entangled: they can remain classically simulable even in regimes with extensive entanglement (Collura et al., 2024).

2. Canonical forms, complex structures, and bipartite normal forms

Any covariance matrix can be brought by an orthogonal transformation to canonical $2n$3 blocks. In one standard form,

$2n$4

with canonical values $2n$5. For mixed states these values encode occupation probabilities and entropies; for pure states they lie at the extremal boundary (Zhang et al., 2022).

A basis-independent formulation is provided by the Kähler description. There the state is encoded by a triple $2n$6, where $2n$7 is the canonical fermionic metric, $2n$8 is the state-dependent antisymmetric form, and $2n$9 is a linear complex structure. Pure fermionic Gaussian states satisfy

$$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$0

and are uniquely characterized by $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$1 up to phase; mixed Gaussian states correspond to $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$2, with eigenvalues $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$3, $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$4 (Hackl et al., 2020). This formulation recasts covariance-matrix manipulations as algebra on the classical phase space and places fermionic Gaussian states on the symmetric space $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$5 (Hackl et al., 2020).

For a bipartition $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$6 with subsystem dimensions $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$7, the reduced covariance matrix on $$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$8 can be brought to the block-diagonal form

$$\gamma_{2j-1}=\hat c_j+\hat c_j^\dagger,\qquad \gamma_{2j}=i(\hat c_j-\hat c_j^\dagger),$$9

with $\{\gamma_p,\gamma_q\}=2\delta_{pq}.$0 and entanglement-spectrum variables

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$1

These $\{\gamma_p,\gamma_q\}=2\delta_{pq}.$2 control all spectral entanglement observables of the bipartition (Huang et al., 2021).

Up to local Gaussian unitaries, a pure bipartite fermionic Gaussian state may also be written in a canonical BCS-like product form,

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$3

which is the fermionic analogue of a Schmidt decomposition (Iannotti et al., 29 Apr 2026). A closely related exact Schmidt decomposition in correlation-matrix language underlies efficient conversion of Gaussian states to matrix product states (Liu et al., 2024).

These canonical forms show that bipartite structure is controlled by a single-particle spectral problem. Entropy, negativity in special cases, non-local magic, and several random-matrix statistics are all functions of the same reduced covariance spectrum.

3. Entanglement observables and random-matrix statistics

For bipartite pure fermionic Gaussian states, the von Neumann entropy of the reduced state is a linear spectral statistic,

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$4

The capacity of entanglement is likewise a spectral sum,

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$5

In the random-state setting, the variables $\{\gamma_p,\gamma_q\}=2\delta_{pq}.$6 are distributed according to the non-standard Jacobi unitary ensemble

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$7

called the fermionic Gaussian ensemble (Huang et al., 2021).

For this ensemble, exact formulas are known for the mean entropy and, now, for its variance. In the unconstrained case the exact variance is

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$8

resolving the conjecture posed earlier for arbitrary subsystem sizes (Huang et al., 2022). In the large-dimension limit this reproduces

$\{\gamma_p,\gamma_q\}=2\delta_{pq}.$9

while the fixed-particle-number ensemble admits an exact variance formula with coefficients $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$0 depending on $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$1 (Huang et al., 2022).

The same random-matrix framework yields second-order information beyond entropy. For equal subsystem dimensions, the average entanglement capacity grows linearly,

$$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$2

whereas the entropy variance saturates to a ratio-dependent limit rather than growing extensively (Huang et al., 2021).

A further statistical phenomenon is Gaussianity of the standardized entropy

$$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$3

Numerically, $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$4 approaches a standard normal law at fixed aspect ratio $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$5, and the asymptotic conjecture

$$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$6

places fermionic Gaussian entropy within the class of linear spectral statistics displaying central-limit behavior (Huang et al., 2021). This suggests that entanglement fluctuations in random fermionic Gaussian states exhibit the same type of spectral universality familiar from other random-matrix ensembles.

4. Partial transpose, negativity, and local equivalence

The mixed-state entanglement theory of fermionic Gaussian states differs sharply from the bosonic case because the fermionic partial transpose does not preserve Gaussianity. For a bipartition $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$7, however, the partial transpose of a fermionic Gaussian state admits a rigid decomposition: $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$8 where $$\rho=C\,\exp\!\left(\frac{i}{2}\gamma^T h\gamma\right), \qquad \rho=\frac{e^{-\frac{i}{4}\bm{\omega}^T\Omega\,\bm{\omega}}}{Z}, \qquad \rho=\frac{1}{Z}e^{-K},$$9 are Gaussian operators determined uniquely by modified covariance matrices

$h$0

This is the basic structural fact behind fermionic negativity calculations (Eisler et al., 2015).

For pure bipartite states, the logarithmic negativity becomes exact and reduces to a spectral expression. If the reduced covariance matrix of $h$1 has eigenvalues $h$2, then

$h$3

so that $h$4 in the pure-state case (Eisler et al., 2015). For mixed states, exact covariance-matrix formulas are generally unavailable, and rigorous upper and lower bounds are obtained instead through pinching constructions, semidefinite programming, and Gaussian-operator product formulas (Eisert et al., 2016).

Local equivalence and separability are likewise constrained. For $h$5-mode $h$6-partite Gaussian fermionic states, the physically appropriate Gaussian separability notion is the one in which the covariance matrix is a direct sum of local covariance matrices. Under Gaussian local unitaries there is a unique mixed-state standard form, and two Gaussian fermionic states are GLU-equivalent if and only if their standard forms coincide (Spee et al., 2017). Deterministic Gaussian LOCC is maximally rigid: for pure fully entangled Gaussian fermionic states, any Gaussian LOCC transformation can already be implemented by Gaussian local unitaries, so there are no non-trivial GLOCC transformations (Spee et al., 2017).

This combination of results establishes an important point of principle. The covariance matrix still organizes the theory, but entanglement monotones and local convertibility are more singular than in bosonic Gaussian theory because partial transpose, separability, and parity constraints do not align with a single Gaussian manifold.

5. Extendibility, de Finetti structure, and Gaussian channels

Extendibility provides a second route to the correlation structure of fermionic Gaussian states. For a bipartite covariance matrix

$h$7

a central theorem states that a fermionic Gaussian state is $h$8-extendible if and only if it admits a fermionic Gaussian extension (Negari et al., 25 Aug 2025). The existence of such an extension is equivalent to a covariance-matrix feasibility condition involving auxiliary antisymmetric matrices $h$9, which yields a semidefinite program whose size scales linearly with the number of modes (Negari et al., 25 Aug 2025).

This has several structural consequences. First,

$\Omega$0

so large extendibility suppresses cross-correlations (Negari et al., 25 Aug 2025). Second, in the Gaussian sector, arbitrary extendibility from either side forces $\Omega$1, hence block-diagonal covariance and factorization. The same work identifies a “classical–quantum binding” phenomenon: $\Omega$2 so Gaussian separability, arbitrary extendibility, and product structure coincide (Negari et al., 25 Aug 2025). A finite de Finetti-type theorem then bounds the trace distance to separable states by

$\Omega$3

with companion bounds for relative entropy of entanglement and squashed entanglement (Negari et al., 25 Aug 2025).

A related graph-theoretic extendibility program applies beyond the Gaussian sector to interacting fermionic systems. There, parity projection and monogamy of odd fermionic correlations lead to non-symmetric de Finetti statements and graph-dependent approximation guarantees for fermionic Gaussian product states as ground-state ansätze (Krumnow et al., 2024). This suggests that Gaussian product structure is not merely a convenient variational class but, on sufficiently high-degree interaction graphs, an approximation class with rigorous error control.

Fermionic Gaussian channels are naturally described in covariance language. A channel specified by matrices $\Omega$4 acts as

$\Omega$5

Its Choi state is Gaussian, anti-degradability becomes an SDP feasibility problem, and entanglement-breaking channels collapse to replacement channels: $\Omega$6 Hence there are no nontrivial entanglement-breaking fermionic Gaussian channels (Negari et al., 25 Aug 2025). In a different physical setting, acceleration of localized Dirac modes is itself described by a fermionic Gaussian channel

$\Omega$7

through which vacuum entanglement is found to increase with acceleration while Bell-state entanglement decreases (Richter et al., 2017).

6. Numerical methods, tensor-network representations, and metrological formulas

The practical power of fermionic Gaussian states lies in exact polynomial-time linear algebra. A quadratic Hamiltonian can be written in Majorana form as

$\Omega$8

with $\Omega$9 real skew-symmetric, and diagonalized by an orthogonal transformation into canonical $K$0 blocks. For Gaussian states, the Dirac correlation matrix satisfies the compact thermal relation

$K$1

and real-time evolution remains closed: $K$2 Reduced states are obtained by restriction of the correlation matrix, and products of Gaussian states remain Gaussian under explicit covariance update rules (Surace et al., 2021).

Several recent algorithms exploit this structure. One is a truncation method for the trace distance between two fermionic Gaussian states. It selects effective modes either by entropy contribution, by large canonical values of $K$3, or by a mixed strategy, and computes the exact trace distance in the reduced subspace. The method is effective for small-entropy nearly commuting states and for nearly orthogonal states, extending subsystem trace-distance calculations in Ising and XX chains from subsystems of roughly ten sites to several hundred sites (Zhang et al., 2022).

Another algorithmic development concerns amplitudes in non-computational bases. Gaussian pure states of the form

$K$4

admit exact Pfaffian amplitudes in the computational basis, and, after a domain-wall mapping, in $K$5, $K$6, and more general $K$7 bases through transformed antisymmetric matrices $K$8 (Tarighi et al., 2024). This makes formation probabilities and Shannon-type observables accessible in large systems.

A complementary tensor-network direction converts fermionic Gaussian states to matrix product states. Using exact Gaussian Schmidt decompositions and virtual-mode encodings, local tensors are constructed recursively from subsystem correlation matrices, and a mode-decimation scheme reduces local tensor generation from $K$9 to

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$0

with $$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$1 the retained bond dimension (Liu et al., 2024). In infinite-cylinder geometries this facilitates extraction of minimally entangled states, entanglement spectra, and anyon sectors for topological parton states (Liu et al., 2024).

Metrological applications are likewise covariance based. The symmetric logarithmic derivative of a fermionic Gaussian state has the quadratic form

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$2

where $$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$3 is determined by the discrete Lyapunov equation

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$4

This gives a direct route to the quantum Fisher information and multiparameter quantum Fisher information matrix for thermal and non-equilibrium fermionic Gaussian probes (Carollo et al., 2019).

7. Non-Gaussianity, magic, and emerging directions

Because fermionic Gaussian states are the free, efficiently simulable sector, several recent developments define resources by distance from the Gaussian manifold. One route introduces fermionic convolution,

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$5

with a cumulant rule that suppresses all higher-order cumulants under repeated self-convolution. The associated Gaussification map $$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$6 converges to the Gaussian state with the same covariance matrix and satisfies

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$7

For pure even states this leads to a three-copy Gaussianity test and to the non-Gaussian entropy

$$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$8

which vanishes exactly on pure Gaussian states (Lyu et al., 2024).

A second line relates non-Gaussianity to particle-number asymmetry. If $$M_{pq}=\operatorname{Tr}(i\gamma_q\gamma_p\,\rho)\quad (p\neq q),\qquad M_{pp}=0,$$9 is the Gaussianification of $$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$0, then the relative entropy of non-Gaussianity is

$$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$1

Using concentration bounds for the particle-number distribution of mixed fermionic Gaussian states, one obtains lower bounds on $$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$2 in terms of the Shannon entropy of the charge distribution. For maximal asymmetry, the bound scales linearly with system size, which shows that large particle-number asymmetry can force extensive non-Gaussianity (Ares et al., 17 Mar 2026).

Magic theory produces a more surprising picture. Perfect sampling from the Majorana characteristic distribution of a fermionic Gaussian state reduces to exact sampling from a determinantal point process with kernel given by the covariance matrix, at cost

$$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$3

per sample (Collura et al., 2024). This enables computation of stabilizer Rényi entropies for systems with hundreds of qubits and reveals that random fermionic Gaussian states have filtered stabilizer Rényi entropies close to the maximal value $$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$4, with logarithmic corrections, despite remaining classically simulable (Collura et al., 2024).

The non-local component of magic can also be extracted exactly from reduced covariance spectra alone. For a bipartition with reduced covariance eigenvalues $$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$5, the fermionic non-local magic is given by

$$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$6

For $$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$7,

$$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$8

This yields an exact Page-like curve for random fermionic Gaussian states, logarithmic scaling at the XY critical point, and a quasiparticle picture for non-local-magic growth after quenches (Iannotti et al., 29 Apr 2026).

Two additional directions emphasize the breadth of the framework. First, projective Bell measurements on half the rungs of two identical half-filled Gaussian layers produce a postselected state that factorizes universally into Bell pairs,

$$\Gamma_{jk}:=\frac12\,\operatorname{Tr}\!\big(\rho[\omega_j,\omega_k]\big).$$9

independently of the microscopic details of the initial Gaussian state (Fang et al., 17 Dec 2025). Second, sufficiently high-temperature Gibbs states of bounded-degree local fermionic Hamiltonians become convex mixtures of fermionic Gaussian states for

$2n$00

and admit efficient classical sampling in the slightly stronger regime

$2n$01

which places a sharp structural boundary on thermal complexity in local fermionic systems (Ramkumar et al., 14 May 2025).

Taken together, these developments present fermionic Gaussian states as more than a solvable subclass. They are simultaneously a structural language for free fermions, a boundary object for resource theories of non-Gaussianity and magic, a testing ground for mixed-state entanglement theory, and a scalable computational interface between many-body physics, random-matrix theory, metrology, and tensor networks.