From dbd4bf70554a6c3681c4d87dd19aa7f856f6d193 Mon Sep 17 00:00:00 2001
From: Adam <aglos@iitis.pl>
Date: Mon, 27 Jul 2026 17:04:06 +0000
Subject: [PATCH 1/6] fix: :bug: Making number of qubits obligatory in order to
 avoid confusion in the interface

---
 benches/_builders.py                |  1 +
 src/monoprop/circuit.py             | 66 +++++++++++++++++++++++------
 src/monoprop/fermi.py               | 25 ++++++-----
 src/monoprop/majorana.py            | 20 +++++++--
 src/monoprop/monomial_propagator.py | 20 +++++++--
 src/monoprop/pauli.py               | 33 ++++++---------
 src/monoprop/qiskit_conversion.py   | 10 ++++-
 src/monoprop/utils.py               | 23 ++++++++++
 tests/cases.py                      |  3 ++
 tests/test_circuit.py               |  4 +-
 tests/test_coeff_trunc.py           |  3 ++
 tests/test_majorana.py              |  1 +
 tests/test_nonfermi.py              |  1 +
 tests/test_only_rotate_k.py         |  1 +
 tests/test_update_methods.py        |  1 +
 15 files changed, 157 insertions(+), 55 deletions(-)

diff --git a/benches/_builders.py b/benches/_builders.py
index 22850d1f..762dca75 100644
--- a/benches/_builders.py
+++ b/benches/_builders.py
@@ -167,6 +167,7 @@ def make_random_problem(
         majoranas=gen_majoranas,
         gen_coeffs=gen_coeffs,
         param_inds=param_inds,
+        num_modes=num_modes,
         parameters=parameters,
         initial_state=[],
     )
diff --git a/src/monoprop/circuit.py b/src/monoprop/circuit.py
index 10cda2e7..4b6e92cc 100644
--- a/src/monoprop/circuit.py
+++ b/src/monoprop/circuit.py
@@ -57,6 +57,7 @@
 from .conversion_utils import _extend_pauli_string, _pauli_to_majorana
 from .majorana import MajoranaOperator
 from .pauli import Pauli, PauliOperator
+from .utils import _validate_system_size
 
 if TYPE_CHECKING:
     from collections.abc import Iterator, Sequence
@@ -255,33 +256,42 @@ class Circuit:
         gates: The ordered exponential gates.
         parameters: The angle values, or empty for an unbound circuit.
         initial_state: The reference state (occupied mode / qubit indices).
+        num_modes: System width (number of fermionic modes / qubits).
         family: The gate family -- ``"pauli"``, ``"majorana"``, or ``"empty"`` -- computed at
             construction; the propagators dispatch on it.
     """
 
-    def __init__(
+    def __init__(  # noqa:C901
         self,
-        gates: Sequence[ExpGate] = (),
+        gates: Sequence[ExpGate],
+        initial_state: Sequence[int],
+        num_modes: int,
         parameters: Sequence[float] = (),
-        initial_state: Sequence[int] = (),
     ) -> None:
         """Build the circuit, dropping identity gates and validating family/mapping/params.
 
         Args:
             gates: The ordered exponential gates.
-            parameters: The angle values, or empty for an unbound circuit.
+            num_modes: Number of modes/qubits for the circuit.
             initial_state: The reference state (occupied mode / qubit indices).
+            parameters: The angle values, or empty for an unbound circuit.
 
         Raises:
             ValueError: On duplicate initial-state indices, a bad parameter mapping, or a
-                bound circuit whose parameter count does not match [n_parameters][].
+                bound circuit whose parameter count does not match [n_parameters][]; also if
+                a gate's operator width differs from ``num_modes``.
             TypeError: On a non-[ExpGate][] gate or a mix of qubit and Majorana gate families.
         """
         gates = tuple(gates)
         parameters = tuple(float(v) for v in parameters)
         initial_state = tuple(int(i) for i in initial_state)
+        num_modes = _validate_system_size(num_modes, argument_name="num_modes")
         if len(set(initial_state)) != len(initial_state):
             raise ValueError("Duplicate indices in initial state")
+        if any(i < 0 or i >= num_modes for i in initial_state):
+            raise ValueError(
+                f"initial_state entries must be in 0..{num_modes - 1}; got {list(initial_state)}."
+            )
 
         # Validate gate types up front: the identity-drop below reads gate.index/.generator,
         # so a non-ExpGate gate must be rejected with a clear TypeError first rather than crashing
@@ -291,6 +301,15 @@ def __init__(
                 raise TypeError(
                     f"Circuit gates must be ExpGate; got {type(gate).__name__}."
                 )
+            gate_size = (
+                gate.generator.num_qubits
+                if isinstance(gate.generator, PauliOperator)
+                else gate.generator.num_modes
+            )
+            if gate_size != num_modes:
+                raise ValueError(
+                    f"Gate generator width {gate_size} does not match circuit num_modes={num_modes}."
+                )
 
         def _is_identity_gate(gate: ExpGate) -> bool:
             return all(coeff == 0 for coeff in gate.generator.terms.values())
@@ -310,6 +329,7 @@ def _is_identity_gate(gate: ExpGate) -> bool:
         self.gates = gates
         self.parameters = parameters
         self.initial_state = initial_state
+        self.num_modes = num_modes
         #: The gate family, computed from the (validated) gates; the propagators dispatch on it.
         self.family = self._resolve_family(gates)
 
@@ -328,6 +348,7 @@ def __eq__(self, other: object) -> bool:
             self.gates == other.gates
             and self.parameters == other.parameters
             and self.initial_state == other.initial_state
+            and self.num_modes == other.num_modes
         )
 
     __hash__ = None  # type: ignore[assignment]  # value-equal but not hashable (mutable gates)
@@ -336,7 +357,8 @@ def __repr__(self) -> str:
         """Return a string representation listing the gates, parameters, and initial state."""
         return (
             f"{self.__class__.__name__}(gates={self.gates!r}, "
-            f"parameters={self.parameters!r}, initial_state={self.initial_state!r})"
+            f"parameters={self.parameters!r}, initial_state={self.initial_state!r}, "
+            f"num_modes={self.num_modes!r})"
         )
 
     @staticmethod
@@ -420,6 +442,11 @@ def __add__(self, other: Circuit) -> Circuit:
             raise ValueError(
                 "Cannot concatenate circuits with different initial states."
             )
+        if self.num_modes != other.num_modes:
+            raise ValueError(
+                f"Cannot concatenate circuits with different num_modes: "
+                f"{self.num_modes} != {other.num_modes}."
+            )
         offset = self.n_parameters
         # Preserve each gate's _structural flag: a dense (wire-format) gate carries structural
         # coefficients that must not be antihermitian-normalized again.
@@ -435,6 +462,7 @@ def __add__(self, other: Circuit) -> Circuit:
             gates=left + right,
             parameters=tuple(self.parameters) + tuple(other.parameters),
             initial_state=self.initial_state or other.initial_state,
+            num_modes=self.num_modes,
         )
 
     @classmethod
@@ -443,6 +471,7 @@ def from_dense_arrays(
         majoranas: Sequence[Sequence[int]],
         gen_coeffs: Sequence[float],
         param_inds: Sequence[int],
+        num_modes: int,
         parameters: Sequence[float] = (),
         initial_state: Sequence[int] = (),
     ) -> Circuit:
@@ -460,6 +489,7 @@ def from_dense_arrays(
             gen_coeffs: Generator coefficient per monomial.
             param_inds: Variational-angle index per monomial (contiguous runs group into
                 gates).
+            num_modes: Number of fermionic modes in the system.
             parameters: Optional angle values.
             initial_state: Optional reference state (occupied mode indices).
 
@@ -467,6 +497,7 @@ def from_dense_arrays(
             A [Circuit][] carrying the grouped gates, angle values, and initial state.
         """
         indices = [int(p) for p in param_inds]
+        num_modes = _validate_system_size(num_modes, argument_name="num_modes")
         gates: list[ExpGate] = []
         current_index: int | None = None
         current_majoranas: list[tuple[int, ...]] = []
@@ -478,7 +509,7 @@ def _flush() -> None:
             gates.append(
                 ExpGate._structural_gate(
                     MajoranaOperator._from_terms(
-                        current_majoranas, current_coeffs, num_modes=None
+                        current_majoranas, current_coeffs, num_modes=num_modes
                     ),
                     index=current_index,
                 )
@@ -490,7 +521,16 @@ def _flush() -> None:
                 current_majoranas = []
                 current_coeffs = []
             current_index = pidx
-            current_majoranas.append(tuple(int(i) for i in maj))
+            majorana = tuple(int(i) for i in maj)
+            if any(i < 0 for i in majorana):
+                raise ValueError(
+                    f"Majorana indices must be non-negative; got {majorana}."
+                )
+            if majorana and max(majorana) >= 2 * num_modes:
+                raise ValueError(
+                    f"Majorana term {majorana} acts on an index >= 2*num_modes={2 * num_modes}."
+                )
+            current_majoranas.append(majorana)
             current_coeffs.append(complex(float(coeff)))
         if current_majoranas:
             _flush()
@@ -499,6 +539,7 @@ def _flush() -> None:
             gates=tuple(gates),
             parameters=tuple(float(p) for p in parameters),
             initial_state=tuple(int(i) for i in initial_state),
+            num_modes=num_modes,
         )
 
 
@@ -632,9 +673,7 @@ def _validate_commuting_majorana_generator(generator: MajoranaOperator) -> None:
             )
 
 
-def _gate_layers(
-    gate: ExpGate, num_qubits: int | None
-) -> list[tuple[tuple[int, ...], float]]:
+def _gate_layers(gate: ExpGate, num_qubits: int) -> list[tuple[tuple[int, ...], float]]:
     r"""Expand one gate into ``(majorana, gen_coeff)`` layers, in application order.
 
     A ``"pauli"``-family [ExpGate][] places each [Pauli][monoprop.pauli.Pauli] term on
@@ -649,8 +688,6 @@ def _gate_layers(
     # ``isinstance`` narrows the fall-through arm to MajoranaOperator).
     generator = gate.generator
     if isinstance(generator, PauliOperator):
-        if num_qubits is None:
-            raise ValueError("num_qubits is required to expand a Pauli gate.")
         layers: list[tuple[tuple[int, ...], float]] = []
         for pauli, coeff in generator.terms.items():
             extended = _extend_pauli_string(pauli.string, pauli.qubits, num_qubits)
@@ -673,7 +710,7 @@ def _gate_layers(
 def expand_monomials(
     gates: Sequence[ExpGate],
     mapping: Sequence[int],
-    num_qubits: int | None = None,
+    num_qubits: int,
 ) -> tuple[list[tuple[int, ...]], list[float], list[int], list[int]]:
     """Flatten gates + an already-resolved per-gate mapping into per-monomial arrays.
 
