API Reference

This page contains the complete API reference for the Free Fermion Library.

Core Module (ff_lib)

Free Fermion Library - Basic Functions For Handling Fermions

This module contains the core quantum physics and linear algebra functions for free fermion systems, including Jordan-Wigner transformations, symplectic diagonalization, Gaussian state constructions, and correlation matrix computations.

Key functionality:

Jordan-Wigner fermionic operators (Dirac and Majorana)
Symplectic free-fermion diagonalization
Gaussian state generation and manipulation
Fermionic correlation matrix computations
Block matrix construction for coefficient matrices
Wick’s theorem implementation

ff.ff_lib.permutation_to_matrix(permutation)[source]

Transform a permutation (list) into a permutation matrix using NumPy.

Parameters:: permutation – A list representing a permutation
Returns:: A NumPy array representing the permutation matrix. Returns None if the input is invalid.

ff.ff_lib.pauli_matrices()[source]: Define the Pauli matrices as numpy arrays.

ff.ff_lib.generate_pauli_group(n, verbose=False)[source]

Generate all 4^n elements of the Pauli group for n qubits.

This function generates all possible tensor products of the Pauli matrices (I, X, Y, Z) for n qubits, creating the complete Pauli group of order 4^n. Each element is a 2^n x 2^n matrix representing a Pauli operator acting on the n-qubit Hilbert space.

The Pauli group is fundamental in quantum error correction, quantum computing, and stabilizer formalism. It consists of all n-fold tensor products of single-qubit Pauli operators.

Parameters:

n (int) – The number of qubits. Must be a positive integer.
verbose (bool) – If True, print the Pauli string name for each operator as it’s generated (default: False).

Returns:

A list of 4^n numpy arrays, each of shape (2^n, 2^n),: representing all Pauli operators for n qubits. The operators are ordered lexicographically by their tensor product structure.

Return type:

list

Raises:

TypeError – If n is not an integer.
ValueError – If n is not positive.

Examples

>>> # Single qubit case (n=1)
>>> pauli_1 = generate_pauli_group(1)
>>> len(pauli_1)
4
>>> pauli_1[0]  # Identity
array([[1, 0],
       [0, 1]])
>>> pauli_1[1]  # Pauli-X
array([[0, 1],
       [1, 0]])

>>> # Two qubit case (n=2) with verbose output
>>> pauli_2 = generate_pauli_group(2, verbose=True)
II
IX
IY
IZ
XI
XX
XY
XZ
YI
YX
YY
YZ
ZI
ZX
ZY
ZZ
>>> len(pauli_2)
16
>>> pauli_2[0].shape
(4, 4)

Notes

The function uses the existing pauli_matrices() function for
consistency
Memory usage scales as O(4^n * 4^n) = O(16^n), so use with caution
for large n
For n > 4, consider using sparse representations or symbolic methods
The ordering follows itertools.product convention: rightmost index
varies fastest
Verbose output shows Pauli strings in standard notation (I, X, Y, Z)

ff.ff_lib.jordan_wigner_lowering(n_sites)[source]

Define annihilation operators using the Jordan-Wigner transform.

Parameters:: n_sites – The number of lattice sites
Returns:: A list of the annihilation operators

ff.ff_lib.jordan_wigner_alphas(n_sites)[source]

Define raising and lowering operators using the Jordan-Wigner transform. With alphas = raising + lowering

Parameters:: n_sites – The number of lattice sites
Returns:: A list of the [raising, lowering] operators

ff.ff_lib.jordan_wigner_majoranas(n_sites)[source]

Generate a set of Majorana operators under Jordan-Wigner.

\[ \begin{align}\begin{aligned}\gamma_j = (a_j + a_j^\dagger)/\sqrt{2}, \qquad j < n_{sites}\\\gamma_k = i(a_k - a_k^\dagger)/\sqrt{2},\qquad k < 2n_{sites}\end{aligned}\end{align} \]

Parameters:: n_sites – The number of lattice sites
Returns:: A list of the Majorana operators

ff.ff_lib.rotate_operators(C, alphas, Left=False, verbose=False)[source]

Transform set of operators using a matrix C per:

\[\beta_j = \sum_i C_{ji} \alpha_i\]

Use Left=True if instead you want

\[\beta_j = \sum_i \alpha_i C_{ij}^*\]

Parameters:

C – A [len(alpha) x len(alpha)] coefficient matrix
alphas – A list of operators (NumPy matrices)
Left – Do the transformation left-handed per above (default=False)
verbose – Print more information (default=False)

Returns:

The list of rotated operators

ff.ff_lib.build_V(n_sites, A, Z=None)[source]

Build the generator V matrix of dimension 2*N for N sites.

\[\begin{split}V=\begin{bmatrix} -Z^* & A \\ -A^* & Z \end{bmatrix}\end{split}\]

Parameters:

n_sites – The number of sites
A – The A coefficient matrix, [N x N]
Z – The Z coefficient matrix, [N x N] (optional)

Returns:

A 2*N x 2*N numpy array representing a generator matrix

ff.ff_lib.build_H(n_sites, A, B=None)[source]

Build the H block coefficient matrix of dimension 2*N for N sites.

\[\begin{split}H=\begin{bmatrix} -A^* & B \\ -B^* & A \end{bmatrix}\end{split}\]

If A = A.conj().T and B = -B.T then the output is compatible with alphas = jordan_wigner_alphas(n_sites)

Parameters:

n_sites – The number of sites, N
A – The A coefficient matrix, [N x N]
B – The B coefficient matrix, [N x N] (optional)

Returns:

A [2*N x 2*N] numpy array representing the H matrix

ff.ff_lib.random_FF_rotation(n_sites, seed=None, returnH=False, returnOrb=False)[source]

Generate a random free fermion rotation matrix

Parameters:

n_sites – The number of sites
seed – Random seed for reproducibility (default: None)
returnH – If True, return the generator matrix instead (default: False)
returnOrb – If True, return the orbital rotation matrix C (default: False)

Returns:

A random free fermion rotation matrix of dimension 2*N for N sites

ff.ff_lib.random_FF_state(n_sites, fixedN=False, seed=None, returnH=False, pure=False)[source]

Generate Haar random free fermion state using symplectic rotations.

This function generates free-fermion states that are uniformly distributed over the space of free-fermion states using Haar random symplectic transformations.

Parameters:

n_sites – The number of sites
fixedN – If True, generator commutes with N_op (default: False)
seed – Random seed for reproducibility (optional)
returnH – If True, return the generator matrix along with rho (default: False)
pure – If True, return a pure state (default: False)

Returns:

A normalized Haar random free fermion state, rho. If returnH is True, also returns the generator matrix H.

ff.ff_lib.random_H_generator(n_sites, fixedN=False, seed=None)[source]

Generate a random generator matrix for free fermions. Entries are drawn from a standard normal distribution for both real and imaginary parts.

Parameters:

n_sites – The number of sites
fixedN – If True, generator will preserve N (default: False)
seed – Random seed for reproducibility (optional)

Returns:

A random generator matrix of dimension 2*N for N sites

ff.ff_lib.kitaev_chain(n_sites, mu, t, delta)[source]: Create Kitaev chain Hamiltonian.

ff.ff_lib.build_Omega(N)[source]

Build the Omega matrix of dimension 2*N for N sites.

\[\begin{split}\Omega = \frac{1}{\sqrt{2}} \begin{bmatrix} \mathbf{1} & \mathbf{1} \\ i \mathbf{1} & -i \mathbf{1} \end{bmatrix}\end{split}\]

Parameters:: N – The number of sites
Returns:: A 2*N x 2*N numpy array representing the Omega matrix

ff.ff_lib.build_reordering_xx_to_xp(n_sites)[source]

A permutation matrix that transforms ordering [x_1, x_2… p_1, p_2…] to [x_1, p_1, x_2, p_2, …]

Parameters:: n_sites – Number of sites
Returns:: Permutation matrix for reordering

ff.ff_lib.build_K(n_sites)[source]: Build the symplectic transform from standard to canonical eigenstruct

ff.ff_lib.is_symp(U)[source]

Function checks if U is symplectic

U is symplectic if U S U^T = S

A known canonical form is:

\[\begin{split}U = \begin{bmatrix} s & t^* \\ t & s^* \end{bmatrix}\end{split}\]

where s and t are square matrices of dimension N.

Parameters:: U – Matrix to check for symplectic property
Returns:: True if U is symplectic, False otherwise

ff.ff_lib.check_canonical_form(L)[source]

Check if a matrix is in standard block diagonal form.

Parameters:: L – A NumPy array representing the matrix
Returns:: True if the matrix is in standard block diagonal form, False otherwise

ff.ff_lib.generate_gaussian_state(n_sites, H, alphas=None)[source]

Generate a Gaussian state using Hamiltonian H.

Parameters:

n_sites – The number of orbitals
H – Quadratic parent Hamiltonian ([2n x 2n] or [2**n x 2**n])
alphas – The list of 2n fermionic operators (optional)

Returns:

The [2**n x 2**n] free fermion state

