diff --git a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index b3f04b4b3e33..e3459bae173e 100644
--- a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
+++ b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
@@ -46,6 +46,31 @@ class SparsePauliOp(LinearOp):
     using the :attr:`~SparsePauliOp.paulis` attribute. The coefficients
     are stored as a complex Numpy array vector and can be accessed using
     the :attr:`~SparsePauliOp.coeffs` attribute.
+
+    .. rubric:: Data type of coefficients
+
+    The default ``dtype`` of the internal ``coeffs`` Numpy array is ``complex128``.  Users can
+    configure this by passing ``np.ndarray`` with a different dtype.  For example, a parameterized
+    :class:`SparsePauliOp` can be made as follows:
+
+    .. code-block:: python
+
+        >>> import numpy as np
+        >>> from qiskit.circuit import ParameterVector
+        >>> from qiskit.quantum_info import SparsePauliOp
+
+        >>> SparsePauliOp(["II", "XZ"], np.array(ParameterVector("a", 2)))
+        SparsePauliOp(['II', 'XZ'],
+              coeffs=[ParameterExpression(1.0*a[0]), ParameterExpression(1.0*a[1])])
+
+    .. note::
+
+      Parameterized :class:`SparsePauliOp` does not support the following methods:
+
+      - ``to_matrix(sparse=True)`` since ``scipy.sparse`` cannot have objects as elements.
+      - ``to_operator()`` since :class:`~.quantum_info.Operator` does not support objects.
+      - ``sort``, ``argsort`` since :class:`.ParameterExpression` does not support comparison.
+      - ``equiv`` since :class:`.ParameterExpression`. cannot be converted into complex.
     """
 
     def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
@@ -86,10 +111,12 @@ def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
 
         pauli_list = PauliList(data.copy() if copy and hasattr(data, "copy") else data)
 
+        dtype = coeffs.dtype if isinstance(coeffs, np.ndarray) else complex
+
         if coeffs is None:
-            coeffs = np.ones(pauli_list.size, dtype=complex)
+            coeffs = np.ones(pauli_list.size, dtype=dtype)
         else:
-            coeffs = np.array(coeffs, copy=copy, dtype=complex)
+            coeffs = np.array(coeffs, copy=copy, dtype=dtype)
 
         if ignore_pauli_phase:
             # Fast path used in copy operations, where the phase of the PauliList is already known
@@ -101,7 +128,7 @@ def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
             # move the phase of `pauli_list` to `self._coeffs`
             phase = pauli_list._phase
             count_y = pauli_list._count_y()
-            self._coeffs = np.asarray((-1j) ** (phase - count_y) * coeffs, dtype=complex)
+            self._coeffs = np.asarray((-1j) ** (phase - count_y) * coeffs, dtype=coeffs.dtype)
             pauli_list._phase = np.mod(count_y, 4)
             self._pauli_list = pauli_list
 
@@ -132,9 +159,14 @@ def __eq__(self, other):
         """Entrywise comparison of two SparsePauliOp operators"""
         return (
             super().__eq__(other)
+            and self.coeffs.dtype == other.coeffs.dtype
             and self.coeffs.shape == other.coeffs.shape
-            and np.allclose(self.coeffs, other.coeffs)
             and self.paulis == other.paulis
+            and (
+                np.allclose(self.coeffs, other.coeffs)
+                if self.coeffs.dtype != object
+                else (self.coeffs == other.coeffs).all()
+            )
         )
 
     def equiv(self, other, atol: Optional[float] = None):
@@ -385,7 +417,19 @@ def simplify(self, atol=None, rtol=None):
             rtol = self.rtol
 
         # Filter non-zero coefficients
-        non_zero = np.logical_not(np.isclose(self.coeffs, 0, atol=atol, rtol=rtol))
+        if self.coeffs.dtype == object:
+
+            def to_complex(coeff):
+                if not hasattr(coeff, "sympify"):
+                    return coeff
+                sympified = coeff.sympify()
+                return complex(sympified) if sympified.is_Number else np.nan
+
+            non_zero = np.logical_not(
+                np.isclose([to_complex(x) for x in self.coeffs], 0, atol=atol, rtol=rtol)
+            )
+        else:
+            non_zero = np.logical_not(np.isclose(self.coeffs, 0, atol=atol, rtol=rtol))
         paulis_x = self.paulis.x[non_zero]
         paulis_z = self.paulis.z[non_zero]
         nz_coeffs = self.coeffs[non_zero]
@@ -398,16 +442,21 @@ def simplify(self, atol=None, rtol=None):
             # No zero operator or duplicate operator
             return self.copy()
 
-        coeffs = np.zeros(indexes.shape[0], dtype=complex)
+        coeffs = np.zeros(indexes.shape[0], dtype=self.coeffs.dtype)
         np.add.at(coeffs, inverses, nz_coeffs)
         # Delete zero coefficient rows
-        is_zero = np.isclose(coeffs, 0, atol=atol, rtol=rtol)
+        if self.coeffs.dtype == object:
+            is_zero = np.array(
+                [np.isclose(to_complex(coeff), 0, atol=atol, rtol=rtol) for coeff in coeffs]
+            )
+        else:
+            is_zero = np.isclose(coeffs, 0, atol=atol, rtol=rtol)
         # Check edge case that we deleted all Paulis
         # In this case we return an identity Pauli with a zero coefficient
         if np.all(is_zero):
             x = np.zeros((1, self.num_qubits), dtype=bool)
             z = np.zeros((1, self.num_qubits), dtype=bool)
-            coeffs = np.array([0j], dtype=complex)
+            coeffs = np.array([0j], dtype=self.coeffs.dtype)
         else:
             non_zero = np.logical_not(is_zero)
             non_zero_indexes = indexes[non_zero]
@@ -659,7 +708,7 @@ def from_operator(obj, atol=None, rtol=None):
         return SparsePauliOp(paulis, coeffs, copy=False)
 
     @staticmethod
-    def from_list(obj):
+    def from_list(obj, dtype=complex):
         """Construct from a list of Pauli strings and coefficients.
 
         For example, the 5-qubit Hamiltonian
