Repository: PennyLaneAI/pennylane
Pull request #7034 (2025-03-04 12:49:56)

Adding Realspace Hamiltonian classes to Labs
This PR introduces the Realspace classes for Trotter error estimation to labs.

[sc-87023]

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7034?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.60%. Comparing base [(`494d18a`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/494d18aedac8e3b5442e38e9437e23e6329211bb?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`bd17ef3`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/bd17ef3dfe88b838939edec97fa56d046eef5ba6?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> Report is 1 commits behind head on master.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##           master    #7034   +/-   ##
=======================================
  Coverage   99.60%   99.60%           
=======================================
  Files         488      488           
  Lines       46459    46459           
=======================================
  Hits        46276    46276           
  Misses        183      183           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7034?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 01418fea7..41f2cd157 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -634,6 +634,52 @@
   [(#7026)](https://github.com/PennyLaneAI/pennylane/pull/7026)
   [(#7054)](https://github.com/PennyLaneAI/pennylane/pull/7054)
 
+* Adding `HOState` and `VibronicHO` classes for representing harmonic oscillator states.
+  [(#7035)](https://github.com/PennyLaneAI/pennylane/pull/7035)
+
+* Adding base classes for Trotter error estimation on Realspace Hamiltonians: ``RealspaceOperator``, ``RealspaceSum``, ``RealspaceCoeffs``, and ``RealspaceMatrix``
+  [(#7034)](https://github.com/PennyLaneAI/pennylane/pull/7034)
+
+* Adding functions for Trotter error estimation and Hamiltonian fragment generation: ``trotter_error``, ``perturbation_error``, ``vibrational_fragments``, ``vibronic_fragments``, and ``generic_fragments``.
+  [(#7036)](https://github.com/PennyLaneAI/pennylane/pull/7036)
+
+  As an example we compute the peruturbation error of a vibrational Hamiltonian.
+  First we generate random harmonic frequences and Taylor coefficients to iniitialize the vibrational Hamiltonian.
+
+  ```pycon
+  >>> from pennylane.labs.trotter_error import HOState, vibrational_fragments, perturbation_error
+  >>> import numpy as np
+  >>> n_modes = 2
+  >>> r_state = np.random.RandomState(42)
+  >>> freqs = r_state.random(n_modes)
+  >>> taylor_coeffs = [
+  >>>     np.array(0),
+  >>>     r_state.random(size=(n_modes, )),
+  >>>     r_state.random(size=(n_modes, n_modes)),
+  >>>     r_state.random(size=(n_modes, n_modes, n_modes))
+  >>> ]
+  ```
+    
+  We call ``vibrational_fragments`` to get the harmonic and anharmonic fragments of the vibrational Hamiltonian.
+  ```pycon
+  >>> frags = vibrational_fragments(n_modes, freqs, taylor_coeffs)
+  ```
+
+  We build state vectors in the harmonic oscilator basis with the ``HOState`` class. 
+
+  ```pycon
+  >>> gridpoints = 5
+  >>> state1 = HOState(n_modes, gridpoints, {(0, 0): 1})
+  >>> state2 = HOState(n_modes, gridpoints, {(1, 1): 1})
+  ```
+
+  Finally, we compute the error by calling ``perturbation_error``.
+
+  ```pycon
+  >>> perturbation_error(frags, [state1, state2])
+  [(-0.9189251160920879+0j), (-4.797716682426851+0j)]
+  ```
+
 <h3>Breaking changes 💔</h3>
 
 * `num_diagonalizing_gates` is no longer accessible in `qml.specs` or `QuantumScript.specs`. The calculation of
diff --git a/pennylane/labs/__init__.py b/pennylane/labs/__init__.py
index 4b1e8da42..1ef3a9e89 100644
--- a/pennylane/labs/__init__.py
+++ b/pennylane/labs/__init__.py
@@ -27,11 +27,13 @@ Modules
 
     dla
     resource_estimation
+    trotter_error
 
 """
 
 from pennylane.labs import dla
 from pennylane.labs import resource_estimation
 from pennylane.labs import vibrational
+from pennylane.labs import trotter_error
 
 __all__ = []
diff --git a/pennylane/labs/trotter_error/__init__.py b/pennylane/labs/trotter_error/__init__.py
new file mode 100644
index 000000000..2c75ca3b4
--- /dev/null
+++ b/pennylane/labs/trotter_error/__init__.py
@@ -0,0 +1,96 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Module containing functionality for computing Trotter error.
+
+.. warning::
+
+    This module is experimental. Frequent changes will occur,
+    with no guarantees of stability or backwards compatibility.
+
+.. currentmodule:: pennylane.labs.trotter_error
+
+Trotter Base Classes
+~~~~~~~~~~~~~~~~~~~~
+Abstract classes specifying the methods needed for implementing fragments and states in the Trotter error workflow.
+
+.. autosummary::
+    :toctree: api
+
+    ~AbstractState
+    ~Fragment
+
+Fragment Classes
+~~~~~~~~~~~~~~~~
+Classes representing fragments of different types of Hamiltonians.
+
+.. autosummary::
+    :toctree: api
+
+    ~GenericFragment
+    ~RealspaceMatrix
+    ~RealspaceSum
+
+Realspace Hamiltonian Classes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Classes used to store representations of realspace Hamiltonians.
+
+.. autosummary::
+    :toctree: api
+
+    ~RealspaceCoeffs
+    ~RealspaceOperator
+
+Fragment Functions
+~~~~~~~~~~~~~~~~~~
+Functions used to retrieve fragments of various Hamiltonians.
+
+.. autosummary::
+    :toctree: api
+
+    ~generic_fragments
+    ~vibrational_fragments
+    ~vibronic_fragments
+
+Harmonic Oscillator Classes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Classes used to construct the representation of harmonic oscillator states.
+
+.. autosummary::
+    :toctree: api
+
+    ~HOState
+    ~VibronicHO
+
+Error Estimation Functions
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+Functions for computing Trotter error estimates.
+
+.. autosummary::
+    :toctree: api
+
+    ~perturbation_error
+    ~trotter_error
+"""
+
+from .abstract import AbstractState, Fragment
+from .fragments import GenericFragment, generic_fragments, vibrational_fragments, vibronic_fragments
+from .product_formulas import perturbation_error, trotter_error
+from .realspace import (
+    RealspaceCoeffs,
+    RealspaceMatrix,
+    RealspaceOperator,
+    RealspaceSum,
+    HOState,
+    VibronicHO,
+)
diff --git a/pennylane/labs/trotter_error/abstract.py b/pennylane/labs/trotter_error/abstract.py
new file mode 100644
index 000000000..82e4903e2
--- /dev/null
+++ b/pennylane/labs/trotter_error/abstract.py
@@ -0,0 +1,173 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""The Fragment class"""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from typing import Dict, Sequence
+
+
+class Fragment(ABC):
+    """Abstract class used to define a fragment object for product formula error estimation.
+
+    A :class:`~.Fragment` is an object that has a well-defined notion of a commutator. To ensure
+    the existence of commutators, the implementation requires the following arithmetic dunder
+    methods:
+
+    * :meth:`~.__add__`: implements addition
+
+    * :meth:`~.__mul__`: implements multiplication
+
+    * :meth:`~.__matmul__`: implements matrix multiplication
+
+    In addition to the arithmetic operators, a ``norm`` method should be defined. The norm is
+    required to compute error estimates of Trotter error operators.
+    """
+
+    @abstractmethod
+    def __add__(self, other: Fragment) -> Fragment:
+        raise NotImplementedError
+
+    def __sub__(self, other: Fragment) -> Fragment:
+        return self + (-1) * other
+
+    @abstractmethod
+    def __mul__(self, scalar: float) -> Fragment:
+        raise NotImplementedError
+
+    @abstractmethod
+    def __matmul__(self, other: Fragment) -> Fragment:
+        raise NotImplementedError
+
+    @abstractmethod
+    def norm(self, params: Dict) -> float:
+        """Compute the norm of the fragment.
+
+        Args:
+            params (Dict): A dictionary of parameters needed to compute the norm. It should be
+                specified for each class inheriting from :class:`~.Fragment`.
+
+        Returns:
+            float: the norm of the :class:`~.Fragment`
+
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def apply(self, state: AbstractState) -> AbstractState:
+        """Apply the Fragment to a state on the right. The type of ``state`` is determined by each class inheriting from ``Fragment``.
+
+        Args:
+            state (AbstractState): an object representing a quantum state
+
+        Returns:
+            AbstractState: the result of applying the ``Fragment`` to ``state``
+
+        """
+        raise NotImplementedError
+
+    def expectation(self, left: AbstractState, right: AbstractState) -> float:
+        """Return the expectation value of a state. The type of ``state`` is determined by each class inheriting from ``Fragment``.
+
+        Args:
+            left (AbstractState): the state to be multiplied on the left of the ``Fragment``
+            right (AbstractState): the state to be multiplied on the right of the ``Fragment``
+
+        Returns:
+            float: the expectation value obtained by applying ``Fragment`` to the given states
+        """
+        return left.dot(self.apply(right))
+
+
+def commutator(a: Fragment, b: Fragment) -> Fragment:
+    """Return the commutator of two :class:`~.Fragment` objects
+
+    Args:
+        a (Fragment): the :class:`~.Fragment` on the left side of the commutator
+        b (Fragment): the :class:`~.Fragment` on the right side of the commutator
+
+    Returns:
+        Fragment: the commutator ``[a, b]``
+    """
+    return a @ b - b @ a
+
+
+def nested_commutator(fragments: Sequence[Fragment]) -> Fragment:
+    """Return the nested commutator of a sequence of :class:`~.Fragment` objects
+
+    Args:
+        fragments (Sequence[Fragment]): a sequence of fragments
+
+    Returns:
+        Fragment: the nested commutator of the fragments
+    """
+
+    if len(fragments) < 2:
+        raise ValueError("Need at least two fragments to commute.")
+
+    if len(fragments) == 2:
+        return commutator(*fragments)
+
+    head, *tail = fragments
+
+    return commutator(head, nested_commutator(tail))
+
+
+class AbstractState(ABC):
+    """Abstract class used to define a state object for product formula error estimation.
+
+    A class inheriting from ``AbstractState`` must implement the following dunder methods.
+
+    * ``__add__``: implements addition
+    * ``__mul__``: implements multiplication
+
+    Additionally, it requires the following methods.
+
+    * ``zero_state``: returns a representation of the zero state
+    * ``dot``: implments the dot product of two states
+    """
+
+    @abstractmethod
+    def __add__(self, other: AbstractState) -> AbstractState:
+        raise NotImplementedError
+
+    def __sub__(self, other: AbstractState) -> AbstractState:
+        return self + (-1) * other
+
+    @abstractmethod
+    def __mul__(self, scalar: float) -> AbstractState:
+        raise NotImplementedError
+
+    @classmethod
+    @abstractmethod
+    def zero_state(cls) -> AbstractState:
+        """Return a representation of the zero state.
+
+        Returns:
+            AbstractState: an ``AbstractState`` representation of the zero state
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def dot(self, other: AbstractState) -> float:
+        """Compute the dot product of two states.
+
+        Args:
+            other (AbstractState): the state to take the dot product with
+
+        Returns:
+           float: the dot product of self and other
+        """
+        raise NotImplementedError
diff --git a/pennylane/labs/trotter_error/fragments/__init__.py b/pennylane/labs/trotter_error/fragments/__init__.py
new file mode 100644
index 000000000..270d576ec
--- /dev/null
+++ b/pennylane/labs/trotter_error/fragments/__init__.py
@@ -0,0 +1,18 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Module containing functions to generate fragments for various Hamiltonians"""
+
+from .generic_fragments import GenericFragment, generic_fragments
+from .vibrational_fragments import vibrational_fragments
+from .vibronic_fragments import vibronic_fragments
diff --git a/pennylane/labs/trotter_error/fragments/generic_fragments.py b/pennylane/labs/trotter_error/fragments/generic_fragments.py
new file mode 100644
index 000000000..f6a208b43
--- /dev/null
+++ b/pennylane/labs/trotter_error/fragments/generic_fragments.py
@@ -0,0 +1,133 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Wrapper class for generic fragment objects"""
+
+from __future__ import annotations
+
+from typing import Any, Callable, Dict, List, Sequence
+
+from pennylane.labs.trotter_error import Fragment
+
+
+def generic_fragments(fragments: Sequence[Any], norm_fn: Callable = None) -> List[GenericFragment]:
+    """Instantiates :class:`~.pennylane.labs.trotter_error.GenericFragment` objects.
+
+    Args:
+        fragments (Sequence[Any]): A sequence of objects of the same type. The type is assumed to implement ``__add__``, ``__mul__``, and ``__matmul__``.
+        norm_fn (Callable): A function that computes the norm of the fragments.
+
+    Returns:
+        List[GenericFragment]: A list of :class:`~.pennylane.labs.trotter_error.GenericFragment` objects instantiated from `fragments`.
+
+
+    **Example**
+
+    This code example demonstrates building fragments from numpy matrices.
+
+    >>> from pennylane.labs.trotter_error import generic_fragments
+    >>> import numpy as np
+    >>> matrices = [np.array([[1, 0], [0, 1]]), np.array([[0, 1], [1, 0]])]
+    >>> fragments = generic_fragments(matrices, norm_fn=np.linalg.norm)
+    >>> fragments
+    [GenericFragment(type=<class 'numpy.ndarray'>), GenericFragment(type=<class 'numpy.ndarray'>)]
+    >>> fragments[0].norm()
+    1.4142135623730951
+    """
+
+    if len(fragments) == 0:
+        return []
+
+    frag_type = type(fragments[0])
+
+    if not all(isinstance(fragment, frag_type) for fragment in fragments):
+        raise TypeError("All fragments must be of the same type.")
+
+    if not hasattr(frag_type, "__add__"):
+        raise TypeError(f"Fragment of type {frag_type} does not implement __add__.")
+
+    if not hasattr(frag_type, "__mul__"):
+        raise TypeError(f"Fragment of type {frag_type} does not implement __mul__.")
+
+    if not hasattr(frag_type, "__matmul__"):
+        raise TypeError(f"Fragment of type {frag_type} does not implement __matmul__.")
+
+    return [GenericFragment(fragment, norm_fn=norm_fn) for fragment in fragments]
+
+
+class GenericFragment(Fragment):
+    """Abstract class used to define a generic fragment object for product formula error estimation.
+
+    This class allows using any object implementing arithmetic dunder methods to be used
+    for product formula error estimation.
+
+    Args:
+        fragment (Any): An object that implements the following arithmetic methods:
+            ``__add__``, ``__sub__``, ``__mul__``, and ``__matmul__``.
+        norm_fn (optional, Callable): A function used to compute the norm of ``fragment``.
+
+    .. note:: :class:`~.pennylane.labs.trotter_error.GenericFragment` objects should be instantated through the ``generic_fragments`` function.
+
+    **Example**
+
+    >>> from pennylane.labs.trotter_error import generic_fragments
+    >>> import numpy as np
+    >>> matrices = [np.array([[1, 0], [0, 1]]), np.array([[0, 1], [1, 0]])]
+    >>> generic_fragments(matrices)
+    [GenericFragment(type=<class 'numpy.ndarray'>), GenericFragment(type=<class 'numpy.ndarray'>)]
+    """
+
+    def __init__(self, fragment: Any, norm_fn: Callable = None):
+        self.fragment = fragment
+        self.norm_fn = norm_fn
+
+    def __add__(self, other: GenericFragment):
+        return GenericFragment(self.fragment + other.fragment, norm_fn=self.norm_fn)
+
+    def __sub__(self, other: GenericFragment):
+        return GenericFragment(self.fragment + (-1) * other.fragment, norm_fn=self.norm_fn)
+
+    def __mul__(self, scalar: float):
+        return GenericFragment(scalar * self.fragment, norm_fn=self.norm_fn)
+
+    def __eq__(self, other: GenericFragment):
+        if not isinstance(self.fragment, type(other)):
+            return False
+
+        return self.fragment == other.fragment
+
+    __rmul__ = __mul__
+
+    def __matmul__(self, other: GenericFragment):
+        return GenericFragment(self.fragment @ other.fragment, norm_fn=self.norm_fn)
+
+    def apply(self, state: Any) -> Any:
+        """Apply the fragment to a state using the underlying object's ``__matmul__`` method."""
+        return self.fragment @ state
+
+    def expectation(self, left: Any, right: Any) -> float:
+        """Compute the expectation value using the underlying object's ``__matmul__`` method."""
+        return left @ self.fragment @ right
+
+    def norm(self, params: Dict = None) -> float:
+        """Compute the norm of the fragment."""
+        if self.norm_fn:
+            params = params or {}
+            return self.norm_fn(self.fragment, **params)
+
+        raise NotImplementedError(
+            "GenericFragment was constructed without specifying the norm function."
+        )
+
+    def __repr__(self):
+        return f"GenericFragment(type={type(self.fragment)})"
diff --git a/pennylane/labs/trotter_error/fragments/vibrational_fragments.py b/pennylane/labs/trotter_error/fragments/vibrational_fragments.py
new file mode 100644
index 000000000..381885dc3
--- /dev/null
+++ b/pennylane/labs/trotter_error/fragments/vibrational_fragments.py
@@ -0,0 +1,159 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""The realspace vibrational Hamiltonian"""
+
+from typing import List, Sequence
+
+import numpy as np
+
+from pennylane.labs.trotter_error.realspace import RealspaceCoeffs, RealspaceOperator, RealspaceSum
+
+
+def vibrational_fragments(
+    modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray], frag_method="harmonic"
+) -> List[RealspaceSum]:
+    """Returns a list of fragments summing to a vibrational Hamiltonian.
+
+    Args:
+        modes (int): the number of vibrational modes
+        freqs (ndarray): the harmonic frequences
+        taylor_coeffs (Sequence[ndarray]): a sequence containing the tensors of coefficients in the Taylor expansion
+        frag_method (string): the fragmentation method, valid options are ``harmonic``, ``kinetic``, and ``position``
+
+    Returns:
+        List[RealspaceSum]: a list of ``RealspaceSum`` objects representing the fragments of the vibrational Hamiltonian
+
+    **Example**
+
+    >>> from pennylane.labs.trotter_error import vibrational_fragments
+    >>> import numpy as np
+    >>> n_modes = 4
+    >>> r_state = np.random.RandomState(42)
+    >>> freqs = r_state.random(4)
+    >>> taylor_coeffs = [np.array(0), r_state.random(size=(n_modes, )), r_state.random(size=(n_modes, n_modes))]
+    >>> fragments = vibrational_fragments(n_modes, freqs, taylor_coeffs)
+    >>> for fragment in fragments:
+    >>>     print(fragment)
+    RealspaceSum((RealspaceOperator(4, ('PP',), omega[idx0]), RealspaceOperator(4, ('QQ',), omega[idx0])))
+    RealspaceSum((RealspaceOperator(4, ('Q',), phi[1][idx0]), RealspaceOperator(4, ('Q', 'Q'), phi[2][idx0,idx1])))
+    """
+
+    if frag_method == "harmonic":
+        return [_harmonic_fragment(modes, freqs), _anharmonic_fragment(modes, taylor_coeffs)]
+
+    if frag_method == "kinetic":
+        return [_kinetic_fragment(modes, freqs), _potential_fragment(modes, freqs, taylor_coeffs)]
+
+    if frag_method == "position":
+        return [_position_fragment(modes, freqs, taylor_coeffs), _momentum_fragment(modes, freqs)]
+
+    raise ValueError(f"{frag_method} is not a valid fragmentation scheme.")
+
+
+def _harmonic_fragment(modes: int, freqs: np.ndarray) -> RealspaceSum:
+    """Returns the fragment of the Hamiltonian corresponding to the harmonic part."""
+    _validate_freqs(modes, freqs)
+
+    coeffs = RealspaceCoeffs(freqs / 2, label="omega")
+    momentum = RealspaceOperator(modes, ("PP",), coeffs)
+    position = RealspaceOperator(modes, ("QQ",), coeffs)
+
+    return RealspaceSum(modes, [momentum, position])
+
+
+def _anharmonic_fragment(modes: int, taylor_coeffs: Sequence[np.ndarray]) -> RealspaceSum:
+    """Returns the fragment of the Hamiltonian corresponding to the anharmonic part."""
+    _validate_taylor_coeffs(modes, taylor_coeffs)
+
+    ops = []
+    for i, phi in enumerate(taylor_coeffs):
+        op = ("Q",) * i
+        coeffs = RealspaceCoeffs(phi, label=f"phi[{i}]")
+        realspace_op = RealspaceOperator(modes, op, coeffs)
+        ops.append(realspace_op)
+
+    return RealspaceSum(modes, ops)
+
+
+def _kinetic_fragment(modes: int, freqs: np.ndarray) -> RealspaceSum:
+    """Returns the fragment of the Hamiltonian corresponding to the kinetic part"""
+    _validate_freqs(modes, freqs)
+
+    coeffs = RealspaceCoeffs(freqs / 2, label="omega")
+    kinetic = RealspaceOperator(modes, ("PP",), coeffs)
+
+    return RealspaceSum(modes, [kinetic])
+
+
+def _potential_fragment(
+    modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]
+) -> RealspaceSum:
+    """Returns the fragment of the Hamiltonian corresponding to the potential part"""
+    _validate_input(modes, freqs, taylor_coeffs)
+
+    ops = []
+    for i, phi in enumerate(taylor_coeffs):
+        op = ("Q",) * i
+        coeffs = RealspaceCoeffs(phi, label=f"phi[{i}]")
+        realspace_op = RealspaceOperator(modes, op, coeffs)
+        ops.append(realspace_op)
+
+    op = ("Q", "Q")
+    diag = np.diag(freqs / 2)
+    coeffs = RealspaceCoeffs(diag, label="omega")
+    realspace_op = RealspaceOperator(modes, op, coeffs)
+    ops.append(realspace_op)
+
+    return RealspaceSum(modes, ops)
+
+
+def _position_fragment(
+    modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]
+) -> RealspaceSum:
+    """Return the position term of the Hamiltonian"""
+    coeffs = RealspaceCoeffs(freqs / 2, label="omega")
+    position = RealspaceOperator(modes, ("QQ",), coeffs)
+
+    return RealspaceSum(modes, [position]) + _anharmonic_fragment(modes, taylor_coeffs)
+
+
+def _momentum_fragment(modes: int, freqs: np.ndarray) -> RealspaceSum:
+    """Return the momentum term of the Hamiltonian"""
+    _validate_freqs(modes, freqs)
+
+    coeffs = RealspaceCoeffs(freqs / 2, label="omega")
+    momentum = RealspaceOperator(modes, ("PP",), coeffs)
+
+    return RealspaceSum(modes, [momentum])
+
+
+def _validate_taylor_coeffs(modes: int, taylor_coeffs: Sequence[np.ndarray]) -> None:
+    """Validate that the Taylor coefficients have the correct shape."""
+    for i, phi in enumerate(taylor_coeffs):
+        shape = (modes,) * i
+
+        if not phi.shape == shape:
+            raise ValueError(f"Expected order {i} to be of shape {shape}, got {phi.shape}.")
+
+
+def _validate_freqs(modes: int, freqs: np.ndarray) -> None:
+    """Validate that the harmonic frequencies have the correct shape."""
+    if not len(freqs) == modes:
+        raise ValueError(f"Expected freqs to be of length {modes}, got {len(freqs)}.")
+
+
+def _validate_input(modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]) -> None:
+    """Validate that the Taylor coefficients and the harmonic frequencies have the correct shape."""
+    _validate_taylor_coeffs(modes, taylor_coeffs)
+    _validate_freqs(modes, freqs)
diff --git a/pennylane/labs/trotter_error/fragments/vibronic_fragments.py b/pennylane/labs/trotter_error/fragments/vibronic_fragments.py
new file mode 100644
index 000000000..94fd4fb75
--- /dev/null
+++ b/pennylane/labs/trotter_error/fragments/vibronic_fragments.py
@@ -0,0 +1,142 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Base class for Vibronic Hamiltonian"""
+
+from __future__ import annotations
+
+from typing import List, Sequence
+
+import numpy as np
+
+from pennylane.labs.trotter_error.realspace import (
+    RealspaceCoeffs,
+    RealspaceMatrix,
+    RealspaceOperator,
+    RealspaceSum,
+)
+from pennylane.labs.trotter_error.realspace.realspace_matrix import _next_pow_2
+
+
+def vibronic_fragments(
+    states: int, modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]
+) -> List[RealspaceMatrix]:
+    """Returns a list of fragments summing to a vibronic Hamiltonian.
+
+    Args:
+        states (int): the number of electronic states
+        modes (int): the number of vibrational modes
+        freqs (ndarray): the harmonic frequences
+        taylor_coeffs (Sequence[ndarray]): a sequence containing the tensors of coefficients in the Taylor expansion
+
+    Returns:
+        List[RealspaceMatrix]: a list of ``RealspaceMatrix`` objects representing the fragments of the vibronic Hamiltonian
+
+    **Example**
+
+    >>> from pennylane.labs.trotter_error import vibronic_fragments
+    >>> import numpy as np
+    >>> n_modes = 4
+    >>> n_states = 2
+    >>> r_state = np.random.RandomState(42)
+    >>> freqs = r_state.random(4)
+    >>> taylor_coeffs = [r_state.random(size=(n_states, n_states, )), r_state.random(size=(n_states, n_states, n_modes))]
+    >>> fragments = vibronic_fragments(n_states, n_modes, freqs, taylor_coeffs)
+    >>> for fragment in fragments:
+    >>>     fragment
+    RealspaceMatrix({(0, 0): RealspaceSum((RealspaceOperator(4, (), 0.15601864044243652), RealspaceOperator(4, ('Q',), phi[1][0, 0][idx0]), RealspaceOperator(4, ('Q', 'Q'), omega[idx0,idx1]))), (1, 1): RealspaceSum((RealspaceOperator(4, (), 0.8661761457749352), RealspaceOperator(4, ('Q',), phi[1][1, 1][idx0]), RealspaceOperator(4, ('Q', 'Q'), omega[idx0,idx1])))})
+    RealspaceMatrix({(0, 1): RealspaceSum((RealspaceOperator(4, (), 0.15599452033620265), RealspaceOperator(4, ('Q',), phi[1][0, 1][idx0]))), (1, 0): RealspaceSum((RealspaceOperator(4, (), 0.05808361216819946), RealspaceOperator(4, ('Q',), phi[1][1, 0][idx0])))})
+    RealspaceMatrix({(0, 0): RealspaceSum((RealspaceOperator(4, ('P', 'P'), omega[idx0,idx1]),)), (1, 1): RealspaceSum((RealspaceOperator(4, ('P', 'P'), omega[idx0,idx1]),))})
+    """
+    _validate_input(states, modes, freqs, taylor_coeffs)
+
+    frags = [
+        _position_fragment(i, states, modes, freqs, taylor_coeffs)
+        for i in range(_next_pow_2(states))
+    ]
+    frags.append(_momentum_fragment(states, modes, freqs))
+
+    return frags
+
+
+def _position_fragment(
+    i: int, states: int, modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]
+) -> RealspaceMatrix:
+    """Return the ``i``th position fragment"""
+    pow2 = _next_pow_2(states)
+    blocks = {
+        (j, i ^ j): _realspace_sum(j, i ^ j, states, modes, freqs, taylor_coeffs)
+        for j in range(pow2)
+    }
+    return RealspaceMatrix(pow2, modes, blocks)
+
+
+def _momentum_fragment(states: int, modes: int, freqs: np.ndarray) -> RealspaceMatrix:
+    """Return the fragment consisting only of momentum operators."""
+    pow2 = _next_pow_2(states)
+    term = RealspaceOperator(
+        modes,
+        ("P", "P"),
+        RealspaceCoeffs(np.diag(freqs) / 2, label="omega"),
+    )
+    word = RealspaceSum(modes, (term,))
+    blocks = {(i, i): word for i in range(states)}
+
+    return RealspaceMatrix(pow2, modes, blocks)
+
+
+# pylint: disable=too-many-arguments,too-many-positional-arguments
+def _realspace_sum(
+    i: int, j: int, states: int, modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]
+) -> RealspaceSum:
+    """Return a RealspaceSum representation of the ``(i, j)`` block in the RealspaceMatrix"""
+    if i > states - 1 or j > states - 1:
+        return RealspaceSum.zero(modes)
+
+    realspace_ops = []
+    for k, phi in enumerate(taylor_coeffs):
+        op = ("Q",) * k
+        realspace_op = RealspaceOperator(
+            modes,
+            op,
+            RealspaceCoeffs(phi[i, j], label=f"phi[{k}][{i}, {j}]"),
+        )
+        realspace_ops.append(realspace_op)
+
+    if i == j:
+        op = ("Q", "Q")
+        coeffs = RealspaceCoeffs(np.diag(freqs) / 2, label="omega")
+        assert coeffs is not None
+        realspace_ops.append(RealspaceOperator(modes, op, coeffs))
+
+    ret_val = RealspaceSum(modes, realspace_ops)
+    for op in ret_val.ops:
+        assert op.coeffs is not None
+
+    return RealspaceSum(modes, realspace_ops)
+
+
+def _validate_input(
+    states: int, modes: int, freqs: np.ndarray, taylor_coeffs: Sequence[np.ndarray]
+) -> None:
+    """Validate that the shapes of the harmonic frequencies and the Taylor coefficients are correct."""
+    for i, phi in enumerate(taylor_coeffs):
+        shape = (states, states) + (modes,) * i
+
+        if phi.shape != shape:
+            raise ValueError(
+                f"{i}th order coefficient tensor must have shape {shape}, got shape {phi.shape}"
+            )
+
+    if freqs.shape != (modes,):
+        raise TypeError(f"Frequencies must have shape {(modes,)}, got shape {freqs.shape}.")
diff --git a/pennylane/labs/trotter_error/product_formulas/__init__.py b/pennylane/labs/trotter_error/product_formulas/__init__.py
new file mode 100644
index 000000000..3e9a0b097
--- /dev/null
+++ b/pennylane/labs/trotter_error/product_formulas/__init__.py
@@ -0,0 +1,16 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Module containing functions for computing product formulas"""
+
+from .error import trotter_error, perturbation_error
diff --git a/pennylane/labs/trotter_error/product_formulas/error.py b/pennylane/labs/trotter_error/product_formulas/error.py
new file mode 100644
index 000000000..cf31b1284
--- /dev/null
+++ b/pennylane/labs/trotter_error/product_formulas/error.py
@@ -0,0 +1,133 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Functions for retreiving effective error from fragments"""
+
+from typing import List, Sequence
+
+from pennylane.labs.trotter_error import AbstractState, Fragment
+from pennylane.labs.trotter_error.abstract import nested_commutator
+
+
+class _AdditiveIdentity:
+    """Only used to initialize accumulators for summing Fragments"""
+
+    def __add__(self, other):
+        return other
+
+    def __radd__(self, other):
+        return other
+
+
+def trotter_error(fragments: Sequence[Fragment], delta: float) -> Fragment:
+    r"""Compute the second-order Trotter error operator.
+
+    For a Hamiltonian :math:`H` expressed as a sum of
+    fragments :math:`\sum_{m=1}^L H_m`, the second order Trotter formula is given by
+
+    .. math:: e^{iH\Delta t} \approx \prod_{m=1}^L e^{iH_m\Delta t / 2} \prod_{m=L}^1 e^{iH_m \Delta t / 2} = e^{i \tilde{H} \Delta t},
+
+    where :math:`\tilde{H} = H + \hat{\epsilon}`. The leading term of the error operator :math:`\hat{\epsilon}` is given by
+
+    .. math:: \hat{\epsilon} = \frac{- \Delta t^2}{24} \sum_{i=1}^{L-1} \sum_{j = i + 1}^L \left[ H_i + 2 \sum_{k = j + 1}^L H_k, \left[ H_i, H_j \right] \right].
+
+    Args:
+        fragments (Sequence[Fragments]): the set of :class:`~.pennylane.labs.trotter_error.Fragment`
+            objects to compute Trotter error from
+        delta (float): time step for the trotter formula.
+
+    Returns:
+        Fragment: the Trotter error operator
+
+    **Example**
+
+    >>> import numpy as np
+    >>> from pennylane.labs.trotter_error.fragments import vibrational_fragments
+    >>> from pennylane.labs.trotter_error.product_formulas import trotter_error
+    >>> n_modes = 4
+    >>> r_state = np.random.RandomState(42)
+    >>> freqs = r_state.random(4)
+    >>> taylor_coeffs = [
+    >>>     np.array(0),
+    >>>     r_state.random(size=(n_modes, )),
+    >>>     r_state.random(size=(n_modes, n_modes)),
+    >>>     r_state.random(size=(n_modes, n_modes, n_modes))
+    >>> ]
+    >>> frags = vibrational_fragments(n_modes, freqs, taylor_coeffs)
+    >>> delta = 0.001
+    >>> type(trotter_error(frags, delta))
+    <class 'pennylane.labs.trotter_error.realspace.realspace_operator.RealspaceSum'>
+
+    """
+
+    if len(fragments) == 0:
+        return fragments
+
+    eff = _AdditiveIdentity()
+    n_frags = len(fragments)
+    scalar = -(delta**2) / 24
+
+    for i in range(n_frags - 1):
+        for j in range(i + 1, n_frags):
+            eff += nested_commutator([fragments[i], fragments[i], fragments[j]])
+            for k in range(i + 1, n_frags):
+                eff += 2 * nested_commutator([fragments[k], fragments[i], fragments[j]])
+
+    eff = scalar * eff
+    return eff
+
+
+def perturbation_error(
+    fragments: Sequence[Fragment], states: Sequence[AbstractState], delta: float = 1
+) -> List[float]:
+    r"""Computes the perturbation theory error using the second-order Trotter error operator.
+
+    The second-order Trotter error operator, :math:`\hat{\epsilon}`, is given by the expression
+
+    .. math:: \hat{\epsilon} = \frac{- \Delta t^2}{24} \sum_{i=1}^{L-1} \sum_{j = i + 1}^L \left[ H_i + 2 \sum_{k = j + 1}^L H_k, \left[ H_i, H_j \right] \right].
+
+    For a state :math:`\left| \psi \right\rangle` the perturbation theory error is given by the expectation value :math:`\left\langle \psi \right| \hat{\epsilon} \left| \psi \right\rangle`.
+
+    Args:
+        fragments (Sequence[Fragments]): the set of :class:`~.pennylane.labs.trotter_error.Fragment`
+            objects to compute Trotter error from
+        states: (Sequence[AbstractState]): the states to compute expectation values from
+        delta (float): time step for the trotter error operator.
+
+    Returns:
+        List[float]: the list of expectation values computed from the Trotter error operator and the input states
+
+    **Example**
+
+    >>> from pennylane.labs.trotter_error import HOState, vibrational_fragments, perturbation_error
+    >>> import numpy as np
+    >>> n_modes = 2
+    >>> r_state = np.random.RandomState(42)
+    >>> freqs = r_state.random(n_modes)
+    >>> taylor_coeffs = [
+    >>>     np.array(0),
+    >>>     r_state.random(size=(n_modes, )),
+    >>>     r_state.random(size=(n_modes, n_modes)),
+    >>>     r_state.random(size=(n_modes, n_modes, n_modes))
+    >>> ]
+    >>> frags = vibrational_fragments(n_modes, freqs, taylor_coeffs)
+    >>> gridpoints = 5
+    >>> state1 = HOState(n_modes, gridpoints, {(0, 0): 1})
+    >>> state2 = HOState(n_modes, gridpoints, {(1, 1): 1})
+    >>> perturbation_error(frags, [state1, state2])
+    [(-0.9189251160920879+0j), (-4.797716682426851+0j)]
+    """
+
+    error = trotter_error(fragments, delta)
+
+    return [error.expectation(state, state) for state in states]
diff --git a/pennylane/labs/trotter_error/realspace/__init__.py b/pennylane/labs/trotter_error/realspace/__init__.py
new file mode 100644
index 000000000..7cc74ce4f
--- /dev/null
+++ b/pennylane/labs/trotter_error/realspace/__init__.py
@@ -0,0 +1,19 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Module for realspace operators"""
+
+from .ho_state import HOState, VibronicHO
+from .realspace_coefficients import RealspaceCoeffs
+from .realspace_operator import RealspaceOperator, RealspaceSum
+from .realspace_matrix import RealspaceMatrix
diff --git a/pennylane/labs/trotter_error/realspace/ho_state.py b/pennylane/labs/trotter_error/realspace/ho_state.py
new file mode 100644
index 000000000..a2933d9de
--- /dev/null
+++ b/pennylane/labs/trotter_error/realspace/ho_state.py
@@ -0,0 +1,337 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Contains the HOState class which represents a wavefunction in the harmonic oscillator basis"""
+
+from __future__ import annotations
+
+from typing import Dict, Sequence, Tuple
+
+import numpy as np
+from scipy.sparse import csr_array, identity, kron
+
+
+class HOState:
+    """Represent a wavefunction in the harmonic oscillator basis.
+
+    Args:
+        modes (int): the number of vibrational modes
+        gridpoints (int): the number of gridpoints used to discretize the state
+        state: (Union[scipy.sparse.csr_array, Dict[Tuple[int], float]]): a sparse state vector for the full wavefunction or a dictionary containing the interacting modes and their non-zero coefficients
+
+
+    **Examples**
+
+    Building an :class:`~.pennylane.labs.trotter_error.HOState` from a dictionary
+
+    >>> from pennylane.labs.trotter_error import HOState
+    >>> n_modes = 3
+    >>> gridpoints = 5
+    >>> state_dict = {(1, 2, 3): 1, (0, 3, 2): 1}
+    >>> HOState(n_modes, gridpoints, state_dict)
+    HOState(modes=3, gridpoints=5, <Compressed Sparse Row sparse array of dtype 'int64'
+        with 2 stored elements and shape (125, 1)>
+      Coords	Values
+      (17, 0)	1
+      (38, 0)	1)
+
+    Building an :class:`~.pennylane.labs.trotter_error.HOState` from a ``scipy.sparse.csr_array``
+
+    >>> from scipy.sparse import csr_array
+    >>> import numpy as np
+    >>> gridpoints = 2
+    >>> n_modes = 2
+    >>> state_vector = csr_array(np.array([0, 1, 0, 0]))
+    >>> HOState(n_modes, gridpoints, state_vector)
+    HOState(modes=2, gridpoints=2, <COOrdinate sparse array of dtype 'int64'
+        with 1 stored elements and shape (4, 1)>
+      Coords	Values
+      (1, 0)	1)
+    """
+
+    def __init__(self, modes: int, gridpoints: int, state: csr_array):
+        if isinstance(state, csr_array):
+            if state.shape == (gridpoints**modes,):
+                state = state.reshape((gridpoints**modes, 1))
+
+            if state.shape != (gridpoints**modes, 1):
+                raise ValueError(
+                    f"Dimension mismatch. Expected vector of shape {(gridpoints ** modes, 1)} but got shape {state.shape}."
+                )
+
+            self.vector = state
+        elif isinstance(state, dict):
+            self.vector = self._vector_from_dict(modes, gridpoints, state)
+        else:
+            raise TypeError(f"state must be of type csr_array or dict, got {type(state)}.")
+
+        self.gridpoints = gridpoints
+        self.modes = modes
+        self.dim = gridpoints**modes
+
+    @classmethod
+    def _vector_from_dict(
+        cls, modes: int, gridpoints: int, coeffs: Dict[Tuple[int], float]
+    ) -> csr_array:
+        """Construct an :class:`~.pennylane.labs.trotter_error.HOState` from a dictionary.
+
+        Args:
+            modes (int): the number of vibrational modes
+            gridpoints (int): the number of gridpoints used to discretize the state
+            coeffs (Dict[Tuple[int]]): a dictionary representation of the state
+
+        Returns:
+            HOState: an :class:`~.pennylane.labs.trotter_error.HOState` representation of the state vector
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import HOState
+        >>> n_modes = 3
+        >>> gridpoints = 5
+        >>> state_dict = {(1, 2, 3): 1, (0, 3, 2): 1}
+        >>> HOState.from_dict(n_modes, gridpoints, state_dict)
+        HOState(modes=3, gridpoints=5, <Compressed Sparse Row sparse array of dtype 'int64'
+            with 2 stored elements and shape (125, 1)>
+          Coords	Values
+          (17, 0)	1
+          (38, 0)	1)
+        """
+        rows, cols, vals = [], [], []
+        for index, val in coeffs.items():
+            if len(index) != modes:
+                raise ValueError(
+                    f"Number of modes given was {modes}, but {index} contains {len(index)} modes."
+                )
+
+            row = sum(i * (gridpoints**exp) for exp, i in enumerate(reversed(index)))
+            rows.append(row)
+            vals.append(val)
+
+        rows = np.array(rows)
+        cols = np.zeros_like(rows)
+        vals = np.array(vals)
+        return csr_array((vals, (rows, cols)), shape=(gridpoints**modes, 1))
+
+    @classmethod
+    def zero_state(cls, modes: int, gridpoints: int) -> HOState:
+        """Construct an :class:`~.pennylane.labs.trotter_error.HOState` whose vector is zero.
+
+        Args:
+            modes (int): the number of vibrational modes
+            gridpoints(int): the number of gridpoints used to discretize the state
+
+        Returns:
+            HOState: an :class:`~.pennylane.labs.trotter_error.HOState` representing the zero state
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import HOState
+        >>> HOState.zero_state(5, 10)
+        HOState(modes=5, gridpoints=10, <Compressed Sparse Row sparse array of dtype 'float64'
+            with 0 stored elements and shape (100000, 1)>)
+        """
+        return cls(modes, gridpoints, csr_array((gridpoints**modes, 1)))
+
+    def __add__(self, other: HOState) -> HOState:
+        if not isinstance(other, HOState):
+            raise TypeError(f"Can only add HOState with another HOState, got {type(other)}.")
+
+        return HOState(self.modes, self.gridpoints, self.vector + other.vector)
+
+    def __mul__(self, scalar: float) -> HOState:
+        return HOState(self.modes, self.gridpoints, scalar * self.vector)
+
+    __rmul__ = __mul__
+
+    def dot(self, other: HOState) -> float:
+        """Return the dot product of two :class:`~.pennylane.labs.trotter_error.HOState` objects.
+
+        Args:
+            other (HOState): the state to take the dot product with
+
+        Returns:
+            float: the dot product of the two states
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import HOState
+        >>> n_modes = 3
+        >>> gridpoints = 5
+        >>> state_dict = {(1, 2, 3): 1, (0, 3, 2): 1}
+        >>> state1 = HOState(n_modes, gridpoints, state_dict)
+        >>> state1.dot(state1)
+        2
+        """
+        if self.dim != other.dim:
+            raise ValueError(
+                f"Dimension mismatch. Attempting to dot product vectors of dimension {self.dim} and {other.dim}."
+            )
+
+        return ((self.vector.T).dot(other.vector))[0, 0]
+
+    def __repr__(self):
+        return f"HOState(modes={self.modes}, gridpoints={self.gridpoints}, {self.vector})"
+
+
+class VibronicHO:
+    """Represent the tensor product of harmonic oscillator states.
+
+    Args:
+        states (int): the number of electronic states
+        modes (int): the number of vibrational modes
+        gridpoints (int): the number of gridpoints used to discretize the state
+        ho_states (Sequence[HOState]): a sequence of :class:`~.pennylane.labs.trotter_error.HOState` objects representing the harmonic oscillator states
+
+    **Example**
+
+    >>> from pennylane.labs.trotter_error import HOState, VibronicHO
+    >>> n_modes = 3
+    >>> n_states = 2
+    >>> gridpoints = 5
+    >>> state_dict = {(1, 2, 3): 1, (0, 3, 2): 1}
+    >>> state = HOState(n_modes, gridpoints, state_dict)
+    >>> VibronicHO(n_states, n_modes, gridpoints, [state, state])
+    VibronicHO([HOState(modes=3, gridpoints=5, <Compressed Sparse Row sparse array of dtype 'int64'
+        with 2 stored elements and shape (125, 1)>
+      Coords	Values
+      (17, 0)	1
+      (38, 0)	1), HOState(modes=3, gridpoints=5, <Compressed Sparse Row sparse array of dtype 'int64'
+        with 2 stored elements and shape (125, 1)>
+      Coords	Values
+      (17, 0)	1
+      (38, 0)	1)])
+    """
+
+    def __init__(self, states: int, modes: int, gridpoints: int, ho_states: Sequence[HOState]):
+
+        if len(ho_states) != states:
+            raise ValueError(
+                f"Got {len(ho_states)} harmonic oscillator states, but expected {states}."
+            )
+
+        for ho_state in ho_states:
+            if ho_state.modes != modes:
+                raise ValueError(
+                    f"Mode mismatch: given {modes} modes, but found an HOState on {ho_state.modes} modes."
+                )
+
+            if ho_state.gridpoints != gridpoints:
+                raise ValueError(
+                    f"Gridpoint mismatch: given {gridpoints} gridpoints, but found an HOState on {ho_state.gridpoints} gridpoints."
+                )
+
+        self.states = states
+        self.modes = modes
+        self.gridpoints = gridpoints
+        self.ho_states = ho_states
+
+    def __repr__(self):
+        return f"VibronicHO({self.ho_states})"
+
+    def __add__(self, other: VibronicHO) -> VibronicHO:
+        if self.states != other.states:
+            raise ValueError(
+                f"Cannot add VibronicHO on {self.states} states with VibronicHO on {other.states} states."
+            )
+
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot add VibronicHO on {self.modes} modes with VibronicHO on {other.modes} modes."
+            )
+
+        if self.gridpoints != other.gridpoints:
+            raise ValueError(
+                f"Cannot add VibronicHO on {self.gridpoints} gridpoints with VibronicHO on {other.gridpoints} gridpoints."
+            )
+
+        return VibronicHO(
+            states=self.states,
+            modes=self.modes,
+            gridpoints=self.gridpoints,
+            ho_states=[x + y for x, y in zip(self.ho_states, other.ho_states)],
+        )
+
+    def __mul__(self, scalar: float) -> VibronicHO:
+        return VibronicHO(
+            states=self.states,
+            modes=self.modes,
+            gridpoints=self.gridpoints,
+            ho_states=[scalar * ho for ho in self.ho_states],
+        )
+
+    __rmul__ = __mul__
+
+    @classmethod
+    def zero_state(cls, states: int, modes: int, gridpoints: int) -> VibronicHO:
+        """Construct a :class:`~.pennylane.labs.trotter_error.VibronicHO` representing the zero state.
+
+        Args:
+            states (int): the number of electronic states
+            modes (int): the number of vibrational modes
+            gridpoints(int): the number of gridpoints used to discretize the state
+
+        Returns:
+            VibronicHO: a :class:`~.pennylane.labs.trotter_error.VibronicHO` representing the zero state
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import VibronicHO
+        >>> VibronicHO.zero_state(2, 3, 5)
+        VibronicHO([HOState(modes=3, gridpoints=5, <Compressed Sparse Row sparse array of dtype 'float64'
+            with 0 stored elements and shape (125, 1)>), HOState(modes=3, gridpoints=5, <Compressed Sparse Row sparse array of dtype 'float64'
+            with 0 stored elements and shape (125, 1)>)])
+        """
+        return cls(
+            states=states,
+            modes=modes,
+            gridpoints=gridpoints,
+            ho_states=[HOState.zero_state(modes, gridpoints)] * states,
+        )
+
+    def dot(self, other: VibronicHO):
+        """Return the dot product of two :class:`~.pennylane.labs.trotter_error.VibronicHO` objects.
+
+        Args:
+            other (VibronicHO): the state to take the dot product with
+
+        Returns:
+            float: the dot product of the two states
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import HOState, VibronicHO
+        >>> n_modes = 3
+        >>> n_states = 2
+        >>> gridpoints = 5
+        >>> state_dict = {(1, 2, 3): 1, (0, 3, 2): 1}
+        >>> state = HOState(n_modes, gridpoints, state_dict)
+        >>> vo_state = VibronicHO(n_states, n_modes, gridpoints, [state, state])
+        >>> vo_state.dot(vo_state)
+        4
+        """
+
+        return np.real(sum(x.dot(y) for x, y in zip(self.ho_states, other.ho_states)))
+
+
+def _tensor_with_identity(op: csr_array, gridpoints: int, n_modes: int, mode: int) -> csr_array:
+    """Return the tensor product of ``op`` with the ``gridpoints**mode`` dimensional identity matrix on the left, and the ``gridpoints ** (n_modes - mode - 1)`` dimensional identity matrix on the right."""
+    if mode == 0:
+        return kron(op, identity(gridpoints ** (n_modes - 1)))
+
+    if mode == n_modes - 1:
+        return kron(identity(gridpoints ** (n_modes - 1)), op)
+
+    id_left = identity(gridpoints**mode)
+    id_right = identity(gridpoints ** (n_modes - mode - 1))
+
+    return kron(id_left, kron(op, id_right))
diff --git a/pennylane/labs/trotter_error/realspace/matrix.py b/pennylane/labs/trotter_error/realspace/matrix.py
new file mode 100644
index 000000000..11b10327b
--- /dev/null
+++ b/pennylane/labs/trotter_error/realspace/matrix.py
@@ -0,0 +1,163 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Private helper functions for converting RealspaceOperator objects to matrices"""
+
+from __future__ import annotations
+
+from functools import reduce
+from typing import TYPE_CHECKING, Tuple, Union
+
+import numpy as np
+import scipy as sp
+
+if TYPE_CHECKING:
+    from pennylane.labs.trotter_error.realspace import RealspaceOperator, RealspaceSum
+
+
+def _position_operator(
+    gridpoints: int, sparse: bool = False, basis: str = "realspace"
+) -> Union[np.ndarray, sp.sparse.csr_array]:
+    """Returns a matrix representation of the position operator"""
+
+    if basis == "realspace":
+        matrix = _realspace_position(gridpoints)
+        return sp.sparse.csr_array(matrix) if sparse else matrix
+
+    if basis == "harmonic":
+        matrix = _harmonic_position(gridpoints)
+        return sp.sparse.csr_array(matrix) if sparse else matrix
+
+    raise ValueError(f'"{basis}" is not a valid basis')
+
+
+def _realspace_position(gridpoints: int) -> np.ndarray:
+    """Return a matrix representation of the position operator in the realspace basis"""
+    values = (np.arange(gridpoints) - gridpoints / 2) * (np.sqrt(2 * np.pi / gridpoints))
+    return np.diag(values)
+
+
+def _harmonic_position(gridpoints: int) -> np.ndarray:
+    """Return a matrix representation of the position operator in the harmonic basis"""
+    rows = np.array(list(range(1, gridpoints)) + list(range(0, gridpoints - 1)))
+    cols = np.array(list(range(0, gridpoints - 1)) + list(range(1, gridpoints)))
+    vals = np.array([np.sqrt(i) for i in range(1, gridpoints)] * 2)
+
+    matrix = np.zeros(shape=(gridpoints, gridpoints))
+    matrix[rows, cols] = vals
+
+    return 1 / np.sqrt(2) * matrix
+
+
+def _momentum_operator(
+    gridpoints: int, sparse: bool = False, basis: str = "realspace"
+) -> Union[np.ndarray, sp.sparse.csr_array]:
+    """Returns a matrix representation of the momentum operator"""
+
+    if basis == "realspace":
+        matrix = _realspace_momentum(gridpoints)
+        return sp.sparse.csr_array(matrix) if sparse else matrix
+
+    if basis == "harmonic":
+        matrix = _harmonic_momentum(gridpoints)
+        return sp.sparse.csr_array(matrix) if sparse else matrix
+
+    raise ValueError(f'"{basis}" is not a valid basis')
+
+
+def _realspace_momentum(gridpoints: int) -> np.ndarray:
+    """Returns a matrix representation of the momenumt operator in the realspace basis"""
+    values = np.arange(gridpoints)
+    values[gridpoints // 2 :] -= gridpoints
+    values = values * (np.sqrt(2 * np.pi / gridpoints))
+    dft = sp.linalg.dft(gridpoints, scale="sqrtn")
+    matrix = dft @ np.diag(values) @ dft.conj().T
+
+    return matrix
+
+
+def _harmonic_momentum(gridpoints: int) -> np.ndarray:
+    """Returns a matrix representation of the momenumt operator in the harmonic basis"""
+    rows = np.array(list(range(1, gridpoints)) + list(range(0, gridpoints - 1)))
+    cols = np.array(list(range(0, gridpoints - 1)) + list(range(1, gridpoints)))
+    vals = np.array(
+        [np.sqrt(i) for i in range(1, gridpoints)] + [-np.sqrt(i) for i in range(1, gridpoints)]
+    )
+
+    matrix = np.zeros(shape=(gridpoints, gridpoints))
+    matrix[rows, cols] = vals
+
+    return (1j / np.sqrt(2)) * matrix
+
+
+def _string_to_matrix(
+    op: str, gridpoints: int, sparse: bool = False, basis: str = "realspace"
+) -> Union[np.ndarray, sp.sparse.csr_array]:
+    """Transforms a string of P's and Q's into a matrix representing the product of position and momenutm operators"""
+
+    matrix = _identity(gridpoints, sparse=sparse)
+    p = _momentum_operator(gridpoints, basis=basis, sparse=sparse)
+    q = _position_operator(gridpoints, basis=basis, sparse=sparse)
+    mat_dict = {"P": p, "Q": q}
+
+    return reduce(lambda x, y: x @ mat_dict[y], op, matrix)
+
+
+def _tensor_with_identity(
+    modes: int,
+    gridpoints: int,
+    index: Tuple[int],
+    ops: Tuple[Union[np.ndarray, sp.sparse.csr_array]],
+    sparse: bool = False,
+) -> Union[np.ndarray, sp.sparse.csr_array]:
+    """Tensor the input matrices with the identity"""
+    lookup = [_identity(gridpoints, sparse=sparse)] * modes
+
+    for mode in range(modes):
+        for count, i in enumerate(index):
+            if i == mode:
+                lookup[mode] = lookup[mode] @ ops[count]
+
+    return reduce(_kron, lookup[1:], lookup[0])
+
+
+def _identity(dim: int, sparse: bool) -> Union[np.ndarray, sp.sparse.csr_array]:
+    """Return a matrix representation of the identity operator"""
+    if sparse:
+        return sp.sparse.eye_array(dim, format="csr")
+
+    return np.eye(dim)
+
+
+def _kron(a: Union[np.ndarray, sp.sparse.csr_array], b: Union[np.ndarray, sp.sparse.csr_array]):
+    """Return the Kronecker product of two matrices"""
+    if isinstance(a, np.ndarray) and isinstance(b, np.ndarray):
+        return np.kron(a, b)
+
+    if isinstance(a, sp.sparse.csr_array) and isinstance(b, sp.sparse.csr_array):
+        return sp.sparse.kron(a, b, format="csr")
+
+    raise TypeError(f"Matrices must both be ndarray or csr_array. Got {type(a)} and {type(b)}.")
+
+
+def _zeros(shape: Tuple[int], sparse: bool = False) -> Union[np.ndarray, sp.sparse.csr_array]:
+    """Return a matrix representation of the zero operator"""
+    if sparse:
+        return sp.sparse.csr_array(shape)
+
+    return np.zeros(shape)
+
+
+def _op_norm(gridpoints: int) -> float:
+    """The norm of P and Q"""
+    return np.sqrt(gridpoints * np.pi / 2)
diff --git a/pennylane/labs/trotter_error/realspace/realspace_coefficients.py b/pennylane/labs/trotter_error/realspace/realspace_coefficients.py
new file mode 100644
index 000000000..3fad86054
--- /dev/null
+++ b/pennylane/labs/trotter_error/realspace/realspace_coefficients.py
@@ -0,0 +1,620 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tree representation of coefficients of a realspace operator"""
+
+from __future__ import annotations
+
+from enum import Enum
+from itertools import product
+from typing import Dict, Tuple
+
+import numpy as np
+from numpy import allclose, isclose, ndarray, zeros
+
+
+class RealspaceCoeffs:
+    """Lightweight representation of a tensor of coefficients.
+
+    The :class:`~.pennylane.labs.trotter_error.RealspaceCoeffs` object is initialized with an array
+    and can be used to represent coefficients of a real space operator. A real space operator
+    is constrcuted from position and momentum operators, e.g., Eq. 4
+    of `arXiv:1703.09313 <https://arxiv.org/abs/1703.09313>`_ which represents a vibrational
+    Hamiltonian.
+
+    Args:
+        tensor (ndarray): a numpy tensor
+        label (string): name of the tensor
+
+    **Examples**
+
+    >>> import numpy as np
+    >>> from pennylane.labs.trotter_error import RealspaceCoeffs
+    >>> coeffs = np.array([[1, 0], [0, 1]])
+    >>> rs_coeffs = RealspaceCoeffs(coeffs, label="alpha")
+    >>> rs_coeffs.shape
+    (2, 2)
+
+    .. details::
+         :title: Usage Details
+
+         The :class:`~.pennylane.labs.trotter_error.RealspaceCoeffs` object allows arithmetic
+         operations such as addition, subtraction, multiplication and matrix multiplication.
+         Printing the resulting objects displays the expression that is used to compute each entry
+         of the tensor.
+
+         >>> coeffs1 = RealspaceCoeffs(np.array([[1, 0], [0, 1]]), label="alpha")
+         >>> coeffs2 = RealspaceCoeffs(np.array([[2, 1], [1, 3]]), label="beta")
+         >>> expr1 = coeffs1 + 2 * coeffs2
+         >>> coeffs3 = RealspaceCoeffs(np.array([3, 2]), label="omega")
+         >>> expr2 = expr1 @ coeffs3
+         >>> expr2
+         ((alpha[idx0,idx1]) + (2 * (beta[idx0,idx1]))) * (omega[idx2])
+    """
+
+    def __init__(self, tensor: np.ndarray, label: str = None):
+        if label is None and tensor.shape != ():
+            raise ValueError(f"A label is required for a tensor of shape {tensor.shape}.")
+
+        self._tree = _RealspaceTree.tensor_node(tensor, label)
+
+    @classmethod
+    def _from_tree(cls, tree: _RealspaceTree):
+        """Initialize directly from a ``_RealspaceTree`` object."""
+        rs_coeffs = cls.__new__(cls)
+        rs_coeffs._tree = tree  # pylint: disable=protected-access
+        return rs_coeffs
+
+    def __add__(self, other: RealspaceCoeffs) -> RealspaceCoeffs:
+        tree = _RealspaceTree.sum_node(self._tree, other._tree)
+        return RealspaceCoeffs._from_tree(tree)
+
+    def __sub__(self, other: RealspaceCoeffs) -> RealspaceCoeffs:
+        tree = _RealspaceTree.sum_node(self._tree, _RealspaceTree.scalar_node(-1, other._tree))
+        return RealspaceCoeffs._from_tree(tree)
+
+    def __mul__(self, scalar: float) -> RealspaceCoeffs:
+        tree = _RealspaceTree.scalar_node(scalar, self._tree)
+        return RealspaceCoeffs._from_tree(tree)
+
+    __rmul__ = __mul__
+
+    def __matmul__(self, other: RealspaceCoeffs) -> RealspaceCoeffs:
+        tree = _RealspaceTree.outer_node(self._tree, other._tree)
+        return RealspaceCoeffs._from_tree(tree)
+
+    def __repr__(self) -> str:
+        return self._tree.__repr__()
+
+    def __getitem__(self, index) -> float:
+        return self._tree.compute(index)
+
+    def __eq__(self, other: RealspaceCoeffs) -> bool:
+        return self._tree == other._tree
+
+    @property
+    def is_zero(self) -> bool:
+        """Determine if the :class:`~.pennylane.labs.trotter_error.RealspaceCoeffs` objects
+        represents the zero tensor.
+
+        Returns:
+            bool: returns ``True`` when the tensor is zero, otherwise returns ``False``
+        """
+        return self._tree.is_zero
+
+    @property
+    def shape(self) -> Tuple[int]:
+        """Return the shape of the tensor."""
+        return self._tree.shape
+
+    def nonzero(self, threshold: float = 0.0):
+        """Return the nonzero coefficients in a dictionary.
+
+        Args:
+            threshold (float): tolerance to return coefficients with magnitude greater than ``threshold``
+
+        Returns:
+            dict: a dictionary representation of the coefficient tensor
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceCoeffs
+        >>> import numpy as np
+        >>> node = RealspaceCoeffs(np.array([[1, 0, 0, 1], [0, 0, 1, 1]]), label="alpha")
+        >>> node.nonzero()
+        {(0, 0): 1, (0, 3): 1, (1, 2): 1, (1, 3): 1}
+        """
+
+        return self._tree.nonzero(threshold)
+
+
+class _NodeType(Enum):
+    """Enum containing the types of nodes"""
+
+    SUM = 1
+    OUTER = 2
+    SCALAR = 3
+    TENSOR = 4
+    FLOAT = 5
+
+
+class _RealspaceTree:  # pylint: disable=too-many-instance-attributes
+    """
+     A tree representing an expression that computes the coefficients of a :class:`~.pennylane.labs.trotter_error.RealspaceOperator`.
+     This class should be instantiated from the following class methods:
+
+        * ``tensor_node(tensor)``: a leaf node containing the coefficients as a tensor
+        * ``outer_node(l_child, r_child)``: a node representing the outer product of two ``RealspaceCoeffs`` objects
+        * ``sum_node(l_child, r_child)``: a node representing the sum of two ``RealspaceCoeffs`` objects
+        * ``scalar_node(scalar, child)``: a node representing the product of a ``RealspaceCoeffs`` object by a scalar
+
+    **Examples**
+
+    Building a node representing a tensor of coefficients
+
+    >>> from pennylane.labs.trotter_error import _RealspaceTree
+    >>> import numpy as np
+    >>> node = _RealspaceTree.tensor_node(np.array([[1, 2, 3], [4, 5, 6]]), label="alpha")
+    >>> node
+    alpha[idx0,idx1]
+
+    Building a node representing the outer product of its children
+
+    >>> from pennylane.labs.trotter_error import _RealspaceTree
+    >>> import numpy as np
+    >>> left_child = _RealspaceTree.tensor_node(np.array([1, 2, 3]), label="alpha")
+    >>> right_child = _RealspaceTree.tensor_node(np.array([[1, 3, 4], [4, 5, 6]]), label="beta")
+    >>> parent = _RealspaceTree.outer_node(left_child, right_child)
+    >>> parent
+    (alpha[idx0]) * (beta[idx1,idx2])
+
+    Building a node representing the sum of its children
+
+    >>> from pennylane.labs.trotter_error import _RealspaceTree
+    >>> import numpy as np
+    >>> left_child = _RealspaceTree.tensor_node(np.array([1, 2, 3]), label="alpha")
+    >>> right_child = _RealspaceTree.tensor_node(np.array([4, 5, 6]), label="beta")
+    >>> parent = _RealspaceTree.sum_node(left_child, right_child)
+    >>> parent
+    (alpha[idx0]) + (beta[idx0])
+
+    Building a node representing the multiplication of its child by a scalar
+
+    >>> from pennylane.labs.trotter_error import _RealspaceTree
+    >>> import numpy as np
+    >>> child = _RealspaceTree.tensor_node(np.array([[1, 2, 3], [4, 5, 6]]), label="alpha")
+    >>> parent = _RealspaceTree.scalar_node(5, child)
+    >>> parent
+    5 * (alpha[idx0,idx1])
+    """
+
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
+    def __init__(
+        self,
+        node_type: _NodeType,
+        l_child: _RealspaceTree = None,
+        r_child: _RealspaceTree = None,
+        tensor: ndarray = None,
+        scalar: float = None,
+        value: float = None,
+        label: str = None,
+    ) -> _RealspaceTree:
+
+        self.node_type = node_type
+        self.l_child = l_child
+        self.r_child = r_child
+        self.tensor = tensor
+        self.scalar = scalar
+        self.value = value
+        self.label = label
+
+        if node_type == _NodeType.SUM:
+            self.shape = l_child.shape
+            self.is_zero = l_child.is_zero and r_child.is_zero
+        elif node_type == _NodeType.OUTER:
+            self.shape = l_child.shape + r_child.shape
+            self.is_zero = l_child.is_zero or r_child.is_zero
+        elif node_type == _NodeType.SCALAR:
+            self.shape = l_child.shape
+            self.is_zero = l_child.is_zero or isclose(scalar, 0)
+        elif node_type == _NodeType.TENSOR:
+            self.shape = tensor.shape
+            self.is_zero = allclose(tensor, zeros(tensor.shape))
+        elif node_type == _NodeType.FLOAT:
+            self.shape = ()
+            self.is_zero = isclose(value, 0)
+        else:
+            raise ValueError(f"Got invalid node type {node_type}.")
+
+    @classmethod
+    def sum_node(cls, l_child: _RealspaceTree, r_child: _RealspaceTree) -> _RealspaceTree:
+        """Returns a ``_RealspaceTree`` representing the sum of the two children nodes.
+
+        Args:
+            l_child (_RealspaceTree): the left child
+            r_child (_RealspaceTree): the right child
+
+        Returns:
+            _RealspaceTree: a `RealspaceCoeff` object representing the sum of `l_child` and `r_child`
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import _RealspaceTree
+        >>> import numpy as np
+        >>> left_child = _RealspaceTree.tensor_node(np.array([1, 2, 3]), label="alpha")
+        >>> right_child = _RealspaceTree.tensor_node(np.array([4, 5, 6]), label="beta")
+        >>> parent = _RealspaceTree.sum_node(left_child, right_child)
+        >>> parent
+        >>> (alpha[idx0]) + (beta[idx0])
+
+        """
+
+        if l_child.shape != r_child.shape:
+            raise ValueError(
+                f"Cannot add _RealspaceTree of shape {l_child.shape} with _RealspaceTree of shape {r_child.shape}."
+            )
+
+        return cls(
+            node_type=_NodeType.SUM,
+            l_child=l_child,
+            r_child=r_child,
+        )
+
+    @classmethod
+    def outer_node(
+        cls,
+        l_child: _RealspaceTree,
+        r_child: _RealspaceTree,
+    ) -> _RealspaceTree:
+        """Returns a ``_RealspaceTree`` representing the outer product of the two children nodes.
+
+        Args:
+            l_child (_RealspaceTree): the left child
+            r_child (RealspaceCOeffs): the right child
+
+        Returns:
+            _RealspaceTree: a ``_RealspaceTree`` object representing the outer product of ``l_child`` and ``r_child``
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import _RealspaceTree
+        >>> import numpy as np
+        >>> left_child = _RealspaceTree.tensor_node(np.array([1, 2, 3]), label="alpha")
+        >>> right_child = _RealspaceTree.tensor_node(np.array([[1, 3, 4], [4, 5, 6]]), label="beta")
+        >>> parent = _RealspaceTree.outer_node(left_child, right_child)
+        >>> parent
+        (alpha[idx0]) * (beta[idx1,idx2])
+        """
+
+        return cls(
+            node_type=_NodeType.OUTER,
+            l_child=l_child,
+            r_child=r_child,
+        )
+
+    @classmethod
+    def tensor_node(cls, tensor: ndarray, label: str = None) -> _RealspaceTree:
+        """Returns a ``_RealspaceTree`` leaf node storing a tensor of coefficients.
+
+        Args:
+            tensor (ndarray): a tensor of coefficients
+            label (string): a label for the tensor to be used when displaying the ``_RealspaceTree`` object as an expression
+
+        Returns:
+            _RealspaceTree: a ``_RealspaceTree`` object representing containing the tensor
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import _RealspaceTree
+        >>> import numpy as np
+        >>> node = _RealspaceTree.tensor_node(np.array([[1, 2, 3], [4, 5, 6]]), label="alpha")
+        >>> node
+        alpha[idx0,idx1]
+        """
+
+        if len(tensor.shape):
+            return cls(
+                node_type=_NodeType.TENSOR,
+                tensor=tensor,
+                label=label,
+            )
+
+        return cls(
+            node_type=_NodeType.FLOAT,
+            value=tensor,
+        )
+
+    @classmethod
+    def scalar_node(cls, scalar: float, child: _RealspaceTree) -> _RealspaceTree:
+        """Returns a ``_RealspaceTree`` representing the scalar product of ``scalar`` and ``child``.
+
+        Args:
+            scalar (float): a scalar to multiply ``child`` by
+            child (_RealspaceTree): the ``_RealspaceTree`` object to be multiplied by ``scalar``
+
+        Returns:
+            _RealspaceTree: a ``_RealspaceTree`` object representing the coefficients of ``child`` multiplied by ``scalar``
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import _RealspaceTree
+        >>> import numpy as np
+        >>> child = _RealspaceTree.tensor_node(np.array([[1, 2, 3], [4, 5, 6]]), label="alpha")
+        >>> parent = _RealspaceTree.scalar_node(5, child)
+        >>> parent
+        5 * (alpha[idx0,idx1])
+        """
+
+        return cls(
+            node_type=_NodeType.SCALAR,
+            l_child=child,
+            scalar=scalar,
+        )
+
+    def __add__(self, other: _RealspaceTree) -> _RealspaceTree:
+        return self.__class__.sum_node(self, other)
+
+    def __mul__(self, scalar: float) -> _RealspaceTree:
+        return self.__class__.scalar_node(scalar, self)
+
+    __rmul__ = __mul__
+
+    def __matmul__(self, other: _RealspaceTree) -> _RealspaceTree:
+        return self.__class__.outer_node(self, other)
+
+    def __eq__(self, other: _RealspaceTree) -> bool:  # pylint: disable=too-many-return-statements
+        if self.node_type != other.node_type:
+            return False
+
+        if self.shape != other.shape:
+            return False
+
+        if self.node_type == _NodeType.OUTER:
+            if self.l_child.shape != other.l_child.shape:
+                return False
+
+            if self.r_child.shape != other.r_child.shape:
+                return False
+
+            return (self.l_child == other.l_child) and (self.r_child == other.r_child)
+
+        if self.node_type == _NodeType.SCALAR:
+            if self.scalar != other.scalar:
+                return False
+
+            return self.l_child == other.l_child
+
+        if self.node_type == _NodeType.SUM:
+            return (self.l_child == other.l_child) and (self.r_child == other.r_child)
+
+        if self.node_type == _NodeType.TENSOR:
+            return allclose(self.tensor, other.tensor)
+
+        if self.node_type == _NodeType.FLOAT:
+            return self.value == other.value
+
+        raise ValueError(f"_RealspaceTree was constructed with invalid _NodeType {self.node_type}.")
+
+    def __repr__(self) -> str:
+        return self.__str__()
+
+    def __str__(self) -> str:
+        indices = [f"idx{i}" for i in range(len(self.shape))]
+
+        return str(self._str(indices))
+
+    # pylint: disable=protected-access
+    def _str(self, indices) -> str:
+
+        if self.node_type == _NodeType.TENSOR:
+            return f"{self.label}[{','.join(indices)}]"
+        if self.node_type == _NodeType.FLOAT:
+            return f"{self.value}"
+        if self.node_type == _NodeType.SCALAR:
+            return f"{self.scalar} * ({self.l_child._str(indices)})"
+        if self.node_type == _NodeType.OUTER:
+            l_indices = indices[: len(self.l_child.shape)]
+            r_indices = indices[len(self.l_child.shape) :]
+            return f"({self.l_child._str(l_indices)}) * ({self.r_child._str(r_indices)})"
+        if self.node_type == _NodeType.SUM:
+            return f"({self.l_child._str(indices)}) + ({self.r_child._str(indices)})"
+
+        raise ValueError(f"_RealspaceTree was constructed with invalid _NodeType {self.node_type}.")
+
+    def compute(self, index: Tuple[int]) -> float:
+        """Evaluate the tree on a given ``index``.
+
+        Args:
+            index (Tuple[int]): the index of the coefficient to be computed
+
+        Returns:
+            float: the coefficient at the given index
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import _RealspaceTree
+        >>> import numpy as np
+        >>> left_child = _RealspaceTree.tensor_node(np.array([1, 2, 3]), label="alpha")
+        >>> right_child = _RealspaceTree.tensor_node(np.array([[1, 3, 4], [4, 5, 6]]), label="beta")
+        >>> parent = _RealspaceTree.outer_node(left_child, right_child)
+        >>> parent.compute((1, 1, 2))
+        12
+        """
+
+        if not self._validate_index(index):
+            raise ValueError(f"Given index {index} is not compatible with shape {self.shape}")
+
+        if self.node_type == _NodeType.TENSOR:
+            return self.tensor[index]
+        if self.node_type == _NodeType.FLOAT:
+            return self.value
+        if self.node_type == _NodeType.SCALAR:
+            return self.scalar * self.l_child.compute(index)
+        if self.node_type == _NodeType.SUM:
+            return self.l_child.compute(index) + self.r_child.compute(index)
+        if self.node_type == _NodeType.OUTER:
+            l_index = index[: len(self.l_child.shape)]
+            r_index = index[len(self.l_child.shape) :]
+            return self.l_child.compute(l_index) * self.r_child.compute(r_index)
+
+        raise ValueError(f"_RealspaceTree was constructed with invalid _NodeType {self.node_type}.")
+
+    def _validate_index(self, index: Tuple[int]) -> bool:
+        """Validate the shape of an index.
+
+        Args:
+            index (Tuple[int]): an index
+
+        Returns:
+            bool: True if ``index`` corresponds to a valid index of the tensor, False otherwise
+        """
+
+        if len(index) != len(self.shape):
+            return False
+
+        for x, y in zip(index, self.shape):
+            if x < 0:
+                return False
+
+            if x >= y:
+                return False
+
+        return True
+
+    def nonzero(self, threshold: float = 0.0):
+        """Return the nonzero coefficients in a dictionary.
+
+        Args:
+            threshold (float): only return coefficients with magnitude greater than ``threshold``
+
+        Returns:
+            dict: a dictionary representation of the coefficient tensor
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import _RealspaceTree
+        >>> import numpy as np
+        >>> node = _RealspaceTree.tensor_node(np.array([[1, 0, 0, 1], [0, 0, 1, 1]]), label="alpha")
+        >>> node.nonzero()
+        {(0, 0): 1, (0, 3): 1, (1, 2): 1, (1, 3): 1}
+        """
+
+        if self.node_type == _NodeType.TENSOR:
+            return _numpy_to_dict(self.tensor, threshold)
+        if self.node_type == _NodeType.FLOAT:
+            return {(): self.value}
+        if self.node_type == _NodeType.SCALAR:
+            return _scale_dict(self.scalar, self.l_child.nonzero(threshold), threshold)
+        if self.node_type == _NodeType.SUM:
+            return _add_dicts(
+                self.l_child.nonzero(threshold),
+                self.r_child.nonzero(threshold),
+                threshold,
+            )
+        if self.node_type == _NodeType.OUTER:
+            return _mul_dicts(
+                self.l_child.nonzero(threshold),
+                self.r_child.nonzero(threshold),
+                threshold,
+            )
+
+        raise ValueError(f"_RealspaceTree was constructed with invalid _NodeType {self.node_type}.")
+
+
+def _add_dicts(d1: Dict, d2: Dict, threshold: float):
+    """Add two coefficient dictionaries
+
+    Args:
+        d1 (dict): the first dictionary to be added
+        d2 (dict): the second dictionary to be added
+        threshold (float): only return coefficients with magnitude greater than ``threshold``
+
+    Returns:
+        dict: the sum of ``d1`` and ``d2``
+    """
+    add_dict = {}
+
+    d1_keys = set(d1.keys())
+    d2_keys = set(d2.keys())
+
+    for key in d1_keys.intersection(d2_keys):
+        if abs(d1[key] + d2[key]) > threshold:
+            add_dict[key] = d1[key] + d2[key]
+
+    for key in d1_keys.difference(d2_keys):
+        add_dict[key] = d1[key]
+
+    for key in d2_keys.difference(d1_keys):
+        add_dict[key] = d2[key]
+
+    return add_dict
+
+
+def _mul_dicts(d1, d2, threshold):
+    """Multiply two coefficient dictionaries
+
+    Args:
+        d1 (dict): the first dictionary to be multiplied
+        d2 (dict): the second dictionary to be multiplied
+        threshold (float): only return coefficients with magnitude greater than ``threshold``
+
+    Returns:
+        dict: the outer product of ``d1`` and ``d2``
+    """
+    mul_dict = {}
+
+    for key1, key2 in product(d1.keys(), d2.keys()):
+        if abs(d1[key1] * d2[key2]) > threshold:
+            mul_dict[key1 + key2] = d1[key1] * d2[key2]
+
+    return mul_dict
+
+
+def _scale_dict(scalar, d, threshold):
+    """Multiply a coefficient dictionary by a scalar
+
+    Args:
+        d (dict): the dictionary to be scaled
+        scalar (float): the scalar to multiply ``d`` by
+        threshold (float): only return coefficients with magnitude greater than ``threshold``
+
+    Returns:
+        dict: the product of ``d`` and ``scalar``
+    """
+    scaled = {}
+    for key in d.keys():
+        if abs(scalar * d[key]) > threshold:
+            scaled[key] = scalar * d[key]
+
+    return scaled
+
+
+def _numpy_to_dict(arr, threshold):
+    """Returns a dictionary representation of a numpy array
+
+    Args:
+        arr (ndarray): a numpy array
+        threshold (float): only return coefficients with magnitude greater than ``threshold``
+
+    Returns:
+        dict: a dictionary representation of the numpy array
+    """
+    nz = arr.nonzero()
+    d = {}
+
+    for index in zip(*nz):
+        if abs(arr[index]) > threshold:
+            index = tuple(map(int, index))
+            d[index] = float(arr[index])
+
+    return d
diff --git a/pennylane/labs/trotter_error/realspace/realspace_matrix.py b/pennylane/labs/trotter_error/realspace/realspace_matrix.py
new file mode 100644
index 000000000..fb394ef41
--- /dev/null
+++ b/pennylane/labs/trotter_error/realspace/realspace_matrix.py
@@ -0,0 +1,491 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""The RealspaceMatrix class"""
+
+from __future__ import annotations
+
+import math
+from itertools import product
+from typing import Dict, Tuple, Union
+
+import numpy as np
+import scipy as sp
+
+from pennylane.labs.trotter_error import Fragment
+from pennylane.labs.trotter_error.realspace import HOState, RealspaceSum, VibronicHO
+from pennylane.labs.trotter_error.realspace.matrix import _kron, _zeros
+
+# pylint: disable=protected-access
+
+
+class RealspaceMatrix(Fragment):
+    r"""Implements a dictionary of :class:`~.pennylane.labs.trotter_error.RealspaceSum` objects.
+
+    This can be used to represent the fragments of a vibronic Hamiltonian given by, Eq. 3
+    of `arXiv:2411.13669 <https://arxiv.org/abs/2411.13669v1>`_,
+
+    .. math:: V_{i,j} = \lambda_{i,j} + \sum_{r} \phi^{(1)}_{i,j,r} Q_r + \sum_{r,s} \phi^{(2)}_{i,j,r,s} Q_r Q_s + \sum_{r,s,t} \phi^{(3)}_{i,j,r,s,t} Q_r Q_s Q_t + \dots,
+
+    where the dictionary is indexed by tuples :math:`(i, j)` and the values are
+    :class:`~.RealspaceSum` objects representing the operator :math:`V_{i,j}`.
+
+    Args:
+        states (int): the number of electronic states
+        modes (int): the number of vibrational modes
+        blocks (Dict[Tuple[int, int], RealspaceSum): a dictionary representation of the block matrix
+
+    **Example**
+
+    >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+    >>> import numpy as np
+    >>> n_states = 1
+    >>> n_modes = 5
+    >>> op1 = RealspaceOperator(n_modes, (), RealspaceCoeffs(np.array(1)))
+    >>> op2 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+    >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+    >>> RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum})
+    RealspaceMatrix({(0, 0): RealspaceSum((RealspaceOperator(5, (), 1), RealspaceOperator(5, 'Q', phi[idx0])))})
+    """
+
+    def __init__(
+        self,
+        states: int,
+        modes: int,
+        blocks: Dict[Tuple[int, int], RealspaceSum] = None,
+    ) -> RealspaceMatrix:
+
+        if blocks is None:
+            blocks = {}
+
+        self._blocks = blocks
+        self.states = states
+        self.modes = modes
+
+    def block(self, row: int, col: int) -> RealspaceSum:
+        """Return the :class:`~.pennylane.labs.trotter_error.RealspaceSum` object located at the
+        ``(row, col)`` entry of the :class:`~.pennylane.labs.trotter_error.RealspaceMatrix`.
+
+        Args:
+            row (int): the row of the index
+            col (int): the column of the index
+
+        Returns:
+            RealspaceSum: the :class:`~.pennylane.labs.trotter_error.RealspaceSum` object indexed
+                at ``(row, col)``
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+        >>> import numpy as np
+        >>> n_states = 1
+        >>> n_modes = 5
+        >>> op1 = RealspaceOperator(n_modes, (), RealspaceCoeffs(np.array(1)))
+        >>> op2 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+        >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+        >>> RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum}).block(0, 0)
+        RealspaceSum((RealspaceOperator(5, (), 1), RealspaceOperator(5, 'Q', phi[idx0])))
+        """
+        if row < 0 or col < 0:
+            raise IndexError(f"Index cannot be negative, got {(row, col)}.")
+        if row >= self.states or col >= self.states:
+            raise IndexError(
+                f"Index out of bounds. Got {(row, col)} but there are only {self.states} states."
+            )
+
+        return self._blocks.get((row, col), RealspaceSum.zero(self.modes))
+
+    def set_block(self, row: int, col: int, rs_sum: RealspaceSum) -> None:
+        """Set the value of the block indexed at ``(row, col)``.
+
+        Args:
+            row (int): the row of the index
+            col (int): the column of the index
+            rs_sum (RealspaceSum): the :class:`~.pennylane.labs.trotter_error.RealspaceSum` object to stored in index ``(row, col)``
+
+        Returns:
+            None
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+        >>> import numpy as np
+        >>> n_states = 2
+        >>> n_modes = 5
+        >>> op1 = RealspaceOperator(n_modes, (), RealspaceCoeffs(np.array(1)))
+        >>> op2 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+        >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+        >>> vib = RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum})
+        >>> vib
+        RealspaceMatrix({(0, 0): RealspaceSum((RealspaceOperator(5, (), 1), RealspaceOperator(5, 'Q', phi[idx0])))})
+        >>> vib.set_block(1, 1, rs_sum)
+        >>> vib
+        RealspaceMatrix({(0, 0): RealspaceSum((RealspaceOperator(5, (), 1), RealspaceOperator(5, 'Q', phi[idx0]))), (1, 1): RealspaceSum((RealspaceOperator(5, (), 1), RealspaceOperator(5, 'Q', phi[idx0])))})
+        """
+
+        if not isinstance(rs_sum, RealspaceSum):
+            raise TypeError(f"Block value must be RealspaceSum. Got {type(rs_sum)}.")
+        if row < 0 or col < 0:
+            raise IndexError(f"Index cannot be negative, got {(row, col)}.")
+        if row >= self.states or col >= self.states:
+            raise IndexError(
+                f"Index out of bounds. Got {(row, col)} but there are only {self.states} states."
+            )
+
+        if rs_sum._is_zero:
+            return
+
+        self._blocks[(row, col)] = rs_sum
+
+    def matrix(
+        self, gridpoints: int, sparse: bool = False, basis: str = "realspace"
+    ) -> Union[np.ndarray, sp.sparse.csr_matrix]:
+        """Return a matrix representation of the operator.
+
+        Args:
+            gridpoints (int): the number of gridpoints used to discretize the position or momentum operators
+            basis (str): the basis of the matrix, available options are ``realspace`` and ``harmonic``
+            sparse (bool): if ``True`` returns a sparse matrix, otherwise returns a dense matrix
+
+        Returns:
+            Union[ndarray, scipy.sparse.csr_array]: the matrix representation of the :class:`~.pennylane.labs.trotter_error.RealspaceOperator`
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+        >>> import numpy as np
+        >>> n_states = 1
+        >>> n_modes = 5
+        >>> op1 = RealspaceOperator(n_modes, (), RealspaceCoeffs(np.array(1)))
+        >>> op2 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+        >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+        >>> RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum}).matrix(2)
+        [[-25.58680776   0.           0.         ...   0.           0.
+            0.        ]
+         [  0.         -16.72453851   0.         ...   0.           0.
+            0.        ]
+         [  0.           0.         -18.49699236 ...   0.           0.
+            0.        ]
+         ...
+         [  0.           0.           0.         ...  -6.0898154    0.
+            0.        ]
+         [  0.           0.           0.         ...   0.          -7.86226925
+            0.        ]
+         [  0.           0.           0.         ...   0.           0.
+            1.        ]]
+        """
+        pow2 = _next_pow_2(self.states)
+        dim = pow2 * gridpoints**self.modes
+        shape = (pow2, pow2)
+        matrix = _zeros((dim, dim), sparse=sparse)
+
+        for index, rs_sum in self._blocks.items():
+            if sparse:
+                data = np.array([1])
+                indices = (np.array([index[0]]), np.array([index[1]]))
+                indicator = sp.sparse.csr_array((data, indices), shape=shape)
+            else:
+                indicator = np.zeros(shape=shape)
+                indicator[index] = 1
+
+            block = rs_sum.matrix(gridpoints, basis=basis, sparse=sparse)
+            matrix = matrix + _kron(indicator, block)
+
+        return matrix
+
+    def norm(self, params: Dict) -> float:
+        """Returns an upper bound on the spectral norm of the operator.
+
+        Args:
+            params (dict[str, Union[int, bool]]): The dictionary of parameters. The supported parameters are
+
+                * ``gridpoints`` (int): the number of gridpoints used to discretize the operator
+                * ``sparse`` (bool): If ``True``, use optimizations for sparse operators. Defaults to ``False``.
+
+        Returns:
+            float: an upper bound on the spectral norm of the operator
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+        >>> import numpy as np
+        >>> n_states = 1
+        >>> n_modes = 5
+        >>> op1 = RealspaceOperator(n_modes, (), RealspaceCoeffs(np.array(1)))
+        >>> op2 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+        >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+        >>> params = {"gridpoints": 2, "sparse": True}
+        >>> RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum}).norm(params)
+        27.586807763582737
+        """
+        try:
+            gridpoints = params["gridpoints"]
+        except KeyError as e:
+            raise KeyError("Need to specify the number of gridpoints") from e
+
+        if not _is_pow_2(gridpoints) or gridpoints <= 0:
+            raise ValueError(
+                f"Number of gridpoints must be a positive power of 2, got {gridpoints}."
+            )
+
+        padded = RealspaceMatrix(_next_pow_2(self.states), self.modes, self._blocks)
+
+        return padded._norm(params)
+
+    def _norm(self, params: Dict) -> float:
+        """Returns an upper bound on the spectral norm. This method assumes ``self.states`` is a power of two."""
+        if self.states == 1:
+            return self.block(0, 0).norm(params)
+
+        top_left, top_right, bottom_left, bottom_right = self._partition_into_quadrants()
+
+        norm1 = max(top_left._norm(params), bottom_right._norm(params))
+        norm2 = math.sqrt(top_right._norm(params) * bottom_left._norm(params))
+
+        return norm1 + norm2
+
+    def __add__(self, other: RealspaceMatrix) -> RealspaceMatrix:
+        if self.states != other.states:
+            raise ValueError(
+                f"Cannot add RealspaceMatrix on {self.states} states with RealspaceMatrix on {other.states} states."
+            )
+
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot add RealspaceMatrix on {self.modes} states with RealspaceMatrix on {other.modes} states."
+            )
+
+        new_blocks = {}
+        l_keys = set(self._blocks.keys())
+        r_keys = set(other._blocks.keys())
+
+        for key in l_keys.intersection(r_keys):
+            new_blocks[key] = self._blocks[key] + other._blocks[key]
+
+        for key in l_keys.difference(r_keys):
+            new_blocks[key] = self._blocks[key]
+
+        for key in r_keys.difference(l_keys):
+            new_blocks[key] = other._blocks[key]
+
+        return RealspaceMatrix(
+            self.states,
+            self.modes,
+            new_blocks,
+        )
+
+    def __sub__(self, other: RealspaceMatrix) -> RealspaceMatrix:
+        if self.states != other.states:
+            raise ValueError(
+                f"Cannot subtract RealspaceMatrix on {self.states} states with RealspaceMatrix on {other.states} states."
+            )
+
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot subtract RealspaceMatrix on {self.modes} states with RealspaceMatrix on {other.modes} states."
+            )
+
+        new_blocks = {}
+
+        l_keys = set(self._blocks.keys())
+        r_keys = set(other._blocks.keys())
+
+        for key in l_keys.intersection(r_keys):
+            new_blocks[key] = self._blocks[key] - other._blocks[key]
+
+        for key in l_keys.difference(r_keys):
+            new_blocks[key] = self._blocks[key]
+
+        for key in r_keys.difference(l_keys):
+            new_blocks[key] = (-1) * other._blocks[key]
+
+        return RealspaceMatrix(
+            self.states,
+            self.modes,
+            new_blocks,
+        )
+
+    def __mul__(self, scalar: float) -> RealspaceMatrix:
+        new_blocks = {}
+        for key in self._blocks.keys():
+            new_blocks[key] = scalar * self._blocks[key]
+
+        return RealspaceMatrix(self.states, self.modes, new_blocks)
+
+    __rmul__ = __mul__
+
+    def __imul__(self, scalar: float) -> RealspaceMatrix:
+        for key in self._blocks.keys():
+            self._blocks[key] *= scalar
+
+        return self
+
+    def __matmul__(self, other: RealspaceMatrix) -> RealspaceMatrix:
+        if self.states != other.states:
+            raise ValueError(
+                f"Cannot multiply RealspaceMatrix on {self.states} states with RealspaceMatrix on {other.states} states."
+            )
+
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot multiply RealspaceMatrix on {self.states} states with RealspaceMatrix on {other.states} states."
+            )
+
+        product_matrix = RealspaceMatrix(self.states, self.modes)
+
+        for i, j in product(range(self.states), repeat=2):
+            for op in self.block(i, j).ops:
+                assert op.coeffs is not None
+            for op in other.block(i, j).ops:
+                assert op.coeffs is not None
+            block_products = [self.block(i, k) @ other.block(k, j) for k in range(self.states)]
+            block_sum = sum(block_products, RealspaceSum.zero(self.modes))
+            product_matrix.set_block(i, j, block_sum)
+
+        return product_matrix
+
+    def __eq__(self, other: RealspaceMatrix):
+        if self.states != other.states:
+            return False
+
+        if self.modes != other.modes:
+            return False
+
+        if self._blocks != other._blocks:
+            return False
+
+        return True
+
+    def _partition_into_quadrants(self) -> Tuple[RealspaceMatrix]:
+        """Partitions the :class:`~.pennylane.labs.trotter_error.RealspaceMatrix` into four :class:`~.pennylane.labs.trotter_error.RealspaceMatrix` objects on ``self.states // 2`` states. This method assumes ``self.states`` is a power of two."""
+        # pylint: disable=chained-comparison
+        half = self.states // 2
+
+        top_left = RealspaceMatrix(half, self.modes, {})
+        top_right = RealspaceMatrix(half, self.modes, {})
+        bottom_left = RealspaceMatrix(half, self.modes, {})
+        bottom_right = RealspaceMatrix(half, self.modes, {})
+
+        for index, word in self._blocks.items():
+            x, y = index
+
+            if x < half and y < half:
+                top_left.set_block(x, y, word)
+            if x < half and y >= half:
+                top_right.set_block(x, y - half, word)
+            if x >= half and y < half:
+                bottom_left.set_block(x - half, y, word)
+            if x >= half and y >= half:
+                bottom_right.set_block(x - half, y - half, word)
+
+        return top_left, top_right, bottom_left, bottom_right
+
+    def apply(self, state: VibronicHO) -> VibronicHO:
+        """Apply the :class:`~.pennylane.labs.trotter_error.RealspaceMatrix` to an input :class:`~.pennylane.labs.trotter_error.VibronicHO` on the right.
+
+        Args:
+            state (VibronicHO): a vibronic wavefunction
+
+        Returns:
+            VibronicHO: the result of applying the :class:`~.pennylane.labs.trotter_error.RealspaceMatrix` to ``state``
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+        >>> from pennylane.labs.trotter_error import HOState, VibronicHO
+        >>> import numpy as np
+        >>> n_states = 1
+        >>> n_modes = 3
+        >>> gridpoints = 2
+        >>> op1 = RealspaceOperator(n_modes, (), RealspaceCoeffs.coeffs(np.array(1), label="lambda"))
+        >>> op2 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs.coeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+        >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+        >>> vib_matrix = RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum})
+        >>> state_dict = {(1, 0, 0): 1, (0, 1, 1): 1}
+        >>> state = HOState.from_dict(n_modes, gridpoints, state_dict)
+        >>> VibronicHO(n_states, n_modes, gridpoints, [state])
+        VibronicHO([HOState(modes=3, gridpoints=2, <Compressed Sparse Row sparse array of dtype 'int64'
+            with 2 stored elements and shape (8, 1)>
+          Coords	Values
+          (3, 0)	1
+          (4, 0)	1)])
+        """
+
+        if self.states != len(state.ho_states):
+            raise ValueError(
+                f"Cannot apply RealspaceMatrix on {self.states} states to VibronicHO on {len(state.ho_states)} states."
+            )
+
+        if self.modes != state.modes:
+            raise ValueError(
+                f"Cannot apply RealspaceMatrix on {self.modes} modes to VibronicHO on {state.modes} modes."
+            )
+
+        ho_states = []
+        for i in range(self.states):
+            ho = sum(
+                (self.block(i, j).apply(ho_state) for j, ho_state in enumerate(state.ho_states)),
+                HOState.zero_state(state.modes, state.gridpoints),
+            )
+            ho_states.append(ho)
+
+        return VibronicHO(
+            states=state.states,
+            modes=state.modes,
+            gridpoints=state.gridpoints,
+            ho_states=ho_states,
+        )
+
+    def get_coefficients(self, threshold: float = 0.0) -> Dict[Tuple[int, int], Dict]:
+        """Return a dictionary containing the coefficients of the :class:`~.pennylane.labs.trotter_error.RealspaceSum`
+
+        Args:
+            threshold (float): tolerance to return coefficients whose magnitude is greater than ``threshold``
+
+        Returns:
+            Dict: a dictionary whose keys are the indices of the :class:`~.pennylane.labs.trotter_error.RealspaceMatrix` and whose values are dictionaries obtained by :func:`~pennylane.labs.trotter_error.RealspaceSum.get_coefficients`
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceSum, RealspaceCoeffs, RealspaceMatrix
+        >>> import numpy as np
+        >>> n_states = 1
+        >>> n_modes = 5
+        >>> op1 = RealspaceOperator(n_modes, ("Q"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="phi"))
+        >>> op2 = RealspaceOperator(n_modes, ("P"), RealspaceCoeffs(np.array([1, 2, 3, 4, 5]), label="chi"))
+        >>> rs_sum = RealspaceSum(n_modes, [op1, op2])
+        >>> RealspaceMatrix(n_states, n_modes, {(0, 0): rs_sum}).get_coefficients()
+        {(0, 0): {'Q': {(0,): 1.0, (1,): 2.0, (2,): 3.0, (3,): 4.0, (4,): 5.0},
+        'P': {(0,): 1.0, (1,): 2.0, (2,): 3.0, (3,): 4.0, (4,): 5.0}}}
+        """
+        d = {}
+        for i, j in product(range(self.states), repeat=2):
+            d[(i, j)] = self.block(i, j).get_coefficients(threshold)
+
+        return d
+
+    @classmethod
+    def zero(cls, states: int, modes: int) -> RealspaceMatrix:
+        """Return a RealspaceMatrix representation of the zero operator"""
+        return cls(states, modes, {})
+
+    def __repr__(self):
+        return f"RealspaceMatrix({self._blocks})"
+
+
+def _is_pow_2(k: int) -> bool:
+    """Test if k is a power of two"""
+    return k & (k - 1) == 0
+
+
+def _next_pow_2(k: int) -> int:
+    """Return the smallest power of 2 greater than or equal to k"""
+    return 2 ** (k - 1).bit_length()
diff --git a/pennylane/labs/trotter_error/realspace/realspace_operator.py b/pennylane/labs/trotter_error/realspace/realspace_operator.py
new file mode 100644
index 000000000..4ea16e67c
--- /dev/null
+++ b/pennylane/labs/trotter_error/realspace/realspace_operator.py
@@ -0,0 +1,492 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""The RealspaceOperator class"""
+
+from __future__ import annotations
+
+import math
+from collections import defaultdict
+from itertools import product
+from typing import Dict, Sequence, Tuple, Union
+
+import numpy as np
+import scipy as sp
+
+from pennylane.labs.trotter_error import Fragment
+from pennylane.labs.trotter_error.realspace import HOState
+from pennylane.labs.trotter_error.realspace.matrix import (
+    _op_norm,
+    _string_to_matrix,
+    _tensor_with_identity,
+    _zeros,
+)
+
+from .realspace_coefficients import RealspaceCoeffs, _RealspaceTree
+
+
+class RealspaceOperator:
+    r"""Represents a linear combination of a product of position and momentum operators.
+    The ``RealspaceOperator`` class can be used to represent components of a vibrational
+    Hamiltonian, e.g., the following sum over a product of two position operators :math:`Q`:
+
+    .. math:: \sum_{i,j=1}^n \phi_{i,j}Q_i Q_j,
+
+    where :math:`\phi_{i, j}` represents the coefficient and is a constant.
+
+    Args:
+        modes (int): the number of vibrational modes
+        ops (Sequence[str]): a sequence representation of the position and momentum operators
+        coeffs (``RealspaceCoeffs``): an expression tree which evaluates the entries of the coefficient tensor
+
+    **Example**
+
+    This example uses :class:`~.pennylane.labs.trotter_error.RealspaceOperator` to build the
+    operator :math:`\sum_{i,j=1}^2 \phi_{i,j}Q_i Q_j`. The operator represents a sum over 2 modes
+    for the position operators :math:`Q_iQ_j`.
+
+    >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs
+    >>> import numpy as np
+    >>> n_modes = 2
+    >>> ops = ("Q", "Q")
+    >>> coeffs = RealspaceCoeffs(np.array([[1, 0], [0, 1]]), label="phi")
+    >>> RealspaceOperator(n_modes, ops, coeffs)
+    RealspaceOperator(5, ('Q', 'Q'), phi[idx0,idx1])
+    """
+
+    def __init__(
+        self, modes: int, ops: Sequence[str], coeffs: Union[RealspaceCoeffs, np.ndarray, float]
+    ) -> RealspaceOperator:
+
+        if coeffs.shape != (modes,) * len(ops):
+            raise ValueError(
+                f"coeffs has shape {coeffs.shape}, but shape {(modes, ) * len(ops)} was expected."
+            )
+
+        self.modes = modes
+        self.ops = ops
+        self.coeffs = coeffs
+
+    def matrix(
+        self, gridpoints: int, basis: str = "realspace", sparse: bool = False
+    ) -> Union[np.ndarray, sp.sparse.csr_array]:
+        """Return a matrix representation of the operator.
+
+        Args:
+            gridpoints (int): the number of gridpoints used to discretize the position or momentum operators
+            basis (str): the basis of the matrix, available options are ``realspace`` and ``harmonic``
+            sparse (bool): if ``True`` returns a sparse matrix, otherwise returns a dense matrix
+
+        Returns:
+            Union[ndarray, scipy.sparse.csr_array]: the matrix representation of the :class:`~.pennylane.labs.trotter_error.RealspaceOperator`
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs
+        >>> import numpy as np
+        >>> n_modes = 2
+        >>> ops = ("Q", "Q")
+        >>> coeffs = RealspaceCoeffs(np.array([[1, 0], [0, 1]]), label="phi")
+        >>> RealspaceOperator(n_modes, ops, coeffs).matrix(2)
+        [[6.28318531 0.         0.         0.        ]
+         [0.         3.14159265 0.         0.        ]
+         [0.         0.         3.14159265 0.        ]
+         [0.         0.         0.         0.        ]]
+        """
+
+        matrices = [
+            _string_to_matrix(op, gridpoints, basis=basis, sparse=sparse) for op in self.ops
+        ]
+        final_matrix = _zeros(shape=(gridpoints**self.modes, gridpoints**self.modes), sparse=sparse)
+
+        if sparse:
+            indices = list(self.coeffs.nonzero().keys())
+        else:
+            indices = product(range(self.modes), repeat=len(self.ops))
+
+        for index in indices:
+            matrix = self.coeffs[index] * _tensor_with_identity(
+                self.modes, gridpoints, index, matrices, sparse=sparse
+            )
+            final_matrix = final_matrix + matrix
+
+        return final_matrix
+
+    def __add__(self, other: RealspaceOperator) -> RealspaceOperator:
+        if self._is_zero:
+            return other
+
+        if other._is_zero:
+            return self
+
+        if self.ops != other.ops:
+            raise ValueError(f"Cannot add term {self.ops} with term {other.ops}.")
+
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot add RealspaceOperator on {self.modes} modes with RealspaceOperator on {other.modes} modes."
+            )
+
+        return RealspaceOperator(self.modes, self.ops, self.coeffs + other.coeffs)
+
+    def __sub__(self, other: RealspaceOperator) -> RealspaceOperator:
+        if self._is_zero:
+            return (-1) * other
+
+        if other._is_zero:
+            return self
+
+        if self.ops != other.ops:
+            raise ValueError(f"Cannot subtract term {self.ops} with term {other.ops}.")
+
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot subtract RealspaceOperator on {self.modes} modes with RealspaceOperator on {other.modes} modes."
+            )
+
+        return RealspaceOperator(self.modes, self.ops, self.coeffs - other.coeffs)
+
+    def __mul__(self, scalar: float) -> RealspaceOperator:
+        if np.isclose(scalar, 0):
+            return RealspaceOperator.zero(self.modes)
+
+        return RealspaceOperator(self.modes, self.ops, scalar * self.coeffs)
+
+    __rmul__ = __mul__
+
+    def __imul__(self, scalar: float) -> RealspaceOperator:
+        if np.isclose(scalar, 0):
+            return RealspaceOperator.zero(self.modes)
+
+        self.coeffs = _RealspaceTree.scalar_node(scalar, self.coeffs)
+        return self
+
+    def __matmul__(self, other: RealspaceOperator) -> RealspaceOperator:
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot multiply RealspaceOperator on {self.modes} modes with RealspaceOperator on {other.modes} modes."
+            )
+
+        return RealspaceOperator(self.modes, self.ops + other.ops, self.coeffs @ other.coeffs)
+
+    def __repr__(self) -> str:
+        return f"RealspaceOperator({self.modes}, {self.ops.__repr__()}, {self.coeffs.__repr__()})"
+
+    def __eq__(self, other: RealspaceOperator) -> bool:
+        if self.ops != other.ops:
+            return False
+
+        return self.coeffs == other.coeffs
+
+    @property
+    def _is_zero(self) -> bool:
+        """Always returns true when the operator is zero, but with false positives
+
+        Returns:
+            bool: if False, the operator is guarnateed to be non-zero, if True the operator is likely non-zero, but with some edge cases
+        """
+        return self.coeffs.is_zero
+
+    @classmethod
+    def zero(cls, modes) -> RealspaceOperator:
+        """Returns a ``RealspaceOperator`` representing the zero operator.
+
+        Args:
+            modes (int): the number of vibrational modes
+
+        Returns:
+            RealspaceOperator: a representation of the zero operator
+
+        """
+        return RealspaceOperator(modes, tuple(), RealspaceCoeffs(np.array(0)))
+
+    def get_coefficients(self, threshold: float = 0.0) -> Dict[Tuple[int], float]:
+        """Return the non-zero coefficients in a dictionary.
+
+        Args:
+            threshold (float): tolerance to return coefficients whose magnitude is greater than ``threshold``
+
+        Returns:
+            Dict[Tuple[int], float]: a dictionary whose keys are the nonzero indices, and values are the coefficients
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs
+        >>> import numpy as np
+        >>> n_modes = 2
+        >>> ops = ("Q", "Q")
+        >>> coeffs = RealspaceCoeffs(np.array([[1, 0], [0, 1]]), label="phi")
+        >>> RealspaceOperator(n_modes, ops, coeffs).get_coefficients()
+        {(0, 0): 1, (1, 1): 1}
+        """
+        return self.coeffs.nonzero(threshold)
+
+
+class RealspaceSum(Fragment):
+    r"""Represents a linear combination of :class:`~.pennylane.labs.trotter_error.RealspaceOperator` objects.
+
+    The :class:`~pennylane.labs.trotter_error.RealspaceSum` class can be used to represent a
+    Hamiltonian that is built from a sum of
+    :class:`~.pennylane.labs.trotter_error.RealspaceOperator` objects. For example, the vibrational
+    hamiltonian, adapted from Eq. 4 of `arXiv:1703.09313 <https://arxiv.org/abs/1703.09313>`_,
+
+    .. math:: \sum_i \frac{\omega_i}{2} P_i^2 + \sum_i \frac{\omega_i}{2} Q_i^2 + \sum_i \phi^{(1)}_i Q_i + \sum_{i,j} \phi^{(2)}_{ij} Q_i Q_j + \dots,
+
+    is a sum of terms where each term can be expressed by a :class:`~.pennylane.labs.trotter_error.RealspaceOperator`.
+
+    Args:
+        modes (int): the number of vibrational modes
+        ops (Sequence[RealspaceOperator]): a sequence containing :class:`~.pennylane.labs.trotter_error.RealspaceOperator` objects
+
+    **Example**
+
+    We can build the harmonic part of a vibrational Hamiltonian,
+    :math:`\sum_i \frac{\omega_i}{2} P_i^2 + \sum_i \frac{\omega_i}{2} Q_i^2`, with the following
+    code.
+
+    >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs, RealspaceSum
+    >>> import numpy as np
+    >>> n_modes = 2
+    >>> freqs = np.array([1.23, 3.45]) / 2
+    >>> coeffs = RealspaceCoeffs(freqs, label="omega")
+    >>> rs_op1 = RealspaceOperator(n_modes, ("PP",), coeffs)
+    >>> rs_op2 = RealspaceOperator(n_modes, ("QQ",), coeffs)
+    >>> RealspaceSum(n_modes, [rs_op1, rs_op2])
+    RealspaceSum((RealspaceOperator(2, ('PP',), omega[idx0]), RealspaceOperator(2, ('QQ',), omega[idx0])))
+    """
+
+    def __init__(self, modes: int, ops: Sequence[RealspaceOperator]):
+        # pylint: disable=unnecessary-lambda, protected-access
+        for op in ops:
+            if op.modes != modes:
+                raise ValueError(
+                    f"RealspaceSum on {modes} modes can only contain RealspaceOperators on {modes}. Found a RealspaceOperator on {op.modes} modes."
+                )
+
+        for op in ops:
+            assert op.coeffs is not None
+
+        ops = tuple(filter(lambda op: not op._is_zero, ops))
+        self._is_zero = len(ops) == 0
+
+        self.modes = modes
+        self._lookup = defaultdict(lambda: RealspaceOperator.zero(self.modes))
+
+        for op in ops:
+            self._lookup[op.ops] += op
+
+        for op in ops:
+            assert self._lookup[op.ops].coeffs is not None
+
+        self.ops = tuple(self._lookup.values())
+
+    def __add__(self, other: RealspaceSum) -> RealspaceSum:
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot add RealspaceSum on {self.modes} modes with RealspaceSum on {other.modes} modes."
+            )
+
+        l_ops = {term.ops for term in self.ops}
+        r_ops = {term.ops for term in other.ops}
+
+        new_ops = []
+
+        for op in l_ops.union(r_ops):
+            new_ops.append(self._lookup[op] + other._lookup[op])
+
+        return RealspaceSum(self.modes, new_ops)
+
+    def __sub__(self, other: RealspaceSum) -> RealspaceSum:
+        if self.modes != other.modes:
+            raise ValueError(
+                f"Cannot subtract RealspaceSum on {self.modes} modes with RealspaceSum on {other.modes} modes."
+            )
+
+        l_ops = {term.ops for term in self.ops}
+        r_ops = {term.ops for term in other.ops}
+
+        new_terms = []
+
+        for op in l_ops.union(r_ops):
+            new_terms.append(self._lookup[op] - other._lookup[op])
+
+        return RealspaceSum(self.modes, new_terms)
+
+    def __mul__(self, scalar: float) -> RealspaceSum:
+        return RealspaceSum(self.modes, [scalar * term for term in self.ops])
+
+    __rmul__ = __mul__
+
+    def __imul__(self, scalar: float) -> RealspaceSum:
+        for term in self.ops:
+            term *= scalar
+
+        return self
+
+    def __matmul__(self, other: RealspaceSum) -> RealspaceSum:
+        return RealspaceSum(
+            self.modes,
+            [
+                RealspaceOperator(
+                    self.modes, l_term.ops + r_term.ops, l_term.coeffs @ r_term.coeffs
+                )
+                for l_term, r_term in product(self.ops, other.ops)
+            ],
+        )
+
+    def __repr__(self) -> str:
+        return f"RealspaceSum({self.ops.__repr__()})"
+
+    def __eq__(self, other: RealspaceSum) -> bool:
+        return self._lookup == other._lookup
+
+    @classmethod
+    def zero(cls, modes: int) -> RealspaceSum:
+        """Returns a :class:`~.pennylane.labs.trotter_error.RealspaceOperator` representing the zero operator
+
+        Args:
+            modes (int): the number of vibrational modes (needed for consistency with arithmetic operations)
+
+        Returns:
+            RealspaceOperator: a representation of the zero operator
+
+        """
+        return RealspaceSum(modes, [RealspaceOperator.zero(modes)])
+
+    def matrix(
+        self, gridpoints: int, basis: str = "realspace", sparse: bool = False
+    ) -> Union[np.ndarray, sp.sparse.cs_array]:
+        """Return a matrix representation of the :class:`~pennylane.labs.trotter_error.RealspaceSum`.
+
+        Args:
+            gridpoints (int): the number of gridpoints used to discretize the position/momentum operators
+            basis (str): the basis of the matrix, available options are ``realspace`` and ``harmonic``
+            sparse (bool): if ``True`` returns a sparse matrix, otherwise a dense matrix
+
+        Returns:
+            Union[ndarray, scipy.sparse.csr_array]: the matrix representation of the :class:`~.pennylane.labs.trotter_error.RealspaceOperator`
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs, RealspaceSum
+        >>> import numpy as np
+        >>> n_modes = 2
+        >>> freqs = np.array([1.23, 3.45])
+        >>> coeffs = RealspaceCoeffs(freqs, label="omega")
+        >>> rs_op1 = RealspaceOperator(n_modes, ("PP",), coeffs)
+        >>> rs_op2 = RealspaceOperator(n_modes, ("QQ",), coeffs)
+        >>> RealspaceSum(n_modes, [rs_op1, rs_op2]).matrix(2)
+        [[22.05398043+0.00000000e+00j -5.41924733+6.63666389e-16j
+          -1.93207948+2.36611495e-16j  0.        +0.00000000e+00j]
+         [-5.41924733-6.63666389e-16j 11.21548577+0.00000000e+00j
+           0.        +0.00000000e+00j -1.93207948+2.36611495e-16j]
+         [-1.93207948-2.36611495e-16j  0.        +0.00000000e+00j
+          18.18982146+0.00000000e+00j -5.41924733+6.63666389e-16j]
+         [ 0.        +0.00000000e+00j -1.93207948-2.36611495e-16j
+          -5.41924733-6.63666389e-16j  7.35132681+0.00000000e+00j]]
+        """
+
+        final_matrix = _zeros(shape=(gridpoints**self.modes, gridpoints**self.modes), sparse=sparse)
+        for op in self.ops:
+            final_matrix = final_matrix + op.matrix(gridpoints, basis=basis, sparse=sparse)
+
+        return final_matrix
+
+    def norm(self, params: Dict) -> float:
+        """Returns an upper bound on the spectral norm of the operator.
+
+        Args:
+            params (Dict): The dictionary of parameters. The supported parameters are
+
+                * ``gridpoints`` (int): the number of gridpoints used to discretize the operator
+                * ``sparse`` (bool): If ``True``, use optimizations for sparse operators. Defaults to ``False``.
+
+        Returns:
+            float: an upper bound on the spectral norm of the operator
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs, RealspaceSum
+        >>> import numpy as np
+        >>> n_modes = 2
+        >>> freqs = np.array([1.23, 3.45])
+        >>> coeffs = RealspaceCoeffs(freqs, label="omega")
+        >>> rs_op1 = RealspaceOperator(n_modes, ("PP",), coeffs)
+        >>> rs_op2 = RealspaceOperator(n_modes, ("QQ",), coeffs)
+        >>> params = {"gridpoints": 2, "sparse": True}
+        >>> RealspaceSum(n_modes, [rs_op1, rs_op2]).norm(params)
+        29.405307237600457
+        """
+
+        sparse = params.get("sparse", False)
+        gridpoints = params.get("gridpoints", None)
+        if gridpoints is None:
+            raise KeyError("Need to specify the number of gridpoints")
+
+        norm = 0
+
+        for op in self.ops:
+            term_op_norm = math.prod(map(lambda op: _op_norm(gridpoints) ** len(op), op.ops))
+
+            if sparse:
+                coeffs = op.coeffs.nonzero()
+                coeff_sum = sum(abs(val) for val in coeffs.values())
+            else:
+                indices = product(range(self.modes), repeat=len(op.ops))
+                coeff_sum = sum(abs(op.coeffs[index]) for index in indices)
+
+            norm += coeff_sum * term_op_norm
+
+        return norm
+
+    def apply(self, state: HOState) -> HOState:
+        """Apply the :class:`~.pennylane.labs.trotter_error.RealspaceSum` to an input :class:`~.pennylane.labs.trotter_error.HOState` object."""
+        if not isinstance(state, HOState):
+            raise TypeError
+
+        mat = self.matrix(state.gridpoints, basis="harmonic", sparse=True)
+
+        return HOState(
+            state.modes,
+            state.gridpoints,
+            mat @ state.vector,
+        )
+
+    def get_coefficients(self, threshold: float = 0.0) -> Dict[Tuple[str], Dict]:
+        """Return a dictionary containing the non-zero coefficients of the :class:`~pennylane.labs.trotter_error.RealspaceSum`.
+
+        Args:
+            threshold (float): tolerance to return coefficients whose magnitude is greater than ``threshold``
+
+        Returns:
+            Dict: a dictionary whose keys correspond to the RealspaceOperators in the sum, and whose
+                values are dictionaries obtained by :func:`~.pennylane.labs.trotter_error.RealspaceOperator.get_coefficients`
+
+        **Example**
+
+        >>> from pennylane.labs.trotter_error import RealspaceOperator, RealspaceCoeffs, RealspaceSum
+        >>> import numpy as np
+        >>> n_modes = 2
+        >>> freqs = np.array([1.23, 3.45])
+        >>> coeffs = RealspaceCoeffs(freqs, label="omega")
+        >>> rs_op1 = RealspaceOperator(n_modes, ("PP",), coeffs)
+        >>> rs_op2 = RealspaceOperator(n_modes, ("QQ",), coeffs)
+        >>> RealspaceSum(n_modes, [rs_op1, rs_op2]).get_coefficients()
+        {('PP',): {(0,): 1.23, (1,): 3.45}, ('QQ',): {(0,): 1.23, (1,): 3.45}}
+        """
+
+        coeffs = {}
+        for op in self.ops:
+            coeffs[op.ops] = op.get_coefficients(threshold)
+
+        return coeffs


Test patch:
diff --git a/pennylane/labs/tests/trotter_error/fragments/test_trotter.py b/pennylane/labs/tests/trotter_error/fragments/test_trotter.py
new file mode 100644
index 000000000..a5020dd84
--- /dev/null
+++ b/pennylane/labs/tests/trotter_error/fragments/test_trotter.py
@@ -0,0 +1,71 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the second order trotter error calculation"""
+
+import numpy as np
+import pytest
+
+from pennylane.labs.trotter_error.abstract import nested_commutator
+from pennylane.labs.trotter_error.fragments import vibronic_fragments
+from pennylane.labs.trotter_error.product_formulas import trotter_error
+from pennylane.labs.trotter_error.realspace import RealspaceMatrix
+
+
+def _coeffs(states: int, modes: int, order: int):
+    """Produce random coefficients used to construct a Vibronic Hamiltonian. The coefficients are modified to reflect the symmetries obtained by commuting operators."""
+
+    phis = []
+    symmetric_phis = []
+    for i in range(order + 1):
+        shape = (states, states) + (modes,) * i
+        phi = np.random.random(size=shape)
+        phis.append(phi)
+        symmetric_phis.append(np.zeros(shape))
+
+    for phi in phis:
+        for i in range(states):
+            for j in range(states):
+                phi[i, j] = (phi[i, j] + phi[i, j].T) / 2
+
+    for phi, symmetric_phi in zip(phis, symmetric_phis):
+        for i in range(states):
+            for j in range(states):
+                symmetric_phi[i, j] = (phi[i, j] + phi[j, i]) / 2
+
+    return np.random.random(size=(modes,)), symmetric_phis
+
+
+@pytest.mark.parametrize("modes", range(5))
+def test_second_order_trotter_error_operator(modes):
+    """Test that the second order Trotter error operator is correct for a 2 state example"""
+    states = 2
+    delta = 0.72
+    scalar = -(delta**2) / 24
+    fragments = vibronic_fragments(states, modes, *_coeffs(states, modes, order=2))
+
+    terms = [
+        nested_commutator([fragments[0], fragments[0], fragments[1]]),
+        2 * nested_commutator([fragments[1], fragments[0], fragments[1]]),
+        2 * nested_commutator([fragments[2], fragments[0], fragments[1]]),
+        nested_commutator([fragments[0], fragments[0], fragments[2]]),
+        2 * nested_commutator([fragments[1], fragments[0], fragments[2]]),
+        2 * nested_commutator([fragments[2], fragments[0], fragments[2]]),
+        nested_commutator([fragments[1], fragments[1], fragments[2]]),
+        2 * nested_commutator([fragments[2], fragments[1], fragments[2]]),
+    ]
+
+    actual = trotter_error(fragments, delta)
+    expected = scalar * sum(terms, RealspaceMatrix(states, modes))
+
+    assert actual == expected
diff --git a/pennylane/labs/tests/trotter_error/fragments/test_vibrational.py b/pennylane/labs/tests/trotter_error/fragments/test_vibrational.py
new file mode 100644
index 000000000..0a61784d8
--- /dev/null
+++ b/pennylane/labs/tests/trotter_error/fragments/test_vibrational.py
@@ -0,0 +1,318 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Test the VibrationalHamiltonian class"""
+
+import itertools
+
+import numpy as np
+import pytest
+import scipy as sp
+from scipy.sparse import csr_array
+
+from pennylane.labs.trotter_error.fragments import vibrational_fragments
+from pennylane.labs.trotter_error.realspace import HOState, RealspaceSum
+
+# pylint: disable=no-self-use
+
+
+def _vibrational_hamiltonian(modes, freqs, taylor_coeffs):
+    frags = vibrational_fragments(modes, freqs, taylor_coeffs)
+    return sum(frags, RealspaceSum.zero(modes))
+
+
+@pytest.mark.parametrize(
+    "n_modes, omegas, phis, gridpoints",
+    [
+        (5, np.random.random(5), [], 2),
+        (5, np.random.random(5), [np.random.random(size=(5,) * i) for i in range(4)], 2),
+    ],
+)
+def test_fragementation_schemes_equal(n_modes, omegas, phis, gridpoints):
+    """Test that harmonic + anharmonic = kinetic + potential"""
+
+    harmonic, anharmonic = vibrational_fragments(n_modes, omegas, phis, frag_method="harmonic")
+    kinetic, potential = vibrational_fragments(n_modes, omegas, phis, frag_method="kinetic")
+
+    mat1 = (harmonic + anharmonic).matrix(gridpoints)
+    mat2 = (kinetic + potential).matrix(gridpoints)
+
+    assert np.allclose(mat1, mat2)
+
+
+class TestHarmonic1Mode:
+    """Test the vibrational Hamiltonian on a single mode"""
+
+    freq = 1.2345
+    n_states = 5
+    omegas = np.array([freq])
+    ham = _vibrational_hamiltonian(1, omegas, [])
+    states = [HOState(1, 10, {(i,): 1}) for i in range(n_states)]
+
+    @pytest.mark.parametrize("n_states, freq, ham, states", [(n_states, freq, ham, states)])
+    def test_expectation_1_mode(self, n_states, freq, ham, states):
+        """Test the expectation computation against known values"""
+
+        expected = np.diag(np.arange(n_states) + 0.5) * freq
+        actual = np.zeros(shape=(n_states, n_states), dtype=np.complex128)
+
+        for i, state1 in enumerate(states):
+            for j, state2 in enumerate(states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+    @pytest.mark.parametrize("n_states, freq, ham, states", [(n_states, freq, ham, states)])
+    def test_linear_combination_harmonic(self, n_states, freq, ham, states):
+        """Test the expectation of a linear combination of harmonics"""
+
+        rot_log = np.zeros(shape=(n_states, n_states))
+
+        for i in range(n_states):
+            for j in range(i):
+                rot_log[i, j] = (i + 1) * (j + 1)
+                rot_log[j, i] = -rot_log[i, j]
+
+        rot = sp.linalg.expm(rot_log)
+
+        comb_states = []
+        for i in range(n_states):
+            state = sum((states[j] * rot[j, i] for j in range(n_states)), HOState.zero_state(1, 10))
+            comb_states.append(state)
+
+        expected = rot.T @ (np.diag(np.arange(n_states) + 0.5) * freq) @ rot
+        actual = np.zeros(shape=(n_states, n_states), dtype=np.complex128)
+
+        for i, state1 in enumerate(comb_states):
+            for j, state2 in enumerate(comb_states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+
+class TestHarmonicMultiMode:
+    """Test the vibrational Hamiltonian with multiple modes"""
+
+    n_states = 3
+    omegas = np.array([1, 2.3])
+    ham = _vibrational_hamiltonian(2, omegas, [])
+    states = [HOState(2, 10, {(i, j): 1}) for i, j in itertools.product(range(n_states), repeat=2)]
+    excitations = list(itertools.product(range(n_states), repeat=2))
+
+    @pytest.mark.parametrize(
+        "omegas, ham, states, excitations", [(omegas, ham, states, excitations)]
+    )
+    def test_harmonic(self, omegas, ham, states, excitations):
+        """Test the expectation value of a harmonic"""
+
+        expected = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i in range(len(states)):
+            expected[i, i] = (0.5 + excitations[i][0]) * omegas[0] + (
+                0.5 + excitations[i][1]
+            ) * omegas[1]
+
+        actual = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i, state1 in enumerate(states):
+            for j, state2 in enumerate(states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+    @pytest.mark.parametrize(
+        "omegas, ham, states, excitations", [(omegas, ham, states, excitations)]
+    )
+    def test_linear_combination_harmonic(self, omegas, ham, states, excitations):
+        """Test the expectation value of a linear combintaion of harmonics"""
+
+        rot_log = np.zeros((len(states), len(states)))
+        for i in range(len(states)):
+            for j in range(i):
+                rot_log[i, j] = (i + 1) * (j + 1)
+                rot_log[j, i] = -rot_log[i, j]
+
+        rot = sp.linalg.expm(rot_log)
+
+        comb_states = []
+        for i in range(len(states)):
+            state = sum(
+                (states[j] * rot[j, i] for j in range(len(states))), HOState.zero_state(2, 10)
+            )
+            comb_states.append(state)
+
+        expected = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i in range(len(states)):
+            expected[i, i] = (0.5 + excitations[i][0]) * omegas[0] + (
+                0.5 + excitations[i][1]
+            ) * omegas[1]
+
+        expected = rot.T @ expected @ rot
+
+        actual = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i, state1 in enumerate(comb_states):
+            for j, state2 in enumerate(comb_states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+
+class TestExpectation:
+    """Test the expectation values against precomputed expected results"""
+
+    phis1 = [np.array(0), np.array([0.0]), np.array([[-0.00306787]]), np.array([[[0.00051144]]])]
+    omegas1 = np.array([0.0249722])
+    mat1 = np.array(
+        [[0.01095216 + 0.0j, 0.00018082 + 0.0j], [0.00018082 + 0.0j, 0.01095216 + 0.0j]]
+    )
+
+    phis2 = [
+        np.array(0),
+        np.array([0.0, 0.0, 0.0]),
+        np.array(
+            [
+                [-6.21741719e-14, 0.00000000e00, 0.00000000e00],
+                [1.99211595e-13, 9.46571153e-04, 0.00000000e00],
+                [-5.71982634e-14, 8.73074134e-04, -1.31910690e-03],
+            ]
+        ),
+        np.array(
+            [
+                [
+                    [5.53904193e-05, 0.00000000e00, 0.00000000e00],
+                    [0.00000000e00, 0.00000000e00, 0.00000000e00],
+                    [0.00000000e00, 0.00000000e00, 0.00000000e00],
+                ],
+                [
+                    [-4.98906472e-15, 0.00000000e00, 0.00000000e00],
+                    [-3.00847220e-03, 7.84674613e-05, 0.00000000e00],
+                    [0.00000000e00, 0.00000000e00, 0.00000000e00],
+                ],
+                [
+                    [5.97655606e-15, 0.00000000e00, 0.00000000e00],
+                    [3.88965240e-14, -9.17507861e-05, 0.00000000e00],
+                    [-2.99010685e-03, -3.09965385e-03, 1.61908332e-04],
+                ],
+            ]
+        ),
+    ]
+    omegas2 = np.array([0.00978462, 0.00978489, 0.01663723])
+    mat2 = np.array(
+        [
+            [
+                1.79171056e-02,
+                2.48044382e-05,
+                -1.06815069e-03,
+                4.36537067e-04,
+                -2.10123449e-03,
+                -2.85917137e-14,
+                9.96152552e-14,
+                1.37476535e-14,
+            ],
+            [
+                2.48044382e-05,
+                1.79171056e-02,
+                4.36537067e-04,
+                -1.06815069e-03,
+                -2.85917137e-14,
+                -2.10123449e-03,
+                1.37476535e-14,
+                9.96152552e-14,
+            ],
+            [
+                -1.06815069e-03,
+                4.36537067e-04,
+                1.79171056e-02,
+                2.48044382e-05,
+                9.96152552e-14,
+                1.37476535e-14,
+                -2.10123449e-03,
+                -2.85917137e-14,
+            ],
+            [
+                4.36537067e-04,
+                -1.06815069e-03,
+                2.48044382e-05,
+                1.79171056e-02,
+                1.37476535e-14,
+                9.96152552e-14,
+                -2.85917137e-14,
+                -2.10123449e-03,
+            ],
+            [
+                -2.10123449e-03,
+                -2.85917137e-14,
+                9.96152552e-14,
+                1.37476535e-14,
+                1.79171056e-02,
+                2.48044382e-05,
+                -1.06815069e-03,
+                4.36537067e-04,
+            ],
+            [
+                -2.85917137e-14,
+                -2.10123449e-03,
+                1.37476535e-14,
+                9.96152552e-14,
+                2.48044382e-05,
+                1.79171056e-02,
+                4.36537067e-04,
+                -1.06815069e-03,
+            ],
+            [
+                9.96152552e-14,
+                1.37476535e-14,
+                -2.10123449e-03,
+                -2.85917137e-14,
+                -1.06815069e-03,
+                4.36537067e-04,
+                1.79171056e-02,
+                2.48044382e-05,
+            ],
+            [
+                1.37476535e-14,
+                9.96152552e-14,
+                -2.85917137e-14,
+                -2.10123449e-03,
+                4.36537067e-04,
+                -1.06815069e-03,
+                2.48044382e-05,
+                1.79171056e-02,
+            ],
+        ]
+    )
+
+    @pytest.mark.parametrize(
+        "phis, omegas, expected", [(phis1, omegas1, mat1), (phis2, omegas2, mat2)]
+    )
+    def test_matrix(self, phis, omegas, expected):
+        """Test the matrix"""
+
+        ham = _vibrational_hamiltonian(len(omegas), omegas, phis)
+        actual = ham.matrix(2, basis="harmonic")
+
+        assert np.allclose(actual, expected)
+
+    @pytest.mark.parametrize(
+        "phis, omegas, expected", [(phis1, omegas1, mat1), (phis2, omegas2, mat2)]
+    )
+    def test_expectation(self, phis, omegas, expected):
+        """Test the expectation values"""
+
+        n_modes = len(omegas)
+        eigvals, eigvecs = np.linalg.eig(expected)
+
+        ham = _vibrational_hamiltonian(n_modes, omegas, phis)
+
+        for i, eigval in enumerate(eigvals):
+            eigvec = eigvecs[:, i]
+            ho_state = HOState(n_modes, 2, csr_array(eigvec.reshape(2**n_modes, 1)))
+            assert np.isclose(ham.expectation(ho_state, ho_state), eigval)
diff --git a/pennylane/labs/tests/trotter_error/fragments/test_vibronic.py b/pennylane/labs/tests/trotter_error/fragments/test_vibronic.py
new file mode 100644
index 000000000..bb27d355c
--- /dev/null
+++ b/pennylane/labs/tests/trotter_error/fragments/test_vibronic.py
@@ -0,0 +1,151 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for Vibronic Hamiltonian"""
+
+from itertools import product
+
+import numpy as np
+import pytest
+import scipy as sp
+
+from pennylane.labs.trotter_error.fragments import vibronic_fragments
+from pennylane.labs.trotter_error.realspace import HOState, RealspaceMatrix, VibronicHO
+
+# pylint: disable=no-self-use
+
+
+def _vibronic_hamiltonian(states, modes, freqs, taylor_coeffs):
+    frags = vibronic_fragments(states, modes, freqs, taylor_coeffs)
+    return sum(frags, RealspaceMatrix.zero(states, modes))
+
+
+class Test1Mode:
+    """Test a simple one mode, one state vibronic Hamiltonian"""
+
+    freq = 1.2345
+    n_states = 5
+    omegas = np.array([freq])
+    ham = _vibronic_hamiltonian(1, 1, omegas, [])
+    states = [VibronicHO(1, 1, 10, [HOState(1, 10, {(i,): 1})]) for i in range(n_states)]
+
+    @pytest.mark.parametrize("n_states, freq, ham, states", [(n_states, freq, ham, states)])
+    def test_expectation_1_mode(self, n_states, freq, ham, states):
+        """Test the expectation computation against known values"""
+
+        expected = np.diag(np.arange(n_states) + 0.5) * freq
+        actual = np.zeros(shape=(n_states, n_states), dtype=np.complex128)
+
+        for i, state1 in enumerate(states):
+            for j, state2 in enumerate(states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+    @pytest.mark.parametrize("n_states, freq, ham, states", [(n_states, freq, ham, states)])
+    def test_linear_combination_1_mode(self, n_states, freq, ham, states):
+        """Test the expectation of a linear combination of harmonics"""
+
+        rot_log = np.zeros(shape=(n_states, n_states))
+
+        for i in range(n_states):
+            for j in range(i):
+                rot_log[i, j] = (i + 1) * (j + 1)
+                rot_log[j, i] = -(i + 1) * (j + 1)
+
+        rot = sp.linalg.expm(rot_log)
+
+        comb_states = []
+        for i in range(n_states):
+            state = sum(
+                (states[j] * rot[j, i] for j in range(n_states)), VibronicHO.zero_state(1, 1, 10)
+            )
+            comb_states.append(state)
+
+        expected = rot.T @ (np.diag(np.arange(n_states) + 0.5) * freq) @ rot
+        actual = np.zeros(shape=(n_states, n_states), dtype=np.complex128)
+
+        for i, state1 in enumerate(comb_states):
+            for j, state2 in enumerate(comb_states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+
+class TestHarmonic:
+    """Test a simple 1 state, 2 mode vibronic Hamiltonian"""
+
+    n_states = 3
+    omegas = np.array([1, 2.3])
+    ham = _vibronic_hamiltonian(1, 2, omegas, [])
+    states = [
+        VibronicHO(1, 2, 10, [HOState(2, 10, {(i, j): 1})])
+        for i, j in product(range(n_states), repeat=2)
+    ]
+
+    excitations = list(product(range(n_states), repeat=2))
+
+    @pytest.mark.parametrize(
+        "omegas, ham, states, excitations", [(omegas, ham, states, excitations)]
+    )
+    def test_harmonic(self, omegas, ham, states, excitations):
+        """Test the expectation value of a harmonic"""
+
+        expected = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i in range(len(states)):
+            expected[i, i] = (0.5 + excitations[i][0]) * omegas[0] + (
+                0.5 + excitations[i][1]
+            ) * omegas[1]
+
+        actual = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i, state1 in enumerate(states):
+            for j, state2 in enumerate(states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
+
+    @pytest.mark.parametrize(
+        "omegas, ham, states, excitations", [(omegas, ham, states, excitations)]
+    )
+    def test_linear_combination_harmonic(self, omegas, ham, states, excitations):
+        """Test the expectation value of a linear combintaion of harmonics"""
+
+        rot_log = np.zeros((len(states), len(states)))
+        for i in range(len(states)):
+            for j in range(i):
+                rot_log[i, j] = (i + 1) * (j + 1)
+                rot_log[j, i] = -(i + 1) * (j + 1)
+
+        rot = sp.linalg.expm(rot_log)
+
+        comb_states = []
+        for i in range(len(states)):
+            state = sum(
+                (states[j] * rot[j, i] for j in range(len(states))), VibronicHO.zero_state(1, 2, 10)
+            )
+            comb_states.append(state)
+
+        expected = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i in range(len(states)):
+            expected[i, i] = (0.5 + excitations[i][0]) * omegas[0] + (
+                0.5 + excitations[i][1]
+            ) * omegas[1]
+
+        expected = rot.T @ expected @ rot
+
+        actual = np.zeros((len(states), len(states)), dtype=np.complex128)
+        for i, state1 in enumerate(comb_states):
+            for j, state2 in enumerate(comb_states):
+                actual[i, j] = ham.expectation(state1, state2)
+
+        assert np.allclose(actual, expected)
diff --git a/pennylane/labs/tests/trotter_error/realspace/test_matrix.py b/pennylane/labs/tests/trotter_error/realspace/test_matrix.py
new file mode 100644
index 000000000..7bdcd9ea6
--- /dev/null
+++ b/pennylane/labs/tests/trotter_error/realspace/test_matrix.py
@@ -0,0 +1,151 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Test matrix operations"""
+
+from typing import Dict, Tuple
+
+import numpy as np
+import pytest
+
+from pennylane.labs.trotter_error.realspace import (
+    RealspaceCoeffs,
+    RealspaceMatrix,
+    RealspaceOperator,
+    RealspaceSum,
+)
+from pennylane.labs.trotter_error.realspace.matrix import _op_norm
+
+# pylint: disable=too-many-arguments,too-many-positional-arguments,no-self-use
+
+
+def _coeffs(states: int, modes: int, order: int):
+    """Produce random coefficients for input"""
+
+    phis = []
+    symmetric_phis = []
+    for i in range(order + 1):
+        shape = (states, states) + (modes,) * i
+        phi = np.random.random(size=shape)
+        phis.append(phi)
+        symmetric_phis.append(np.zeros(shape))
+
+    for phi in phis:
+        for i in range(states):
+            for j in range(states):
+                phi[i, j] = (phi[i, j] + phi[i, j].T) / 2
+
+    for phi, symmetric_phi in zip(phis, symmetric_phis):
+        for i in range(states):
+            for j in range(states):
+                symmetric_phi[i, j] = (phi[i, j] + phi[j, i]) / 2
+
+    return np.random.random(size=modes), symmetric_phis
+
+
+vword0 = RealspaceSum(2, [RealspaceOperator(2, tuple(), RealspaceCoeffs(np.array(0.5)))])
+blocks0 = {(0, 0): vword0}
+vmat0 = RealspaceMatrix(1, 2, blocks0)
+
+vword1 = RealspaceSum(
+    1, [RealspaceOperator(1, ("P",), RealspaceCoeffs(np.array([1]), label="omega"))]
+)
+blocks1 = {(0, 0): vword1}
+vmat1 = RealspaceMatrix(1, 1, blocks1)
+
+vword2a = RealspaceSum(
+    2,
+    [RealspaceOperator(2, ("P", "P"), RealspaceCoeffs(np.array([[0, 1], [2, 3]]), label="omega"))],
+)
+vword2b = RealspaceSum(
+    2, [RealspaceOperator(2, ("P",), RealspaceCoeffs(np.array([1, 2]), label="omega"))]
+)
+blocks2 = {(0, 0): vword2a, (1, 1): vword2b}
+vmat2 = RealspaceMatrix(2, 2, blocks2)
+
+blocks3 = {(0, 1): vword2a, (1, 0): vword2b}
+vmat3 = RealspaceMatrix(2, 2, blocks3)
+
+blocks4 = {(0, 0): vword2a, (0, 1): vword2a, (1, 0): vword2b, (1, 1): vword2b}
+vmat4 = RealspaceMatrix(2, 2, blocks4)
+
+
+class TestMatrix:
+    """Test properties of the RealspaceMatrix class"""
+
+    params = [
+        (blocks2, 2, 2, 2, False, 6 * _op_norm(2) ** 2),
+        (blocks2, 4, 2, 2, False, 6 * _op_norm(4) ** 2),
+        (blocks2, 2, 2, 2, True, 6 * _op_norm(2) ** 2),
+        (blocks2, 4, 2, 2, True, 6 * _op_norm(4) ** 2),
+        (blocks3, 2, 2, 2, False, np.sqrt(18 * _op_norm(2) ** 3)),
+        (blocks3, 4, 2, 2, False, np.sqrt(18 * _op_norm(4) ** 3)),
+        (blocks3, 2, 2, 2, True, np.sqrt(18 * _op_norm(2) ** 3)),
+        (blocks3, 4, 2, 2, True, np.sqrt(18 * _op_norm(4) ** 3)),
+        (blocks4, 2, 2, 2, False, 6 * _op_norm(2) ** 2 + np.sqrt(18 * _op_norm(2) ** 3)),
+        (blocks4, 4, 2, 2, False, 6 * _op_norm(4) ** 2 + np.sqrt(18 * _op_norm(4) ** 3)),
+        (blocks4, 2, 2, 2, True, 6 * _op_norm(2) ** 2 + np.sqrt(18 * _op_norm(2) ** 3)),
+        (blocks4, 4, 2, 2, True, 6 * _op_norm(4) ** 2 + np.sqrt(18 * _op_norm(4) ** 3)),
+    ]
+
+    @pytest.mark.parametrize("blocks, gridpoints, states, modes, sparse, expected", params)
+    def test_norm(
+        self,
+        blocks: Dict[Tuple[int], RealspaceSum],
+        gridpoints: int,
+        states: int,
+        modes: int,
+        sparse: bool,
+        expected: float,
+    ):
+        """Test that the norm is correct"""
+
+        vmatrix = RealspaceMatrix(states, modes, blocks)
+
+        params = {"gridpoints": gridpoints, "sparse": sparse}
+
+        assert np.isclose(vmatrix.norm(params), expected)
+
+    params = [
+        (blocks2, 2, 2, 2, False),
+        (blocks2, 4, 2, 2, False),
+        (blocks2, 2, 2, 2, True),
+        (blocks2, 4, 2, 2, True),
+        (blocks3, 2, 2, 2, False),
+        (blocks3, 4, 2, 2, False),
+        (blocks3, 2, 2, 2, True),
+        (blocks3, 4, 2, 2, True),
+        (blocks4, 2, 2, 2, False),
+        (blocks4, 4, 2, 2, False),
+        (blocks4, 2, 2, 2, True),
+        (blocks4, 4, 2, 2, True),
+    ]
+
+    @pytest.mark.parametrize("blocks, gridpoints, states, modes, sparse", params)
+    def test_norm_against_numpy(
+        self,
+        blocks: Dict[Tuple[int], RealspaceSum],
+        gridpoints: int,
+        states: int,
+        modes: int,
+        sparse: bool,
+    ):
+        """Test that .norm is an upper bound on the true norm"""
+
+        params = {"gridpoints": gridpoints, "sparse": sparse}
+
+        vmatrix = RealspaceMatrix(states, modes, blocks)
+        upper_bound = vmatrix.norm(params)
+        norm = np.abs(np.max(np.linalg.eigvals(vmatrix.matrix(gridpoints))))
+
+        assert np.isclose(norm, upper_bound) or norm < upper_bound