diff --git a/qiskit/circuit/library/pauli_evolution.py b/qiskit/circuit/library/pauli_evolution.py
index 2711d81cbacf..dad68239f090 100644
--- a/qiskit/circuit/library/pauli_evolution.py
+++ b/qiskit/circuit/library/pauli_evolution.py
@@ -51,7 +51,7 @@ class PauliEvolutionGate(Gate):
     def __init__(
         self,
         operator,
-        time: Union[float, ParameterExpression] = 1.0,
+        time: Union[int, float, ParameterExpression] = 1.0,
         label: Optional[str] = None,
         synthesis: Optional[EvolutionSynthesis] = None,
     ) -> None:
@@ -79,10 +79,27 @@ def __init__(
         num_qubits = operator[0].num_qubits if isinstance(operator, list) else operator.num_qubits
         super().__init__(name=name, num_qubits=num_qubits, params=[time], label=label)
 
-        self.time = time
         self.operator = operator
         self.synthesis = synthesis
 
+    @property
+    def time(self) -> Union[float, ParameterExpression]:
+        """Return the evolution time as stored in the gate parameters.
+
+        Returns:
+            The evolution time.
+        """
+        return self.params[0]
+
+    @time.setter
+    def time(self, time: Union[float, ParameterExpression]) -> None:
+        """Set the evolution time.
+
+        Args:
+            time: The evolution time.
+        """
+        self.params = [time]
+
     def _define(self):
         """Unroll, where the default synthesis is matrix based."""
         self.definition = self.synthesis.synthesize(self)
@@ -90,6 +107,15 @@ def _define(self):
     def inverse(self) -> "PauliEvolutionGate":
         return PauliEvolutionGate(operator=self.operator, time=-self.time, synthesis=self.synthesis)
 
+    def validate_parameter(
+        self, parameter: Union[int, float, ParameterExpression]
+    ) -> Union[float, ParameterExpression]:
+        """Gate parameters should be int, float, or ParameterExpression"""
+        if isinstance(parameter, int):
+            parameter = float(parameter)
+
+        return super().validate_parameter(parameter)
+
 
 def _to_sparse_pauli_op(operator):
     """Cast the operator to a SparsePauliOp.
diff --git a/qiskit/circuit/qpy_serialization.py b/qiskit/circuit/qpy_serialization.py
index 3785f4d433d6..0e63e46a37ec 100644
--- a/qiskit/circuit/qpy_serialization.py
+++ b/qiskit/circuit/qpy_serialization.py
@@ -100,6 +100,66 @@
 by ``num_circuits`` in the file header). There is no padding between the
 circuits in the data.
 
+.. _version_3:
+
+Version 3
+=========
+
+Version 3 of the QPY format is identical to :ref:`version_2` except that it defines
+a struct format to represent a :class:`~qiskit.circuit.library.PauliEvolutionGate`
+natively in QPY. To accomplish this the :ref:`custom_definition` struct now supports
+a new type value ``'p'`` to represent a :class:`~qiskit.circuit.library.PauliEvolutionGate`.
+Enties in the custom instructions tables have unique name generated that start with the
+string ``"###PauliEvolutionGate_"`` followed by a uuid string. This gate name is reservered
+in QPY and if you have a custom :class:`~qiskit.circuit.Instruction` object with a definition
+set and that name prefix it will error. If it's of type ``'p'`` the data payload is defined
+as follows:
+
+.. _pauli_evo_qpy:
+
+PAULI_EVOLUTION
+---------------
+
+This represents the high level :class:`~qiskit.circuit.library.PauliEvolutionGate`
+
+.. code-block:: c
+
+    struct {
+        uint64_t operator_count;
+        _Bool standalone_op;
+        char time_type;
+        uint64_t time_size;
+        uint64_t synthesis_size;
+    }
+
+This is immediately followed by ``operator_count`` elements defined by the :ref:`pauli_sum_op`
+payload.  Following that we have ``time_size`` bytes representing the ``time`` attribute. If
+``standalone_op`` is ``True`` then there must only be a single operator. The
+encoding of these bytes is determined by the value of ``time_type``. Possible values of
+``time_type`` are ``'f'``, ``'p'``, and ``'e'``. If ``time_type`` is ``'f'`` it's a double,
+``'p'`` defines a :class:`~qiskit.circuit.Parameter` object  which is represented by a
+:ref:`param_struct`, ``e`` defines a :class:`~qiskit.circuit.ParameterExpression` object
+(that's not a :class:`~qiskit.circuit.Parameter`) which is represented by a :ref:`param_expr`.
+Following that is ``synthesis_size`` bytes which is a utf8 encoded json payload representing
+the :class:`.EvolutionSynthesis` class used by the gate.
+
+.. _pauli_sum_op:
+
+SPARSE_PAULI_OP_LIST_ELEM
+-------------------------
+
+This represents an instance of :class:`.PauliSumOp`.
+
+
+.. code-block:: c
+
+    struct {
+        uint32_t pauli_op_size;
+    }
+
+which is immediately followed by ``pauli_op_size`` bytes which are .npy format [#f2]_
+data which represents the :class:`~qiskit.quantum_info.SparsePauliOp`.
+
 .. _version_2:
 
 Version 2
@@ -213,6 +273,8 @@
 ``qr`` would have ``standalone`` set to ``False``.
 
 
+.. _custom_definition:
+
 CUSTOM_DEFINITIONS
 ------------------
 
@@ -244,7 +306,10 @@
 If ``custom_definition`` is ``True`` that means that the immediately following
 ``size`` bytes contains a QPY circuit data which can be used for the custom
 definition of that gate. If ``custom_definition`` is ``False`` then the
-instruction can be considered opaque (ie no definition).
+instruction can be considered opaque (ie no definition). The ``type`` field
+determines what type of object will get created with the custom definition.
+If it's ``'g'`` it will be a :class:`~qiskit.circuit.Gate` object, ``'i'``
+it will be a :class:`~qiskit.circuit.Instruction` object.
 
 INSTRUCTIONS
 ------------
@@ -316,6 +381,8 @@ class if it's defined in Qiskit. Otherwise it falls back to the custom
 (see below), and ``'n'`` represents an object from numpy (either an ``ndarray``
 or a numpy type) which means the data is .npy format [#f2]_ data.
 
+.. _param_struct:
+
 PARAMETER
 ---------
 
@@ -407,6 +474,8 @@ class if it's defined in Qiskit. Otherwise it falls back to the custom
 from qiskit.extensions import quantum_initializer
 from qiskit.version import __version__
 from qiskit.exceptions import QiskitError
+from qiskit.quantum_info.operators import SparsePauliOp
+from qiskit.synthesis import evolution as evo_synth
 
 try:
     import symengine
@@ -518,6 +587,17 @@ class if it's defined in Qiskit. Otherwise it falls back to the custom
 COMPLEX = namedtuple("COMPLEX", ["real", "imag"])
 COMPLEX_PACK = "!dd"
 COMPLEX_SIZE = struct.calcsize(COMPLEX_PACK)
+# Pauli Evolution Gate
+PAULI_EVOLUTION_DEF = namedtuple(
+    "PAULI_EVOLUTION_DEF",
+    ["operator_size", "standalone_op", "time_type", "time_size", "synth_method_size"],
+)
+PAULI_EVOLUTION_DEF_PACK = "!Q?1cQQ"
+PAULI_EVOLUTION_DEF_SIZE = struct.calcsize(PAULI_EVOLUTION_DEF_PACK)
+# SparsePauliOp List
+SPARSE_PAULI_OP_LIST_ELEM = namedtuple("SPARSE_PAULI_OP_LIST_ELEMENT", ["size"])
+SPARSE_PAULI_OP_LIST_ELEM_PACK = "!Q"
+SPARSE_PAULI_OP_LIST_ELEM_SIZE = struct.calcsize(SPARSE_PAULI_OP_LIST_ELEM_PACK)
 
 
 def _read_header_v2(file_obj):
@@ -751,6 +831,8 @@ def _parse_custom_instruction(custom_instructions, gate_name, params):
     elif type_str == "g":
         inst_obj = Gate(gate_name, num_qubits, params)
         inst_obj.definition = definition
+    elif type_str == "p":
+        inst_obj = definition
     else:
         raise ValueError("Invalid custom instruction type '%s'" % type_str)
     return inst_obj
@@ -779,7 +861,10 @@ def _read_custom_instructions(file_obj, version):
             definition_circuit = None
             if has_custom_definition:
                 definition_buffer = io.BytesIO(file_obj.read(size))
-                definition_circuit = _read_circuit(definition_buffer, version)
+                if version < 3 or not name.startswith(r"###PauliEvolutionGate_"):
+                    definition_circuit = _read_circuit(definition_buffer, version)
+                elif name.startswith(r"###PauliEvolutionGate_"):
+                    definition_circuit = _read_pauli_evolution_gate(definition_buffer)
             custom_instructions[name] = (type_str, num_qubits, num_clbits, definition_circuit)
     return custom_instructions
 
@@ -842,12 +927,16 @@ def _write_instruction(file_obj, instruction_tuple, custom_instructions, index_m
         )
         or gate_class_name == "Gate"
         or gate_class_name == "Instruction"
-        or isinstance(instruction_tuple[0], (library.BlueprintCircuit, library.PauliEvolutionGate))
+        or isinstance(instruction_tuple[0], library.BlueprintCircuit)
     ):
         if instruction_tuple[0].name not in custom_instructions:
             custom_instructions[instruction_tuple[0].name] = instruction_tuple[0]
         gate_class_name = instruction_tuple[0].name
 
+    elif isinstance(instruction_tuple[0], library.PauliEvolutionGate):
+        gate_class_name = r"###PauliEvolutionGate_" + str(uuid.uuid4())
+        custom_instructions[gate_class_name] = instruction_tuple[0]
+
     has_condition = False
     condition_register = b""
     condition_value = 0
@@ -936,8 +1025,97 @@ def _write_instruction(file_obj, instruction_tuple, custom_instructions, index_m
         container.close()
 
 
+def _write_pauli_evolution_gate(file_obj, evolution_gate):
+    operator_list = evolution_gate.operator
+    standalone = False
+    if not isinstance(operator_list, list):
+        operator_list = [operator_list]
+        standalone = True
+    num_operators = len(operator_list)
+    pauli_data_buf = io.BytesIO()
+    for operator in operator_list:
+        with io.BytesIO() as element_buf:
+            with io.BytesIO() as buf:
+                pauli_list = operator.to_list(array=True)
+                np.save(buf, pauli_list)
+                data = buf.getvalue()
+            element_metadata = struct.pack(SPARSE_PAULI_OP_LIST_ELEM_PACK, len(data))
+            element_buf.write(element_metadata)
+            element_buf.write(data)
+            pauli_data_buf.write(element_buf.getvalue())
+    time = evolution_gate.time
+    if isinstance(time, float):
+        time_type = b"f"
+        time_data = struct.pack("!d", time)
+        time_size = struct.calcsize("!d")
+    elif isinstance(time, Parameter):
+        time_type = b"p"
+        with io.BytesIO() as buf:
+            _write_parameter(buf, time)
+            time_data = buf.getvalue()
+            time_size = len(time_data)
+    elif isinstance(time, ParameterExpression):
+        time_type = b"e"
+        with io.BytesIO() as buf:
+            _write_parameter_expression(buf, time)
+            time_data = buf.getvalue()
+            time_size = len(time_data)
+    else:
+        raise TypeError(f"Invalid time type {time} for PauliEvolutionGate")
+
+    synth_class = str(type(evolution_gate.synthesis).__name__)
+    settings_dict = evolution_gate.synthesis.settings
+    synth_data = json.dumps({"class": synth_class, "settings": settings_dict}).encode("utf8")
+    synth_size = len(synth_data)
+    pauli_evolution_raw = struct.pack(
+        PAULI_EVOLUTION_DEF_PACK, num_operators, standalone, time_type, time_size, synth_size
+    )
+    file_obj.write(pauli_evolution_raw)
+    file_obj.write(pauli_data_buf.getvalue())
+    pauli_data_buf.close()
+    file_obj.write(time_data)
+    file_obj.write(synth_data)
+
+
+def _read_pauli_evolution_gate(file_obj):
+    pauli_evolution_raw = struct.unpack(
+        PAULI_EVOLUTION_DEF_PACK, file_obj.read(PAULI_EVOLUTION_DEF_SIZE)
+    )
+    if pauli_evolution_raw[0] != 1 and pauli_evolution_raw[1]:
+        raise ValueError(
+            "Can't have a standalone operator with {pauli_evolution_raw[0]} operators in the payload"
+        )
+    operator_list = []
+    for _ in range(pauli_evolution_raw[0]):
+        op_size = struct.unpack(
+            SPARSE_PAULI_OP_LIST_ELEM_PACK, file_obj.read(SPARSE_PAULI_OP_LIST_ELEM_SIZE)
+        )[0]
+        operator_list.append(SparsePauliOp.from_list(np.load(io.BytesIO(file_obj.read(op_size)))))
+    if pauli_evolution_raw[1]:
+        pauli_op = operator_list[0]
+    else:
+        pauli_op = operator_list
+
+    time_type = pauli_evolution_raw[2]
+    time_data = file_obj.read(pauli_evolution_raw[3])
+    if time_type == b"f":
+        time = struct.unpack("!d", time_data)[0]
+    elif time_type == b"p":
+        with io.BytesIO(time_data) as buf:
+            time = _read_parameter(buf)
+    elif time_type == b"e":
+        with io.BytesIO(time_data) as buf:
+            time = _read_parameter_expression(buf)
+    synth_data = json.loads(file_obj.read(pauli_evolution_raw[4]))
+    synthesis = getattr(evo_synth, synth_data["class"])(**synth_data["settings"])
+    return_gate = library.PauliEvolutionGate(pauli_op, time=time, synthesis=synthesis)
+    return return_gate
+
+
 def _write_custom_instruction(file_obj, name, instruction):
-    if isinstance(instruction, Gate):
+    if isinstance(instruction, library.PauliEvolutionGate):
+        type_str = b"p"
+    elif isinstance(instruction, Gate):
         type_str = b"g"
     else:
         type_str = b"i"
@@ -946,10 +1124,13 @@ def _write_custom_instruction(file_obj, name, instruction):
     data = None
     num_qubits = instruction.num_qubits
     num_clbits = instruction.num_clbits
-    if instruction.definition:
+    if instruction.definition or type_str == b"p":
         has_definition = True
         definition_buffer = io.BytesIO()
-        _write_circuit(definition_buffer, instruction.definition)
+        if type_str == b"p":
+            _write_pauli_evolution_gate(definition_buffer, instruction)
+        else:
+            _write_circuit(definition_buffer, instruction.definition)
         definition_buffer.seek(0)
         data = definition_buffer.read()
         definition_buffer.close()
@@ -1019,7 +1200,7 @@ def dump(circuits, file_obj):
     header = struct.pack(
         FILE_HEADER_PACK,
         b"QISKIT",
-        2,
+        3,
         version_parts[0],
         version_parts[1],
         version_parts[2],
diff --git a/qiskit/opflow/evolutions/pauli_trotter_evolution.py b/qiskit/opflow/evolutions/pauli_trotter_evolution.py
index 45f7b2d8f502..5521ba3f86a3 100644
--- a/qiskit/opflow/evolutions/pauli_trotter_evolution.py
+++ b/qiskit/opflow/evolutions/pauli_trotter_evolution.py
@@ -119,7 +119,7 @@ def _recursive_convert(self, operator: OperatorBase) -> OperatorBase:
                 evo = PauliEvolutionGate(
                     pauli, time=time, synthesis=self._get_evolution_synthesis()
                 )
-                return CircuitOp(evo.definition)
+                return CircuitOp(evo)
                 # operator = EvolvedOp(operator.primitive.to_pauli_op(), coeff=operator.coeff)
             if not {"Pauli"} == operator.primitive_strings():
                 logger.warning(
diff --git a/qiskit/synthesis/evolution/evolution_synthesis.py b/qiskit/synthesis/evolution/evolution_synthesis.py
index b97be353c2ca..f904f457e49b 100644
--- a/qiskit/synthesis/evolution/evolution_synthesis.py
+++ b/qiskit/synthesis/evolution/evolution_synthesis.py
@@ -13,6 +13,7 @@
 """Evolution synthesis."""
 
 from abc import ABC, abstractmethod
+from typing import Any, Dict
 
 
 class EvolutionSynthesis(ABC):
@@ -29,3 +30,17 @@ def synthesize(self, evolution):
             QuantumCircuit: A circuit implementing the evolution.
         """
         raise NotImplementedError
+
+    @property
+    def settings(self) -> Dict[str, Any]:
+        """Return the settings in a dictionary, which can be used to reconstruct the object.
+
+        Returns:
+            A dictionary containing the settings of this product formula.
+
+        Raises:
+            NotImplementedError: The interface does not implement this method.
+        """
+        raise NotImplementedError(
+            "The settings property is not implemented for the base interface."
+        )
diff --git a/qiskit/synthesis/evolution/lie_trotter.py b/qiskit/synthesis/evolution/lie_trotter.py
index ddcd11b8f2bf..c71932dcaa57 100644
--- a/qiskit/synthesis/evolution/lie_trotter.py
+++ b/qiskit/synthesis/evolution/lie_trotter.py
@@ -12,7 +12,7 @@
 
 """The Lie-Trotter product formula."""
 
-from typing import Callable, Optional, Union
+from typing import Callable, Optional, Union, Dict, Any
 import numpy as np
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
@@ -97,3 +97,24 @@ def synthesize(self, evolution):
                 )
 
         return evolution_circuit
