diff --git a/crates/accelerate/src/circuit_library/mod.rs b/crates/accelerate/src/circuit_library/mod.rs
index 5fa711070b48..d2e8755b2d90 100644
--- a/crates/accelerate/src/circuit_library/mod.rs
+++ b/crates/accelerate/src/circuit_library/mod.rs
@@ -14,10 +14,12 @@ use pyo3::prelude::*;
 
 mod entanglement;
 mod iqp;
+mod pauli_evolution;
 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!(iqp::py_iqp))?;
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..be41b51347bf
--- /dev/null
+++ b/crates/accelerate/src/circuit_library/pauli_evolution.rs
@@ -0,0 +1,356 @@
+// 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 pyo3::types::{PyList, PyString, PyTuple};
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::{multiply_param, radd_param, Param, PyInstruction, StandardGate};
+use qiskit_circuit::packed_instruction::PackedOperation;
+use qiskit_circuit::{imports, Clbit, Qubit};
+use smallvec::{smallvec, SmallVec};
+
+// custom types for a more readable code
+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).
+///
+/// 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.
+///     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.
+///
+/// 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,
+    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()
+        .zip(indices)
+        .filter(|(pauli, _)| *pauli != 'i');
+    let (paulis, indices): (Vec<char>, Vec<u32>) = active.unzip();
+
+    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, indices, time),
+        _ => Box::new(multi_qubit_evolution(
+            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];
+
+    std::iter::once(match pauli {
+        'x' => (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."),
+    })
+}
+
+/// 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>,
+    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))),
+        // 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)),
+    }
+}
+
+/// 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()
+        .zip(indices.into_iter().map(Qubit))
+        .collect();
+
+    // get the basis change: x -> HGate, y -> SXdgGate, z -> nothing
+    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::SXGate, smallvec![], smallvec![*q]),
+            _ => unreachable!("Invalid Pauli string."), // "z" and "i" have been filtered out
+        })
+        .collect();
+
+    // get the inverse basis change
+    let inverse_basis_change: Vec<StandardInstruction> = basis_change
+        .iter()
+        .map(|(gate, _, qubit)| match gate {
+            StandardGate::HGate => (StandardGate::HGate, smallvec![], qubit.clone()),
+            StandardGate::SXGate => (StandardGate::SXdgGate, 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 {
+        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 first qubit
+    let first_qubit = active_paulis.first().unwrap().1;
+    let z_rotation = std::iter::once((
+        if phase_gate {
+            StandardGate::PhaseGate
+        } else {
+            StandardGate::RZGate
+        },
+        smallvec![time],
+        smallvec![first_qubit],
+    ));
+
+    // and finally chain everything together
+    basis_change
+        .into_iter()
+        .chain(chain_down)
+        .chain(z_rotation)
+        .chain(chain_up)
+        .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: ┤ √Xdg ├┤ X ├──■─────────────■──┤ X ├┤ √X ├
+///    └──────┘└─┬─┘                   └─┬─┘└────┘
+/// 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.
+///     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.
+///
+/// Returns:
+///     Circuit data for to implement the evolution.
+#[pyfunction]
+#[pyo3(name = "pauli_evolution", signature = (num_qubits, sparse_paulis, insert_barriers=false, do_fountain=false))]
+pub fn py_pauli_evolution(
+    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>()?;
+        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);
+            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)?.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(
+                |(gate, params, qubits)| -> PyResult<Instruction> {
+                    Ok((
+                        gate.into(),
+                        params,
+                        Vec::from_iter(qubits.into_iter()),
+                        Vec::new(),
+                    ))
+                },
+            );
+
+            // 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)
+        },
+    );
+
+    // 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)
+}
+
+/// 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>> {
+    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])),
+    )
+}
+
+/// 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>> {
+    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],
+        )
+    }))
+}
+
+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/pauli_feature_map.rs b/crates/accelerate/src/circuit_library/pauli_feature_map.rs
index 6fa88d187182..bb9f8c25eb24 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,17 @@ 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.into_iter().rev().collect(),
+                multiply_param(&angle, alpha, py),
+                true,
+                false,
+            )
+            .map(|(gate, params, qargs)| {
+                (gate.into(), params, qargs.to_vec(), vec![] as Vec<Clbit>)
+            })
+            .collect::<Vec<Instruction>>();
             insts.extend(evo);
         }
     }