ff.ff_lib.build_op(n_sites, R, alphas, verbose=None, direct=False)[source]

Build the N-body lift of a quadratic coefficient matrix

R̂ = Σ_ij R_ij α_i† α_j

Parameters:

n_sites – The number of sites
R – Coefficient matrix [2N x 2N]
alphas – A list of 2N fermionic operators
verbose – Increases output if set to True or an integer (default: None)
direct – Do R_ij α_i α_j (no adjoint) (default: False)

Returns:

The N-body Hamiltonian operator as a matrix

ff.ff_lib.compute_cov_matrix(rho, n_sites=None, alphas=None)[source]

Calculate the fermionic covariance matrix of rho

K_ij = Tr[rho α_i α_j] - Tr[rho α_j α_i]

Parameters:

rho – should be a [2**n_sites x 2**n_sites] matrix
n_sites – number of sites
alphas – a set of 2 n_sites of fermionic operators (optional)

Returns:

An [2 n_sites x 2 n_sites] covariance matrix

ff.ff_lib.correlation_matrix(rho)[source]

Calculates the following two-point correlation matrix

\[\Gamma = \langle \vec \alpha \vec \alpha ^t \rangle\]

\[\Gamma_{ij} = Tr[\rho \alpha_i \alpha_j]\]

for JW fermionic operators in the [a^+ a] ordering

Parameters:: rho – should be a [2**n_sites x 2**n_sites] matrix
Returns:: An [2 n_sites x 2 n_sites] correlation matrix

ff.ff_lib.compute_2corr_matrix(rho, n_sites, alphas=None, conjugation=None)[source]

Calculate the two-point correlation matrix

\[\Gamma = \langle \vec \alpha \vec \alpha ^t \rangle\]

\[\Gamma_{ij} = Tr[\rho \alpha_i \alpha_j]\]

By changing the input conjugation parameter, this function also computes P = ⟨vec alpha vec α†⟩ (conjugation == T) η = ⟨vec alpha^dagger vec α^t⟩ (conjugation < 0)

By default the operators are not conjugated but if conjugation is True or positive:

P_ij = Tr[ρ α_i α_j†]

And if conjugation is negative:: η_ij = Tr[ρ α_i† α_j]

Parameters:

rho – should be a [2**n_sites x 2**n_sites] matrix
n_sites – number of sites
alphas – a set of 2 n_sites of fermionic operators (optional)
conjugation – Indicates if the operators should be conjugated (optional)

Returns:

An [2 n_sites x 2 n_sites] correlation matrix

ff.ff_lib.compute_algebra_S(alphas, verbose=False)[source]

Given operators α_i, compute the S matrix

α_i α_j = S_ij 𝟙 - α_j α_i

Parameters:

alphas – A list of 2N operators
verbose – Output more information (default: False)

Returns:

The algebraic S matrix

Return type:

S

ff.ff_lib.is_matchgate(M, verbose=False)[source]

Check the matchgate condition for 4x4 input matrices.

A 4x4 matrix satisfies the matchgate condition if the determinant of its inner 2x2 submatrix equals the determinant of its corner 2x2 submatrix.

Parameters:

M – 4x4 numpy array to test
verbose – Print detailed information (default: False)

Returns:

True if M satisfies the matchgate condition, False otherwise

ff.ff_lib.eigh_sp(H)[source]

Return the eigenvectors in symplectic form and the eigenvalues of a complex Hermitian (conjugate symmetric) FF coefficient array.

Returns two objects, a 1-D array containing the eigenvalues of H and a 2-D array containing a symplectic unitary that diagonalizes H.

This routine transforms to Majorana forms, gets the canonicalizing orthogonal matrix for skew-symmetric real matrices. This is transformed to symplectic form using

\[U_{symp} = \Omega^\dagger O \Omega\]

Then a symplectic transformation is constructed to transform from canonical eigenstructure

\[\begin{split}L_{canonical} = \oplus_k \begin{pmatrix} 0 & -\lambda_k\\ \lambda_k & 0 \end{pmatrix}\end{split}\]

to standard eigenstructure

\[\begin{split}L_{standard} = \begin{pmatrix} -\Sigma & 0\\ 0 & \Sigma \end{pmatrix}\end{split}\]

Parameters:

H –

(2N, 2N) array Should be of the form

[[-A.conj(), B],
[-B.conj(), A]]

with \(A = A^\dagger\) and \(B = -B^T\)

Returns:

A tuple with eigenvalues: (2N) ndarray eigenvectors: (2N,2N) ndarray

Ortho-normal eigenvectors arranged in the form \(\begin{pmatrix} s & t\\t^* & s^* \end{pmatrix}\)

ff.ff_lib.eigv_sp(V)[source]

Returns the left eigenvectors in symplectic form and eigenvalues in Y-form

\[\begin{split}L_{Y} = \begin{pmatrix} 0 & -\Sigma\\ \Sigma & 0 \end{pmatrix}\end{split}\]

of the coefficient matrix in V-form whereby \(\hat V = \hat V^\dagger = \sum V_{ij}\alpha_i \alpha_j\)

Returns two objects, a 1-D array containing the eigenvalues of V in Y-form and a 2-D array containing a symplectic unitary that diagonalizes V.

\[\begin{split}L_{standard} = \begin{pmatrix} -\Sigma & 0\\ 0 & \Sigma \end{pmatrix}\end{split}\]

Parameters:

V – (2N, 2N) array Should be of the form \(V = \begin{pmatrix}-B^* & A\\ -A^* & B\end{pmatrix}\) with A = A† and B = -B^T

Returns:

A tuple with eigenvalues: (2N) ndarray in \(L_Y\) form eigenvectors: (2N,2N) ndarray

Ortho-normal eigenvectors arranged in the form \(U = \begin{pmatrix} s & t\\ t^* & s^* \end{pmatrix}\) such that U.T @ V @ U = L_Y

ff.ff_lib.eigm_sp_can(G)[source]

Returns the orthogonal eigenvectors and the eigenvalues in direct sum canonical form for the Majorana coefficient matrix G.

G should be of form .. math:

G =  i\begin{bmatrix}Im(A+Z) & Re(A+Z)\\
Re(Z-A)&  Im(A-Z)
\end{bmatrix}

where Z = -Z^T and A = A^dagger.

Returns two objects, the eigenvalues of G in canonical form and the sympletic orthogonal matrix that transforms G to canonical eigenstructure form .. math:

L_{canonical} = \oplus_k
\begin{matrix} 0 & -\lambda_k\\
\lambda_k & 0 \end{matrix}

Parameters:

G (( 2N, 2N ) array properly formatted)

Returns:

A tuple with
eigenvalues ((2N) ndarray) –

Lc = oplus_k [[ 0 , -ev_k ],
[ ev_k , 0 ]]
eigenvectors ((2N,2N) ndarray) – Ortho-normal eigenvectors arranged in the form U = [[ s , t ]

[ t^* , s^* ]]

ff.ff_lib.eigm_sp(G)[source]

Returns the orthogonal eigenvectors and the eigenvalues in L_Y form for the Majorana coefficient matrix G. G should be of form

\[\begin{split}G = i\begin{bmatrix}Im(A+Z) & Re(A+Z)\\ Re(Z-A)& Im(A-Z) \end{bmatrix}\end{split}\]

where \(Z = -Z^T\) and \(A = A^\dagger\).

Returns two objects, the eigenvalues of G in canonical form and the sympletic orthogonal matrix that transforms G to the Y eigenstructure form

\[\begin{split}L_{Y} = \begin{pmatrix} 0 & -\Sigma\\ \Sigma & 0 \end{pmatrix}\end{split}\]

Parameters:

G (( 2N, 2N ) array properly formatted)

Returns:

A tuple with
eigenvalues (( 2N ) ndarray) – i*L in L_Y form
eigenvectors (( 2N, 2N ) ndarray) – Orthogonal matrix such that O^T O = Id

Combinatorics Module (ff_combinatorics)

Free Fermion Combinatorial Functions Module

Core combinatorial matrix functions including determinants, pfaffians, permanents, and hafnians.

ff.ff_combinatorics.sgn(permutation)[source]

Computes the sign of a permutation by counting inversions.

Parameters:: permutation – List or array representing a permutation
Returns:: 1 if even number of inversions, -1 if odd

ff.ff_combinatorics.pf(A)[source]

Computes pfaffian via combinatorial formula:

Given \(A\) of dimension \(N = 2 n_e\) we have the pfaffian as:

\[pf(A)= \sum_{\sigma\in S_N}^{N!} sgn(\sigma) W_\sigma(A)/(2^{n_e} n_e!)\]

Where the weight of matching \(\sigma\) is given by

\[W_\sigma(A) = \prod_{i=0}^{n_e} A_{\sigma(2i), \sigma(2i+1)}\]

Parameters:: A – Square matrix (numpy array or array-like)
Returns:: Pfaffian value (complex or real)

ff.ff_combinatorics.hf(A)[source]

Computes hafnian via combinatorial formula.

