The ffsim simulation backend

Important

The concepts in this guide are currently available only in the Python API. Equivalent functionality will be made available through the C API in a future release.

The getting-started guides simulate fermionic circuits and Hamiltonians with ffsim – calling ffsim.apply_unitary(), ffsim.linear_operator(), or ffsim.sample_state_vector() directly on FermionicCircuits and FermionOperators. This guide explains why that works and what is actually happening underneath: this package couples deliberately with ffsim’s simulation protocols rather than building its own simulation API, and it does so in a way that keeps a native (scipy-only) simulation path available for users who cannot or do not want to install ffsim.

Why couple with ffsim

ffsim is a high-performance simulator for fermionic quantum circuits that exploits particle-number and spin-Z conservation to represent state vectors far more compactly than a generic \(2^n\)-dimensional qubit statevector. It defines two small protocols that any object can implement to participate in its simulation machinery (see qiskit_fermions.protocols for this package’s own protocols, which follow the same design):

Rather than invent a separate simulation interface, every fermionic gate in qiskit_fermions.circuit.library and FermionOperator implement these exact protocols. The direct payoff is that ffsim’s own tools work natively on this package’s objects, with no conversion step: ffsim.apply_unitary() can simulate a FermionicCircuit end to end, ffsim.linear_operator() (or plain scipy.sparse.linalg.eigsh()) can diagonalize a FermionOperator, and sampling utilities like ffsim.sample_state_vector() – used in the SKQD guide to turn a simulated statevector into measurement counts – work out of the box.

>>> import numpy as np
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution
>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>> from qiskit_fermions.utils.optionals import HAS_FFSIM
>>>
>>> if HAS_FFSIM:
...     import ffsim
>>>
>>> norb, nelec = 2, (1, 1)
>>> hamiltonian = FermionOperator.from_terms([
...     ([cre(0), ann(1)], 0.5),
...     ([cre(1), ann(0)], 0.5),
... ])
>>>
>>> circuit = FermionicCircuit(2 * norb)
>>> circuit.append(Evolution(2 * norb, hamiltonian, time=1.0), circuit.modes)
>>>
>>> if HAS_FFSIM:
...     reference = ffsim.hartree_fock_state(norb, nelec)
...     state = ffsim.apply_unitary(reference, circuit, norb=norb, nelec=nelec)  # native ffsim call

This is the concrete reason for coupling with ffsim’s protocols rather than, say, only offering a bespoke simulate() method: it lets ffsim’s existing, actively developed ecosystem of simulation and sampling utilities apply to this package’s circuits and operators unchanged, and it lets users already working with ffsim mix in this package’s gates and operators without learning a second simulation API.

A native path when ffsim is unavailable

Coupling with ffsim’s protocols does not make ffsim a hard dependency. ffsim transitively depends on PySCF, which does not support Windows – so ffsim is declared as an optional extra (pip install "qiskit-fermions[simulation]", or transitively via [all]), guarded at runtime by HAS_FFSIM (as you already saw above). On Windows, that extra simply resolves to nothing (a silent no-op via a sys_platform marker), rather than an install failure.

Simulation must still work without ffsim, so SupportsLinearOperator is backed by an independent native Rust FCI (full configuration interaction) kernel – not a wrapper around ffsim’s own linear-algebra routines. It compiles an operator’s terms once into a scatter map over the fixed-particle-number determinant basis, then reuses that compiled form across repeated matrix-vector products. This is exactly the protocol method ffsim itself expects (ffsim.SupportsLinearOperator): it happens to be implemented without ffsim, and it is cross-checked against ffsim’s own matrix elements in this package’s test suite, but it does not call into ffsim at all. Handing a FermionOperator to scipy.sparse.linalg.eigsh() or scipy.sparse.linalg.expm_multiply() – as Evolution._apply_unitary_placed_() does internally – therefore works identically whether or not ffsim is installed:

>>> import scipy.sparse.linalg
>>> from qiskit_fermions.linalg import linear_operator
>>>
>>> linop = linear_operator(hamiltonian, norb, nelec)  # pure scipy + native Rust kernel, no ffsim
>>> energy, _ = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> print(f"ground-state energy: {energy[0]:.6f}")
ground-state energy: -0.500000

Some gates go further and branch on HAS_FFSIM at simulate time purely for performance: OrbitalRotation is one example, whose _apply_unitary_ implementation delegates to ffsim’s dedicated Givens-rotation kernel (ffsim.apply_orbital_rotation()) when ffsim is installed, and otherwise falls back to expressing the rotation as \(\exp(G)\) for a one-body generator \(G\) and applying that through the same native-kernel-plus-expm_multiply() path used throughout this section. Both paths implement the same protocol method and produce the same result – ffsim is a performance choice here, not a correctness dependency. Gates without such a fast path (Evolution, and everything built out of it, such as UCJ) simply use the ffsim-independent path unconditionally.