@@ -689,6 +726,7 @@ def expand_monomials(
         the authoring gate monomial ``i`` came from, so the engine can recover gate
         boundaries; monomials from a multi-term gate share one gate index.
     """
+    num_qubits = _validate_system_size(num_qubits, argument_name="num_qubits")
     majoranas: list[tuple[int, ...]] = []
     gen_coeffs: list[float] = []
     per_monomial: list[int] = []
diff --git a/src/monoprop/fermi.py b/src/monoprop/fermi.py
index 79a01c75..9058d92f 100644
--- a/src/monoprop/fermi.py
+++ b/src/monoprop/fermi.py
@@ -22,6 +22,7 @@
 
 from .conversion_utils import _n_product
 from .majorana import MajoranaOperator
+from .utils import _validate_system_size
 
 if TYPE_CHECKING:
     from collections.abc import Sequence
@@ -109,31 +110,35 @@ def __init__(
         self,
         terms: Sequence[FermiString] | Sequence[Sequence[tuple[int, str]]],
         coefficients: Sequence[complex],
-        num_modes: int | None = None,
+        num_modes: int,
     ) -> None:
         """Initialize the fermi operator.
 
         Args:
             terms: List of FermiString objects representing the operator.
             coefficients: List of coefficients corresponding to the terms.
-            num_modes: Optional number of modes. If not provided, it will be inferred from the terms.
+            num_modes: Number of modes in the system.
 
         Raises:
-            ValueError: If any index is out of bounds or if there are duplicate indices.
+            TypeError: If ``num_modes`` is not an integer.
+            ValueError: If ``num_modes`` is negative or any index is out of bounds.
         """
         self.terms = [
             t if isinstance(t, FermiString) else FermiString(t) for t in terms
         ]
         self.coefficients = list(coefficients)
-        self.num_modes = (
-            num_modes
-            if num_modes is not None
-            else max((idx for f in self.terms for idx, _ in f.expression)) + 1
-        )
+        self.num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        for term in self.terms:
+            for idx, _ in term.expression:
+                if idx >= self.num_modes:
+                    raise ValueError(
+                        "Fermi term index out of bounds: "
+                        f"{idx} >= num_modes={self.num_modes}."
+                    )
 
     @classmethod
     def from_dict(
-        cls, terms_dict: dict[tuple[tuple[int, str], ...], complex]
+        cls, terms_dict: dict[tuple[tuple[int, str], ...], complex], num_modes: int
     ) -> FermiOperator:
         """Construct a FermiOperator from a dictionary."""
         terms = []
@@ -141,7 +146,7 @@ def from_dict(
         for key, value in terms_dict.items():
             terms.append(FermiString(key))
             coefficients.append(value)
-        return cls(terms=terms, coefficients=coefficients)
+        return cls(terms=terms, coefficients=coefficients, num_modes=num_modes)
 
     def __len__(self) -> int:
         """Number of terms in the operator."""
diff --git a/src/monoprop/majorana.py b/src/monoprop/majorana.py
index 4291b550..839cacfd 100644
--- a/src/monoprop/majorana.py
+++ b/src/monoprop/majorana.py
@@ -21,6 +21,8 @@
 
 import numpy as np
 
+from .utils import _validate_system_size
+
 if TYPE_CHECKING:
     from collections.abc import Mapping, Sequence
 
@@ -99,6 +101,10 @@ def __init__(
                 also authored as a [MajoranaOperator][] (wrapped in
                 [ExpGate][monoprop.circuit.ExpGate]) -- bare [Majorana][] terms are not accepted
                 by ``ExpGate``, since the operator is what carries the mode count.
+
+        Raises:
+            TypeError: If ``num_modes`` is not an integer.
+            ValueError: If ``num_modes`` is negative or a term index is out of range.
         """
         # Route raw index tuples through Majorana so they get the same non-negative/distinct
         # validation a Majorana key already carries (a bare tuple would otherwise slip past it).
@@ -106,7 +112,13 @@ def __init__(
             (key if isinstance(key, Majorana) else Majorana(*key)).indices
             for key in terms
         ]
-        self.num_modes = num_modes
+        self.num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        for majorana in majoranas:
+            # majoranas are sorted in here
+            if majorana and majorana[-1] >= 2 * self.num_modes:
+                raise ValueError(
+                    f"Majorana term {majorana} acts on an index >= num_modes={self.num_modes}."
+                )
         self.terms = self._accumulate(majoranas, list(terms.values()))
 
     @classmethod
@@ -114,7 +126,7 @@ def _from_terms(
         cls,
         majoranas: Sequence[Sequence[int]],
         coefficients: Sequence[complex],
-        num_modes: int | None = None,
+        num_modes: int,
     ) -> MajoranaOperator:
         """Build from parallel ``majoranas``/``coefficients`` lists (internal).
 
@@ -122,7 +134,9 @@ def _from_terms(
         Jordan-Wigner and fermionic conversions ([get_majorana_operator][]) rely on.
         """
         obj = cls.__new__(cls)
-        obj.num_modes = num_modes
+        obj.num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        # _accumulate is taking care of the sorting/removing repeating terms
+        # so no need to check in here
         obj.terms = cls._accumulate(majoranas, coefficients)
         return obj
 
diff --git a/src/monoprop/monomial_propagator.py b/src/monoprop/monomial_propagator.py
index bb2c2a94..35f29a23 100644
--- a/src/monoprop/monomial_propagator.py
+++ b/src/monoprop/monomial_propagator.py
@@ -78,7 +78,7 @@ class MonomialPropagator(ABC):
 
     _comm: MPI.Comm | None
     _n_params: int
-    _num_qubits: int | None
+    _num_modes: int
     _initial_state: list[int]
     _simulator: object
 
@@ -115,6 +115,7 @@ def _init_simulator(
 
         self._comm = comm
         self._n_params = 0
+        self._num_modes = num_modes
         # System qubit count for expanding Pauli gates; set by PauliPropagator from the
         # observable. None for a native Majorana propagator (its gates need no qubit count).
         self._num_qubits = None
@@ -195,6 +196,15 @@ def _check_initial_state(self, circuit: Circuit) -> None:
                 "propagator with this circuit's initial state (or via from_circuit)."
             )
 
+    def _check_circuit_width(self, circuit: Circuit) -> None:
+        """Reject a circuit with a system width that disagrees with the propagator."""
+        expected = self._num_qubits if self._num_qubits is not None else self._num_modes
+        if circuit.num_modes != expected:
+            raise ValueError(
+                f"Circuit num_modes={circuit.num_modes} does not match propagator width "
+                f"{expected}."
+            )
+
     def _validate_and_correct_only_rotate_len_k(
         self, only_rotate_len_k: int | None
     ) -> int:
@@ -255,6 +265,7 @@ def build_graph(
                 expectation-value estimation in Schrodinger-picture simulations.
         """
         self._check_initial_state(circuit)
+        self._check_circuit_width(circuit)
         # Resolve the coefficient seed handed to the engine (the operator coefficients the new
         # layers are contracted against while the graph is built, informing coefficient
         # truncation). The engine validates its length against the accumulated parameter axis.
@@ -275,7 +286,7 @@ def build_graph(
         else:
             seed = None
         gates = self._circuit_gates(circuit)
-        num_qubits = self._num_qubits
+        num_qubits = 0 if self._num_qubits is None else self._num_qubits
         # Shift the circuit's local 0-based angle indices onto the accumulated axis.
         mapping = [self._n_params + m for m in circuit.resolved_mapping]
         self._n_params += circuit.n_parameters
@@ -312,8 +323,9 @@ def propagate(
             only_rotate_len_k
         )
         self._check_initial_state(circuit)
+        self._check_circuit_width(circuit)
         gates = self._circuit_gates(circuit)
-        num_qubits = self._num_qubits
+        num_qubits = 0 if self._num_qubits is None else self._num_qubits
         majoranas, gen_coeffs, mapping, _gate_indices = expand_monomials(
             gates, circuit.resolved_mapping, num_qubits
         )
@@ -585,7 +597,7 @@ def graph_size(self) -> tuple[int, int]:
     @property
     def num_modes(self) -> int:
         """Number of fermionic modes for the simulator."""
-        return self._simulator.num_modes
+        return self._num_modes
 
     @property
     def graph_layers(self) -> int:
diff --git a/src/monoprop/pauli.py b/src/monoprop/pauli.py
index ae136fab..ccb0bd92 100644
--- a/src/monoprop/pauli.py
+++ b/src/monoprop/pauli.py
@@ -23,6 +23,7 @@
 
 from .conversion_utils import _extend_pauli_string, _pauli_to_majorana
 from .majorana import MajoranaOperator
+from .utils import _validate_system_size
 
 if TYPE_CHECKING:
     from collections.abc import Mapping, Sequence
@@ -116,7 +117,7 @@ class PauliOperator:
     def __init__(
         self,
         terms: Mapping[Pauli | str, complex],
-        num_qubits: int | None,
+        num_qubits: int,
     ) -> None:
         """Initialize the Pauli operator from a term mapping.
 
@@ -125,11 +126,10 @@ def __init__(
                 coefficients.
             num_qubits: Total number of qubits the operator acts on. An operator carries its
                 own qubit count so a propagator can be built from it directly; every term must
-                act within ``0..num_qubits-1``. ``None`` defers the qubit count (only reachable
-                via `_from_terms`, e.g. while building a generator whose width is not yet
-                known); [get_majorana_operator][] then raises.
+                act within ``0..num_qubits-1``.
 
         Raises:
+            TypeError: If ``num_qubits`` is not an integer.
             ValueError: If a term acts on a qubit index ``>= num_qubits``.
         """
         accumulated: dict[Pauli, complex] = defaultdict(complex)
@@ -137,21 +137,20 @@ def __init__(
             pauli = key if isinstance(key, Pauli) else Pauli(key)
             accumulated[pauli] += coeff
         self.terms: dict[Pauli, complex] = dict(accumulated)
-        self.num_qubits = num_qubits
-        if num_qubits is not None:
-            for pauli in self.terms:
-                if pauli.qubits and pauli.qubits[-1] >= num_qubits:
-                    raise ValueError(
-                        f"Pauli term {pauli} acts on a qubit index >= num_qubits="
-                        f"{num_qubits}."
-                    )
+        self.num_qubits = _validate_system_size(num_qubits, argument_name="num_qubits")
+        for pauli in self.terms:
+            if pauli.qubits and pauli.qubits[-1] >= self.num_qubits:
+                raise ValueError(
+                    f"Pauli term {pauli} acts on a qubit index >= num_qubits="
+                    f"{self.num_qubits}."
+                )
 
     @classmethod
     def _from_terms(
         cls,
         strings: Sequence[Pauli | str],
         coefficients: Sequence[complex],
-        num_qubits: int | None = None,
+        num_qubits: int,
     ) -> PauliOperator:
         """Build from parallel ``strings``/``coefficients`` lists (internal)."""
         accumulated: dict[Pauli, complex] = defaultdict(complex)
@@ -217,15 +216,7 @@ def get_majorana_operator(self) -> MajoranaOperator:
 
         Each local term is extended to the full ``num_qubits`` width (identities filled in)
         before the Jordan-Wigner map, so the resulting Majorana indices are global.
-
-        Raises:
-            ValueError: If ``num_qubits`` is unset.
         """
-        if self.num_qubits is None:
-            raise ValueError(
-                "PauliOperator.get_majorana_operator() needs num_qubits; construct the "
-                "operator with an explicit num_qubits."
-            )
         majoranas: list[Sequence[int]] = []
         coefficients: list[complex] = []
         for pauli, coeff in self.terms.items():
diff --git a/src/monoprop/qiskit_conversion.py b/src/monoprop/qiskit_conversion.py
index 14501dba..5156c171 100644
--- a/src/monoprop/qiskit_conversion.py
+++ b/src/monoprop/qiskit_conversion.py
@@ -31,6 +31,7 @@
 
 from monoprop.circuit import Circuit, ExpGate
 from monoprop.pauli import Pauli, PauliOperator
+from monoprop.utils import _validate_system_size
 
 PAULI_EVOLUTION_EQUIVALENT = {
     "rx",
@@ -168,6 +169,7 @@ def from_qiskit_circuit(
         gates=tuple(gates),
         parameters=tuple(parameters),
         initial_state=tuple(initial_state),
+        num_modes=num_qubits,
     )
 
 
@@ -197,12 +199,16 @@ def to_qiskit_circuit(circuit: Circuit, num_qubits: int) -> QuantumCircuit:
 
     Args:
         circuit: A [Circuit][monoprop.circuit.Circuit] representing the given circuit.
-        num_qubits: Total number of qubits (the circuit no longer carries it; supply the
-            observable's ``num_qubits``).
+        num_qubits: Total number of qubits. Must match ``circuit.num_modes``.
 
     Returns:
         A qiskit quantum circuit.
     """
+    num_qubits = _validate_system_size(num_qubits, argument_name="num_qubits")
+    if num_qubits != circuit.num_modes:
+        raise ValueError(
+            f"num_qubits={num_qubits} does not match circuit.num_modes={circuit.num_modes}."
+        )
     if len(circuit.parameters) != circuit.n_parameters:
         raise ValueError(
             f"to_qiskit_circuit needs a bound circuit: it has {circuit.n_parameters} "
diff --git a/src/monoprop/utils.py b/src/monoprop/utils.py
index e5a8c1c5..820e1b5a 100644
--- a/src/monoprop/utils.py
+++ b/src/monoprop/utils.py
@@ -17,6 +17,29 @@
 from __future__ import annotations
 
 
+def _validate_system_size(size: int, *, argument_name: str) -> int:
+    """Validate and normalize a positive system-size argument.
+
+    Args:
+        size: Number of modes/qubits.
+        argument_name: Public argument name for error messages.
+
+    Returns:
+        The normalized integer value.
+
+    Raises:
+        TypeError: If ``size`` is not an integer (or is ``bool``).
+        ValueError: If ``size`` is not positive.
+    """
+    if isinstance(size, bool) or not isinstance(size, int):
+        raise TypeError(
+            f"{argument_name} must be an integer (not {type(size).__name__})."
+        )
+    if size <= 0:
+        raise ValueError(f"{argument_name} must be positive; got {size}.")
+    return int(size)
+
+
 def jordan_wigner_basis_change(n_qubits: int) -> list[list[int]]:
     """Generate a basis change for Jordan-Wigner representation.
 
diff --git a/tests/cases.py b/tests/cases.py
index 5dd8f1a9..b4406431 100644
--- a/tests/cases.py
+++ b/tests/cases.py
@@ -44,6 +44,7 @@ class DenseMajoranaArrays:
     parameters: list[float] | ndarray
     gen_coeffs: list[float] | ndarray
     param_inds: list[int] | ndarray
+    num_modes: int
 
     def to_circuit(self) -> Circuit:
         """Group the dense arrays into a :class:`~monoprop.circuit.Circuit`."""
@@ -51,6 +52,7 @@ def to_circuit(self) -> Circuit:
             majoranas=self.majoranas,
             gen_coeffs=self.gen_coeffs,
             param_inds=self.param_inds,
+            num_modes=self.num_modes,
             parameters=self.parameters,
             initial_state=self.initial_state,
         )
@@ -95,6 +97,7 @@ def load_problem(path: Path) -> FermionicProblem:
         gen_coeffs=np.asarray(data["gen_coeffs"]),
         param_inds=np.asarray(data["param_inds"], dtype=int),
         parameters=np.asarray(data["parameters"]),
+        num_modes=int(data["num_modes"]),
     )
 
     ham = data["hamiltonian"]
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index 59eed9c6..839d7c60 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -320,7 +320,9 @@ def test_hermitian_majorana_generator_matches_structural() -> None:
     hermitian = Circuit(
         gates=(ExpGate(MajoranaOperator({(4, 5): 1j}, num_modes=8)),), parameters=(0.5,)
     )
-    structural = Circuit.from_dense_arrays([[4, 5]], [-1.0], [0], parameters=[0.5])
+    structural = Circuit.from_dense_arrays(
+        [[4, 5]], [-1.0], [0], num_modes=8, parameters=[0.5]
+    )
 
     from_hermitian = MajoranaPropagator.from_circuit(
         hermitian, obs, cutoff=16
diff --git a/tests/test_coeff_trunc.py b/tests/test_coeff_trunc.py
index b5b87b52..3cf50378 100644
--- a/tests/test_coeff_trunc.py
+++ b/tests/test_coeff_trunc.py
@@ -62,6 +62,7 @@ def test_coeff_trunc(serial_comm):
         parameters=[np.pi / 6],
         gen_coeffs=[1.0],
         param_inds=[0],
+        num_modes=n_modes,
     )
     circuit = sequence
 
@@ -146,6 +147,7 @@ def test_evolution_coeff_trunc_no_atols(serial_comm):
         parameters=[p],
         gen_coeffs=[1.0],
         param_inds=[0],
+        num_modes=n_modes,
     )
     circuit = sequence
 
@@ -183,6 +185,7 @@ def test_evolution_coeff_trunc_small_coeffs(serial_comm):
         parameters=[p],
         gen_coeffs=[1.0],
         param_inds=[0],
+        num_modes=n_modes,
     )
     circuit = sequence
 
diff --git a/tests/test_majorana.py b/tests/test_majorana.py
index 4eed640d..f3853761 100644
--- a/tests/test_majorana.py
+++ b/tests/test_majorana.py
@@ -64,6 +64,7 @@ def test_from_dense_arrays_groups_by_param_ind():
         majoranas=[(0, 1), (2, 3), (0, 3)],
         gen_coeffs=[0.5, -0.5, 1.0],
         param_inds=[0, 0, 1],
+        num_modes=2,
         parameters=[1.0, 2.0],
         initial_state=[0, 1],
     )
diff --git a/tests/test_nonfermi.py b/tests/test_nonfermi.py
index 098d51db..c534d940 100644
--- a/tests/test_nonfermi.py
+++ b/tests/test_nonfermi.py
@@ -33,6 +33,7 @@ def test_nonfermi(serial_comm):
         majoranas=majoranas,
         gen_coeffs=gen_coeffs,
         param_inds=param_inds,
+        num_modes=num_modes,
         parameters=parameters,
         initial_state=[],
     )
diff --git a/tests/test_only_rotate_k.py b/tests/test_only_rotate_k.py
index 1f34e85e..020344d4 100644
--- a/tests/test_only_rotate_k.py
+++ b/tests/test_only_rotate_k.py
@@ -55,6 +55,7 @@ def test_basic_orbital_rotation(serial_comm):
         parameters=[np.pi / 4],
         gen_coeffs=[1.0],
         param_inds=[0],
+        num_modes=n_modes,
     )
     circuit = sequence
     kwargs = {"cutoff": 6, "schrodinger_cutoff": 8, "comm": serial_comm}
diff --git a/tests/test_update_methods.py b/tests/test_update_methods.py
index 90a88cea..81b4cbf1 100644
--- a/tests/test_update_methods.py
+++ b/tests/test_update_methods.py
@@ -144,6 +144,7 @@ def test_integration(self, serial_comm):
             majoranas=[(0, 2), (1, 3)],
             gen_coeffs=[0.0, 0.0],
             param_inds=[0, 1],