Given \(A\) of even dimension \(N = 2 n_e\) we have the hafnian as:

\[hf(A) = \sum_{\sigma\in S_N}^{N!} W_\sigma(A)/(2^{n_e} n_e!)\]

Where the weight of matching \(\sigma\) is given by

\[W_\sigma(A) = \prod_{i=0}^{n_e} A_{\sigma(2i), \sigma(2i+1)}\]

Parameters:: A – Square matrix (numpy array or array-like)
Returns:: Hafnian value (complex or real)

ff.ff_combinatorics.pt(A)[source]

Computes permanent via combinatorial formula.

Given \(A\) of dimension \(N\) we have the permanent as:

\[pt(A) = \sum_{\sigma\in S_N}^{N!} B_\sigma(A)\]

Where the bipartite matching of \(\sigma\) is given by

\[B_\sigma(A) = \prod_{i=1}^N A_{i,\sigma(i)}\]

Parameters:: A – Square matrix (numpy array or array-like)
Returns:: Permanent value (complex or real)

ff.ff_combinatorics.dt(A)[source]

Computes determinant via combinatorial formula.

Given \(A\) of dimension \(N\) we have the determinant as:

\[dt(A) = \sum_{\sigma\in S_N}^{N!} sgn(\sigma) B_\sigma(A)\]

Where the bipartite matching of \(\sigma\) is given by

\[B_\sigma(A) = \prod_{i=1}^N A_{i,\sigma(i)}\]

Parameters:: A – Square matrix (numpy array or array-like)
Returns:: Determinant value (complex or real)

ff.ff_combinatorics.dt_eigen(A)[source]

Computes the determinant of a matrix using eigenvalue decomposition.

\[dt(A) = \prod_{i=0}^{n-1} \lambda_i\]

with \(\lambda_i\) the eigenvalues of \(A\).

Parameters:: A – A NumPy array representing the square matrix.
Returns:: The determinant of the matrix. Returns None if the input is not a square matrix or if eigenvalue decomposition fails.

Graph Theory Module (ff_graph_theory)

Graph Theory and Visualization Functions for Free Fermion Library

This module contains graph theory algorithms and visualization functions, particularly focused on planar graphs and the Pfaffian ordering algorithm (FKT algorithm).

Key algorithms: - Pfaffian ordering algorithm (FKT algorithm) for planar graphs - Graph visualization and plotting functions - Perfect matching algorithms - Dual graph construction for planar embeddings

ff.ff_graph_theory.plot_graph_with_edge_weights(A, matching=None, title=None)[source]

Plot the graph associated with adjacency matrix A with edge weights.

Given a real [n x n] matrix A, interpret the entries as directed edge weights, then plot the graph with the edge weights.

Parameters:

A – A NumPy array representing the square matrix or NetworkX graph
matching – List of edges to highlight in red (optional)
title – Title for the plot (optional)

Returns:

NetworkX graph object with edge weights

ff.ff_graph_theory.generate_random_planar_graph(n, seed=None)[source]

Generate a random planar graph with n nodes.

Parameters:

n – The number of nodes in the graph
seed – Random seed for reproducibility (optional)

Returns:

A NetworkX graph object representing the planar graph. Returns None if a planar graph with n nodes cannot be generated within a reasonable number of attempts.

ff.ff_graph_theory.plot_planar_embedding(G, pfo=None, title=None)[source]

Plot the planar embedding of a graph with oriented faces.

Parameters:

G – The original NetworkX graph
pfo – Pfaffian ordering matrix (optional)
title – Title for the plot (optional)

ff.ff_graph_theory.dual_graph_H(G, F, T)[source]

Create the dual graph according to the PFO algorithm (Vazirani 1987).

“Construct a new graph H having one vertex corresponding to each face (including the infinite face) of G. Two vertices u and v of H are connected by an edge iff the corresponding faces share an edge not in T. Let r be the vertex in H corresponding to the infinite face of G. Root H at r.”

Parameters:

G – Original planar graph
F – List of faces from planar embedding
T – Spanning tree of G

Returns:

NetworkX graph H representing the dual graph

ff.ff_graph_theory.faces(graph, emb=None)[source]

Return the faces of an embedded graph using the networkx embedding scheme.

ADAPTED HEAVILY FROM SAGE: https://github.com/sagemath/sage/blob/develop/src/sage/graphs/generic_graph.py#L6815

Parameters:

graph – NetworkX graph object
emb – Planar embedding from NetworkX

Returns:

List of faces, where each face is a list of edges

ff.ff_graph_theory.complete_face(pfo, edges, verbose=False)[source]

Update one edge while maintaining the PFO condition.

Parameters:

pfo – Pfaffian ordering matrix
edges – List of edges in the face
verbose – Print debug information (optional)

Returns:

Updated PFO matrix

ff.ff_graph_theory.count_perfect_matchings(graph)[source]

Find all perfect matchings in a graph.

If planar, this method uses a pfo, otherwise it is computed via brute force

Parameters:: graph – A NetworkX graph
Returns:: The number of weighted matching for the given graph.

ff.ff_graph_theory.count_perfect_matchings_planar(graph)[source]

Find all perfect matchings in a planar graph using the PFO algorithm.

Parameters:: graph – A NetworkX graph
Returns:: The number of weighted matching for the given graph

ff.ff_graph_theory.find_perfect_matchings_brute(graph)[source]

Find all perfect matchings in a graph using a brute force enumeration.

Parameters:: graph – A NetworkX graph
Returns:: A list of tuples, where each tuple represents a perfect matching

ff.ff_graph_theory.pfo_algorithm(graph, verbose=False)[source]

Pfaffian ordering algorithm for planar graphs (FKT algorithm).

This algorithm is due to Kasteleyn 1967 by way of Vazirani 1987.

The algorithm constructs a pfaffian ordering of the edges of a planar graph such that the pfaffian of the resulting skew-symmetric matrix equals the number of perfect matchings in the graph.

Parameters:

graph – NetworkX planar graph
verbose – Print debug information and plots (optional)

Returns:

NumPy array representing the pfaffian ordering matrix

ff.ff_graph_theory.compute_tree_depth(graph)[source]

Compute the depth of a tree graph.

Parameters:: graph – NetworkX tree graph
Returns:: Integer representing the maximum depth of the tree

Distance Measures Module (ff_distance_measures)

Free Fermion Distance Measures Module

This module contains core distance measure functions for analyzing quantum states in the context of free fermion systems and stabilizer formalism. These measures quantify various notions of distance from classical or free-fermion behavior.

Key functionality:

Stabilizer distribution and Stabilizer Rényi Entropy (SRE)
General Rényi entropy and linear entropy measures
Fermionic covariance distribution analysis
Fermionic Anti-Flatness (FAF) distance measure
Kullback-Leibler divergence for probability distributions
Jensen-Shannon divergence for probability distributions
Trace distance for density matrices
Total variation distance for probability distributions
Bhattacharyya coefficient for probability distributions

The distance measures implemented here are fundamental tools for characterizing quantum states and their proximity to classical, stabilizer, or free-fermion subspaces in quantum many-body systems.

ff.ff_distance_measures.stabilizer_distribution(rho)[source]

Compute the probability distribution over Pauli operators in the stabilizer representation.

This function calculates the stabilizer distribution ξ_P for a quantum state ρ, where each element represents the probability of finding a specific Pauli operator P in the representation of the state. The distribution is defined as:

\[\xi_P = \frac{|\langle P \rangle_\rho|^2}{d Purity(\rho)}\]

where \(\langle P \rangle_\rho = \text{Tr}[\rho P]\) is the expectation value of Pauli operator P in state ρ, \(d = 2^n\) is the dimension of the Hilbert space for n qubits, and \(Purity(\rho) = \text{Tr}(\rho^2)\) is the purity of the state.

The stabilizer distribution is a fundamental quantity in stabilizer formalism and quantum error correction, providing insight into how “classical” or “stabilizer-like” a quantum state is.

Parameters:

rho (numpy.ndarray) – Input quantum state. Can be: - Density matrix of shape (2^n, 2^n) - Pure state wavefunction of shape (2^n,) or (2^n, 1) The function automatically detects and converts wavefunctions to density matrices.

Returns:

Stabilizer distribution ξ of length 4^n, where each: element ξ_j corresponds to the probability associated with the j-th Pauli operator in the Pauli group. The distribution is normalized such that Σ_j ξ_j = 1.

Return type:

numpy.ndarray

Raises:

AssertionError – If the computed distribution does not sum to 1 (within
numerical precision). –

Examples

>>> # For a single qubit in |0⟩ state
>>> psi = np.array([1, 0])
>>> xi = stabilizer_distribution(psi)
>>> len(xi)
4
>>> np.allclose(np.sum(xi), 1.0)
True

>>> # For a two-qubit maximally mixed state
>>> rho = np.eye(4) / 4
>>> xi = stabilizer_distribution(rho)
>>> len(xi)
16

Notes

For n qubits, the Pauli group has 4^n elements
The function uses the existing generate_pauli_group() for consistency
Numerical precision is handled with tolerance of 1e-9 for imaginary
parts
The distribution provides a complete characterization of the
state in the Pauli basis