diff --git a/crates/circuit/src/operations.rs b/crates/circuit/src/operations.rs
index f8c1e68c0da9..a7b147c7e371 100644
--- a/crates/circuit/src/operations.rs
+++ b/crates/circuit/src/operations.rs
@@ -2305,7 +2305,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."),
     }
 }
 
@@ -2339,7 +2339,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)] => {
diff --git a/pyproject.toml b/pyproject.toml
index 326d28f2a273..91ef14f2e84d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -101,6 +101,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"
+"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/circuit/library/data_preparation/pauli_feature_map.py b/qiskit/circuit/library/data_preparation/pauli_feature_map.py
index 8c04f37ca9b0..1bc14d169273 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.sxdg(i)
+                    else:
+                        circuit.sx(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/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index 8bf7b3848c46..e8e19bd7af56 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..5b4e91721a6e 100644
--- a/qiskit/synthesis/evolution/lie_trotter.py
+++ b/qiskit/synthesis/evolution/lie_trotter.py
@@ -14,18 +14,15 @@
 
 from __future__ import annotations
 
-import inspect
 from collections.abc import Callable
 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
 
-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
@@ -40,7 +37,7 @@ class LieTrotter(ProductFormula):
 
     .. 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:
 
@@ -52,21 +49,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,26 +82,6 @@ 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
-
-        # 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()]
-        else:
-            pauli_list = [(op, 1) for op in operators]
-
-        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()
-
     @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 28bbc5711bf8..4de838215ff1 100644
--- a/qiskit/synthesis/evolution/product_formula.py
+++ b/qiskit/synthesis/evolution/product_formula.py
@@ -16,16 +16,19 @@
 
 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
+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 pauli_evolution
 
 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.
@@ -78,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
@@ -88,14 +92,16 @@ def __init__(
 
         # user-provided atomic evolution, stored for serialization
         self._atomic_evolution = atomic_evolution
+
+        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 = partial(
-                _default_atomic_evolution, cx_structure=cx_structure, wrap=wrap
-            )
+            self.atomic_evolution = None
 
         elif len(inspect.signature(atomic_evolution).parameters) == 2:
 
@@ -108,8 +114,50 @@ def wrap_atomic_evolution(output, operator, time):
         else:
             self.atomic_evolution = atomic_evolution
 
+    def expand(
+        self, evolution: PauliEvolutionGate
+    ) -> list[tuple[str, tuple[int], ParameterValueType]]:
+        """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", [1], 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)
+        num_qubits = evolution.num_qubits
+
+        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
+            cx_fountain = self._cx_structure == "fountain"
+            data = pauli_evolution(num_qubits, pauli_rotations, self.insert_barriers, cx_fountain)
+            circuit = QuantumCircuit._from_circuit_data(data, add_regs=True)
+
+        return circuit
+
     @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:
@@ -131,254 +179,63 @@ def settings(self) -> dict[str, 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,
-    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
+    def _custom_evolution(self, num_qubits, pauli_rotations):
+        """Implement the evolution for the non-standard path.
 
-    .. math::
+        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)
+        cx_fountain = self._cx_structure == "fountain"
 
-        U(t) = e^{-itP}.
+        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
 
-    Since only a single Pauli string is evolved the circuit decomposition is exact.
+            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)
 
-    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)
+                # build the circuit Rust-side
+                data = pauli_evolution(
+                    len(qubits),
+                    [local_pauli],
+                    False,
+                    cx_fountain,
+                )
+                evo = QuantumCircuit._from_circuit_data(data)
 
-    return cliff
+                # 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()
 
-def cnot_chain(pauli: Pauli) -> QuantumCircuit:
-    """CX chain.
+        return circuit
 
-    For example, for the Pauli with the label 'XYZIX'.
 
-    .. code-block:: text
+def real_or_fail(value, tol=100):
+    """Return real if close, otherwise fail. Unbound parameters are left unchanged.
 
-                       ┌───┐
-        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.
+    Based on NumPy's ``real_if_close``, i.e. ``tol`` is in terms of machine precision for float.
     """
+    if isinstance(value, ParameterExpression):
+        return value
 
-    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'.
-
-    .. code-block:: text
-
-             ┌───┐┌───┐┌───┐
-        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.
-    """
+    abstol = tol * np.finfo(float).eps
+    if abs(np.imag(value)) < abstol:
+        return np.real(value)
 
-    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)
+    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 7c68f66dc8ea..64ae5fba637b 100644
--- a/qiskit/synthesis/evolution/qdrift.py
+++ b/qiskit/synthesis/evolution/qdrift.py
@@ -16,14 +16,19 @@
 
 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
 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