@@ -677,6 +726,7 @@ def from_list(obj):
 
         Args:
             obj (Iterable[Tuple[str, complex]]): The list of 2-tuples specifying the Pauli terms.
+            dtype (type): The dtype of coeffs (Default complex).
 
         Returns:
             SparsePauliOp: The SparsePauliOp representation of the Pauli terms.
@@ -693,7 +743,7 @@ def from_list(obj):
         # determine the number of qubits
         num_qubits = len(obj[0][0])
 
-        coeffs = np.zeros(size, dtype=complex)
+        coeffs = np.zeros(size, dtype=dtype)
         labels = np.zeros(size, dtype=f"<U{num_qubits}")
         for i, item in enumerate(obj):
             labels[i] = item[0]
@@ -703,7 +753,7 @@ def from_list(obj):
         return SparsePauliOp(paulis, coeffs, copy=False)
 
     @staticmethod
-    def from_sparse_list(obj, num_qubits, do_checks=True):
+    def from_sparse_list(obj, num_qubits, do_checks=True, dtype=complex):
         """Construct from a list of local Pauli strings and coefficients.
 
         Each list element is a 3-tuple of a local Pauli string, indices where to apply it,
@@ -729,6 +779,7 @@ def from_sparse_list(obj, num_qubits, do_checks=True):
             obj (Iterable[Tuple[str, List[int], complex]]): The list 3-tuples specifying the Paulis.
             num_qubits (int): The number of qubits of the operator.
             do_checks (bool): The flag of checking if the input indices are not duplicated.
+            dtype (type): The dtype of coeffs (Default complex).
 
         Returns:
             SparsePauliOp: The SparsePauliOp representation of the Pauli terms.
@@ -744,7 +795,7 @@ def from_sparse_list(obj, num_qubits, do_checks=True):
         if size == 0:
             raise QiskitError("Input Pauli list is empty.")
 
-        coeffs = np.zeros(size, dtype=complex)
+        coeffs = np.zeros(size, dtype=dtype)
         labels = np.zeros(size, dtype=f"<U{num_qubits}")
 
         for i, (paulis, indices, coeff) in enumerate(obj):
@@ -781,7 +832,9 @@ def to_list(self, array=False):
         """
         # Dtype for a structured array with string labels and complex coeffs
         pauli_labels = self.paulis.to_labels(array=True)