References

Stabilizer Renyi Entropy paper by Leone, Oliviero, and Hamma Phys. Rev. Lett. 128, 050402 (2022). arXiv:2106.12587

ff.ff_distance_measures.SRE(rho, a=2)[source]

Compute the Stabilizer Rényi Entropy (SRE) of a quantum state.

The Stabilizer Rényi Entropy is a measure of how far a quantum state is from being a stabilizer state. It is defined as:

\[\text{SRE}_\alpha(\rho) = S_\alpha(\xi) - \log(d)\]

where \(S_\alpha(\xi)\) is the α-Rényi entropy of the stabilizer distribution ξ, and \(d = 2^n\) is the dimension of the Hilbert space.

For stabilizer states, SRE = 0, while for maximally mixed states or Haar-random states, SRE approaches its maximum value. This makes SRE a useful measure for quantifying “magic” or non-stabilizer resources in quantum computation.

Parameters:

rho (numpy.ndarray) – Input quantum state as density matrix of shape (2^n, 2^n) or wavefunction that will be converted to density matrix.
a (float, optional) – Rényi parameter α. Default is 2 (quadratic Rényi entropy). For α=1, this reduces to the von Neumann entropy case.

Returns:

The Stabilizer Rényi Entropy SRE_α(ρ).

SRE = 0 for stabilizer states
SRE > 0 indicates deviation from stabilizer behavior
Higher values indicate more “magic” or quantum resources

Return type:

float

Examples

>>> # Stabilizer state (should give SRE ≈ 0)
>>> psi_stab = np.array([1, 1, 1, 1]) / 2  # |++⟩ state
>>> rho_stab = np.outer(psi_stab, psi_stab.conj())
>>> sre_val = SRE(rho_stab)
>>> np.allclose(sre_val, 0, atol=1e-10)
True

>>> # Maximally mixed state (should give positive SRE)
>>> rho_mixed = np.eye(4) / 4
>>> sre_mixed = SRE(rho_mixed)
>>> sre_mixed > 0
True

Notes

The function uses the stabilizer_distribution() to compute ξ
The logarithmic offset -log(d) ensures SRE=0 for stabilizer states
For α=2, this is sometimes called the “linear entropy” variant
The measure is basis-independent and unitarily invariant

References

Leone, L., Oliviero, S. F., & Hamma, A. (2021). Stabilizer Rényi entropy. Physical Review Letters, 128(5), 050402. arXiv:2106.12587

ff.ff_distance_measures.renyi_entropy(p, a)[source]

Compute the Rényi entropy of a probability distribution.

The Rényi entropy is a generalization of the Shannon entropy and is defined as:

\[\]

S_alpha(p) = frac{1}{1-alpha} logleft(sum_i p_i^alpharight)

For α=1, this reduces to the Shannon entropy:

\[S_1(p) = -\sum_i p_i \log(p_i)\]

For α=2, this gives the collision entropy:

\[S_2(p) = -\log\left(\sum_i p_i^2\right)\]

The Rényi entropy provides a family of entropy measures that capture different aspects of the probability distribution’s structure and are widely used in information theory, quantum information, and statistical physics.

Parameters:

p (numpy.ndarray) – Probability distribution as a 1D array. Must be normalized such that Σ_i p_i = 1 and all elements p_i ≥ 0.
a (float) – Rényi parameter α. Must be non-negative and α ≠ 1. - α → 0: Max entropy (log of support size) - α = 1: Shannon entropy (handled as special case) - α = 2: Collision entropy - α → ∞: Min entropy (-log(max(p_i)))

Returns:

The α-Rényi entropy S_α(p).

For α < 1: S_α ≥ S_1 (Rényi entropy is larger)
For α > 1: S_α ≤ S_1 (Rényi entropy is smaller)
All Rényi entropies are non-negative for normalized distributions

Return type:

float

Examples

>>> # Uniform distribution
>>> p_uniform = np.ones(4) / 4
>>> s1 = renyi_entropy(p_uniform, 1)  # Shannon entropy
>>> s2 = renyi_entropy(p_uniform, 2)  # Collision entropy
>>> np.allclose(s1, np.log(4))  # Should equal log(4) for uniform
True
>>> s2 < s1  # Collision entropy is smaller
True

>>> # Delta distribution (pure state)
>>> p_delta = np.array([1, 0, 0, 0])
>>> renyi_entropy(p_delta, 2)
0.0

Notes

For α=1, the function uses scipy.stats.entropy for numerical
stability
The function handles the limit α→1 by using the Shannon entropy
formula
Numerical precision may affect results for very small probabilities
The function assumes input is a valid probability distribution

References

Rényi, A. (1961). On measures of entropy and information.
Cover, T., & Thomas, J. (2006). Elements of information theory.

ff.ff_distance_measures.linear_entropy(p)[source]

Compute the linear entropy of a probability distribution.

The linear entropy is a measure of mixedness defined as:

\[S_{\text{lin}}(p) = 1 - d \|p\|_2^2\]

where \(d\) is the dimension (length of the probability vector) and \(\|p\|_2^2\) is the squared 2-norm of the distribution. This is closely related to the 2-Rényi entropy and provides a simple measure of how far a distribution is from being uniform.

For quantum states, the linear entropy measures the degree of mixedness:

\[S_{\text{lin}}(\rho) = 1 - \text{Tr}(\rho^2)\]

Parameters:

p (numpy.ndarray) – Probability distribution as a 1D array. Must be normalized such that Σ_i p_i = 1 and all elements p_i ≥ 0. The length must be a perfect square (d = 2^n for some integer n).

Returns:

The linear entropy S_lin(p).

S_lin = 0 for pure states (delta distributions)
S_lin = (d-1)/d for maximally mixed states
Higher values indicate more mixedness/entropy

Return type:

float

Raises:

AssertionError – If the length of p is not a perfect square.

Examples

>>> # Pure state (delta distribution)
>>> p_pure = np.array([1, 0, 0, 0])
>>> linear_entropy(p_pure)
0.0

>>> # Maximally mixed state
>>> p_mixed = np.ones(4) / 4
>>> lin_ent = linear_entropy(p_mixed)
>>> np.allclose(lin_ent, 3/4)  # (d-1)/d = 3/4 for d=4
True

Notes

The dimension d must be a perfect square for consistency with qubit
systems
This measure is computationally efficient compared to von Neumann
entropy
Linear entropy is closely related to purity: Purity = 1 - S_lin
For quantum states: 0 ≤ S_lin ≤ (d-1)/d

ff.ff_distance_measures.cov_distribution(rho)[source]

Compute the fermionic covariance distribution for a quantum state.

This function computes the probability distribution over fermionic covariance matrix eigenvalues, which characterizes the fermionic correlations in the state. The distribution is obtained from the eigenvalues of the fermionic covariance matrix:

\[M_{ij} = -i \langle \alpha_i \alpha_j \rangle_c\]

where \(\langle \alpha_i \alpha_j \rangle_c\) is the connected correlation function and \(\alpha_i\) are the Majorana operators under Jordan-Wigner transformation.

The covariance matrix M is real and antisymmetric, and its eigenvalues come in pairs ±λ_k. The distribution is formed from the squared eigenvalues.

Parameters:

rho (numpy.ndarray) – Input quantum state. Can be: - Density matrix of shape (2^n, 2^n) - Pure state wavefunction of shape (2^n,) or (2^n, 1) The function automatically detects and converts wavefunctions to density matrices.

Returns:

Distribution over squared eigenvalues of the covariance: matrix. The length depends on the number of non-zero eigenvalue pairs.

Return type:

numpy.ndarray

Examples

>>> # Two-site system
>>> rho = np.eye(4) / 4  # Maximally mixed state
>>> pm = cov_distribution(rho)
>>> len(pm) <= 4  # At most 2 sites worth of eigenvalues
True

Notes

Uses Jordan-Wigner Majorana operators from ff_lib
The covariance matrix is computed using compute_cov_matrix
from ff_lib
Eigenvalues are extracted using Schur decomposition for numerical
stability
Only eigenvalues with magnitude > 1e-10 are included in the
distribution

References

Fermionic Magic Resources of Quantum Many-Body Systems. Sierant, Stornati, and Turkeshi (arXiv:2506.00116)

ff.ff_distance_measures.total_variation_distance(p, q)[source]

Compute the total variation distance between two probability vectors.

The trace distance is a measure of distinguishability between two probability distributions. It is defined as: .. math:

D_{\text{trace}}(p, q) = \frac{1}{2} \sum_i |p_i - q_i|

where \(p\) and \(q\) are the two probability distributions. :param p: First probability distribution as a 1D array.

Must be normalized such that Σ_i p_i = 1 and all elements p_i ≥ 0.

Parameters:

q (numpy.ndarray) – Second probability distribution as a 1D array. Must be normalized such that Σ_i q_i = 1 and all elements q_i ≥ 0.

Returns:

The trace distance D_trace(p, q).