-from .lie_trotter import LieTrotter
+
+if typing.TYPE_CHECKING:
+    from qiskit.circuit.library import PauliEvolutionGate
 
 
 class QDrift(ProductFormula):
@@ -88,44 +93,36 @@ 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 determine 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()
+
+        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)
         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 = 2 * lambd / num_gates * time
+        sampled_paulis = [
+            (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 e03fd27e26d4..4dfa92579bcb 100644
--- a/qiskit/synthesis/evolution/suzuki_trotter.py
+++ b/qiskit/synthesis/evolution/suzuki_trotter.py
@@ -15,17 +15,20 @@
 from __future__ import annotations
 
 import inspect
+import typing
 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
 from qiskit.utils.deprecation import deprecate_arg
 
-
 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.
@@ -44,7 +47,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,
@@ -105,51 +108,91 @@ 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."
+                f"for when the order is 1 or even, not {order}."
             )
         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
+    def expand(
+        self, evolution: PauliEvolutionGate
+    ) -> list[tuple[str, list[int], ParameterValueType]]:
+        """Expand the Hamiltonian into a Suzuki-Trotter sequence of sparse gates.
 
-        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]
+        For example, the Hamiltonian ``H = IX + ZZ`` for an evolution time ``t`` and
+        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), ("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)``.
 
-        ops_to_evolve = self._recurse(self.order, time / self.reps, pauli_list)
+        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
 
         # construct the evolution circuit
-        single_rep = QuantumCircuit(operators[0].num_qubits)
+        if isinstance(operators, list):  # already sorted into commuting bits
+            non_commuting = [
+                (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 (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)
 
-        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()
+        # 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 single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
+        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/qiskit/transpiler/passes/synthesis/hls_plugins.py b/qiskit/transpiler/passes/synthesis/hls_plugins.py
index fa827a20ed48..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
 '''''''''''''''''''''
@@ -251,6 +264,8 @@
 
 """
 
+from __future__ import annotations
+
 import numpy as np
 import rustworkx as rx
 
@@ -288,8 +303,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 +1044,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)
diff --git a/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
new file mode 100644
index 000000000000..e4693aa4f3f0
--- /dev/null
+++ b/releasenotes/notes/pauli-evo-plugins-612850146c3f7d49.yaml
@@ -0,0 +1,35 @@
+---
+features_quantum_info:
+  - |
+    Added :meth:`.SparsePauliOperator.to_sparse_list` to convert an operator into
+    a sparse list format. This works inversely 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. 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:
+  - | 
+    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 a7c9a73abb5f..7d1740b67f30 100644
--- a/test/python/circuit/library/test_evolution_gate.py
+++ b/test/python/circuit/library/test_evolution_gate.py
@@ -18,9 +18,8 @@
 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.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
@@ -40,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)
@@ -59,7 +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 ^ Y  # use a string for which we do not have a basis gate
+        time = 0.321
+        evo_gate = PauliEvolutionGate(op, time)
+        self.assertSuzukiTrotterIsCorrect(evo_gate)
 
     def test_rzx_order(self):
         """Test ZX and XZ is mapped onto the correct qubits."""
@@ -105,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."""
@@ -132,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)]),
@@ -185,6 +208,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)
@@ -251,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"),
@@ -275,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."""
@@ -286,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.
@@ -294,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."""
@@ -336,6 +362,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)
 
@@ -365,5 +397,65 @@ 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.sx(i)
+        elif pauli_i == "X":
+            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()
diff --git a/test/python/circuit/library/test_evolved_op_ansatz.py b/test/python/circuit/library/test_evolved_op_ansatz.py
index 8ea66b7012fe..2dbcec86a120 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.sx(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..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: ┤ Rx(π/2) ├┤ X ├──■──────────────■──┤ X ├┤ Rx(-π/2) ├
-            #      └──┬───┬──┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└──────────┘
-            # 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.rx(np.pi / 2, 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.rx(-np.pi / 2, 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: ┤ Rx(π/2) ├┤ X ├──■──────────────■──┤ X ├┤ Rx(-π/2) ├
-        #      └──┬───┬──┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤└──────────┘
-        # 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.rx(np.pi / 2, 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.rx(-np.pi / 2, 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.rx(np.pi / 2, 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.rx(-np.pi / 2, range(3))
+        ref.sxdg(range(3))
 
         self.assertEqual(ref, circuit)