-        labels = np.zeros(self.size, dtype=[("labels", pauli_labels.dtype), ("coeffs", "c16")])
+        labels = np.zeros(
+            self.size, dtype=[("labels", pauli_labels.dtype), ("coeffs", self.coeffs.dtype)]
+        )
         labels["labels"] = pauli_labels
         labels["coeffs"] = self.coeffs
         if array:
diff --git a/releasenotes/notes/operator-parameters-c81b7c05bffb740b.yaml b/releasenotes/notes/operator-parameters-c81b7c05bffb740b.yaml
new file mode 100644
index 000000000000..5d6edf00bad7
--- /dev/null
+++ b/releasenotes/notes/operator-parameters-c81b7c05bffb740b.yaml
@@ -0,0 +1,25 @@
+---
+features:
+  - |
+    :class:`.SparsePauliOp`\ s can now be constructed with coefficient arrays
+    that are general Python objects.  This is purely intended for use with Terra's
+    :class:`.ParameterExpression` objects; other objects may work, but do not
+    have first-class support.  Some :class:`.SparsePauliOp` methods (such as
+    conversion to other class representations) may not work when using ``object``
+    arrays, if the desired target cannot represent these general arrays.
+
+    For example, a :class:`.ParameterExpression` :class:`.SparsePauliOp` could
+    be constructed by::
+
+        import numpy as np
+        from qiskit.circuit import Parameter
+        from qiskit.quantum_info import SparsePauliOp
+
+        print(SparsePauliOp(["II", "XZ"], np.array([Parameter("a"), Parameter("b")])))
+
+    which gives
+
+    .. code-block:: text
+
+        SparsePauliOp(['II', 'XZ'],
+              coeffs=[ParameterExpression(1.0*a), ParameterExpression(1.0*b)])
diff --git a/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py b/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py
index e73777b4815d..bf58220d1d63 100644
--- a/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py
+++ b/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py
@@ -20,6 +20,7 @@
 from ddt import ddt
 
 from qiskit import QiskitError
+from qiskit.circuit import Parameter, ParameterVector
 from qiskit.quantum_info.operators import Operator, Pauli, PauliList, PauliTable, SparsePauliOp
 from qiskit.test import QiskitTestCase
 
@@ -83,6 +84,18 @@ def test_pauli_list_init(self):
             np.testing.assert_array_equal(spp_op.coeffs, [1, 2, 3, -4, 5j, -6j])
             paulis.phase = 0
             self.assertEqual(spp_op.paulis, paulis)
+        paulis = PauliList(labels)
+        with self.subTest(msg="with Parameterized coeffs"):
+            params = ParameterVector("params", 6)
+            coeffs = np.array(params)
+            spp_op = SparsePauliOp(paulis, coeffs)
+            target = coeffs.copy()
+            target[3] *= -1
+            target[4] *= 1j
+            target[5] *= -1j
+            np.testing.assert_array_equal(spp_op.coeffs, target)
+            paulis.phase = 0
+            self.assertEqual(spp_op.paulis, paulis)
 
     def test_sparse_pauli_op_init(self):
         """Test SparsePauliOp initialization."""
@@ -156,6 +169,14 @@ def test_from_list(self):
         np.testing.assert_array_equal(spp_op.coeffs, coeffs)
         self.assertEqual(spp_op.paulis, PauliList(labels))
 
+    def test_from_list_parameters(self):
+        """Test from_list method with parameters."""
+        labels = ["XXZ", "IXI", "YZZ", "III"]
+        coeffs = ParameterVector("a", 4)
+        spp_op = SparsePauliOp.from_list(zip(labels, coeffs), dtype=object)
+        np.testing.assert_array_equal(spp_op.coeffs, coeffs)
+        self.assertEqual(spp_op.paulis, PauliList(labels))
+
     def test_from_index_list(self):
         """Test from_list method specifying the Paulis via indices."""
         expected_labels = ["XXZ", "IXI", "YIZ", "III"]