0 ≤ D_trace ≤ 1
D_trace = 0 if and only if p = q

Return type:

float

Examples

>>> # Simple distributions
>>> p = np.array([0.5, 0.5])
>>> q = np.array([0.9, 0.1])
>>> dtrace = trace_distance(p, q)
>>> dtrace > 0
True
>>> dtrace_zero = trace_distance(p, p)
>>> np.allclose(dtrace_zero, 0)
True

Notes

The function coverts input wave functions and density matrices to
probability distributions

References

Nielsen, M. A., & Chuang, I. L. (2010). Quantum computation and quantum information.

ff.ff_distance_measures.trace_distance(A, B)[source]

Computes the trace distance between two matrices A and B.

The trace distance is a quantum measure used to quantify the distinguishability between two quantum states represented by density matrices A and B. It is defined as: .. math:: D_{trace}(A, B) = frac{1}{2} text{Tr} |A - B|

where \(|X| = \sqrt{X^\dagger X}\) is the positive square root of the operator X. :param A: First density matrix (numpy array or array-like) :param B: Second density matrix (numpy array or array-like)

Returns:: Trace distance value (float)

Notes

The function assumes A and B are matrices and casts them to density
matrices as needed
The trace distance ranges from 0 (identical states) to 1 (orthogonal
states)
This measure is widely used in quantum information theory to assess
state distinguishability

References

Nielsen, M. A., & Chuang, I. L. (2010). Quantum computation and quantum information.

ff.ff_distance_measures.relative_entropy(p, q)[source]

Compute the Kullback-Leibler divergence between two probability vectors.

The KL divergence is a measure of how one probability distribution diverges from a second, expected probability distribution. It is defined as:

\[D_{KL}(p || q) = \sum_i p_i \log\left(\frac{p_i}{q_i}\right)\]

where \(p\) is the true distribution and \(q\) is the reference distribution.

Parameters:

p (numpy.ndarray) – True probability distribution as a 1D array. Must be normalized such that Σ_i p_i = 1 and all elements p_i ≥ 0.
q (numpy.ndarray) – Reference probability distribution as a 1D array. Must be normalized such that Σ_i q_i = 1 and all elements q_i ≥ 0.

Returns:

The KL divergence D_KL(p || q).

D_KL ≥ 0, with equality if and only if p = q

Return type:

float

Examples

>>> # Simple distributions
>>> p = np.array([0.5, 0.5])
>>> q = np.array([0.9, 0.1])
>>> dkl = relative_entropy(p, q)
>>> dkl > 0
True
>>> dkl_zero = relative_entropy(p, p)
>>> np.allclose(dkl_zero, 0)
True

Notes

The function uses scipy.stats.entropy for numerical stability
KL divergence is not symmetric: D_KL(p || q) ≠ D_KL(q || p)
The function coverts input wave functions and density matrices to
probability distributions

References

Kullback, S., & Leibler, R. A. (1951).
On information and sufficiency.
Cover, T. M., & Thomas, J. A. (2006).
Elements of information theory.

ff.ff_distance_measures.bhattacharyya_coeff(p, q)[source]

Compute Bhattacharyya coefficient between two probability distributions.

The Bhattacharyya coefficient is a measure of similarity between two probability distributions. It is defined as:

\[BC(p, q) = \sum_i \sqrt{p_i q_i}\]

where \(p\) and \(q\) are the two probability distributions.

Parameters:

p (numpy.ndarray) – First probability distribution as a 1D array. Must be normalized such that Σ_i p_i = 1 and all elements p_i ≥ 0.
q (numpy.ndarray) – Second probability distribution as a 1D array. Must be normalized such that Σ_i q_i = 1 and all elements q_i ≥ 0.

Returns:

The Bhattacharyya coefficient BC(p, q).

Return type:

float

Notes

The function coverts input wave functions and density matrices to
probability distributions

ff.ff_distance_measures.jensen_shannon_divergence(p, q)[source]

Compute the Jensen-Shannon (JS) divergence between two probability distributions.

The JS divergence is a symmetric and bounded measure of similarity between two probability distributions. It is defined as:

\[\]

D_{JS}(p || q) = frac{1}{2} D_{KL}(p || m) + frac{1}{2} D_{KL}(q || m)

where \(m = \frac{1}{2}(p + q)\) is the average distribution, and \(D_{KL}\) is the Kullback-Leibler divergence.

Parameters:

p (numpy.ndarray) – First probability distribution as a 1D array. Must be normalized such that Σ_i p_i = 1 and all elements p_i ≥ 0.
q (numpy.ndarray) – Second probability distribution as a 1D array. Must be normalized such that Σ_i q_i = 1 and all elements q_i ≥ 0.

Returns:

The JS divergence D_JS(p || q).

0 ≤ D_JS ≤ log(2)
D_JS = 0 if and only if p = q

Return type:

float

Examples

>>> # Simple distributions
>>> p = np.array([0.5, 0.5])
>>> q = np.array([0.9, 0.1])
>>> djs = jensen_shannon_divergence(p, q)
>>> djs > 0
True
>>> djs_zero = jensen_shannon_divergence(p, p)
>>> np.allclose(djs_zero, 0)
True

Notes

The function uses scipy.stats.entropy for numerical stability
JS divergence is symmetric: D_JS(p || q) = D_JS(q || p)
The function coverts input wave functions and density matrices to
probability distributions

References

Lin, J. (1991). Divergence measures based on the Shannon entropy.
Cover, T. M., & Thomas, J. A. (2006). Elements of information theory.

ff.ff_distance_measures.FAF(rho, k=2)[source]

Compute the Fermionic Anti-Flatness (FAF) distance measure.

The Fermionic Anti-Flatness is a distance measure that quantifies how far a fermionic state is from being a free-fermion (Gaussian) state. It is defined as:

\[\text{FAF}_k(\rho) = n - \|\mathbf{p}\|_k^k\]

where \(n\) is the number of fermionic sites, \(\mathbf{p}\) is the distribution of squared eigenvalues of the fermionic covariance matrix, and \(\|\cdot\|_k\) is the k-norm.

Equivalently, this can be expressed as:

\[\text{FAF}_k(\rho) = n - \frac{1}{2} \text{Tr}[(M^T M)^k]\]

where \(M\) is the fermionic covariance matrix.

For free-fermion states, FAF = 0, while for highly correlated fermionic states, FAF approaches its maximum value. This makes FAF a useful measure for quantifying fermionic correlations and deviations from free-fermion behavior.

Parameters:

rho (numpy.ndarray) – Input quantum state as density matrix of shape (2^n, 2^n) or wavefunction that will be converted to density matrix.
k (int, optional) – Norm parameter for the FAF measure. Default is 2. Higher values of k make the measure more sensitive to large eigenvalues in the covariance spectrum.

Returns:

The Fermionic Anti-Flatness FAF_k(ρ).

FAF = 0 for free-fermion (Gaussian) states
FAF > 0 indicates deviation from free-fermion behavior
Higher values indicate stronger fermionic correlations

Return type:

float

Examples

>>> # Free-fermion state (should give FAF ≈ 0)
>>> # This would require constructing a proper Gaussian state
>>> rho_ff = np.eye(4) / 4  # Placeholder - not actually free-fermion
>>> faf_val = FAF(rho_ff)
>>> isinstance(faf_val, float)
True

Notes

Uses the fermionic covariance matrix computed via Jordan-Wigner
transformation
The function verifies consistency between different computational
approaches
For k=2, this reduces to a quadratic measure of correlations
The measure is basis-independent within the fermionic representation

References

Sierant, P., Stornati, G., & Turkeshi, X. (2024). Fermionic Magic Resources of Quantum Many-Body Systems. arXiv:2506.00116

The distance measures module provides fundamental tools for quantifying quantum states and their proximity to classical, stabilizer, or free-fermion subspaces. Key measures include:

Stabilizer Rényi Entropy (SRE): Quantifies deviation from stabilizer states
Fermionic Anti-Flatness (FAF): Measures distance from free-fermion behavior
General entropy measures: Rényi entropy, linear entropy for probability distributions
Stabilizer and covariance distributions: Core probability distributions for analysis
Distance measures: Total variation distance, trace distance, relative entropy
Divergence measures: Jensen-Shannon divergence, Bhattacharyya coefficient

Distance Measures Functions

`ff.ff_distance_measures.stabilizer_distribution`(rho)	Compute the probability distribution over Pauli operators in the stabilizer representation.
`ff.ff_distance_measures.SRE`(rho[, a])	Compute the Stabilizer Rényi Entropy (SRE) of a quantum state.
`ff.ff_distance_measures.FAF`(rho[, k])	Compute the Fermionic Anti-Flatness (FAF) distance measure.
`ff.ff_distance_measures.renyi_entropy`(p, a)	Compute the Rényi entropy of a probability distribution.
`ff.ff_distance_measures.linear_entropy`(p)	Compute the linear entropy of a probability distribution.
`ff.ff_distance_measures.cov_distribution`(rho)	Compute the fermionic covariance distribution for a quantum state.
`ff.ff_distance_measures.total_variation_distance`(p, q)	Compute the total variation distance between two probability vectors.
`ff.ff_distance_measures.trace_distance`(A, B)	Computes the trace distance between two matrices A and B.
`ff.ff_distance_measures.relative_entropy`(p, q)	Compute the Kullback-Leibler divergence between two probability vectors.
`ff.ff_distance_measures.jensen_shannon_divergence`(p, q)	Compute the Jensen-Shannon (JS) divergence between two probability distributions.
`ff.ff_distance_measures.bhattacharyya_coeff`(p, q)	Compute Bhattacharyya coefficient between two probability distributions.

