From 3a155896966cab361d9a5bfb311dd458f9253eca Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Thu, 3 Oct 2024 16:11:27 +0200
Subject: [PATCH 01/25] py version for expand

---
 .../operators/symplectic/sparse_pauli_op.py   | 16 +++++
 qiskit/synthesis/evolution/lie_trotter.py     | 27 ++++----
 qiskit/synthesis/evolution/product_formula.py | 49 ++++++++++++++-
 qiskit/synthesis/evolution/qdrift.py          | 51 +++++++--------
 qiskit/synthesis/evolution/suzuki_trotter.py  | 62 +++++++++++--------
 .../circuit/library/test_evolution_gate.py    |  1 +
 6 files changed, 132 insertions(+), 74 deletions(-)

diff --git a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index 5699944788be..18f3af600e80 100644
--- a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
+++ b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
@@ -935,6 +935,14 @@ def to_list(self, array: bool = False):
             return labels
         return labels.tolist()
 
+    def to_sparse_list(self):
+        """Convert to a sparse Pauli list format with elements (pauli, qubits, coefficient)."""
+        pauli_labels = self.paulis.to_labels()
+        sparse_list = [
+            (*sparsify_label(label), coeff) for label, coeff in zip(pauli_labels, self.coeffs)
+        ]
+        return sparse_list
+
     def to_matrix(self, sparse: bool = False, force_serial: bool = False) -> np.ndarray:
         """Convert to a dense or sparse matrix.
 
@@ -1176,5 +1184,13 @@ def apply_layout(
         return new_op.compose(self, qargs=layout)
 
 
+def sparsify_label(pauli_string):
+    """Return a sparse format of a Pauli string, e.g. "XIIIZ" -> ("XZ", [0, 4])."""
+    qubits = [i for i, label in enumerate(reversed(pauli_string)) if label != "I"]
+    sparse_label = "".join(pauli_string[~i] for i in qubits)
+
+    return sparse_label, qubits
+
+
 # Update docstrings for API docs
 generate_apidocs(SparsePauliOp)
diff --git a/qiskit/synthesis/evolution/lie_trotter.py b/qiskit/synthesis/evolution/lie_trotter.py
index 1a01675a6782..954102f85d06 100644
--- a/qiskit/synthesis/evolution/lie_trotter.py
+++ b/qiskit/synthesis/evolution/lie_trotter.py
@@ -16,8 +16,8 @@
 
 import inspect
 from collections.abc import Callable
+from itertools import chain
 from typing import Any
-import numpy as np
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
@@ -100,25 +100,22 @@ def __init__(
         """
         super().__init__(1, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
-    def synthesize(self, evolution):
-        # get operators and time to evolve
-        operators = evolution.operator
-        time = evolution.time
+    def expand(self, evolution):
+        operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
 
         # construct the evolution circuit
-        single_rep = QuantumCircuit(operators[0].num_qubits)
 
-        if not isinstance(operators, list):
-            pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operators.to_list()]
+        if isinstance(operators, list):
+            # the expansion formula is the same for commuting and non-commuting bits, so
+            # we'll just concatenate all Paulis
+            non_commuting = list(chain.from_iterable([op.to_sparse_list() for op in operators]))
         else:
-            pauli_list = [(op, 1) for op in operators]
+            # Assume no commutativity here. If we were to group commuting Paulis,
+            # here would be the location to do so.
+            non_commuting = operators.to_sparse_list()
 
-        for i, (op, coeff) in enumerate(pauli_list):
-            self.atomic_evolution(single_rep, op, coeff * time / self.reps)
-            if self.insert_barriers and i != len(pauli_list) - 1:
-                single_rep.barrier()
-
-        return single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
+        # we're already done here since Lie Trotter does not do any operator repetition
+        return non_commuting
 
     @property
     def settings(self) -> dict[str, Any]:
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index df38a2a541ab..cb0f75c9334a 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -6,7 +6,7 @@
 # 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
+# typing.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.
 
@@ -16,8 +16,8 @@
 
 import inspect
 from collections.abc import Callable
-from typing import Any
 from functools import partial
+import typing
 import numpy as np
 from qiskit.circuit.parameterexpression import ParameterExpression
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -26,6 +26,9 @@
 
 from .evolution_synthesis import EvolutionSynthesis
 
+if typing.TYPE_CHECKING:
+    from qiskit.circuit.library import PauliEvolutionGate
+
 
 class ProductFormula(EvolutionSynthesis):
     """Product formula base class for the decomposition of non-commuting operator exponentials.
@@ -108,8 +111,48 @@ def wrap_atomic_evolution(output, operator, time):
         else:
             self.atomic_evolution = atomic_evolution
 
+    def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], float]]:
+        """Apply the product formula to expand the Hamiltonian in the evolution gate.
+
+        Args:
+            evolution: The :class:`.PauliEvolutionGate`, whose Hamiltonian we expand.
+
+        Returns:
+            A list of Pauli rotations in a sparse format, where each element is
+            ``(paulistring, qubits, coefficient)``. For example, the Lie-Trotter expansion
+            of ``H = XI + ZZ`` would return ``[("X", [0], 1), ("ZZ", [0, 1], 1)]``.
+        """
+        raise NotImplementedError(
+            f"The method ``expand`` is not implemented for {self.__class__}. Implement it to "
+            f"automatically enable the call to {self.__class__}.synthesize."
+        )
+
+    def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
+        """Synthesize a :class:`.PauliEvolutionGate`.
+
+        Args:
+            evolution: The evolution gate to synthesize.
+
+        Returns:
+            QuantumCircuit: A circuit implementing the evolution.
+        """
+        pauli_rotations = self.expand(evolution)
+
+        time = evolution.time
+        num_qubits = evolution.num_qubits
+
+        single_rep = QuantumCircuit(num_qubits)
+        for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
+            op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
+            angle = np.real_if_close(coeff * time / self.reps)
+            self.atomic_evolution(single_rep, op, angle)
+            if self.insert_barriers and i != len(pauli_rotations) - 1:
+                single_rep.barrier()
+
+        return single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
+
     @property
-    def settings(self) -> dict[str, Any]:
+    def settings(self) -> dict[str, typing.Any]:
         """Return the settings in a dictionary, which can be used to reconstruct the object.
 
         Returns:
diff --git a/qiskit/synthesis/evolution/qdrift.py b/qiskit/synthesis/evolution/qdrift.py
index 7c68f66dc8ea..653516dbd6c1 100644
--- a/qiskit/synthesis/evolution/qdrift.py
+++ b/qiskit/synthesis/evolution/qdrift.py
@@ -16,6 +16,8 @@
 
 import inspect
 import math
+import typing
+from itertools import chain
 from collections.abc import Callable
 import numpy as np
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -23,7 +25,9 @@
 from qiskit.utils.deprecation import deprecate_arg
 
 from .product_formula import ProductFormula
-from .lie_trotter import LieTrotter
+
+if typing.TYPE_CHECKING:
+    from qiskit.circuit.library import PauliEvolutionGate
 
 
 class QDrift(ProductFormula):
@@ -88,44 +92,33 @@ def __init__(
         self.sampled_ops = None
         self.rng = np.random.default_rng(seed)
 
-    def synthesize(self, evolution):
-        # get operators and time to evolve
+    def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], float]]:
         operators = evolution.operator
-        time = evolution.time
+        time = evolution.time  # used to dtermine the number of gates
 
-        if not isinstance(operators, list):
-            pauli_list = [(Pauli(op), coeff) for op, coeff in operators.to_list()]
-            coeffs = [np.real(coeff) for op, coeff in operators.to_list()]
+        # QDrift is based on first-order Lie-Trotter, hence we can just concatenate all
+        # Pauli terms and ignore commutations
+        if isinstance(operators, list):
+            paulis = list(chain.from_iterable([op.to_sparse_list() for op in operators]))
         else:
-            pauli_list = [(op, 1) for op in operators]
-            coeffs = [1 for op in operators]
+            paulis = operators.to_sparse_list()
+
+        coeffs = [np.real(coeff) for _, _, coeff in paulis]
 
         # We artificially make the weights positive
         weights = np.abs(coeffs)
         lambd = np.sum(weights)
 
         num_gates = math.ceil(2 * (lambd**2) * (time**2) * self.reps)
+
         # The protocol calls for the removal of the individual coefficients,
         # and multiplication by a constant evolution time.
-        evolution_time = lambd * time / num_gates
-
-        self.sampled_ops = self.rng.choice(
-            np.array(pauli_list, dtype=object),
-            size=(num_gates,),
-            p=weights / lambd,
-        )
-
-        # pylint: disable=cyclic-import
-        from qiskit.circuit.library.pauli_evolution import PauliEvolutionGate
-
-        # Build the evolution circuit using the LieTrotter synthesis with the sampled operators
-        lie_trotter = LieTrotter(
-            insert_barriers=self.insert_barriers, atomic_evolution=self.atomic_evolution
+        sampled = self.rng.choice(
+            np.array(paulis, dtype=object), size=(num_gates,), p=weights / lambd
         )
-        evolution_circuit = PauliEvolutionGate(
-            sum(SparsePauliOp(np.sign(coeff) * op) for op, coeff in self.sampled_ops),
-            time=evolution_time,
-            synthesis=lie_trotter,
-        ).definition
 
-        return evolution_circuit
+        rescaled_time = lambd / num_gates
+        sampled_paulis = [
+            (pauli[0], pauli[1], np.sign(pauli[2]) * rescaled_time) for pauli in sampled
+        ]
+        return sampled_paulis
diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index e03fd27e26d4..0ae0b1da9636 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -16,8 +16,7 @@
 
 import inspect
 from collections.abc import Callable
-
-import numpy as np
+from itertools import chain
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
@@ -105,51 +104,60 @@ def __init__(
             ValueError: If order is not even
         """
 
-        if order % 2 == 1:
+        if order > 1 and order % 2 == 1:
             raise ValueError(
                 "Suzuki product formulae are symmetric and therefore only defined "
-                "for even orders."
+                "for even orders (or order==1)."
             )
         super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
-    def synthesize(self, evolution):
-        # get operators and time to evolve
-        operators = evolution.operator
-        time = evolution.time
-
-        if not isinstance(operators, list):
-            pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operators.to_list()]
-        else:
-            pauli_list = [(op, 1) for op in operators]
-
-        ops_to_evolve = self._recurse(self.order, time / self.reps, pauli_list)
+    def expand(self, evolution):
+        operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
 
         # construct the evolution circuit
-        single_rep = QuantumCircuit(operators[0].num_qubits)
 
-        for i, (op, coeff) in enumerate(ops_to_evolve):
-            self.atomic_evolution(single_rep, op, coeff)
-            if self.insert_barriers and i != len(ops_to_evolve) - 1:
-                single_rep.barrier()
+        if isinstance(operators, list):  # already sorted into commuting bits
+            non_commuting = [operator.to_sparse_list() for operator in operators]
+        else:
+            # Assume no commutativity here. If we were to group commuting Paulis,
+            # here would be the location to do so.
+            non_commuting = [[op] for op in operators.to_sparse_list()]
 
-        return single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
+        # we're already done here since Lie Trotter does not do any operator repetition
+        product_formula = self._recurse(self.order, non_commuting)
+        flattened = list(chain.from_iterable(product_formula))
+        return flattened
 
     @staticmethod
-    def _recurse(order, time, pauli_list):
+    def _recurse(order, grouped_paulis):
         if order == 1:
-            return pauli_list
+            return grouped_paulis
 
         elif order == 2:
-            halves = [(op, coeff * time / 2) for op, coeff in pauli_list[:-1]]
-            full = [(pauli_list[-1][0], time * pauli_list[-1][1])]
+            halves = [
+                [(label, qubits, coeff / 2) for label, qubits, coeff in paulis]
+                for paulis in grouped_paulis[:-1]
+            ]
+            full = [grouped_paulis[-1]]
             return halves + full + list(reversed(halves))
 
         else:
             reduction = 1 / (4 - 4 ** (1 / (order - 1)))
             outer = 2 * SuzukiTrotter._recurse(
-                order - 2, time=reduction * time, pauli_list=pauli_list
+                order - 2,
+                [
+                    [(label, qubits, coeff * reduction) for label, qubits, coeff in paulis]
+                    for paulis in grouped_paulis
+                ],
             )
             inner = SuzukiTrotter._recurse(
-                order - 2, time=(1 - 4 * reduction) * time, pauli_list=pauli_list
+                order - 2,
+                [
+                    [
+                        (label, qubits, coeff * (1 - 4 * reduction))
+                        for label, qubits, coeff in paulis
+                    ]
+                    for paulis in grouped_paulis
+                ],
             )
             return outer + inner + outer
diff --git a/test/python/circuit/library/test_evolution_gate.py b/test/python/circuit/library/test_evolution_gate.py
index a7c9a73abb5f..94825729a70b 100644
--- a/test/python/circuit/library/test_evolution_gate.py
+++ b/test/python/circuit/library/test_evolution_gate.py
@@ -185,6 +185,7 @@ def test_passing_grouped_paulis(self, wrap):
             decomposed = evo_gate.definition.decompose()
         else:
             decomposed = evo_gate.definition
+
         self.assertEqual(decomposed.count_ops()["rz"], 4)
         self.assertEqual(decomposed.count_ops()["rzz"], 1)
         self.assertEqual(decomposed.count_ops()["rxx"], 1)

From 72f26897707438af03bd097a7f2968f526cda248 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 7 Oct 2024 17:43:11 +0200
Subject: [PATCH 02/25] expand fully & simplify lie trotter

---
 qiskit/synthesis/evolution/lie_trotter.py     | 36 ++-----------------
 qiskit/synthesis/evolution/product_formula.py | 14 ++++----
 qiskit/synthesis/evolution/qdrift.py          |  4 +--
 qiskit/synthesis/evolution/suzuki_trotter.py  | 15 +++++---
 4 files changed, 23 insertions(+), 46 deletions(-)

diff --git a/qiskit/synthesis/evolution/lie_trotter.py b/qiskit/synthesis/evolution/lie_trotter.py
index 954102f85d06..89ba03590326 100644
--- a/qiskit/synthesis/evolution/lie_trotter.py
+++ b/qiskit/synthesis/evolution/lie_trotter.py
@@ -22,10 +22,10 @@
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
 
-from .product_formula import ProductFormula
+from .suzuki_trotter import SuzukiTrotter
 
 
-class LieTrotter(ProductFormula):
+class LieTrotter(SuzukiTrotter):
     r"""The Lie-Trotter product formula.
 
     The Lie-Trotter formula approximates the exponential of two non-commuting operators
@@ -52,21 +52,6 @@ class LieTrotter(ProductFormula):
         `arXiv:math-ph/0506007 <https://arxiv.org/pdf/math-ph/0506007.pdf>`_
     """
 
-    @deprecate_arg(
-        name="atomic_evolution",
-        since="1.2",
-        predicate=lambda callable: callable is not None
-        and len(inspect.signature(callable).parameters) == 2,
-        deprecation_description=(
-            "The 'Callable[[Pauli | SparsePauliOp, float], QuantumCircuit]' signature of the "
-            "'atomic_evolution' argument"
-        ),
-        additional_msg=(
-            "Instead you should update your 'atomic_evolution' function to be of the following "
-            "type: 'Callable[[QuantumCircuit, Pauli | SparsePauliOp, float], None]'."
-        ),
-        pending=True,
-    )
     def __init__(
         self,
         reps: int = 1,
@@ -100,23 +85,6 @@ def __init__(
         """
         super().__init__(1, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
-    def expand(self, evolution):
-        operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
-
-        # construct the evolution circuit
-
-        if isinstance(operators, list):
-            # the expansion formula is the same for commuting and non-commuting bits, so
-            # we'll just concatenate all Paulis
-            non_commuting = list(chain.from_iterable([op.to_sparse_list() for op in operators]))
-        else:
-            # Assume no commutativity here. If we were to group commuting Paulis,
-            # here would be the location to do so.
-            non_commuting = operators.to_sparse_list()
-
-        # we're already done here since Lie Trotter does not do any operator repetition
-        return non_commuting
-
     @property
     def settings(self) -> dict[str, Any]:
         """Return the settings in a dictionary, which can be used to reconstruct the object.
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index cb0f75c9334a..46e7d8067854 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -120,7 +120,7 @@ def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], f
         Returns:
             A list of Pauli rotations in a sparse format, where each element is
             ``(paulistring, qubits, coefficient)``. For example, the Lie-Trotter expansion
-            of ``H = XI + ZZ`` would return ``[("X", [0], 1), ("ZZ", [0, 1], 1)]``.
+            of ``H = XI + ZZ`` would return ``[("X", [1], 1), ("ZZ", [0, 1], 1)]``.
         """
         raise NotImplementedError(
             f"The method ``expand`` is not implemented for {self.__class__}. Implement it to "
@@ -141,15 +141,17 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
         time = evolution.time
         num_qubits = evolution.num_qubits
 
-        single_rep = QuantumCircuit(num_qubits)
+        circuit = QuantumCircuit(num_qubits)
+
+        # we could cache the circuit decomposition for the circuits we already saw
         for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
             op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
-            angle = np.real_if_close(coeff * time / self.reps)
-            self.atomic_evolution(single_rep, op, angle)
+
+            self.atomic_evolution(circuit, op, np.real_if_close(coeff))
             if self.insert_barriers and i != len(pauli_rotations) - 1:
-                single_rep.barrier()
+                circuit.barrier()
 
-        return single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
+        return circuit
 
     @property
     def settings(self) -> dict[str, typing.Any]:
diff --git a/qiskit/synthesis/evolution/qdrift.py b/qiskit/synthesis/evolution/qdrift.py
index 653516dbd6c1..d3882066204f 100644
--- a/qiskit/synthesis/evolution/qdrift.py
+++ b/qiskit/synthesis/evolution/qdrift.py
@@ -94,7 +94,7 @@ def __init__(
 
     def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], float]]:
         operators = evolution.operator
-        time = evolution.time  # used to dtermine the number of gates
+        time = evolution.time  # used to determine the number of gates
 
         # QDrift is based on first-order Lie-Trotter, hence we can just concatenate all
         # Pauli terms and ignore commutations
@@ -117,7 +117,7 @@ def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], f
             np.array(paulis, dtype=object), size=(num_gates,), p=weights / lambd
         )
 
-        rescaled_time = lambd / num_gates
+        rescaled_time = lambd * time / num_gates
         sampled_paulis = [
             (pauli[0], pauli[1], np.sign(pauli[2]) * rescaled_time) for pauli in sampled
         ]
diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index 0ae0b1da9636..0a18eefd3e67 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -112,20 +112,27 @@ def __init__(
         super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
     def expand(self, evolution):
+        """
+        H = ZZ + IX --> ("X", [0], 1/2), ("ZZ", [0, 1], 1), ("X", [0], 1/2)
+
+        ("X", [0], 1/2), ("ZZ", [0, 1], 1), ("X", [0], 1), ("ZZ", [0, 1], 1), ("X", [0], 1/2)
+        """
         operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
+        time = evolution.time
 
         # construct the evolution circuit
-
         if isinstance(operators, list):  # already sorted into commuting bits
-            non_commuting = [operator.to_sparse_list() for operator in operators]
+            non_commuting = [
+                (time / self.reps * operator).to_sparse_list() for operator in operators
+            ]
         else:
             # Assume no commutativity here. If we were to group commuting Paulis,
             # here would be the location to do so.
-            non_commuting = [[op] for op in operators.to_sparse_list()]
+            non_commuting = [[op] for op in (time / self.reps * operators).to_sparse_list()]
 
         # we're already done here since Lie Trotter does not do any operator repetition
         product_formula = self._recurse(self.order, non_commuting)
-        flattened = list(chain.from_iterable(product_formula))
+        flattened = self.reps * list(chain.from_iterable(product_formula))
         return flattened
 
     @staticmethod

From 8cd2c25d7eb21952897d5c5909149c92cd673e7e Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 8 Oct 2024 08:28:49 +0200
Subject: [PATCH 03/25] use examples that actually do not commute

---
 qiskit/synthesis/evolution/lie_trotter.py    | 5 +----
 qiskit/synthesis/evolution/suzuki_trotter.py | 2 +-
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/qiskit/synthesis/evolution/lie_trotter.py b/qiskit/synthesis/evolution/lie_trotter.py
index 89ba03590326..5b4e91721a6e 100644
--- a/qiskit/synthesis/evolution/lie_trotter.py
+++ b/qiskit/synthesis/evolution/lie_trotter.py
@@ -14,13 +14,10 @@
 
 from __future__ import annotations
 
-import inspect
 from collections.abc import Callable
-from itertools import chain
 from typing import Any
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
-from qiskit.utils.deprecation import deprecate_arg
 
 from .suzuki_trotter import SuzukiTrotter
 
@@ -40,7 +37,7 @@ class LieTrotter(SuzukiTrotter):
 
     .. math::
 
-        e^{-it(XX + ZZ)} = e^{-it XX}e^{-it ZZ} + \mathcal{O}(t^2).
+        e^{-it(XI + ZZ)} = e^{-it XI}e^{-it ZZ} + \mathcal{O}(t^2).
 
     References:
 
diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index 0a18eefd3e67..264591430f0b 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -43,7 +43,7 @@ class SuzukiTrotter(ProductFormula):
 
     .. math::
 
-        e^{-it(XX + ZZ)} = e^{-it/2 ZZ}e^{-it XX}e^{-it/2 ZZ} + \mathcal{O}(t^3).
+        e^{-it(XI + ZZ)} = e^{-it/2 XI}e^{-it ZZ}e^{-it/2 XI} + \mathcal{O}(t^3).
 
     References:
         [1]: D. Berry, G. Ahokas, R. Cleve and B. Sanders,

From 681105d634063843af120c8a9a11245a0bf1a91b Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 8 Oct 2024 14:07:30 +0200
Subject: [PATCH 04/25] add plugin structure

---
 pyproject.toml                                    |  1 +
 qiskit/synthesis/evolution/product_formula.py     |  1 -
 qiskit/transpiler/passes/synthesis/hls_plugins.py | 12 +++++++++++-
 3 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 70a6e7c9d0cf..939c17acb9dd 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -100,6 +100,7 @@ sk = "qiskit.transpiler.passes.synthesis.solovay_kitaev_synthesis:SolovayKitaevS
 "qft.line" = "qiskit.transpiler.passes.synthesis.hls_plugins:QFTSynthesisLine"
 "qft.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:QFTSynthesisFull"
 "permutation.token_swapper" = "qiskit.transpiler.passes.synthesis.hls_plugins:TokenSwapperSynthesisPermutation"
+"evolution.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:PauliEvolutionSynthesisDefault"
 
 [project.entry-points."qiskit.transpiler.init"]
 default = "qiskit.transpiler.preset_passmanagers.builtin_plugins:DefaultInitPassManager"
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index 46e7d8067854..839940151aa1 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -138,7 +138,6 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
         """
         pauli_rotations = self.expand(evolution)
 
-        time = evolution.time
         num_qubits = evolution.num_qubits
 
         circuit = QuantumCircuit(num_qubits)
diff --git a/qiskit/transpiler/passes/synthesis/hls_plugins.py b/qiskit/transpiler/passes/synthesis/hls_plugins.py
index fa827a20ed48..3f563274f871 100644
--- a/qiskit/transpiler/passes/synthesis/hls_plugins.py
+++ b/qiskit/transpiler/passes/synthesis/hls_plugins.py
@@ -251,6 +251,8 @@
 
 """
 
+from __future__ import annotations
+
 import numpy as np
 import rustworkx as rx
 
@@ -288,8 +290,8 @@
     synth_mcx_1_clean_b95,
     synth_mcx_gray_code,
     synth_mcx_noaux_v24,
+    synth_mcmt_vchain,
 )
-from qiskit.synthesis.multi_controlled import synth_mcmt_vchain
 from qiskit.transpiler.passes.routing.algorithms import ApproximateTokenSwapper
 from .plugin import HighLevelSynthesisPlugin
 
@@ -1029,3 +1031,11 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **
             high_level_object.num_target_qubits,
             ctrl_state,
         )
+
+
+class PauliEvolutionSynthesisDefault(HighLevelSynthesisPlugin):
+    """The default implementation calling the attached synthesis algorithm."""
+
+    def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **options):
+        algo = high_level_object.synthesis
+        return algo.synthesize(high_level_object)

From 58e4a03c5c99d1667c28fba0c6088f514492219a Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Fri, 11 Oct 2024 10:13:18 +0200
Subject: [PATCH 05/25] fix plugin name, rm complex from expand

--> expand should return float | ParameterExpression
---
 pyproject.toml                                |  2 +-
 qiskit/synthesis/evolution/product_formula.py | 29 +++++++++++++++++--
 qiskit/synthesis/evolution/qdrift.py          |  8 +++--
 qiskit/synthesis/evolution/suzuki_trotter.py  |  4 ++-
 4 files changed, 36 insertions(+), 7 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 939c17acb9dd..72f3b4c14d54 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -100,7 +100,7 @@ sk = "qiskit.transpiler.passes.synthesis.solovay_kitaev_synthesis:SolovayKitaevS
 "qft.line" = "qiskit.transpiler.passes.synthesis.hls_plugins:QFTSynthesisLine"
 "qft.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:QFTSynthesisFull"
 "permutation.token_swapper" = "qiskit.transpiler.passes.synthesis.hls_plugins:TokenSwapperSynthesisPermutation"
-"evolution.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:PauliEvolutionSynthesisDefault"
+"PauliEvolution.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:PauliEvolutionSynthesisDefault"
 
 [project.entry-points."qiskit.transpiler.init"]
 default = "qiskit.transpiler.preset_passmanagers.builtin_plugins:DefaultInitPassManager"
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index 839940151aa1..484a29295d81 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -20,7 +20,7 @@
 import typing
 import numpy as np
 from qiskit.circuit.parameterexpression import ParameterExpression
-from qiskit.circuit.quantumcircuit import QuantumCircuit
+from qiskit.circuit.quantumcircuit import QuantumCircuit, ParameterValueType
 from qiskit.quantum_info import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
 
@@ -111,7 +111,9 @@ def wrap_atomic_evolution(output, operator, time):
         else:
             self.atomic_evolution = atomic_evolution
 
