diff --git a/.pylintdict b/.pylintdict
index d1ef95662..08d3581c5 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -61,6 +61,7 @@ bopes
 boson
 bosons
 bosonic
+bosoniclinearmapper
 bosonicop
 bravais
 bravyi
diff --git a/qiskit_nature/second_q/mappers/__init__.py b/qiskit_nature/second_q/mappers/__init__.py
index 9d98b9515..6b33d3ee9 100644
--- a/qiskit_nature/second_q/mappers/__init__.py
+++ b/qiskit_nature/second_q/mappers/__init__.py
@@ -90,6 +90,16 @@
 
    TaperedQubitMapper
 
+   
+MixedOp Mappers
++++++++++++++++
+
+.. autosummary::
+   :toctree: ../stubs/
+   :nosignatures:
+
+   MixedMapper
+
 """
 
 from .bksf import BravyiKitaevSuperFastMapper
@@ -103,6 +113,7 @@
 from .qubit_mapper import QubitMapper
 from .interleaved_qubit_mapper import InterleavedQubitMapper
 from .tapered_qubit_mapper import TaperedQubitMapper
+from .mixed_mapper import MixedMapper
 
 __all__ = [
     "BravyiKitaevMapper",
@@ -113,6 +124,7 @@
     "LinearMapper",
     "BosonicLinearMapper",
     "LogarithmicMapper",
+    "MixedMapper",
     "QubitMapper",
     "InterleavedQubitMapper",
     "TaperedQubitMapper",
diff --git a/qiskit_nature/second_q/mappers/mixed_mapper.py b/qiskit_nature/second_q/mappers/mixed_mapper.py
new file mode 100644
index 000000000..442e6951b
--- /dev/null
+++ b/qiskit_nature/second_q/mappers/mixed_mapper.py
@@ -0,0 +1,175 @@
+# This code is part of a Qiskit project.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""The Mixed Mapper class."""
+
+from __future__ import annotations
+
+import logging
+from abc import ABC
+from functools import reduce
+
+from qiskit.quantum_info import SparsePauliOp
+
+from qiskit_algorithms.list_or_dict import ListOrDict as ListOrDictType
+
+from qiskit_nature.second_q.operators import MixedOp, SparseLabelOp
+
+from .qubit_mapper import QubitMapper, _ListOrDict
+
+LOGGER = logging.getLogger(__name__)
+
+
+class MixedMapper(ABC):
+    """Mapper of a Mixed Operator to a Qubit Operator.
+
+    This class is intended to be used for handling the mapping of composed fermionic and (or) bosonic
+    systems, defined as :class:`~qiskit_nature.second_q.operators.MixedOp`, into qubit operators.
+
+    Please note that the creation and usage of this class requires the precise definition of the
+    composite Hilbert size corresponding to the problem.
+    The ordering of the qubit registers associated to the bosonic and fermionic degrees of freedom
+    (for example) must be provided by the user through the definition of the hilbert space
+    registers dictionary. This ordering corresponds to a specific way to take the tensor product
+    of the fermionic and bosonic operators.
+
+    .. note::
+
+      This class is limited to one instance of a Fermionic Hilbert space, to ensure the anticommutation
+      relations of the fermions.
+
+    .. note::
+
+      This class enforces the register lengths to the mappers. Note that for the bosonic mappers, the
+      register lengths is not directly equal to the qubit register length but to the number of modes.
+      See the documentation of the class :class:``~.BosonicLinearMapper``.
+
+    The following attributes can be read and updated once the ``MixedMapper`` object has been
+    constructed.
+
+    Attributes:
+        mappers: Dictionary of mappers corresponding to "local" Hilbert spaces of the global problem.
+        hilbert_space_register_lengths: Ordered dictionary of local registers and their respective sizes.
+    """
+
+    def __init__(self, mappers: dict[str, QubitMapper], hilbert_space_register_lengths: dict):
+        """
+        Args:
+            mappers: Dictionary of mappers corresponding to the "local" Hilbert spaces.
+            hilbert_space_register_lengths: Ordered dictionary of local registers with their sizes.
+        """
+        super().__init__()
+        self.mappers = mappers
+        self.hilbert_space_register_lengths = hilbert_space_register_lengths
+
+    def _map_tuple_product(
+        self, active_indices: tuple[str], active_operators: tuple[SparseLabelOp]
+    ) -> SparsePauliOp:
+        """Maps a product of operators defined on the local Hilbert spaces defined by the active
+        indices. Note that the order of the active indices is not relevant. The only relevant ordering
+        is that given by :attr:`~MixedMapper.hilbert_space_register_lengths` at initialization.
+
+        When the operator is not present in the tuple, we use a padding operator with identities.
+
+        Args:
+            active_indices: Reference names of the Hilbert spaces on which the operator acts.
+            active_operators: List of operators to compose..
+        """
+
+        product_op_dict = {
+            index: SparsePauliOp("I" * value)
+            for index, value in self.hilbert_space_register_lengths.items()
+        }
+
+        for active_index, active_op in zip(active_indices, active_operators):
+            register_length = self.hilbert_space_register_lengths[active_index]
+            product_op_dict[active_index] = self.mappers[active_index].map(
+                active_op, register_length=register_length
+            )
+
+        product_op = reduce(SparsePauliOp.tensor, list(product_op_dict.values()))
+
+        return product_op
+
+    def _distribute_map(
+        self, operator_dict: dict[tuple[str], list[tuple[float, SparseLabelOp]]]
+    ) -> SparsePauliOp:
+        """Distributes the mapping of operators to each of the terms defined across specific
+        Hilbert spaces.
+
+        Args:
+            operator_dict: Dictionary of (key, operator list) pairs where the key specify which
+                local "Hilbert space" the operators act on. Note that the first element of the
+                operator list is the coefficient of this operator product across these Hilbert spaces.
+        """
+
+        mapped_op: SparsePauliOp = 0
+        for active_indices, operator_list in operator_dict.items():
+            for coef_and_operators in operator_list:
+                coef: float = coef_and_operators[0]
+                active_operators: tuple[SparseLabelOp] = coef_and_operators[1:]
+                mapped_op += coef * self._map_tuple_product(active_indices, active_operators)
+
+        return mapped_op.simplify()
+
+    def _map_single(
+        self,
+        mixed_op: MixedOp,
+        *,
+        register_length: int | None = None,
+    ) -> SparsePauliOp:
+        """Map the :class:`~qiskit_nature.second_q.operators.MixedOp` into a qubit operator.
+
+        The ``MixedOp`` is a representation of sums of products of operators corresponding to different
+        Hilbert spaces. The mapping procedure first runs through all of the terms to be summed,
+        and then maps the operator product by tensoring the individually mapped operators.
+
+        Args:
+            mixed_op: Operator to map.
+            register_length: UNUSED.
+        """
+
+        if register_length is not None:
+            LOGGER.info("Argument register length = %s was ignored.", register_length)
+        mapped_op: SparsePauliOp = self._distribute_map(mixed_op.data)
+
+        return mapped_op
+
+    def map(
+        self,
+        mixed_ops: MixedOp | ListOrDictType[MixedOp],
+        *,
+        register_length: int | None = None,
+    ) -> SparsePauliOp | ListOrDictType[SparsePauliOp]:
+        """Maps a second quantized operator or a list, dict of second quantized operators based on
+        the current mapper.
+
+        Args:
+            mixed_ops: A second quantized operator, or list thereof.
+            register_length: when provided, this will be used to overwrite the ``register_length``
+                attribute of the ``SparseLabelOp`` being mapped. This is possible because the
+                ``register_length`` is considered a lower bound in a ``SparseLabelOp``.
+
+        Returns:
+            A qubit operator in the form of a ``SparsePauliOp``, or list (resp. dict) thereof if a
+            list (resp. dict) of second quantized operators was supplied.
+        """
+        wrapped_second_q_ops, wrapped_type = _ListOrDict.wrap(mixed_ops)
+
+        qubit_ops: _ListOrDict = _ListOrDict()
+        for name, second_q_op in iter(wrapped_second_q_ops):
+            qubit_ops[name] = self._map_single(second_q_op, register_length=register_length)
+
+        returned_ops = qubit_ops.unwrap(wrapped_type)
+        # Note the output of the mapping will never be None for standard mappers other than the
+        # TaperedQubitMapper.
+        return returned_ops
diff --git a/qiskit_nature/second_q/operators/__init__.py b/qiskit_nature/second_q/operators/__init__.py
index 5ac9b844e..a66d7af1d 100644
--- a/qiskit_nature/second_q/operators/__init__.py
+++ b/qiskit_nature/second_q/operators/__init__.py
@@ -30,6 +30,7 @@
    VibrationalIntegrals
    PolynomialTensor
    Tensor