Random States Module (ff_random_states)

Free Fermion Random States Module

This module contains random state generators and path construction functions for free fermion systems, including Haar random states, Clifford states, and various free fermion state generation methods.

Key functionality:

Random qubit pure states (Haar random)
Random Clifford states using Stim package
Random free fermion mixed states from random Hamiltonians
Random free fermion mixed states from rotated probability distributions
Random free fermion pure states with various constraints
Path construction functions for state interpolation

ff.ff_random_states.random_qubit_state(n, seed=None)[source]

Generate a Haar random qubit state.

This function generates a uniformly distributed mixed state over the n-qubit Hilbert space using the Haar measure on the unitary group. The state is created by applying a Haar random unitary to a random pdf sampled according to the Dirichlet distribution.

Parameters:

n (int) – The number of qubits
seed (int, optional) – Random seed for reproducibility

Returns:

A random qubit state of dimension (2^n, 2^n) as a: normalized density matrix

Return type:

numpy.ndarray

Raises:

ValueError – If n is not a positive integer
TypeError – If n is not an integer

Examples

>>> # Generate a random single qubit state
>>> rho = random_qubit_state(1, seed=42)
>>> print(rho.shape)
(2, 2)
>>> print(np.allclose(np.trace(rho), 1.0))
True

>>> # Generate a random two-qubit state
>>> rho = random_qubit_state(2, seed=123)
>>> print(rho.shape)
(4, 4)

ff.ff_random_states.random_qubit_pure_state(n, seed=None)[source]

Generate a Haar random qubit state.

This function generates a uniformly distributed pure state over the n-qubit Hilbert space using the Haar measure on the unitary group. The state is created by applying a Haar random unitary to the computational basis state |00…0⟩.

Parameters:

n (int) – The number of qubits
seed (int, optional) – Random seed for reproducibility

Returns:

A random qubit state of dimension (2^n, 1) as a: normalized complex vector

Return type:

numpy.ndarray

Raises:

ValueError – If n is not a positive integer
TypeError – If n is not an integer

Examples

>>> # Generate a random single qubit state
>>> psi = random_qubit_pure_state(1, seed=42)
>>> print(psi.shape)
(2, 1)
>>> print(np.allclose(np.linalg.norm(psi), 1.0))
True

>>> # Generate a random two-qubit state
>>> psi = random_qubit_pure_state(2, seed=123)
>>> print(psi.shape)
(4, 1)

Notes

The generated states are uniformly distributed according to the Haar measure on the space of pure quantum states
Memory usage scales as O(4^n), so use with caution for large n
The state is always normalized to unit length

ff.ff_random_states.random_CHP_state(n_qubits)[source]

Generate a random Clifford state using Stim package.

This function generates a random stabilizer state (Clifford state) using the Stim quantum circuit simulator package. The state is uniformly distributed over the space of n-qubit stabilizer states.

Parameters:

n_qubits (int) – The number of qubits for the CHP state

Returns:

A random Clifford state of dimension (2^n_qubits, 1): as a normalized complex vector

Return type:

numpy.ndarray

Raises:

ImportError – If the stim package is not available
ValueError – If n_qubits is not a positive integer
TypeError – If n_qubits is not an integer

Examples

>>> # Generate a random single qubit Clifford state
>>> psi = random_CHP_state(1)
>>> print(psi.shape)
(2, 1)

>>> # Generate a random two-qubit Clifford state
>>> psi = random_CHP_state(2)
>>> print(psi.shape)
(4, 1)

Notes

Requires the ‘stim’ package to be installed
Clifford states form a subset of all quantum states that can be efficiently simulated classically
The generated states are uniformly distributed over the stabilizer state space

ff.ff_random_states.random_FF_state_randH(n_sites, seed=None)[source]

Generate a random free fermion mixed state from a random Hamiltonian.

This function generates a random free fermion thermal state by constructing a random quadratic Hamiltonian and computing the corresponding FF state via exponentiation. The resulting state is a density matrix.

NOTE: This function does not generate Haar random states since: exponentiating a random Hamiltonian does not generate correctly distribute eigenvalues of the resulting random matrix.

Parameters:

n_sites (int) – The number of fermionic sites
seed (int, optional) – Random seed for reproducibility

Returns:

A random free fermion state of dimension: (2^n_sites, 2^n_sites) as a normalized density matrix

Return type:

numpy.ndarray

Raises:

ValueError – If n_sites is not a positive integer
TypeError – If n_sites is not an integer

Examples

>>> # Generate a random 2-site FF mixed state
>>> rho = random_FF_state_randH(2, seed=42)
>>> print(rho.shape)
(4, 4)
>>> print(np.allclose(np.trace(rho), 1.0))
True

Notes

The Hamiltonian is constructed from random Gaussian matrices
The resulting state is a thermal state at infinite temperature
Memory usage scales as O(4^n_sites)

ff.ff_random_states.random_FF_state_rotPDF(n_sites, returnS=False, seed=None)[source]

Generate random free fermion state from rotated probability vector.

This function generates a random free fermion state by applying a random free fermion rotation to a random probability distribution (Dirichlet distributed eigenvalues).

Parameters:

n_sites (int) – The number of fermionic sites
returnS (bool, optional) – If True, also return the probability distribution (default: False)
seed (int, optional) – Random seed for reproducibility

Returns:

If returnS is False, returns a random free: fermion state of dimension (2^n_sites, 2^n_sites). If returnS is True, returns (rho, s) where s is the probability distribution used.

Return type:

numpy.ndarray or tuple

Raises:

ValueError – If n_sites is not a positive integer
TypeError – If n_sites is not an integer

Examples

>>> # Generate a random 2-site FF mixed state
>>> rho = random_FF_state_rotPDF(2, seed=42)
>>> print(rho.shape)
(4, 4)

>>> # Generate state and return probability distribution
>>> rho, s = random_FF_state_rotPDF(2, returnS=True, seed=42)
>>> print(len(s))
4

Notes

Uses Dirichlet distribution for random eigenvalues
The rotation preserves the free fermion structure
Memory usage scales as O(4^n_sites)

ff.ff_random_states.random_FF_pure_state_W0(n_sites, seed=None)[source]

Generate a Haar random free fermion pure state from vacuum rotation.

This function generates free-fermion pure states that are uniformly distributed over the space of free-fermion states using Haar random symplectic transformations applied to the vacuum state |0⟩.

Parameters:

n_sites (int) – The number of fermionic sites
seed (int, optional) – Random seed for reproducibility

Returns:

A normalized Haar random free fermion pure: state of dimension (2^n_sites, 1). If returnH is True, also returns the generator matrix H.

Return type:

numpy.ndarray or tuple

Raises:

ValueError – If n_sites is not a positive integer
TypeError – If n_sites is not an integer

Examples

>>> # Generate a random 2-site FF pure state
>>> psi = random_FF_pure_state_W0(2, seed=42)
>>> print(psi.shape)
(4, 1)
>>> print(np.allclose(np.linalg.norm(psi), 1.0))
True

Notes

Uses Haar random symplectic rotations
The state is created by rotating the vacuum state
Memory usage scales as O(4^n_sites)

ff.ff_random_states.random_FF_pure_state_WN(n_sites, N=None, seed=None)[source]

Generates Haar random free fermion pure state with fixed particle number

This function generates free-fermion pure states with a specified particle number N, uniformly distributed over the space of N-particle free-fermion states using Haar random symplectic transformations.

NOTE: This function is most-likely generating states that will not have: fixed particle number. Further testing and validation is required to ensure that the generated states are in a fixed-N subspace (also which fermionic basis should be used).

Parameters:

n_sites (int) – The number of fermionic sites
N (int, optional) – The particle number. If None, chosen randomly (default: None)
seed (int, optional) – Random seed for reproducibility

Returns:

A normalized Haar random free fermion pure: state of dimension (2^n_sites, 1) with N particles. If returnH is True, also returns the generator matrix H.

Return type:

numpy.ndarray or tuple

Raises:

ValueError – If n_sites is not a positive integer or N is invalid
TypeError – If n_sites is not an integer

Examples

>>> # Generate a random 2-site FF pure state with 1 particle
>>> psi = random_FF_pure_state_WN(2, N=1, seed=42)
>>> print(psi.shape)
(4, 1)

>>> # Generate with random particle number
>>> psi = random_FF_pure_state_WN(3, seed=42)
>>> print(psi.shape)
(8, 1)