+            num_modes=4,
             parameters=[1.0, 1.0],
         )
         mp = MajoranaPropagator(

From afd910682e95cb4421561057ddad15b3bed29153 Mon Sep 17 00:00:00 2001
From: Adam <aglos@iitis.pl>
Date: Mon, 27 Jul 2026 20:58:53 +0000
Subject: [PATCH 2/6] test: :white_check_mark: tests fixed

---
 src/monoprop/circuit.py             |   4 +-
 src/monoprop/integral_conversion.py |   3 +-
 src/monoprop/utils.py               |   8 +-
 tests/test_circuit.py               | 165 ++++++++++++++++++++++------
 tests/test_fermi.py                 |  38 ++++---
 tests/test_integral_conversion.py   |   5 +-
 tests/test_monoprop_trivial.py      |  10 +-
 tests/test_only_rotate_k.py         |  14 ++-
 tests/test_parameter_validation.py  |  17 ++-
 tests/test_pauli.py                 |  25 +++--
 tests/test_qiskit_conversion.py     |  24 ++--
 tests/test_update_methods.py        |   2 +-
 12 files changed, 228 insertions(+), 87 deletions(-)

diff --git a/src/monoprop/circuit.py b/src/monoprop/circuit.py
index 4b6e92cc..e026fcdc 100644
--- a/src/monoprop/circuit.py
+++ b/src/monoprop/circuit.py
@@ -264,16 +264,16 @@ class Circuit:
     def __init__(  # noqa:C901
         self,
         gates: Sequence[ExpGate],
-        initial_state: Sequence[int],
         num_modes: int,
+        initial_state: Sequence[int] = (),
         parameters: Sequence[float] = (),
     ) -> None:
         """Build the circuit, dropping identity gates and validating family/mapping/params.
 
         Args:
             gates: The ordered exponential gates.
-            num_modes: Number of modes/qubits for the circuit.
             initial_state: The reference state (occupied mode / qubit indices).
+            num_modes: Number of modes/qubits for the circuit.
             parameters: The angle values, or empty for an unbound circuit.
 
         Raises:
diff --git a/src/monoprop/integral_conversion.py b/src/monoprop/integral_conversion.py
index 906ce583..5261e31e 100644
--- a/src/monoprop/integral_conversion.py
+++ b/src/monoprop/integral_conversion.py
@@ -104,4 +104,5 @@ def integrals_to_fermion(
             continue
         terms[ind] += coeff
 
-    return FermiOperator.from_dict(terms)
+    # 2* because we have alpha and beta orbitals
+    return FermiOperator.from_dict(terms, 2 * hamiltonian[1].shape[1])
diff --git a/src/monoprop/utils.py b/src/monoprop/utils.py
index 820e1b5a..0db69fb9 100644
--- a/src/monoprop/utils.py
+++ b/src/monoprop/utils.py
@@ -18,7 +18,7 @@
 
 
 def _validate_system_size(size: int, *, argument_name: str) -> int:
-    """Validate and normalize a positive system-size argument.
+    """Validate and normalize a non-negative system-size argument.
 
     Args:
         size: Number of modes/qubits.
@@ -29,14 +29,14 @@ def _validate_system_size(size: int, *, argument_name: str) -> int:
 
     Raises:
         TypeError: If ``size`` is not an integer (or is ``bool``).
-        ValueError: If ``size`` is not positive.
+        ValueError: If ``size`` is negative.
     """
     if isinstance(size, bool) or not isinstance(size, int):
         raise TypeError(
             f"{argument_name} must be an integer (not {type(size).__name__})."
         )
-    if size <= 0:
-        raise ValueError(f"{argument_name} must be positive; got {size}.")
+    if size < 0:
+        raise ValueError(f"{argument_name} must be non-negative; got {size}.")
     return int(size)
 
 
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index 839d7c60..9ffefc56 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -162,17 +162,19 @@ def test_exp_gate_applies_atol_truncation(
 def test_circuit_equality() -> None:
     """Circuits are equal on gates/parameters/initial_state; family is derived, not compared."""
     gen = MajoranaOperator({(0, 1): 1.0j}, num_modes=2)
-    a = Circuit((ExpGate(gen),), parameters=(0.3,), initial_state=(0,))
-    b = Circuit((ExpGate(gen),), parameters=(0.3,), initial_state=(0,))
+    a = Circuit((ExpGate(gen),), initial_state=(0,), num_modes=2, parameters=(0.3,))
+    b = Circuit((ExpGate(gen),), initial_state=(0,), num_modes=2, parameters=(0.3,))
     assert a == b
-    assert a != Circuit((ExpGate(gen),), parameters=(0.9,), initial_state=(0,))
+    assert a != Circuit(
+        (ExpGate(gen),), initial_state=(0,), num_modes=2, parameters=(0.9,)
+    )
     assert a != "not a circuit"
 
 
 def test_circuit_rejects_non_exp_gate() -> None:
     """A gate that is not an ExpGate is rejected with a clear TypeError."""
     with pytest.raises(TypeError, match="Circuit gates must be ExpGate"):
-        Circuit(("not a gate",))  # type: ignore[arg-type]
+        Circuit(("not a gate",), initial_state=(), num_modes=0)  # type: ignore[arg-type]
 
 
 def test_to_circuit_round_trips_sequence() -> None:
@@ -197,7 +199,9 @@ def test_default_mapping_is_identity() -> None:
         (
             ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),
             ExpGate(MajoranaOperator({(2, 3): 1.0j}, num_modes=2)),
-        )
+        ),
+        initial_state=(),
+        num_modes=2,
     )
     assert circuit.resolved_mapping == (0, 1)
     assert circuit.n_parameters == 2
@@ -212,7 +216,7 @@ def test_shared_mapping_index_ties_gates() -> None:
             MajoranaOperator({(0, 3): 1.0j}, num_modes=2), index=0
         ),  # ties to the first
     )
-    circuit = Circuit(gates)
+    circuit = Circuit(gates, initial_state=(), num_modes=2)
     assert circuit.resolved_mapping == (0, 1, 0)
     assert circuit.n_parameters == 2
 
@@ -224,7 +228,7 @@ def test_circuit_rejects_non_contiguous_mapping() -> None:
         ExpGate(MajoranaOperator({(2, 3): 1.0j}, num_modes=2), index=2),
     )
     with pytest.raises(ValueError, match="contiguous"):
-        Circuit(gates)
+        Circuit(gates, initial_state=(), num_modes=2)
 
 
 def test_circuit_rejects_mixed_param_scheme() -> None:
@@ -234,14 +238,14 @@ def test_circuit_rejects_mixed_param_scheme() -> None:
         ExpGate(MajoranaOperator({(2, 3): 1.0j}, num_modes=2)),
     )
     with pytest.raises(ValueError, match="every gate must set"):
-        Circuit(gates)
+        Circuit(gates, initial_state=(), num_modes=2)
 
 
 def test_circuit_rejects_wrong_parameter_length() -> None:
     """A bound circuit must supply exactly one value per distinct angle."""
     gates = (ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),)
     with pytest.raises(ValueError, match="1 parameters"):
-        Circuit(gates, parameters=(0.1, 0.2))
+        Circuit(gates, initial_state=(), num_modes=2, parameters=(0.1, 0.2))
 
 
 def test_circuit_add_offsets_second_axis() -> None:
@@ -251,10 +255,15 @@ def test_circuit_add_offsets_second_axis() -> None:
             ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),
             ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),
         ),
+        initial_state=(),
+        num_modes=2,
         parameters=(0.1, 0.2),
     )
     b = Circuit(
-        (ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2)),), parameters=(0.3,)
+        (ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2)),),
+        initial_state=(),
+        num_modes=2,
+        parameters=(0.3,),
     )
     combined = a + b
     assert combined.resolved_mapping == (0, 1, 2)
@@ -264,8 +273,16 @@ def test_circuit_add_offsets_second_axis() -> None:
 
 def test_circuit_add_rejects_mixed_families() -> None:
     """Concatenating a Majorana circuit with a qubit circuit raises a clear TypeError."""
-    maj = Circuit((ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),))
-    qubit = Circuit((ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=2)),))
+    maj = Circuit(
+        (ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),),
+        initial_state=(),
+        num_modes=2,
+    )
+    qubit = Circuit(
+        (ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=2)),),
+        initial_state=(),
+        num_modes=2,
+    )
     with pytest.raises(TypeError, match="gate families differ"):
         _ = maj + qubit
 
@@ -273,10 +290,14 @@ def test_circuit_add_rejects_mixed_families() -> None:
 def test_circuit_add_rejects_different_initial_states() -> None:
     """Concatenating circuits with conflicting initial states raises a ValueError."""
     a = Circuit(
-        (ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),), initial_state=(0,)
+        (ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),),
+        initial_state=(0,),
+        num_modes=2,
     )
     b = Circuit(
-        (ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),), initial_state=(1,)
+        (ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),),
+        initial_state=(1,),
+        num_modes=2,
     )
     with pytest.raises(ValueError, match="different initial states"):
         _ = a + b
@@ -289,7 +310,12 @@ def test_bound_circuit_with_identity_gate_wrong_param_count_raises() -> None:
         ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),
     )
     with pytest.raises(ValueError, match="2 gates"):
-        Circuit(gates, parameters=(0.1,))  # 1 value, but 2 gates before the drop
+        Circuit(
+            gates,
+            initial_state=(),
+            num_modes=2,
+            parameters=(0.1,),
+        )  # 1 value, but 2 gates before the drop
 
 
 def test_non_hermitian_majorana_generator_rejected() -> None:
@@ -302,7 +328,10 @@ def test_non_hermitian_majorana_generator_rejected() -> None:
     obs = MajoranaOperator({(0, 1): 1.0j}, num_modes=2)
     prop = MajoranaPropagator(obs, [0, 1], cutoff=4)
     bad = Circuit(
-        (ExpGate(MajoranaOperator({(0, 1): 1.0}, num_modes=2)),), parameters=(0.3,)
+        (ExpGate(MajoranaOperator({(0, 1): 1.0}, num_modes=2)),),
+        initial_state=(),
+        num_modes=2,
+        parameters=(0.3,),
     )
     with pytest.raises(ValueError, match="not Hermitian"):
         prop.propagate(bad)
@@ -318,7 +347,10 @@ def test_hermitian_majorana_generator_matches_structural() -> None:
     """
     obs = MajoranaOperator({(0, 1, 2, 4): 1.0}, 8)
     hermitian = Circuit(
-        gates=(ExpGate(MajoranaOperator({(4, 5): 1j}, num_modes=8)),), parameters=(0.5,)
+        gates=(ExpGate(MajoranaOperator({(4, 5): 1j}, num_modes=8)),),
+        initial_state=(),
+        num_modes=8,
+        parameters=(0.5,),
     )
     structural = Circuit.from_dense_arrays(
         [[4, 5]], [-1.0], [0], num_modes=8, parameters=[0.5]
@@ -458,7 +490,7 @@ def _multi_term_gate_propagator():
     # two monomials -> two layers
     g0 = ExpGate(MajoranaOperator({(0, 2): 1.0j, (1, 3): 1.0j}, num_modes=2))
     g1 = ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2))
-    prop.build_graph(Circuit((g0, g1)))
+    prop.build_graph(Circuit((g0, g1), initial_state=(), num_modes=2))
     return prop
 
 
@@ -473,9 +505,21 @@ def test_n_gates_accumulates_across_builds() -> None:
     """Each build_graph call extends the gate count on the accumulated axis."""
     op = MajoranaOperator({(0, 1): 1.0j}, num_modes=2)
     prop = MajoranaPropagator(op, [0, 1], cutoff=4)
-    prop.build_graph(Circuit((ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),)))
+    prop.build_graph(
+        Circuit(
+            (ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),),
+            initial_state=(),
+            num_modes=2,
+        )
+    )
     assert prop.n_gates == 1
-    prop.build_graph(Circuit((ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),)))
+    prop.build_graph(
+        Circuit(
+            (ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),),
+            initial_state=(),
+            num_modes=2,
+        )
+    )
     assert prop.n_gates == 2
 
 
@@ -499,7 +543,9 @@ def test_majorana_propagator_rejects_pauli_circuit() -> None:
     problem = load_problem(DATA / "rx_rz_ry_rz_exact.msgpack")
     prop = _propagator(problem)
     circuit = Circuit(
-        gates=(ExpGate(PauliOperator({"Z": 1.0}, num_qubits=1)),),
+        gates=(ExpGate(PauliOperator({"Z": 1.0}, num_qubits=problem.n_modes)),),
+        initial_state=(),
+        num_modes=problem.n_modes,
     )
 
     with pytest.raises(TypeError, match="qubit"):
@@ -510,8 +556,13 @@ def test_propagate_rejects_mismatched_initial_state() -> None:
     """A circuit whose initial_state disagrees with the propagator's raises ValueError."""
     problem = load_problem(DATA / "rx_rz_ry_rz_exact.msgpack")
     prop = _propagator(problem)  # built with the fixture's initial state ([])
-    gate = ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2))
-    circuit = Circuit((gate,), parameters=(0.1,), initial_state=(0, 1))
+    gate = ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=problem.n_modes))
+    circuit = Circuit(
+        (gate,),
+        initial_state=(0, 1),
+        num_modes=problem.n_modes,
+        parameters=(0.1,),
+    )
 
     with pytest.raises(ValueError, match="initial_state"):
         prop.propagate(circuit)
@@ -521,8 +572,10 @@ def test_propagate_accepts_empty_initial_state() -> None:
     """An empty circuit.initial_state defers to the propagator's reference state."""
     problem = load_problem(DATA / "rx_rz_ry_rz_exact.msgpack")
     prop = _propagator(problem)
-    gate = ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2))
-    circuit = Circuit((gate,), parameters=(0.1,))  # empty initial_state
+    gate = ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=problem.n_modes))
+    circuit = Circuit(
+        (gate,), initial_state=(), num_modes=problem.n_modes, parameters=(0.1,)
+    )  # empty initial_state
 
     prop.propagate(circuit)  # does not raise
 
@@ -551,9 +604,13 @@ def test_build_graph_accumulates_layers_and_parameters(fixture: str) -> None:
     single.build_graph(circuit)
 
     twice = _propagator(problem)
-    twice.build_graph(Circuit(_rebase(gates[:split])))
+    twice.build_graph(
+        Circuit(_rebase(gates[:split]), initial_state=(), num_modes=problem.n_modes)
+    )
     layers_after_first = twice.graph_layers
-    twice.build_graph(Circuit(_rebase(gates[split:])))
+    twice.build_graph(
+        Circuit(_rebase(gates[split:]), initial_state=(), num_modes=problem.n_modes)
+    )
 
     assert 0 < layers_after_first < twice.graph_layers
     assert twice.graph_layers == single.graph_layers
@@ -573,8 +630,18 @@ def test_compose_then_single_build_matches_single_call(fixture: str) -> None:
     gates = circuit.gates
     split = len(gates) // 2
     params = list(map(float, problem.monomial_circuit.parameters))
-    a = Circuit(_rebase(gates[:split]), parameters=tuple(params[:split]))
-    b = Circuit(_rebase(gates[split:]), parameters=tuple(params[split:]))
+    a = Circuit(
+        _rebase(gates[:split]),
+        initial_state=(),
+        num_modes=problem.n_modes,
+        parameters=tuple(params[:split]),
+    )
+    b = Circuit(
+        _rebase(gates[split:]),
+        initial_state=(),
+        num_modes=problem.n_modes,
+        parameters=tuple(params[split:]),
+    )
     composed = a + b
 
     single = _propagator(problem)