@@ -166,6 +187,18 @@ def test_from_index_list(self):
         np.testing.assert_array_equal(spp_op.coeffs, coeffs)
         self.assertEqual(spp_op.paulis, PauliList(expected_labels))
 
+    def test_from_index_list_parameters(self):
+        """Test from_list method specifying the Paulis via indices with paramteres."""
+        expected_labels = ["XXZ", "IXI", "YIZ", "III"]
+        paulis = ["XXZ", "X", "YZ", ""]
+        indices = [[2, 1, 0], [1], [2, 0], []]
+        coeffs = ParameterVector("a", 4)
+        spp_op = SparsePauliOp.from_sparse_list(
+            zip(paulis, indices, coeffs), num_qubits=3, dtype=object
+        )
+        np.testing.assert_array_equal(spp_op.coeffs, coeffs)
+        self.assertEqual(spp_op.paulis, PauliList(expected_labels))
+
     def test_from_index_list_endianness(self):
         """Test the construction from index list has the right endianness."""
         spp_op = SparsePauliOp.from_sparse_list([("ZX", [1, 4], 1)], num_qubits=5)
@@ -204,6 +237,16 @@ def to_matrix(self):
             target += coeff * pauli_mat(label)
         np.testing.assert_array_equal(spp_op.to_matrix(), target)
 
+    def to_matrix_parameters(self):
+        """Test to_matrix method for parameterized SparsePauliOp."""
+        labels = ["XI", "YZ", "YY", "ZZ"]
+        coeffs = ParameterVector("a", 4)
+        spp_op = SparsePauliOp(labels, coeffs)
+        target = np.zeros((4, 4), dtype=object)
+        for coeff, label in zip(coeffs, labels):
+            target += coeff * pauli_mat(label)
+        np.testing.assert_array_equal(spp_op.to_matrix(), target)
+
     def to_operator(self):
         """Test to_operator method."""
         labels = ["XI", "YZ", "YY", "ZZ"]
@@ -222,6 +265,14 @@ def to_list(self):
         target = list(zip(labels, coeffs))
         self.assertEqual(op.to_list(), target)
 
+    def to_list_parameters(self):
+        """Test to_operator method with paramters."""
+        labels = ["XI", "YZ", "YY", "ZZ"]
+        coeffs = np.array(ParameterVector("a", 4))
+        op = SparsePauliOp(labels, coeffs)
+        target = list(zip(labels, coeffs))
+        self.assertEqual(op.to_list(), target)
+
 
 class TestSparsePauliOpIteration(QiskitTestCase):
     """Tests for SparsePauliOp iterators class."""
@@ -234,6 +285,14 @@ def test_enumerate(self):
         for idx, i in enumerate(op):
             self.assertEqual(i, SparsePauliOp(labels[idx], coeffs[[idx]]))
 
+    def test_enumerate_parameters(self):
+        """Test enumerate with SparsePauliOp with parameters."""
+        labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
+        coeffs = np.array(ParameterVector("a", 6))
+        op = SparsePauliOp(labels, coeffs)
+        for idx, i in enumerate(op):
+            self.assertEqual(i, SparsePauliOp(labels[idx], coeffs[[idx]]))
+
     def test_iter(self):
         """Test iter with SparsePauliOp."""
         labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
@@ -242,6 +301,14 @@ def test_iter(self):
         for idx, i in enumerate(iter(op)):
             self.assertEqual(i, SparsePauliOp(labels[idx], coeffs[[idx]]))
 
+    def test_iter_parameters(self):
+        """Test iter with SparsePauliOp with parameters."""
+        labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
+        coeffs = np.array(ParameterVector("a", 6))
+        op = SparsePauliOp(labels, coeffs)
+        for idx, i in enumerate(iter(op)):
+            self.assertEqual(i, SparsePauliOp(labels[idx], coeffs[[idx]]))
+
     def test_label_iter(self):
         """Test SparsePauliOp label_iter method."""
         labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