Notes

The particle number N must be between 0 and n_sites
Uses Haar random symplectic rotations
Memory usage scales as O(4^n_sites)

ff.ff_random_states.random_FF_pure_state_CN(n_sites, seed=None)[source]

Generate random free fermion pure state using random orbital rotations.

This function generates a random free fermion pure state by first creating a random vacuum state through projections, then applying creation operators according to a random occupation pattern, all with randomly rotated orbitals.

NOTE: This function is most-likely generating states that will not have: fixed particle number. Further testing and validation is required to ensure that the generated states are in a fixed-N subspace (also which fermionic basis should be used).

Parameters:

n_sites (int) – The number of fermionic sites
seed (int, optional) – Random seed for reproducibility

Returns:

A random free fermion pure state of: dimension (2^n_sites, 1)

Return type:

numpy.ndarray

Raises:

ValueError – If n_sites is not a positive integer
TypeError – If n_sites is not an integer

Examples

>>> # Generate a random 2-site FF pure state
>>> psi = random_FF_pure_state_CN(2, seed=42)
>>> print(psi.shape)
(4, 1)

>>> # Generate a random 3-site FF pure state
>>> psi = random_FF_pure_state_CN(3, seed=123)
>>> print(psi.shape)
(8, 1)

Notes

Uses random orbital rotations respecting fermionic antisymmetry
Creates a random vacuum through projections
Applies creation operators according to random occupation
Memory usage scales as O(4^n_sites)

ff.ff_random_states.get_orthogonal_vectors(v)[source]

Get orthogonal vectors to a given vector using QR decomposition.

Given a vector v (n x 1), return n-1 vectors that are orthogonal to v and form an orthonormal basis together with the normalized v.

Parameters:

v (numpy.ndarray) – The input vector (n x 1 or n-dimensional)

Returns:

A matrix of shape (n, n) where the first column is: the normalized input vector and the remaining n-1 columns are orthogonal to it

Return type:

numpy.ndarray

Raises:

ValueError – If the input vector has insufficient dimensions
TypeError – If the input is not a numpy array

Examples

>>> # Get orthogonal vectors to a 3D vector
>>> v = np.array([[1], [2], [3]])
>>> orth_vecs = get_orthogonal_vectors(v)
>>> print(orth_vecs.shape)
(3, 3)

>>> # Check orthogonality
>>> v_norm = v / np.linalg.norm(v)
>>> print(np.allclose(orth_vecs[:, 0:1], v_norm))
True

Notes

Uses QR decomposition for numerical stability
Returns an orthonormal basis including the normalized input vector
For zero vectors, returns standard basis vectors

ff.ff_random_states.build_unitary_path(w, v)[source]

Build a unitary path between two vectors.

This function constructs a smooth unitary path between two normalized vectors w and v, parameterized by t ∈ [0,1], such that the path starts at w (t=0) and ends at v (t=1).

Parameters:

w (numpy.ndarray) – The initial vector
v (numpy.ndarray) – The target vector

Returns:

A function path(t) that returns the interpolated vector: at parameter t ∈ [0,1]

Return type:

callable

Raises:

ValueError – If vectors have incompatible dimensions
TypeError – If inputs are not numpy arrays

Examples

>>> # Create a path between two 3D vectors
>>> w = np.array([[1], [0], [0]])
>>> v = np.array([[0], [1], [0]])
>>> path = build_unitary_path(w, v)
>>>
>>> # Evaluate at t=0 (should be close to w)
>>> start = path(0.0)
>>> print(np.allclose(start, w / np.linalg.norm(w)))
True
>>>
>>> # Evaluate at t=1 (should be close to v)
>>> end = path(1.0)
>>> print(np.allclose(end, v / np.linalg.norm(v)))
True

Notes

The path preserves normalization at all points
Uses matrix exponential for smooth interpolation
The path is the shortest geodesic on the unit sphere

ff.ff_random_states.build_linear_path(w, v)[source]

Build a linear path between two vectors.

This function constructs a linear interpolation path between two vectors w and v, parameterized by s ∈ [0,1], such that the path starts at w (s=0) and ends at v (s=1). The path maintains normalization at each point.

Parameters:

w (numpy.ndarray) – The initial vector
v (numpy.ndarray) – The target vector

Returns:

A function path(s) that returns the normalized interpolated: vector at parameter s ∈ [0,1]

Return type:

callable

Raises:

ValueError – If vectors have incompatible dimensions
TypeError – If inputs are not numpy arrays

Examples

>>> # Create a linear path between two 3D vectors
>>> w = np.array([[1], [0], [0]])
>>> v = np.array([[0], [1], [0]])
>>> path = build_linear_path(w, v)
>>>
>>> # Evaluate at s=0 (should be close to w)
>>> start = path(0.0)
>>> print(np.allclose(start, w / np.linalg.norm(w)))
True
>>>
>>> # Evaluate at s=1 (should be close to v)
>>> end = path(1.0)
>>> print(np.allclose(end, v / np.linalg.norm(v)))
True

Notes

The path is a straight line in the ambient space, then normalized
Simpler than unitary path but may not be the shortest geodesic
Maintains normalization at all points

The random states module provides comprehensive random state generation capabilities for quantum systems, including:

Haar random states: Uniformly distributed pure states over qubit Hilbert spaces
Clifford states: Random stabilizer states using the Stim package (optional dependency)
Free fermion states: Various methods for generating random FF states with different constraints
Path construction: Unitary and linear interpolation between quantum states

Note: Some functions require the optional stim package for Clifford state generation.

Utilities Module (ff_utils)

Free Fermion Utilities Module

Common utility functions for the free fermion codebase.

ff.ff_utils.print_custom(obj, k=9)[source]

Custom print function with small number suppression

Parameters:

obj – Any object to be printed
k – The number of decimal places to print

Returns:

None

ff.ff_utils.clean(obj, threshold=1e-06)[source]

Clean small numerical values from arrays or matrices.

Parameters:

obj – array, scalar, NumPy array or matrix
threshold – Values below this threshold are set to zero

Note: if threshold is an integer, it will be converted to 10^-threshold

Returns:: Cleaned obj with rounded values and small values set to zero.

ff.ff_utils.formatted_output(obj, precision=6)[source]

Format numerical output with specified precision.

Parameters:

obj – Object to format
precision – Number of decimal places

Returns:

Formatted string representation

ff.ff_utils.generate_random_bitstring(n, k)[source]

Generates a random bit string of length n with Hamming weight k.

Based on np.random.choice

Parameters:

n – The length of the bit string.
k – The Hamming weight (number of 1s).

Returns:

A NumPy array representing the bit string, or None if k is invalid.

ff.ff_utils.kron_plus(a, b)[source]

Computes the direct sum of two matrices

Parameters:

a – First matrix
b – Second matrix

Returns:

Direct sum matrix [[a, 0], [0, b]]

ff.ff_utils.cast_to_pdf(rho)[source]: Casts input to a probability distribution.

ff.ff_utils.cast_to_density_matrix(rho)[source]

Casts input to a density matrix.

Parameters:: rho – Input: can be a density matrix, wavefunction, or prob. distro.
Returns:: Density matrix as a NumPy array.

Notes

Threshold for positivity checks is 1e-12

ff.ff_utils.analyze_pdf(rho, name=None, stem=True)[source]

Population analysis and visualization for probability distributions.

Analyzes and visualizes probability distributions from density matrices or wavefunctions. Computes both diagonal elements and eigenvalues of the density matrix and displays them with optional logarithmic scaling.

Parameters:

rho (numpy.ndarray) –
Density matrix, wavefunction, or prob. distro. - If 2D square matrix: treated as density matrix - If 1D array or column vector: interpreted based on normalization
- L2 norm ≈ 1: treated as normalized wavefunction
- L1 norm ≈ 1: treated as normalized probability distribution
name (str, optional) – Name for the plot title. If provided, will appear as “Population analysis of {name}”. Default is None.
stem (bool, optional) – Whether to use stem (True) or line plot (False). Default is True.

Returns:

Function creates a matplotlib plot but does not return values.

Return type:

None

Raises:

ValueError – If input array has invalid dimensions or normalization.
TypeError – If input is not a numpy array or array-like object.

Notes

The function automatically converts wavefunctions to density matrices
Diagonal elements are plotted in red, eigenvalues in black dashed
Y-axis uses logarithmic scaling to highlight small probabilities
Small numerical values cleaned using the FF library clean() function
Plot is created but not displayed; use plt.show() separately

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # Analyze a simple 2x2 density matrix
>>> rho = np.array([[0.7, 0.1], [0.1, 0.3]])
>>> analyze_pdf(rho, name="2-level system")
>>> plt.show()

>>> # Analyze a wavefunction
>>> psi = np.array([1/np.sqrt(2), 1/np.sqrt(2)])
>>> analyze_pdf(psi, name="Bell state", stem=False)
>>> plt.show()

ff.ff_utils.partial_trace_diagblocksum(AB, d)[source]