@@ -606,8 +673,12 @@ def test_build_graph_in_two_calls_schrodinger(fixture: str) -> None:
     single.build_graph(circuit)
 
     twice = _schrodinger_propagator(problem)
-    twice.build_graph(Circuit(_rebase(gates[:split])))
-    twice.build_graph(Circuit(_rebase(gates[split:])))
+    twice.build_graph(
+        Circuit(_rebase(gates[:split]), initial_state=(), num_modes=problem.n_modes)
+    )
+    twice.build_graph(
+        Circuit(_rebase(gates[split:]), initial_state=(), num_modes=problem.n_modes)
+    )
 
     np.testing.assert_allclose(
         twice.expectation_value(params), single.expectation_value(params)
@@ -625,10 +696,15 @@ def test_build_graph_twice_with_seed_regeneration(fixture: str) -> None:
     split = len(gates) // 2
 
     prop = _schrodinger_propagator(problem)
-    prop.build_graph(Circuit(_rebase(gates[:split])))
+    prop.build_graph(
+        Circuit(_rebase(gates[:split]), initial_state=(), num_modes=problem.n_modes)
+    )
     # seed_parameters on the second call exercises the internal seed regeneration
     # (the former operator_coeffs round-trip) used for coefficient-informed truncation.
-    prop.build_graph(Circuit(_rebase(gates[split:])), seed_parameters=params)
+    prop.build_graph(
+        Circuit(_rebase(gates[split:]), initial_state=(), num_modes=problem.n_modes),
+        seed_parameters=params,
+    )
 
     np.testing.assert_allclose(prop.expectation_value(params), problem.exact_expval)
 
@@ -652,6 +728,8 @@ def test_empty_default_mapping_gate_dropped_and_evaluable() -> None:
                 MajoranaOperator({(6, 7): 0.0}, num_modes=8)
             ),  # identity generator: dropped
         ),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.5, 0.3),
     )
     assert len(circuit.gates) == 1
@@ -674,6 +752,8 @@ def test_empty_gate_in_middle_builds_contiguously() -> None:
             ),  # identity in the middle
             ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),
         ),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.5, 0.3, 0.2),
     )
     assert len(circuit.gates) == 2
@@ -688,7 +768,12 @@ def test_surplus_parameters_raise_not_truncated() -> None:
         [[(0, "+"), (1, "-")], [(1, "+"), (0, "-")]], [1.0, 1.0], num_modes=4
     )
     with pytest.raises(ValueError, match="1 parameter"):
-        Circuit(gates=(ExpGate(generator),), parameters=(1.0, 2.0))
+        Circuit(
+            gates=(ExpGate(generator),),
+            initial_state=(),
+            num_modes=4,
+            parameters=(1.0, 2.0),
+        )
 
 
 def test_non_commuting_pauli_generator_rejected() -> None:
@@ -742,6 +827,8 @@ def test_build_graph_seed_parameters_accepts_numpy() -> None:
             ExpGate(MajoranaOperator({(4, 5): -1.0j}, num_modes=8)),
             ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),
         ),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.5, 0.3),
     )
     prop = _small_propagator(lower_atol=1e-12)
@@ -755,6 +842,8 @@ def test_build_graph_rejects_too_short_seed() -> None:
             ExpGate(MajoranaOperator({(4, 5): -1.0j}, num_modes=8)),
             ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),
         ),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.5, 0.3),
     )
     prop = _small_propagator()
@@ -768,10 +857,14 @@ def test_extend_without_seed_builds_structurally() -> None:
     full-axis seed is still accepted."""
     c1 = Circuit(
         gates=(ExpGate(MajoranaOperator({(4, 5): -1.0j}, num_modes=8)),),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.3,),
     )
     c2 = Circuit(
         gates=(ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.4,),
     )
     params = [0.3, 0.4]
@@ -796,10 +889,14 @@ def test_propagate_after_build_graph_rejected() -> None:
     """propagate() on top of a build_graph() graph raises rather than corrupting it."""
     c1 = Circuit(
         gates=(ExpGate(MajoranaOperator({(4, 5): -1.0j}, num_modes=8)),),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.3,),
     )
     c2 = Circuit(
         gates=(ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),),
+        initial_state=(),
+        num_modes=8,
         parameters=(0.4,),
     )
     prop = _small_propagator()
diff --git a/tests/test_fermi.py b/tests/test_fermi.py
index a86c3cec..e6568dd7 100644
--- a/tests/test_fermi.py
+++ b/tests/test_fermi.py
@@ -64,13 +64,13 @@ def test_repr_empty(self):
 class TestFermiOperator:
     def test_valid_creation(self):
         terms = [FermiString([(0, "+")]), FermiString([(1, "-")])]
-        op = FermiOperator(terms, [1.0, -1.0])
+        op = FermiOperator(terms, [1.0, -1.0], num_modes=2)
         assert op.terms == terms
         assert op.coefficients == [1.0, -1.0]
 
     def test_num_modes_max_index(self):
         terms = [FermiString([(0, "+"), (3, "-")]), FermiString([(5, "+")])]
-        op = FermiOperator(terms, [1.0, 1.0])
+        op = FermiOperator(terms, [1.0, 1.0], num_modes=6)
         assert op.num_modes == 6
 
     def test_num_modes_empty_terms(self):
@@ -78,24 +78,24 @@ def test_num_modes_empty_terms(self):
         assert op.num_modes == 10
 
     def test_num_modes_single_term(self):
-        op = FermiOperator([FermiString([(2, "+"), (7, "-")])], [1.0])
+        op = FermiOperator([FermiString([(2, "+"), (7, "-")])], [1.0], num_modes=8)
         assert op.num_modes == 8
 
     def test_terms_is_copy(self):
         terms = [FermiString([(0, "+")])]
-        op = FermiOperator(terms, [1.0])
+        op = FermiOperator(terms, [1.0], num_modes=1)
         terms.append(FermiString([(1, "-")]))
         assert len(op.terms) == 1
 
     def test_coefficients_is_copy(self):
         coeffs = [1.0]
-        op = FermiOperator([FermiString([(0, "+")])], coeffs)
+        op = FermiOperator([FermiString([(0, "+")])], coeffs, num_modes=1)
         coeffs.append(2.0)
         assert len(op.coefficients) == 1
 
     def test_str(self):
         terms = [FermiString([(0, "+")])]
-        op = FermiOperator(terms, [2.0])
+        op = FermiOperator(terms, [2.0], num_modes=1)
         assert str(op) == "FermiOperator(1 terms, 1 modes: 2.0*FermiString(c_0^+))"
 
     def test_num_modes_explicit_override(self):
@@ -104,11 +104,8 @@ def test_num_modes_explicit_override(self):
         assert op.num_modes == 12
 
     def test_empty_terms_without_num_modes_raises(self):
-        with pytest.raises(
-            ValueError,
-            match=r"max\(\) (arg is an empty sequence|iterable argument is empty)",
-        ):
-            FermiOperator([], [])
+        with pytest.raises(TypeError):
+            FermiOperator([], [])  # type: ignore[call-arg]
 
     @pytest.mark.parametrize(
         ("left", "right", "expected"),
@@ -323,7 +320,12 @@ class TestCircuit:
     def test_len(self):
         gates = [ExpGate(_number_op()), ExpGate(_number_op())]
 
-        circuit = Circuit(gates=gates, parameters=[0.1, 0.2], initial_state=[0])
+        circuit = Circuit(
+            gates=gates,
+            initial_state=[0],
+            num_modes=1,
+            parameters=[0.1, 0.2],
+        )
 
         assert len(circuit) == 2
 
@@ -336,7 +338,10 @@ def test_converts_fermi_gates_to_majorana(self):
         gate_1 = ExpGate(hop)
 
         circuit = Circuit(
-            gates=[gate_0, gate_1], parameters=[0.3, -0.7], initial_state=[0, 1]
+            gates=[gate_0, gate_1],
+            initial_state=[0, 1],
+            num_modes=2,
+            parameters=[0.3, -0.7],
         )
 
         np.testing.assert_array_equal(circuit.initial_state, np.array([0, 1]))
@@ -347,7 +352,7 @@ def test_converts_fermi_gates_to_majorana(self):
         assert all(g.family == "majorana" for g in circuit.gates)
 
         majoranas, gen_coeffs, per_monomial_mapping, gate_indices = expand_monomials(
-            circuit.gates, circuit.resolved_mapping
+            circuit.gates, circuit.resolved_mapping, circuit.num_modes
         )
         # One gate per fermi gate; each generator here has two monomials.
         n_terms = len(gate_0.generator.terms)
@@ -364,12 +369,13 @@ def test_drops_identity_generators_and_aligned_parameters(self):
         identity = FermiOperator([], [], num_modes=1)
         circuit = Circuit(
             gates=[ExpGate(_number_op()), ExpGate(identity), ExpGate(_number_op())],
-            parameters=[0.1, 0.2, 0.3],
             initial_state=[0],
+            num_modes=1,
+            parameters=[0.1, 0.2, 0.3],
         )
         assert len(circuit) == 2  # the identity gate is dropped
         assert circuit.parameters == (0.1, 0.3)  # its aligned parameter goes with it
 
     def test_validate_inputs_duplicate_initial_state_raises(self):
         with pytest.raises(ValueError, match="Duplicate indices in initial state"):
-            Circuit(gates=[ExpGate(_number_op())], initial_state=[0, 0])
+            Circuit(gates=[ExpGate(_number_op())], initial_state=[0, 0], num_modes=1)
diff --git a/tests/test_integral_conversion.py b/tests/test_integral_conversion.py
index 9d0210b9..72f7f902 100644
--- a/tests/test_integral_conversion.py
+++ b/tests/test_integral_conversion.py
@@ -28,12 +28,15 @@ def _read_openfermion(path: Path) -> FermiOperator:
         terms = json.load(f)
     coeffs = []
     fermi_strings = []
+    num_modes = 0
     for key, value in terms.items():
         fixed_key = ast.literal_eval(key) if key != "()" else ()
         coeffs.append(value)
         fermi_terms = [(el[0], "-" if el[1] == 0 else "+") for el in fixed_key]
         fermi_strings.append(FermiString(fermi_terms))
-    return FermiOperator(fermi_strings, coeffs)
+        if fermi_terms:
+            num_modes = max(num_modes, max(el[0] for el in fermi_terms) + 1)
+    return FermiOperator(fermi_strings, coeffs, num_modes)
 
 
 @pytest.fixture
diff --git a/tests/test_monoprop_trivial.py b/tests/test_monoprop_trivial.py
index f9fc659a..634a1800 100644
--- a/tests/test_monoprop_trivial.py
+++ b/tests/test_monoprop_trivial.py
@@ -39,7 +39,9 @@ def test_trivial_evolved_operator_cases(
 ):
     """Test trivial evolved operator dict for various initial conditions."""
     kwargs = {"schrodinger_cutoff": schrodinger_cutoff} if schrodinger_cutoff else {}
-    quantum_circuit = Circuit(initial_state=[], gates=[])
+    quantum_circuit = Circuit(
+        initial_state=[], num_modes=initial_op.num_modes, gates=[]
+    )
     mp = MajoranaPropagator(
         initial_op,
         quantum_circuit.initial_state,
@@ -53,7 +55,9 @@ def test_trivial_evolved_operator_cases(
 
 def test_trivial_evolved_operator(serial_comm):
     initial_op = MajoranaOperator({(0, 1, 2, 4): 1}, 8)
-    quantum_circuit = Circuit(initial_state=[], gates=[])
+    quantum_circuit = Circuit(
+        initial_state=[], num_modes=initial_op.num_modes, gates=[]
+    )
     mp = MajoranaPropagator(
         initial_op, quantum_circuit.initial_state, cutoff=16, comm=serial_comm
     )
@@ -102,7 +106,7 @@ def test_update_initial_operator(
 ):
     """Test updating coefficients in both regular and Schrodinger pictures."""
     kwargs = {"schrodinger_cutoff": schrodinger_cutoff} if schrodinger_cutoff else {}
-    quantum_circuit = Circuit(initial_state=[], gates=[])
+    quantum_circuit = Circuit(initial_state=[], num_modes=init_op.num_modes, gates=[])
     mp = MajoranaPropagator(
         init_op,
         quantum_circuit.initial_state,
diff --git a/tests/test_only_rotate_k.py b/tests/test_only_rotate_k.py
index 020344d4..06ce66fc 100644
--- a/tests/test_only_rotate_k.py
+++ b/tests/test_only_rotate_k.py
@@ -96,10 +96,14 @@ def test_only_rotate_len_k(problem, inplace, serial_mp_kwargs):
     # circuit, so a plain ExpGate(gate.generator) would re-normalize their structural coefficients).
     non_orbital = Circuit(
         tuple(ExpGate._with_index(gate, None) for gate in non_orbital_gates),
+        initial_state=(),
+        num_modes=problem.n_modes,
         parameters=tuple(parameters[:split]),
     )
     orbital = Circuit(
         tuple(ExpGate._with_index(gate, None) for gate in orbital_gates),
+        initial_state=(),
+        num_modes=problem.n_modes,
         parameters=tuple(parameters[split:]),
     )
 
@@ -153,7 +157,10 @@ def test_only_rotate_len_k_errors_majorana(only_rotate_len_k, err, method_name):
     """Test that invalid only_rotate_len_k raises ValueError."""
     mp = MajoranaPropagator(MajoranaOperator({}, 4), [], cutoff=6, schrodinger_cutoff=8)
     with err:
-        getattr(mp, method_name)(Circuit(), only_rotate_len_k=only_rotate_len_k)
+        getattr(mp, method_name)(
+            Circuit((), initial_state=(), num_modes=4),
+            only_rotate_len_k=only_rotate_len_k,
+        )
 
 
 @pytest.mark.parametrize(
@@ -188,4 +195,7 @@ def test_only_rotate_len_k_errors_pauli(only_rotate_len_k, err, method_name):
     """Test that invalid only_rotate_len_k raises ValueError."""
     mp = PauliPropagator(PauliOperator({}, 4), [], cutoff=6, schrodinger_cutoff=8)
     with err:
-        getattr(mp, method_name)(Circuit(), only_rotate_len_k=only_rotate_len_k)
+        getattr(mp, method_name)(
+            Circuit((), initial_state=(), num_modes=4),
+            only_rotate_len_k=only_rotate_len_k,
+        )
diff --git a/tests/test_parameter_validation.py b/tests/test_parameter_validation.py
index 182cef90..c638f537 100644
--- a/tests/test_parameter_validation.py
+++ b/tests/test_parameter_validation.py
@@ -30,7 +30,8 @@ def _two_gate_graph(serial_comm):
         (
             ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),
             ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),
-        )
+        ),
+        2,
     )
     mp.build_graph(circuit)  # identity mapping -> two distinct angles
     return mp, circuit
@@ -76,7 +77,7 @@ def test_non_contiguous_mapping_raises(self):
             ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2), index=2),
         )
         with pytest.raises(ValueError, match="contiguous"):
-            Circuit(gates)
+            Circuit(gates, 2)
 
     def test_mixed_param_scheme_rejected(self):
         """Setting `index` on some gates but not others is rejected as ambiguous."""
@@ -85,7 +86,7 @@ def test_mixed_param_scheme_rejected(self):
             ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),
         )
         with pytest.raises(ValueError, match="every gate must set"):
-            Circuit(gates)
+            Circuit(gates, 2)
 
     def test_shared_mapping_index_ties_gates(self, serial_comm):
         """Repeating an index in the mapping ties gates to one angle (one parameter)."""
@@ -96,6 +97,7 @@ def test_shared_mapping_index_ties_gates(self, serial_comm):
                 ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2), index=0),
                 ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2), index=0),
             ),
+            2,
         )
         mp.build_graph(circuit)
         assert mp.graph_layers == 2
@@ -109,13 +111,16 @@ def test_functional_invalidated_after_graph_mutation(self, serial_comm):
                 (
                     ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),
                     ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),
-                )
-            )
+                ),
+                2,
+            ),
         )
         functional = mp.expectation_value_functional()
         # Appending another layer mutates the graph, so the previously-built functional
         # must reject being called against the stale plan.
-        mp.build_graph(Circuit((ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2)),)))
+        mp.build_graph(
+            Circuit((ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2)),), 2)
+        )
         # Call with the parameter count the functional was built with (2), so the
         # stale-graph guard fires rather than the parameter-length check.
         with pytest.raises(RuntimeError, match=r"MP object has been modified"):
diff --git a/tests/test_pauli.py b/tests/test_pauli.py
index 9c029536..a6aff845 100644
--- a/tests/test_pauli.py
+++ b/tests/test_pauli.py
@@ -49,7 +49,9 @@ def test_num_qubits(self, serial_comm):
     def test_non_hermitian_pauli_gate_rejected(self, serial_comm):
         """An ExpGate with a complex (non-Hermitian) Pauli coefficient is rejected."""
         circuit = Circuit(
-            (ExpGate(PauliOperator({Pauli("X", 0): 1.0j}, num_qubits=1)),),
+            (ExpGate(PauliOperator({Pauli("X", 0): 1.0j}, num_qubits=2)),),
+            initial_state=(),
+            num_modes=2,
             parameters=(0.3,),
         )
         with pytest.raises(ValueError, match="not Hermitian"):
@@ -183,10 +185,9 @@ def test_all_valid_pauli_chars(self):
         assert set(op.terms) == {Pauli("XYIZ")}
 
     def test_get_majorana_operator_requires_num_qubits(self):
-        """Converting to Majorana without a qubit count raises a clear ValueError."""
-        op = PauliOperator._from_terms(["X"], [1.0], num_qubits=None)
-        with pytest.raises(ValueError, match="needs num_qubits"):
-            op.get_majorana_operator()
+        """Constructing with no qubit count raises a clear TypeError."""
+        with pytest.raises(TypeError, match="num_qubits must be an integer"):
+            PauliOperator._from_terms(["X"], [1.0], num_qubits=None)
 
     def test_str_few_terms(self):
         op = PauliOperator({"XY": 1.0}, num_qubits=2)
@@ -289,16 +290,18 @@ def _make_gate(self):
 
     def test_basic_construction(self):
         gates = (self._make_gate(), self._make_gate())
-        circuit = Circuit(gates, parameters=(0.5, 0.5), initial_state=(0,))
+        circuit = Circuit(gates, initial_state=(0,), num_modes=1, parameters=(0.5, 0.5))
         assert len(circuit) == 2
         assert circuit.initial_state == (0,)
 
     def test_empty_gates(self):
-        circuit = Circuit((), initial_state=(0,))
+        circuit = Circuit((), initial_state=(0,), num_modes=1)
         assert len(circuit) == 0
 
     def test_default_mapping_is_identity(self):
-        circuit = Circuit((self._make_gate(), self._make_gate()))
+        circuit = Circuit(
+            (self._make_gate(), self._make_gate()), initial_state=(), num_modes=1
+        )
         assert list(circuit.resolved_mapping) == [0, 1]
         assert circuit.n_parameters == 2
 
@@ -307,9 +310,11 @@ def test_rejects_mixed_gate_families(self):
         with pytest.raises(TypeError, match="mix"):
             Circuit(
                 (
-                    ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=1)),
+                    ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=2)),
                     ExpGate(MajoranaOperator({(0, 1): 1.0}, num_modes=2)),
-                )
+                ),
+                initial_state=(),
+                num_modes=2,
             )
 
     def test_pauli_gate_equality(self):
diff --git a/tests/test_qiskit_conversion.py b/tests/test_qiskit_conversion.py
index 68444a96..81ef4038 100644
--- a/tests/test_qiskit_conversion.py
+++ b/tests/test_qiskit_conversion.py
@@ -190,8 +190,9 @@ class ToQiskitCircuitCases:
     def case_single_gate(self):
         circuit = Circuit(
             (ExpGate(PauliOperator({Pauli("Z", 0): 1.0}, num_qubits=1)),),
-            parameters=(0.7,),
             initial_state=(),
+            num_modes=1,
+            parameters=(0.7,),
         )
         expected_circuit = QuantumCircuit(1)
         expected_circuit.append(
@@ -204,8 +205,9 @@ def case_single_gate(self):
     def case_local_gate(self):
         circuit = Circuit(
             (ExpGate(PauliOperator({Pauli("ZXY", (3, 1, 2)): 1.5}, num_qubits=5)),),
-            parameters=(0.7,),
             initial_state=(),
+            num_modes=5,
+            parameters=(0.7,),
         )
         expected_circuit = QuantumCircuit(5)
         # qiskit uses reversed Pauli string order and sorted gate qubit indices.
@@ -233,8 +235,9 @@ def case_single_pauli_evolution_gate(self):
         circuit.append(PauliEvolutionGate(operator, time=0.7), [0])
         expected = Circuit(
             gates=(ExpGate(PauliOperator({Pauli("Z", 0): 1.0}, num_qubits=1)),),
-            parameters=(0.7,),
             initial_state=(),
+            num_modes=1,
+            parameters=(0.7,),
         )
         return circuit, expected
 
@@ -249,8 +252,9 @@ def case_multiple_pauli_evolution_gates(self):
                 ExpGate(PauliOperator({Pauli("ZX", (0, 1)): 1.0}, num_qubits=2)),
                 ExpGate(PauliOperator({Pauli("Y", 0): 0.5}, num_qubits=2)),
             ),
-            parameters=(0.3, 0.5),
             initial_state=(),
+            num_modes=2,
+            parameters=(0.3, 0.5),
         )
         return circuit, expected
 
@@ -265,8 +269,9 @@ def case_rotation_gates_equivalent_to_pauli_evolution(self):
                 ExpGate(PauliOperator({Pauli("Y", 0): 0.5}, num_qubits=1)),
                 ExpGate(PauliOperator({Pauli("Z", 0): 0.5}, num_qubits=1)),
             ),
-            parameters=(0.5, 0.3, 0.7),
             initial_state=(),
+            num_modes=1,
+            parameters=(0.5, 0.3, 0.7),
         )
         return circuit, expected
 
@@ -277,8 +282,9 @@ def case_barrier_ignored(self):
         circuit.barrier()
         expected = Circuit(
             gates=(ExpGate(PauliOperator({Pauli("Z", 0): 1.0}, num_qubits=1)),),
-            parameters=(0.7,),
             initial_state=(),
+            num_modes=1,
+            parameters=(0.7,),
         )
         return circuit, expected
 
@@ -316,7 +322,9 @@ def test_non_commuting_generator_rejected(self):
 def test_to_qiskit_circuit_rejects_unbound() -> None:
     """to_qiskit_circuit on an unbound circuit raises a clear error, not an IndexError."""
     circuit = Circuit(
-        gates=(ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=1)),)
+        gates=(ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=1)),),
+        initial_state=(),
+        num_modes=1,
     )  # no parameter values
     with pytest.raises(ValueError, match="bound circuit"):
         to_qiskit_circuit(circuit, num_qubits=1)
@@ -333,6 +341,8 @@ def test_from_to_qiskit_circuit_roundtrip() -> None:
     """
     circuit = Circuit(
         gates=(ExpGate(PauliOperator({Pauli("XYZ", (3, 1, 2)): 1.0}, num_qubits=4)),),
+        initial_state=[],
+        num_modes=4,
         parameters=[-1.2],
     )  # no parameter values
     qcirc = to_qiskit_circuit(circuit, num_qubits=4)