+
+    @property
+    def settings(self) -> Dict[str, Any]:
+        """Return the settings in a dictionary, which can be used to reconstruct the object.
+
+        Returns:
+            A dictionary containing the settings of this product formula.
+
+        Raises:
+            NotImplementedError: If a custom atomic evolution is set, which cannot be serialized.
+        """
+        if self._atomic_evolution is not None:
+            raise NotImplementedError(
+                "Cannot serialize a product formula with a custom atomic evolution."
+            )
+
+        return {
+            "reps": self.reps,
+            "insert_barriers": self.insert_barriers,
+            "cx_structure": self._cx_structure,
+        }
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index db5007f8fa2b..eda3ff938ffc 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -12,7 +12,7 @@
 
 """A product formula base for decomposing non-commuting operator exponentials."""
 
-from typing import Callable, Optional, Union
+from typing import Callable, Optional, Union, Any, Dict
 from functools import partial
 import numpy as np
 from qiskit.circuit.parameterexpression import ParameterExpression
@@ -55,11 +55,38 @@ def __init__(
         self.reps = reps
         self.insert_barriers = insert_barriers
 
+        # user-provided atomic evolution, stored for serialization
+        self._atomic_evolution = atomic_evolution
+        self._cx_structure = cx_structure
+
+        # if atomic evolution is not provided, set a default
         if atomic_evolution is None:
             atomic_evolution = partial(_default_atomic_evolution, cx_structure=cx_structure)
 
         self.atomic_evolution = atomic_evolution
 
+    @property
+    def settings(self) -> Dict[str, Any]:
+        """Return the settings in a dictionary, which can be used to reconstruct the object.
+
+        Returns:
+            A dictionary containing the settings of this product formula.
+
+        Raises:
+            NotImplementedError: If a custom atomic evolution is set, which cannot be serialized.
+        """
+        if self._atomic_evolution is not None:
+            raise NotImplementedError(
+                "Cannot serialize a product formula with a custom atomic evolution."
+            )
+
+        return {
+            "order": self.order,
+            "reps": self.reps,
+            "insert_barriers": self.insert_barriers,
+            "cx_structure": self._cx_structure,
+        }
+
 
 def evolve_pauli(
     pauli: Pauli,
diff --git a/releasenotes/notes/fix-pauli-evolution-gate-bf85592f0f8f0ba7.yaml b/releasenotes/notes/fix-pauli-evolution-gate-bf85592f0f8f0ba7.yaml
new file mode 100644
index 000000000000..79b27f4a06c0
--- /dev/null
+++ b/releasenotes/notes/fix-pauli-evolution-gate-bf85592f0f8f0ba7.yaml
@@ -0,0 +1,24 @@
+---
+fixes:
+  - |
+    Fixed the :mod:`~qiskit.circuit.qpy_serialization` support for serializing
+    a :class:`~qiskit.circuit.library.PauliEvolutionGate` object. Previously,
+    the :class:`~qiskit.circuit.library.PauliEvolutionGate` was treated as
+    a custom gate for serialization and would be deserialized as a
+    :class:`~qiskit.circuit.Gate` object that had the same definition and
+    name as the original :class:`~qiskit.circuit.library.PauliEvolutionGate`.
+    However, this would lose the original state from the
+    :class:`~qiskit.circuit.library.PauliEvolutionGate`. This has been fixed
+    so that starting in this release a
+    :class:`~qiskit.circuit.library.PauliEvolutionGate` in the circuit will
+    be preserved 1:1 across QPY serialization now. The only limitation with
+    this is that it does not support custom
+    :class:`~qiskit.synthesis.EvolutionSynthesis` classes. Only the classes
+    available from :mod:`qiskit.synthesis` can be used with a
+    :class:`~qiskit.circuit.library.PauliEvolutionGate` for qpy serialization.
+
+    To fix this issue a new QPY format version, :ref:`version_3`, was required.
+    This new format version includes a representation of the
+    :class:`~qiskit.circuit.library.PauliEvolutionGate` class which is
+    described in the :mod:`~qiskit.circuit.qpy_serialization` documentation at
+    :ref:`pauli_evo_qpy`.
diff --git a/test/python/circuit/library/test_evolved_op_ansatz.py b/test/python/circuit/library/test_evolved_op_ansatz.py
index 4d14b9c03a02..93b489c7dc82 100644
--- a/test/python/circuit/library/test_evolved_op_ansatz.py
+++ b/test/python/circuit/library/test_evolved_op_ansatz.py
@@ -37,7 +37,7 @@ def test_evolved_op_ansatz(self):
         for string, time in zip(strings, parameters):
             reference.compose(evolve(string, time), inplace=True)
 
-        self.assertEqual(evo.decompose(), reference)
+        self.assertEqual(evo.decompose().decompose(), reference)
 
     def test_custom_evolution(self):
         """Test using another evolution than the default (e.g. matrix evolution)."""
diff --git a/test/python/circuit/library/test_qaoa_ansatz.py b/test/python/circuit/library/test_qaoa_ansatz.py
index 840c691f8ded..564f329e0db0 100644
--- a/test/python/circuit/library/test_qaoa_ansatz.py
+++ b/test/python/circuit/library/test_qaoa_ansatz.py
@@ -35,7 +35,7 @@ def test_default_qaoa(self):
         circuit = circuit.decompose()
         self.assertEqual(1, len(parameters))
         self.assertIsInstance(circuit.data[0][0], HGate)
-        self.assertIsInstance(circuit.data[1][0], RXGate)
+        self.assertIsInstance(circuit.decompose().data[1][0], RXGate)
 
     def test_custom_initial_state(self):
         """Test circuit with a custom initial state."""
@@ -47,7 +47,7 @@ def test_custom_initial_state(self):
         circuit = circuit.decompose()
         self.assertEqual(1, len(parameters))
         self.assertIsInstance(circuit.data[0][0], YGate)
-        self.assertIsInstance(circuit.data[1][0], RXGate)
+        self.assertIsInstance(circuit.decompose().data[1][0], RXGate)
 
     def test_invalid_reps(self):
         """Test negative reps."""
diff --git a/test/python/circuit/test_circuit_load_from_qpy.py b/test/python/circuit/test_circuit_load_from_qpy.py
index 53dbe21e45db..9e98ee23564d 100644
--- a/test/python/circuit/test_circuit_load_from_qpy.py
+++ b/test/python/circuit/test_circuit_load_from_qpy.py
@@ -26,7 +26,7 @@
 from qiskit.circuit.library import XGate, QFT, QAOAAnsatz, PauliEvolutionGate
 from qiskit.circuit.instruction import Instruction
 from qiskit.circuit.parameter import Parameter
-from qiskit.synthesis import LieTrotter
+from qiskit.synthesis import LieTrotter, SuzukiTrotter
 from qiskit.extensions import UnitaryGate
 from qiskit.opflow import I, X, Y, Z
 from qiskit.test import QiskitTestCase
@@ -523,21 +523,70 @@ def test_single_bit_teleportation(self):
     def test_qaoa(self):
         """Test loading a QAOA circuit works."""
         cost_operator = Z ^ I ^ I ^ Z
-        qaoa = QAOAAnsatz(cost_operator)
+        qaoa = QAOAAnsatz(cost_operator, reps=2)
         qpy_file = io.BytesIO()
         dump(qaoa, qpy_file)
         qpy_file.seek(0)
         new_circ = load(qpy_file)[0]
-        # decompose to check the circuits, not their labels
-        qaoa = qaoa.decompose().decompose()
-        new_circ = new_circ.decompose().decompose()
         self.assertEqual(qaoa, new_circ)
         self.assertEqual([x[0].label for x in qaoa.data], [x[0].label for x in new_circ.data])
 
     def test_evolutiongate(self):
         """Test loading a circuit with evolution gate works."""
         synthesis = LieTrotter(reps=2)
-        evo = PauliEvolutionGate((Z ^ I) + (I ^ Z), time=0.2, synthesis=synthesis)
+        evo = PauliEvolutionGate((Z ^ I) + (I ^ Z), time=2, synthesis=synthesis)
+        qc = QuantumCircuit(2)
+        qc.append(evo, range(2))
+        qpy_file = io.BytesIO()
+        dump(qc, qpy_file)
+        qpy_file.seek(0)
+        new_circ = load(qpy_file)[0]
+
+        self.assertEqual(qc, new_circ)
+        self.assertEqual([x[0].label for x in qc.data], [x[0].label for x in new_circ.data])
+
+        new_evo = new_circ.data[0][0]
+        self.assertIsInstance(new_evo, PauliEvolutionGate)
+
+    def test_evolutiongate_param_time(self):
+        """Test loading a circuit with an evolution gate that has a parameter for time."""
+        synthesis = LieTrotter(reps=2)
+        time = Parameter("t")
+        evo = PauliEvolutionGate((Z ^ I) + (I ^ Z), time=time, synthesis=synthesis)
+        qc = QuantumCircuit(2)
+        qc.append(evo, range(2))
+        qpy_file = io.BytesIO()
+        dump(qc, qpy_file)
+        qpy_file.seek(0)
+        new_circ = load(qpy_file)[0]
+
+        self.assertEqual(qc, new_circ)
+        self.assertEqual([x[0].label for x in qc.data], [x[0].label for x in new_circ.data])
+
+        new_evo = new_circ.data[0][0]
+        self.assertIsInstance(new_evo, PauliEvolutionGate)
+
+    def test_evolutiongate_param_expr_time(self):
+        """Test loading a circuit with an evolution gate that has a parameter for time."""
+        synthesis = LieTrotter(reps=2)
+        time = Parameter("t")
+        evo = PauliEvolutionGate((Z ^ I) + (I ^ Z), time=time * time, synthesis=synthesis)
+        qc = QuantumCircuit(2)
+        qc.append(evo, range(2))
+        qpy_file = io.BytesIO()
+        dump(qc, qpy_file)
+        qpy_file.seek(0)
+        new_circ = load(qpy_file)[0]
+
+        self.assertEqual(qc, new_circ)
+        self.assertEqual([x[0].label for x in qc.data], [x[0].label for x in new_circ.data])
+
+        new_evo = new_circ.data[0][0]
+        self.assertIsInstance(new_evo, PauliEvolutionGate)
+
+    def test_op_list_evolutiongate(self):
+        """Test loading a circuit with evolution gate works."""
+        evo = PauliEvolutionGate([(Z ^ I) + (I ^ Z)] * 5, time=0.2, synthesis=None)
         qc = QuantumCircuit(2)
         qc.append(evo, range(2))
         qpy_file = io.BytesIO()
@@ -545,18 +594,28 @@ def test_evolutiongate(self):
         qpy_file.seek(0)
         new_circ = load(qpy_file)[0]
 
-        # remove wrapping of instructions
-        qc = qc.decompose().decompose()
-        new_circ = new_circ.decompose().decompose()
+        self.assertEqual(qc, new_circ)
+        self.assertEqual([x[0].label for x in qc.data], [x[0].label for x in new_circ.data])
+
+        new_evo = new_circ.data[0][0]
+        self.assertIsInstance(new_evo, PauliEvolutionGate)
+
+    def test_op_evolution_gate_suzuki_trotter(self):
+        """Test qpy path with a suzuki trotter synthesis method on an evolution gate."""
+        synthesis = SuzukiTrotter()
+        evo = PauliEvolutionGate((Z ^ I) + (I ^ Z), time=0.2, synthesis=synthesis)
+        qc = QuantumCircuit(2)
+        qc.append(evo, range(2))
+        qpy_file = io.BytesIO()
+        dump(qc, qpy_file)
+        qpy_file.seek(0)
+        new_circ = load(qpy_file)[0]
 
         self.assertEqual(qc, new_circ)
         self.assertEqual([x[0].label for x in qc.data], [x[0].label for x in new_circ.data])
 
-        # enable these tests once we can can serialize allPauliEvolutionGate parameters such as
-        # new_evo = new_circ.data[0][0]
-        # SparsePauliOp and EvolutionSynthesis
-        # self.assertIsInstance(new_evo,PauliEvolutionGate)
-        # self.assertIsInstance(new_evo.synthesis, LieTrotter)
+        new_evo = new_circ.data[0][0]
+        self.assertIsInstance(new_evo, PauliEvolutionGate)
 
     def test_parameter_expression_global_phase(self):
         """Test a circuit with a parameter expression global_phase."""
diff --git a/test/qpy_compat/test_qpy.py b/test/qpy_compat/test_qpy.py
index fba61507f302..11b58c5d2f90 100755
--- a/test/qpy_compat/test_qpy.py
+++ b/test/qpy_compat/test_qpy.py
@@ -24,7 +24,7 @@
 from qiskit.circuit.quantumregister import Qubit
 from qiskit.circuit.parameter import Parameter
 from qiskit.circuit.qpy_serialization import dump, load
-from qiskit.opflow import X, Y, Z
+from qiskit.opflow import X, Y, Z, I
 from qiskit.quantum_info.random import random_unitary
 from qiskit.circuit.library import U1Gate, U2Gate, U3Gate, QFT
 
@@ -230,7 +230,8 @@ def generate_param_phase():
     return output_circuits
 
 
-def generate_single_clbit_condition_teleportation():
+def generate_single_clbit_condition_teleportation():  # pylint: disable=invalid-name
+    """Generate single clbit condition teleportation circuit."""
     qr = QuantumRegister(1)
     cr = ClassicalRegister(2, name="name")
     teleport_qc = QuantumCircuit(qr, cr, name="Reset Test")
@@ -241,6 +242,19 @@ def generate_single_clbit_condition_teleportation():
     return teleport_qc
 
 
+def generate_evolution_gate():
+    """Generate a circuit with a pauli evolution gate."""
+    # Runtime import since this only exists in terra 0.19.0
+    from qiskit.circuit.library import PauliEvolutionGate
+    from qiskit.synthesis import SuzukiTrotter
+
+    synthesis = SuzukiTrotter()
+    evo = PauliEvolutionGate([(Z ^ I) + (I ^ Z)] * 5, time=2.0, synthesis=synthesis)
+    qc = QuantumCircuit(2)
+    qc.append(evo, range(2))
+    return qc
+
+
 def generate_circuits(version_str=None):
     """Generate reference circuits."""
     version_parts = None
@@ -265,6 +279,9 @@ def generate_circuits(version_str=None):
     if version_parts >= (0, 19, 0):
         output_circuits["param_phase.qpy"] = generate_param_phase()
 
+    if version_parts >= (0, 19, 1):
+        output_circuits["pauli_evo.qpy"] = [generate_evolution_gate()]
+
     return output_circuits