<a id="operators-explanation"></a>

# Design principles of operator representations

This guide explains the common design principles and core concepts shared across
all operator representations in the [`operators`](../pydoc/qiskit_fermions.operators.html.md#module-qiskit_fermions.operators) module.

## Overview

The operator representations provided by this module share several fundamental
design principles:

- [**Sparse data structure**]():
  Operators encode only non-identity operations. The internal data layout is
  generally inspired by sparse matrix data formats, enabling efficient storage
  and computation for systems with many modes but relatively few significant
  contributions.
- [**Term iteration and reconstruction**]():
  Irrespective of the internal sparse storage, operators provide a consistent
  iteration interface, so you can inspect, filter, and transform terms without
  understanding the underlying data structure, then reconstruct new operators
  from modified terms.
- [**Mode-based indexing**]():
  Operators use abstract mode indices to label fermionic degrees of freedom,
  enabling flexible mapping from physical systems to the operator representation.
- [**Term grouping and commutation relationships**]():
  Operators natively support grouping information that associates terms with
  group indices. This enables optimizations and physical structure preservation
  without requiring separate data structures. See the [grouping
  guide](grouping.html.md#grouping-explanation) for practical usage.
- [**Arithmetic and mathematical operations**]():
  All operators implement a consistent set of arithmetic operations (such as
  addition, multiplication, and composition) and mathematical functions by using
  the [`OperatorTrait`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait) protocol, enabling uniform code across different
  operator types.
- [**Operator term ordering and normal forms**]():
  Mathematically equivalent operators can have very different representations and
  behavior in quantum algorithms. All operator representations support various
  normal forms (based on algebra-specific commutation relations) to achieve
  canonical, predictable, and optimizable operator representations.

<a id="sparse-term-representation"></a>

## Sparse data structure

All operators use a **sparse representation** where only non-identity operations
are tracked. Each non-identity operation consists of a **coefficient** (complex
number) and a **sequence of actions** on specific modes.

This approach dramatically reduces memory usage and computation time, especially for
systems with many modes but relatively few significant contributions. By encoding only
non-identity operations, operations focus only on what matters, enabling work with large
systems that would be infeasible with dense representations. Additionally, operators
naturally scale to any number of modes; an operator acting on modes `{0, 1}` works
unchanged in systems with many more modes since unaffected modes are implicitly identity.

Identical terms are preserved separately during arithmetic operations
and must be explicitly combined if needed.

### Internal storage format

Internally, operators are stored in arrays inspired by sparse matrix data formats:

- **Coefficients array**: The complex coefficient for each term
- **Mode indices array**: The fermionic modes that each action acts upon
- **Boundaries array**: Indices marking where each term’s modes begin and end in the mode array

#### IMPORTANT
Additional arrays might be present depending on the operator type. For example,
[`FermionOperator`](../stubs/qiskit_fermions.operators.FermionOperator.html.md#qiskit_fermions.operators.FermionOperator) instances include an **actions array** of booleans
that specifies the type of fermionic action acting on the respective mode
index. In contrast, the `MajoranaOperator` class does not require
this distinction since it encodes that information in the parity of the mode
index. See the API documentation for your operator type to
understand the full storage format.

The following examples show how these arrays are organized. First is a direct
construction using the sparse arrays:

### PYTHON

```python
>>> from qiskit_fermions.operators import FermionOperator
>>>
>>> # Construct operators directly using sparse arrays
>>> # First operator: 1.0 * +0 -1
>>> op1 = FermionOperator(
...     coeffs=[1.0],
...     actions=[True, False],
...     modes=[0, 1],
...     boundaries=[0, 2],
... )
>>>
>>> # Second operator: 1.0 * +2 -3
>>> op2 = FermionOperator(
...     coeffs=[1.0],
...     actions=[True, False],
...     modes=[2, 3],
...     boundaries=[0, 2],
... )
>>>
>>> # Combine the sparse operators
>>> op1 += op2
>>> print(format(op1))
  1.000000e0 +0.000000e0j * (+0 -1)
  1.000000e0 +0.000000e0j * (+2 -3)
```

### C

```c
#include <qiskit_fermions.h>

// Construct first operator: 1.0 * +0 -1
QkComplex64 coeff1[1] = {{1.0, 0.0}};
bool actions1[2] = {true, false};
uint32_t modes1[2] = {0, 1};
uint32_t boundaries1[2] = {0, 2};
QfFermionOperator *op1 = qf_ferm_op_new(1, 2, coeff1, actions1, modes1, boundaries1);

// Construct second operator: 1.0 * +2 -3
QkComplex64 coeff2[1] = {{1.0, 0.0}};
bool actions2[2] = {true, false};
uint32_t modes2[2] = {2, 3};
uint32_t boundaries2[2] = {0, 2};
QfFermionOperator *op2 = qf_ferm_op_new(1, 2, coeff2, actions2, modes2, boundaries2);

// Combine the sparse operators. There is no in-place add in the C API, so
// the sum is a new operator that the caller owns.
QfFermionOperator *sum = qf_ferm_op_add(op1, op2);

qf_ferm_op_free(op1);
qf_ferm_op_free(op2);
qf_ferm_op_free(sum);
```

### Convenient construction methods

For Python developers, several convenient construction methods are available that
abstract away the sparse storage details. These make it easier to build operators
without worrying about managing coefficient, mode, and boundary arrays:

### PYTHON

```python
>>> from qiskit_fermions.operators import FermionOperator, cre, ann
>>>
>>> # Construct operators using operator algebra notation
>>> op1 = FermionOperator.from_dict({(cre(0), ann(1)): 1.0})
>>> op2 = FermionOperator.from_dict({(cre(2), ann(3)): 1.0})
>>>
>>> # The result is sparse even when combining them
>>> op1 += op2
>>> print(format(op1))
1.000000e0 +0.000000e0j * (+0 -1)
1.000000e0 +0.000000e0j * (+2 -3)
```

### C

```c
// The C API uses direct array construction; convenience methods are not available.
```

#### HINT
Individual operator implementations might support additional construction methods
suited to their specific use case. Check the API documentation for your operator
type to see all available construction options.

<a id="term-iteration-and-reconstruction"></a>

## Term iteration and reconstruction

Operators provide a consistent iteration interface by using [`OperatorTrait.iter_terms()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.iter_terms)
irrespective of their internal sparse representation. You can therefore inspect, filter,
or transform terms without needing to understand the underlying data structure.
You can then reconstruct a new operator from the transformed terms by using
[`OperatorTrait.from_terms()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.from_terms).

### PYTHON

```python
>>> from qiskit_fermions.operators import FermionOperator, cre, ann
>>>
>>> # Construct an operator with terms of different orders
>>> op = FermionOperator.from_dict({
...     (): 0.5,  # constant term (order 0)
...     (cre(0), ann(1)): 1.0,  # two-body term (order 2)
...     (cre(0), cre(1), ann(1), ann(0)): 0.25  # four-body term (order 4)
... })
>>>
>>> # Filter to keep only terms of order 2
>>> order_two_terms = [
...     (term, coeff) for term, coeff in op.iter_terms()
...     if len(term) == 2
... ]
>>>
>>> # Reconstruct operator from filtered terms
>>> filtered_op = FermionOperator.from_terms(order_two_terms)
>>> print(f"Original operator has {len(op)} terms")
Original operator has 3 terms
>>> print(f"Filtered operator has {len(filtered_op)} term")
Filtered operator has 1 term
```

### C

```c
// WARNING: Term iteration and filtering are not yet available in the C API.
```

<a id="mode-based-indexing"></a>

## Mode-based indexing

All operator representations refer to the indices that their terms act upon as
**modes**. A mode is simply an index identifying a fermionic degree of freedom
in your system. The mapping from physical degrees of freedom (such as spatial
orbitals, spin states, or other quantum numbers) to mode indices is left to the
user, enabling maximum flexibility.

This abstraction is also present in the [`qiskit_fermions.circuit`](../pydoc/qiskit_fermions.circuit.html.md#module-qiskit_fermions.circuit) module,
where the [`FermionicCircuit`](../stubs/qiskit_fermions.circuit.FermionicCircuit.html.md#qiskit_fermions.circuit.FermionicCircuit) operates on a register of fermionic modes.
In both the operator and circuit representations, modes provide a consistent,
algebra-independent way to specify which degrees of freedom participate in a
given operation.

#### IMPORTANT
**Current implementations use spinless modes**: All operator representations
currently provided by this module treat modes as spinless fermionic degrees
of freedom. This means if your system has both spin-up and spin-down electrons
or fermions, you must explicitly map them to distinct modes (for example, modes 0-3
for spin-up of four spatial orbitals, modes 4-7 for spin-down, or any other
convention you choose).

This design keeps the core representations simple and general while avoiding
imposing a specific spin-ordering convention. Utility modules like
[`qiskit_fermions.operators.library`](../pydoc/qiskit_fermions.operators.library.html.md#module-qiskit_fermions.operators.library) provide convenience functions
(for example, [`FCIDump.from_file()`](../stubs/qiskit_fermions.operators.library.FCIDump.html.md#qiskit_fermions.operators.library.FCIDump.from_file)) that handle such mappings for you when
loading electronic structure data.

#### HINT
As the package evolves, spinful operator representations might be added that
natively support spin degrees of freedom in their data model. These will be
clearly distinguished from the current spinless implementations and will
coexist with them in the module.

<a id="term-grouping"></a>

## Term grouping and commutation relationships

Like the coefficients and mode indices, operators can optionally store a **groups array**
that associates each term with a group index. By integrating grouping directly into the
operator representation as part of the sparse data structure, grouping information
naturally travels with the operator through transformations. This enables systematic
structure usage, whether from physical properties, algebraic relationships, or
problem-specific symmetries. The structured information can then be used in downstream
operations like circuit synthesis and decomposition by using methods like
[`OperatorTrait.split_out_groups()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.split_out_groups).

For detailed guidance on how to group operator terms in your workflows, see
the [grouping](grouping.html.md#grouping-explanation) guide.

### PYTHON

```python
>>> from qiskit_fermions.operators import MajoranaOperator, gamma
>>> op = MajoranaOperator.from_dict({
...     (gamma(0, False),): 1.0,
...     (gamma(1, False),): 1.0,
...     (gamma(2, False), gamma(3, False)): 1.0
... })
>>> # Assign group indices to terms
>>> op.groups = [0, 0, 1]
>>> # Partition operator by groups
>>> grouped_ops = op.split_out_groups()
```

### C

```c
#include <qiskit_fermions.h>

// Create operator with 3 terms
QkComplex64 coeffs[3] = {{1.0, 0.0}, {1.0, 0.0}, {1.0, 0.0}};
uint32_t modes[4] = {0, 1, 2, 3};
uint32_t boundaries[4] = {0, 1, 2, 4};
QfMajoranaOperator *op = qf_maj_op_new(3, 4, coeffs, modes, boundaries);

// Assign grouping information
uint32_t groups[3] = {0, 0, 1};
qf_maj_op_set_groups(op, groups, 3);

// Partition operator by groups
QfMajoranaOperator *grouped_ops[2];
qf_maj_op_split_out_groups(op, NULL, 0, grouped_ops);
```

<a id="arithmetic-and-mathematical-operations"></a>

## Arithmetic and mathematical operations

All operators implement the [`OperatorTrait`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait) protocol, which provides
a unified set of operations across different operator types. This ensures that code
written for one operator representation works uniformly with others. See
[`qiskit_fermions.protocols`](../pydoc/qiskit_fermions.protocols.html.md#module-qiskit_fermions.protocols) for how this protocol relates to the others this
package defines.

The protocol includes arithmetic operations (such as addition, multiplication,
and composition), structural operations (term iteration, mode support analysis,
relabeling), mathematical functions (normal ordering, simplification, equivalence
checking), and more. See the [`OperatorTrait`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait) documentation for a complete
reference of all available operations.

### PYTHON

```python
>>> from qiskit_fermions.operators import FermionOperator, cre, ann
>>>
>>> # Construct a Hermitian operator: H = +0 -1 + +1 -0
>>> op = FermionOperator.from_dict({
...     (cre(0), ann(1)): 1.0,
...     (cre(1), ann(0)): 1.0
... })
>>>
>>> # Check if the operator is Hermitian
>>> print(f"Operator is Hermitian: {op.is_hermitian(atol=1e-10)}")
Operator is Hermitian: True
```

### C

```c
#include <qiskit_fermions.h>

// Construct a Hermitian operator: H = +0 -1 + +1 -0
QkComplex64 coeffs[2] = {{1.0, 0.0}, {1.0, 0.0}};
bool actions[4] = {true, false, true, false};
uint32_t modes[4] = {0, 1, 1, 0};
uint32_t boundaries[3] = {0, 2, 4};
QfFermionOperator *op = qf_ferm_op_new(2, 4, coeffs, actions, modes, boundaries);

// Check if the operator is Hermitian
bool is_hermitian = qf_ferm_op_is_hermitian(op, 1e-10);
printf("Operator is Hermitian: %s\n", is_hermitian ? "true" : "false");

// Clean up
qf_ferm_op_free(op);
```

#### NOTE
[`is_hermitian()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.is_hermitian) verifies that $H - H^\dagger = 0$ by reducing that
difference to a normal form and comparing it against [`zero()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.zero). Spelled out
in terms of the other protocol operations, it performs the following:

```python
>>> difference = (op - op.adjoint()).normal_ordered().simplify(atol=1e-10)
>>> difference.equiv(FermionOperator.zero(), atol=1e-10)
True
```

The strength of this check depends on the operator type, because it depends on how completely
that type’s normal form contracts adjacent generators. A `True` result is always reliable, but
for some types a `False` result is conservative: see
[`is_hermitian()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.is_hermitian) for details.

#### IMPORTANT
The example uses `atol=1e-10`. The `atol` (absolute tolerance) parameter specifies a
threshold. Coefficients with magnitude smaller than `atol` are treated as zero and discarded.
This is essential for numerical stability when comparing operators, since floating-point
arithmetic can introduce small rounding errors that would otherwise prevent equivalent operators
from being recognized as such.

#### HINT
While the [`OperatorTrait`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait) protocol provides a common interface,
individual operator implementations might offer additional convenience methods
not part of the protocol. For example, [`FermionOperator`](../stubs/qiskit_fermions.operators.FermionOperator.html.md#qiskit_fermions.operators.FermionOperator) and
[`MajoranaOperator`](../stubs/qiskit_fermions.operators.MajoranaOperator.html.md#qiskit_fermions.operators.MajoranaOperator) provide a `max_rank()` method that reports the
length of their longest term. Always consult the API documentation for your
operator type to discover all available functionality.

<a id="operator-term-ordering"></a>

## Operator term ordering and normal forms

A fundamental challenge in quantum operator algebra is that mathematically
equivalent operators can be represented in many different ways, each with
different implications for quantum algorithms. The same operator can be written
in algebraically equivalent forms (for example, $a^\dagger b$ can be
expressed as $ba^\dagger + [a^\dagger,b]$), yet these representations lead
to different behavior in circuit synthesis, simplification, and numerical
algorithms.

The operator representations support **normal ordering** operations that
transform operators into algebra-specific canonical forms using commutation
relations. This enables reliable comparisons (two equivalent operators have
identical normal-ordered forms), reveals simplifications (commutation relations
cause terms to cancel or combine), and supports algorithms that require specific
operator forms for correctness and efficiency.

#### IMPORTANT
Different operator representations might implement normal ordering based on
different commutation relations appropriate to their algebra. For example,
fermionic normal ordering uses anticommutation relations ($\{c_i, c_j^\dagger\} = \delta_{ij}$),
while Majorana normal ordering uses different algebra conventions
($\{\gamma_i, \gamma_j\} = 2\delta_{ij}$). Always consult the
documentation for your operator type to understand how normal
ordering is implemented.

### PYTHON

```python
>>> from qiskit_fermions.operators import FermionOperator, cre, ann
>>>
>>> # Two different representations of the same operator
>>> op1 = FermionOperator.from_dict({(ann(0), cre(0)): 1.0})
>>> op2 = FermionOperator.from_dict({(): 1.0, (cre(0), ann(0)): -1.0})
>>>
>>> # Direct comparison fails due to different forms
>>> op1.equiv(op2, atol=1e-10)
False
>>>
>>> # Normal-order both and compare again
>>> op1_normal = op1.normal_ordered()
>>> op2_normal = op2.normal_ordered()
>>> op1_normal.equiv(op2_normal, atol=1e-10)
True
```

### C

```c
#include <qiskit_fermions.h>
#include <stdbool.h>

// Two different representations of the same operator
QkComplex64 coeff1[1] = {{1.0, 0.0}};
bool actions1[2] = {false, true};
uint32_t modes1[2] = {0, 0};
uint32_t boundaries1[2] = {0, 2};
QfFermionOperator *op1 = qf_ferm_op_new(1, 2, coeff1, actions1, modes1, boundaries1);

QkComplex64 coeff2[2] = {{1.0, 0.0}, {-1.0, 0.0}};
bool actions2[2] = {true, false};
uint32_t modes2[2] = {0, 0};
uint32_t boundaries2[3] = {0, 0, 2};
QfFermionOperator *op2 = qf_ferm_op_new(2, 2, coeff2, actions2, modes2, boundaries2);

// Direct comparison fails due to different forms
bool equiv_before = qf_ferm_op_equiv(op1, op2, 1e-10);
printf("Equivalent before normal ordering: %s\n", equiv_before ? "true" : "false");

// Normal-order both and compare again
QfFermionOperator *op1_normal = qf_ferm_op_normal_ordered(op1, NULL);
QfFermionOperator *op2_normal = qf_ferm_op_normal_ordered(op2, NULL);
bool equiv_after = qf_ferm_op_equiv(op1_normal, op2_normal, 1e-10);
printf("Equivalent after normal ordering: %s\n", equiv_after ? "true" : "false");

// Clean up
qf_ferm_op_free(op1);
qf_ferm_op_free(op2);
qf_ferm_op_free(op1_normal);
qf_ferm_op_free(op2_normal);
```

#### HINT
The [`OperatorTrait.normal_ordered()`](../stubs/qiskit_fermions.operators.OperatorTrait.html.md#qiskit_fermions.operators.OperatorTrait.normal_ordered) protocol method leaves its positional
and keyword arguments deliberately unspecified, allowing concrete operator
implementations to define adjustable parameters that control the exact canonical
form produced. This enables operator-specific optimizations and normal form
variants tailored to their algebra or use case. Check the API documentation for
your operator type to see what parameters are available.