diff --git a/tests/test_update_methods.py b/tests/test_update_methods.py
index 81b4cbf1..e3711f9c 100644
--- a/tests/test_update_methods.py
+++ b/tests/test_update_methods.py
@@ -107,7 +107,7 @@ def test_update_cutoff_valid(self, mp):
         mp.cutoff = 6
         # A weight-6 monomial takes an imaginary Hermitian coefficient (like weight-2).
         gate = ExpGate(MajoranaOperator({(0, 1, 2, 3, 4, 5): 1.0j}, num_modes=4))
-        mp.build_graph(Circuit((gate,)))
+        mp.build_graph(Circuit((gate,), initial_state=(), num_modes=4))
         assert mp.size() > 0
 
     def test_update_cutoff_invalid(self, mp):

From 9592bc87e7709a8b5d934aef27b003e5ad79e485 Mon Sep 17 00:00:00 2001
From: Adam <aglos@iitis.pl>
Date: Mon, 27 Jul 2026 21:27:46 +0000
Subject: [PATCH 3/6] refactor: :truck: num_modes rename to system_size for
 circuit

---
 benches/_builders.py                          | 10 ++-
 .../monoprop_hubbard1d_benchmark.py           |  1 +
 src/monoprop/circuit.py                       | 46 +++++-----
 src/monoprop/monomial_propagator.py           |  4 +-
 src/monoprop/qiskit_conversion.py             |  8 +-
 tests/cases.py                                |  2 +-
 tests/test_circuit.py                         | 84 +++++++++----------
 tests/test_coeff_trunc.py                     |  6 +-
 tests/test_fermi.py                           | 10 +--
 tests/test_majorana.py                        |  2 +-
 tests/test_monoprop_trivial.py                |  6 +-
 tests/test_nonfermi.py                        |  2 +-
 tests/test_only_rotate_k.py                   | 10 +--
 tests/test_pauli.py                           | 12 +--
 tests/test_qiskit_conversion.py               | 16 ++--
 tests/test_update_methods.py                  |  4 +-
 16 files changed, 116 insertions(+), 107 deletions(-)

diff --git a/benches/_builders.py b/benches/_builders.py
index 762dca75..9657ef7e 100644
--- a/benches/_builders.py
+++ b/benches/_builders.py
@@ -167,7 +167,7 @@ def make_random_problem(
         majoranas=gen_majoranas,
         gen_coeffs=gen_coeffs,
         param_inds=param_inds,
-        num_modes=num_modes,
+        system_size=num_modes,
         parameters=parameters,
         initial_state=[],
     )