@@ -250,6 +317,14 @@ def test_label_iter(self):
         for idx, i in enumerate(op.label_iter()):
             self.assertEqual(i, (labels[idx], coeffs[idx]))
 
+    def test_label_iter_parameters(self):
+        """Test SparsePauliOp label_iter method with parameters."""
+        labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
+        coeffs = np.array(ParameterVector("a", 6))
+        op = SparsePauliOp(labels, coeffs)
+        for idx, i in enumerate(op.label_iter()):
+            self.assertEqual(i, (labels[idx], coeffs[idx]))
+
     def test_matrix_iter(self):
         """Test SparsePauliOp dense matrix_iter method."""
         labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
@@ -258,6 +333,14 @@ def test_matrix_iter(self):
         for idx, i in enumerate(op.matrix_iter()):
             np.testing.assert_array_equal(i, coeffs[idx] * pauli_mat(labels[idx]))
 
+    def test_matrix_iter_parameters(self):
+        """Test SparsePauliOp dense matrix_iter method. with parameters"""
+        labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
+        coeffs = np.array(ParameterVector("a", 6))
+        op = SparsePauliOp(labels, coeffs)
+        for idx, i in enumerate(op.matrix_iter()):
+            np.testing.assert_array_equal(i, coeffs[idx] * pauli_mat(labels[idx]))
+
     def test_matrix_iter_sparse(self):
         """Test SparsePauliOp sparse matrix_iter method."""
         labels = ["III", "IXI", "IYY", "YIZ", "XYZ", "III"]
@@ -267,85 +350,116 @@ def test_matrix_iter_sparse(self):
             np.testing.assert_array_equal(i.toarray(), coeffs[idx] * pauli_mat(labels[idx]))
 
 
+def bind_parameters_to_one(array):
+    """Bind parameters to one. The purpose of using this method is to bind some value and
+    use ``assert_allclose``, since it is impossible to verify equivalence in the case of
+    numerical errors with parameters existing.
+    """
+
+    def bind_one(a):
+        parameters = a.parameters
+        return complex(a.bind(dict(zip(parameters, [1] * len(parameters)))))
+
+    return np.vectorize(bind_one, otypes=[complex])(array)
+
+
 @ddt
 class TestSparsePauliOpMethods(QiskitTestCase):
     """Tests for SparsePauliOp operator methods."""
 
     RNG = np.random.default_rng(1994)
 
-    def random_spp_op(self, num_qubits, num_terms):
+    def setUp(self):
+        super().setUp()
+
+        self.parameter_names = (f"param_{x}" for x in it.count())
+
+    def random_spp_op(self, num_qubits, num_terms, use_parameters=False):
         """Generate a pseudo-random SparsePauliOp"""
-        coeffs = self.RNG.uniform(-1, 1, size=num_terms) + 1j * self.RNG.uniform(
-            -1, 1, size=num_terms
-        )
+        if use_parameters:
+            coeffs = np.array(ParameterVector(next(self.parameter_names), num_terms))
+        else:
+            coeffs = self.RNG.uniform(-1, 1, size=num_terms) + 1j * self.RNG.uniform(
+                -1, 1, size=num_terms
+            )
         labels = [
             "".join(self.RNG.choice(["I", "X", "Y", "Z"], size=num_qubits))
             for _ in range(num_terms)
         ]
         return SparsePauliOp(labels, coeffs)
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_conjugate(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_conjugate(self, num_qubits, use_parameters):
         """Test conjugate method for {num_qubits}-qubits."""
-        spp_op = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op).conjugate()
+        spp_op = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op.to_matrix().conjugate()
         op = spp_op.conjugate()
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        np.testing.assert_array_equal(value, target)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_transpose(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_transpose(self, num_qubits, use_parameters):
         """Test transpose method for {num_qubits}-qubits."""
-        spp_op = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op).transpose()
+        spp_op = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op.to_matrix().transpose()
         op = spp_op.transpose()
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        np.testing.assert_array_equal(value, target)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_adjoint(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_adjoint(self, num_qubits, use_parameters):
         """Test adjoint method for {num_qubits}-qubits."""
-        spp_op = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op).adjoint()
+        spp_op = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op.to_matrix().transpose().conjugate()
         op = spp_op.adjoint()
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        np.testing.assert_array_equal(value, target)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_compose(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_compose(self, num_qubits, use_parameters):
         """Test {num_qubits}-qubit compose methods."""
-        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits)
-        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op1).compose(Operator(spp_op2))
+        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op2.to_matrix() @ spp_op1.to_matrix()
 
         op = spp_op1.compose(spp_op2)
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
         op = spp_op1 & spp_op2
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_dot(self, num_qubits):
-        """Test {num_qubits}-qubit compose methods."""
-        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits)
-        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op1).dot(Operator(spp_op2))
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_dot(self, num_qubits, use_parameters):
+        """Test {num_qubits}-qubit dot methods."""
+        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op1.to_matrix() @ spp_op2.to_matrix()
 
         op = spp_op1.dot(spp_op2)
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
         op = spp_op1 @ spp_op2
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
     @combine(num_qubits=[1, 2, 3])