-    def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], float]]:
+    def expand(
+        self, evolution: PauliEvolutionGate
+    ) -> list[tuple[str, tuple[int], ParameterValueType]]:
         """Apply the product formula to expand the Hamiltonian in the evolution gate.
 
         Args:
@@ -146,7 +148,7 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
         for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
             op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
 
-            self.atomic_evolution(circuit, op, np.real_if_close(coeff))
+            self.atomic_evolution(circuit, op, coeff)
             if self.insert_barriers and i != len(pauli_rotations) - 1:
                 circuit.barrier()
 
@@ -175,6 +177,12 @@ def settings(self) -> dict[str, typing.Any]:
             "wrap": self._wrap,
         }
 
+    def _normalize_coefficients(
+        self, paulis: list[str | list[int], float | complex | ParameterExpression]
+    ) -> list[str | list[int] | ParameterValueType]:
+        """Ensure the coefficients are real (or parameter expressions)."""
+        return [[(op, qubits, real_or_fail(coeff)) for op, qubits, coeff in ops] for ops in paulis]
+
 
 def evolve_pauli(
     output: QuantumCircuit,
@@ -426,3 +434,18 @@ def _default_atomic_evolution(output, operator, time, cx_structure, wrap):
         pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operator.to_list()]
         for pauli, coeff in pauli_list:
             evolve_pauli(output, pauli, coeff * time, cx_structure, wrap)
+
+
+def real_or_fail(value, tol=100):
+    """Return real if close, otherwise fail. Unbound parameters are left unchanged.
+
+    Based on NumPy's ``real_if_close``, i.e. ``tol`` is in terms of machine precision for float.
+    """
+    if isinstance(value, ParameterExpression):
+        return value
+
+    abstol = tol * np.finfo(float).eps
+    if abs(np.imag(value)) < abstol:
+        return np.real(value)
+
+    raise ValueError(f"Encountered complex value {value}, but expected real.")
diff --git a/qiskit/synthesis/evolution/qdrift.py b/qiskit/synthesis/evolution/qdrift.py
index d3882066204f..e214d4c113dc 100644
--- a/qiskit/synthesis/evolution/qdrift.py
+++ b/qiskit/synthesis/evolution/qdrift.py
@@ -23,6 +23,7 @@
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
+from qiskit.exceptions import QiskitError
 
 from .product_formula import ProductFormula
 
@@ -103,7 +104,10 @@ def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], f
         else:
             paulis = operators.to_sparse_list()
 
-        coeffs = [np.real(coeff) for _, _, coeff in paulis]
+        try:
+            coeffs = [float(np.real_if_close(coeff)) for _, _, coeff in paulis]
+        except TypeError as exc:
+            raise QiskitError("QDrift requires bound, real coefficients.") from exc
 
         # We artificially make the weights positive
         weights = np.abs(coeffs)
@@ -119,6 +123,6 @@ def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], f
 
         rescaled_time = lambd * time / num_gates
         sampled_paulis = [
-            (pauli[0], pauli[1], np.sign(pauli[2]) * rescaled_time) for pauli in sampled
+            (pauli[0], pauli[1], np.real(np.sign(pauli[2])) * rescaled_time) for pauli in sampled
         ]
         return sampled_paulis
diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index 264591430f0b..44f66f3d5f3b 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -22,7 +22,6 @@
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
 
-
 from .product_formula import ProductFormula
 
 
@@ -130,6 +129,9 @@ def expand(self, evolution):
             # here would be the location to do so.
             non_commuting = [[op] for op in (time / self.reps * operators).to_sparse_list()]
 
+        # normalize coefficients, i.e. ensure they are float or ParameterExpression
+        non_commuting = self._normalize_coefficients(non_commuting)
+
         # we're already done here since Lie Trotter does not do any operator repetition
         product_formula = self._recurse(self.order, non_commuting)
         flattened = self.reps * list(chain.from_iterable(product_formula))

From 7f2ec3dbfeec2d70aee9a38f2b75ff1883d445e6 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Fri, 11 Oct 2024 16:11:53 +0200
Subject: [PATCH 06/25] paulievo in Rust

---
 crates/accelerate/src/circuit_library/mod.rs  |   1 +
 .../src/circuit_library/pauli_evolution.rs    | 209 ++++++++++++++++++
 .../src/circuit_library/pauli_feature_map.rs  | 109 +--------
 3 files changed, 221 insertions(+), 98 deletions(-)
 create mode 100644 crates/accelerate/src/circuit_library/pauli_evolution.rs

diff --git a/crates/accelerate/src/circuit_library/mod.rs b/crates/accelerate/src/circuit_library/mod.rs
index bc960bab2aab..09692d049f12 100644
--- a/crates/accelerate/src/circuit_library/mod.rs
+++ b/crates/accelerate/src/circuit_library/mod.rs
@@ -13,6 +13,7 @@
 use pyo3::prelude::*;
 
 mod entanglement;
+mod pauli_evolution;
 mod pauli_feature_map;
 mod quantum_volume;
 
diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
new file mode 100644
index 000000000000..7a8db8380258
--- /dev/null
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -0,0 +1,209 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2024
+//
+// 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.
+
+use itertools::Itertools;
+use pyo3::prelude::*;
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::Qubit;
+use smallvec::{smallvec, SmallVec};
+use std::f64::consts::PI;
+
+// custom math and types for a more readable code
+const PI2: f64 = PI / 2.;
+type StandardInstruction = (StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>);
+
+/// Return instructions (using only StandardGate operations) to implement a Pauli evolution
+/// of a given Pauli string over a given time (as Param).
+///
+/// The Pauli evolution is implemented as a basis transformation to the Pauli-Z basis,
+/// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
+/// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
+/// Pauli string XIYZ we have the circuit
+///                     ┌───┐┌───────┐┌───┐
+/// 0: ─────────────────┤ X ├┤ Rz(2) ├┤ X ├──────────────────
+///    ┌──────────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌───────────┐
+/// 1: ┤ Rx(pi/2) ├┤ X ├──■─────────────■──┤ X ├┤ Rx(-pi/2) ├
+///    └──────────┘└─┬─┘                   └─┬─┘└───────────┘
+/// 2: ──────────────┼───────────────────────┼───────────────
+///     ┌───┐        │                       │  ┌───┐
+/// 3: ─┤ H ├────────■───────────────────────■──┤ H ├────────
+///     └───┘                                   └───┘
+///
+/// If ``phase_gate`` is ``true``, use the ``PhaseGate`` instead of ``RZGate`` as single-qubit
+/// rotation.
+pub fn pauli_evolution<'a>(
+    pauli: &'a str,
+    indices: Vec<u32>,
+    time: Param,
+    phase_gate: bool,
+    // check_sparse: bool
+) -> Box<dyn Iterator<Item = StandardInstruction> + 'a> {
+    // ensure the Pauli has no identity terms
+    let binding = pauli.to_lowercase(); // lowercase for convenience
+    let active = binding
+        .as_str()
+        .chars()
+        .rev()
+        .zip(indices)
+        .filter(|(pauli, _)| *pauli != 'i');
+    let (paulis, indices): (Vec<char>, Vec<u32>) = active.unzip();
+
+    match (phase_gate, indices.len()) {
+        (false, 0) => Box::new(std::iter::empty()),
+        (false, 1) => Box::new(single_qubit_evolution(paulis[0], indices[0], time)),
+        (false, 2) => two_qubit_evolution(paulis.clone(), indices.clone(), time),
+        _ => Box::new(multi_qubit_evolution(
+            paulis.clone(),
+            indices.clone(),
+            time,
+            phase_gate,
+        )),
+    }
+}
+
+fn single_qubit_evolution(
+    pauli: char,
+    index: u32,
+    time: Param,
+) -> impl Iterator<Item = StandardInstruction> {
+    let qubit: SmallVec<[Qubit; 2]> = smallvec![Qubit(index)];
+    let param: SmallVec<[Param; 3]> = smallvec![time.clone()];
+
+    std::iter::once(match pauli {
+        'x' => (StandardGate::RXGate, param, qubit),
+        'y' => (StandardGate::RXGate, param, qubit),
+        'z' => (StandardGate::RXGate, param, qubit),
+        _ => unreachable!("Unsupported Pauli, at this point we expected one of x, y, z."),
+    })
+}
+
+fn two_qubit_evolution<'a>(
+    pauli: Vec<char>,
+    indices: Vec<u32>,
+    time: Param,
+) -> Box<dyn Iterator<Item = StandardInstruction> + 'a> {
+    let qubits: SmallVec<[Qubit; 2]> = smallvec![Qubit(indices[0]), Qubit(indices[1])];
+    let param: SmallVec<[Param; 3]> = smallvec![time.clone()];
+    let paulistring: String = pauli.iter().collect();
+
+    match paulistring.as_str() {
+        "xx" => Box::new(std::iter::once((StandardGate::RXXGate, param, qubits))),
+        "zx" => Box::new(std::iter::once((StandardGate::RZXGate, param, qubits))),
+        "yy" => Box::new(std::iter::once((StandardGate::RYYGate, param, qubits))),
+        "zz" => Box::new(std::iter::once((StandardGate::RZZGate, param, qubits))),
+        _ => Box::new(multi_qubit_evolution(pauli, indices, time, false)),
+    }
+}
+
+fn multi_qubit_evolution(
+    pauli: Vec<char>,
+    indices: Vec<u32>,
+    time: Param,
+    phase_gate: bool,
+) -> impl Iterator<Item = StandardInstruction> {
+    // Get pairs of (pauli, qubit) that are active, i.e. that are not the identity. Note that
+    // the rest of the code also works if there are only identities, in which case we will
+    // effectively return an empty iterator.
+    // let qubits: Vec<Qubit> = indices.iter().map(|i| Qubit(*i)).collect();
+    // let binding = pauli.to_lowercase(); // lowercase for convenience
+    // let active_paulis = binding
+    //     .as_str()
+    //     .chars()
+    //     .rev() // reverse due to Qiskit's bit ordering convention
+    //     .zip(qubits)
+    //     .filter(|(p, _)| *p != 'i')
+    //     .collect_vec();
+    let active_paulis: Vec<(char, Qubit)> = pauli
+        .into_iter()
+        .zip(indices.into_iter().map(Qubit))
+        .collect();
+
+    // get the basis change: x -> HGate, y -> RXGate(pi/2), z -> nothing
+    let basis_change = active_paulis
+        .clone()
+        .into_iter()
+        .filter(|(p, _)| *p != 'z')
+        .map(|(p, q)| match p {
+            'x' => (StandardGate::HGate, smallvec![], smallvec![q]),
+            'y' => (
+                StandardGate::RXGate,
+                smallvec![Param::Float(PI2)],
+                smallvec![q],
+            ),
+            _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
+        });
+
+    // get the inverse basis change
+    let inverse_basis_change = basis_change.clone().map(|(gate, _, qubit)| match gate {
+        StandardGate::HGate => (gate, smallvec![], qubit),
+        StandardGate::RXGate => (gate, smallvec![Param::Float(-PI2)], qubit),
+        _ => unreachable!(),
+    });
+
+    // get the CX chain down to the target rotation qubit
+    let chain_down = active_paulis
+        .clone()
+        .into_iter()
+        .map(|(_, q)| q)
+        .tuple_windows() // iterates over (q[i], q[i+1]) windows
+        .map(|(ctrl, target)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target]));
+
+    // get the CX chain up (cannot use chain_down.rev since tuple_windows is not double ended)
+    let chain_up = active_paulis
+        .clone()
+        .into_iter()
+        .rev()
+        .map(|(_, q)| q)
+        .tuple_windows()
+        .map(|(target, ctrl)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target]));
+
+    // get the RZ gate on the last qubit
+    let last_qubit = active_paulis.last().unwrap().1;
+    let z_rotation = std::iter::once((
+        if phase_gate {
+            StandardGate::PhaseGate
+        } else {
+            StandardGate::RZGate
+        },
+        smallvec![time],
+        smallvec![last_qubit],
+    ));
+
+    // and finally chain everything together
+    basis_change
+        .chain(chain_down)
+        .chain(z_rotation)
+        .chain(chain_up)
+        .chain(inverse_basis_change)
+}
+
+#[pyfunction]
+#[pyo3(signature = (num_qubits, pauli, qubits, time))]
+pub fn py_pauli_evolution(
+    py: Python,
+    num_qubits: i64,
+    pauli: &Bound<PyString>,
+    qubits: &Bound<PyList>,
+    time: &Bound<PyAny>,
+) -> PyResult<CircuitData> {
+    let pauli = pauli.to_string();
+    let indices = qubits
+        .iter()
+        .map(|el| el.extract::<u32>())
+        .collect::<PyResult<_>>()?;
+    let time = Param::extract_no_coerce(time)?;
+    let instructions = pauli_evolution(pauli.as_str(), indices, time, false);
+    let data =
+        CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0));
+    Ok(data)
+}
diff --git a/crates/accelerate/src/circuit_library/pauli_feature_map.rs b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
index 6fa88d187182..39374d4c6d9d 100644
--- a/crates/accelerate/src/circuit_library/pauli_feature_map.rs
+++ b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
@@ -10,7 +10,6 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use itertools::Itertools;
 use pyo3::prelude::*;
 use pyo3::types::PySequence;
 use pyo3::types::PyString;
@@ -24,106 +23,15 @@ use smallvec::{smallvec, SmallVec};
 use std::f64::consts::PI;
 
 use crate::circuit_library::entanglement;
+use crate::circuit_library::pauli_evolution;
 use crate::QiskitError;
 
-// custom math and types for a more readable code
-const PI2: f64 = PI / 2.;
 type Instruction = (
     PackedOperation,
     SmallVec<[Param; 3]>,
     Vec<Qubit>,
     Vec<Clbit>,
 );
-type StandardInstruction = (StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>);
-
-/// Return instructions (using only StandardGate operations) to implement a Pauli evolution
-/// of a given Pauli string over a given time (as Param).
-///
-/// The Pauli evolution is implemented as a basis transformation to the Pauli-Z basis,
-/// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
-/// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
-/// Pauli string XIYZ we have the circuit
-///                     ┌───┐┌───────┐┌───┐
-/// 0: ─────────────────┤ X ├┤ Rz(2) ├┤ X ├──────────────────
-///    ┌──────────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌───────────┐
-/// 1: ┤ Rx(pi/2) ├┤ X ├──■─────────────■──┤ X ├┤ Rx(-pi/2) ├
-///    └──────────┘└─┬─┘                   └─┬─┘└───────────┘
-/// 2: ──────────────┼───────────────────────┼───────────────
-///     ┌───┐        │                       │  ┌───┐
-/// 3: ─┤ H ├────────■───────────────────────■──┤ H ├────────
-///     └───┘                                   └───┘
-fn pauli_evolution(
-    pauli: &str,
-    indices: Vec<u32>,
-    time: Param,
-) -> impl Iterator<Item = StandardInstruction> + '_ {
-    // Get pairs of (pauli, qubit) that are active, i.e. that are not the identity. Note that
-    // the rest of the code also works if there are only identities, in which case we will
-    // effectively return an empty iterator.
-    let qubits = indices.iter().map(|i| Qubit(*i)).collect_vec();
-    let binding = pauli.to_lowercase(); // lowercase for convenience
-    let active_paulis = binding
-        .as_str()
-        .chars()
-        .rev() // reverse due to Qiskit's bit ordering convention
-        .zip(qubits)
-        .filter(|(p, _)| *p != 'i')
-        .collect_vec();
-
-    // get the basis change: x -> HGate, y -> RXGate(pi/2), z -> nothing
-    let basis_change = active_paulis
-        .clone()
-        .into_iter()
-        .filter(|(p, _)| *p != 'z')
-        .map(|(p, q)| match p {
-            'x' => (StandardGate::HGate, smallvec![], smallvec![q]),
-            'y' => (
-                StandardGate::RXGate,
-                smallvec![Param::Float(PI2)],
-                smallvec![q],
-            ),
-            _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
-        });
-
-    // get the inverse basis change
-    let inverse_basis_change = basis_change.clone().map(|(gate, _, qubit)| match gate {
-        StandardGate::HGate => (gate, smallvec![], qubit),
-        StandardGate::RXGate => (gate, smallvec![Param::Float(-PI2)], qubit),
-        _ => unreachable!(),
-    });
-
-    // get the CX chain down to the target rotation qubit
-    let chain_down = active_paulis
-        .clone()
-        .into_iter()
-        .map(|(_, q)| q)
-        .tuple_windows() // iterates over (q[i], q[i+1]) windows
-        .map(|(ctrl, target)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target]));
-
-    // get the CX chain up (cannot use chain_down.rev since tuple_windows is not double ended)
-    let chain_up = active_paulis
-        .clone()
-        .into_iter()
-        .rev()
-        .map(|(_, q)| q)
-        .tuple_windows()
-        .map(|(target, ctrl)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target]));
-
-    // get the RZ gate on the last qubit
-    let last_qubit = active_paulis.last().unwrap().1;
-    let z_rotation = std::iter::once((
-        StandardGate::PhaseGate,
-        smallvec![time],
-        smallvec![last_qubit],
-    ));
-
-    // and finally chain everything together
-    basis_change
-        .chain(chain_down)
-        .chain(z_rotation)
-        .chain(chain_up)
-        .chain(inverse_basis_change)
-}
 
 /// Build a Pauli feature map circuit.
 ///
@@ -263,11 +171,16 @@ fn _get_evolution_layer<'a>(
             // to call CircuitData::from_packed_operations. This is needed since we might
             // have to interject barriers, which are not a standard gate and prevents us
             // from using CircuitData::from_standard_gates.
-            let evo = pauli_evolution(pauli, indices.clone(), multiply_param(&angle, alpha, py))
-                .map(|(gate, params, qargs)| {
-                    (gate.into(), params, qargs.to_vec(), vec![] as Vec<Clbit>)
-                })
-                .collect::<Vec<Instruction>>();
+            let evo = pauli_evolution::pauli_evolution(
+                pauli,
+                indices.clone(),
+                multiply_param(&angle, alpha, py),
+                true,
+            )
+            .map(|(gate, params, qargs)| {
+                (gate.into(), params, qargs.to_vec(), vec![] as Vec<Clbit>)
+            })
+            .collect::<Vec<Instruction>>();
             insts.extend(evo);
         }
     }

From f9b72cf4ec4f3d63c2960fe0f31c38d1f898bf7d Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Fri, 11 Oct 2024 17:32:13 +0200
Subject: [PATCH 07/25] take care of global phase

---
 crates/accelerate/src/circuit_library/mod.rs  |  1 +
 .../src/circuit_library/pauli_evolution.rs    | 54 +++++++++++++------
 crates/circuit/src/operations.rs              |  2 +-
 3 files changed, 41 insertions(+), 16 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/mod.rs b/crates/accelerate/src/circuit_library/mod.rs
index 09692d049f12..3897167f9489 100644
--- a/crates/accelerate/src/circuit_library/mod.rs
+++ b/crates/accelerate/src/circuit_library/mod.rs
@@ -18,6 +18,7 @@ mod pauli_feature_map;
 mod quantum_volume;
 
 pub fn circuit_library(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(pauli_evolution::py_pauli_evolution))?;
     m.add_wrapped(wrap_pyfunction!(pauli_feature_map::pauli_feature_map))?;
     m.add_wrapped(wrap_pyfunction!(entanglement::get_entangler_map))?;
     m.add_wrapped(wrap_pyfunction!(quantum_volume::quantum_volume))?;
diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 7a8db8380258..0f49b917ddac 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -12,8 +12,9 @@
 
 use itertools::Itertools;
 use pyo3::prelude::*;
+use pyo3::types::{PyList, PyString, PyTuple};
 use qiskit_circuit::circuit_data::CircuitData;
-use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::operations::{radd_param, Param, StandardGate};
 use qiskit_circuit::Qubit;
 use smallvec::{smallvec, SmallVec};
 use std::f64::consts::PI;
@@ -59,7 +60,7 @@ pub fn pauli_evolution<'a>(
     let (paulis, indices): (Vec<char>, Vec<u32>) = active.unzip();
 
     match (phase_gate, indices.len()) {
-        (false, 0) => Box::new(std::iter::empty()),
+        (_, 0) => Box::new(std::iter::empty()),
         (false, 1) => Box::new(single_qubit_evolution(paulis[0], indices[0], time)),
         (false, 2) => two_qubit_evolution(paulis.clone(), indices.clone(), time),
         _ => Box::new(multi_qubit_evolution(
@@ -188,22 +189,45 @@ fn multi_qubit_evolution(
 }
 
 #[pyfunction]
-#[pyo3(signature = (num_qubits, pauli, qubits, time))]
+#[pyo3(signature = (num_qubits, sparse_paulis))]
 pub fn py_pauli_evolution(
     py: Python,
     num_qubits: i64,
-    pauli: &Bound<PyString>,
-    qubits: &Bound<PyList>,
-    time: &Bound<PyAny>,
+    sparse_paulis: &Bound<PyList>,
 ) -> PyResult<CircuitData> {
-    let pauli = pauli.to_string();
-    let indices = qubits
+    let num_paulis = sparse_paulis.len();
+    let mut paulis: Vec<String> = Vec::with_capacity(num_paulis);
+    let mut indices: Vec<Vec<u32>> = Vec::with_capacity(num_paulis);
+    let mut times: Vec<Param> = Vec::with_capacity(num_paulis);
+    let mut global_phase = Param::Float(0.0);
+
+    for el in sparse_paulis.iter() {
+        let tuple = el.downcast::<PyTuple>()?;
+        let pauli = tuple.get_item(0)?.downcast::<PyString>()?.to_string();
+        let time = Param::extract_no_coerce(&tuple.get_item(2)?)?;
+
+        if pauli.as_str().chars().all(|p| p == 'i') {
+            global_phase = radd_param(global_phase, time, py); // apply factor -1 at the end
+            continue;
+        }
+
+        paulis.push(pauli);
+        times.push(time);
+        indices.push(
+            tuple
+                .get_item(1)?
+                .downcast::<PyList>()?
+                .iter()
+                .map(|index| index.extract::<u32>())
+                .collect::<PyResult<_>>()?,
+        );
+    }
+
+    let evos = paulis
         .iter()
-        .map(|el| el.extract::<u32>())
-        .collect::<PyResult<_>>()?;
-    let time = Param::extract_no_coerce(time)?;
-    let instructions = pauli_evolution(pauli.as_str(), indices, time, false);
-    let data =
-        CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0));
-    Ok(data)
+        .zip(indices)
+        .zip(times)
+        .flat_map(|((pauli, qubits), time)| pauli_evolution(pauli, qubits, time, false));
+
+    CircuitData::from_standard_gates(py, num_qubits as u32, evos, global_phase)
 }
diff --git a/crates/circuit/src/operations.rs b/crates/circuit/src/operations.rs
index a73d81ad7180..205446e60077 100644
--- a/crates/circuit/src/operations.rs
+++ b/crates/circuit/src/operations.rs
@@ -2116,7 +2116,7 @@ pub fn add_param(param: &Param, summand: f64, py: Python) -> Param {
     }
 }
 
-fn radd_param(param1: Param, param2: Param, py: Python) -> Param {
+pub fn radd_param(param1: Param, param2: Param, py: Python) -> Param {
     match [param1, param2] {
         [Param::Float(theta), Param::Float(lambda)] => Param::Float(theta + lambda),
         [Param::ParameterExpression(theta), Param::ParameterExpression(lambda)] => {

From 39853940308df35ed8fda369e26862690815ab2b Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Fri, 11 Oct 2024 17:52:33 +0200
Subject: [PATCH 08/25] support barriers

---
 crates/accelerate/src/circuit_library/mod.rs  |  1 +
 .../src/circuit_library/pauli_evolution.rs    | 30 ++++++++--
 .../accelerate/src/circuit_library/utils.rs   | 58 +++++++++++++++++++
 qiskit/synthesis/evolution/product_formula.py | 18 +++---
 4 files changed, 96 insertions(+), 11 deletions(-)
 create mode 100644 crates/accelerate/src/circuit_library/utils.rs

diff --git a/crates/accelerate/src/circuit_library/mod.rs b/crates/accelerate/src/circuit_library/mod.rs
index 3897167f9489..b943881e4a14 100644
--- a/crates/accelerate/src/circuit_library/mod.rs
+++ b/crates/accelerate/src/circuit_library/mod.rs
@@ -16,6 +16,7 @@ mod entanglement;
 mod pauli_evolution;
 mod pauli_feature_map;
 mod quantum_volume;
+mod utils;
 
 pub fn circuit_library(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(pauli_evolution::py_pauli_evolution))?;
diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 0f49b917ddac..03b6b94f09f6 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -15,13 +15,22 @@ use pyo3::prelude::*;
 use pyo3::types::{PyList, PyString, PyTuple};
 use qiskit_circuit::circuit_data::CircuitData;
 use qiskit_circuit::operations::{radd_param, Param, StandardGate};
-use qiskit_circuit::Qubit;
+use qiskit_circuit::packed_instruction::PackedOperation;
+use qiskit_circuit::{Clbit, Qubit};
 use smallvec::{smallvec, SmallVec};
 use std::f64::consts::PI;
 
+use crate::circuit_library::utils;
+
 // custom math and types for a more readable code
 const PI2: f64 = PI / 2.;
 type StandardInstruction = (StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>);
+type Instruction = (
+    PackedOperation,
+    SmallVec<[Param; 3]>,
+    Vec<Qubit>,
+    Vec<Clbit>,
+);
 
 /// Return instructions (using only StandardGate operations) to implement a Pauli evolution
 /// of a given Pauli string over a given time (as Param).
@@ -189,11 +198,12 @@ fn multi_qubit_evolution(
 }
 
 #[pyfunction]
-#[pyo3(signature = (num_qubits, sparse_paulis))]
+#[pyo3(signature = (num_qubits, sparse_paulis, insert_barriers=false))]
 pub fn py_pauli_evolution(
     py: Python,
     num_qubits: i64,
     sparse_paulis: &Bound<PyList>,
+    insert_barriers: bool,
 ) -> PyResult<CircuitData> {
     let num_paulis = sparse_paulis.len();
     let mut paulis: Vec<String> = Vec::with_capacity(num_paulis);
@@ -227,7 +237,19 @@ pub fn py_pauli_evolution(
         .iter()
         .zip(indices)
         .zip(times)
-        .flat_map(|((pauli, qubits), time)| pauli_evolution(pauli, qubits, time, false));
+        .flat_map(|((pauli, qubits), time)| {
+            let as_packed = pauli_evolution(pauli, qubits, time, false).map(
+                |(gate, params, qubits)| -> PyResult<Instruction> {
+                    Ok((
+                        gate.into(),
+                        params.clone(),
+                        Vec::from_iter(qubits.into_iter()),
+                        Vec::new(),
+                    ))
+                },
+            );
+            as_packed.chain(utils::maybe_barrier(py, num_qubits as u32, insert_barriers))
+        });
 
-    CircuitData::from_standard_gates(py, num_qubits as u32, evos, global_phase)
+    CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)
 }
diff --git a/crates/accelerate/src/circuit_library/utils.rs b/crates/accelerate/src/circuit_library/utils.rs
new file mode 100644
index 000000000000..4df8d64865f0
--- /dev/null
+++ b/crates/accelerate/src/circuit_library/utils.rs
@@ -0,0 +1,58 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2024
+//
+// 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.
+
+use pyo3::prelude::*;
+use qiskit_circuit::{
+    imports,
+    operations::{Param, PyInstruction},
+    packed_instruction::PackedOperation,
+    Clbit, Qubit,
+};
+use smallvec::{smallvec, SmallVec};
+
+type Instruction = (
+    PackedOperation,
+    SmallVec<[Param; 3]>,
+    Vec<Qubit>,
+    Vec<Clbit>,
+);
+
+/// Get an iterator that returns a barrier or an empty element.
+pub fn maybe_barrier(
+    py: Python,
+    num_qubits: u32,
+    insert_barriers: bool,
+) -> Box<dyn Iterator<Item = PyResult<Instruction>>> {
+    // TODO could speed this up by only defining the barrier class once
+    if !insert_barriers {
+        Box::new(std::iter::empty())
+    } else {
+        let barrier_cls = imports::BARRIER.get_bound(py);
+        let barrier = barrier_cls
+            .call1((num_qubits,))
+            .expect("Could not create Barrier Python-side");
+        let barrier_inst = PyInstruction {
+            qubits: num_qubits,
+            clbits: 0,
+            params: 0,
+            op_name: "barrier".to_string(),
+            control_flow: false,
+            instruction: barrier.into(),
+        };
+        Box::new(std::iter::once(Ok((
+            barrier_inst.into(),
+            smallvec![],
+            (0..num_qubits).map(Qubit).collect(),
+            vec![] as Vec<Clbit>,
+        ))))
+    }
+}
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index 484a29295d81..4f96e461defe 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -23,6 +23,7 @@
 from qiskit.circuit.quantumcircuit import QuantumCircuit, ParameterValueType
 from qiskit.quantum_info import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
+from qiskit._accelerate.circuit_library import py_pauli_evolution
 
 from .evolution_synthesis import EvolutionSynthesis
 
@@ -142,15 +143,18 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
 
         num_qubits = evolution.num_qubits
 
-        circuit = QuantumCircuit(num_qubits)
+        data = py_pauli_evolution(num_qubits, pauli_rotations, self.insert_barriers)
+        circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
 
-        # we could cache the circuit decomposition for the circuits we already saw
-        for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
-            op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
+        # circuit = QuantumCircuit(num_qubits)
 
-            self.atomic_evolution(circuit, op, coeff)
-            if self.insert_barriers and i != len(pauli_rotations) - 1:
-                circuit.barrier()
+        # # we could cache the circuit decomposition for the circuits we already saw
+        # for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
+        #     op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
+
+        #     self.atomic_evolution(circuit, op, coeff)
+        #     if self.insert_barriers and i != len(pauli_rotations) - 1:
+        #         circuit.barrier()
 
         return circuit
 

From 99a115319d23b7fddb92b4e55fd5b2f4960699d3 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Fri, 11 Oct 2024 18:04:55 +0200
Subject: [PATCH 09/25] fix time parameter

some things are still upside down and it seems like it would be cleaner to do the time multiplication inside pauli_evolution
---
 .../src/circuit_library/pauli_evolution.rs    | 22 +++++--------------
 1 file changed, 6 insertions(+), 16 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 03b6b94f09f6..22e7845c2a18 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -14,7 +14,7 @@ use itertools::Itertools;
 use pyo3::prelude::*;
 use pyo3::types::{PyList, PyString, PyTuple};
 use qiskit_circuit::circuit_data::CircuitData;
-use qiskit_circuit::operations::{radd_param, Param, StandardGate};
+use qiskit_circuit::operations::{multiply_param, radd_param, Param, StandardGate};
 use qiskit_circuit::packed_instruction::PackedOperation;
 use qiskit_circuit::{Clbit, Qubit};
 use smallvec::{smallvec, SmallVec};
@@ -91,8 +91,8 @@ fn single_qubit_evolution(
 
     std::iter::once(match pauli {
         'x' => (StandardGate::RXGate, param, qubit),
-        'y' => (StandardGate::RXGate, param, qubit),
-        'z' => (StandardGate::RXGate, param, qubit),
+        'y' => (StandardGate::RYGate, param, qubit),
+        'z' => (StandardGate::RZGate, param, qubit),
         _ => unreachable!("Unsupported Pauli, at this point we expected one of x, y, z."),
     })
 }
@@ -121,18 +121,6 @@ fn multi_qubit_evolution(
     time: Param,
     phase_gate: bool,
 ) -> impl Iterator<Item = StandardInstruction> {
-    // Get pairs of (pauli, qubit) that are active, i.e. that are not the identity. Note that
-    // the rest of the code also works if there are only identities, in which case we will
-    // effectively return an empty iterator.
-    // let qubits: Vec<Qubit> = indices.iter().map(|i| Qubit(*i)).collect();
-    // let binding = pauli.to_lowercase(); // lowercase for convenience
-    // let active_paulis = binding
-    //     .as_str()
-    //     .chars()
-    //     .rev() // reverse due to Qiskit's bit ordering convention
-    //     .zip(qubits)
-    //     .filter(|(p, _)| *p != 'i')
-    //     .collect_vec();
     let active_paulis: Vec<(char, Qubit)> = pauli
         .into_iter()
         .zip(indices.into_iter().map(Qubit))
@@ -222,7 +210,7 @@ pub fn py_pauli_evolution(
         }
 
         paulis.push(pauli);
-        times.push(time);
+        times.push(multiply_param(&time, 2.0, py));
         indices.push(
             tuple
                 .get_item(1)?
@@ -251,5 +239,7 @@ pub fn py_pauli_evolution(
             as_packed.chain(utils::maybe_barrier(py, num_qubits as u32, insert_barriers))
         });
 
+    // apply factor -1 for global phase
+    global_phase = multiply_param(&global_phase, -1.0, py);
     CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)
 }

From 214b96961bbce585dd6eb662608ac2dbdf07b2f9 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 14 Oct 2024 14:20:25 +0200
Subject: [PATCH 10/25] fix lint

---
 qiskit/synthesis/evolution/product_formula.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index 484a29295d81..a378d6b275f4 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -6,7 +6,7 @@
 # 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.
 #
-# typing.Any modifications or derivative works of this code must retain this
+# 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.
 

From e8bcc4f5bab6d593ca10fca2cef8e9f1a43eecda Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 14 Oct 2024 14:27:11 +0200
Subject: [PATCH 11/25] fix final barrier

---
 .../src/circuit_library/pauli_evolution.rs       | 16 +++++++++-------
 qiskit/synthesis/evolution/product_formula.py    | 10 ----------
 2 files changed, 9 insertions(+), 17 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 22e7845c2a18..2db1d6592382 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -221,11 +221,8 @@ pub fn py_pauli_evolution(
         );
     }
 
-    let evos = paulis
-        .iter()
-        .zip(indices)
-        .zip(times)
-        .flat_map(|((pauli, qubits), time)| {
+    let evos = paulis.iter().enumerate().zip(indices).zip(times).flat_map(
+        |(((i, pauli), qubits), time)| {
             let as_packed = pauli_evolution(pauli, qubits, time, false).map(
                 |(gate, params, qubits)| -> PyResult<Instruction> {
                     Ok((
@@ -236,8 +233,13 @@ pub fn py_pauli_evolution(
                     ))
                 },
             );
-            as_packed.chain(utils::maybe_barrier(py, num_qubits as u32, insert_barriers))
-        });
+            as_packed.chain(utils::maybe_barrier(
+                py,
+                num_qubits as u32,
+                insert_barriers && i < num_paulis, // do not add barrier on final block
+            ))
+        },
+    );
 
     // apply factor -1 for global phase
     global_phase = multiply_param(&global_phase, -1.0, py);
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index c0473d8f32ca..663fc7f3046d 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -146,16 +146,6 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
         data = py_pauli_evolution(num_qubits, pauli_rotations, self.insert_barriers)
         circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
 
-        # circuit = QuantumCircuit(num_qubits)
-
-        # # we could cache the circuit decomposition for the circuits we already saw
-        # for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
-        #     op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
-
-        #     self.atomic_evolution(circuit, op, coeff)
-        #     if self.insert_barriers and i != len(pauli_rotations) - 1:
-        #         circuit.barrier()
-
         return circuit
 
     @property

From aba92e279968a338ec2f133db60bae6b8439af89 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 15 Oct 2024 13:00:19 +0200
Subject: [PATCH 12/25] fix wrapping

---
 .../src/circuit_library/pauli_evolution.rs    | 157 ++++++---
 .../src/circuit_library/pauli_feature_map.rs  |   1 +
 qiskit/synthesis/evolution/product_formula.py | 316 ++++--------------
 qiskit/synthesis/evolution/suzuki_trotter.py  |  29 +-
 .../circuit/library/test_evolution_gate.py    |  59 +++-
 5 files changed, 252 insertions(+), 310 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 2db1d6592382..15c725de218e 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -10,7 +10,6 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use itertools::Itertools;
 use pyo3::prelude::*;
 use pyo3::types::{PyList, PyString, PyTuple};
 use qiskit_circuit::circuit_data::CircuitData;
@@ -35,35 +34,32 @@ type Instruction = (
 /// Return instructions (using only StandardGate operations) to implement a Pauli evolution
 /// of a given Pauli string over a given time (as Param).
 ///
-/// The Pauli evolution is implemented as a basis transformation to the Pauli-Z basis,
-/// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
-/// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
-/// Pauli string XIYZ we have the circuit
-///                     ┌───┐┌───────┐┌───┐
-/// 0: ─────────────────┤ X ├┤ Rz(2) ├┤ X ├──────────────────
-///    ┌──────────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌───────────┐
-/// 1: ┤ Rx(pi/2) ├┤ X ├──■─────────────■──┤ X ├┤ Rx(-pi/2) ├
-///    └──────────┘└─┬─┘                   └─┬─┘└───────────┘
-/// 2: ──────────────┼───────────────────────┼───────────────
-///     ┌───┐        │                       │  ┌───┐
-/// 3: ─┤ H ├────────■───────────────────────■──┤ H ├────────
-///     └───┘                                   └───┘
+/// Args:
+///     pauli: The Pauli string, e.g. "IXYZ".
+///     indices: The qubit indices the Pauli acts on, e.g. if given as [0, 1, 2, 3] with the
+///         Pauli "IXYZ", then the correspondence is I_0 X_1 Y_2 Z_3.
+///     time: The rotation angle. Note that this will directly be used as input of the
+///         rotation gate and not be multiplied by a factor of 2 (that should be done before so
+///         that this function can remain Rust-only).
+///     phase_gate: If ``true``, use the ``PhaseGate`` instead of ``RZGate`` as single-qubit rotation.
+///     cx_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
+///         CX uses the top qubit as target. If ``false``, uses a "chain" shape, where CX in between
+///         neighboring qubits are used.
 ///
-/// If ``phase_gate`` is ``true``, use the ``PhaseGate`` instead of ``RZGate`` as single-qubit
-/// rotation.
+/// Returns:
+///     A pointer to an iterator over standard instructions.
 pub fn pauli_evolution<'a>(
     pauli: &'a str,
     indices: Vec<u32>,
     time: Param,
     phase_gate: bool,
-    // check_sparse: bool
+    do_fountain: bool,
 ) -> Box<dyn Iterator<Item = StandardInstruction> + 'a> {
     // ensure the Pauli has no identity terms
     let binding = pauli.to_lowercase(); // lowercase for convenience
     let active = binding
         .as_str()
         .chars()
-        .rev()
         .zip(indices)
         .filter(|(pauli, _)| *pauli != 'i');
     let (paulis, indices): (Vec<char>, Vec<u32>) = active.unzip();
@@ -71,23 +67,27 @@ pub fn pauli_evolution<'a>(
     match (phase_gate, indices.len()) {
         (_, 0) => Box::new(std::iter::empty()),
         (false, 1) => Box::new(single_qubit_evolution(paulis[0], indices[0], time)),
-        (false, 2) => two_qubit_evolution(paulis.clone(), indices.clone(), time),
+        (false, 2) => two_qubit_evolution(paulis, indices, time),
         _ => Box::new(multi_qubit_evolution(
-            paulis.clone(),
-            indices.clone(),
+            paulis,
+            indices,
             time,
             phase_gate,
+            do_fountain,
         )),
     }
 }
 
+/// Implement a single-qubit Pauli evolution of a Pauli given as char, on a given index and
+/// for given time. Note that the time here equals the angle of the rotation and is not
+/// multiplied by a factor of 2.
 fn single_qubit_evolution(
     pauli: char,
     index: u32,
     time: Param,
 ) -> impl Iterator<Item = StandardInstruction> {
     let qubit: SmallVec<[Qubit; 2]> = smallvec![Qubit(index)];
-    let param: SmallVec<[Param; 3]> = smallvec![time.clone()];
+    let param: SmallVec<[Param; 3]> = smallvec![time];
 
     std::iter::once(match pauli {
         'x' => (StandardGate::RXGate, param, qubit),
@@ -97,6 +97,12 @@ fn single_qubit_evolution(
     })
 }
 
+/// Implement a 2-qubit Pauli evolution of a Pauli string, on a given indices and
+/// for given time. Note that the time here equals the angle of the rotation and is not
+/// multiplied by a factor of 2.
+///
+/// If possible, Qiskit's native 2-qubit Pauli rotations are used. Otherwise, the general
+/// multi-qubit evolution is called.
 fn two_qubit_evolution<'a>(
     pauli: Vec<char>,
     indices: Vec<u32>,
@@ -111,15 +117,19 @@ fn two_qubit_evolution<'a>(
         "zx" => Box::new(std::iter::once((StandardGate::RZXGate, param, qubits))),
         "yy" => Box::new(std::iter::once((StandardGate::RYYGate, param, qubits))),
         "zz" => Box::new(std::iter::once((StandardGate::RZZGate, param, qubits))),
-        _ => Box::new(multi_qubit_evolution(pauli, indices, time, false)),
+        // note that the CX modes (do_fountain=true/false) give the same circuit for a 2-qubit
+        // Pauli, so we just set it to false here
+        _ => Box::new(multi_qubit_evolution(pauli, indices, time, false, false)),
     }
 }
 
+/// Implement a multi-qubit Pauli evolution. See ``pauli_evolution`` detailed docs.
 fn multi_qubit_evolution(
     pauli: Vec<char>,
     indices: Vec<u32>,
     time: Param,
     phase_gate: bool,
+    do_fountain: bool,
 ) -> impl Iterator<Item = StandardInstruction> {
     let active_paulis: Vec<(char, Qubit)> = pauli
         .into_iter()
@@ -147,26 +157,22 @@ fn multi_qubit_evolution(
         StandardGate::RXGate => (gate, smallvec![Param::Float(-PI2)], qubit),
         _ => unreachable!(),
     });
+    println!("Fountain: {}", do_fountain);
 
-    // get the CX chain down to the target rotation qubit
-    let chain_down = active_paulis
-        .clone()
-        .into_iter()
-        .map(|(_, q)| q)
-        .tuple_windows() // iterates over (q[i], q[i+1]) windows
-        .map(|(ctrl, target)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target]));
-
-    // get the CX chain up (cannot use chain_down.rev since tuple_windows is not double ended)
-    let chain_up = active_paulis
-        .clone()
-        .into_iter()
-        .rev()
-        .map(|(_, q)| q)
-        .tuple_windows()
-        .map(|(target, ctrl)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target]));
+    // get the CX propagation up to the first qubit, and down
+    let (chain_up, chain_down) = match do_fountain {
+        true => (
+            cx_fountain(active_paulis.clone()),
+            cx_fountain(active_paulis.clone()).rev(),
+        ),
+        false => (
+            cx_chain(active_paulis.clone()),
+            cx_chain(active_paulis.clone()).rev(),
+        ),
+    };
 
-    // get the RZ gate on the last qubit
-    let last_qubit = active_paulis.last().unwrap().1;
+    // get the RZ gate on the first qubit
+    let first_qubit = active_paulis.first().unwrap().1;
     let z_rotation = std::iter::once((
         if phase_gate {
             StandardGate::PhaseGate
@@ -174,9 +180,11 @@ fn multi_qubit_evolution(
             StandardGate::RZGate
         },
         smallvec![time],
-        smallvec![last_qubit],
+        smallvec![first_qubit],
     ));
 
+    println!("first qubit {:?}", first_qubit);
+
     // and finally chain everything together
     basis_change
         .chain(chain_down)
@@ -185,13 +193,44 @@ fn multi_qubit_evolution(
         .chain(inverse_basis_change)
 }
 
+/// Implement a Pauli evolution circuit.
+///
+/// The Pauli evolution is implemented as a basis transformation to the Pauli-Z basis,
+/// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
+/// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
+/// Pauli string XIYZ we have the circuit
+///                     ┌───┐┌───────┐┌───┐
+/// 0: ─────────────────┤ X ├┤ Rz(2) ├┤ X ├──────────────────
+///    ┌──────────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌───────────┐
+/// 1: ┤ Rx(pi/2) ├┤ X ├──■─────────────■──┤ X ├┤ Rx(-pi/2) ├
+///    └──────────┘└─┬─┘                   └─┬─┘└───────────┘
+/// 2: ──────────────┼───────────────────────┼───────────────
+///     ┌───┐        │                       │  ┌───┐
+/// 3: ─┤ H ├────────■───────────────────────■──┤ H ├────────
+///     └───┘                                   └───┘
+///
+/// Args:
+///     num_qubits: The number of qubits in the Hamiltonian.
+///     sparse_paulis: The Paulis to implement. Given in a sparse-list format with elements
+///         ``(pauli_string, qubit_indices, coefficient)``. An element of the form
+///         ``("IXYZ", [0,1,2,3], 0.2)``, for example, is interpreted in terms of qubit indices as
+///          I_q0 X_q1 Y_q2 Z_q3 and will use a RZ rotation angle of 0.4.
+///     insert_barriers: If ``true``, insert a barrier in between the evolution of individual
+///         Pauli terms.
+///     cx_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
+///         CX uses the top qubit as target. If ``false``, uses a "chain" shape, where CX in between
+///         neighboring qubits are used.
+///
+/// Returns:
+///     Circuit data for to implement the evolution.
 #[pyfunction]
-#[pyo3(signature = (num_qubits, sparse_paulis, insert_barriers=false))]
+#[pyo3(signature = (num_qubits, sparse_paulis, insert_barriers=false, do_fountain=false))]
 pub fn py_pauli_evolution(
     py: Python,
     num_qubits: i64,
     sparse_paulis: &Bound<PyList>,
     insert_barriers: bool,
+    do_fountain: bool,
 ) -> PyResult<CircuitData> {
     let num_paulis = sparse_paulis.len();
     let mut paulis: Vec<String> = Vec::with_capacity(num_paulis);
@@ -223,11 +262,11 @@ pub fn py_pauli_evolution(
 
     let evos = paulis.iter().enumerate().zip(indices).zip(times).flat_map(
         |(((i, pauli), qubits), time)| {
-            let as_packed = pauli_evolution(pauli, qubits, time, false).map(
+            let as_packed = pauli_evolution(pauli, qubits, time, false, do_fountain).map(
                 |(gate, params, qubits)| -> PyResult<Instruction> {
                     Ok((
                         gate.into(),
-                        params.clone(),
+                        params,
                         Vec::from_iter(qubits.into_iter()),
                         Vec::new(),
                     ))
@@ -236,7 +275,7 @@ pub fn py_pauli_evolution(
             as_packed.chain(utils::maybe_barrier(
                 py,
                 num_qubits as u32,
-                insert_barriers && i < num_paulis, // do not add barrier on final block
+                insert_barriers && i < (num_paulis - 1), // do not add barrier on final block
             ))
         },
     );
@@ -245,3 +284,29 @@ pub fn py_pauli_evolution(
     global_phase = multiply_param(&global_phase, -1.0, py);
     CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)
 }
+
+fn cx_chain(
+    active_paulis: Vec<(char, Qubit)>,
+) -> Box<dyn DoubleEndedIterator<Item = StandardInstruction>> {
+    let num_terms = active_paulis.len();
+    Box::new(
+        (0..num_terms - 1)
+            .map(move |i| (active_paulis[i].1, active_paulis[i + 1].1))
+            .map(|(target, ctrl)| (StandardGate::CXGate, smallvec![], smallvec![ctrl, target])),
+    )
+}
+
+fn cx_fountain(
+    active_paulis: Vec<(char, Qubit)>,
+) -> Box<dyn DoubleEndedIterator<Item = StandardInstruction>> {
+    let num_terms = active_paulis.len();
+    let first_qubit = active_paulis[0].1;
+    Box::new((1..num_terms).rev().map(move |i| {
+        let ctrl = active_paulis[i].1;
+        (
+            StandardGate::CXGate,
+            smallvec![],
+            smallvec![ctrl, first_qubit],
+        )
+    }))
+}
diff --git a/crates/accelerate/src/circuit_library/pauli_feature_map.rs b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
index 39374d4c6d9d..b6b56fd55f29 100644
--- a/crates/accelerate/src/circuit_library/pauli_feature_map.rs
+++ b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
@@ -176,6 +176,7 @@ fn _get_evolution_layer<'a>(
                 indices.clone(),
                 multiply_param(&angle, alpha, py),
                 true,
+                false,
             )
             .map(|(gate, params, qargs)| {
                 (gate.into(), params, qargs.to_vec(), vec![] as Vec<Clbit>)
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index 663fc7f3046d..aadea5d6de00 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -16,7 +16,6 @@
 
 import inspect
 from collections.abc import Callable
-from functools import partial
 import typing
 import numpy as np
 from qiskit.circuit.parameterexpression import ParameterExpression
@@ -82,8 +81,9 @@ def __init__(
                 Alternatively, the function can also take Pauli operator and evolution time as
                 inputs and returns the circuit that will be appended to the overall circuit being
                 built.
-            wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
-                effect when ``atomic_evolution is None``.
+            wrap: Whether to wrap the atomic evolutions into custom gate objects. Note that setting
+                this to ``True`` is slower than ``False``. This only takes effect when
+                ``atomic_evolution is None``.
         """
         super().__init__()
         self.order = order