@@ -294,7 +294,12 @@ def build_hubbard_problem(
     fermi_gates = [ExpGate(term) for term in _hubbard_fermion_terms(config)]
     parameters = [config.trotter_dt] * len(fermi_gates)
     occupied = _neel_occupied_modes(config.num_sites, config.neel_start_spin)
-    circuit = Circuit(gates=fermi_gates, parameters=parameters, initial_state=occupied)
+    circuit = Circuit(
+        gates=fermi_gates,
+        parameters=parameters,
+        initial_state=occupied,
+        system_size=config.num_qubits,
+    )
 
     observable = FermiOperator(
         terms=[
@@ -526,6 +531,7 @@ def build_kicked_ising_problem(
         gates=tuple(gate for gate, _ in gate_angles),
         parameters=tuple(angle for _, angle in gate_angles),
         initial_state=[],
+        system_size=config.num_qubits,
     )
 
     obs_str = (
diff --git a/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py b/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py
index 0ea13926..638cd167 100644
--- a/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py
+++ b/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py
@@ -157,6 +157,7 @@ def main():
         gates=trotter_gates,
         parameters=trotter_parameters,
         initial_state=intial_state,
+        system_size=num_qubits,
     )
 
     simulator = MajoranaPropagator(
diff --git a/src/monoprop/circuit.py b/src/monoprop/circuit.py
index e026fcdc..23a0ce0c 100644
--- a/src/monoprop/circuit.py
+++ b/src/monoprop/circuit.py
@@ -256,7 +256,7 @@ class Circuit:
         gates: The ordered exponential gates.
         parameters: The angle values, or empty for an unbound circuit.
         initial_state: The reference state (occupied mode / qubit indices).
-        num_modes: System width (number of fermionic modes / qubits).
+        system_size: System width (number of fermionic modes / qubits).
         family: The gate family -- ``"pauli"``, ``"majorana"``, or ``"empty"`` -- computed at
             construction; the propagators dispatch on it.
     """
@@ -264,7 +264,7 @@ class Circuit:
     def __init__(  # noqa:C901
         self,
         gates: Sequence[ExpGate],
-        num_modes: int,
+        system_size: int,
         initial_state: Sequence[int] = (),
         parameters: Sequence[float] = (),
     ) -> None:
@@ -273,24 +273,24 @@ def __init__(  # noqa:C901
         Args:
             gates: The ordered exponential gates.
             initial_state: The reference state (occupied mode / qubit indices).
-            num_modes: Number of modes/qubits for the circuit.
+            system_size: Number of modes/qubits for the circuit.
             parameters: The angle values, or empty for an unbound circuit.
 
         Raises:
             ValueError: On duplicate initial-state indices, a bad parameter mapping, or a
                 bound circuit whose parameter count does not match [n_parameters][]; also if
-                a gate's operator width differs from ``num_modes``.
+                a gate's operator width differs from ``system_size``.
             TypeError: On a non-[ExpGate][] gate or a mix of qubit and Majorana gate families.
         """
         gates = tuple(gates)
         parameters = tuple(float(v) for v in parameters)
         initial_state = tuple(int(i) for i in initial_state)
-        num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        system_size = _validate_system_size(system_size, argument_name="system_size")
         if len(set(initial_state)) != len(initial_state):
             raise ValueError("Duplicate indices in initial state")
-        if any(i < 0 or i >= num_modes for i in initial_state):
+        if any(i < 0 or i >= system_size for i in initial_state):
             raise ValueError(
-                f"initial_state entries must be in 0..{num_modes - 1}; got {list(initial_state)}."
+                f"initial_state entries must be in 0..{system_size - 1}; got {list(initial_state)}."
             )
 
         # Validate gate types up front: the identity-drop below reads gate.index/.generator,
@@ -306,9 +306,9 @@ def __init__(  # noqa:C901
                 if isinstance(gate.generator, PauliOperator)
                 else gate.generator.num_modes
             )
-            if gate_size != num_modes:
+            if gate_size != system_size:
                 raise ValueError(
-                    f"Gate generator width {gate_size} does not match circuit num_modes={num_modes}."
+                    f"Gate generator width {gate_size} does not match circuit system_size={system_size}."
                 )
 
         def _is_identity_gate(gate: ExpGate) -> bool:
@@ -329,7 +329,7 @@ def _is_identity_gate(gate: ExpGate) -> bool:
         self.gates = gates
         self.parameters = parameters
         self.initial_state = initial_state
-        self.num_modes = num_modes
+        self.system_size = system_size
         #: The gate family, computed from the (validated) gates; the propagators dispatch on it.
         self.family = self._resolve_family(gates)
 
@@ -348,7 +348,7 @@ def __eq__(self, other: object) -> bool:
             self.gates == other.gates
             and self.parameters == other.parameters
             and self.initial_state == other.initial_state
-            and self.num_modes == other.num_modes
+            and self.system_size == other.system_size
         )
 
     __hash__ = None  # type: ignore[assignment]  # value-equal but not hashable (mutable gates)
@@ -358,7 +358,7 @@ def __repr__(self) -> str:
         return (
             f"{self.__class__.__name__}(gates={self.gates!r}, "
             f"parameters={self.parameters!r}, initial_state={self.initial_state!r}, "
-            f"num_modes={self.num_modes!r})"
+            f"system_size={self.system_size!r})"
         )
 
     @staticmethod
@@ -442,10 +442,10 @@ def __add__(self, other: Circuit) -> Circuit:
             raise ValueError(
                 "Cannot concatenate circuits with different initial states."
             )
-        if self.num_modes != other.num_modes:
+        if self.system_size != other.system_size:
             raise ValueError(
-                f"Cannot concatenate circuits with different num_modes: "
-                f"{self.num_modes} != {other.num_modes}."
+                f"Cannot concatenate circuits with different system_size: "
+                f"{self.system_size} != {other.system_size}."
             )
         offset = self.n_parameters
         # Preserve each gate's _structural flag: a dense (wire-format) gate carries structural
@@ -462,7 +462,7 @@ def __add__(self, other: Circuit) -> Circuit:
             gates=left + right,
             parameters=tuple(self.parameters) + tuple(other.parameters),
             initial_state=self.initial_state or other.initial_state,
-            num_modes=self.num_modes,
+            system_size=self.system_size,
         )
 
     @classmethod
@@ -471,7 +471,7 @@ def from_dense_arrays(
         majoranas: Sequence[Sequence[int]],
         gen_coeffs: Sequence[float],
         param_inds: Sequence[int],
-        num_modes: int,
+        system_size: int,
         parameters: Sequence[float] = (),
         initial_state: Sequence[int] = (),
     ) -> Circuit:
@@ -489,7 +489,7 @@ def from_dense_arrays(
             gen_coeffs: Generator coefficient per monomial.
             param_inds: Variational-angle index per monomial (contiguous runs group into
                 gates).
-            num_modes: Number of fermionic modes in the system.
+            system_size: Number of fermionic modes in the system.
             parameters: Optional angle values.
             initial_state: Optional reference state (occupied mode indices).
 
@@ -497,7 +497,7 @@ def from_dense_arrays(
             A [Circuit][] carrying the grouped gates, angle values, and initial state.
         """
         indices = [int(p) for p in param_inds]
-        num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        system_size = _validate_system_size(system_size, argument_name="system_size")
         gates: list[ExpGate] = []
         current_index: int | None = None
         current_majoranas: list[tuple[int, ...]] = []
@@ -509,7 +509,7 @@ def _flush() -> None:
             gates.append(
                 ExpGate._structural_gate(
                     MajoranaOperator._from_terms(
-                        current_majoranas, current_coeffs, num_modes=num_modes
+                        current_majoranas, current_coeffs, num_modes=system_size
                     ),
                     index=current_index,
                 )
@@ -526,9 +526,9 @@ def _flush() -> None:
                 raise ValueError(
                     f"Majorana indices must be non-negative; got {majorana}."
                 )
-            if majorana and max(majorana) >= 2 * num_modes:
+            if majorana and max(majorana) >= 2 * system_size:
                 raise ValueError(
-                    f"Majorana term {majorana} acts on an index >= 2*num_modes={2 * num_modes}."
+                    f"Majorana term {majorana} acts on an index >= 2*system_size={2 * system_size}."
                 )
             current_majoranas.append(majorana)
             current_coeffs.append(complex(float(coeff)))
@@ -539,7 +539,7 @@ def _flush() -> None:
             gates=tuple(gates),
             parameters=tuple(float(p) for p in parameters),
             initial_state=tuple(int(i) for i in initial_state),
-            num_modes=num_modes,
+            system_size=system_size,
         )
 
 
diff --git a/src/monoprop/monomial_propagator.py b/src/monoprop/monomial_propagator.py
index 35f29a23..44c6c405 100644
--- a/src/monoprop/monomial_propagator.py
+++ b/src/monoprop/monomial_propagator.py
@@ -199,9 +199,9 @@ def _check_initial_state(self, circuit: Circuit) -> None:
     def _check_circuit_width(self, circuit: Circuit) -> None:
         """Reject a circuit with a system width that disagrees with the propagator."""
         expected = self._num_qubits if self._num_qubits is not None else self._num_modes
-        if circuit.num_modes != expected:
+        if circuit.system_size != expected:
             raise ValueError(
-                f"Circuit num_modes={circuit.num_modes} does not match propagator width "
+                f"Circuit system_size={circuit.system_size} does not match propagator width "
                 f"{expected}."
             )
 
diff --git a/src/monoprop/qiskit_conversion.py b/src/monoprop/qiskit_conversion.py
index 5156c171..b660f8a1 100644
--- a/src/monoprop/qiskit_conversion.py
+++ b/src/monoprop/qiskit_conversion.py
@@ -169,7 +169,7 @@ def from_qiskit_circuit(
         gates=tuple(gates),
         parameters=tuple(parameters),
         initial_state=tuple(initial_state),
-        num_modes=num_qubits,
+        system_size=num_qubits,
     )
 
 
@@ -199,15 +199,15 @@ def to_qiskit_circuit(circuit: Circuit, num_qubits: int) -> QuantumCircuit:
 
     Args:
         circuit: A [Circuit][monoprop.circuit.Circuit] representing the given circuit.
-        num_qubits: Total number of qubits. Must match ``circuit.num_modes``.
+        num_qubits: Total number of qubits. Must match ``circuit.system_size``.
 
     Returns:
         A qiskit quantum circuit.
     """
     num_qubits = _validate_system_size(num_qubits, argument_name="num_qubits")
-    if num_qubits != circuit.num_modes:
+    if num_qubits != circuit.system_size:
         raise ValueError(
-            f"num_qubits={num_qubits} does not match circuit.num_modes={circuit.num_modes}."
+            f"num_qubits={num_qubits} does not match circuit.system_size={circuit.system_size}."
         )
     if len(circuit.parameters) != circuit.n_parameters:
         raise ValueError(
diff --git a/tests/cases.py b/tests/cases.py
index b4406431..b049b475 100644
--- a/tests/cases.py
+++ b/tests/cases.py
@@ -52,7 +52,7 @@ def to_circuit(self) -> Circuit:
             majoranas=self.majoranas,
             gen_coeffs=self.gen_coeffs,
             param_inds=self.param_inds,
-            num_modes=self.num_modes,
+            system_size=self.num_modes,
             parameters=self.parameters,
             initial_state=self.initial_state,
         )
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index 9ffefc56..4d57b94b 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -162,11 +162,11 @@ def test_exp_gate_applies_atol_truncation(
 def test_circuit_equality() -> None:
     """Circuits are equal on gates/parameters/initial_state; family is derived, not compared."""
     gen = MajoranaOperator({(0, 1): 1.0j}, num_modes=2)
-    a = Circuit((ExpGate(gen),), initial_state=(0,), num_modes=2, parameters=(0.3,))
-    b = Circuit((ExpGate(gen),), initial_state=(0,), num_modes=2, parameters=(0.3,))
+    a = Circuit((ExpGate(gen),), initial_state=(0,), system_size=2, parameters=(0.3,))
+    b = Circuit((ExpGate(gen),), initial_state=(0,), system_size=2, parameters=(0.3,))
     assert a == b
     assert a != Circuit(
-        (ExpGate(gen),), initial_state=(0,), num_modes=2, parameters=(0.9,)
+        (ExpGate(gen),), initial_state=(0,), system_size=2, parameters=(0.9,)
     )
     assert a != "not a circuit"
 
@@ -174,7 +174,7 @@ def test_circuit_equality() -> None:
 def test_circuit_rejects_non_exp_gate() -> None:
     """A gate that is not an ExpGate is rejected with a clear TypeError."""
     with pytest.raises(TypeError, match="Circuit gates must be ExpGate"):
-        Circuit(("not a gate",), initial_state=(), num_modes=0)  # type: ignore[arg-type]
+        Circuit(("not a gate",), initial_state=(), system_size=0)  # type: ignore[arg-type]
 
 
 def test_to_circuit_round_trips_sequence() -> None:
@@ -201,7 +201,7 @@ def test_default_mapping_is_identity() -> None:
             ExpGate(MajoranaOperator({(2, 3): 1.0j}, num_modes=2)),
         ),
         initial_state=(),
-        num_modes=2,
+        system_size=2,
     )
     assert circuit.resolved_mapping == (0, 1)
     assert circuit.n_parameters == 2
@@ -216,7 +216,7 @@ def test_shared_mapping_index_ties_gates() -> None:
             MajoranaOperator({(0, 3): 1.0j}, num_modes=2), index=0
         ),  # ties to the first
     )
-    circuit = Circuit(gates, initial_state=(), num_modes=2)
+    circuit = Circuit(gates, initial_state=(), system_size=2)
     assert circuit.resolved_mapping == (0, 1, 0)
     assert circuit.n_parameters == 2
 
@@ -228,7 +228,7 @@ def test_circuit_rejects_non_contiguous_mapping() -> None:
         ExpGate(MajoranaOperator({(2, 3): 1.0j}, num_modes=2), index=2),
     )
     with pytest.raises(ValueError, match="contiguous"):
-        Circuit(gates, initial_state=(), num_modes=2)
+        Circuit(gates, initial_state=(), system_size=2)
 
 
 def test_circuit_rejects_mixed_param_scheme() -> None:
@@ -238,14 +238,14 @@ def test_circuit_rejects_mixed_param_scheme() -> None:
         ExpGate(MajoranaOperator({(2, 3): 1.0j}, num_modes=2)),
     )
     with pytest.raises(ValueError, match="every gate must set"):
-        Circuit(gates, initial_state=(), num_modes=2)
+        Circuit(gates, initial_state=(), system_size=2)
 
 
 def test_circuit_rejects_wrong_parameter_length() -> None:
     """A bound circuit must supply exactly one value per distinct angle."""
     gates = (ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),)
     with pytest.raises(ValueError, match="1 parameters"):
-        Circuit(gates, initial_state=(), num_modes=2, parameters=(0.1, 0.2))
+        Circuit(gates, initial_state=(), system_size=2, parameters=(0.1, 0.2))
 
 
 def test_circuit_add_offsets_second_axis() -> None:
@@ -256,13 +256,13 @@ def test_circuit_add_offsets_second_axis() -> None:
             ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),
         ),
         initial_state=(),
-        num_modes=2,
+        system_size=2,
         parameters=(0.1, 0.2),
     )
     b = Circuit(
         (ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2)),),
         initial_state=(),
-        num_modes=2,
+        system_size=2,
         parameters=(0.3,),
     )
     combined = a + b
@@ -276,12 +276,12 @@ def test_circuit_add_rejects_mixed_families() -> None:
     maj = Circuit(
         (ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=2)),),
         initial_state=(),
-        num_modes=2,
+        system_size=2,
     )
     qubit = Circuit(
         (ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=2)),),
         initial_state=(),
-        num_modes=2,
+        system_size=2,
     )
     with pytest.raises(TypeError, match="gate families differ"):
         _ = maj + qubit
@@ -292,12 +292,12 @@ def test_circuit_add_rejects_different_initial_states() -> None:
     a = Circuit(
         (ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),),
         initial_state=(0,),
-        num_modes=2,
+        system_size=2,
     )
     b = Circuit(
         (ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),),
         initial_state=(1,),
-        num_modes=2,
+        system_size=2,
     )
     with pytest.raises(ValueError, match="different initial states"):
         _ = a + b
@@ -313,7 +313,7 @@ def test_bound_circuit_with_identity_gate_wrong_param_count_raises() -> None:
         Circuit(
             gates,
             initial_state=(),
-            num_modes=2,
+            system_size=2,
             parameters=(0.1,),
         )  # 1 value, but 2 gates before the drop
 
@@ -330,7 +330,7 @@ def test_non_hermitian_majorana_generator_rejected() -> None:
     bad = Circuit(
         (ExpGate(MajoranaOperator({(0, 1): 1.0}, num_modes=2)),),
         initial_state=(),
-        num_modes=2,
+        system_size=2,
         parameters=(0.3,),
     )
     with pytest.raises(ValueError, match="not Hermitian"):
@@ -349,11 +349,11 @@ def test_hermitian_majorana_generator_matches_structural() -> None:
     hermitian = Circuit(
         gates=(ExpGate(MajoranaOperator({(4, 5): 1j}, num_modes=8)),),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.5,),
     )
     structural = Circuit.from_dense_arrays(
-        [[4, 5]], [-1.0], [0], num_modes=8, parameters=[0.5]
+        [[4, 5]], [-1.0], [0], system_size=8, parameters=[0.5]
     )
 
     from_hermitian = MajoranaPropagator.from_circuit(
@@ -490,7 +490,7 @@ def _multi_term_gate_propagator():
     # two monomials -> two layers
     g0 = ExpGate(MajoranaOperator({(0, 2): 1.0j, (1, 3): 1.0j}, num_modes=2))
     g1 = ExpGate(MajoranaOperator({(2,): 1.0}, num_modes=2))
-    prop.build_graph(Circuit((g0, g1), initial_state=(), num_modes=2))
+    prop.build_graph(Circuit((g0, g1), initial_state=(), system_size=2))
     return prop
 
 
@@ -509,7 +509,7 @@ def test_n_gates_accumulates_across_builds() -> None:
         Circuit(
             (ExpGate(MajoranaOperator({(0,): 1.0}, num_modes=2)),),
             initial_state=(),
-            num_modes=2,
+            system_size=2,
         )
     )
     assert prop.n_gates == 1
@@ -517,7 +517,7 @@ def test_n_gates_accumulates_across_builds() -> None:
         Circuit(
             (ExpGate(MajoranaOperator({(1,): 1.0}, num_modes=2)),),
             initial_state=(),
-            num_modes=2,
+            system_size=2,
         )
     )
     assert prop.n_gates == 2
@@ -545,7 +545,7 @@ def test_majorana_propagator_rejects_pauli_circuit() -> None:
     circuit = Circuit(
         gates=(ExpGate(PauliOperator({"Z": 1.0}, num_qubits=problem.n_modes)),),
         initial_state=(),
-        num_modes=problem.n_modes,
+        system_size=problem.n_modes,
     )
 
     with pytest.raises(TypeError, match="qubit"):