@@ -379,48 +493,60 @@ def test_qargs_dot(self, num_qubits):
         self.assertEqual(value, target)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits1=[1, 2, 3], num_qubits2=[1, 2, 3])
-    def test_tensor(self, num_qubits1, num_qubits2):
+    @combine(num_qubits1=[1, 2, 3], num_qubits2=[1, 2, 3], use_parameters=[True, False])
+    def test_tensor(self, num_qubits1, num_qubits2, use_parameters):
         """Test tensor method for {num_qubits1} and {num_qubits2} qubits."""
-        spp_op1 = self.random_spp_op(num_qubits1, 2**num_qubits1)
-        spp_op2 = self.random_spp_op(num_qubits2, 2**num_qubits2)
-        target = Operator(spp_op1).tensor(Operator(spp_op2))
+        spp_op1 = self.random_spp_op(num_qubits1, 2**num_qubits1, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits2, 2**num_qubits2, use_parameters)
+        target = np.kron(spp_op1.to_matrix(), spp_op2.to_matrix())
         op = spp_op1.tensor(spp_op2)
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits1=[1, 2, 3], num_qubits2=[1, 2, 3])
-    def test_expand(self, num_qubits1, num_qubits2):
+    @combine(num_qubits1=[1, 2, 3], num_qubits2=[1, 2, 3], use_parameters=[True, False])
+    def test_expand(self, num_qubits1, num_qubits2, use_parameters):
         """Test expand method for {num_qubits1} and {num_qubits2} qubits."""
-        spp_op1 = self.random_spp_op(num_qubits1, 2**num_qubits1)
-        spp_op2 = self.random_spp_op(num_qubits2, 2**num_qubits2)
-        target = Operator(spp_op1).expand(Operator(spp_op2))
+        spp_op1 = self.random_spp_op(num_qubits1, 2**num_qubits1, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits2, 2**num_qubits2, use_parameters)
+        target = np.kron(spp_op2.to_matrix(), spp_op1.to_matrix())
         op = spp_op1.expand(spp_op2)
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_add(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_add(self, num_qubits, use_parameters):
         """Test + method for {num_qubits} qubits."""
-        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits)
-        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op1) + Operator(spp_op2)
+        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op1.to_matrix() + spp_op2.to_matrix()
         op = spp_op1 + spp_op2
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_sub(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_sub(self, num_qubits, use_parameters):
         """Test + method for {num_qubits} qubits."""
-        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits)
-        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op1) - Operator(spp_op2)
+        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        target = spp_op1.to_matrix() - spp_op2.to_matrix()
         op = spp_op1 - spp_op2
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value = op.to_matrix()
+        if use_parameters:
+            value = bind_parameters_to_one(value)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
     @combine(num_qubits=[1, 2, 3])
@@ -447,29 +573,44 @@ def test_sub_qargs(self, num_qubits):
         self.assertEqual(value, target)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3], value=[0, 1, 1j, -3 + 4.4j, np.int64(2)])