@@ -92,14 +92,18 @@ def __init__(
 
         # user-provided atomic evolution, stored for serialization
         self._atomic_evolution = atomic_evolution
-        self._cx_structure = cx_structure
         self._wrap = wrap
 
+        if cx_structure == "chain":
+            self._cx_fountain = False
+        elif cx_structure == "fountain":
+            self._cx_fountain = True
+        else:
+            raise ValueError(f"Unsupported CX structure: {cx_structure}")
+
         # if atomic evolution is not provided, set a default
         if atomic_evolution is None:
-            self.atomic_evolution = partial(
-                _default_atomic_evolution, cx_structure=cx_structure, wrap=wrap
-            )
+            self.atomic_evolution = None
 
         elif len(inspect.signature(atomic_evolution).parameters) == 2:
 
@@ -140,11 +144,18 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
             QuantumCircuit: A circuit implementing the evolution.
         """
         pauli_rotations = self.expand(evolution)
-
         num_qubits = evolution.num_qubits
 
-        data = py_pauli_evolution(num_qubits, pauli_rotations, self.insert_barriers)
-        circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
+        if self._wrap or self._atomic_evolution is not None:
+            # this is the slow path, where each Pauli evolution is constructed in Rust
+            # separately and then wrapped into a gate object
+            circuit = self._custom_evolution(num_qubits, pauli_rotations)
+        else:
+            # this is the fast path, where the whole evolution is constructed Rust-side
+            data = py_pauli_evolution(
+                num_qubits, pauli_rotations, self.insert_barriers, self._cx_fountain
+            )
+            circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
 
         return circuit
 
@@ -177,257 +188,44 @@ def _normalize_coefficients(
         """Ensure the coefficients are real (or parameter expressions)."""
         return [[(op, qubits, real_or_fail(coeff)) for op, qubits, coeff in ops] for ops in paulis]
 
+    def _custom_evolution(self, num_qubits, pauli_rotations):
+        """Implement the evolution for the non-standard path.
 
-def evolve_pauli(
-    output: QuantumCircuit,
-    pauli: Pauli,
-    time: float | ParameterExpression = 1.0,
-    cx_structure: str = "chain",
-    wrap: bool = False,
-    label: str | None = None,
-) -> None:
-    r"""Construct a circuit implementing the time evolution of a single Pauli string.
-
-    For a Pauli string :math:`P = \{I, X, Y, Z\}^{\otimes n}` on :math:`n` qubits and an
-    evolution time :math:`t`, the returned circuit implements the unitary operation
-
-    .. math::
-
-        U(t) = e^{-itP}.
-
-    Since only a single Pauli string is evolved the circuit decomposition is exact.
-
-    Args:
-        output: The circuit object to which to append the evolved Pauli.
-        pauli: The Pauli to evolve.
-        time: The evolution time.
-        cx_structure: Determine the structure of CX gates, can be either ``"chain"`` for
-            next-neighbor connections or ``"fountain"`` to connect directly to the top qubit.
-        wrap: Whether to wrap the single Pauli evolutions into custom gate objects.
-        label: A label for the gate.
-    """
-    num_non_identity = len([label for label in pauli.to_label() if label != "I"])
-
-    # first check, if the Pauli is only the identity, in which case the evolution only
-    # adds a global phase
-    if num_non_identity == 0:
-        output.global_phase -= time
-    # if we evolve on a single qubit, if yes use the corresponding qubit rotation
-    elif num_non_identity == 1:
-        _single_qubit_evolution(output, pauli, time, wrap)
-    # same for two qubits, use Qiskit's native rotations
-    elif num_non_identity == 2:
-        _two_qubit_evolution(output, pauli, time, cx_structure, wrap)
-    # otherwise do basis transformation and CX chains
-    else:
-        _multi_qubit_evolution(output, pauli, time, cx_structure, wrap)
-
-
-def _single_qubit_evolution(output, pauli, time, wrap):
-    dest = QuantumCircuit(1) if wrap else output
-    # Note that all phases are removed from the pauli label and are only in the coefficients.
-    # That's because the operators we evolved have all been translated to a SparsePauliOp.
-    qubits = []
-    label = ""
-    for i, pauli_i in enumerate(reversed(pauli.to_label())):
-        idx = 0 if wrap else i
-        if pauli_i == "X":
-            dest.rx(2 * time, idx)
-            qubits.append(i)
-            label += "X"
-        elif pauli_i == "Y":
-            dest.ry(2 * time, idx)
-            qubits.append(i)
-            label += "Y"
-        elif pauli_i == "Z":
-            dest.rz(2 * time, idx)
-            qubits.append(i)
-            label += "Z"
-
-    if wrap:
-        gate = dest.to_gate(label=f"exp(it {label})")
-        qubits = [output.qubits[q] for q in qubits]
-        output.append(gate, qargs=qubits, copy=False)
-
-
-def _two_qubit_evolution(output, pauli, time, cx_structure, wrap):
-    # Get the Paulis and the qubits they act on.
-    # Note that all phases are removed from the pauli label and are only in the coefficients.
-    # That's because the operators we evolved have all been translated to a SparsePauliOp.
-    labels_as_array = np.array(list(reversed(pauli.to_label())))
-    qubits = np.where(labels_as_array != "I")[0]
-    indices = [0, 1] if wrap else qubits
-    labels = np.array([labels_as_array[idx] for idx in qubits])
-
-    dest = QuantumCircuit(2) if wrap else output
-
-    # go through all cases we have implemented in Qiskit
-    if all(labels == "X"):  # RXX
-        dest.rxx(2 * time, indices[0], indices[1])
-    elif all(labels == "Y"):  # RYY
-        dest.ryy(2 * time, indices[0], indices[1])
-    elif all(labels == "Z"):  # RZZ
-        dest.rzz(2 * time, indices[0], indices[1])
-    elif labels[0] == "Z" and labels[1] == "X":  # RZX
-        dest.rzx(2 * time, indices[0], indices[1])
-    elif labels[0] == "X" and labels[1] == "Z":  # RXZ
-        dest.rzx(2 * time, indices[1], indices[0])
-    else:  # all the others are not native in Qiskit, so use default the decomposition
-        _multi_qubit_evolution(output, pauli, time, cx_structure, wrap)
-        return
-
-    if wrap:
-        gate = dest.to_gate(label=f"exp(it {''.join(labels)})")
-        qubits = [output.qubits[q] for q in qubits]
-        output.append(gate, qargs=qubits, copy=False)
-
-
-def _multi_qubit_evolution(output, pauli, time, cx_structure, wrap):
-    # get diagonalizing clifford
-    cliff = diagonalizing_clifford(pauli)
-
-    # get CX chain to reduce the evolution to the top qubit
-    if cx_structure == "chain":
-        chain = cnot_chain(pauli)
-    else:
-        chain = cnot_fountain(pauli)
-
-    # determine qubit to do the rotation on
-    target = None
-    # Note that all phases are removed from the pauli label and are only in the coefficients.
-    # That's because the operators we evolved have all been translated to a SparsePauliOp.
-    for i, pauli_i in enumerate(reversed(pauli.to_label())):
-        if pauli_i != "I":
-            target = i
-            break
-
-    # build the evolution as: diagonalization, reduction, 1q evolution, followed by inverses
-    dest = QuantumCircuit(pauli.num_qubits) if wrap else output
-    dest.compose(cliff, inplace=True)
-    dest.compose(chain, inplace=True)
-    dest.rz(2 * time, target)
-    dest.compose(chain.inverse(), inplace=True)
-    dest.compose(cliff.inverse(), inplace=True)
-
-    if wrap:
-        gate = dest.to_gate(label=f"exp(it {pauli.to_label()})")
-        output.append(gate, qargs=output.qubits, copy=False)
-
-
-def diagonalizing_clifford(pauli: Pauli) -> QuantumCircuit:
-    """Get the clifford circuit to diagonalize the Pauli operator.
-
-    Args:
-        pauli: The Pauli to diagonalize.
-
-    Returns:
-        A circuit to diagonalize.
-    """
-    cliff = QuantumCircuit(pauli.num_qubits)
-    for i, pauli_i in enumerate(reversed(pauli.to_label())):
-        if pauli_i == "Y":
-            cliff.sdg(i)
-        if pauli_i in ["X", "Y"]:
-            cliff.h(i)
-
-    return cliff
-
-
-def cnot_chain(pauli: Pauli) -> QuantumCircuit:
-    """CX chain.
-
-    For example, for the Pauli with the label 'XYZIX'.
-
-    .. parsed-literal::
-
-                       ┌───┐
-        q_0: ──────────┤ X ├
-                       └─┬─┘
-        q_1: ────────────┼──
-                  ┌───┐  │
-        q_2: ─────┤ X ├──■──
-             ┌───┐└─┬─┘
-        q_3: ┤ X ├──■───────
-             └─┬─┘
-        q_4: ──■────────────
-
-    Args:
-        pauli: The Pauli for which to construct the CX chain.
-
-    Returns:
-        A circuit implementing the CX chain.
-    """
-
-    chain = QuantumCircuit(pauli.num_qubits)
-    control, target = None, None
-
-    # iterate over the Pauli's and add CNOTs
-    for i, pauli_i in enumerate(pauli.to_label()):
-        i = pauli.num_qubits - i - 1
-        if pauli_i != "I":
-            if control is None:
-                control = i
-            else:
-                target = i
-
-        if control is not None and target is not None:
-            chain.cx(control, target)
-            control = i
-            target = None
-
-    return chain
-
-
-def cnot_fountain(pauli: Pauli) -> QuantumCircuit:
-    """CX chain in the fountain shape.
-
-    For example, for the Pauli with the label 'XYZIX'.
-
-    .. parsed-literal::
-
-             ┌───┐┌───┐┌───┐
-        q_0: ┤ X ├┤ X ├┤ X ├
-             └─┬─┘└─┬─┘└─┬─┘
-        q_1: ──┼────┼────┼──
-               │    │    │
-        q_2: ──■────┼────┼──
-                    │    │
-        q_3: ───────■────┼──
-                         │
-        q_4: ────────────■──
-
-    Args:
-        pauli: The Pauli for which to construct the CX chain.
-
-    Returns:
-        A circuit implementing the CX chain.
-    """
+        This is either because a user-defined atomic evolution is given, or because the evolution
+        of individual Paulis needs to be wrapped in gates.
+        """
+        circuit = QuantumCircuit(num_qubits)
+
+        num_paulis = len(pauli_rotations)
+        for i, pauli_rotation in enumerate(pauli_rotations):
+            if self._atomic_evolution is not None:
+                # use the user-provided evolution with a global operator
+                operator = SparsePauliOp.from_sparse_list([pauli_rotation], num_qubits)
+                self.atomic_evolution(circuit, operator, time=1)  # time is inside the Pauli coeff
+
+            else:  # this means self._wrap is True
+                # we create a local sparse Pauli representation such that the operator
+                # does not span over all qubits of the circuit
+                pauli_string, qubits, coeff = pauli_rotation
+                local_pauli = (pauli_string, list(range(len(qubits))), coeff)
+
+                # build the circuit Rust-side
+                data = py_pauli_evolution(
+                    len(qubits),
+                    [local_pauli],
+                    False,
+                    self._cx_fountain,
+                )
+                evo = QuantumCircuit._from_circuit_data(data)
+
+                # and append it to the circuit with the correct label
+                gate = evo.to_gate(label=f"exp(it {pauli_string})")
+                circuit.append(gate, qubits)
+
+            if self.insert_barriers and i < num_paulis - 1:
+                circuit.barrier()
 
-    chain = QuantumCircuit(pauli.num_qubits)
-    control, target = None, None
-    for i, pauli_i in enumerate(reversed(pauli.to_label())):
-        if pauli_i != "I":
-            if target is None:
-                target = i
-            else:
-                control = i
-
-        if control is not None and target is not None:
-            chain.cx(control, target)
-            control = None
-
-    return chain
-
-
-def _default_atomic_evolution(output, operator, time, cx_structure, wrap):
-    if isinstance(operator, Pauli):
-        # single Pauli operator: just exponentiate it
-        evolve_pauli(output, operator, time, cx_structure, wrap)
-    else:
-        # sum of Pauli operators: exponentiate each term (this assumes they commute)
-        pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operator.to_list()]
-        for pauli, coeff in pauli_list:
-            evolve_pauli(output, pauli, coeff * time, cx_structure, wrap)
+        return circuit
 
 
 def real_or_fail(value, tol=100):
diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index 44f66f3d5f3b..0201c5367dbe 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -15,6 +15,7 @@
 from __future__ import annotations
 
 import inspect
+import typing
 from collections.abc import Callable
 from itertools import chain
 
@@ -24,6 +25,10 @@
 
 from .product_formula import ProductFormula
 
+if typing.TYPE_CHECKING:
+    from qiskit.circuit.quantumcircuit import ParameterValueType
+    from qiskit.circuit.library.pauli_evolution import PauliEvolutionGate
+
 
 class SuzukiTrotter(ProductFormula):
     r"""The (higher order) Suzuki-Trotter product formula.
@@ -110,11 +115,26 @@ def __init__(
             )
         super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
-    def expand(self, evolution):
-        """
-        H = ZZ + IX --> ("X", [0], 1/2), ("ZZ", [0, 1], 1), ("X", [0], 1/2)
+    def expand(
+        self, evolution: PauliEvolutionGate
+    ) -> list[tuple[str, list[int], ParameterValueType]]:
+        """Expand the Hamiltonian into a Suzuki-Trotter sequence of sparse gates.
+
+        For example, the Hamiltonian ``H = IX + ZZ`` for an evolution time ``t`` and
+        1 repetition for an order 2 formula would get decomposed into
+
+        .. code-block:: text
+
+            ("X", [0], t/2), ("ZZ", [0, 1], t), ("X", [0], t/2)
 
-        ("X", [0], 1/2), ("ZZ", [0, 1], 1), ("X", [0], 1), ("ZZ", [0, 1], 1), ("X", [0], 1/2)
+        For ``N`` repetitions, this sequence would be repeated ``N`` times and the coefficients
+        divided by ``N``.
+
+        Args:
+            evolution: The evolution gate to expand.
+
+        Returns:
+            The Pauli network implementing the Trotter expansion.
         """
         operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
         time = evolution.time
@@ -135,6 +155,7 @@ def expand(self, evolution):
         # we're already done here since Lie Trotter does not do any operator repetition
         product_formula = self._recurse(self.order, non_commuting)
         flattened = self.reps * list(chain.from_iterable(product_formula))
+
         return flattened
 
     @staticmethod
diff --git a/test/python/circuit/library/test_evolution_gate.py b/test/python/circuit/library/test_evolution_gate.py
index 94825729a70b..3b94003af35b 100644
--- a/test/python/circuit/library/test_evolution_gate.py
+++ b/test/python/circuit/library/test_evolution_gate.py
@@ -20,7 +20,6 @@
 from qiskit.circuit import QuantumCircuit, Parameter
 from qiskit.circuit.library import PauliEvolutionGate
 from qiskit.synthesis import LieTrotter, SuzukiTrotter, MatrixExponential, QDrift
-from qiskit.synthesis.evolution.product_formula import cnot_chain, diagonalizing_clifford
 from qiskit.converters import circuit_to_dag
 from qiskit.quantum_info import Operator, SparsePauliOp, Pauli, Statevector
 from test import QiskitTestCase  # pylint: disable=wrong-import-order
@@ -337,6 +336,12 @@ def test_atomic_evolution(self):
         """Test a custom atomic_evolution."""
 
         def atomic_evolution(pauli, time):
+            if isinstance(pauli, SparsePauliOp):
+                if len(pauli.paulis) != 1:
+                    raise ValueError("Unsupported input.")
+                time *= np.real(pauli.coeffs[0])
+                pauli = pauli.paulis[0]
+
             cliff = diagonalizing_clifford(pauli)
             chain = cnot_chain(pauli)
 
@@ -366,5 +371,57 @@ def atomic_evolution(pauli, time):
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)
 
 
+def diagonalizing_clifford(pauli: Pauli) -> QuantumCircuit:
+    """Get the clifford circuit to diagonalize the Pauli operator."""
+    cliff = QuantumCircuit(pauli.num_qubits)
+    for i, pauli_i in enumerate(reversed(pauli.to_label())):
+        if pauli_i == "Y":
+            cliff.sdg(i)
+        if pauli_i in ["X", "Y"]:
+            cliff.h(i)
+
+    return cliff
+
+
+def cnot_chain(pauli: Pauli) -> QuantumCircuit:
+    """CX chain.
+
+    For example, for the Pauli with the label 'XYZIX'.
+
+    .. parsed-literal::
+
+                       ┌───┐
+        q_0: ──────────┤ X ├
+                       └─┬─┘
+        q_1: ────────────┼──
+                  ┌───┐  │
+        q_2: ─────┤ X ├──■──
+             ┌───┐└─┬─┘
+        q_3: ┤ X ├──■───────
+             └─┬─┘
+        q_4: ──■────────────
+
+    """
+
+    chain = QuantumCircuit(pauli.num_qubits)
+    control, target = None, None
+
+    # iterate over the Pauli's and add CNOTs
+    for i, pauli_i in enumerate(pauli.to_label()):
+        i = pauli.num_qubits - i - 1
+        if pauli_i != "I":
+            if control is None:
+                control = i
+            else:
+                target = i
+
+        if control is not None and target is not None:
+            chain.cx(control, target)
+            control = i
+            target = None
+
+    return chain
+
+
 if __name__ == "__main__":
     unittest.main()

From 4fcf6396396f479aea70c0de34ad584e1fd5c93f Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 15 Oct 2024 13:34:59 +0200
Subject: [PATCH 13/25] add reno and HLS docs

---
 .../operators/symplectic/sparse_pauli_op.py   |  3 ++-
 .../passes/synthesis/hls_plugins.py           | 13 ++++++++++
 .../pauli-evo-plugins-612850146c3f7d49.yaml   | 24 +++++++++++++++++++
 3 files changed, 39 insertions(+), 1 deletion(-)
 create mode 100644 releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml

diff --git a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index 18f3af600e80..31faaeb5ebaf 100644
--- a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
+++ b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
@@ -939,7 +939,8 @@ def to_sparse_list(self):
         """Convert to a sparse Pauli list format with elements (pauli, qubits, coefficient)."""
         pauli_labels = self.paulis.to_labels()
         sparse_list = [
-            (*sparsify_label(label), coeff) for label, coeff in zip(pauli_labels, self.coeffs)
+            (*sparsify_label(label), np.real_if_close(coeff))
+            for label, coeff in zip(pauli_labels, self.coeffs)
         ]
         return sparse_list
 
diff --git a/qiskit/transpiler/passes/synthesis/hls_plugins.py b/qiskit/transpiler/passes/synthesis/hls_plugins.py
index 3f563274f871..3e26809a55f7 100644
--- a/qiskit/transpiler/passes/synthesis/hls_plugins.py
+++ b/qiskit/transpiler/passes/synthesis/hls_plugins.py
@@ -92,6 +92,19 @@
    PMHSynthesisLinearFunction
    DefaultSynthesisLinearFunction
 
+   
+Pauli Evolution Synthesis
+'''''''''''''''''''''''''
+
+.. list-table:: Plugins for :class:`.PauliEvolutionGate` (key = ``"PauliEvolution"``)
+    :header-rows: 1
+
+    * - Plugin name
+      - Plugin class
+      - Description
+    * - ``"default"``
+      - :class:`.PauliEvolutionSynthesisDefault`
+      - use a diagonalizing Clifford per Pauli term
 
 Permutation Synthesis
 '''''''''''''''''''''
diff --git a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
new file mode 100644
index 000000000000..cd56af8ffc7f
--- /dev/null
+++ b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
@@ -0,0 +1,24 @@
+---
+features_quantum_info:
+  - |
+    Added :meth:`.SparsePauliOperator.to_sparse_list` to convert an operator into
+    a sparse list format. This works analogous to :meth:`.SparsePauliOperator.from_sparse_list`.
+    For example::
+
+      from qiskit.quantum_info import SparsePauliOp
+
+      op = SparsePauliOp(["XIII", "IZZI"], coeffs=[1, 2])
+      sparse = op.to_sparse_list()  # [("X", [3], 1), ("ZZ", [1, 2], 2)]
+
+      other = SparsePauliOp.from_sparse_list(sparse, op.num_qubits)
+      print(other == op)  # True
+
+features_synthesis:
+  - |
+    Added :meth:`.ProductFormula.expand` which allows to view the expansion of a product formula
+    in a sparse Pauli format.
+  - |
+    Added the plugin structure for the :class:`.PauliEvolutionGate`. The default plugin,
+    :class:`.PauliEvolutionSynthesisDefault`, constructs circuit as before (but faster as it
+    internally uses Rust).
+

From 5c95b0dc090e52f3fbf600c45905b281813ae566 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 15 Oct 2024 13:49:32 +0200
Subject: [PATCH 14/25] fix unreachable

---
 crates/accelerate/src/circuit_library/pauli_evolution.rs    | 3 +--
 crates/circuit/src/operations.rs                            | 2 +-
 qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py | 3 +--
 3 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 15c725de218e..e2630c08e2d8 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -155,9 +155,8 @@ fn multi_qubit_evolution(
     let inverse_basis_change = basis_change.clone().map(|(gate, _, qubit)| match gate {
         StandardGate::HGate => (gate, smallvec![], qubit),
         StandardGate::RXGate => (gate, smallvec![Param::Float(-PI2)], qubit),
-        _ => unreachable!(),
+        _ => unreachable!("Invalid basis-changing Clifford."),
     });
-    println!("Fountain: {}", do_fountain);
 
     // get the CX propagation up to the first qubit, and down
     let (chain_up, chain_down) = match do_fountain {
diff --git a/crates/circuit/src/operations.rs b/crates/circuit/src/operations.rs
index 38f36118cdc2..e43fa548d245 100644
--- a/crates/circuit/src/operations.rs
+++ b/crates/circuit/src/operations.rs
@@ -2282,7 +2282,7 @@ pub fn multiply_param(param: &Param, mult: f64, py: Python) -> Param {
                 .call_method1(py, intern!(py, "__rmul__"), (mult,))
                 .expect("Multiplication of Parameter expression by float failed."),
         ),
-        Param::Obj(_) => unreachable!(),
+        Param::Obj(_) => unreachable!("Unsupported multiplication of a Param::Obj."),
     }
 }
 
diff --git a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index 31faaeb5ebaf..18f3af600e80 100644
--- a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
+++ b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
@@ -939,8 +939,7 @@ def to_sparse_list(self):
         """Convert to a sparse Pauli list format with elements (pauli, qubits, coefficient)."""
         pauli_labels = self.paulis.to_labels()
         sparse_list = [
-            (*sparsify_label(label), np.real_if_close(coeff))
-            for label, coeff in zip(pauli_labels, self.coeffs)
+            (*sparsify_label(label), coeff) for label, coeff in zip(pauli_labels, self.coeffs)
         ]
         return sparse_list
 

From 2eb66f113b3e50f4efa90b7f75a74aed912f69ed Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 15 Oct 2024 16:26:45 +0200
Subject: [PATCH 15/25] fix QPY test & pauli feature map

---
 .../src/circuit_library/pauli_feature_map.rs     |  2 +-
 qiskit/synthesis/evolution/product_formula.py    | 16 ++++++++--------
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_feature_map.rs b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
index b6b56fd55f29..bb9f8c25eb24 100644
--- a/crates/accelerate/src/circuit_library/pauli_feature_map.rs
+++ b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
@@ -173,7 +173,7 @@ fn _get_evolution_layer<'a>(
             // from using CircuitData::from_standard_gates.
             let evo = pauli_evolution::pauli_evolution(
                 pauli,
-                indices.clone(),
+                indices.into_iter().rev().collect(),
                 multiply_param(&angle, alpha, py),
                 true,
                 false,
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index aadea5d6de00..2ac4f151feda 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -92,15 +92,13 @@ def __init__(
 
         # user-provided atomic evolution, stored for serialization
         self._atomic_evolution = atomic_evolution
-        self._wrap = wrap
 
-        if cx_structure == "chain":
-            self._cx_fountain = False
-        elif cx_structure == "fountain":
-            self._cx_fountain = True
-        else:
+        if cx_structure not in ["chain", "fountain"]:
             raise ValueError(f"Unsupported CX structure: {cx_structure}")
 
+        self._cx_structure = cx_structure
+        self._wrap = wrap
+
         # if atomic evolution is not provided, set a default
         if atomic_evolution is None:
             self.atomic_evolution = None
@@ -152,8 +150,9 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
             circuit = self._custom_evolution(num_qubits, pauli_rotations)
         else:
             # this is the fast path, where the whole evolution is constructed Rust-side
+            cx_fountain = self._cx_structure == "fountain"
             data = py_pauli_evolution(
-                num_qubits, pauli_rotations, self.insert_barriers, self._cx_fountain
+                num_qubits, pauli_rotations, self.insert_barriers, cx_fountain
             )
             circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
 
@@ -195,6 +194,7 @@ def _custom_evolution(self, num_qubits, pauli_rotations):
         of individual Paulis needs to be wrapped in gates.
         """
         circuit = QuantumCircuit(num_qubits)
+        cx_fountain = self._cx_structure == "fountain"
 
         num_paulis = len(pauli_rotations)
         for i, pauli_rotation in enumerate(pauli_rotations):
@@ -214,7 +214,7 @@ def _custom_evolution(self, num_qubits, pauli_rotations):
                     len(qubits),
                     [local_pauli],
                     False,
-                    self._cx_fountain,
+                    cx_fountain,
                 )
                 evo = QuantumCircuit._from_circuit_data(data)
 