@@ -560,7 +560,7 @@ def test_propagate_rejects_mismatched_initial_state() -> None:
     circuit = Circuit(
         (gate,),
         initial_state=(0, 1),
-        num_modes=problem.n_modes,
+        system_size=problem.n_modes,
         parameters=(0.1,),
     )
 
@@ -574,7 +574,7 @@ def test_propagate_accepts_empty_initial_state() -> None:
     prop = _propagator(problem)
     gate = ExpGate(MajoranaOperator({(0, 1): 1.0j}, num_modes=problem.n_modes))
     circuit = Circuit(
-        (gate,), initial_state=(), num_modes=problem.n_modes, parameters=(0.1,)
+        (gate,), initial_state=(), system_size=problem.n_modes, parameters=(0.1,)
     )  # empty initial_state
 
     prop.propagate(circuit)  # does not raise
@@ -605,11 +605,11 @@ def test_build_graph_accumulates_layers_and_parameters(fixture: str) -> None:
 
     twice = _propagator(problem)
     twice.build_graph(
-        Circuit(_rebase(gates[:split]), initial_state=(), num_modes=problem.n_modes)
+        Circuit(_rebase(gates[:split]), initial_state=(), system_size=problem.n_modes)
     )
     layers_after_first = twice.graph_layers
     twice.build_graph(
-        Circuit(_rebase(gates[split:]), initial_state=(), num_modes=problem.n_modes)
+        Circuit(_rebase(gates[split:]), initial_state=(), system_size=problem.n_modes)
     )
 
     assert 0 < layers_after_first < twice.graph_layers
@@ -633,13 +633,13 @@ def test_compose_then_single_build_matches_single_call(fixture: str) -> None:
     a = Circuit(
         _rebase(gates[:split]),
         initial_state=(),
-        num_modes=problem.n_modes,
+        system_size=problem.n_modes,
         parameters=tuple(params[:split]),
     )
     b = Circuit(
         _rebase(gates[split:]),
         initial_state=(),
-        num_modes=problem.n_modes,
+        system_size=problem.n_modes,
         parameters=tuple(params[split:]),
     )
     composed = a + b
@@ -674,10 +674,10 @@ def test_build_graph_in_two_calls_schrodinger(fixture: str) -> None:
 
     twice = _schrodinger_propagator(problem)
     twice.build_graph(
-        Circuit(_rebase(gates[:split]), initial_state=(), num_modes=problem.n_modes)
+        Circuit(_rebase(gates[:split]), initial_state=(), system_size=problem.n_modes)
     )
     twice.build_graph(
-        Circuit(_rebase(gates[split:]), initial_state=(), num_modes=problem.n_modes)
+        Circuit(_rebase(gates[split:]), initial_state=(), system_size=problem.n_modes)
     )
 
     np.testing.assert_allclose(
@@ -697,12 +697,12 @@ def test_build_graph_twice_with_seed_regeneration(fixture: str) -> None:
 
     prop = _schrodinger_propagator(problem)
     prop.build_graph(
-        Circuit(_rebase(gates[:split]), initial_state=(), num_modes=problem.n_modes)
+        Circuit(_rebase(gates[:split]), initial_state=(), system_size=problem.n_modes)
     )
     # seed_parameters on the second call exercises the internal seed regeneration
     # (the former operator_coeffs round-trip) used for coefficient-informed truncation.
     prop.build_graph(
-        Circuit(_rebase(gates[split:]), initial_state=(), num_modes=problem.n_modes),
+        Circuit(_rebase(gates[split:]), initial_state=(), system_size=problem.n_modes),
         seed_parameters=params,
     )
 
@@ -729,7 +729,7 @@ def test_empty_default_mapping_gate_dropped_and_evaluable() -> None:
             ),  # identity generator: dropped
         ),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.5, 0.3),
     )
     assert len(circuit.gates) == 1
@@ -753,7 +753,7 @@ def test_empty_gate_in_middle_builds_contiguously() -> None:
             ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),
         ),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.5, 0.3, 0.2),
     )
     assert len(circuit.gates) == 2
@@ -771,7 +771,7 @@ def test_surplus_parameters_raise_not_truncated() -> None:
         Circuit(
             gates=(ExpGate(generator),),
             initial_state=(),
-            num_modes=4,
+            system_size=4,
             parameters=(1.0, 2.0),
         )
 
@@ -828,7 +828,7 @@ def test_build_graph_seed_parameters_accepts_numpy() -> None:
             ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),
         ),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.5, 0.3),
     )
     prop = _small_propagator(lower_atol=1e-12)
@@ -843,7 +843,7 @@ def test_build_graph_rejects_too_short_seed() -> None:
             ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),
         ),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.5, 0.3),
     )
     prop = _small_propagator()
@@ -858,13 +858,13 @@ def test_extend_without_seed_builds_structurally() -> None:
     c1 = Circuit(
         gates=(ExpGate(MajoranaOperator({(4, 5): -1.0j}, num_modes=8)),),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.3,),
     )
     c2 = Circuit(
         gates=(ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.4,),
     )
     params = [0.3, 0.4]
@@ -890,13 +890,13 @@ def test_propagate_after_build_graph_rejected() -> None:
     c1 = Circuit(
         gates=(ExpGate(MajoranaOperator({(4, 5): -1.0j}, num_modes=8)),),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.3,),
     )
     c2 = Circuit(
         gates=(ExpGate(MajoranaOperator({(2, 3): -1.0j}, num_modes=8)),),
         initial_state=(),
-        num_modes=8,
+        system_size=8,
         parameters=(0.4,),
     )
     prop = _small_propagator()
diff --git a/tests/test_coeff_trunc.py b/tests/test_coeff_trunc.py
index 3cf50378..ed58ee20 100644
--- a/tests/test_coeff_trunc.py
+++ b/tests/test_coeff_trunc.py
@@ -62,7 +62,7 @@ def test_coeff_trunc(serial_comm):
         parameters=[np.pi / 6],
         gen_coeffs=[1.0],
         param_inds=[0],
-        num_modes=n_modes,
+        system_size=n_modes,
     )
     circuit = sequence
 
@@ -147,7 +147,7 @@ def test_evolution_coeff_trunc_no_atols(serial_comm):
         parameters=[p],
         gen_coeffs=[1.0],
         param_inds=[0],
-        num_modes=n_modes,
+        system_size=n_modes,
     )
     circuit = sequence
 
@@ -185,7 +185,7 @@ def test_evolution_coeff_trunc_small_coeffs(serial_comm):
         parameters=[p],
         gen_coeffs=[1.0],
         param_inds=[0],
-        num_modes=n_modes,
+        system_size=n_modes,
     )
     circuit = sequence
 
diff --git a/tests/test_fermi.py b/tests/test_fermi.py
index e6568dd7..629007bd 100644
--- a/tests/test_fermi.py
+++ b/tests/test_fermi.py
@@ -323,7 +323,7 @@ def test_len(self):
         circuit = Circuit(
             gates=gates,
             initial_state=[0],
-            num_modes=1,
+            system_size=1,
             parameters=[0.1, 0.2],
         )
 
@@ -340,7 +340,7 @@ def test_converts_fermi_gates_to_majorana(self):
         circuit = Circuit(
             gates=[gate_0, gate_1],
             initial_state=[0, 1],
-            num_modes=2,
+            system_size=2,
             parameters=[0.3, -0.7],
         )
 
@@ -352,7 +352,7 @@ def test_converts_fermi_gates_to_majorana(self):
         assert all(g.family == "majorana" for g in circuit.gates)
 
         majoranas, gen_coeffs, per_monomial_mapping, gate_indices = expand_monomials(
-            circuit.gates, circuit.resolved_mapping, circuit.num_modes
+            circuit.gates, circuit.resolved_mapping, circuit.system_size
         )
         # One gate per fermi gate; each generator here has two monomials.
         n_terms = len(gate_0.generator.terms)
@@ -370,7 +370,7 @@ def test_drops_identity_generators_and_aligned_parameters(self):
         circuit = Circuit(
             gates=[ExpGate(_number_op()), ExpGate(identity), ExpGate(_number_op())],
             initial_state=[0],
-            num_modes=1,
+            system_size=1,
             parameters=[0.1, 0.2, 0.3],
         )
         assert len(circuit) == 2  # the identity gate is dropped
@@ -378,4 +378,4 @@ def test_drops_identity_generators_and_aligned_parameters(self):
 
     def test_validate_inputs_duplicate_initial_state_raises(self):
         with pytest.raises(ValueError, match="Duplicate indices in initial state"):
-            Circuit(gates=[ExpGate(_number_op())], initial_state=[0, 0], num_modes=1)
+            Circuit(gates=[ExpGate(_number_op())], initial_state=[0, 0], system_size=1)
diff --git a/tests/test_majorana.py b/tests/test_majorana.py
index f3853761..efbbe76b 100644
--- a/tests/test_majorana.py
+++ b/tests/test_majorana.py
@@ -64,7 +64,7 @@ def test_from_dense_arrays_groups_by_param_ind():
         majoranas=[(0, 1), (2, 3), (0, 3)],
         gen_coeffs=[0.5, -0.5, 1.0],
         param_inds=[0, 0, 1],
-        num_modes=2,
+        system_size=2,
         parameters=[1.0, 2.0],
         initial_state=[0, 1],
     )
diff --git a/tests/test_monoprop_trivial.py b/tests/test_monoprop_trivial.py
index 634a1800..86e608eb 100644
--- a/tests/test_monoprop_trivial.py
+++ b/tests/test_monoprop_trivial.py
@@ -40,7 +40,7 @@ def test_trivial_evolved_operator_cases(
     """Test trivial evolved operator dict for various initial conditions."""
     kwargs = {"schrodinger_cutoff": schrodinger_cutoff} if schrodinger_cutoff else {}
     quantum_circuit = Circuit(
-        initial_state=[], num_modes=initial_op.num_modes, gates=[]
+        initial_state=[], system_size=initial_op.num_modes, gates=[]
     )
     mp = MajoranaPropagator(
         initial_op,
@@ -56,7 +56,7 @@ def test_trivial_evolved_operator_cases(
 def test_trivial_evolved_operator(serial_comm):
     initial_op = MajoranaOperator({(0, 1, 2, 4): 1}, 8)
     quantum_circuit = Circuit(
-        initial_state=[], num_modes=initial_op.num_modes, gates=[]
+        initial_state=[], system_size=initial_op.num_modes, gates=[]
     )
     mp = MajoranaPropagator(
         initial_op, quantum_circuit.initial_state, cutoff=16, comm=serial_comm
@@ -106,7 +106,7 @@ def test_update_initial_operator(
 ):
     """Test updating coefficients in both regular and Schrodinger pictures."""
     kwargs = {"schrodinger_cutoff": schrodinger_cutoff} if schrodinger_cutoff else {}
-    quantum_circuit = Circuit(initial_state=[], num_modes=init_op.num_modes, gates=[])
+    quantum_circuit = Circuit(initial_state=[], system_size=init_op.num_modes, gates=[])
     mp = MajoranaPropagator(
         init_op,
         quantum_circuit.initial_state,
diff --git a/tests/test_nonfermi.py b/tests/test_nonfermi.py
index c534d940..9a81e29b 100644
--- a/tests/test_nonfermi.py
+++ b/tests/test_nonfermi.py
@@ -33,7 +33,7 @@ def test_nonfermi(serial_comm):
         majoranas=majoranas,
         gen_coeffs=gen_coeffs,
         param_inds=param_inds,
-        num_modes=num_modes,
+        system_size=num_modes,
         parameters=parameters,
         initial_state=[],
     )
diff --git a/tests/test_only_rotate_k.py b/tests/test_only_rotate_k.py
index 06ce66fc..8a7109fb 100644
--- a/tests/test_only_rotate_k.py
+++ b/tests/test_only_rotate_k.py
@@ -55,7 +55,7 @@ def test_basic_orbital_rotation(serial_comm):
         parameters=[np.pi / 4],
         gen_coeffs=[1.0],
         param_inds=[0],
-        num_modes=n_modes,
+        system_size=n_modes,
     )
     circuit = sequence
     kwargs = {"cutoff": 6, "schrodinger_cutoff": 8, "comm": serial_comm}
@@ -97,13 +97,13 @@ def test_only_rotate_len_k(problem, inplace, serial_mp_kwargs):
     non_orbital = Circuit(
         tuple(ExpGate._with_index(gate, None) for gate in non_orbital_gates),
         initial_state=(),
-        num_modes=problem.n_modes,
+        system_size=problem.n_modes,
         parameters=tuple(parameters[:split]),
     )
     orbital = Circuit(
         tuple(ExpGate._with_index(gate, None) for gate in orbital_gates),
         initial_state=(),
-        num_modes=problem.n_modes,
+        system_size=problem.n_modes,
         parameters=tuple(parameters[split:]),
     )
 
@@ -158,7 +158,7 @@ def test_only_rotate_len_k_errors_majorana(only_rotate_len_k, err, method_name):
     mp = MajoranaPropagator(MajoranaOperator({}, 4), [], cutoff=6, schrodinger_cutoff=8)
     with err:
         getattr(mp, method_name)(
-            Circuit((), initial_state=(), num_modes=4),
+            Circuit((), initial_state=(), system_size=4),
             only_rotate_len_k=only_rotate_len_k,
         )
 
@@ -196,6 +196,6 @@ def test_only_rotate_len_k_errors_pauli(only_rotate_len_k, err, method_name):
     mp = PauliPropagator(PauliOperator({}, 4), [], cutoff=6, schrodinger_cutoff=8)
     with err:
         getattr(mp, method_name)(
-            Circuit((), initial_state=(), num_modes=4),
+            Circuit((), initial_state=(), system_size=4),
             only_rotate_len_k=only_rotate_len_k,
         )
diff --git a/tests/test_pauli.py b/tests/test_pauli.py
index a6aff845..0a1df2ad 100644
--- a/tests/test_pauli.py
+++ b/tests/test_pauli.py
@@ -51,7 +51,7 @@ def test_non_hermitian_pauli_gate_rejected(self, serial_comm):
         circuit = Circuit(
             (ExpGate(PauliOperator({Pauli("X", 0): 1.0j}, num_qubits=2)),),
             initial_state=(),
-            num_modes=2,
+            system_size=2,
             parameters=(0.3,),
         )
         with pytest.raises(ValueError, match="not Hermitian"):
@@ -290,17 +290,19 @@ def _make_gate(self):
 
     def test_basic_construction(self):
         gates = (self._make_gate(), self._make_gate())
-        circuit = Circuit(gates, initial_state=(0,), num_modes=1, parameters=(0.5, 0.5))
+        circuit = Circuit(
+            gates, initial_state=(0,), system_size=1, parameters=(0.5, 0.5)
+        )
         assert len(circuit) == 2
         assert circuit.initial_state == (0,)
 
     def test_empty_gates(self):
-        circuit = Circuit((), initial_state=(0,), num_modes=1)
+        circuit = Circuit((), initial_state=(0,), system_size=1)
         assert len(circuit) == 0
 
     def test_default_mapping_is_identity(self):
         circuit = Circuit(
-            (self._make_gate(), self._make_gate()), initial_state=(), num_modes=1
+            (self._make_gate(), self._make_gate()), initial_state=(), system_size=1
         )
         assert list(circuit.resolved_mapping) == [0, 1]
         assert circuit.n_parameters == 2
@@ -314,7 +316,7 @@ def test_rejects_mixed_gate_families(self):
                     ExpGate(MajoranaOperator({(0, 1): 1.0}, num_modes=2)),
                 ),
                 initial_state=(),
-                num_modes=2,
+                system_size=2,
             )
 
     def test_pauli_gate_equality(self):
diff --git a/tests/test_qiskit_conversion.py b/tests/test_qiskit_conversion.py
index 81ef4038..3407fdb7 100644
--- a/tests/test_qiskit_conversion.py
+++ b/tests/test_qiskit_conversion.py
@@ -191,7 +191,7 @@ def case_single_gate(self):
         circuit = Circuit(
             (ExpGate(PauliOperator({Pauli("Z", 0): 1.0}, num_qubits=1)),),
             initial_state=(),
-            num_modes=1,
+            system_size=1,
             parameters=(0.7,),
         )
         expected_circuit = QuantumCircuit(1)
@@ -206,7 +206,7 @@ def case_local_gate(self):
         circuit = Circuit(
             (ExpGate(PauliOperator({Pauli("ZXY", (3, 1, 2)): 1.5}, num_qubits=5)),),
             initial_state=(),
-            num_modes=5,
+            system_size=5,
             parameters=(0.7,),
         )
         expected_circuit = QuantumCircuit(5)
@@ -236,7 +236,7 @@ def case_single_pauli_evolution_gate(self):
         expected = Circuit(
             gates=(ExpGate(PauliOperator({Pauli("Z", 0): 1.0}, num_qubits=1)),),
             initial_state=(),
-            num_modes=1,
+            system_size=1,
             parameters=(0.7,),
         )
         return circuit, expected
@@ -253,7 +253,7 @@ def case_multiple_pauli_evolution_gates(self):
                 ExpGate(PauliOperator({Pauli("Y", 0): 0.5}, num_qubits=2)),
             ),
             initial_state=(),
-            num_modes=2,
+            system_size=2,
             parameters=(0.3, 0.5),
         )
         return circuit, expected
@@ -270,7 +270,7 @@ def case_rotation_gates_equivalent_to_pauli_evolution(self):
                 ExpGate(PauliOperator({Pauli("Z", 0): 0.5}, num_qubits=1)),
             ),
             initial_state=(),
-            num_modes=1,
+            system_size=1,
             parameters=(0.5, 0.3, 0.7),
         )
         return circuit, expected
@@ -283,7 +283,7 @@ def case_barrier_ignored(self):
         expected = Circuit(
             gates=(ExpGate(PauliOperator({Pauli("Z", 0): 1.0}, num_qubits=1)),),
             initial_state=(),
-            num_modes=1,
+            system_size=1,
             parameters=(0.7,),
         )
         return circuit, expected
@@ -324,7 +324,7 @@ def test_to_qiskit_circuit_rejects_unbound() -> None:
     circuit = Circuit(
         gates=(ExpGate(PauliOperator({Pauli("X", 0): 1.0}, num_qubits=1)),),
         initial_state=(),
-        num_modes=1,
+        system_size=1,
     )  # no parameter values
     with pytest.raises(ValueError, match="bound circuit"):
         to_qiskit_circuit(circuit, num_qubits=1)
@@ -342,7 +342,7 @@ def test_from_to_qiskit_circuit_roundtrip() -> None:
     circuit = Circuit(
         gates=(ExpGate(PauliOperator({Pauli("XYZ", (3, 1, 2)): 1.0}, num_qubits=4)),),
         initial_state=[],
-        num_modes=4,
+        system_size=4,
         parameters=[-1.2],
     )  # no parameter values
     qcirc = to_qiskit_circuit(circuit, num_qubits=4)
diff --git a/tests/test_update_methods.py b/tests/test_update_methods.py
index e3711f9c..cbd9056e 100644
--- a/tests/test_update_methods.py
+++ b/tests/test_update_methods.py
@@ -107,7 +107,7 @@ def test_update_cutoff_valid(self, mp):
         mp.cutoff = 6
         # A weight-6 monomial takes an imaginary Hermitian coefficient (like weight-2).
         gate = ExpGate(MajoranaOperator({(0, 1, 2, 3, 4, 5): 1.0j}, num_modes=4))
-        mp.build_graph(Circuit((gate,), initial_state=(), num_modes=4))
+        mp.build_graph(Circuit((gate,), initial_state=(), system_size=4))
         assert mp.size() > 0
 
     def test_update_cutoff_invalid(self, mp):
@@ -144,7 +144,7 @@ def test_integration(self, serial_comm):
             majoranas=[(0, 2), (1, 3)],
             gen_coeffs=[0.0, 0.0],
             param_inds=[0, 1],
-            num_modes=4,
+            system_size=4,
             parameters=[1.0, 1.0],
         )
         mp = MajoranaPropagator(

From 28b64e1efe8e8e6cfe1650224e73d0e4de646d98 Mon Sep 17 00:00:00 2001
From: Adam <aglos@iitis.pl>
Date: Mon, 27 Jul 2026 21:31:19 +0000
Subject: [PATCH 4/6] feat: :sparkles: expgate has system size method now

---
 src/monoprop/circuit.py | 15 +++++++++++++++
 tests/test_circuit.py   | 26 ++++++++++++++++++++++++++
 2 files changed, 41 insertions(+)

diff --git a/src/monoprop/circuit.py b/src/monoprop/circuit.py
index 23a0ce0c..266c69d4 100644
--- a/src/monoprop/circuit.py
+++ b/src/monoprop/circuit.py
@@ -224,6 +224,21 @@ def __repr__(self) -> str:
         """Return a string representation such as ``ExpGate(<generator>, index=0)``."""
         return f"{self.__class__.__name__}({self.generator}, index={self.index})"
 
+    @property
+    def system_size(self) -> int:
+        """Number of modes/qubits the generator acts on.
+
+        Reads ``num_modes`` off a [MajoranaOperator][monoprop.majorana.MajoranaOperator]
+        generator, or ``num_qubits`` off a [PauliOperator][monoprop.pauli.PauliOperator]
+        generator.
+
+        Returns:
+            Number of modes/qubits
+        """
+        if isinstance(self.generator, PauliOperator):
+            return self.generator.num_qubits
+        return self.generator.num_modes
+
 
 class Circuit:
     """A variational circuit: an ordered sequence of exponential gates, angles, and a state.
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index 4d57b94b..70a0cea5 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -88,6 +88,32 @@ def test_exp_from_fermi_generator_becomes_majorana() -> None:
     assert isinstance(gate.generator, MajoranaOperator)
 
 
+@pytest.mark.parametrize(
+    ("generator", "expected"),
+    [
+        pytest.param(
+            MajoranaOperator({(0, 1): 1.0j}, num_modes=3), 3, id="majorana_operator"
+        ),
+        pytest.param(
+            PauliOperator({Pauli("XX"): 1.0}, num_qubits=2), 2, id="pauli_operator"
+        ),
+        pytest.param(
+            # A FermiOperator generator is converted to Majorana; system_size follows suit.
+            FermiOperator(
+                [[(0, "+"), (1, "-")], [(1, "+"), (0, "-")]], [1.0, -1.0], num_modes=4
+            ),
+            4,
+            id="fermi_operator",
+        ),
+    ],
+)
+def test_exp_gate_system_size_reads_from_generator(
+    generator: MajoranaOperator | PauliOperator | FermiOperator, expected: int
+) -> None:
+    """system_size reads num_modes off a Majorana generator, num_qubits off a Pauli one."""
+    assert ExpGate(generator).system_size == expected
+
+
 class ExpGateAtolCases:
     @case(id="single_excitation")
     def case_single_excitation(self):

From a7734896b310baedf3e2de591fe45e82035479de Mon Sep 17 00:00:00 2001
From: Adam <aglos@iitis.pl>
Date: Mon, 27 Jul 2026 21:32:56 +0000
Subject: [PATCH 5/6] refactor: :recycle: expgate verification in circuit
 simplified

---
 src/monoprop/circuit.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/monoprop/circuit.py b/src/monoprop/circuit.py
index 266c69d4..4e08a921 100644
--- a/src/monoprop/circuit.py
+++ b/src/monoprop/circuit.py
@@ -316,11 +316,7 @@ def __init__(  # noqa:C901
                 raise TypeError(
                     f"Circuit gates must be ExpGate; got {type(gate).__name__}."
                 )
-            gate_size = (
-                gate.generator.num_qubits
-                if isinstance(gate.generator, PauliOperator)
-                else gate.generator.num_modes
-            )
+            gate_size = gate.system_size
             if gate_size != system_size:
                 raise ValueError(
                     f"Gate generator width {gate_size} does not match circuit system_size={system_size}."

From 4ddcc1ad77e5861050d3d14dcfe766d66eb073b0 Mon Sep 17 00:00:00 2001
From: Adam <aglos@iitis.pl>
Date: Mon, 27 Jul 2026 21:50:04 +0000
Subject: [PATCH 6/6] refactor: :recycle: propagators refactoring

---
 src/monoprop/monomial_propagator.py | 23 ++++++++---------------
 src/monoprop/pauli_propagator.py    | 11 ++---------
 tests/test_only_rotate_k.py         |  8 +++++++-
 3 files changed, 17 insertions(+), 25 deletions(-)

diff --git a/src/monoprop/monomial_propagator.py b/src/monoprop/monomial_propagator.py
index 44c6c405..5222924c 100644
--- a/src/monoprop/monomial_propagator.py
+++ b/src/monoprop/monomial_propagator.py
@@ -78,7 +78,7 @@ class MonomialPropagator(ABC):
 
     _comm: MPI.Comm | None
     _n_params: int
-    _num_modes: int
+    _system_size: int
     _initial_state: list[int]
     _simulator: object
 
@@ -115,10 +115,7 @@ def _init_simulator(
 
         self._comm = comm
         self._n_params = 0
-        self._num_modes = num_modes
-        # System qubit count for expanding Pauli gates; set by PauliPropagator from the
-        # observable. None for a native Majorana propagator (its gates need no qubit count).
-        self._num_qubits = None
+        self._system_size = num_modes
         self._initial_state = list(initial_state)
         # dispatch() is typed to return the base `type[_SimulatorAdapter]`, whose __init__ takes
         # extra positional args the generated per-mode subclasses fill in; the kwargs below match
@@ -198,11 +195,10 @@ def _check_initial_state(self, circuit: Circuit) -> None:
 
     def _check_circuit_width(self, circuit: Circuit) -> None:
         """Reject a circuit with a system width that disagrees with the propagator."""
-        expected = self._num_qubits if self._num_qubits is not None else self._num_modes
-        if circuit.system_size != expected:
+        if circuit.system_size != self._system_size:
             raise ValueError(
                 f"Circuit system_size={circuit.system_size} does not match propagator width "
-                f"{expected}."
+                f"{self._system_size}."
             )
 
     def _validate_and_correct_only_rotate_len_k(
@@ -224,10 +220,7 @@ def _validate_and_correct_only_rotate_len_k(
         """
         if only_rotate_len_k is None:
             return 0
-        if only_rotate_len_k <= 0 or (
-            isinstance(self._num_qubits, int)
-            and only_rotate_len_k > 2 * self._num_qubits
-        ):
+        if only_rotate_len_k <= 0 or only_rotate_len_k > 2 * self._system_size:
             raise ValueError(
                 f"only_rotate_len_k={only_rotate_len_k} is out of range; must be 0 < k <= 2*num_qubits "
             )
@@ -286,7 +279,7 @@ def build_graph(
         else:
             seed = None
         gates = self._circuit_gates(circuit)
-        num_qubits = 0 if self._num_qubits is None else self._num_qubits
+        num_qubits = self._system_size
         # Shift the circuit's local 0-based angle indices onto the accumulated axis.
         mapping = [self._n_params + m for m in circuit.resolved_mapping]
         self._n_params += circuit.n_parameters
@@ -325,7 +318,7 @@ def propagate(
         self._check_initial_state(circuit)
         self._check_circuit_width(circuit)
         gates = self._circuit_gates(circuit)
-        num_qubits = 0 if self._num_qubits is None else self._num_qubits
+        num_qubits = self._system_size
         majoranas, gen_coeffs, mapping, _gate_indices = expand_monomials(
             gates, circuit.resolved_mapping, num_qubits
         )
@@ -597,7 +590,7 @@ def graph_size(self) -> tuple[int, int]:
     @property
     def num_modes(self) -> int:
         """Number of fermionic modes for the simulator."""
-        return self._num_modes
+        return self._system_size
 
     @property
     def graph_layers(self) -> int:
diff --git a/src/monoprop/pauli_propagator.py b/src/monoprop/pauli_propagator.py
index 1da8d11d..e7b25f85 100644
--- a/src/monoprop/pauli_propagator.py
+++ b/src/monoprop/pauli_propagator.py
@@ -102,18 +102,11 @@ def __init__(
             basis_change=jordan_wigner_basis_change(num_qubits),
             comm=comm,
         )
-        # The qubit count comes from the observable and is carried into Pauli gate expansion
-        # via build_graph (_init_simulator initializes it to None).
-        self._num_qubits = num_qubits
 
     @property
     def num_qubits(self) -> int:
         """Number of qubits the propagator acts on."""
-        # Always set in __init__ (which raises if the observable has no qubit count); the base
-        # declares it Optional for the native Majorana propagator.
-        if self._num_qubits is None:
-            raise RuntimeError("PauliPropagator has no qubit count set.")
-        return self._num_qubits
+        return self._system_size
 
     def _circuit_gates(self, circuit: Circuit) -> Sequence[ExpGate]:
         """Accept a qubit circuit; its gates are expanded by the shared pipeline.
@@ -121,7 +114,7 @@ def _circuit_gates(self, circuit: Circuit) -> Sequence[ExpGate]:
         A ``PauliPropagator`` rejects a Majorana/fermionic circuit. The Jordan-Wigner mapping
         and antihermitian normalization live in [expand_monomials][monoprop.circuit.expand_monomials];
         the propagator's ``num_qubits`` (from the observable) reaches the expander via
-        ``self._num_qubits``.
+        ``self._system_size``.
         """
         if circuit.family == "majorana":
             raise TypeError(
diff --git a/tests/test_only_rotate_k.py b/tests/test_only_rotate_k.py
index 8a7109fb..1c0ac309 100644
--- a/tests/test_only_rotate_k.py
+++ b/tests/test_only_rotate_k.py
@@ -148,7 +148,13 @@ def test_only_rotate_len_k(problem, inplace, serial_mp_kwargs):
                 match=r"only_rotate_len_k=0 is out of range; must be 0 < k <= 2\*num_qubits",
             ),
         ),
-        (9, does_not_raise()),
+        (
+            9,
+            pytest.raises(
+                ValueError,
+                match=r"only_rotate_len_k=9 is out of range; must be 0 < k <= 2\*num_qubits",
+            ),
+        ),
         (8, does_not_raise()),
     ],
 )