-    def test_mul(self, num_qubits, value):
+    @combine(num_qubits=[1, 2, 3], value=[0, 1, 1j, -3 + 4.4j, np.int64(2)], param=[None, "a"])
+    def test_mul(self, num_qubits, value, param):
         """Test * method for {num_qubits} qubits and value {value}."""
-        spp_op = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = value * Operator(spp_op)
+        spp_op = self.random_spp_op(num_qubits, 2**num_qubits, param)
+        target = value * spp_op.to_matrix()
         op = value * spp_op
-        value_mat = op.to_operator()
-        self.assertEqual(value_mat, target)
+        value_mat = op.to_matrix()
+        if value != 0 and param is not None:
+            value_mat = bind_parameters_to_one(value_mat)
+            target = bind_parameters_to_one(target)
+        if value == 0:
+            np.testing.assert_array_equal(value_mat, target.astype(complex))
+        else:
+            np.testing.assert_allclose(value_mat, target)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
-        target = Operator(spp_op) * value
+        target = spp_op.to_matrix() * value
         op = spp_op * value
-        value_mat = op.to_operator()
-        self.assertEqual(value_mat, target)
+        value_mat = op.to_matrix()
+        if value != 0 and param is not None:
+            value_mat = bind_parameters_to_one(value_mat)
+            target = bind_parameters_to_one(target)
+        if value == 0:
+            np.testing.assert_array_equal(value_mat, target.astype(complex))
+        else:
+            np.testing.assert_allclose(value_mat, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
-    @combine(num_qubits=[1, 2, 3], value=[1, 1j, -3 + 4.4j])
-    def test_div(self, num_qubits, value):
+    @combine(num_qubits=[1, 2, 3], value=[1, 1j, -3 + 4.4j], param=[None, "a"])
+    def test_div(self, num_qubits, value, param):
         """Test / method for {num_qubits} qubits and value {value}."""
-        spp_op = self.random_spp_op(num_qubits, 2**num_qubits)
-        target = Operator(spp_op) / value
+        spp_op = self.random_spp_op(num_qubits, 2**num_qubits, param)
+        target = spp_op.to_matrix() / value
         op = spp_op / value
-        value = op.to_operator()
-        self.assertEqual(value, target)
+        value_mat = op.to_matrix()
+        if param is not None:
+            value_mat = bind_parameters_to_one(value_mat)
+            target = bind_parameters_to_one(target)
+        np.testing.assert_allclose(value_mat, target, atol=1e-8)
         np.testing.assert_array_equal(op.paulis.phase, np.zeros(op.size))
 
     def test_simplify(self):
@@ -510,6 +651,19 @@ def test_simplify_zero(self, num_qubits):
         np.testing.assert_array_equal(zero_op.paulis.phase, np.zeros(zero_op.size))
         np.testing.assert_array_equal(simplified_op.paulis.phase, np.zeros(simplified_op.size))
 
+    def test_simplify_parameters(self):
+        """Test simplify methods for parameterized SparsePauliOp."""
+        a = Parameter("a")
+        coeffs = np.array([a, -a, 0, a, a, a, 2 * a])
+        labels = ["IXI", "IXI", "ZZZ", "III", "III", "XXX", "XXX"]
+        spp_op = SparsePauliOp(labels, coeffs)
+        simplified_op = spp_op.simplify()
+        target_coeffs = np.array([2 * a, 3 * a])
+        target_labels = ["III", "XXX"]
+        target_op = SparsePauliOp(target_labels, target_coeffs)
+        self.assertEqual(simplified_op, target_op)
+        np.testing.assert_array_equal(simplified_op.paulis.phase, np.zeros(simplified_op.size))
+
     def test_sort(self):
         """Test sort method."""
         with self.assertRaises(QiskitError):
@@ -656,14 +810,22 @@ def test_chop_all(self):
         expected = SparsePauliOp(["I"], coeffs=[0.0])
         self.assertEqual(simplified, expected)
 
-    @combine(num_qubits=[1, 2, 3, 4], num_ops=[1, 2, 3, 4])
-    def test_sum(self, num_qubits, num_ops):
+    @combine(num_qubits=[1, 2, 3, 4], num_ops=[1, 2, 3, 4], param=[None, "a"])
+    def test_sum(self, num_qubits, num_ops, param):
         """Test sum method for {num_qubits} qubits with {num_ops} operators."""
-        ops = [self.random_spp_op(num_qubits, 2**num_qubits) for _ in range(num_ops)]
+        ops = [
+            self.random_spp_op(
+                num_qubits, 2**num_qubits, param if param is None else f"{param}_{i}"
+            )
+            for i in range(num_ops)
+        ]
         sum_op = SparsePauliOp.sum(ops)
-        value = Operator(sum_op)
-        target_operator = sum((Operator(op) for op in ops[1:]), Operator(ops[0]))
-        self.assertEqual(value, target_operator)
+        value = sum_op.to_matrix()
+        target_operator = sum((op.to_matrix() for op in ops[1:]), ops[0].to_matrix())
+        if param is not None:
+            value = bind_parameters_to_one(value)
+            target_operator = bind_parameters_to_one(target_operator)
+        np.testing.assert_allclose(value, target_operator, atol=1e-8)
         target_spp_op = sum((op for op in ops[1:]), ops[0])
         self.assertEqual(sum_op, target_spp_op)
         np.testing.assert_array_equal(sum_op.paulis.phase, np.zeros(sum_op.size))
@@ -678,12 +840,12 @@ def test_sum_error(self):
         with self.assertRaises(QiskitError):
             SparsePauliOp.sum([1, 2])
 
-    @combine(num_qubits=[1, 2, 3, 4])
-    def test_eq(self, num_qubits):
+    @combine(num_qubits=[1, 2, 3, 4], use_parameters=[True, False])
+    def test_eq(self, num_qubits, use_parameters):
         """Test __eq__ method for {num_qubits} qubits."""
-        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits)
-        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits)
-        spp_op3 = self.random_spp_op(num_qubits, 2**num_qubits)
+        spp_op1 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        spp_op2 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
+        spp_op3 = self.random_spp_op(num_qubits, 2**num_qubits, use_parameters)
         zero = spp_op3 - spp_op3
         self.assertEqual(spp_op1, spp_op1)
         self.assertEqual(spp_op2, spp_op2)