From 87496d29c23c8d8151493349866dfaa8db17ad59 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 15 Oct 2024 16:45:22 +0200
Subject: [PATCH 16/25] use SX as basis tranfo

---
 .../src/circuit_library/pauli_evolution.rs    | 34 +++++++----------
 .../data_preparation/pauli_feature_map.py     |  5 ++-
 qiskit/circuit/library/pauli_evolution.py     | 19 +++++-----
 .../pauli-evo-plugins-612850146c3f7d49.yaml   |  8 ++++
 .../circuit/library/test_evolution_gate.py    |  4 +-
 .../circuit/library/test_evolved_op_ansatz.py |  3 +-
 .../circuit/library/test_pauli_feature_map.py | 38 ++++++++++---------
 7 files changed, 57 insertions(+), 54 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index e2630c08e2d8..4bf4432d9576 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -17,12 +17,10 @@ use qiskit_circuit::operations::{multiply_param, radd_param, Param, StandardGate
 use qiskit_circuit::packed_instruction::PackedOperation;
 use qiskit_circuit::{Clbit, Qubit};
 use smallvec::{smallvec, SmallVec};
-use std::f64::consts::PI;
 
 use crate::circuit_library::utils;
 
-// custom math and types for a more readable code
-const PI2: f64 = PI / 2.;
+// custom types for a more readable code
 type StandardInstruction = (StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>);
 type Instruction = (
     PackedOperation,
@@ -136,25 +134,21 @@ fn multi_qubit_evolution(
         .zip(indices.into_iter().map(Qubit))
         .collect();
 
-    // get the basis change: x -> HGate, y -> RXGate(pi/2), z -> nothing
+    // get the basis change: x -> HGate, y -> SXdgGate, z -> nothing
     let basis_change = active_paulis
         .clone()
         .into_iter()
         .filter(|(p, _)| *p != 'z')
         .map(|(p, q)| match p {
             'x' => (StandardGate::HGate, smallvec![], smallvec![q]),
-            'y' => (
-                StandardGate::RXGate,
-                smallvec![Param::Float(PI2)],
-                smallvec![q],
-            ),
+            'y' => (StandardGate::SXdgGate, smallvec![], smallvec![q]),
             _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
         });
 
     // get the inverse basis change
     let inverse_basis_change = basis_change.clone().map(|(gate, _, qubit)| match gate {
         StandardGate::HGate => (gate, smallvec![], qubit),
-        StandardGate::RXGate => (gate, smallvec![Param::Float(-PI2)], qubit),
+        StandardGate::SXdgGate => (StandardGate::SXGate, smallvec![], qubit),
         _ => unreachable!("Invalid basis-changing Clifford."),
     });
 
@@ -182,8 +176,6 @@ fn multi_qubit_evolution(
         smallvec![first_qubit],
     ));
 
-    println!("first qubit {:?}", first_qubit);
-
     // and finally chain everything together
     basis_change
         .chain(chain_down)
@@ -198,15 +190,15 @@ fn multi_qubit_evolution(
 /// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
 /// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
 /// Pauli string XIYZ we have the circuit
-///                     ┌───┐┌───────┐┌───┐
-/// 0: ─────────────────┤ X ├┤ Rz(2) ├┤ X ├──────────────────
-///    ┌──────────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌───────────┐
-/// 1: ┤ Rx(pi/2) ├┤ X ├──■─────────────■──┤ X ├┤ Rx(-pi/2) ├
-///    └──────────┘└─┬─┘                   └─┬─┘└───────────┘
-/// 2: ──────────────┼───────────────────────┼───────────────
-///     ┌───┐        │                       │  ┌───┐
-/// 3: ─┤ H ├────────■───────────────────────■──┤ H ├────────
-///     └───┘                                   └───┘
+///                 ┌───┐┌───────┐┌───┐
+/// 0: ─────────────┤ X ├┤ Rz(2) ├┤ X ├───────────
+///    ┌──────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌────┐
+/// 1: ┤ √Xdg ├┤ X ├──■─────────────■──┤ X ├┤ √X ├
+///    └──────┘└─┬─┘                   └─┬─┘└────┘
+/// 2: ──────────┼───────────────────────┼────────
+///     ┌───┐    │                       │  ┌───┐
+/// 3: ─┤ H ├────■───────────────────────■──┤ H ├─
+///     └───┘                               └───┘
 ///
 /// Args:
 ///     num_qubits: The number of qubits in the Hamiltonian.
diff --git a/qiskit/circuit/library/data_preparation/pauli_feature_map.py b/qiskit/circuit/library/data_preparation/pauli_feature_map.py
index 75107ee5546f..87d9b709bc92 100644
--- a/qiskit/circuit/library/data_preparation/pauli_feature_map.py
+++ b/qiskit/circuit/library/data_preparation/pauli_feature_map.py
@@ -585,7 +585,10 @@ def basis_change(circuit, inverse=False):
                 if pauli == "X":
                     circuit.h(i)
                 elif pauli == "Y":
-                    circuit.rx(-np.pi / 2 if inverse else np.pi / 2, i)
+                    if inverse:
+                        circuit.sx(i)
+                    else:
+                        circuit.sxdg(i)
 
         def cx_chain(circuit, inverse=False):
             num_cx = len(indices) - 1
diff --git a/qiskit/circuit/library/pauli_evolution.py b/qiskit/circuit/library/pauli_evolution.py
index 23ccd065aad7..f7ab0393f32c 100644
--- a/qiskit/circuit/library/pauli_evolution.py
+++ b/qiskit/circuit/library/pauli_evolution.py
@@ -14,10 +14,11 @@
 
 from __future__ import annotations
 
-from typing import Union, Optional, TYPE_CHECKING
+from typing import TYPE_CHECKING
 import numpy as np
 
 from qiskit.circuit.gate import Gate
+from qiskit.circuit.quantumcircuit import ParameterValueType
 from qiskit.circuit.parameterexpression import ParameterExpression
 from qiskit.quantum_info import Pauli, SparsePauliOp
 
@@ -85,10 +86,10 @@ class PauliEvolutionGate(Gate):
 
     def __init__(
         self,
-        operator,
-        time: Union[int, float, ParameterExpression] = 1.0,
-        label: Optional[str] = None,
-        synthesis: Optional[EvolutionSynthesis] = None,
+        operator: Pauli | SparsePauliOp | list[Pauli | SparsePauliOp],
+        time: ParameterValueType = 1.0,
+        label: str | None = None,
+        synthesis: EvolutionSynthesis | None = None,
     ) -> None:
         """
         Args:
@@ -123,7 +124,7 @@ class docstring for an example.
         self.synthesis = synthesis
 
     @property
-    def time(self) -> Union[float, ParameterExpression]:
+    def time(self) -> ParameterValueType:
         """Return the evolution time as stored in the gate parameters.
 
         Returns:
@@ -132,7 +133,7 @@ def time(self) -> Union[float, ParameterExpression]:
         return self.params[0]
 
     @time.setter
-    def time(self, time: Union[float, ParameterExpression]) -> None:
+    def time(self, time: ParameterValueType) -> None:
         """Set the evolution time.
 
         Args:
@@ -144,9 +145,7 @@ def _define(self):
         """Unroll, where the default synthesis is matrix based."""
         self.definition = self.synthesis.synthesize(self)
 
-    def validate_parameter(
-        self, parameter: Union[int, float, ParameterExpression]
-    ) -> Union[float, ParameterExpression]:
+    def validate_parameter(self, parameter: ParameterValueType) -> ParameterValueType:
         """Gate parameters should be int, float, or ParameterExpression"""
         if isinstance(parameter, int):
             parameter = float(parameter)
diff --git a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
index cd56af8ffc7f..4a01f7c75888 100644
--- a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
+++ b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
@@ -22,3 +22,11 @@ features_synthesis:
     :class:`.PauliEvolutionSynthesisDefault`, constructs circuit as before (but faster as it
     internally uses Rust).
 
+upgrade:
+  - | 
+    The following classes now use the :math:`\sqrt{X}` operation to diagonalize the Pauli-Y
+    operator: :class:`.PauliEvolutionGate`, :class:`.EvolvedOperatorAnsatz`, 
+    :class:`.PauliFeatureMap`. Previously, these classes used either :math:`H S` or
+    :math:`R_X(-\pi/2)` as basis transformation. Using the :math:`\sqrt{X}` operation, 
+    represented by the :class:`.SXGate` is more efficient as it uses only a single gate
+    implemented as singleton.
diff --git a/test/python/circuit/library/test_evolution_gate.py b/test/python/circuit/library/test_evolution_gate.py
index 3b94003af35b..9bcc6565a6e9 100644
--- a/test/python/circuit/library/test_evolution_gate.py
+++ b/test/python/circuit/library/test_evolution_gate.py
@@ -376,8 +376,8 @@ def diagonalizing_clifford(pauli: Pauli) -> QuantumCircuit:
     cliff = QuantumCircuit(pauli.num_qubits)
     for i, pauli_i in enumerate(reversed(pauli.to_label())):
         if pauli_i == "Y":
-            cliff.sdg(i)
-        if pauli_i in ["X", "Y"]:
+            cliff.sxdg(i)
+        elif pauli_i == "X":
             cliff.h(i)
 
     return cliff
diff --git a/test/python/circuit/library/test_evolved_op_ansatz.py b/test/python/circuit/library/test_evolved_op_ansatz.py
index 8ea66b7012fe..37528e00857e 100644
--- a/test/python/circuit/library/test_evolved_op_ansatz.py
+++ b/test/python/circuit/library/test_evolved_op_ansatz.py
@@ -119,8 +119,7 @@ def evolve(pauli_string, time):
         if pauli == "x":
             forward.h(i)
         elif pauli == "y":
-            forward.sdg(i)
-            forward.h(i)
+            forward.sxdg(i)
 
     for i in range(1, num_qubits):
         forward.cx(num_qubits - i, num_qubits - i - 1)
diff --git a/test/python/circuit/library/test_pauli_feature_map.py b/test/python/circuit/library/test_pauli_feature_map.py
index 5d61944d4032..26961a7d8d06 100644
--- a/test/python/circuit/library/test_pauli_feature_map.py
+++ b/test/python/circuit/library/test_pauli_feature_map.py
@@ -70,24 +70,26 @@ def test_pauli_evolution(self):
             self.assertTrue(Operator(pauli).equiv(evo))
 
         with self.subTest(pauli_string="XYZ"):
-            # q_0: ─────────────■────────────────────────■──────────────
-            #      ┌─────────┐┌─┴─┐                    ┌─┴─┐┌──────────┐
-            # q_1: ┤ Rx(π/2) ├┤ X ├──■──────────────■──┤ X ├┤ Rx(-π/2) ├
-            #      └──┬───┬──┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└──────────┘
-            # q_2: ───┤ H ├────────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├────────────
-            #         └───┘        └───┘└────────┘└───┘└───┘
+            # q_0: ──────────■────────────────────────■───────
+            #      ┌──────┐┌─┴─┐                    ┌─┴─┐┌────┐
+            # q_1: ┤ √Xdg ├┤ X ├──■──────────────■──┤ X ├┤ √X ├
+            #      └─┬───┬┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└────┘
+            # q_2: ──┤ H ├──────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├──────
+            #        └───┘      └───┘└────────┘└───┘└───┘
             evo = QuantumCircuit(3)
             # X on the most-significant, bottom qubit, Z on the top
             evo.h(2)
-            evo.rx(np.pi / 2, 1)
+            evo.sxdg(1)
             evo.cx(0, 1)
             evo.cx(1, 2)
             evo.p(2 * time, 2)
             evo.cx(1, 2)
             evo.cx(0, 1)
-            evo.rx(-np.pi / 2, 1)
+            evo.sx(1)
             evo.h(2)
 
+            print(evo.draw())
+
             pauli = encoding.pauli_evolution("XYZ", time)
             self.assertTrue(Operator(pauli).equiv(evo))
 
@@ -347,23 +349,23 @@ def test_pauli_xyz(self):
 
         params = encoding.parameters
 
-        # q_0: ─────────────■────────────────────────■──────────────
-        #      ┌─────────┐┌─┴─┐                    ┌─┴─┐┌──────────┐
-        # q_1: ┤ Rx(π/2) ├┤ X ├──■──────────────■──┤ X ├┤ Rx(-π/2) ├
-        #      └──┬───┬──┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└──────────┘
-        # q_2: ───┤ H ├────────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├────────────
-        #         └───┘        └───┘└────────┘└───┘└───┘
+        # q_0: ──────────■────────────────────────■────────
+        #      ┌──────┐┌─┴─┐                    ┌─┴─┐┌────┐
+        # q_1: ┤ √Xdg ├┤ X ├──■──────────────■──┤ X ├┤ √X ├
+        #      └─┬───┬┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└────┘
+        # q_2: ──┤ H ├──────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├──────
+        #        └───┘      └───┘└────────┘└───┘└───┘
         # X on the most-significant, bottom qubit, Z on the top
         ref = QuantumCircuit(3)
         ref.h(range(3))
         ref.h(2)
-        ref.rx(np.pi / 2, 1)
+        ref.sxdg(1)
         ref.cx(0, 1)
         ref.cx(1, 2)
         ref.p(2 * np.prod([np.pi - p for p in params]), 2)
         ref.cx(1, 2)
         ref.cx(0, 1)
-        ref.rx(-np.pi / 2, 1)
+        ref.sx(1)
         ref.h(2)
 
         self.assertEqual(ref, encoding)
@@ -483,13 +485,13 @@ def test_dict_entanglement(self):
         ref.cx(1, 2)
         ref.h([1, 2])
 
-        ref.rx(np.pi / 2, range(3))
+        ref.sxdg(range(3))
         ref.cx(0, 1)
         ref.cx(1, 2)
         ref.p(2 * np.prod([np.pi - xi for xi in x]), 2)
         ref.cx(1, 2)
         ref.cx(0, 1)
-        ref.rx(-np.pi / 2, range(3))
+        ref.sx(range(3))
 
         self.assertEqual(ref, circuit)
 

From 546cc5fc946c8571f8da49129ed76987bfb36ebd Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Wed, 23 Oct 2024 16:55:10 +0200
Subject: [PATCH 17/25] slight performance tweak

move the x2 multiplication of the rotation angle to Python, where we can pull it together with other floats and minimize the float * Param multiplications
---
 .../src/circuit_library/pauli_evolution.rs          |  2 +-
 qiskit/synthesis/evolution/qdrift.py                |  2 +-
 qiskit/synthesis/evolution/suzuki_trotter.py        | 13 +++++++++----
 3 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 4bf4432d9576..243b70e6edee 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -240,7 +240,7 @@ pub fn py_pauli_evolution(
         }
 
         paulis.push(pauli);
-        times.push(multiply_param(&time, 2.0, py));
+        times.push(time); // note we do not multiply by 2 here, this is done Python side!
         indices.push(
             tuple
                 .get_item(1)?
diff --git a/qiskit/synthesis/evolution/qdrift.py b/qiskit/synthesis/evolution/qdrift.py
index e214d4c113dc..64ae5fba637b 100644
--- a/qiskit/synthesis/evolution/qdrift.py
+++ b/qiskit/synthesis/evolution/qdrift.py
@@ -121,7 +121,7 @@ def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], f
             np.array(paulis, dtype=object), size=(num_gates,), p=weights / lambd
         )
 
-        rescaled_time = lambd * time / num_gates
+        rescaled_time = 2 * lambd / num_gates * time
         sampled_paulis = [
             (pauli[0], pauli[1], np.real(np.sign(pauli[2])) * rescaled_time) for pauli in sampled
         ]
diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index 0201c5367dbe..6bbc27ad3b60 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -121,11 +121,16 @@ def expand(
         """Expand the Hamiltonian into a Suzuki-Trotter sequence of sparse gates.
 
         For example, the Hamiltonian ``H = IX + ZZ`` for an evolution time ``t`` and
-        1 repetition for an order 2 formula would get decomposed into
+        1 repetition for an order 2 formula would get decomposed into a list of 3-tuples
+        containing ``(pauli, indices, rz_rotation_angle)``, that is:
 
         .. code-block:: text
 
-            ("X", [0], t/2), ("ZZ", [0, 1], t), ("X", [0], t/2)
+            ("X", [0], t), ("ZZ", [0, 1], 2t), ("X", [0], 2)
+
+        Note that the rotation angle contains a factor of 2, such that that evolution
+        of a Pauli :math:`P` over time :math:`t`, which is :math:`e^{itP}`, is represented
+        by ``(P, indices, 2 * t)``.
 
         For ``N`` repetitions, this sequence would be repeated ``N`` times and the coefficients
         divided by ``N``.
@@ -142,12 +147,12 @@ def expand(
         # construct the evolution circuit
         if isinstance(operators, list):  # already sorted into commuting bits
             non_commuting = [
-                (time / self.reps * operator).to_sparse_list() for operator in operators
+                (2 / self.reps * time * operator).to_sparse_list() for operator in operators
             ]
         else:
             # Assume no commutativity here. If we were to group commuting Paulis,
             # here would be the location to do so.
-            non_commuting = [[op] for op in (time / self.reps * operators).to_sparse_list()]
+            non_commuting = [[op] for op in (2 / self.reps * time * operators).to_sparse_list()]
 
         # normalize coefficients, i.e. ensure they are float or ParameterExpression
         non_commuting = self._normalize_coefficients(non_commuting)

From cf7d8b09681c400c2213c6c35722f85a7f7c63b4 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 28 Oct 2024 13:23:00 +0100
Subject: [PATCH 18/25] implement review comments

---
 crates/accelerate/src/circuit_library/mod.rs  |  1 -
 .../src/circuit_library/pauli_evolution.rs    | 95 ++++++++++++-------
 .../accelerate/src/circuit_library/utils.rs   | 58 -----------
 qiskit/synthesis/evolution/product_formula.py |  8 +-
 .../pauli-evo-plugins-612850146c3f7d49.yaml   |  2 +-
 5 files changed, 65 insertions(+), 99 deletions(-)
 delete mode 100644 crates/accelerate/src/circuit_library/utils.rs

diff --git a/crates/accelerate/src/circuit_library/mod.rs b/crates/accelerate/src/circuit_library/mod.rs
index b943881e4a14..3897167f9489 100644
--- a/crates/accelerate/src/circuit_library/mod.rs
+++ b/crates/accelerate/src/circuit_library/mod.rs
@@ -16,7 +16,6 @@ mod entanglement;
 mod pauli_evolution;
 mod pauli_feature_map;
 mod quantum_volume;
-mod utils;
 
 pub fn circuit_library(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(pauli_evolution::py_pauli_evolution))?;
diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 243b70e6edee..9b0a43cc04b1 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -13,13 +13,11 @@
 use pyo3::prelude::*;
 use pyo3::types::{PyList, PyString, PyTuple};
 use qiskit_circuit::circuit_data::CircuitData;
-use qiskit_circuit::operations::{multiply_param, radd_param, Param, StandardGate};
+use qiskit_circuit::operations::{multiply_param, radd_param, Param, PyInstruction, StandardGate};
 use qiskit_circuit::packed_instruction::PackedOperation;
-use qiskit_circuit::{Clbit, Qubit};
+use qiskit_circuit::{imports, Clbit, Qubit};
 use smallvec::{smallvec, SmallVec};
 
-use crate::circuit_library::utils;
-
 // custom types for a more readable code
 type StandardInstruction = (StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>);
 type Instruction = (
@@ -115,7 +113,7 @@ fn two_qubit_evolution<'a>(
         "zx" => Box::new(std::iter::once((StandardGate::RZXGate, param, qubits))),
         "yy" => Box::new(std::iter::once((StandardGate::RYYGate, param, qubits))),
         "zz" => Box::new(std::iter::once((StandardGate::RZZGate, param, qubits))),
-        // note that the CX modes (do_fountain=true/false) give the same circuit for a 2-qubit
+        // Note: the CX modes (do_fountain=true/false) give the same circuit for a 2-qubit
         // Pauli, so we just set it to false here
         _ => Box::new(multi_qubit_evolution(pauli, indices, time, false, false)),
     }
@@ -135,22 +133,25 @@ fn multi_qubit_evolution(
         .collect();
 
     // get the basis change: x -> HGate, y -> SXdgGate, z -> nothing
-    let basis_change = active_paulis
-        .clone()
-        .into_iter()
+    let basis_change: Vec<StandardInstruction> = active_paulis
+        .iter()
         .filter(|(p, _)| *p != 'z')
         .map(|(p, q)| match p {
-            'x' => (StandardGate::HGate, smallvec![], smallvec![q]),
-            'y' => (StandardGate::SXdgGate, smallvec![], smallvec![q]),
+            'x' => (StandardGate::HGate, smallvec![], smallvec![q.clone()]),
+            'y' => (StandardGate::SXdgGate, smallvec![], smallvec![q.clone()]),
             _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
-        });
+        })
+        .collect();
 
     // get the inverse basis change
-    let inverse_basis_change = basis_change.clone().map(|(gate, _, qubit)| match gate {
-        StandardGate::HGate => (gate, smallvec![], qubit),
-        StandardGate::SXdgGate => (StandardGate::SXGate, smallvec![], qubit),
-        _ => unreachable!("Invalid basis-changing Clifford."),
-    });
+    let inverse_basis_change: Vec<StandardInstruction> = basis_change
+        .iter()
+        .map(|(gate, _, qubit)| match gate {
+            StandardGate::HGate => (StandardGate::HGate, smallvec![], qubit.clone()),
+            StandardGate::SXdgGate => (StandardGate::SXGate, smallvec![], qubit.clone()),
+            _ => unreachable!("Invalid basis-changing Clifford."),
+        })
+        .collect();
 
     // get the CX propagation up to the first qubit, and down
     let (chain_up, chain_down) = match do_fountain {
@@ -178,6 +179,7 @@ fn multi_qubit_evolution(
 
     // and finally chain everything together
     basis_change
+        .into_iter()
         .chain(chain_down)
         .chain(z_rotation)
         .chain(chain_up)
@@ -215,19 +217,20 @@ fn multi_qubit_evolution(
 /// Returns:
 ///     Circuit data for to implement the evolution.
 #[pyfunction]
-#[pyo3(signature = (num_qubits, sparse_paulis, insert_barriers=false, do_fountain=false))]
+#[pyo3(name = "pauli_evolution", signature = (num_qubits, sparse_paulis, insert_barriers=false, do_fountain=false))]
 pub fn py_pauli_evolution(
-    py: Python,
     num_qubits: i64,
     sparse_paulis: &Bound<PyList>,
     insert_barriers: bool,
     do_fountain: bool,
 ) -> PyResult<CircuitData> {
+    let py = sparse_paulis.py();
     let num_paulis = sparse_paulis.len();
     let mut paulis: Vec<String> = Vec::with_capacity(num_paulis);
     let mut indices: Vec<Vec<u32>> = Vec::with_capacity(num_paulis);
     let mut times: Vec<Param> = Vec::with_capacity(num_paulis);
     let mut global_phase = Param::Float(0.0);
+    let mut modified_phase = false; // keep track of whether we modified the phase
 
     for el in sparse_paulis.iter() {
         let tuple = el.downcast::<PyTuple>()?;
@@ -235,22 +238,18 @@ pub fn py_pauli_evolution(
         let time = Param::extract_no_coerce(&tuple.get_item(2)?)?;
 
         if pauli.as_str().chars().all(|p| p == 'i') {
-            global_phase = radd_param(global_phase, time, py); // apply factor -1 at the end
+            global_phase = radd_param(global_phase, time, py);
+            modified_phase = true;
             continue;
         }
 
         paulis.push(pauli);
         times.push(time); // note we do not multiply by 2 here, this is done Python side!
-        indices.push(
-            tuple
-                .get_item(1)?
-                .downcast::<PyList>()?
-                .iter()
-                .map(|index| index.extract::<u32>())
-                .collect::<PyResult<_>>()?,
-        );
+        indices.push(tuple.get_item(1)?.extract::<Vec<u32>>()?)
     }
 
+    let barrier = get_barrier(py, num_qubits as u32);
+
     let evos = paulis.iter().enumerate().zip(indices).zip(times).flat_map(
         |(((i, pauli), qubits), time)| {
             let as_packed = pauli_evolution(pauli, qubits, time, false, do_fountain).map(
@@ -263,16 +262,23 @@ pub fn py_pauli_evolution(
                     ))
                 },
             );
-            as_packed.chain(utils::maybe_barrier(
-                py,
-                num_qubits as u32,
-                insert_barriers && i < (num_paulis - 1), // do not add barrier on final block
-            ))
+
+            // this creates an iterator containing a barrier only if required, otherwise it is empty
+            let maybe_barrier = (insert_barriers && i < (num_paulis - 1))
+                .then_some(Ok(barrier.clone()))
+                .into_iter();
+            as_packed.chain(maybe_barrier)
         },
     );
 
-    // apply factor -1 for global phase
-    global_phase = multiply_param(&global_phase, -1.0, py);
+    // When handling all-identity Paulis above, we added the time as global phase.
+    // However, the all-identity Paulis should add a negative phase, as they implement
+    // exp(-i t I). We apply the negative sign here, to only do a single (-1) multiplication,
+    // instead of doing it every time we find an all-identity Pauli.
+    if modified_phase {
+        global_phase = multiply_param(&global_phase, -1.0, py);
+    }
+
     CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)
 }
 
@@ -301,3 +307,24 @@ fn cx_fountain(
         )
     }))
 }