+   MixedOp
 
 Modules
 -------
@@ -51,6 +52,7 @@
 from .polynomial_tensor import PolynomialTensor
 from .sparse_label_op import SparseLabelOp
 from .tensor import Tensor
+from .mixed_op import MixedOp
 
 __all__ = [
     "ElectronicIntegrals",
@@ -62,4 +64,5 @@
     "PolynomialTensor",
     "SparseLabelOp",
     "Tensor",
+    "MixedOp",
 ]
diff --git a/qiskit_nature/second_q/operators/mixed_op.py b/qiskit_nature/second_q/operators/mixed_op.py
new file mode 100644
index 000000000..dbf70f8f9
--- /dev/null
+++ b/qiskit_nature/second_q/operators/mixed_op.py
@@ -0,0 +1,201 @@
+# This code is part of a Qiskit project.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""The Mixed Operator class."""
+
+from __future__ import annotations
+
+import itertools
+from copy import deepcopy
+from typing import Callable
+
+from qiskit.quantum_info.operators.mixins import LinearMixin
+
+from .sparse_label_op import SparseLabelOp
+
+
+class MixedOp(LinearMixin):
+    """Mixed operator.
+
+    A ``MixedOp`` represents a weighted sum of products of fermionic/bosonic operators potentially
+    acting on different "local" Hilbert spaces. The terms to be summed are encoded in a dictionary,
+    where each operator product is identified by its key, a tuple of string specifying the names of the
+    local Hilbert spaces on which it acts, and by its value, a list of tuple (corresponding to a sum of
+    operators acting on the same composite Hilbert space) where each tuple encodes the coupling
+    coefficient and the operators themselves (that might also have coefficients associated with them).
+
+
+    **Initialization**
+
+    A ``MixedOp`` is initialized with a dictionary, mapping terms to their respective
+    coefficients:
+
+    .. code-block:: python
+
+        from qiskit_nature.second_q.operators import FermionicOp, SpinOp, MixedOp
+
+        fop1 = FermionicOp({"+_0 -_0": 1}) # Acting on Hilbert space "h1"
+        sop1 = SpinOp({"X_0 Y_0": 1}, num_spins=1) # Acting on Hilbert space "s1"
+
+        mop1 = MixedOp({("h1",): [(5.0, fop1)]}) # 5.0 * fop1
+        mop2 = MixedOp(
+            {
+                ("h1", "s1"): [(3, fop1, sop1)],
+                ("s1",): [(2, sop1)],
+            }
+        ) # 3*(fop1 @ sop1) + 2*(sop1)
+
+
+    **Algebra**
+
+    This class supports the following basic arithmetic operations: addition, subtraction, scalar
+    multiplication, operator multiplication.
+    For example,
+
+    Addition
+
+    .. code-block:: python
+
+        fop1 = FermionicOp({"+_0 -_0": 1}) # Acting on Hilbert space "h1"
+        sop1 = SpinOp({"X_0 Y_0": 1}, num_spins=1) # Acting on Hilbert space "s1"
+        MixedOp({("h1",): [(5.0, fop1)]}) + MixedOp({("s1",): [(6.0, sop1)]})
+        # MixedOp({("h1",): [(5.0, fop1)], ("s1",): [(6.0, sop1)]})
+
+    Scalar multiplication
+
+    .. code-block:: python
+
+        fop1 = FermionicOp({"+_0 -_0": 1})
+        0.5 * MixedOp({("h1",): [(5.0, fop1)]})
+
+    Operator multiplication
+
+    .. code-block:: python
+
+        fop1 = FermionicOp({"+_0 -_0": 1}) # Acting on Hilbert space "h1"
+        sop1 = SpinOp({"X_0 Y_0": 1}, num_spins=1) # Acting on Hilbert space "s1"
+        MixedOp({("h1",): [(5.0, fop1)]}) @ MixedOp({("s1",): [(6.0, sop1)]})
+        # MixedOp({("h1", "s1"): [(30.0, fop1, sop1)]})
+
+    """
+
+    def __init__(self, data: dict[tuple[str], list[tuple[float, SparseLabelOp]]]):
+        self.data = deepcopy(data)
+
+    def __repr__(self, indentation_level=0) -> str:
+        out_str = "Mixed Op\n"
+        out_str += f"Nb terms = {len(self.data)}\n"
+        for active_indices, oplist in self.data.items():
+            for op_tuple in oplist:
+                coef, active_operators = op_tuple[0], op_tuple[1:]
+                out_str += f"- Coefficient: {coef:.02f}\n"
+
+                for index, op in zip(active_indices, active_operators):
+                    out_str += "" * indentation_level + f"{index}: {repr(op)}\n"
+
+                out_str += "\n"
+
+        return out_str
+
+    @staticmethod
+    def _tuple_prod(tup1: tuple[int, ...], tup2: tuple) -> tuple[float, ...]:
+        """Implements the composition of operator tuples representing tensor products of operators."""
+        new_coeff = tup1[0] * tup2[0]
+        new_op_tuple = tup1[1:] + tup2[1:]
+        return (new_coeff,) + new_op_tuple
+
+    @staticmethod
+    def _tuple_multiply(tup: tuple[int, ...], coef: float) -> tuple[float, ...]:
+        """Implements the dilation by a coefficient of an operator tuple representing a tensor product
+        of operators."""
+        new_coeff = tup[0] * coef
+        return (new_coeff,) + tup[1:]
+
+    @classmethod
+    def _distribute_on_tuples(
+        cls, method: Callable, op_left: MixedOp, op_right: MixedOp = None, **kwargs
+    ) -> MixedOp:
+        """Implements the distributions of a method to the tuples of operators representing the product
+        of operators."""
+        new_op_data: dict = {}
+        if op_right is None:
+            # Distribute method over all tuples.
+            for key, op_tuple_list in op_left.data.items():
+                new_op_data[key] = [method(op_tuple, **kwargs) for op_tuple in op_tuple_list]
+        else:
+            # Distribute method over all combinations of tuples from the first and second operators.
+            for (key1, op_tuple_list1), (key2, op_tuple_list2) in itertools.product(
+                op_left.data.items(), op_right.data.items()
+            ):
+                new_op_data[key1 + key2] = [
+                    method(op_tuple1, op_tuple2, **kwargs)
+                    for (op_tuple1, op_tuple2) in itertools.product(op_tuple_list1, op_tuple_list2)
+                ]
+
+        return MixedOp(new_op_data)
+
+    def _multiply(self, other: float) -> MixedOp:
+        """Return Operator multiplication of self and other.
+
+        Args:
+            other: the second ``MixedOp`` to multiply to the first.
+            qargs: UNUSED.
+
+        Returns:
+            The new multiplied ``MixedOp``.
+        """
+        return MixedOp._distribute_on_tuples(MixedOp._tuple_multiply, op_left=self, coef=other)
+
+    def _add(self, other: MixedOp, qargs: None = None) -> MixedOp:
+        """Return Operator addition of self and other.
+
+        Args:
+            other: the second ``MixedOp`` to add to the first.
+            qargs: UNUSED.
+
+        Returns:
+            The new summed ``MixedOp``.
+        """
+
+        sum_op = MixedOp(self.data)  # deepcopy
+        for key in other.data.keys():
+            # If the key for the composite Hilbert space already exists in the dictionary, then the
+            # addition is performed by appending the new operator to the corresponding list.
+            # Otherwise, the addition is performed by adding a new pair key, value to the dictionary.
+            if key in sum_op.data.keys():
+                sum_op.data[key] += other.data[key]
+            else:
+                sum_op.data[key] = other.data[key]
+        return sum_op
+
+    @classmethod
+    def compose(cls, op_left: MixedOp, op_right: MixedOp) -> MixedOp:
+        """Returns Operator composition of self and other.
+
+        Args:
+            op_left: left MixedOp to tensor.
+            op_right: right MixedOp to tensor.
+
+        Returns:
+            The tensor product of left with right.
+        """
+
+        # Lazy composition without applying the products.
+        return MixedOp._distribute_on_tuples(
+            MixedOp._tuple_prod, op_left=op_left, op_right=op_right
+        )
+
+    def __eq__(self, other):
+        if not isinstance(other, self.__class__):
+            return False
+
+        return self.data == other.data
diff --git a/test/second_q/mappers/test_mixed_mapper.py b/test/second_q/mappers/test_mixed_mapper.py
new file mode 100644
index 000000000..e3ec5a5c0
--- /dev/null
+++ b/test/second_q/mappers/test_mixed_mapper.py
@@ -0,0 +1,108 @@
+# This code is part of a Qiskit project.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+""" Test Mixed Mapper """
+
+import unittest
+
+
+from test import QiskitNatureTestCase
+
+from ddt import ddt, data, unpack
+import numpy as np
+
+from qiskit.quantum_info import SparsePauliOp
+from qiskit_nature.second_q.operators import BosonicOp, FermionicOp, MixedOp
+from qiskit_nature.second_q.mappers import (
+    BosonicLinearMapper,
+    JordanWignerMapper,
+    MixedMapper,
+)
+
+
+@ddt
+class TestMixedMapper(QiskitNatureTestCase):
+    """Test Mixed Mapper"""
+
+    # Define some useful coefficients
+    sq_2 = np.sqrt(2)
+
+    bos_op1 = BosonicOp({"+_0": 1})
+    mapped_bos_op1 = SparsePauliOp(["XX", "YY", "YX", "XY"], coeffs=[0.25, 0.25, -0.25j, 0.25j])
+
+    bos_op2 = BosonicOp({"-_0": 1})
+    mapped_bos_op2 = SparsePauliOp(["XX", "YY", "YX", "XY"], coeffs=[0.25, 0.25, 0.25j, -0.25j])
+
+    fer_op1 = FermionicOp({"+_0": 1}, num_spin_orbitals=1)
+    mapped_fer_op1 = SparsePauliOp.from_list([("X", 0.5), ("Y", -0.5j)])
+
+    fer_op2 = FermionicOp({"-_0": 1}, num_spin_orbitals=1)
+    mapped_fer_op2 = SparsePauliOp.from_list([("X", 0.5), ("Y", 0.5j)])
+
+    bos_op5 = BosonicOp({"+_0 -_0": 1})
+    bos_op6 = BosonicOp({"-_0 +_0": 1})
+
+    bos_mapper = BosonicLinearMapper(max_occupation=1)
+    fer_mapper = JordanWignerMapper()
+    mappers = {"b1": bos_mapper, "f1": fer_mapper}
+    hilbert_space_register_lengths = {"b1": 1, "f1": 1}
+    mix_mapper = MixedMapper(
+        mappers=mappers, hilbert_space_register_lengths=hilbert_space_register_lengths
+    )
+
+    @data(
+        (bos_op1, fer_op1, 2.0 + 1.0j),
+        (bos_op2, fer_op2, 3.0 + 2.0j),
+        (bos_op5, fer_op1, -4.0 - 3.0j),
+        (bos_op6, fer_op2, -5.0 + 4j),
+    )
+    @unpack
+    def test_relative_mapping(self, bos_op, fer_op, coef):
+        """Test the ``MixedOp`` mapping and compare to the composition of the mapped operators."""
+
+        mop = MixedOp({("b1", "f1"): [(coef, bos_op, fer_op)]})
+
+        target = coef * self.bos_mapper.map(bos_op).tensor(self.fer_mapper.map(fer_op))
+        test = self.mix_mapper.map(mop)
+        self.assertTrue(target.equiv(test))
+
+    @data(
+        (bos_op1, fer_op1, 2.0 + 1.0j),
+        (bos_op2, fer_op2, 3.0 + 2.0j),
+        (bos_op5, fer_op1, -4.0 - 3.0j),
+        (bos_op6, fer_op2, -5.0 + 4j),
+    )
+    @unpack
+    def test_map_list_or_dict(self, bos_op, fer_op, coef):
+        """Test the ``MixedOp`` mapping on list and dictionaries."""
+
+        mop_1 = MixedOp({("b1", "f1"): [(coef, bos_op, fer_op)]})
+        mop_2 = MixedOp({("b1", "f1"): [(coef, bos_op, fer_op)]})
+
+        calc_op_dict = self.mix_mapper.map({"A": mop_1, "B": mop_2})
+        calc_op_list = self.mix_mapper.map([mop_1, mop_2])
+
+        with self.subTest("Type dict"):
+            self.assertTrue(isinstance(calc_op_dict, dict))
+
+        with self.subTest("Type list"):
+            self.assertTrue(isinstance(calc_op_list, list))
+
+        with self.subTest("Value"):
+            self.assertTrue(
+                calc_op_list[0].equiv(calc_op_dict["A"])
+                and calc_op_list[1].equiv(calc_op_dict["B"])
+            )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/second_q/operators/test_mixed_op.py b/test/second_q/operators/test_mixed_op.py
new file mode 100644
index 000000000..45803d04b
--- /dev/null
+++ b/test/second_q/operators/test_mixed_op.py
@@ -0,0 +1,93 @@
+# This code is part of a Qiskit project.
+#
+# (C) Copyright IBM 2021, 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Test for MixedOp"""
+
+import unittest
+from test import QiskitNatureTestCase
+
+from qiskit_nature.second_q.operators import FermionicOp, MixedOp
+
+
+class TestFermionicOp(QiskitNatureTestCase):
+    """FermionicOp tests."""
+
+    op1 = FermionicOp({"+_0 -_0": 1})
+    op2 = FermionicOp({"-_0 +_0": 2})
+    mop1_h1 = MixedOp({("h1",): [(2.0, op1)]})
+    mop2_h1 = MixedOp({("h1",): [(3.0, op2)]})
+    mop2_h2 = MixedOp({("h2",): [(3.0, op2)]})
+    sumop2 = MixedOp({("h1",): [(1, FermionicOp({"+_0 -_0": 2, "-_0 +_0": 3}))]})
+
+    def test_neg(self):
+        """Test __neg__"""
+        minus_mop1 = -self.mop1_h1
+        target_v1 = MixedOp({("h1",): [(-2.0, self.op1)]})
+        self.assertEqual(minus_mop1, target_v1)
+
+    def test_mul(self):
+        """Test __mul__, and __rmul__"""
+        with self.subTest("rightmul"):
+            minus_mop1 = self.mop1_h1 * 2.0
+            target_v1 = MixedOp({("h1",): [(4.0, self.op1)]})
+            self.assertEqual(minus_mop1, target_v1)
+
+        with self.subTest("leftmul"):
+            minus_mop1 = (2.0 + 1.0j) * self.mop1_h1
+            target_v1 = MixedOp({("h1",): [((4.0 + 2.0j), self.op1)]})
+            self.assertEqual(minus_mop1, target_v1)
+
+    def test_div(self):
+        """Test __truediv__"""
+        fer_op = self.op1 / 2
+        target = FermionicOp({"+_0 -_0": 0.5}, num_spin_orbitals=1)
+        self.assertEqual(fer_op, target)
+
+    def test_add(self):
+        """Test __add__"""
+        with self.subTest("same hilbert space"):
+            sum_mop = self.mop1_h1 + self.mop2_h1
+            target = MixedOp({("h1",): [(2, self.op1), (3, self.op2)]})
+            self.assertEqual(sum_mop, target)
+
+        with self.subTest("different hilbert space"):
+            sum_mop = self.mop1_h1 + self.mop2_h2
+            target = MixedOp({("h1",): [(2.0, self.op1)], ("h2",): [(3.0, self.op2)]})
+            self.assertEqual(sum_mop, target)
+
+    def test_sub(self):
+        """Test __sub__"""
+        with self.subTest("same hilbert space"):
+            sum_mop = self.mop1_h1 - self.mop2_h1
+            target = MixedOp({("h1",): [(2, self.op1), (-3, self.op2)]})
+            self.assertEqual(sum_mop, target)
+
+        with self.subTest("different hilbert space"):
+            sum_mop = self.mop1_h1 - self.mop2_h2
+            target = MixedOp({("h1",): [(2.0, self.op1)], ("h2",): [(-3.0, self.op2)]})
+            self.assertEqual(sum_mop, target)
+
+    def test_compose(self):
+        """Test operator composition"""
+        with self.subTest("same hilbert spaces"):
+            composed_op = MixedOp.compose(self.mop1_h1, self.mop2_h1)
+            target = MixedOp({("h1", "h1"): [(6.0, self.op1, self.op2)]})
+            self.assertEqual(composed_op, target)
+
+        with self.subTest("different hilbert spaces"):
+            composed_op = MixedOp.compose(self.mop1_h1, self.mop2_h2)
+            target = MixedOp({("h1", "h2"): [(6.0, self.op1, self.op2)]})
+            self.assertEqual(composed_op, target)
+
+
+if __name__ == "__main__":
+    unittest.main()