This function implements the partial trace using a block matrix approach, where the composite matrix AB is viewed as a \(d \times d\) array of \(d_2 \times d_2\) blocks, and sum of the diagonal blocks gives the partial trace result.

Parameters:

AB – Any composite matrix over space \(H_1 \otimes H_2\).
d – The dimension of the reduced state. Must be a positive integer such

:param that D is divisible by \(d\): :param where D is the dimension of the: :param composite system AB.:

Returns:: the reduced matrix of subsystem B with shape :math:` dtimes d `.
Return type:: B

Notes

Dimension of subsystem A is automatically computed as \(d_1 = D/d\)
The composite system dimension D must be exactly divisible by \(d\)
This is the partial trace over subsystem A.

Mathematical Details

The reduced density matrix A is computed as:

\[A = \sum_j AB_{j,j}\]

where \(AB_{i,j}\) is the \((i,j)\)-th block of size \(d \times d\).

ff.ff_utils.partial_trace_blockTr(AB, d)[source]

This function implements the partial trace using a block matrix approach, where the composite matrix AB is viewed as a \(d \times d\) array of \(d_1 \times d_1\) blocks, and the trace of each block gives the corresponding matrix element of the result matrix.

Parameters:

AB – Any composite matrix over space \(H_1 \otimes H_2\).
d – The dimension of the reduced state. Must be a positive integer such

:param that D is divisible by \(d\): :param where D is the dimension of the: :param composite system AB.:

Returns:: the reduced matrix of subsystem A with shape :math:` dtimes d `.
Return type:: A

Notes

Dimension of subsystem B is automatically computed as \(d_2 = D/d\)
The composite system dimension D must be exactly divisible by \(d\)
This is the partial trace over subsystem B.

Mathematical Details

The reduced density matrix element \((i,j)\) is computed as:

\[B[i,j] = \text{Tr}(AB_{i,j})\]

where \(AB_{i,j}\) is the \((i,j)\)-th block of size \(d_1 \times d_1\).

Function Index

Jordan-Wigner Transformations

`ff.ff_lib.jordan_wigner_lowering`(n_sites)	Define annihilation operators using the Jordan-Wigner transform.
`ff.ff_lib.jordan_wigner_alphas`(n_sites)	Define raising and lowering operators using the Jordan-Wigner transform.
`ff.ff_lib.jordan_wigner_majoranas`(n_sites)	Generate a set of Majorana operators under Jordan-Wigner.
`ff.ff_lib.rotate_operators`(C, alphas[, ...])	Transform set of operators using a matrix C per:
`ff.ff_lib._perform_rotation`(C, alphas[, verbose])	Helper function to perform the actual rotation of operators.

Matrix Construction

`ff.ff_lib.build_V`(n_sites, A[, Z])	Build the generator V matrix of dimension 2*N for N sites.
`ff.ff_lib.build_H`(n_sites, A[, B])	Build the H block coefficient matrix of dimension 2*N for N sites.
`ff.ff_lib.build_Omega`(N)	Build the Omega matrix of dimension 2*N for N sites.
`ff.ff_lib.build_K`(n_sites)	Build the symplectic transform from standard to canonical eigenstruct
`ff.ff_lib.build_reordering_xx_to_xp`(n_sites)	A permutation matrix that transforms ordering [x_1, x_2.
`ff.ff_lib.permutation_to_matrix`(permutation)	Transform a permutation (list) into a permutation matrix using NumPy.
`ff.ff_lib.pauli_matrices`()	Define the Pauli matrices as numpy arrays.
`ff.ff_lib.build_op`(n_sites, R, alphas[, ...])	Build the N-body lift of a quadratic coefficient matrix

Gaussian States and Correlation Matrices

`ff.ff_lib.generate_gaussian_state`(n_sites, H)	Generate a Gaussian state using Hamiltonian H.
`ff.ff_lib.compute_cov_matrix`(rho[, n_sites, ...])	Calculate the fermionic covariance matrix of rho
`ff.ff_lib.compute_2corr_matrix`(rho, n_sites)	Calculate the two-point correlation matrix
`ff.ff_lib.compute_algebra_S`(alphas[, verbose])	Given operators α_i, compute the S matrix
`ff.ff_lib.correlation_matrix`(rho)	Calculates the following two-point correlation matrix
`ff.ff_lib.random_FF_rotation`(n_sites[, ...])	Generate a random free fermion rotation matrix
`ff.ff_lib.random_FF_state`(n_sites[, fixedN, ...])	Generate Haar random free fermion state using symplectic rotations.
`ff.ff_lib.random_H_generator`(n_sites[, ...])	Generate a random generator matrix for free fermions.
`ff.ff_lib.kitaev_chain`(n_sites, mu, t, delta)	Create Kitaev chain Hamiltonian.

Symplectic Operations

`ff.ff_lib.eigh_sp`(H)	Return the eigenvectors in symplectic form and the eigenvalues of a complex Hermitian (conjugate symmetric) FF coefficient array.
`ff.ff_lib.eigv_sp`(V)	Returns the left eigenvectors in symplectic form and eigenvalues in Y-form
`ff.ff_lib.eigm_sp_can`(G)	Returns the orthogonal eigenvectors and the eigenvalues in direct sum canonical form for the Majorana coefficient matrix G.
`ff.ff_lib.eigm_sp`(G)	Returns the orthogonal eigenvectors and the eigenvalues in L_Y form for the Majorana coefficient matrix G.
`ff.ff_lib.is_symp`(U)	Function checks if U is symplectic
`ff.ff_lib.check_canonical_form`(L)	Check if a matrix is in standard block diagonal form.

Quantum Physics Analysis Functions

ff.ff_lib.is_matchgate(M[, verbose])

Check the matchgate condition for 4x4 input matrices.

Combinatorial Functions

`ff.ff_combinatorics.sgn`(permutation)	Computes the sign of a permutation by counting inversions.
`ff.ff_combinatorics.pf`(A)	Computes pfaffian via combinatorial formula:
`ff.ff_combinatorics.hf`(A)	Computes hafnian via combinatorial formula.
`ff.ff_combinatorics.pt`(A)	Computes permanent via combinatorial formula.
`ff.ff_combinatorics.dt`(A)	Computes determinant via combinatorial formula.
`ff.ff_combinatorics.dt_eigen`(A)	Computes the determinant of a matrix using eigenvalue decomposition.

Graph Theory Functions

`ff.ff_graph_theory.pfo_algorithm`(graph[, ...])	Pfaffian ordering algorithm for planar graphs (FKT algorithm).
`ff.ff_graph_theory.plot_graph_with_edge_weights`(A)	Plot the graph associated with adjacency matrix A with edge weights.
`ff.ff_graph_theory.generate_random_planar_graph`(n)	Generate a random planar graph with n nodes.
`ff.ff_graph_theory.plot_planar_embedding`(G)	Plot the planar embedding of a graph with oriented faces.
`ff.ff_graph_theory.dual_graph_H`(G, F, T)	Create the dual graph according to the PFO algorithm (Vazirani 1987).
`ff.ff_graph_theory.faces`(graph[, emb])	Return the faces of an embedded graph using the networkx embedding scheme.
`ff.ff_graph_theory.complete_face`(pfo, edges)	Update one edge while maintaining the PFO condition.

Graph Visualization Functions

`ff.ff_graph_theory._draw_labeled_multigraph`(G)	Draw a labeled multigraph with curved edges for multiple connections.
`ff.ff_graph_theory._draw_labeled_graph`(G[, ...])	Draw a labeled graph using positions (if given).

Perfect Matching Functions

`ff.ff_graph_theory.count_perfect_matchings`(graph)	Find all perfect matchings in a graph.
`ff.ff_graph_theory.count_perfect_matchings_planar`(graph)	Find all perfect matchings in a planar graph using the PFO algorithm.
`ff.ff_graph_theory.find_perfect_matchings_brute`(graph)	Find all perfect matchings in a graph using a brute force enumeration.

Tree Analysis Functions

ff.ff_graph_theory.compute_tree_depth(graph)

Compute the depth of a tree graph.

Utility Functions

`ff.ff_utils._print`(obj[, k])	Printing with small number suppression (using numpy printoptions)
`ff.ff_utils.clean`(obj[, threshold])	Clean small numerical values from arrays or matrices.
`ff.ff_utils.formatted_output`(obj[, precision])	Format numerical output with specified precision.
`ff.ff_utils.generate_random_bitstring`(n, k)	Generates a random bit string of length n with Hamming weight k.
`ff.ff_utils.kron_plus`(a, b)	Computes the direct sum of two matrices
`ff.partial_trace_over_1`(AB, d)	This function implements the partial trace using a block matrix approach, where the composite matrix AB is viewed as a \(d \times d\) array of \(d_2 \times d_2\) blocks, and sum of the diagonal blocks gives the partial trace result.
`ff.partial_trace_over_2`(AB, d)	This function implements the partial trace using a block matrix approach, where the composite matrix AB is viewed as a \(d \times d\) array of \(d_1 \times d_1\) blocks, and the trace of each block gives the corresponding matrix element of the result matrix.