+
+fn get_barrier(py: Python, num_qubits: u32) -> Instruction {
+    let barrier_cls = imports::BARRIER.get_bound(py);
+    let barrier = barrier_cls
+        .call1((num_qubits,))
+        .expect("Could not create Barrier Python-side");
+    let barrier_inst = PyInstruction {
+        qubits: num_qubits,
+        clbits: 0,
+        params: 0,
+        op_name: "barrier".to_string(),
+        control_flow: false,
+        instruction: barrier.into(),
+    };
+    (
+        barrier_inst.into(),
+        smallvec![],
+        (0..num_qubits).map(Qubit).collect(),
+        vec![],
+    )
+}
diff --git a/crates/accelerate/src/circuit_library/utils.rs b/crates/accelerate/src/circuit_library/utils.rs
deleted file mode 100644
index 4df8d64865f0..000000000000
--- a/crates/accelerate/src/circuit_library/utils.rs
+++ /dev/null
@@ -1,58 +0,0 @@
-// This code is part of Qiskit.
-//
-// (C) Copyright IBM 2024
-//
-// 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.
-
-use pyo3::prelude::*;
-use qiskit_circuit::{
-    imports,
-    operations::{Param, PyInstruction},
-    packed_instruction::PackedOperation,
-    Clbit, Qubit,
-};
-use smallvec::{smallvec, SmallVec};
-
-type Instruction = (
-    PackedOperation,
-    SmallVec<[Param; 3]>,
-    Vec<Qubit>,
-    Vec<Clbit>,
-);
-
-/// Get an iterator that returns a barrier or an empty element.
-pub fn maybe_barrier(
-    py: Python,
-    num_qubits: u32,
-    insert_barriers: bool,
-) -> Box<dyn Iterator<Item = PyResult<Instruction>>> {
-    // TODO could speed this up by only defining the barrier class once
-    if !insert_barriers {
-        Box::new(std::iter::empty())
-    } else {
-        let barrier_cls = imports::BARRIER.get_bound(py);
-        let barrier = barrier_cls
-            .call1((num_qubits,))
-            .expect("Could not create Barrier Python-side");
-        let barrier_inst = PyInstruction {
-            qubits: num_qubits,
-            clbits: 0,
-            params: 0,
-            op_name: "barrier".to_string(),
-            control_flow: false,
-            instruction: barrier.into(),
-        };
-        Box::new(std::iter::once(Ok((
-            barrier_inst.into(),
-            smallvec![],
-            (0..num_qubits).map(Qubit).collect(),
-            vec![] as Vec<Clbit>,
-        ))))
-    }
-}
diff --git a/qiskit/synthesis/evolution/product_formula.py b/qiskit/synthesis/evolution/product_formula.py
index 2ac4f151feda..4de838215ff1 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -22,7 +22,7 @@
 from qiskit.circuit.quantumcircuit import QuantumCircuit, ParameterValueType
 from qiskit.quantum_info import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