@@ -739,7 +901,8 @@ def test_eq_equiv(self):
             self.assertNotEqual(spp_op1, spp_op2)
             self.assertTrue(spp_op1.equiv(spp_op2))
 
-    def test_group_commuting(self):
+    @combine(parameterized=[True, False])
+    def test_group_commuting(self, parameterized):
         """Test general grouping commuting operators"""
 
         def commutes(left: Pauli, right: Pauli) -> bool:
@@ -747,8 +910,11 @@ def commutes(left: Pauli, right: Pauli) -> bool:
 
         input_labels = ["IX", "IY", "IZ", "XX", "YY", "ZZ", "XY", "YX", "ZX", "ZY", "XZ", "YZ"]
         np.random.shuffle(input_labels)
-        coefs = np.random.random(len(input_labels)) + np.random.random(len(input_labels)) * 1j
-        sparse_pauli_list = SparsePauliOp(input_labels, coefs)
+        if parameterized:
+            coeffs = np.array(ParameterVector("a", len(input_labels)))
+        else:
+            coeffs = np.random.random(len(input_labels)) + np.random.random(len(input_labels)) * 1j
+        sparse_pauli_list = SparsePauliOp(input_labels, coeffs)
         groups = sparse_pauli_list.group_commuting()
         # checking that every input Pauli in sparse_pauli_list is in a group in the ouput
         output_labels = [pauli.to_label() for group in groups for pauli in group.paulis]
@@ -757,7 +923,7 @@ def commutes(left: Pauli, right: Pauli) -> bool:
         paulis_coeff_dict = dict(
             sum([list(zip(group.paulis.to_labels(), group.coeffs)) for group in groups], [])
         )
-        self.assertDictEqual(dict(zip(input_labels, coefs)), paulis_coeff_dict)
+        self.assertDictEqual(dict(zip(input_labels, coeffs)), paulis_coeff_dict)
 
         # Within each group, every operator commutes with every other operator.
         for group in groups: