From f873d594e799254a648d85fea281922eb8d471da Mon Sep 17 00:00:00 2001
From: linusschulte <linus.schulte@gmx.com>
Date: Mon, 6 Jul 2026 19:38:15 +0200
Subject: [PATCH 1/9] Added non-qubit local observables to facilitate
 trapped-ion position measurement

---
 CHANGELOG.md                                  |  1 +
 docs/examples/trapped_ion.md                  | 56 +++++++++--------
 src/mqt/yaqs/core/data_structures/mps.py      | 62 ++++++++++++-------
 .../data_structures/simulation_parameters.py  | 30 +++++----
 src/mqt/yaqs/core/libraries/gate_library.py   | 33 ++++++++++
 tests/core/data_structures/test_mpo.py        | 41 ++++++++++++
 tests/core/data_structures/test_mps.py        | 33 ++++++++++
 .../test_simulation_parameters.py             | 19 ++++++
 8 files changed, 214 insertions(+), 61 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 158313bf9..22a119a8e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,7 @@ This project adheres to [Semantic Versioning], with the exception that minor rel
 
 ### Added
 
+- added custom one-site matrix observables for non-qubit local dimensions
 - added analytical noise characterization module and digital twin pipeline ([#288]) ([**@aleramos119**], [**@aaronleesander**])
 - added memory characterization module for analyzing non-Markovian processes ([#482]) ([**@aaronleesander**])
 - added an MPO constructor for static one- and two-ion trapped-ion Hamiltonians in the position basis ([#476]) ([**@linusschulte**])
diff --git a/docs/examples/trapped_ion.md b/docs/examples/trapped_ion.md
index caf376b7b..f1de7bc7f 100644
--- a/docs/examples/trapped_ion.md
+++ b/docs/examples/trapped_ion.md
@@ -29,7 +29,7 @@ so after half a trap period it reaches the opposite turning point.
 ```{code-cell} ipython3
 import numpy as np
 
-from mqt.yaqs import Hamiltonian, MPO, State
+from mqt.yaqs import Hamiltonian, MPO, Observable, State
 
 omega = 1.0
 initial_displacement = 1.0
@@ -43,6 +43,7 @@ initial_grid_state /= np.linalg.norm(initial_grid_state)
 
 hamiltonian = Hamiltonian.from_mpo(MPO.trapped_ion(positions, masses=[1.0], omega=omega))
 state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim])
+position_observable = Observable(np.diag(positions), 0)
 ```
 
 ## 2. Noiseless evolution to $T/2$
@@ -51,7 +52,7 @@ state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim
 from mqt.yaqs import AnalogSimParams, Simulator
 
 params = AnalogSimParams(
-    observables=[],
+    observables=[position_observable],
     elapsed_time=half_period,
     dt=half_period / 16,
     max_bond_dim=None,
@@ -59,16 +60,18 @@ params = AnalogSimParams(
     krylov_tol=1e-12,
     preset="exact",
     get_state=True,
-    sample_timesteps=False,
+    sample_timesteps=True,
 )
 
 result = Simulator(show_progress=False).run(state, hamiltonian, params)
 final_state = result.output_state.vector
-final_x = float(np.sum(positions * np.abs(final_state) ** 2))
+position_expectation = np.real(result.expectation_values[0])
+final_x = float(position_expectation[-1])
 ```
 
-The final $\langle x\rangle$ is close to $-x_0$ but not exact because the simulation uses a
-finite grid and a finite-difference kinetic operator.
+The position observable is a custom one-site matrix on the grid basis. The final
+$\langle x\rangle$ is close to $-x_0$ but not exact because the simulation uses a finite
+grid and a finite-difference kinetic operator.
 
 ```{code-cell} ipython3
 print(f"Initial <x>       = {initial_displacement:.6f}")
@@ -76,7 +79,7 @@ print(f"Final <x> at T/2  = {final_x:.6f}")
 print(f"Continuum target  = {-initial_displacement:.6f}")
 ```
 
-## 3. Wavepacket at $t=0$ and $t=T/2$
+## 3. Wavepacket over time
 
 ```{code-cell} ipython3
 ---
@@ -87,24 +90,27 @@ mystnb:
 ---
 import matplotlib.pyplot as plt
 
-prob_initial = np.abs(initial_grid_state) ** 2
-prob_final = np.abs(final_state) ** 2
-
-fig, axes = plt.subplots(1, 2, figsize=(8, 3.2), layout="constrained", sharey=True)
-axes[0].fill_between(positions, prob_initial, alpha=0.35, color="tab:blue")
-axes[0].plot(positions, prob_initial, color="tab:blue", lw=1.5)
-axes[0].set_title(r"$t = 0$")
-axes[0].set_xlabel(r"$x$")
-axes[0].set_ylabel(r"$|\psi(x)|^2$")
-axes[0].grid(alpha=0.3)
-
-axes[1].fill_between(positions, prob_final, alpha=0.35, color="tab:orange")
-axes[1].plot(positions, prob_final, color="tab:orange", lw=1.5)
-axes[1].set_title(rf"$t = T/2$")
-axes[1].set_xlabel(r"$x$")
-axes[1].grid(alpha=0.3)
-
-fig.suptitle("Harmonic wavepacket reflection on a position grid")
+dense_hamiltonian = hamiltonian.to_matrix()
+eigenvalues, eigenvectors = np.linalg.eigh(dense_hamiltonian)
+coefficients = eigenvectors.conj().T @ initial_grid_state
+phases = np.exp(-1j * eigenvalues[:, None] * params.times[None, :])
+states = eigenvectors @ (coefficients[:, None] * phases)
+probability_density = np.abs(states) ** 2
+
+fig, ax = plt.subplots(figsize=(7.2, 3.6), layout="constrained")
+image = ax.imshow(
+    probability_density,
+    aspect="auto",
+    origin="lower",
+    extent=(params.times[0], params.times[-1], positions[0], positions[-1]),
+    cmap="viridis",
+)
+ax.plot(params.times, position_expectation, color="white", lw=1.4, label=r"$\langle x\rangle$")
+ax.set_xlabel(r"$t$")
+ax.set_ylabel(r"$x$")
+ax.set_title("Position-grid wavepacket density")
+ax.legend(loc="upper right")
+fig.colorbar(image, ax=ax, label=r"$|\psi(x,t)|^2$")
 plt.show()
 ```
 
diff --git a/src/mqt/yaqs/core/data_structures/mps.py b/src/mqt/yaqs/core/data_structures/mps.py
index bf413bf37..50be8bf11 100644
--- a/src/mqt/yaqs/core/data_structures/mps.py
+++ b/src/mqt/yaqs/core/data_structures/mps.py
@@ -973,32 +973,34 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
         Returns:
             np.complex128: The computed expectation value (typically, its real part is of interest).
 
+        Raises:
+            ValueError: If the observable is not supported or its matrix shape does not match the target site.
+
         Notes:
             A deep copy of the state is used to prevent modifications to the original MPS.
             Requires :meth:`check_covers_sites` to hold for ``sites``; prefer :meth:`expect` for
             gauge-safe evaluation.
         """
         temp_state = copy.deepcopy(self)
-        if operator.gate.matrix.shape[0] == 2:  # Local observable
-            i = None
-            if isinstance(sites, list):
-                i = sites[0]
-            elif isinstance(sites, int):
-                i = sites
-
-            if isinstance(operator.sites, list):
-                assert operator.sites[0] == i, f"Operator sites mismatch {operator.sites[0]}, {i}"
-            elif isinstance(operator.sites, int):
-                assert operator.sites == i, f"Operator sites mismatch {operator.sites}, {i}"
+        sites_list = [sites] if isinstance(sites, int) else list(sites)
+        operator_sites = [operator.sites] if isinstance(operator.sites, int) else list(operator.sites)
 
-            assert i is not None, f"Invalid type for 'sites': expected int or list[int], got {type(sites).__name__}"
+        if operator.gate.interaction == 1:
+            assert len(sites_list) == 1, f"One-site observable requires one site, got {sites_list}."
+            assert operator_sites == sites_list, f"Operator sites mismatch {operator_sites}, {sites_list}"
+            i = sites_list[0]
             a = temp_state.tensors[i]
-            temp_state.tensors[i] = oe.contract("ab, bcd->acd", operator.gate.matrix, a)
+            local_dim = a.shape[0]
+            matrix = np.asarray(operator.gate.matrix, dtype=np.complex128)
+            if matrix.shape != (local_dim, local_dim):
+                msg = f"Local observable matrix shape {matrix.shape} does not match site {i} dimension {local_dim}."
+                raise ValueError(msg)
+            temp_state.tensors[i] = oe.contract("ab, bcd->acd", matrix, a)
 
-        elif operator.gate.matrix.shape[0] == 4:  # Two-site correlator
+        elif operator.gate.interaction == 2:
             assert isinstance(sites, list)
             assert isinstance(operator.sites, list)
-            i, j = sites
+            i, j = sites_list
 
             assert operator.sites[0] == i, "Observable sites mismatch"
             assert operator.sites[1] == j, "Observable sites mismatch"
@@ -1010,6 +1012,10 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
             b = temp_state.tensors[j]
             d_i, left, _ = a.shape
             d_j, _, right = b.shape
+            matrix = np.asarray(operator.gate.matrix, dtype=np.complex128)
+            if matrix.shape != (d_i * d_j, d_i * d_j):
+                msg = f"Two-site observable matrix shape {matrix.shape} does not match site dimensions {d_i} and {d_j}."
+                raise ValueError(msg)
 
             # 1) merge A,B into theta of shape (l, d_i*d_j, r)
             theta = np.tensordot(a, b, axes=(2, 1))  # (d_i, l, d_j, r)
@@ -1017,7 +1023,7 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
             theta = theta.reshape(left, d_i * d_j, right)  # (l, d_i*d_j, r)
 
             # 2) apply operator on the combined phys index
-            theta = oe.contract("ab, cbd->cad", operator.gate.matrix, theta)  # (l, d_i*d_j, r)
+            theta = oe.contract("ab, cbd->cad", matrix, theta)  # (l, d_i*d_j, r)
             theta = theta.reshape(left, d_i, d_j, right)  # back to (l, d_i, d_j, r)
 
             # 3) split via SVD
@@ -1035,6 +1041,9 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
 
             temp_state.tensors[i] = a_new
             temp_state.tensors[j] = b_new
+        else:
+            msg = "Local observable must be one-site or nearest-neighbor two-site."
+            raise ValueError(msg)
 
         return self.scalar_product(temp_state, sites)
 
@@ -1099,26 +1108,33 @@ def bubble_swaps_backward(state: MPS) -> None:
             for i in reversed(range(state.length - 2)):
                 apply_two_site_nn_inplace(state, i, sw)
 
-        sites = [observable.sites] if isinstance(observable.sites, int) else observable.sites
+        sites = [observable.sites] if isinstance(observable.sites, int) else list(observable.sites)
 
-        if observable.gate.matrix.shape[0] == 2:
+        if observable.gate.interaction == 1:
+            if len(sites) != 1:
+                msg = f"One-site local observable requires one site, got {sites}."
+                raise ValueError(msg)
             site = sites[0]
-            self.tensors[site] = oe.contract("ab, bcd->acd", observable.gate.matrix, self.tensors[site])
+            local_dim = self.tensors[site].shape[0]
+            matrix = np.asarray(observable.gate.matrix, dtype=np.complex128)
+            if matrix.shape != (local_dim, local_dim):
+                msg = f"Local observable matrix shape {matrix.shape} does not match site {site} dimension {local_dim}."
+                raise ValueError(msg)
+            self.tensors[site] = oe.contract("ab, bcd->acd", matrix, self.tensors[site])
             return
 
-        if observable.gate.matrix.shape[0] == 4:
+        if observable.gate.interaction == 2:
             i, j = int(sites[0]), int(sites[1])
             length = self.length
+            mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
 
             if length == 2:
                 if i == length - 1 and j == 0:
-                    mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
                     g_merged = permuted_periodic_wrap(mat)
                     apply_two_site_nn_inplace(self, 0, g_merged)
                     return
                 i, j = min(i, j), max(i, j)
             elif (i == length - 1 and j == 0) or (i == 0 and j == length - 1):
-                mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
                 bubble_swaps_forward(self)
                 g_merged = permuted_periodic_wrap(mat)
                 apply_two_site_nn_inplace(self, length - 2, g_merged)
@@ -1129,7 +1145,7 @@ def bubble_swaps_backward(state: MPS) -> None:
                 msg = "Only nearest-neighbor two-site observables are currently implemented."
                 raise ValueError(msg)
 
-            apply_two_site_nn_inplace(self, i, np.asarray(observable.gate.matrix, dtype=np.complex128))
+            apply_two_site_nn_inplace(self, i, mat)
             return
 
         msg = "Local observable must be one-site or nearest-neighbor two-site."
diff --git a/src/mqt/yaqs/core/data_structures/simulation_parameters.py b/src/mqt/yaqs/core/data_structures/simulation_parameters.py
index daf891b69..90f8b2540 100644
--- a/src/mqt/yaqs/core/data_structures/simulation_parameters.py
+++ b/src/mqt/yaqs/core/data_structures/simulation_parameters.py
@@ -23,10 +23,10 @@
 
 import numpy as np
 
-from mqt.yaqs.core.libraries.gate_library import GateLibrary
+from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary
 
 if TYPE_CHECKING:
-    from mqt.yaqs.core.libraries.gate_library import BaseGate
+    from numpy.typing import ArrayLike
 
 SimulationPreset = Literal["fast", "balanced", "accurate", "exact"]
 GateMode = Literal["tdvp", "full-tdvp", "swaps", "mpo"]
@@ -233,31 +233,35 @@ class Observable:
         sites: The site or site indices on which this observable is measured.
     """
 
-    def __init__(self, gate: BaseGate | str, sites: int | list[int] | None = None) -> None:
+    def __init__(self, gate: BaseGate | str | ArrayLike, sites: int | list[int] | None = None) -> None:
         """Initializes an Observable instance.
 
         Args:
-            gate: The gate that will act as the observable.
+            gate: The gate or one-site local matrix that will act as the observable.
             sites: The qubit or site indices on which this observable is measured.
         """
         if isinstance(gate, str):
             if gate == "entropy":
-                gate = GateLibrary.entropy()
+                resolved_gate = GateLibrary.entropy()
             elif gate == "schmidt_spectrum":
-                gate = GateLibrary.schmidt_spectrum()
+                resolved_gate = GateLibrary.schmidt_spectrum()
             elif gate == "pvm":
-                gate = GateLibrary.pvm(gate)
+                resolved_gate = GateLibrary.pvm(gate)
             elif hasattr(GateLibrary, gate):
                 attr = getattr(GateLibrary, gate)
                 try:
-                    gate = attr()
+                    resolved_gate = attr()
                 except TypeError:
-                    gate = GateLibrary.pvm(gate)
+                    resolved_gate = GateLibrary.pvm(gate)
             else:
-                gate = GateLibrary.pvm(gate)
-        assert hasattr(GateLibrary, gate.name), f"Observable {gate.name} not found in GateLibrary."
-        self.gate = copy.deepcopy(gate)
-        if gate.name != "pvm":
+                resolved_gate = GateLibrary.pvm(gate)
+        elif isinstance(gate, BaseGate):
+            resolved_gate = gate
+        else:
+            resolved_gate = GateLibrary.local(gate)
+        assert hasattr(GateLibrary, resolved_gate.name), f"Observable {resolved_gate.name} not found in GateLibrary."
+        self.gate: BaseGate = copy.deepcopy(resolved_gate)
+        if resolved_gate.name != "pvm":
             assert sites is not None
             self.sites = sites
             self.gate.set_sites(self.sites)
diff --git a/src/mqt/yaqs/core/libraries/gate_library.py b/src/mqt/yaqs/core/libraries/gate_library.py
index 0cba89819..20a8a7274 100644
--- a/src/mqt/yaqs/core/libraries/gate_library.py
+++ b/src/mqt/yaqs/core/libraries/gate_library.py
@@ -1650,6 +1650,37 @@ def __init__(self, bitstring: str) -> None:
         super().__init__(mat)
 
 
+class LocalOperator(BaseGate):
+    """Custom one-site operator for arbitrary local Hilbert-space dimensions.
+
+    This gate is intended for observables such as position-grid operators on
+    qudits or oscillator truncations. Unlike :class:`BaseGate`, it does not
+    interpret the matrix dimension as a qubit interaction count.
+    """
+
+    name = "local"
+
+    def __init__(self, matrix: ArrayLike) -> None:
+        """Create a one-site local operator.
+
+        Args:
+            matrix: Square matrix acting on one local site.
+
+        Raises:
+            ValueError: If ``matrix`` is not a square two-dimensional array.
+        """
+        mat = np.asarray(matrix, dtype=np.complex128)
+        if mat.ndim != 2:
+            msg = "Local operator matrix must be a 2-D array."
+            raise ValueError(msg)
+        if mat.shape[0] != mat.shape[1]:
+            msg = "Local operator matrix must be square."
+            raise ValueError(msg)
+        self.matrix = mat
+        self.tensor = mat
+        self.interaction = 1
+
+
 class Entropy(BaseGate):
     """Meta-observable for bipartite entanglement entropy across a cut.
 
@@ -1761,6 +1792,7 @@ class GateLibrary:
         p0: Class for projector ``|0⟩⟨0|``.
         p1: Class for projector ``|1⟩⟨1|``.
         pvm: Class for projection-valued measurement onto a given bitstring.
+        local: Class for arbitrary one-site local operators.
 
         entropy:      Class representing a request for bipartite entanglement entropy across a cut.
         schmidt_spectrum: Class representing a request for the Schmidt spectrum across a cut.
@@ -1811,6 +1843,7 @@ class GateLibrary:
     p0 = P0
     p1 = P1
     pvm = PVM
+    local = LocalOperator
 
     entropy = Entropy
     schmidt_spectrum = SchmidtSpectrum
diff --git a/tests/core/data_structures/test_mpo.py b/tests/core/data_structures/test_mpo.py
index 1b6b33f9e..19d0c1499 100644
--- a/tests/core/data_structures/test_mpo.py
+++ b/tests/core/data_structures/test_mpo.py
@@ -513,6 +513,47 @@ def test_trapped_ion_coulomb_truncation() -> None:
     np.testing.assert_allclose(truncated_coulomb, expected_rank_2, atol=1e-12)
 
 
+def test_trapped_ion_one_ion_position_observable_centers_on_trap() -> None:
+    """The one-ion ground-state position expectation follows the static trap center."""
+    positions = np.linspace(-2.0, 2.0, 9, dtype=np.float64)
+    trap_center = 0.4
+    mpo = MPO.trapped_ion(positions, [1.0], omega=1.0, trap_center=trap_center)
+    _energy, eigenvectors = np.linalg.eigh(mpo.to_matrix())
+    ground_state = eigenvectors[:, 0]
+    mps = MPS(
+        length=1,
+        tensors=[ground_state.reshape(positions.size, 1, 1)],
+        physical_dimensions=[positions.size],
+    )
+
+    position = Observable(np.diag(positions), 0)
+
+    np.testing.assert_allclose(mps.expect(position), trap_center, atol=6e-2)
+
+
+def test_trapped_ion_two_ion_coulomb_increases_ground_state_separation() -> None:
+    """Softened Coulomb repulsion increases the two-ion ground-state separation."""
+    positions = np.linspace(-3.0, 3.0, 9, dtype=np.float64)
+    grid_dim = positions.size
+    separation = np.abs(positions[:, None] - positions[None, :])
+
+    def ground_state_separation(coulomb_strength: float) -> float:
+        mpo = MPO.trapped_ion(
+            positions,
+            [1.0, 1.0],
+            omega=0.6,
+            coulomb_strength=coulomb_strength,
+        )
+        _energy, eigenvectors = np.linalg.eigh(mpo.to_matrix())
+        probabilities = np.abs(eigenvectors[:, 0].reshape(grid_dim, grid_dim)) ** 2
+        return float(np.sum(separation * probabilities))
+
+    uncoupled_separation = ground_state_separation(0.0)
+    repulsive_separation = ground_state_separation(0.8)
+
+    assert repulsive_separation > uncoupled_separation + 0.2
+
+
 @pytest.mark.parametrize(
     ("kwargs", "match"),
     [
diff --git a/tests/core/data_structures/test_mps.py b/tests/core/data_structures/test_mps.py
index 2949f4109..f6728df7e 100644
--- a/tests/core/data_structures/test_mps.py
+++ b/tests/core/data_structures/test_mps.py
@@ -441,6 +441,39 @@ def test_local_expect_x_on_plus_state() -> None:
     np.testing.assert_allclose(val, 1.0, atol=1e-12)
 
 
+def test_non_qubit_local_expectation_from_matrix_observable() -> None:
+    """A matrix observable can be measured on a non-qubit local site."""
+    amplitudes = np.sqrt(np.array([0.2, 0.3, 0.5], dtype=np.float64)).astype(np.complex128)
+    tensor0 = amplitudes.reshape(3, 1, 1)
+    tensor1 = np.array([1.0, 0.0], dtype=np.complex128).reshape(2, 1, 1)
+    psi_mps = MPS(length=2, tensors=[tensor0, tensor1], physical_dimensions=[3, 2])
+    position = np.diag(np.array([-1.0, 0.5, 2.0], dtype=np.float64))
+
+    val = psi_mps.expect(Observable(position, 0))
+
+    expected = 0.2 * -1.0 + 0.3 * 0.5 + 0.5 * 2.0
+    np.testing.assert_allclose(val, expected, atol=1e-12)
+
+
+def test_four_level_local_observable_is_not_treated_as_two_site() -> None:
+    """A ``4 x 4`` matrix observable acts on one four-level site when wrapped as local."""
+    amplitudes = np.array([0.5, 0.5j, -0.5, 0.5], dtype=np.complex128)
+    psi_mps = MPS(length=1, tensors=[amplitudes.reshape(4, 1, 1)], physical_dimensions=[4])
+    matrix = np.diag(np.array([0.0, 1.0, 2.0, 3.0], dtype=np.float64))
+
+    val = psi_mps.expect(Observable(matrix, 0))
+
+    np.testing.assert_allclose(val, 1.5, atol=1e-12)
+
+
+def test_local_observable_dimension_mismatch_raises() -> None:
+    """One-site matrix observables must match the measured site's local dimension."""
+    psi_mps = MPS(length=1, physical_dimensions=[3], state="zeros")
+
+    with pytest.raises(ValueError, match="does not match site 0 dimension 3"):
+        psi_mps.expect(Observable(np.eye(2), 0))
+
+
 def test_mps_apply_local_l2_periodic_wrap_matches_permuted_nn() -> None:
     """For ``L == 2``, wrap-ordered and permuted NN applications must agree."""
     length = 2
diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py
index d5472d821..e4a4d6d7d 100644
--- a/tests/core/data_structures/test_simulation_parameters.py
+++ b/tests/core/data_structures/test_simulation_parameters.py
@@ -62,6 +62,25 @@ def test_observable_creation_valid() -> None:
     assert obs.sites == site
 
 
+def test_observable_accepts_custom_local_matrix() -> None:
+    """Observable accepts square matrices as arbitrary one-site local operators."""
+    matrix = np.diag(np.array([-1.0, 0.25, 2.0]))
+
+    obs = Observable(matrix, 0)
+
+    assert obs.gate.name == "local"
+    assert obs.gate.interaction == 1
+    np.testing.assert_allclose(obs.gate.matrix, matrix)
+    assert obs.sites == 0
+
+
+@pytest.mark.parametrize("matrix", [np.ones(3), np.ones((2, 3))])
+def test_observable_rejects_invalid_custom_local_matrix(matrix: np.ndarray) -> None:
+    """Matrix observables must be two-dimensional and square."""
+    with pytest.raises(ValueError, match="Local operator matrix"):
+        Observable(matrix, 0)
+
+
 def test_analog_simparams_basic() -> None:
     """Test that AnalogSimParams is initialized with correct parameters.
 

From 70009dfeec9d252a724f630304c4543afccfaa73 Mon Sep 17 00:00:00 2001
From: linusschulte <119408387+linusschulte@users.noreply.github.com>
Date: Wed, 8 Jul 2026 09:01:55 +0200
Subject: [PATCH 2/9] adjusted changelog

Signed-off-by: linusschulte <119408387+linusschulte@users.noreply.github.com>
---
 CHANGELOG.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 22a119a8e..dd83dd4fb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,7 +13,7 @@ This project adheres to [Semantic Versioning], with the exception that minor rel
 
 ### Added
 
-- added custom one-site matrix observables for non-qubit local dimensions
+- added support for custom non-qubit local observables ([#497]) ([**@linusschulte**])
 - added analytical noise characterization module and digital twin pipeline ([#288]) ([**@aleramos119**], [**@aaronleesander**])
 - added memory characterization module for analyzing non-Markovian processes ([#482]) ([**@aaronleesander**])
 - added an MPO constructor for static one- and two-ion trapped-ion Hamiltonians in the position basis ([#476]) ([**@linusschulte**])
@@ -151,6 +151,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool
 
 <!-- PR links -->
 
+[#497]: https://github.com/munich-quantum-toolkit/yaqs/pull/497
 [#482]: https://github.com/munich-quantum-toolkit/yaqs/pull/482
 [#481]: https://github.com/munich-quantum-toolkit/yaqs/pull/481
 [#476]: https://github.com/munich-quantum-toolkit/yaqs/pull/476

From efd1d999019a315ec99b5137851b03b7815ab9b5 Mon Sep 17 00:00:00 2001
From: linusschulte <linus.schulte@gmx.com>
Date: Wed, 8 Jul 2026 09:39:26 +0200
Subject: [PATCH 3/9] adjusted changelog

---
 CHANGELOG.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index dd83dd4fb..87ee0268d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,10 @@ This project adheres to [Semantic Versioning], with the exception that minor rel
 
 ## [Unreleased]
 
+### Added
+
+- added custom one-site matrix observables for non-qubit local dimensions ([#497]) ([**@linusschulte**])
+
 ## [0.6.0] - 2026-07-06
 
 ### Added

From d38154d379f13c9463175ea81bb7d6122b2093de Mon Sep 17 00:00:00 2001
From: linusschulte <linus.schulte@gmx.com>
Date: Wed, 29 Jul 2026 16:41:44 +0200
Subject: [PATCH 4/9] added tests to satisfy code coverage

---
 tests/core/data_structures/test_mps.py        | 28 +++++++++++++++++++
 .../test_simulation_parameters.py             | 10 +++++++
 2 files changed, 38 insertions(+)

diff --git a/tests/core/data_structures/test_mps.py b/tests/core/data_structures/test_mps.py
index f6728df7e..40b27d8d1 100644
--- a/tests/core/data_structures/test_mps.py
+++ b/tests/core/data_structures/test_mps.py
@@ -474,6 +474,34 @@ def test_local_observable_dimension_mismatch_raises() -> None:
         psi_mps.expect(Observable(np.eye(2), 0))
 
 
+def test_two_site_local_observable_dimension_mismatch_raises() -> None:
+    """Two-site observables must match the product of both local dimensions."""
+    psi_mps = MPS(length=2, physical_dimensions=[3, 2], state="zeros")
+    observable = Observable(BaseGate(np.eye(4)), [0, 1])
+
+    with pytest.raises(ValueError, match="does not match site dimensions 3 and 2"):
+        psi_mps.local_expect(observable, [0, 1])
+
+
+def test_local_expect_rejects_observables_with_unsupported_interaction() -> None:
+    """Local expectation values support at most two-site observables."""
+    psi_mps = MPS(length=3, state="zeros")
+    observable = Observable(BaseGate(np.eye(8)), [0, 1, 2])
+
+    with pytest.raises(ValueError, match="Local observable must be one-site or nearest-neighbor two-site"):
+        psi_mps.local_expect(observable, [0, 1, 2])
+
+
+def test_apply_local_rejects_one_site_observable_with_multiple_sites() -> None:
+    """One-site observables retain a defensive site-count check during application."""
+    psi_mps = MPS(length=2, state="zeros")
+    observable = Observable(X(), 0)
+    observable.sites = [0, 1]
+
+    with pytest.raises(ValueError, match=r"One-site local observable requires one site, got \[0, 1\]"):
+        psi_mps.apply_local(observable)
+
+
 def test_mps_apply_local_l2_periodic_wrap_matches_permuted_nn() -> None:
     """For ``L == 2``, wrap-ordered and permuted NN applications must agree."""
     length = 2
diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py
index e4a4d6d7d..a543d5c7b 100644
--- a/tests/core/data_structures/test_simulation_parameters.py
+++ b/tests/core/data_structures/test_simulation_parameters.py
@@ -410,6 +410,16 @@ def test_observable_from_string_falls_back_to_pvm() -> None:
     assert np.allclose(obs.gate.matrix, np.eye(2))
 
 
+@pytest.mark.parametrize("gate_name", ["pvm", "rx"])
+def test_observable_from_non_default_constructible_string_falls_back_to_pvm(gate_name: str) -> None:
+    """PVM and parameterized gate names use the string as a measurement bitstring."""
+    obs = Observable(gate_name)
+
+    assert obs.gate.name == "pvm"
+    assert hasattr(obs.gate, "bitstring")
+    assert obs.gate.bitstring == gate_name
+
+
 def test_observable_from_gate_instance_keeps_gate_and_sites_int() -> None:
     """Passing a concrete BaseGate instance should be preserved and sites can be an int."""
     x_gate = GateLibrary.x()

From ad1a58ac57629e05bd8255743880b68a5df5d17c Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 29 Jul 2026 15:14:26 +0000
Subject: [PATCH 5/9] =?UTF-8?q?=F0=9F=8E=A8=20pre-commit=20fixes?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md                 | 3 ++-
 docs/examples/trapped_ion.md | 4 ++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 75347913b..67ac90340 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,7 +12,8 @@ releases may include breaking changes.
 
 ### Added
 
-- added custom one-site matrix observables for non-qubit local dimensions ([#497]) ([**@linusschulte**])
+- added custom one-site matrix observables for non-qubit local dimensions
+  ([#497]) ([**@linusschulte**])
 - added direct MPO Process Tensor construction and temporal entanglement
   ([#508]) ([**@aaronleesander**])
 
diff --git a/docs/examples/trapped_ion.md b/docs/examples/trapped_ion.md
index f4ac908a2..084c7ee01 100644
--- a/docs/examples/trapped_ion.md
+++ b/docs/examples/trapped_ion.md
@@ -71,8 +71,8 @@ final_x = float(position_expectation[-1])
 ```
 
 The position observable is a custom one-site matrix on the grid basis. The final
-$\langle x\rangle$ is close to $-x_0$ but not exact because the simulation uses a finite
-grid and a finite-difference kinetic operator.
+$\langle x\rangle$ is close to $-x_0$ but not exact because the simulation uses
+a finite grid and a finite-difference kinetic operator.
 
 ```{code-cell} ipython3
 print(f"Initial <x>       = {initial_displacement:.6f}")

From d9d1cf65a6fa4f99abd88536d487cae3677fdc26 Mon Sep 17 00:00:00 2001
From: linusschulte <linus.schulte@gmx.com>
Date: Thu, 30 Jul 2026 12:11:11 +0200
Subject: [PATCH 6/9] made positions a named observable type requiring
 appropriate kwargs

---
 CHANGELOG.md                                  |  2 +-
 docs/examples/simulation_parameters.md        |  4 ++
 docs/examples/trapped_ion.md                  |  2 +-
 .../data_structures/simulation_parameters.py  | 34 +++++++----
 src/mqt/yaqs/core/libraries/gate_library.py   | 26 +++++++++
 tests/core/data_structures/test_mpo.py        |  2 +-
 .../test_simulation_parameters.py             | 57 +++++++++++++++++--
 7 files changed, 109 insertions(+), 18 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 67ac90340..211aaa28e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,7 +13,7 @@ releases may include breaking changes.
 ### Added
 
 - added custom one-site matrix observables for non-qubit local dimensions
-  ([#497]) ([**@linusschulte**])
+ , including named position observables ([#497]) ([**@linusschulte**])
 - added direct MPO Process Tensor construction and temporal entanglement
   ([#508]) ([**@aaronleesander**])
 
diff --git a/docs/examples/simulation_parameters.md b/docs/examples/simulation_parameters.md
index 4f9a3a8e6..64542ce45 100644
--- a/docs/examples/simulation_parameters.md
+++ b/docs/examples/simulation_parameters.md
@@ -34,6 +34,7 @@ do not import gate classes for standard measurements.
 | `"x"`, `"y"`, `"z"`            | Single-qubit Pauli operators                                   | `Observable("z", sites=0)`                  |
 | `"h"`, `"s"`, `"t"`, `"rx"`, … | Other single-qubit gates from the built-in library             | `Observable("h", sites=0)`                  |
 | `"xx"`, `"yy"`, `"zz"`         | Two-qubit Pauli strings                                        | `Observable("zz", sites=[0, 1])`            |
+| `"position"`                   | Position operator for a supplied local position basis          | `Observable("position", 0, positions=grid)` |
 | `"entropy"`                    | Bipartite entanglement entropy across a cut                    | `Observable("entropy", sites=cut)`          |
 | `"schmidt_spectrum"`           | Schmidt spectrum across a cut                                  | `Observable("schmidt_spectrum", sites=cut)` |
 | bitstring / `"pvm"`            | Projection-valued measurement onto a computational basis state | see {doc}`strong_simulation`                |
@@ -41,6 +42,9 @@ do not import gate classes for standard measurements.
 For custom unitaries and circuit gates, use {doc}`custom_gates` — those
 workflows still use `GateLibrary` or Qiskit circuits directly.
 
+Named observables that require configuration accept keyword-only factory arguments. Missing or unknown arguments raise
+`TypeError`, so misspelled parameters are not silently ignored.
+
 ## Start with a preset
 
 You do **not** need to tune every numerical knob before running a simulation.
diff --git a/docs/examples/trapped_ion.md b/docs/examples/trapped_ion.md
index 084c7ee01..c4c7c71e2 100644
--- a/docs/examples/trapped_ion.md
+++ b/docs/examples/trapped_ion.md
@@ -44,7 +44,7 @@ initial_grid_state /= np.linalg.norm(initial_grid_state)
 
 hamiltonian = Hamiltonian.from_mpo(MPO.trapped_ion(positions, masses=[1.0], omega=omega))
 state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim])
-position_observable = Observable(np.diag(positions), 0)
+position_observable = Observable("position", 0, positions=positions)
 ```
 
 ## 2. Noiseless evolution to $T/2$
diff --git a/src/mqt/yaqs/core/data_structures/simulation_parameters.py b/src/mqt/yaqs/core/data_structures/simulation_parameters.py
index 90f8b2540..b4fa01e69 100644
--- a/src/mqt/yaqs/core/data_structures/simulation_parameters.py
+++ b/src/mqt/yaqs/core/data_structures/simulation_parameters.py
@@ -233,31 +233,45 @@ class Observable:
         sites: The site or site indices on which this observable is measured.
     """
 
-    def __init__(self, gate: BaseGate | str | ArrayLike, sites: int | list[int] | None = None) -> None:
+    def __init__(
+        self,
+        gate: BaseGate | str | ArrayLike,
+        sites: int | list[int] | None = None,
+        **gate_kwargs: object,
+    ) -> None:
         """Initializes an Observable instance.
 
         Args:
             gate: The gate or one-site local matrix that will act as the observable.
             sites: The qubit or site indices on which this observable is measured.
+            **gate_kwargs: Keyword-only arguments for a named gate or observable factory.
+
+        Raises:
+            TypeError: If factory arguments are missing, unexpected, or supplied for a gate instance or matrix.
         """
         if isinstance(gate, str):
-            if gate == "entropy":
-                resolved_gate = GateLibrary.entropy()
-            elif gate == "schmidt_spectrum":
-                resolved_gate = GateLibrary.schmidt_spectrum()
-            elif gate == "pvm":
+            if gate == "pvm":
+                if gate_kwargs:
+                    msg = "'pvm' does not accept observable parameters."
+                    raise TypeError(msg)
                 resolved_gate = GateLibrary.pvm(gate)
             elif hasattr(GateLibrary, gate):
                 attr = getattr(GateLibrary, gate)
-                try:
-                    resolved_gate = attr()
-                except TypeError:
-                    resolved_gate = GateLibrary.pvm(gate)
+                resolved_gate = attr(**gate_kwargs)
             else:
+                if gate_kwargs:
+                    msg = f"Unknown observable {gate!r} does not accept observable parameters."
+                    raise TypeError(msg)
                 resolved_gate = GateLibrary.pvm(gate)
         elif isinstance(gate, BaseGate):
+            if gate_kwargs:
+                msg = "Observable parameters are only supported for named observables."
+                raise TypeError(msg)
             resolved_gate = gate
         else:
+            if gate_kwargs:
+                msg = "Observable parameters are only supported for named observables."
+                raise TypeError(msg)
             resolved_gate = GateLibrary.local(gate)
         assert hasattr(GateLibrary, resolved_gate.name), f"Observable {resolved_gate.name} not found in GateLibrary."
         self.gate: BaseGate = copy.deepcopy(resolved_gate)
diff --git a/src/mqt/yaqs/core/libraries/gate_library.py b/src/mqt/yaqs/core/libraries/gate_library.py
index 20a8a7274..694ff4d9e 100644
--- a/src/mqt/yaqs/core/libraries/gate_library.py
+++ b/src/mqt/yaqs/core/libraries/gate_library.py
@@ -1681,6 +1681,30 @@ def __init__(self, matrix: ArrayLike) -> None:
         self.interaction = 1
 
 
+class Position(LocalOperator):
+    """One-site position operator for a supplied position basis."""
+
+    name = "position"
+
+    def __init__(self, *, positions: ArrayLike) -> None:
+        """Create a position operator that is diagonal in the supplied basis.
+
+        Args:
+            positions: One-dimensional position values defining the local basis.
+
+        Raises:
+            ValueError: If ``positions`` is not a non-empty, finite one-dimensional array.
+        """
+        position_values = np.asarray(positions, dtype=np.float64)
+        if position_values.ndim != 1 or position_values.size == 0:
+            msg = "positions must be a non-empty one-dimensional array."
+            raise ValueError(msg)
+        if not np.all(np.isfinite(position_values)):
+            msg = "positions must contain only finite values."
+            raise ValueError(msg)
+        super().__init__(np.diag(position_values))
+
+
 class Entropy(BaseGate):
     """Meta-observable for bipartite entanglement entropy across a cut.
 
@@ -1793,6 +1817,7 @@ class GateLibrary:
         p1: Class for projector ``|1⟩⟨1|``.
         pvm: Class for projection-valued measurement onto a given bitstring.
         local: Class for arbitrary one-site local operators.
+        position: Class for a one-site position operator in a supplied position basis.
 
         entropy:      Class representing a request for bipartite entanglement entropy across a cut.
         schmidt_spectrum: Class representing a request for the Schmidt spectrum across a cut.
@@ -1844,6 +1869,7 @@ class GateLibrary:
     p1 = P1
     pvm = PVM
     local = LocalOperator
+    position = Position
 
     entropy = Entropy
     schmidt_spectrum = SchmidtSpectrum
diff --git a/tests/core/data_structures/test_mpo.py b/tests/core/data_structures/test_mpo.py
index c3a3a72d1..f8756b137 100644
--- a/tests/core/data_structures/test_mpo.py
+++ b/tests/core/data_structures/test_mpo.py
@@ -526,7 +526,7 @@ def test_trapped_ion_one_ion_position_observable_centers_on_trap() -> None:
         physical_dimensions=[positions.size],
     )
 
-    position = Observable(np.diag(positions), 0)
+    position = Observable("position", 0, positions=positions)
 
     np.testing.assert_allclose(mps.expect(position), trap_center, atol=6e-2)
 
diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py
index 00452b888..aaac92e35 100644
--- a/tests/core/data_structures/test_simulation_parameters.py
+++ b/tests/core/data_structures/test_simulation_parameters.py
@@ -74,6 +74,54 @@ def test_observable_accepts_custom_local_matrix() -> None:
     assert obs.sites == 0
 
 
+def test_observable_accepts_named_position_operator() -> None:
+    """Position observables build a diagonal local operator from the supplied basis."""
+    positions = np.array([-1.5, 0.0, 2.5])
+
+    obs = Observable("position", 1, positions=positions)
+
+    assert obs.gate.name == "position"
+    assert obs.gate.interaction == 1
+    np.testing.assert_allclose(obs.gate.matrix, np.diag(positions))
+    assert obs.sites == 1
+
+
+def test_position_observable_requires_positions() -> None:
+    """Position observables require their basis values as a keyword argument."""
+    with pytest.raises(TypeError, match="required keyword-only argument: 'positions'"):
+        Observable("position", 0)
+
+
+@pytest.mark.parametrize(
+    ("gate", "kwargs", "match"),
+    [
+        ("position", {"position_values": [0.0, 1.0]}, "unexpected keyword argument 'position_values'"),
+        ("z", {"positions": [0.0, 1.0]}, "unexpected keyword argument 'positions'"),
+    ],
+)
+def test_named_observable_rejects_unexpected_parameters(
+    gate: str,
+    kwargs: dict[str, object],
+    match: str,
+) -> None:
+    """Named observable factories reject misspelled or inapplicable parameters."""
+    with pytest.raises(TypeError, match=match):
+        Observable(gate, 0, **kwargs)
+
+
+def test_matrix_observable_rejects_named_parameters() -> None:
+    """Factory parameters cannot be supplied with a matrix observable."""
+    with pytest.raises(TypeError, match="only supported for named observables"):
+        Observable(np.eye(2), 0, positions=[0.0, 1.0])
+
+
+@pytest.mark.parametrize("positions", [np.array([]), np.array([0.0, np.nan])])
+def test_position_observable_rejects_invalid_positions(positions: np.ndarray) -> None:
+    """Position bases must be non-empty and finite."""
+    with pytest.raises(ValueError, match="positions must"):
+        Observable("position", 0, positions=positions)
+
+
 @pytest.mark.parametrize("matrix", [np.ones(3), np.ones((2, 3))])
 def test_observable_rejects_invalid_custom_local_matrix(matrix: np.ndarray) -> None:
     """Matrix observables must be two-dimensional and square."""
@@ -410,14 +458,13 @@ def test_observable_from_string_falls_back_to_pvm() -> None:
     assert np.allclose(obs.gate.matrix, np.eye(2))
 
 
-@pytest.mark.parametrize("gate_name", ["pvm", "rx"])
-def test_observable_from_non_default_constructible_string_falls_back_to_pvm(gate_name: str) -> None:
-    """PVM and parameterized gate names use the string as a measurement bitstring."""
-    obs = Observable(gate_name)
+def test_observable_from_explicit_pvm_string_falls_back_to_pvm() -> None:
+    """The explicit PVM name retains its historical string-resolution behavior."""
+    obs = Observable("pvm")
 
     assert obs.gate.name == "pvm"
     assert hasattr(obs.gate, "bitstring")
-    assert obs.gate.bitstring == gate_name
+    assert obs.gate.bitstring == "pvm"
 
 
 def test_observable_from_gate_instance_keeps_gate_and_sites_int() -> None:

From 0ec45feeb6dcc18472f4710e7fec9b27bd0cea6f Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 30 Jul 2026 10:25:30 +0000
Subject: [PATCH 7/9] =?UTF-8?q?=F0=9F=8E=A8=20pre-commit=20fixes?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md                           | 4 +++-
 docs/examples/simulation_parameters.md | 5 +++--
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index d1245ab1b..59ec78762 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,7 +13,9 @@ releases may include breaking changes.
 ### Added
 
 - added custom one-site matrix observables for non-qubit local dimensions
- , including named position observables ([#497]) ([**@linusschulte**])
+
+  , including named position observables ([#497]) ([**@linusschulte**])
+
 - added direct MPO Process Tensor construction and temporal entanglement
   ([#508]) ([**@aaronleesander**])
 
diff --git a/docs/examples/simulation_parameters.md b/docs/examples/simulation_parameters.md
index 729277653..8df8c166d 100644
--- a/docs/examples/simulation_parameters.md
+++ b/docs/examples/simulation_parameters.md
@@ -41,8 +41,9 @@ do not import gate classes for standard measurements.
 For custom unitaries and circuit gates, use {doc}`custom_gates` — those
 workflows still use `GateLibrary` or Qiskit circuits directly.
 
-Named observables that require configuration accept keyword-only factory arguments. Missing or unknown arguments raise
-`TypeError`, so misspelled parameters are not silently ignored.
+Named observables that require configuration accept keyword-only factory
+arguments. Missing or unknown arguments raise `TypeError`, so misspelled
+parameters are not silently ignored.
 
 ## Start with a preset
 

From c1adc2e86b8831f52f29e37d1c954e3384661a4f Mon Sep 17 00:00:00 2001
From: linusschulte <linus.schulte@gmx.com>
Date: Thu, 30 Jul 2026 13:49:10 +0200
Subject: [PATCH 8/9] add test coverage

---
 tests/core/data_structures/test_mps.py             |  8 ++++++++
 .../data_structures/test_simulation_parameters.py  | 14 +++++++++++++-
 2 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/tests/core/data_structures/test_mps.py b/tests/core/data_structures/test_mps.py
index 99973bdf0..b73bcc02e 100644
--- a/tests/core/data_structures/test_mps.py
+++ b/tests/core/data_structures/test_mps.py
@@ -502,6 +502,14 @@ def test_apply_local_rejects_one_site_observable_with_multiple_sites() -> None:
         psi_mps.apply_local(observable)
 
 
+def test_apply_local_rejects_mismatched_one_site_dimension() -> None:
+    """One-site observables must match the local dimension during direct application."""
+    psi_mps = MPS(length=1, physical_dimensions=[3], state="zeros")
+
+    with pytest.raises(ValueError, match="does not match site 0 dimension 3"):
+        psi_mps.apply_local(Observable(np.eye(2), 0))
+
+
 def test_mps_apply_local_l2_periodic_wrap_matches_permuted_nn() -> None:
     """For ``L == 2``, wrap-ordered and permuted NN applications must agree."""
     length = 2
diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py
index b7a354918..27ee02c69 100644
--- a/tests/core/data_structures/test_simulation_parameters.py
+++ b/tests/core/data_structures/test_simulation_parameters.py
@@ -35,7 +35,7 @@
     Observable,
     _validate_tdvp_sweeps,
 )
-from mqt.yaqs.core.libraries.gate_library import GateLibrary, X
+from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary, X
 from mqt.yaqs.core.methods.tdvp import primitives as tdvp_primitives
 
 if TYPE_CHECKING:
@@ -114,6 +114,18 @@ def test_matrix_observable_rejects_named_parameters() -> None:
         Observable(np.eye(2), 0, positions=[0.0, 1.0])
 
 
+def test_observable_rejects_parameters_without_a_matching_factory() -> None:
+    """Only recognized named factories accept additional observable parameters."""
+    with pytest.raises(TypeError, match="'pvm' does not accept observable parameters"):
+        Observable("pvm", bitstring="0")
+
+    with pytest.raises(TypeError, match="Unknown observable 'unknown'"):
+        Observable("unknown", parameter=1)
+
+    with pytest.raises(TypeError, match="only supported for named observables"):
+        Observable(BaseGate(np.eye(2)), 0, parameter=1)
+
+
 @pytest.mark.parametrize("positions", [np.array([]), np.array([0.0, np.nan])])
 def test_position_observable_rejects_invalid_positions(positions: np.ndarray) -> None:
     """Position bases must be non-empty and finite."""

From 5c01f4fa2adccf531d263d58df0f02ee2b67e320 Mon Sep 17 00:00:00 2001
From: linusschulte <linus.schulte@gmx.com>
Date: Thu, 30 Jul 2026 14:05:30 +0200
Subject: [PATCH 9/9] address minor stability feedback

---
 src/mqt/yaqs/core/data_structures/mps.py           |  8 ++++++++
 src/mqt/yaqs/core/libraries/gate_library.py        |  8 ++++++--
 tests/core/data_structures/test_mps.py             | 14 ++++++++++++++
 .../data_structures/test_simulation_parameters.py  |  4 ++--
 4 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/src/mqt/yaqs/core/data_structures/mps.py b/src/mqt/yaqs/core/data_structures/mps.py
index 83a024125..43b831291 100644
--- a/src/mqt/yaqs/core/data_structures/mps.py
+++ b/src/mqt/yaqs/core/data_structures/mps.py
@@ -1119,9 +1119,17 @@ def bubble_swaps_backward(state: MPS) -> None:
             return
 
         if observable.gate.interaction == 2:
+            if len(sites) != 2:
+                msg = f"Two-site local observable requires two sites, got {sites}."
+                raise ValueError(msg)
             i, j = int(sites[0]), int(sites[1])
             length = self.length
             mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
+            d_i = self.tensors[i].shape[0]
+            d_j = self.tensors[j].shape[0]
+            if mat.shape != (d_i * d_j, d_i * d_j):
+                msg = f"Two-site observable matrix shape {mat.shape} does not match site dimensions {d_i} and {d_j}."
+                raise ValueError(msg)
 
             if length == 2:
                 if i == length - 1 and j == 0:
diff --git a/src/mqt/yaqs/core/libraries/gate_library.py b/src/mqt/yaqs/core/libraries/gate_library.py
index 694ff4d9e..ea3de5e2b 100644
--- a/src/mqt/yaqs/core/libraries/gate_library.py
+++ b/src/mqt/yaqs/core/libraries/gate_library.py
@@ -1693,9 +1693,13 @@ def __init__(self, *, positions: ArrayLike) -> None:
             positions: One-dimensional position values defining the local basis.
 
         Raises:
-            ValueError: If ``positions`` is not a non-empty, finite one-dimensional array.
+            ValueError: If ``positions`` is complex or not a non-empty, finite one-dimensional array.
         """
-        position_values = np.asarray(positions, dtype=np.float64)
+        position_values = np.asarray(positions)
+        if np.iscomplexobj(position_values):
+            msg = "positions must contain only real values."
+            raise ValueError(msg)
+        position_values = np.asarray(position_values, dtype=np.float64)
         if position_values.ndim != 1 or position_values.size == 0:
             msg = "positions must be a non-empty one-dimensional array."
             raise ValueError(msg)
diff --git a/tests/core/data_structures/test_mps.py b/tests/core/data_structures/test_mps.py
index b73bcc02e..1c3472567 100644
--- a/tests/core/data_structures/test_mps.py
+++ b/tests/core/data_structures/test_mps.py
@@ -510,6 +510,20 @@ def test_apply_local_rejects_mismatched_one_site_dimension() -> None:
         psi_mps.apply_local(Observable(np.eye(2), 0))
 
 
+def test_apply_local_rejects_invalid_two_site_observable_shape() -> None:
+    """Two-site observables must supply two sites and match both local dimensions."""
+    psi_mps = MPS(length=2, physical_dimensions=[3, 2], state="zeros")
+    observable = Observable(BaseGate(np.eye(4)), [0, 1])
+
+    observable.sites = [0]
+    with pytest.raises(ValueError, match=r"requires two sites, got \[0\]"):
+        psi_mps.apply_local(observable)
+
+    observable.sites = [0, 1]
+    with pytest.raises(ValueError, match="does not match site dimensions 3 and 2"):
+        psi_mps.apply_local(observable)
+
+
 def test_mps_apply_local_l2_periodic_wrap_matches_permuted_nn() -> None:
     """For ``L == 2``, wrap-ordered and permuted NN applications must agree."""
     length = 2
diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py
index 27ee02c69..73db2262d 100644
--- a/tests/core/data_structures/test_simulation_parameters.py
+++ b/tests/core/data_structures/test_simulation_parameters.py
@@ -126,9 +126,9 @@ def test_observable_rejects_parameters_without_a_matching_factory() -> None:
         Observable(BaseGate(np.eye(2)), 0, parameter=1)
 
 
-@pytest.mark.parametrize("positions", [np.array([]), np.array([0.0, np.nan])])
+@pytest.mark.parametrize("positions", [np.array([]), np.array([0.0, np.nan]), np.array([0.0, 1.0j])])
 def test_position_observable_rejects_invalid_positions(positions: np.ndarray) -> None:
-    """Position bases must be non-empty and finite."""
+    """Position bases must be non-empty, finite, and real."""
     with pytest.raises(ValueError, match="positions must"):
         Observable("position", 0, positions=positions)