-from qiskit._accelerate.circuit_library import py_pauli_evolution
+from qiskit._accelerate.circuit_library import pauli_evolution
 
 from .evolution_synthesis import EvolutionSynthesis
 
@@ -151,9 +151,7 @@ def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
         else:
             # this is the fast path, where the whole evolution is constructed Rust-side
             cx_fountain = self._cx_structure == "fountain"
-            data = py_pauli_evolution(
-                num_qubits, pauli_rotations, self.insert_barriers, cx_fountain
-            )
+            data = pauli_evolution(num_qubits, pauli_rotations, self.insert_barriers, cx_fountain)
             circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
 
         return circuit
@@ -210,7 +208,7 @@ def _custom_evolution(self, num_qubits, pauli_rotations):
                 local_pauli = (pauli_string, list(range(len(qubits))), coeff)
 
                 # build the circuit Rust-side
-                data = py_pauli_evolution(
+                data = pauli_evolution(
                     len(qubits),
                     [local_pauli],
                     False,
diff --git a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
index 4a01f7c75888..67ccdfbae3b5 100644
--- a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
+++ b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
@@ -2,7 +2,7 @@
 features_quantum_info:
   - |
     Added :meth:`.SparsePauliOperator.to_sparse_list` to convert an operator into
-    a sparse list format. This works analogous to :meth:`.SparsePauliOperator.from_sparse_list`.
+    a sparse list format. This works inversely to :meth:`.SparsePauliOperator.from_sparse_list`.
     For example::
 
       from qiskit.quantum_info import SparsePauliOp

From 6f4fcbad4e283d6db25710524f2039981bd3c2a5 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 28 Oct 2024 13:34:34 +0100
Subject: [PATCH 19/25] do_fountain

---
 crates/accelerate/src/circuit_library/pauli_evolution.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 9b0a43cc04b1..5f0ec3d68e15 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -38,7 +38,7 @@ type Instruction = (
 ///         rotation gate and not be multiplied by a factor of 2 (that should be done before so
 ///         that this function can remain Rust-only).
 ///     phase_gate: If ``true``, use the ``PhaseGate`` instead of ``RZGate`` as single-qubit rotation.
-///     cx_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
+///     do_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
 ///         CX uses the top qubit as target. If ``false``, uses a "chain" shape, where CX in between
 ///         neighboring qubits are used.
 ///
@@ -210,7 +210,7 @@ fn multi_qubit_evolution(
 ///          I_q0 X_q1 Y_q2 Z_q3 and will use a RZ rotation angle of 0.4.
 ///     insert_barriers: If ``true``, insert a barrier in between the evolution of individual
 ///         Pauli terms.
-///     cx_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
+///     do_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
 ///         CX uses the top qubit as target. If ``false``, uses a "chain" shape, where CX in between
 ///         neighboring qubits are used.
 ///

From 6217d0ed0406569aeb1a3b1ddd3fffaf6b861b79 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 28 Oct 2024 13:56:53 +0100
Subject: [PATCH 20/25] clippy

---
 crates/accelerate/src/circuit_library/pauli_evolution.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 5f0ec3d68e15..169cece41810 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -137,8 +137,8 @@ fn multi_qubit_evolution(
         .iter()
         .filter(|(p, _)| *p != 'z')
         .map(|(p, q)| match p {
-            'x' => (StandardGate::HGate, smallvec![], smallvec![q.clone()]),
-            'y' => (StandardGate::SXdgGate, smallvec![], smallvec![q.clone()]),
+            'x' => (StandardGate::HGate, smallvec![], smallvec![*q]),
+            'y' => (StandardGate::SXdgGate, smallvec![], smallvec![*q]),
             _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
         })
         .collect();

From e802c45f85d411b5d9549f9d7266e9f628e720b0 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Thu, 31 Oct 2024 13:38:01 +0100
Subject: [PATCH 21/25] review comments

---
 qiskit/synthesis/evolution/suzuki_trotter.py          | 2 +-
 test/python/circuit/library/test_pauli_feature_map.py | 2 --
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/qiskit/synthesis/evolution/suzuki_trotter.py b/qiskit/synthesis/evolution/suzuki_trotter.py
index 6bbc27ad3b60..4dfa92579bcb 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -111,7 +111,7 @@ def __init__(
         if order > 1 and order % 2 == 1:
             raise ValueError(
                 "Suzuki product formulae are symmetric and therefore only defined "
-                "for even orders (or order==1)."
+                f"for when the order is 1 or even, not {order}."
             )
         super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
diff --git a/test/python/circuit/library/test_pauli_feature_map.py b/test/python/circuit/library/test_pauli_feature_map.py
index 26961a7d8d06..93f230f26ace 100644
--- a/test/python/circuit/library/test_pauli_feature_map.py
+++ b/test/python/circuit/library/test_pauli_feature_map.py
@@ -88,8 +88,6 @@ def test_pauli_evolution(self):
             evo.sx(1)
             evo.h(2)
 
-            print(evo.draw())
-
             pauli = encoding.pauli_evolution("XYZ", time)
             self.assertTrue(Operator(pauli).equiv(evo))
 

From 2b6f4ebc38a0eaa7bc25060985d84d8a06bedfeb Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Thu, 31 Oct 2024 15:36:12 +0100
Subject: [PATCH 22/25] fix the basis trafo!

---
 crates/accelerate/src/circuit_library/pauli_evolution.rs | 4 ++--
 test/python/circuit/library/test_evolution_gate.py       | 8 ++++++++
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index 169cece41810..d110310e9e45 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -138,7 +138,7 @@ fn multi_qubit_evolution(
         .filter(|(p, _)| *p != 'z')
         .map(|(p, q)| match p {
             'x' => (StandardGate::HGate, smallvec![], smallvec![*q]),
-            'y' => (StandardGate::SXdgGate, smallvec![], smallvec![*q]),
+            'y' => (StandardGate::SXGate, smallvec![], smallvec![*q]),
             _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
         })
         .collect();
@@ -148,7 +148,7 @@ fn multi_qubit_evolution(
         .iter()
         .map(|(gate, _, qubit)| match gate {
             StandardGate::HGate => (StandardGate::HGate, smallvec![], qubit.clone()),
-            StandardGate::SXdgGate => (StandardGate::SXGate, smallvec![], qubit.clone()),
+            StandardGate::SXGate => (StandardGate::SXdgGate, smallvec![], qubit.clone()),
             _ => unreachable!("Invalid basis-changing Clifford."),
         })
         .collect();
diff --git a/test/python/circuit/library/test_evolution_gate.py b/test/python/circuit/library/test_evolution_gate.py
index 9bcc6565a6e9..6b623b4f6cb2 100644
--- a/test/python/circuit/library/test_evolution_gate.py
+++ b/test/python/circuit/library/test_evolution_gate.py
@@ -60,6 +60,14 @@ def test_lie_trotter(self):
         decomposed = evo_gate.definition.decompose()
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)
 
+    def test_basis_change(self):
+        """Test the basis change is correctly implemented."""
+        op = I ^ X ^ Y ^ Z  # use a string for which we do not have a basis gate
+        time = 0.321
+        evolved = scipy.linalg.expm(-1j * time * op.to_matrix())
+        evo_gate = PauliEvolutionGate(op, time)
+        self.assertTrue(Operator(evo_gate).equiv(evolved))
+
     def test_rzx_order(self):
         """Test ZX and XZ is mapped onto the correct qubits."""
 

From 1af569bcc97a0a6bdc62269f74789d2bbd817c1f Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 4 Nov 2024 11:49:19 +0100
Subject: [PATCH 23/25] properly test the evolution basis changes

---
 .../data_preparation/pauli_feature_map.py     |  4 +-
 .../circuit/library/test_evolution_gate.py    | 40 +++++++++++++++----
 .../circuit/library/test_evolved_op_ansatz.py |  2 +-
 3 files changed, 36 insertions(+), 10 deletions(-)

diff --git a/qiskit/circuit/library/data_preparation/pauli_feature_map.py b/qiskit/circuit/library/data_preparation/pauli_feature_map.py
index c91bbea6b6a8..1bc14d169273 100644
--- a/qiskit/circuit/library/data_preparation/pauli_feature_map.py
+++ b/qiskit/circuit/library/data_preparation/pauli_feature_map.py
@@ -586,9 +586,9 @@ def basis_change(circuit, inverse=False):
                     circuit.h(i)
                 elif pauli == "Y":
                     if inverse:
-                        circuit.sx(i)
-                    else:
                         circuit.sxdg(i)
+                    else:
+                        circuit.sx(i)
 
         def cx_chain(circuit, inverse=False):
             num_cx = len(indices) - 1
diff --git a/test/python/circuit/library/test_evolution_gate.py b/test/python/circuit/library/test_evolution_gate.py
index 6b623b4f6cb2..7d1740b67f30 100644
--- a/test/python/circuit/library/test_evolution_gate.py
+++ b/test/python/circuit/library/test_evolution_gate.py
@@ -18,7 +18,7 @@
 from ddt import ddt, data, unpack
 
 from qiskit.circuit import QuantumCircuit, Parameter
-from qiskit.circuit.library import PauliEvolutionGate
+from qiskit.circuit.library import PauliEvolutionGate, HamiltonianGate
 from qiskit.synthesis import LieTrotter, SuzukiTrotter, MatrixExponential, QDrift
 from qiskit.converters import circuit_to_dag
 from qiskit.quantum_info import Operator, SparsePauliOp, Pauli, Statevector
@@ -39,6 +39,19 @@ def setUp(self):
         # fix random seed for reproducibility (used in QDrift)
         self.seed = 2
 
+    def assertSuzukiTrotterIsCorrect(self, gate):
+        """Assert the Suzuki Trotter evolution is correct."""
+        op = gate.operator
+        time = gate.time
+        synthesis = gate.synthesis
+
+        exact_suzuki = SuzukiTrotter(
+            reps=synthesis.reps, order=synthesis.order, atomic_evolution=exact_atomic_evolution
+        )
+        exact_gate = PauliEvolutionGate(op, time, synthesis=exact_suzuki)
+
+        self.assertTrue(Operator(gate).equiv(exact_gate))
+
     def test_matrix_decomposition(self):
         """Test the default decomposition."""
         op = (X ^ X ^ X) + (Y ^ Y ^ Y) + (Z ^ Z ^ Z)
@@ -58,15 +71,16 @@ def test_lie_trotter(self):
         reps = 4
         evo_gate = PauliEvolutionGate(op, time, synthesis=LieTrotter(reps=reps))
         decomposed = evo_gate.definition.decompose()
+
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)
+        self.assertSuzukiTrotterIsCorrect(evo_gate)
 
     def test_basis_change(self):
         """Test the basis change is correctly implemented."""
-        op = I ^ X ^ Y ^ Z  # use a string for which we do not have a basis gate
+        op = I ^ Y  # use a string for which we do not have a basis gate
         time = 0.321
-        evolved = scipy.linalg.expm(-1j * time * op.to_matrix())
         evo_gate = PauliEvolutionGate(op, time)
-        self.assertTrue(Operator(evo_gate).equiv(evolved))
+        self.assertSuzukiTrotterIsCorrect(evo_gate)
 
     def test_rzx_order(self):
         """Test ZX and XZ is mapped onto the correct qubits."""
@@ -112,6 +126,7 @@ def test_suzuki_trotter(self):
                 )
                 decomposed = evo_gate.definition.decompose()
                 self.assertEqual(decomposed.count_ops()["cx"], expected_cx)
+                self.assertSuzukiTrotterIsCorrect(evo_gate)
 
     def test_suzuki_trotter_manual(self):
         """Test the evolution circuit of Suzuki Trotter against a manually constructed circuit."""
@@ -139,6 +154,7 @@ def test_suzuki_trotter_manual(self):
             expected.rx(p_4 * time, 0)
 
         self.assertEqual(evo_gate.definition, expected)
+        self.assertSuzukiTrotterIsCorrect(evo_gate)
 
     @data(
         (X + Y, 0.5, 1, [(Pauli("X"), 0.5), (Pauli("X"), 0.5)]),
@@ -259,6 +275,7 @@ def test_cnot_chain_options(self, option):
             expected.cx(1, 0)
 
         self.assertEqual(expected, evo.definition)
+        self.assertSuzukiTrotterIsCorrect(evo)
 
     @data(
         Pauli("XI"),
@@ -283,6 +300,7 @@ def test_pauliop_coefficients_respected(self):
         circuit = evo.definition.decompose()
         rz_angle = circuit.data[0].operation.params[0]
         self.assertEqual(rz_angle, 10)
+        self.assertSuzukiTrotterIsCorrect(evo)
 
     def test_paulisumop_coefficients_respected(self):
         """Test that global ``PauliSumOp`` coefficients are being taken care of."""
@@ -294,6 +312,7 @@ def test_paulisumop_coefficients_respected(self):
             circuit.data[2].operation.params[0],  # Z
         ]
         self.assertListEqual(rz_angles, [20, 30, -10])
+        self.assertSuzukiTrotterIsCorrect(evo)
 
     def test_lie_trotter_two_qubit_correct_order(self):
         """Test that evolutions on two qubit operators are in the right order.
@@ -302,10 +321,9 @@ def test_lie_trotter_two_qubit_correct_order(self):
         """
         operator = I ^ Z ^ Z
         time = 0.5
-        exact = scipy.linalg.expm(-1j * time * operator.to_matrix())
         lie_trotter = PauliEvolutionGate(operator, time, synthesis=LieTrotter())
 
-        self.assertTrue(Operator(lie_trotter).equiv(exact))
+        self.assertSuzukiTrotterIsCorrect(lie_trotter)
 
     def test_complex_op_raises(self):
         """Test an operator with complex coefficient raises an error."""
@@ -379,12 +397,20 @@ def atomic_evolution(pauli, time):
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)
 
 
+def exact_atomic_evolution(circuit, pauli, time):
+    """An exact atomic evolution for Suzuki-Trotter.
+
+    Note that the Pauli has a x2 coefficient already, hence we evolve for time/2.
+    """
+    circuit.append(HamiltonianGate(pauli.to_matrix(), time / 2), circuit.qubits)
+
+
 def diagonalizing_clifford(pauli: Pauli) -> QuantumCircuit:
     """Get the clifford circuit to diagonalize the Pauli operator."""
     cliff = QuantumCircuit(pauli.num_qubits)
     for i, pauli_i in enumerate(reversed(pauli.to_label())):
         if pauli_i == "Y":
-            cliff.sxdg(i)
+            cliff.sx(i)
         elif pauli_i == "X":
             cliff.h(i)
 
diff --git a/test/python/circuit/library/test_evolved_op_ansatz.py b/test/python/circuit/library/test_evolved_op_ansatz.py
index 37528e00857e..2dbcec86a120 100644
--- a/test/python/circuit/library/test_evolved_op_ansatz.py
+++ b/test/python/circuit/library/test_evolved_op_ansatz.py
@@ -119,7 +119,7 @@ def evolve(pauli_string, time):
         if pauli == "x":
             forward.h(i)
         elif pauli == "y":
-            forward.sxdg(i)
+            forward.sx(i)
 
     for i in range(1, num_qubits):
         forward.cx(num_qubits - i, num_qubits - i - 1)

From c478cc162ec46a89f341b857af337cef37d38104 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 4 Nov 2024 13:50:48 +0100
Subject: [PATCH 24/25] fix more tests

---
 .../circuit/library/test_pauli_feature_map.py | 36 +++++++++----------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/test/python/circuit/library/test_pauli_feature_map.py b/test/python/circuit/library/test_pauli_feature_map.py
index 93f230f26ace..d12e89fc556f 100644
--- a/test/python/circuit/library/test_pauli_feature_map.py
+++ b/test/python/circuit/library/test_pauli_feature_map.py
@@ -70,22 +70,22 @@ def test_pauli_evolution(self):
             self.assertTrue(Operator(pauli).equiv(evo))
 
         with self.subTest(pauli_string="XYZ"):
-            # q_0: ──────────■────────────────────────■───────
-            #      ┌──────┐┌─┴─┐                    ┌─┴─┐┌────┐
-            # q_1: ┤ √Xdg ├┤ X ├──■──────────────■──┤ X ├┤ √X ├
-            #      └─┬───┬┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└────┘
-            # q_2: ──┤ H ├──────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├──────
-            #        └───┘      └───┘└────────┘└───┘└───┘
+            # q_0: ────────■────────────────────────■──────────
+            #      ┌────┐┌─┴─┐                    ┌─┴─┐┌──────┐
+            # q_1: ┤ √X ├┤ X ├──■──────────────■──┤ X ├┤ √Xdg ├
+            #      └┬───┤└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└──────┘
+            # q_2: ─┤ H ├─────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├────────
+            #       └───┘     └───┘└────────┘└───┘└───┘
             evo = QuantumCircuit(3)
             # X on the most-significant, bottom qubit, Z on the top
             evo.h(2)
-            evo.sxdg(1)
+            evo.sx(1)
             evo.cx(0, 1)
             evo.cx(1, 2)
             evo.p(2 * time, 2)
             evo.cx(1, 2)
             evo.cx(0, 1)
-            evo.sx(1)
+            evo.sxdg(1)
             evo.h(2)
 
             pauli = encoding.pauli_evolution("XYZ", time)
@@ -347,23 +347,23 @@ def test_pauli_xyz(self):
 
         params = encoding.parameters
 
-        # q_0: ──────────■────────────────────────■────────
-        #      ┌──────┐┌─┴─┐                    ┌─┴─┐┌────┐
-        # q_1: ┤ √Xdg ├┤ X ├──■──────────────■──┤ X ├┤ √X ├
-        #      └─┬───┬┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└────┘
-        # q_2: ──┤ H ├──────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├──────
-        #        └───┘      └───┘└────────┘└───┘└───┘
+        # q_0: ────────■────────────────────────■──────────
+        #      ┌────┐┌─┴─┐                    ┌─┴─┐┌──────┐
+        # q_1: ┤ √X ├┤ X ├──■──────────────■──┤ X ├┤ √Xdg ├
+        #      └┬───┤└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└──────┘
+        # q_2: ─┤ H ├─────┤ X ├┤ P(2.8) ├┤ X ├┤ H ├────────
+        #       └───┘     └───┘└────────┘└───┘└───┘
         # X on the most-significant, bottom qubit, Z on the top
         ref = QuantumCircuit(3)
         ref.h(range(3))
         ref.h(2)
-        ref.sxdg(1)
+        ref.sx(1)
         ref.cx(0, 1)
         ref.cx(1, 2)
         ref.p(2 * np.prod([np.pi - p for p in params]), 2)
         ref.cx(1, 2)
         ref.cx(0, 1)
-        ref.sx(1)
+        ref.sxdg(1)
         ref.h(2)
 
         self.assertEqual(ref, encoding)
@@ -483,13 +483,13 @@ def test_dict_entanglement(self):
         ref.cx(1, 2)
         ref.h([1, 2])
 
-        ref.sxdg(range(3))
+        ref.sx(range(3))
         ref.cx(0, 1)
         ref.cx(1, 2)
         ref.p(2 * np.prod([np.pi - xi for xi in x]), 2)
         ref.cx(1, 2)
         ref.cx(0, 1)
-        ref.sx(range(3))
+        ref.sxdg(range(3))
 
         self.assertEqual(ref, circuit)
 

From 8d5e3f3bb3a651953730ea52f3d2d442e26adec0 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 4 Nov 2024 16:29:42 +0100
Subject: [PATCH 25/25] Shelly's review comments

---
 crates/accelerate/src/circuit_library/mod.rs  |  2 +-
 .../src/circuit_library/pauli_evolution.rs    | 26 +++++++++++++++++++
 .../pauli-evo-plugins-612850146c3f7d49.yaml   |  7 +++--
 3 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/crates/accelerate/src/circuit_library/mod.rs b/crates/accelerate/src/circuit_library/mod.rs
index 3066ea67f1a7..d2e8755b2d90 100644
--- a/crates/accelerate/src/circuit_library/mod.rs
+++ b/crates/accelerate/src/circuit_library/mod.rs
@@ -13,8 +13,8 @@
 use pyo3::prelude::*;
 
 mod entanglement;
-mod pauli_evolution;
 mod iqp;
+mod pauli_evolution;
 mod pauli_feature_map;
 mod quantum_volume;
 
diff --git a/crates/accelerate/src/circuit_library/pauli_evolution.rs b/crates/accelerate/src/circuit_library/pauli_evolution.rs
index d110310e9e45..be41b51347bf 100644
--- a/crates/accelerate/src/circuit_library/pauli_evolution.rs
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -282,6 +282,19 @@ pub fn py_pauli_evolution(
     CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)
 }
 
+/// Build a CX chain over the active qubits. E.g. with q_1 inactive, this would return
+///
+///                    ┌───┐
+///     q_0: ──────────┤ X ├
+///                    └─┬─┘
+///     q_1: ────────────┼──
+///               ┌───┐  │
+///     q_2: ─────┤ X ├──■──
+///          ┌───┐└─┬─┘
+///     q_3: ┤ X ├──■───────
+///          └─┬─┘
+///     q_4: ──■────────────
+///     
 fn cx_chain(
     active_paulis: Vec<(char, Qubit)>,
 ) -> Box<dyn DoubleEndedIterator<Item = StandardInstruction>> {
@@ -293,6 +306,19 @@ fn cx_chain(
     )
 }
 
+/// Build a CX fountain over the active qubits. E.g. with q_1 inactive, this would return
+///
+////         ┌───┐┌───┐┌───┐
+////    q_0: ┤ X ├┤ X ├┤ X ├
+////         └─┬─┘└─┬─┘└─┬─┘
+////    q_1: ──┼────┼────┼──
+////           │    │    │
+////    q_2: ──■────┼────┼──
+////                │    │
+////    q_3: ───────■────┼──
+////                     │
+////    q_4: ────────────■──
+///     
 fn cx_fountain(
     active_paulis: Vec<(char, Qubit)>,
 ) -> Box<dyn DoubleEndedIterator<Item = StandardInstruction>> {
diff --git a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
index 67ccdfbae3b5..e4693aa4f3f0 100644
--- a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
+++ b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
@@ -19,8 +19,11 @@ features_synthesis:
     in a sparse Pauli format.
   - |
     Added the plugin structure for the :class:`.PauliEvolutionGate`. The default plugin,
-    :class:`.PauliEvolutionSynthesisDefault`, constructs circuit as before (but faster as it
-    internally uses Rust).
+    :class:`.PauliEvolutionSynthesisDefault`, constructs circuit as before, but faster as it
+    internally uses Rust. The larger the circuit (e.g. by the Hamiltonian size, the number 
+    of timesteps, or the Suzuki-Trotter order), the higher the speedup. For example,
+    a 100-qubit Heisenberg Hamiltonian with 10 timesteps and a 4th-order Trotter formula is 
+    now constructed ~9.4x faster.
 
 upgrade:
   - |