In short: ffsim’s protocols are the interface; ffsim itself is an accelerator you can uninstall. This is also why the two getting-started guides that use ffsim directly (LUCJ and SKQD) internally guard their ffsim-specific code with HAS_FFSIM the same way this guide does.

Fermionic simulation lives in a fixed particle-number sector

Both _apply_unitary_ and _linear_operator_ take an nelec argument and represent the state vector over the fixed-particle-number determinant basis for that (norb, nelec) sector – mirroring ffsim’s (and, transitively, PySCF’s) FCI space setup, rather than the full \(2^{\text{num\_modes}}\)-dimensional space a general qubit simulator would use. This is a much smaller space (its size is a product of binomial coefficients rather than a power of two), but it comes with a hard restriction: only operators and gates that preserve particle number (and, in the spinful case, each spin species’ particle number individually) can be represented in it. An operator whose action would move amplitude to a different particle number simply has nowhere to go in this fixed-sector picture.

The native kernel resolves this by silently dropping any term that would leave the sector – projecting its contribution to zero rather than raising an error, since the kernel’s contract is “matrix-vector product on this sector,” not “validate this operator.” For most callers this dropping would be the wrong thing: applying \(\exp(-i t H)\) for a Hamiltonian \(H\) with a non-particle-conserving term would silently turn a would-be unitary into a non-unitary, physically-meaningless map. So the higher-level entry points that build a unitary out of the kernel add an explicit guard before ever calling it:

>>> from qiskit_fermions.circuit.library import Evolution
>>> from qiskit_fermions.linalg import apply_unitary
>>>
>>> non_conserving = FermionOperator.from_terms([([cre(0), cre(1)], 1.0)])  # creates 2 particles
>>> gate = Evolution(2 * norb, non_conserving, time=1.0)
>>>
>>> if not HAS_FFSIM:
...     # on Windows we manually define our reference state vector
...     reference = np.asarray([1, 0, 0, 0], dtype=complex)
>>>
>>> try:
...     apply_unitary(reference, gate, norb, nelec, copy=True)
... except ValueError as exc:
...     print("rejected:", exc)
rejected: Evolution requires an operator that conserves the (norb, nelec) sector: every term must preserve the particle number of each spin species (norb=2, nelec=(1, 1)).

Calling SupportsLinearOperator._linear_operator_() directly on a non-conserving operator bypasses this guard – it is a lower-level building block, not a validated simulation entry point – and will return a matrix-vector product that silently zeroes the non-conserving amplitude rather than raising.

The practical consequence is that non-particle-preserving simulation is not possible in fermionic space at all, by construction of the fixed-sector representation – not as a missing feature, but as the price of the compact FCI-space representation that makes fermionic simulation tractable in the first place. If your algorithm genuinely needs particle-number-violating operators (for example, a qubit-native error channel, or an operator built for a mapped Hamiltonian that does not conserve particle number term by term), transpile to qubits first and simulate the resulting QuantumCircuit with a qubit-level simulator instead, which has no such restriction. Any transpilation route works for this; the fermionic circuit and transpilation guides go into more detail, but generate_preset_jw_pass_manager() is a reasonable default choice to reach for.

The block-spin convention for spinful systems

nelec is typed int | tuple[int, int], and its type – not a separate flag – is what selects between the two supported mode layouts:

  • an int selects the spinless interpretation: the norb modes are treated directly as spinless orbitals, and nelec is simply the total particle count.

  • a pair (n_alpha, n_beta) selects the spinful interpretation of 2 * norb modes under a fixed block-spin convention: modes 0 .. norb are the alpha (spin-up) orbitals and modes norb .. 2 * norb are the beta (spin-down) orbitals. Each spin species is conserved independently – an alpha-only term can move an electron between alpha modes but never into a beta mode, and vice versa.

This dispatch happens at every ffsim-protocol entry point in this package – gate constructors, _apply_unitary_placed_ implementations, and the native Rust kernel’s own sector compilation – and it is exactly the convention used throughout the LUCJ and SKQD guides (for example, InitializeModes.from_hartree_fock() fills alpha occupations into modes 0 .. n_alpha and beta occupations into modes norb .. norb + n_beta). It is worth contrasting this with the qiskit_fermions.operators module and FermionicCircuit in general, which use generic mode indices with no inherent spin semantics (see the fermionic circuit guide) – the block-spin meaning is imposed only when a spinful nelec is supplied to a simulation call, not baked into the operator or circuit representation itself.

>>> from qiskit_fermions.circuit.library import InitializeModes
>>>
>>> norb, nelec = 3, (2, 1)
>>> init = InitializeModes.from_hartree_fock(norb, nelec)
>>> print([bool(occ) for occ in init.occupation])  # alpha modes 0,1; beta mode 3 (= norb + 0)
[True, True, False, True, False, False]

Next steps