diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs
index 01725269bb84..24c4f6e87c2a 100644
--- a/crates/accelerate/src/euler_one_qubit_decomposer.rs
+++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs
@@ -18,7 +18,6 @@ use num_complex::{Complex64, ComplexFloat};
 use smallvec::{smallvec, SmallVec};
 use std::cmp::Ordering;
 use std::f64::consts::PI;
-use std::ops::Deref;
 use std::str::FromStr;
 
 use pyo3::exceptions::PyValueError;
@@ -31,8 +30,12 @@ use ndarray::prelude::*;
 use numpy::PyReadonlyArray2;
 use pyo3::pybacked::PyBackedStr;
 
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::dag_node::DAGOpNode;
+use qiskit_circuit::operations::{Operation, Param, StandardGate};
 use qiskit_circuit::slice::{PySequenceIndex, SequenceIndex};
 use qiskit_circuit::util::c64;
+use qiskit_circuit::Qubit;
 
 pub const ANGLE_ZERO_EPSILON: f64 = 1e-12;
 
@@ -68,12 +71,12 @@ impl OneQubitGateErrorMap {
 
 #[pyclass(sequence)]
 pub struct OneQubitGateSequence {
-    pub gates: Vec<(String, SmallVec<[f64; 3]>)>,
+    pub gates: Vec<(StandardGate, SmallVec<[f64; 3]>)>,
     #[pyo3(get)]
     pub global_phase: f64,
 }
 
-type OneQubitGateSequenceState = (Vec<(String, SmallVec<[f64; 3]>)>, f64);
+type OneQubitGateSequenceState = (Vec<(StandardGate, SmallVec<[f64; 3]>)>, f64);
 
 #[pymethods]
 impl OneQubitGateSequence {
@@ -115,15 +118,15 @@ fn circuit_kak(
     phi: f64,
     lam: f64,
     phase: f64,
-    k_gate: &str,
-    a_gate: &str,
+    k_gate: StandardGate,
+    a_gate: StandardGate,
     simplify: bool,
     atol: Option<f64>,
 ) -> OneQubitGateSequence {
     let mut lam = lam;
     let mut theta = theta;
     let mut phi = phi;
-    let mut circuit: Vec<(String, SmallVec<[f64; 3]>)> = Vec::with_capacity(3);
+    let mut circuit: Vec<(StandardGate, SmallVec<[f64; 3]>)> = Vec::with_capacity(3);
     let mut atol = match atol {
         Some(atol) => atol,
         None => ANGLE_ZERO_EPSILON,
@@ -139,7 +142,7 @@ fn circuit_kak(
         // slippage coming from _mod_2pi injecting multiples of 2pi.
         lam = mod_2pi(lam, atol);
         if lam.abs() > atol {
-            circuit.push((String::from(k_gate), smallvec![lam]));
+            circuit.push((k_gate, smallvec![lam]));
             global_phase += lam / 2.;
         }
         return OneQubitGateSequence {
@@ -160,13 +163,13 @@ fn circuit_kak(
     lam = mod_2pi(lam, atol);
     if lam.abs() > atol {
         global_phase += lam / 2.;
-        circuit.push((String::from(k_gate), smallvec![lam]));
+        circuit.push((k_gate, smallvec![lam]));
     }
-    circuit.push((String::from(a_gate), smallvec![theta]));
+    circuit.push((a_gate, smallvec![theta]));
     phi = mod_2pi(phi, atol);
     if phi.abs() > atol {
         global_phase += phi / 2.;
-        circuit.push((String::from(k_gate), smallvec![phi]));
+        circuit.push((k_gate, smallvec![phi]));
     }
     OneQubitGateSequence {
         gates: circuit,
@@ -190,7 +193,7 @@ fn circuit_u3(
     let phi = mod_2pi(phi, atol);
     let lam = mod_2pi(lam, atol);
     if !simplify || theta.abs() > atol || phi.abs() > atol || lam.abs() > atol {
-        circuit.push((String::from("u3"), smallvec![theta, phi, lam]));
+        circuit.push((StandardGate::U3Gate, smallvec![theta, phi, lam]));
     }
     OneQubitGateSequence {
         gates: circuit,
@@ -217,16 +220,16 @@ fn circuit_u321(
     if theta.abs() < atol {
         let tot = mod_2pi(phi + lam, atol);
         if tot.abs() > atol {
-            circuit.push((String::from("u1"), smallvec![tot]));
+            circuit.push((StandardGate::U1Gate, smallvec![tot]));
         }
     } else if (theta - PI / 2.).abs() < atol {
         circuit.push((
-            String::from("u2"),
+            StandardGate::U2Gate,
             smallvec![mod_2pi(phi, atol), mod_2pi(lam, atol)],
         ));
     } else {
         circuit.push((
-            String::from("u3"),
+            StandardGate::U3Gate,
             smallvec![theta, mod_2pi(phi, atol), mod_2pi(lam, atol)],
         ));
     }
@@ -255,7 +258,7 @@ fn circuit_u(
     let phi = mod_2pi(phi, atol);
     let lam = mod_2pi(lam, atol);
     if theta.abs() > atol || phi.abs() > atol || lam.abs() > atol {
-        circuit.push((String::from("u"), smallvec![theta, phi, lam]));
+        circuit.push((StandardGate::UGate, smallvec![theta, phi, lam]));
     }
     OneQubitGateSequence {
         gates: circuit,
@@ -358,7 +361,7 @@ fn circuit_rr(
         // This can be expressed as a single R gate
         if theta.abs() > atol {
             circuit.push((
-                String::from("r"),
+                StandardGate::RGate,
                 smallvec![theta, mod_2pi(PI / 2. + phi, atol)],
             ));
         }
@@ -366,12 +369,12 @@ fn circuit_rr(
         // General case: use two R gates
         if (theta - PI).abs() > atol {
             circuit.push((
-                String::from("r"),
+                StandardGate::RGate,
                 smallvec![theta - PI, mod_2pi(PI / 2. - lam, atol)],
             ));
         }
         circuit.push((
-            String::from("r"),
+            StandardGate::RGate,
             smallvec![PI, mod_2pi(0.5 * (phi - lam + PI), atol)],
         ));
     }
@@ -393,10 +396,46 @@ pub fn generate_circuit(
     atol: Option<f64>,
 ) -> PyResult<OneQubitGateSequence> {
     let res = match target_basis {
-        EulerBasis::ZYZ => circuit_kak(theta, phi, lam, phase, "rz", "ry", simplify, atol),
-        EulerBasis::ZXZ => circuit_kak(theta, phi, lam, phase, "rz", "rx", simplify, atol),
-        EulerBasis::XZX => circuit_kak(theta, phi, lam, phase, "rx", "rz", simplify, atol),
-        EulerBasis::XYX => circuit_kak(theta, phi, lam, phase, "rx", "ry", simplify, atol),
+        EulerBasis::ZYZ => circuit_kak(
+            theta,
+            phi,
+            lam,
+            phase,
+            StandardGate::RZGate,
+            StandardGate::RYGate,
+            simplify,
+            atol,
+        ),
+        EulerBasis::ZXZ => circuit_kak(
+            theta,
+            phi,
+            lam,
+            phase,
+            StandardGate::RZGate,
+            StandardGate::RXGate,
+            simplify,
+            atol,
+        ),
+        EulerBasis::XZX => circuit_kak(
+            theta,
+            phi,
+            lam,
+            phase,
+            StandardGate::RXGate,
+            StandardGate::RZGate,
+            simplify,
+            atol,
+        ),
+        EulerBasis::XYX => circuit_kak(
+            theta,
+            phi,
+            lam,
+            phase,
+            StandardGate::RXGate,
+            StandardGate::RYGate,
+            simplify,
+            atol,
+        ),
         EulerBasis::U3 => circuit_u3(theta, phi, lam, phase, simplify, atol),
         EulerBasis::U321 => circuit_u321(theta, phi, lam, phase, simplify, atol),
         EulerBasis::U => circuit_u(theta, phi, lam, phase, simplify, atol),
@@ -411,11 +450,13 @@ pub fn generate_circuit(
             let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| {
                 let phi = mod_2pi(phi, inner_atol);
                 if phi.abs() > inner_atol {
-                    circuit.gates.push((String::from("p"), smallvec![phi]));
+                    circuit
+                        .gates
+                        .push((StandardGate::PhaseGate, smallvec![phi]));
                 }
             };
             let fnx = |circuit: &mut OneQubitGateSequence| {
-                circuit.gates.push((String::from("sx"), SmallVec::new()));
+                circuit.gates.push((StandardGate::SXGate, SmallVec::new()));
             };
 
             circuit_psx_gen(
@@ -441,12 +482,12 @@ pub fn generate_circuit(
             let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| {
                 let phi = mod_2pi(phi, inner_atol);
                 if phi.abs() > inner_atol {
-                    circuit.gates.push((String::from("rz"), smallvec![phi]));
+                    circuit.gates.push((StandardGate::RZGate, smallvec![phi]));
                     circuit.global_phase += phi / 2.;
                 }
             };
             let fnx = |circuit: &mut OneQubitGateSequence| {
-                circuit.gates.push((String::from("sx"), SmallVec::new()));
+                circuit.gates.push((StandardGate::SXGate, SmallVec::new()));
             };
             circuit_psx_gen(
                 theta,
@@ -471,12 +512,14 @@ pub fn generate_circuit(
             let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| {
                 let phi = mod_2pi(phi, inner_atol);
                 if phi.abs() > inner_atol {
-                    circuit.gates.push((String::from("u1"), smallvec![phi]));
+                    circuit.gates.push((StandardGate::U1Gate, smallvec![phi]));
                 }
             };
             let fnx = |circuit: &mut OneQubitGateSequence| {
                 circuit.global_phase += PI / 4.;
-                circuit.gates.push((String::from("rx"), smallvec![PI / 2.]));
+                circuit
+                    .gates
+                    .push((StandardGate::RXGate, smallvec![PI / 2.]));
             };
             circuit_psx_gen(
                 theta,
@@ -501,15 +544,15 @@ pub fn generate_circuit(
             let fnz = |circuit: &mut OneQubitGateSequence, phi: f64| {
                 let phi = mod_2pi(phi, inner_atol);
                 if phi.abs() > inner_atol {
-                    circuit.gates.push((String::from("rz"), smallvec![phi]));
+                    circuit.gates.push((StandardGate::RZGate, smallvec![phi]));
                     circuit.global_phase += phi / 2.;
                 }
             };
             let fnx = |circuit: &mut OneQubitGateSequence| {
-                circuit.gates.push((String::from("sx"), SmallVec::new()));
+                circuit.gates.push((StandardGate::SXGate, SmallVec::new()));
             };
             let fnxpi = |circuit: &mut OneQubitGateSequence| {
-                circuit.gates.push((String::from("x"), SmallVec::new()));
+                circuit.gates.push((StandardGate::XGate, SmallVec::new()));
             };
             circuit_psx_gen(
                 theta,
@@ -633,7 +676,7 @@ fn compare_error_fn(
             let fidelity_product: f64 = circuit
                 .gates
                 .iter()
-                .map(|x| 1. - err_map.get(&x.0).unwrap_or(&0.))
+                .map(|gate| 1. - err_map.get(gate.0.name()).unwrap_or(&0.))
                 .product();
             (1. - fidelity_product, circuit.gates.len())
         }
@@ -642,6 +685,28 @@ fn compare_error_fn(
 }
 
 fn compute_error(
+    gates: &[(StandardGate, SmallVec<[f64; 3]>)],
+    error_map: Option<&OneQubitGateErrorMap>,
+    qubit: usize,
+) -> (f64, usize) {
+    match error_map {
+        Some(err_map) => {
+            let num_gates = gates.len();
+            let gate_fidelities: f64 = gates
+                .iter()
+                .map(|gate| 1. - err_map.error_map[qubit].get(gate.0.name()).unwrap_or(&0.))
+                .product();
+            (1. - gate_fidelities, num_gates)
+        }
+        None => (gates.len() as f64, gates.len()),
+    }
+}
+
+fn compute_error_term(gate: &str, error_map: &OneQubitGateErrorMap, qubit: usize) -> f64 {
+    1. - error_map.error_map[qubit].get(gate).unwrap_or(&0.)
+}
+
+fn compute_error_str(
     gates: &[(String, SmallVec<[f64; 3]>)],
     error_map: Option<&OneQubitGateErrorMap>,
     qubit: usize,
@@ -651,7 +716,7 @@ fn compute_error(
             let num_gates = gates.len();
             let gate_fidelities: f64 = gates
                 .iter()
-                .map(|x| 1. - err_map.error_map[qubit].get(&x.0).unwrap_or(&0.))
+                .map(|gate| compute_error_term(gate.0.as_str(), err_map, qubit))
                 .product();
             (1. - gate_fidelities, num_gates)
         }
@@ -670,11 +735,20 @@ pub fn compute_error_one_qubit_sequence(
 
 #[pyfunction]
 pub fn compute_error_list(
-    circuit: Vec<(String, SmallVec<[f64; 3]>)>,
+    circuit: Vec<PyRef<DAGOpNode>>,
     qubit: usize,
     error_map: Option<&OneQubitGateErrorMap>,
 ) -> (f64, usize) {
-    compute_error(&circuit, error_map, qubit)
+    let circuit_list: Vec<(String, SmallVec<[f64; 3]>)> = circuit
+        .iter()
+        .map(|node| {
+            (
+                node.instruction.operation.name().to_string(),
+                smallvec![], // Params not needed in this path
+            )
+        })
+        .collect();
+    compute_error_str(&circuit_list, error_map, qubit)
 }
 
 #[pyfunction]
@@ -687,15 +761,13 @@ pub fn unitary_to_gate_sequence(
     simplify: bool,
     atol: Option<f64>,
 ) -> PyResult<Option<OneQubitGateSequence>> {
-    let mut target_basis_vec: Vec<EulerBasis> = Vec::with_capacity(target_basis_list.len());
-    for basis in target_basis_list {
-        let basis_enum = EulerBasis::__new__(basis.deref())?;
-        target_basis_vec.push(basis_enum)
-    }
-    let unitary_mat = unitary.as_array();
+    let target_basis_vec: PyResult<Vec<EulerBasis>> = target_basis_list
+        .iter()
+        .map(|basis| EulerBasis::__new__(basis))
+        .collect();
     Ok(unitary_to_gate_sequence_inner(
-        unitary_mat,
-        &target_basis_vec,
+        unitary.as_array(),
+        &target_basis_vec?,
         qubit,
         error_map,
         simplify,
@@ -725,6 +797,46 @@ pub fn unitary_to_gate_sequence_inner(
         })
 }
 
+#[pyfunction]
+#[pyo3(signature = (unitary, target_basis_list, qubit, error_map=None, simplify=true, atol=None))]
+pub fn unitary_to_circuit(
+    py: Python,
+    unitary: PyReadonlyArray2<Complex64>,
+    target_basis_list: Vec<PyBackedStr>,
+    qubit: usize,
+    error_map: Option<&OneQubitGateErrorMap>,
+    simplify: bool,
+    atol: Option<f64>,
+) -> PyResult<Option<CircuitData>> {
+    let target_basis_vec: PyResult<Vec<EulerBasis>> = target_basis_list
+        .iter()
+        .map(|basis| EulerBasis::__new__(basis))
+        .collect();
+    let circuit_sequence = unitary_to_gate_sequence_inner(
+        unitary.as_array(),
+        &target_basis_vec?,
+        qubit,
+        error_map,
+        simplify,
+        atol,
+    );
+    Ok(circuit_sequence.map(|seq| {
+        CircuitData::from_standard_gates(
+            py,
+            1,
+            seq.gates.into_iter().map(|(gate, params)| {
+                (
+                    gate,
+                    params.into_iter().map(Param::Float).collect(),
+                    smallvec![Qubit(0)],
+                )
+            }),
+            Param::Float(seq.global_phase),
+        )
+        .expect("Unexpected Qiskit python bug")
+    }))
+}
+
 #[inline]
 pub fn det_one_qubit(mat: ArrayView2<Complex64>) -> Complex64 {
     mat[[0, 0]] * mat[[1, 1]] - mat[[0, 1]] * mat[[1, 0]]
@@ -853,6 +965,107 @@ pub fn params_zxz(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
     params_zxz_inner(mat)
 }
 
+type OptimizeDecompositionReturn = Option<((f64, usize), (f64, usize), OneQubitGateSequence)>;
+
+#[pyfunction]
+pub fn optimize_1q_gates_decomposition(
+    runs: Vec<Vec<PyRef<DAGOpNode>>>,
+    qubits: Vec<usize>,
+    bases: Vec<Vec<PyBackedStr>>,
+    simplify: bool,
+    error_map: Option<&OneQubitGateErrorMap>,
+    atol: Option<f64>,
+) -> Vec<OptimizeDecompositionReturn> {
+    runs.iter()
+        .enumerate()
+        .map(|(index, raw_run)| -> OptimizeDecompositionReturn {
+            let mut error = match error_map {
+                Some(_) => 1.,
+                None => raw_run.len() as f64,
+            };
+            let qubit = qubits[index];
+            let operator = &raw_run
+                .iter()
+                .map(|node| {
+                    if let Some(err_map) = error_map {
+                        error *=
+                            compute_error_term(node.instruction.operation.name(), err_map, qubit)
+                    }
+                    node.instruction
+                        .operation
+                        .matrix(&node.instruction.params)
+                        .expect("No matrix defined for operation")
+                })
+                .fold(
+                    [
+                        [Complex64::new(1., 0.), Complex64::new(0., 0.)],
+                        [Complex64::new(0., 0.), Complex64::new(1., 0.)],
+                    ],
+                    |mut operator, node| {
+                        matmul_1q(&mut operator, node);
+                        operator
+                    },
+                );
+            let old_error = if error_map.is_some() {
+                (1. - error, raw_run.len())
+            } else {
+                (error, raw_run.len())
+            };
+            let target_basis_vec: Vec<EulerBasis> = bases[index]
+                .iter()
+                .map(|basis| EulerBasis::__new__(basis).unwrap())
+                .collect();
+            unitary_to_gate_sequence_inner(
+                aview2(operator),
+                &target_basis_vec,
+                qubit,
+                error_map,
+                simplify,
+                atol,
+            )
+            .map(|out_seq| {
+                let new_error = compute_error_one_qubit_sequence(&out_seq, qubit, error_map);
+                (old_error, new_error, out_seq)
+            })
+        })
+        .collect()
+}
+
+fn matmul_1q(operator: &mut [[Complex64; 2]; 2], other: Array2<Complex64>) {
+    *operator = [
+        [
+            other[[0, 0]] * operator[0][0] + other[[0, 1]] * operator[1][0],
+            other[[0, 0]] * operator[0][1] + other[[0, 1]] * operator[1][1],
+        ],
+        [
+            other[[1, 0]] * operator[0][0] + other[[1, 1]] * operator[1][0],
+            other[[1, 0]] * operator[0][1] + other[[1, 1]] * operator[1][1],
+        ],
+    ];
+}
+
+#[pyfunction]
+pub fn collect_1q_runs_filter(node: &Bound<PyAny>) -> bool {
+    let op_node = node.downcast::<DAGOpNode>();
+    match op_node {
+        Ok(bound_node) => {
+            let node = bound_node.borrow();
+            node.instruction.operation.num_qubits() == 1
+                && node.instruction.operation.num_clbits() == 0
+                && node
+                    .instruction
+                    .operation
+                    .matrix(&node.instruction.params)
+                    .is_some()
+                && match &node.instruction.extra_attrs {
+                    None => true,
+                    Some(attrs) => attrs.condition.is_none(),
+                }
+        }
+        Err(_) => false,
+    }
+}
+
 #[pymodule]
 pub fn euler_one_qubit_decomposer(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(params_zyz))?;
@@ -863,8 +1076,11 @@ pub fn euler_one_qubit_decomposer(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(params_u1x))?;
     m.add_wrapped(wrap_pyfunction!(generate_circuit))?;
     m.add_wrapped(wrap_pyfunction!(unitary_to_gate_sequence))?;
+    m.add_wrapped(wrap_pyfunction!(unitary_to_circuit))?;
     m.add_wrapped(wrap_pyfunction!(compute_error_one_qubit_sequence))?;
     m.add_wrapped(wrap_pyfunction!(compute_error_list))?;
+    m.add_wrapped(wrap_pyfunction!(optimize_1q_gates_decomposition))?;
+    m.add_wrapped(wrap_pyfunction!(collect_1q_runs_filter))?;
     m.add_class::<OneQubitGateSequence>()?;
     m.add_class::<OneQubitGateErrorMap>()?;
     m.add_class::<EulerBasis>()?;
diff --git a/crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs b/crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs
new file mode 100644
index 000000000000..e53e25282005
--- /dev/null
+++ b/crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs
@@ -0,0 +1,441 @@
+// 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 indexmap::IndexSet;
+use ndarray::{s, ArrayView2};
+use smallvec::smallvec;
+
+use crate::synthesis::clifford::utils::CliffordGatesVec;
+use crate::synthesis::clifford::utils::{adjust_final_pauli_gates, SymplecticMatrix};
+use qiskit_circuit::operations::StandardGate;
+use qiskit_circuit::Qubit;
+
+/// Converts a pair of Paulis pauli_x and pauli_z acting on a specific qubit
+/// to the corresponding index in [PauliPairsClass] or [SingleQubitGate] classes.
+/// The input is given as a 4-tuple: (pauli_x stabilizer, pauli_x destabilizer,
+/// pauli_z stabilizer, pauli_z destabilizer), and the output is an unsigned
+/// integer from 0 to 15.
+fn pauli_pair_to_index(xs: bool, xd: bool, zs: bool, zd: bool) -> usize {
+    ((xs as usize) << 3) | ((xd as usize) << 2) | ((zs as usize) << 1) | (zd as usize)
+}
+
+/// The five classes of Pauli 2-qubit operators as described in the paper.
+#[derive(Clone, Copy)]
+enum PauliPairsClass {
+    ClassA,
+    ClassB,
+    ClassC,
+    ClassD,
+    ClassE,
+}
+
+/// The 16 Pauli 2-qubit operators are divided into 5 equivalence classes
+/// under the action of single-qubit Cliffords.
+static PAULI_INDEX_TO_CLASS: [PauliPairsClass; 16] = [
+    PauliPairsClass::ClassE, // 'II'
+    PauliPairsClass::ClassD, // 'IX'
+    PauliPairsClass::ClassD, // 'IZ'
+    PauliPairsClass::ClassD, // 'IY'
+    PauliPairsClass::ClassC, // 'XI'
+    PauliPairsClass::ClassB, // 'XX'
+    PauliPairsClass::ClassA, // 'XZ'
+    PauliPairsClass::ClassA, // 'XY'
+    PauliPairsClass::ClassC, // 'ZI'
+    PauliPairsClass::ClassA, // 'ZX'
+    PauliPairsClass::ClassB, // 'ZZ'
+    PauliPairsClass::ClassA, // 'ZY'
+    PauliPairsClass::ClassC, // 'YI'
+    PauliPairsClass::ClassA, // 'YX'
+    PauliPairsClass::ClassA, // 'YZ'
+    PauliPairsClass::ClassB, // 'YY'
+];
+
+/// Single-qubit Clifford gates modulo Paulis.
+#[derive(Clone, Copy)]
+enum SingleQubitGate {
+    GateI,
+    GateS,
+    GateH,
+    GateSH,
+    GateHS,
+    GateSHS,
+}
+
+/// Maps pair of pauli operators to the single-qubit gate required
+/// for the decoupling step.
+static PAULI_INDEX_TO_1Q_GATE: [SingleQubitGate; 16] = [
+    SingleQubitGate::GateI,   // 'II'
+    SingleQubitGate::GateH,   // 'IX'
+    SingleQubitGate::GateI,   // 'IZ'
+    SingleQubitGate::GateSH,  // 'IY'
+    SingleQubitGate::GateI,   // 'XI'
+    SingleQubitGate::GateI,   // 'XX'
+    SingleQubitGate::GateI,   // 'XZ'
+    SingleQubitGate::GateSHS, // 'XY'
+    SingleQubitGate::GateH,   // 'ZI'
+    SingleQubitGate::GateH,   // 'ZX'
+    SingleQubitGate::GateH,   // 'ZZ'
+    SingleQubitGate::GateSH,  // 'ZY'
+    SingleQubitGate::GateS,   // 'YI'
+    SingleQubitGate::GateHS,  // 'YX'
+    SingleQubitGate::GateS,   // 'YZ'
+    SingleQubitGate::GateS,   // 'YY'
+];
+
+pub struct GreedyCliffordSynthesis<'a> {
+    /// The Clifford tableau to be synthesized.
+    tableau: ArrayView2<'a, bool>,
+
+    /// The total number of qubits.
+    num_qubits: usize,
+
+    /// Symplectic matrix being reduced.
+    symplectic_matrix: SymplecticMatrix,
+
+    /// Unprocessed qubits.
+    unprocessed_qubits: IndexSet<usize>,
+}
+
+impl GreedyCliffordSynthesis<'_> {
+    pub(crate) fn new(tableau: ArrayView2<bool>) -> Result<GreedyCliffordSynthesis<'_>, String> {
+        let tableau_shape = tableau.shape();
+        if (tableau_shape[0] % 2 == 1) || (tableau_shape[1] != tableau_shape[0] + 1) {
+            return Err("The shape of the Clifford tableau is invalid".to_string());
+        }
+
+        let num_qubits = tableau_shape[0] / 2;
+
+        // We are going to modify symplectic_matrix in-place until it
+        // becomes the identity.
+        let symplectic_matrix = SymplecticMatrix {
+            num_qubits,
+            smat: tableau.slice(s![.., 0..2 * num_qubits]).to_owned(),
+        };
+
+        let unprocessed_qubits: IndexSet<usize> = (0..num_qubits).collect();
+
+        Ok(GreedyCliffordSynthesis {
+            tableau,
+            num_qubits,
+            symplectic_matrix,
+            unprocessed_qubits,
+        })
+    }
+
+    /// Computes the CX cost of decoupling the symplectic matrix on the
+    /// given qubit.
+    fn compute_cost(&self, qubit: usize) -> Result<usize, String> {
+        let mut a_num = 0;
+        let mut b_num = 0;
+        let mut c_num = 0;
+        let mut d_num = 0;
+
+        let mut qubit_is_in_a = false;
+
+        for q in &self.unprocessed_qubits {
+            let pauli_pair_index = pauli_pair_to_index(
+                self.symplectic_matrix.smat[[*q, qubit + self.num_qubits]],
+                self.symplectic_matrix.smat[[*q + self.num_qubits, qubit + self.num_qubits]],
+                self.symplectic_matrix.smat[[*q, qubit]],
+                self.symplectic_matrix.smat[[*q + self.num_qubits, qubit]],
+            );
+            let pauli_class = PAULI_INDEX_TO_CLASS[pauli_pair_index];
+
+            match pauli_class {
+                PauliPairsClass::ClassA => {
+                    a_num += 1;
+                    if *q == qubit {
+                        qubit_is_in_a = true;
+                    }
+                }
+                PauliPairsClass::ClassB => {
+                    b_num += 1;
+                }
+                PauliPairsClass::ClassC => {
+                    c_num += 1;
+                }
+                PauliPairsClass::ClassD => {
+                    d_num += 1;
+                }
+                PauliPairsClass::ClassE => {}
+            }
+        }
+
+        if a_num % 2 == 0 {
+            return Err("Symplectic Gaussian elimination failed.".to_string());
+        }
+
+        let mut cnot_cost: usize =
+            3 * (a_num - 1) / 2 + (b_num + 1) * ((b_num > 0) as usize) + c_num + d_num;
+
+        if !qubit_is_in_a {
+            cnot_cost += 3;
+        }
+
+        Ok(cnot_cost)
+    }
+
+    /// Calculate a decoupling operator D:
+    /// D^{-1} * Ox * D = x1
+    /// D^{-1} * Oz * D = z1
+    /// and reduces the clifford such that it will act trivially on min_qubit.
+    fn decouple_qubit(
+        &mut self,
+        gate_seq: &mut CliffordGatesVec,
+        min_qubit: usize,
+    ) -> Result<(), String> {
+        let mut a_qubits = IndexSet::new();
+        let mut b_qubits = IndexSet::new();
+        let mut c_qubits = IndexSet::new();
+        let mut d_qubits = IndexSet::new();
+
+        for qubit in &self.unprocessed_qubits {
+            let pauli_pair_index = pauli_pair_to_index(
+                self.symplectic_matrix.smat[[*qubit, min_qubit + self.num_qubits]],
+                self.symplectic_matrix.smat
+                    [[*qubit + self.num_qubits, min_qubit + self.num_qubits]],
+                self.symplectic_matrix.smat[[*qubit, min_qubit]],
+                self.symplectic_matrix.smat[[*qubit + self.num_qubits, min_qubit]],
+            );
+
+            let single_qubit_gate = PAULI_INDEX_TO_1Q_GATE[pauli_pair_index];
+            match single_qubit_gate {
+                SingleQubitGate::GateS => {
+                    gate_seq.push((
+                        StandardGate::SGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    self.symplectic_matrix.prepend_s(*qubit);
+                }
+                SingleQubitGate::GateH => {
+                    gate_seq.push((
+                        StandardGate::HGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    self.symplectic_matrix.prepend_h(*qubit);
+                }
+                SingleQubitGate::GateSH => {
+                    gate_seq.push((
+                        StandardGate::SGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    gate_seq.push((
+                        StandardGate::HGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    self.symplectic_matrix.prepend_s(*qubit);
+                    self.symplectic_matrix.prepend_h(*qubit);
+                }
+                SingleQubitGate::GateHS => {
+                    gate_seq.push((
+                        StandardGate::HGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    gate_seq.push((
+                        StandardGate::SGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    self.symplectic_matrix.prepend_h(*qubit);
+                    self.symplectic_matrix.prepend_s(*qubit);
+                }
+                SingleQubitGate::GateSHS => {
+                    gate_seq.push((
+                        StandardGate::SGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    gate_seq.push((
+                        StandardGate::HGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    gate_seq.push((
+                        StandardGate::SGate,
+                        smallvec![],
+                        smallvec![Qubit(*qubit as u32)],
+                    ));
+                    self.symplectic_matrix.prepend_s(*qubit);
+                    self.symplectic_matrix.prepend_h(*qubit);
+                    self.symplectic_matrix.prepend_s(*qubit);
+                }
+                SingleQubitGate::GateI => {}
+            }
+
+            let pauli_class = PAULI_INDEX_TO_CLASS[pauli_pair_index];
+            match pauli_class {
+                PauliPairsClass::ClassA => {
+                    a_qubits.insert(*qubit);
+                }
+                PauliPairsClass::ClassB => {
+                    b_qubits.insert(*qubit);
+                }
+                PauliPairsClass::ClassC => {
+                    c_qubits.insert(*qubit);
+                }
+                PauliPairsClass::ClassD => {
+                    d_qubits.insert(*qubit);
+                }
+                PauliPairsClass::ClassE => {}
+            }
+        }
+
+        if a_qubits.len() % 2 != 1 {
+            return Err("Symplectic Gaussian elimination failed.".to_string());
+        }
+
+        if !a_qubits.contains(&min_qubit) {
+            let qubit_a = a_qubits[0];
+            gate_seq.push((
+                StandardGate::SwapGate,
+                smallvec![],
+                smallvec![Qubit(min_qubit as u32), Qubit(qubit_a as u32)],
+            ));
+            self.symplectic_matrix.prepend_swap(min_qubit, qubit_a);
+
+            if b_qubits.contains(&min_qubit) {
+                b_qubits.swap_remove(&min_qubit);
+                b_qubits.insert(qubit_a);
+            } else if c_qubits.contains(&min_qubit) {
+                c_qubits.swap_remove(&min_qubit);
+                c_qubits.insert(qubit_a);
+            } else if d_qubits.contains(&min_qubit) {
+                d_qubits.swap_remove(&min_qubit);
+                d_qubits.insert(qubit_a);
+            }
+
+            a_qubits.swap_remove(&qubit_a);
+            a_qubits.insert(min_qubit);
+        }
+
+        for qubit in c_qubits {
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(min_qubit as u32), Qubit(qubit as u32)],
+            ));
+            self.symplectic_matrix.prepend_cx(min_qubit, qubit);
+        }
+
+        for qubit in d_qubits {
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(qubit as u32), Qubit(min_qubit as u32)],
+            ));
+            self.symplectic_matrix.prepend_cx(qubit, min_qubit);
+        }
+
+        if b_qubits.len() > 1 {
+            let qubit_b = b_qubits[0];
+            for qubit in &b_qubits[1..] {
+                gate_seq.push((
+                    StandardGate::CXGate,
+                    smallvec![],
+                    smallvec![Qubit(qubit_b as u32), Qubit(*qubit as u32)],
+                ));
+                self.symplectic_matrix.prepend_cx(qubit_b, *qubit);
+            }
+        }
+
+        if !b_qubits.is_empty() {
+            let qubit_b = b_qubits[0];
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(min_qubit as u32), Qubit(qubit_b as u32)],
+            ));
+            self.symplectic_matrix.prepend_cx(min_qubit, qubit_b);
+
+            gate_seq.push((
+                StandardGate::HGate,
+                smallvec![],
+                smallvec![Qubit(qubit_b as u32)],
+            ));
+            self.symplectic_matrix.prepend_h(qubit_b);
+
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(qubit_b as u32), Qubit(min_qubit as u32)],
+            ));
+            self.symplectic_matrix.prepend_cx(qubit_b, min_qubit);
+        }
+
+        let a_len: usize = (a_qubits.len() - 1) / 2;
+        if a_len > 0 {
+            a_qubits.swap_remove(&min_qubit);
+        }
+
+        for qubit in 0..a_len {
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![
+                    Qubit(a_qubits[2 * qubit + 1] as u32),
+                    Qubit(a_qubits[2 * qubit] as u32)
+                ],
+            ));
+            self.symplectic_matrix
+                .prepend_cx(a_qubits[2 * qubit + 1], a_qubits[2 * qubit]);
+
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(a_qubits[2 * qubit] as u32), Qubit(min_qubit as u32)],
+            ));
+            self.symplectic_matrix
+                .prepend_cx(a_qubits[2 * qubit], min_qubit);
+
+            gate_seq.push((
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![
+                    Qubit(min_qubit as u32),
+                    Qubit(a_qubits[2 * qubit + 1] as u32)
+                ],
+            ));
+            self.symplectic_matrix
+                .prepend_cx(min_qubit, a_qubits[2 * qubit + 1]);
+        }
+
+        Ok(())
+    }
+
+    /// The main synthesis function.
+    pub(crate) fn run(&mut self) -> Result<(usize, CliffordGatesVec), String> {
+        let mut clifford_gates = CliffordGatesVec::new();
+
+        while !self.unprocessed_qubits.is_empty() {
+            let costs: Vec<(usize, usize)> = self
+                .unprocessed_qubits
+                .iter()
+                .map(|q| self.compute_cost(*q).map(|cost| (cost, *q)))
+                .collect::<Result<Vec<_>, _>>()?;
+
+            let min_cost_qubit = costs.iter().min_by_key(|(cost, _)| cost).unwrap().1;
+
+            self.decouple_qubit(&mut clifford_gates, min_cost_qubit)?;
+
+            self.unprocessed_qubits.swap_remove(&min_cost_qubit);
+        }
+
+        adjust_final_pauli_gates(&mut clifford_gates, self.tableau, self.num_qubits)?;
+
+        Ok((self.num_qubits, clifford_gates))
+    }
+}
diff --git a/crates/accelerate/src/synthesis/clifford/mod.rs b/crates/accelerate/src/synthesis/clifford/mod.rs
new file mode 100644
index 000000000000..6772228acf88
--- /dev/null
+++ b/crates/accelerate/src/synthesis/clifford/mod.rs
@@ -0,0 +1,48 @@
+// 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.
+
+mod greedy_synthesis;
+mod utils;
+
+use crate::synthesis::clifford::greedy_synthesis::GreedyCliffordSynthesis;
+use crate::QiskitError;
+use numpy::PyReadonlyArray2;
+use pyo3::prelude::*;
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::Param;
+
+/// Create a circuit that synthesizes a given Clifford operator represented as a tableau.
+///
+/// This is an implementation of the "greedy Clifford compiler" presented in
+/// Appendix A of the paper "Clifford Circuit Optimization with Templates and Symbolic
+/// Pauli Gates" by Bravyi, Shaydulin, Hu, and Maslov (2021), `<https://arxiv.org/abs/2105.02291>`__.
+///
+/// This method typically yields better CX cost compared to the Aaronson-Gottesman method.
+///
+/// Note that this function only implements the greedy Clifford compiler and not the
+/// templates and symbolic Pauli gates optimizations that are also described in the paper.
+#[pyfunction]
+#[pyo3(signature = (clifford))]
+fn synth_clifford_greedy(py: Python, clifford: PyReadonlyArray2<bool>) -> PyResult<CircuitData> {
+    let tableau = clifford.as_array();
+    let mut greedy_synthesis =
+        GreedyCliffordSynthesis::new(tableau.view()).map_err(QiskitError::new_err)?;
+    let (num_qubits, clifford_gates) = greedy_synthesis.run().map_err(QiskitError::new_err)?;
+
+    CircuitData::from_standard_gates(py, num_qubits as u32, clifford_gates, Param::Float(0.0))
+}
+
+#[pymodule]
+pub fn clifford(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(synth_clifford_greedy, m)?)?;
+    Ok(())
+}
diff --git a/crates/accelerate/src/synthesis/clifford/utils.rs b/crates/accelerate/src/synthesis/clifford/utils.rs
new file mode 100644
index 000000000000..766d84ed179d
--- /dev/null
+++ b/crates/accelerate/src/synthesis/clifford/utils.rs
@@ -0,0 +1,289 @@
+// 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 crate::synthesis::linear::utils::calc_inverse_matrix_inner;
+use ndarray::{azip, s, Array1, Array2, ArrayView2};
+use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::Qubit;
+use smallvec::{smallvec, SmallVec};
+
+/// Symplectic matrix.
+/// Currently this class is internal to the synthesis library.
+pub struct SymplecticMatrix {
+    /// Number of qubits.
+    pub num_qubits: usize,
+    /// Matrix with dimensions (2 * num_qubits) x (2 * num_qubits).
+    pub smat: Array2<bool>,
+}
+
+/// Clifford.
+/// Currently this class is internal to the synthesis library and
+/// has a very different functionality from Qiskit's python-based
+/// Clifford class.
+pub struct Clifford {
+    /// Number of qubits.
+    pub num_qubits: usize,
+    /// Matrix with dimensions (2 * num_qubits) x (2 * num_qubits + 1).
+    pub tableau: Array2<bool>,
+}
+
+impl SymplecticMatrix {
+    /// Modifies the matrix in-place by appending S-gate
+    #[allow(dead_code)]
+    pub fn append_s(&mut self, qubit: usize) {
+        let (x, mut z) = self
+            .smat
+            .multi_slice_mut((s![.., qubit], s![.., self.num_qubits + qubit]));
+        azip!((z in &mut z, &x in &x) *z ^= x);
+    }
+
+    /// Modifies the matrix in-place by prepending S-gate
+    pub fn prepend_s(&mut self, qubit: usize) {
+        let (x, mut z) = self
+            .smat
+            .multi_slice_mut((s![self.num_qubits + qubit, ..], s![qubit, ..]));
+        azip!((z in &mut z, &x in &x) *z ^= x);
+    }
+
+    /// Modifies the matrix in-place by appending H-gate
+    #[allow(dead_code)]
+    pub fn append_h(&mut self, qubit: usize) {
+        let (mut x, mut z) = self
+            .smat
+            .multi_slice_mut((s![.., qubit], s![.., self.num_qubits + qubit]));
+        azip!((x in &mut x, z in &mut z)  (*x, *z) = (*z, *x));
+    }
+
+    /// Modifies the matrix in-place by prepending H-gate
+    pub fn prepend_h(&mut self, qubit: usize) {
+        let (mut x, mut z) = self
+            .smat
+            .multi_slice_mut((s![qubit, ..], s![self.num_qubits + qubit, ..]));
+        azip!((x in &mut x, z in &mut z)  (*x, *z) = (*z, *x));
+    }
+
+    /// Modifies the matrix in-place by appending SWAP-gate
+    #[allow(dead_code)]
+    pub fn append_swap(&mut self, qubit0: usize, qubit1: usize) {
+        let (mut x0, mut z0, mut x1, mut z1) = self.smat.multi_slice_mut((
+            s![.., qubit0],
+            s![.., self.num_qubits + qubit0],
+            s![.., qubit1],
+            s![.., self.num_qubits + qubit1],
+        ));
+        azip!((x0 in &mut x0, x1 in &mut x1)  (*x0, *x1) = (*x1, *x0));
+        azip!((z0 in &mut z0, z1 in &mut z1)  (*z0, *z1) = (*z1, *z0));
+    }
+
+    /// Modifies the matrix in-place by prepending SWAP-gate
+    pub fn prepend_swap(&mut self, qubit0: usize, qubit1: usize) {
+        let (mut x0, mut z0, mut x1, mut z1) = self.smat.multi_slice_mut((
+            s![qubit0, ..],
+            s![self.num_qubits + qubit0, ..],
+            s![qubit1, ..],
+            s![self.num_qubits + qubit1, ..],
+        ));
+        azip!((x0 in &mut x0, x1 in &mut x1)  (*x0, *x1) = (*x1, *x0));
+        azip!((z0 in &mut z0, z1 in &mut z1)  (*z0, *z1) = (*z1, *z0));
+    }
+
+    /// Modifies the matrix in-place by appending CX-gate
+    #[allow(dead_code)]
+    pub fn append_cx(&mut self, qubit0: usize, qubit1: usize) {
+        let (x0, mut z0, mut x1, z1) = self.smat.multi_slice_mut((
+            s![.., qubit0],
+            s![.., self.num_qubits + qubit0],
+            s![.., qubit1],
+            s![.., self.num_qubits + qubit1],
+        ));
+        azip!((x1 in &mut x1, &x0 in &x0) *x1 ^= x0);
+        azip!((z0 in &mut z0, &z1 in &z1) *z0 ^= z1);
+    }
+
+    /// Modifies the matrix in-place by prepending CX-gate
+    pub fn prepend_cx(&mut self, qubit0: usize, qubit1: usize) {
+        let (x0, mut z0, mut x1, z1) = self.smat.multi_slice_mut((
+            s![qubit1, ..],
+            s![self.num_qubits + qubit1, ..],
+            s![qubit0, ..],
+            s![self.num_qubits + qubit0, ..],
+        ));
+        azip!((x1 in &mut x1, &x0 in &x0) *x1 ^= x0);
+        azip!((z0 in &mut z0, &z1 in &z1) *z0 ^= z1);
+    }
+}
+
+impl Clifford {
+    /// Modifies the tableau in-place by appending S-gate
+    pub fn append_s(&mut self, qubit: usize) {
+        let (x, mut z, mut p) = self.tableau.multi_slice_mut((
+            s![.., qubit],
+            s![.., self.num_qubits + qubit],
+            s![.., 2 * self.num_qubits],
+        ));
+
+        azip!((p in &mut p, &x in &x, &z in &z)  *p ^= x & z);
+        azip!((z in &mut z, &x in &x) *z ^= x);
+    }
+
+    /// Modifies the tableau in-place by appending Sdg-gate
+    #[allow(dead_code)]
+    pub fn append_sdg(&mut self, qubit: usize) {
+        let (x, mut z, mut p) = self.tableau.multi_slice_mut((
+            s![.., qubit],
+            s![.., self.num_qubits + qubit],
+            s![.., 2 * self.num_qubits],
+        ));
+
+        azip!((p in &mut p, &x in &x, &z in &z)  *p ^= x & !z);
+        azip!((z in &mut z, &x in &x) *z ^= x);
+    }
+
+    /// Modifies the tableau in-place by appending H-gate
+    pub fn append_h(&mut self, qubit: usize) {
+        let (mut x, mut z, mut p) = self.tableau.multi_slice_mut((
+            s![.., qubit],
+            s![.., self.num_qubits + qubit],
+            s![.., 2 * self.num_qubits],
+        ));
+
+        azip!((p in &mut p, &x in &x, &z in &z)  *p ^= x & z);
+        azip!((x in &mut x, z in &mut z)  (*x, *z) = (*z, *x));
+    }
+
+    /// Modifies the tableau in-place by appending SWAP-gate
+    pub fn append_swap(&mut self, qubit0: usize, qubit1: usize) {
+        let (mut x0, mut z0, mut x1, mut z1) = self.tableau.multi_slice_mut((
+            s![.., qubit0],
+            s![.., self.num_qubits + qubit0],
+            s![.., qubit1],
+            s![.., self.num_qubits + qubit1],
+        ));
+        azip!((x0 in &mut x0, x1 in &mut x1)  (*x0, *x1) = (*x1, *x0));
+        azip!((z0 in &mut z0, z1 in &mut z1)  (*z0, *z1) = (*z1, *z0));
+    }
+
+    /// Modifies the tableau in-place by appending CX-gate
+    pub fn append_cx(&mut self, qubit0: usize, qubit1: usize) {
+        let (x0, mut z0, mut x1, z1, mut p) = self.tableau.multi_slice_mut((
+            s![.., qubit0],
+            s![.., self.num_qubits + qubit0],
+            s![.., qubit1],
+            s![.., self.num_qubits + qubit1],
+            s![.., 2 * self.num_qubits],
+        ));
+        azip!((p in &mut p, &x0 in &x0, &z0 in &z0, &x1 in &x1, &z1 in &z1) *p ^= (x1 ^ z0 ^ true) & z1 & x0);
+        azip!((x1 in &mut x1, &x0 in &x0) *x1 ^= x0);
+        azip!((z0 in &mut z0, &z1 in &z1) *z0 ^= z1);
+    }
+
+    /// Creates a Clifford from a given sequence of Clifford gates.
+    /// In essence, starts from the identity tableau and modifies it
+    /// based on the gates in the sequence.
+    pub fn from_gate_sequence(
+        gate_seq: &CliffordGatesVec,
+        num_qubits: usize,
+    ) -> Result<Clifford, String> {
+        // create the identity
+        let mut clifford = Clifford {
+            num_qubits,
+            tableau: Array2::from_shape_fn((2 * num_qubits, 2 * num_qubits + 1), |(i, j)| i == j),
+        };
+
+        gate_seq
+            .iter()
+            .try_for_each(|(gate, _params, qubits)| match *gate {
+                StandardGate::SGate => {
+                    clifford.append_s(qubits[0].0 as usize);
+                    Ok(())
+                }
+                StandardGate::HGate => {
+                    clifford.append_h(qubits[0].0 as usize);
+                    Ok(())
+                }
+                StandardGate::CXGate => {
+                    clifford.append_cx(qubits[0].0 as usize, qubits[1].0 as usize);
+                    Ok(())
+                }
+                StandardGate::SwapGate => {
+                    clifford.append_swap(qubits[0].0 as usize, qubits[1].0 as usize);
+                    Ok(())
+                }
+                _ => Err(format!("Unsupported gate {:?}", gate)),
+            })?;
+        Ok(clifford)
+    }
+}
+
+/// A sequence of Clifford gates.
+/// Represents the return type of Clifford synthesis algorithms.
+pub type CliffordGatesVec = Vec<(StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>)>;
+
+/// Given a sequence of Clifford gates that correctly implements the symplectic matrix
+/// of the target clifford tableau, adds the Pauli gates to also match the phase of
+/// the tableau.
+pub fn adjust_final_pauli_gates(
+    gate_seq: &mut CliffordGatesVec,
+    target_tableau: ArrayView2<bool>,
+    num_qubits: usize,
+) -> Result<(), String> {
+    // simulate the clifford circuit that we have constructed
+    let simulated_clifford = Clifford::from_gate_sequence(gate_seq, num_qubits)?;
+
+    // compute the phase difference
+    let target_phase = target_tableau.column(2 * num_qubits);
+    let sim_phase = simulated_clifford.tableau.column(2 * num_qubits);
+
+    let delta_phase: Vec<bool> = target_phase
+        .iter()
+        .zip(sim_phase.iter())
+        .map(|(&a, &b)| a ^ b)
+        .collect();
+
+    // compute inverse of the symplectic matrix
+    let smat = target_tableau.slice(s![.., ..2 * num_qubits]);
+    let smat_inv = calc_inverse_matrix_inner(smat, false)?;
+
+    // compute smat_inv * delta_phase
+    let arr1 = smat_inv.map(|v| *v as usize);
+    let vec2: Vec<usize> = delta_phase.into_iter().map(|v| v as usize).collect();
+    let arr2 = Array1::from(vec2);
+    let delta_phase_pre = arr1.dot(&arr2).map(|v| v % 2 == 1);
+
+    // add pauli gates
+    for qubit in 0..num_qubits {
+        if delta_phase_pre[qubit] && delta_phase_pre[qubit + num_qubits] {
+            // println!("=> Adding Y-gate on {}", qubit);
+            gate_seq.push((
+                StandardGate::YGate,
+                smallvec![],
+                smallvec![Qubit(qubit as u32)],
+            ));
+        } else if delta_phase_pre[qubit] {
+            // println!("=> Adding Z-gate on {}", qubit);
+            gate_seq.push((
+                StandardGate::ZGate,
+                smallvec![],
+                smallvec![Qubit(qubit as u32)],
+            ));
+        } else if delta_phase_pre[qubit + num_qubits] {
+            // println!("=> Adding X-gate on {}", qubit);
+            gate_seq.push((
+                StandardGate::XGate,
+                smallvec![],
+                smallvec![Qubit(qubit as u32)],
+            ));
+        }
+    }
+
+    Ok(())
+}
diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs
index 2fa158ea761f..49dfeefd3869 100644
--- a/crates/accelerate/src/synthesis/linear/mod.rs
+++ b/crates/accelerate/src/synthesis/linear/mod.rs
@@ -14,7 +14,8 @@ use crate::QiskitError;
 use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2};
 use pyo3::prelude::*;
 
-mod utils;
+mod pmh;
+pub mod utils;
 
 #[pyfunction]
 #[pyo3(signature = (mat, ncols=None, full_elim=false))]
@@ -186,5 +187,6 @@ pub fn linear(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(binary_matmul))?;
     m.add_wrapped(wrap_pyfunction!(random_invertible_binary_matrix))?;
     m.add_wrapped(wrap_pyfunction!(check_invertible_binary_matrix))?;
+    m.add_wrapped(wrap_pyfunction!(pmh::synth_cnot_count_full_pmh))?;
     Ok(())
 }
diff --git a/crates/accelerate/src/synthesis/linear/pmh.rs b/crates/accelerate/src/synthesis/linear/pmh.rs
new file mode 100644
index 000000000000..d7283fd580e4
--- /dev/null
+++ b/crates/accelerate/src/synthesis/linear/pmh.rs
@@ -0,0 +1,195 @@
+// 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 hashbrown::HashMap;
+use ndarray::{s, Array1, Array2, ArrayViewMut2, Axis};
+use numpy::PyReadonlyArray2;
+use smallvec::smallvec;
+use std::cmp;
+
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::Qubit;
+
+use pyo3::prelude::*;
+
+use super::utils::_add_row_or_col;
+
+/// This helper function allows transposed access to a matrix.
+fn _index(transpose: bool, i: usize, j: usize) -> (usize, usize) {
+    if transpose {
+        (j, i)
+    } else {
+        (i, j)
+    }
+}
+
+fn _ceil_fraction(numerator: usize, denominator: usize) -> usize {
+    let mut fraction = numerator / denominator;
+    if numerator % denominator > 0 {
+        fraction += 1;
+    }
+    fraction
+}
+
+/// This function is a helper function of the algorithm for optimal synthesis
+/// of linear reversible circuits (the Patel–Markov–Hayes algorithm). It works
+/// like gaussian elimination, except that it works a lot faster, and requires
+/// fewer steps (and therefore fewer CNOTs). It takes the matrix and
+/// splits it into sections of size section_size. Then it eliminates all non-zero
+/// sub-rows within each section, which are the same as a non-zero sub-row
+/// above. Once this has been done, it continues with normal gaussian elimination.
+/// The benefit is that with small section sizes, most of the sub-rows will
+/// be cleared in the first step, resulting in a factor ``section_size`` fewer row row operations
+/// during Gaussian elimination.
+///
+/// The algorithm is described in detail in the following paper
+/// "Optimal synthesis of linear reversible circuits."
+/// Patel, Ketan N., Igor L. Markov, and John P. Hayes.
+/// Quantum Information & Computation 8.3 (2008): 282-294.
+///
+/// Note:
+/// This implementation tweaks the Patel, Markov, and Hayes algorithm by adding
+/// a "back reduce" which adds rows below the pivot row with a high degree of
+/// overlap back to it. The intuition is to avoid a high-weight pivot row
+/// increasing the weight of lower rows.
+///
+/// Args:
+///     matrix: square matrix, describing a linear quantum circuit
+///     section_size: the section size the matrix columns are divided into
+///
+/// Returns:
+///     A vector of CX locations (control, target) that need to be applied.
+fn lower_cnot_synth(
+    mut matrix: ArrayViewMut2<bool>,
+    section_size: usize,
+    transpose: bool,
+) -> Vec<(usize, usize)> {
+    // The vector of CNOTs to be applied. Called ``circuit`` here for consistency with the paper.
+    let mut circuit: Vec<(usize, usize)> = Vec::new();
+    let cutoff = 1;
+
+    // to apply to the transposed matrix, we can just set axis = 1
+    let row_axis = if transpose { Axis(1) } else { Axis(0) };
+
+    // get number of columns (same as rows) and the number of sections
+    let n = matrix.raw_dim()[0];
+    let num_sections = _ceil_fraction(n, section_size);
+
+    // iterate over the columns
+    for section in 1..num_sections + 1 {
+        // store sub section row patterns here, which we saw already
+        let mut patterns: HashMap<Array1<bool>, usize> = HashMap::new();
+        let section_slice = s![(section - 1) * section_size..cmp::min(section * section_size, n)];
+
+        // iterate over the rows (note we only iterate from the diagonal downwards)
+        for row_idx in (section - 1) * section_size..n {
+            // we need to keep track of the rows we saw already, called ``pattern`` here
+            let pattern: Array1<bool> = matrix
+                .index_axis(row_axis, row_idx)
+                .slice(section_slice)
+                .to_owned();
+
+            // skip if the row is empty (i.e. all elements are false)
+            if pattern.iter().any(|&el| el) {
+                if patterns.contains_key(&pattern) {
+                    // store CX location
+                    circuit.push((patterns[&pattern], row_idx));
+                    // remove the row
+                    _add_row_or_col(matrix.view_mut(), &transpose, patterns[&pattern], row_idx);
+                } else {
+                    // if we have not seen this pattern yet, keep track of it
+                    patterns.insert(pattern, row_idx);
+                }
+            }
+        }
+
+        // gaussian eliminate the remainder of the section
+        for col_idx in (section - 1) * section_size..cmp::min(section * section_size, n) {
+            let mut diag_el = matrix[[col_idx, col_idx]];
+
+            for r in col_idx + 1..n {
+                if matrix[_index(transpose, r, col_idx)] {
+                    if !diag_el {
+                        _add_row_or_col(matrix.view_mut(), &transpose, r, col_idx);
+                        circuit.push((r, col_idx));
+                        diag_el = true
+                    }
+                    _add_row_or_col(matrix.view_mut(), &transpose, col_idx, r);
+                    circuit.push((col_idx, r));
+                }
+
+                // back-reduce to the pivot row: this one-line-magic checks if the logical AND
+                // between the two target rows has more ``true`` elements is larger than the cutoff
+                if matrix
+                    .index_axis(row_axis, col_idx)
+                    .iter()
+                    .zip(matrix.index_axis(row_axis, r).iter())
+                    .map(|(&i, &j)| i & j)
+                    .filter(|&x| x)
+                    .count()
+                    > cutoff
+                {
+                    _add_row_or_col(matrix.view_mut(), &transpose, r, col_idx);
+                    circuit.push((r, col_idx));
+                }
+            }
+        }
+    }
+    circuit
+}
+
+/// Synthesize a linear function, given by a boolean square matrix, into a circuit.
+/// This function uses the Patel-Markov-Hayes algorithm, described in arXiv:quant-ph/0302002,
+/// using section-wise elimination of the rows.
+#[pyfunction]
+#[pyo3(signature = (matrix, section_size=None))]
+pub fn synth_cnot_count_full_pmh(
+    py: Python,
+    matrix: PyReadonlyArray2<bool>,
+    section_size: Option<i64>,
+) -> PyResult<CircuitData> {
+    let arrayview = matrix.as_array();
+    let mut mat: Array2<bool> = arrayview.to_owned();
+    let num_qubits = mat.nrows(); // is a quadratic matrix
+
+    // If given, use the user-specified input size. If None, we default to
+    // alpha * log2(num_qubits), as suggested in the paper, where the coefficient alpha
+    // is calibrated to minimize the upper bound on the number of row operations.
+    // In addition, we at least set a block size of 2, which, in practice, minimizes the bound
+    // until ~100 qubits.
+    let alpha = 0.56;
+    let blocksize = match section_size {
+        Some(section_size) => section_size as usize,
+        None => std::cmp::max(2, (alpha * (num_qubits as f64).log2()).floor() as usize),
+    };
+
+    // compute the synthesis for the lower triangular part of the matrix, and then
+    // apply it on the transposed part for the full synthesis
+    let lower_cnots = lower_cnot_synth(mat.view_mut(), blocksize, false);
+    let upper_cnots = lower_cnot_synth(mat.view_mut(), blocksize, true);
+
+    // iterator over the gates
+    let instructions = upper_cnots
+        .iter()
+        .map(|(i, j)| (*j, *i))
+        .chain(lower_cnots.into_iter().rev())
+        .map(|(ctrl, target)| {
+            (
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(ctrl as u32), Qubit(target as u32)],
+            )
+        });
+
+    CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0))
+}
diff --git a/crates/accelerate/src/synthesis/mod.rs b/crates/accelerate/src/synthesis/mod.rs
index db28751437f6..1b9908ef80cf 100644
--- a/crates/accelerate/src/synthesis/mod.rs
+++ b/crates/accelerate/src/synthesis/mod.rs
@@ -10,7 +10,8 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-mod linear;
+mod clifford;
+pub mod linear;
 mod permutation;
 
 use pyo3::prelude::*;
@@ -18,7 +19,8 @@ use pyo3::wrap_pymodule;
 
 #[pymodule]
 pub fn synthesis(m: &Bound<PyModule>) -> PyResult<()> {
-    m.add_wrapped(wrap_pymodule!(permutation::permutation))?;
     m.add_wrapped(wrap_pymodule!(linear::linear))?;
+    m.add_wrapped(wrap_pymodule!(permutation::permutation))?;
+    m.add_wrapped(wrap_pymodule!(clifford::clifford))?;
     Ok(())
 }
diff --git a/crates/accelerate/src/synthesis/permutation/utils.rs b/crates/accelerate/src/synthesis/permutation/utils.rs
index 620ce4628741..47a9e1c3a7a9 100644
--- a/crates/accelerate/src/synthesis/permutation/utils.rs
+++ b/crates/accelerate/src/synthesis/permutation/utils.rs
@@ -120,8 +120,10 @@ pub fn pattern_to_cycles(pattern: &ArrayView1<usize>) -> Vec<Vec<usize>> {
 /// Periodic (or Python-like) access to a vector.
 /// Util used below in ``decompose_cycles``.
 #[inline]
-fn pget(vec: &Vec<usize>, index: isize) -> Result<usize, PySequenceIndexError> {
-    let SequenceIndex::Int(wrapped) = PySequenceIndex::Int(index).with_len(vec.len())? else {unreachable!()};
+fn pget(vec: &[usize], index: isize) -> Result<usize, PySequenceIndexError> {
+    let SequenceIndex::Int(wrapped) = PySequenceIndex::Int(index).with_len(vec.len())? else {
+        unreachable!()
+    };
     Ok(vec[wrapped])
 }
 
diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs
index 37061d5159f4..568206925c2b 100644
--- a/crates/accelerate/src/two_qubit_decompose.rs
+++ b/crates/accelerate/src/two_qubit_decompose.rs
@@ -52,6 +52,7 @@ use rand_distr::StandardNormal;
 use rand_pcg::Pcg64Mcg;
 
 use qiskit_circuit::gate_matrix::{CX_GATE, H_GATE, ONE_QUBIT_IDENTITY, SX_GATE, X_GATE};
+use qiskit_circuit::operations::Operation;
 use qiskit_circuit::slice::{PySequenceIndex, SequenceIndex};
 use qiskit_circuit::util::{c64, GateArray1Q, GateArray2Q, C_M_ONE, C_ONE, C_ZERO, IM, M_IM};
 
@@ -1045,7 +1046,7 @@ impl TwoQubitWeylDecomposition {
         )
         .unwrap();
         for gate in c2r.gates {
-            gate_sequence.push((gate.0, gate.1, smallvec![0]))
+            gate_sequence.push((gate.0.name().to_string(), gate.1, smallvec![0]))
         }
         global_phase += c2r.global_phase;
         let c2l = unitary_to_gate_sequence_inner(
@@ -1058,7 +1059,7 @@ impl TwoQubitWeylDecomposition {
         )
         .unwrap();
         for gate in c2l.gates {
-            gate_sequence.push((gate.0, gate.1, smallvec![1]))
+            gate_sequence.push((gate.0.name().to_string(), gate.1, smallvec![1]))
         }
         global_phase += c2l.global_phase;
         self.weyl_gate(
@@ -1077,7 +1078,7 @@ impl TwoQubitWeylDecomposition {
         )
         .unwrap();
         for gate in c1r.gates {
-            gate_sequence.push((gate.0, gate.1, smallvec![0]))
+            gate_sequence.push((gate.0.name().to_string(), gate.1, smallvec![0]))
         }
         global_phase += c2r.global_phase;
         let c1l = unitary_to_gate_sequence_inner(
@@ -1090,7 +1091,7 @@ impl TwoQubitWeylDecomposition {
         )
         .unwrap();
         for gate in c1l.gates {
-            gate_sequence.push((gate.0, gate.1, smallvec![1]))
+            gate_sequence.push((gate.0.name().to_string(), gate.1, smallvec![1]))
         }
         Ok(TwoQubitGateSequence {
             gates: gate_sequence,
@@ -1459,7 +1460,7 @@ impl TwoQubitBasisDecomposer {
         if let Some(sequence) = sequence {
             *global_phase += sequence.global_phase;
             for gate in sequence.gates {
-                gates.push((gate.0, gate.1, smallvec![qubit]));
+                gates.push((gate.0.name().to_string(), gate.1, smallvec![qubit]));
             }
         }
     }
@@ -1847,13 +1848,13 @@ impl TwoQubitBasisDecomposer {
         for i in 0..best_nbasis as usize {
             if let Some(euler_decomp) = &euler_decompositions[2 * i] {
                 for gate in &euler_decomp.gates {
-                    gates.push((gate.0.clone(), gate.1.clone(), smallvec![0]));
+                    gates.push((gate.0.name().to_string(), gate.1.clone(), smallvec![0]));
                 }
                 global_phase += euler_decomp.global_phase
             }
             if let Some(euler_decomp) = &euler_decompositions[2 * i + 1] {
                 for gate in &euler_decomp.gates {
-                    gates.push((gate.0.clone(), gate.1.clone(), smallvec![1]));
+                    gates.push((gate.0.name().to_string(), gate.1.clone(), smallvec![1]));
                 }
                 global_phase += euler_decomp.global_phase
             }
@@ -1861,13 +1862,13 @@ impl TwoQubitBasisDecomposer {
         }
         if let Some(euler_decomp) = &euler_decompositions[2 * best_nbasis as usize] {
             for gate in &euler_decomp.gates {
-                gates.push((gate.0.clone(), gate.1.clone(), smallvec![0]));
+                gates.push((gate.0.name().to_string(), gate.1.clone(), smallvec![0]));
             }
             global_phase += euler_decomp.global_phase
         }
         if let Some(euler_decomp) = &euler_decompositions[2 * best_nbasis as usize + 1] {
             for gate in &euler_decomp.gates {
-                gates.push((gate.0.clone(), gate.1.clone(), smallvec![1]));
+                gates.push((gate.0.name().to_string(), gate.1.clone(), smallvec![1]));
             }
             global_phase += euler_decomp.global_phase
         }
diff --git a/crates/circuit/src/circuit_data.rs b/crates/circuit/src/circuit_data.rs
index f10911cc440f..501645415874 100644
--- a/crates/circuit/src/circuit_data.rs
+++ b/crates/circuit/src/circuit_data.rs
@@ -1048,13 +1048,7 @@ impl CircuitData {
         Ok(PySet::new_bound(py, self.param_table.uuid_map.values())?.unbind())
     }
 
-    pub fn pop_param(
-        &mut self,
-        py: Python,
-        uuid: u128,
-        name: String,
-        default: PyObject,
-    ) -> PyObject {
+    pub fn pop_param(&mut self, py: Python, uuid: u128, name: &str, default: PyObject) -> PyObject {
         match self.param_table.pop(uuid, name) {
             Some(res) => res.into_py(py),
             None => default.clone_ref(py),
diff --git a/crates/circuit/src/circuit_instruction.rs b/crates/circuit/src/circuit_instruction.rs
index d6516722fbac..ed1c358cbc5b 100644
--- a/crates/circuit/src/circuit_instruction.rs
+++ b/crates/circuit/src/circuit_instruction.rs
@@ -460,6 +460,12 @@ impl CircuitInstruction {
             .and_then(|attrs| attrs.unit.as_deref())
     }
 
+    pub fn is_parameterized(&self) -> bool {
+        self.params
+            .iter()
+            .any(|x| matches!(x, Param::ParameterExpression(_)))
+    }
+
     /// Creates a shallow copy with the given fields replaced.
     ///
     /// Returns:
diff --git a/crates/circuit/src/dag_node.rs b/crates/circuit/src/dag_node.rs
index ffd7920a36fd..55a40c83dc39 100644
--- a/crates/circuit/src/dag_node.rs
+++ b/crates/circuit/src/dag_node.rs
@@ -14,11 +14,12 @@ use crate::circuit_instruction::{
     convert_py_to_operation_type, operation_type_to_py, CircuitInstruction,
     ExtraInstructionAttributes,
 };
+use crate::imports::QUANTUM_CIRCUIT;
 use crate::operations::Operation;
 use numpy::IntoPyArray;
 use pyo3::prelude::*;
 use pyo3::types::{PyDict, PyList, PySequence, PyString, PyTuple};
-use pyo3::{intern, IntoPy, PyObject, PyResult};
+use pyo3::{intern, IntoPy, PyObject, PyResult, ToPyObject};
 use smallvec::smallvec;
 
 /// Parent class for DAGOpNode, DAGInNode, and DAGOutNode.
@@ -135,6 +136,50 @@ impl DAGOpNode {
         ))
     }
 
+    #[staticmethod]
+    fn from_instruction(
+        py: Python,
+        instruction: CircuitInstruction,
+        dag: Option<&Bound<PyAny>>,
+    ) -> PyResult<PyObject> {
+        let qargs = instruction.qubits.clone_ref(py).into_bound(py);
+        let cargs = instruction.clbits.clone_ref(py).into_bound(py);
+
+        let sort_key = match dag {
+            Some(dag) => {
+                let cache = dag
+                    .getattr(intern!(py, "_key_cache"))?
+                    .downcast_into_exact::<PyDict>()?;
+                let cache_key = PyTuple::new_bound(py, [&qargs, &cargs]);
+                match cache.get_item(&cache_key)? {
+                    Some(key) => key,
+                    None => {
+                        let indices: PyResult<Vec<_>> = qargs
+                            .iter()
+                            .chain(cargs.iter())
+                            .map(|bit| {
+                                dag.call_method1(intern!(py, "find_bit"), (bit,))?
+                                    .getattr(intern!(py, "index"))
+                            })
+                            .collect();
+                        let index_strs: Vec<_> =
+                            indices?.into_iter().map(|i| format!("{:04}", i)).collect();
+                        let key = PyString::new_bound(py, index_strs.join(",").as_str());
+                        cache.set_item(&cache_key, &key)?;
+                        key.into_any()
+                    }
+                }
+            }
+            None => qargs.str()?.into_any(),
+        };
+        let base = PyClassInitializer::from(DAGNode { _node_id: -1 });
+        let sub = base.add_subclass(DAGOpNode {
+            instruction,
+            sort_key: sort_key.unbind(),
+        });
+        Ok(Py::new(py, sub)?.to_object(py))
+    }
+
     fn __reduce__(slf: PyRef<Self>, py: Python) -> PyResult<PyObject> {
         let state = (slf.as_ref()._node_id, &slf.sort_key);
         Ok((
@@ -184,6 +229,16 @@ impl DAGOpNode {
         Ok(())
     }
 
+    #[getter]
+    fn num_qubits(&self) -> u32 {
+        self.instruction.operation.num_qubits()
+    }
+
+    #[getter]
+    fn num_clbits(&self) -> u32 {
+        self.instruction.operation.num_clbits()
+    }
+
     #[getter]
     fn get_qargs(&self, py: Python) -> Py<PyTuple> {
         self.instruction.qubits.clone_ref(py)
@@ -206,8 +261,8 @@ impl DAGOpNode {
 
     /// Returns the Instruction name corresponding to the op for this node
     #[getter]
-    fn get_name(&self, py: Python) -> PyObject {
-        self.instruction.operation.name().to_object(py)
+    fn get_name(&self) -> &str {
+        self.instruction.operation.name()
     }
 
     #[getter]
@@ -215,6 +270,10 @@ impl DAGOpNode {
         self.instruction.params.to_object(py)
     }
 
+    pub fn is_parameterized(&self) -> bool {
+        self.instruction.is_parameterized()
+    }
+
     #[getter]
     fn matrix(&self, py: Python) -> Option<PyObject> {
         let matrix = self.instruction.operation.matrix(&self.instruction.params);
@@ -281,6 +340,21 @@ impl DAGOpNode {
         }
     }
 
+    #[getter]
+    fn definition<'py>(&self, py: Python<'py>) -> PyResult<Option<Bound<'py, PyAny>>> {
+        let definition = self
+            .instruction
+            .operation
+            .definition(&self.instruction.params);
+        definition
+            .map(|data| {
+                QUANTUM_CIRCUIT
+                    .get_bound(py)
+                    .call_method1(intern!(py, "_from_circuit_data"), (data,))
+            })
+            .transpose()
+    }
+
     /// Sets the Instruction name corresponding to the op for this node
     #[setter]
     fn set_name(&mut self, py: Python, new_name: PyObject) -> PyResult<()> {
diff --git a/crates/circuit/src/operations.rs b/crates/circuit/src/operations.rs
index c752c0af494a..6cb7ed7893bf 100644
--- a/crates/circuit/src/operations.rs
+++ b/crates/circuit/src/operations.rs
@@ -279,6 +279,12 @@ pub enum StandardGate {
     RZXGate = 52,
 }
 
+impl ToPyObject for StandardGate {
+    fn to_object(&self, py: Python) -> PyObject {
+        self.into_py(py)
+    }
+}
+
 // TODO: replace all 34s (placeholders) with actual number
 static STANDARD_GATE_NUM_QUBITS: [u32; STANDARD_GATE_SIZE] = [
     1, 1, 1, 2, 2, 2, 3, 1, 1, 1, // 0-9
@@ -746,8 +752,38 @@ impl Operation for StandardGate {
                     .expect("Unexpected Qiskit python bug"),
                 )
             }),
-            Self::RXGate => todo!("Add when we have R"),
-            Self::RYGate => todo!("Add when we have R"),
+            Self::RXGate => Python::with_gil(|py| -> Option<CircuitData> {
+                let theta = &params[0];
+                Some(
+                    CircuitData::from_standard_gates(
+                        py,
+                        1,
+                        [(
+                            Self::RGate,
+                            smallvec![theta.clone(), FLOAT_ZERO],
+                            smallvec![Qubit(0)],
+                        )],
+                        FLOAT_ZERO,
+                    )
+                    .expect("Unexpected Qiskit python bug"),
+                )
+            }),
+            Self::RYGate => Python::with_gil(|py| -> Option<CircuitData> {
+                let theta = &params[0];
+                Some(
+                    CircuitData::from_standard_gates(
+                        py,
+                        1,
+                        [(
+                            Self::RGate,
+                            smallvec![theta.clone(), Param::Float(PI / 2.0)],
+                            smallvec![Qubit(0)],
+                        )],
+                        FLOAT_ZERO,
+                    )
+                    .expect("Unexpected Qiskit python bug"),
+                )
+            }),
             Self::RZGate => Python::with_gil(|py| -> Option<CircuitData> {
                 let theta = &params[0];
                 Some(
diff --git a/crates/circuit/src/parameter_table.rs b/crates/circuit/src/parameter_table.rs
index 48c779eed3a3..9e5b31245391 100644
--- a/crates/circuit/src/parameter_table.rs
+++ b/crates/circuit/src/parameter_table.rs
@@ -153,8 +153,8 @@ impl ParamTable {
         self.uuid_map.clear();
     }
 
-    pub fn pop(&mut self, key: u128, name: String) -> Option<ParamEntry> {
-        self.names.remove(&name);
+    pub fn pop(&mut self, key: u128, name: &str) -> Option<ParamEntry> {
+        self.names.remove(name);
         self.uuid_map.remove(&key);
         self.table.remove(&key)
     }
diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 1913f8bfa8b3..4c0399a2bac8 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -83,6 +83,7 @@
 sys.modules["qiskit._accelerate.vf2_layout"] = _accelerate.vf2_layout
 sys.modules["qiskit._accelerate.synthesis.permutation"] = _accelerate.synthesis.permutation
 sys.modules["qiskit._accelerate.synthesis.linear"] = _accelerate.synthesis.linear
+sys.modules["qiskit._accelerate.synthesis.clifford"] = _accelerate.synthesis.clifford
 
 from qiskit.exceptions import QiskitError, MissingOptionalLibraryError
 
diff --git a/qiskit/circuit/classicalfunction/__init__.py b/qiskit/circuit/classicalfunction/__init__.py
index a072d910f97a..b045227b167e 100644
--- a/qiskit/circuit/classicalfunction/__init__.py
+++ b/qiskit/circuit/classicalfunction/__init__.py
@@ -51,6 +51,14 @@ def grover_oracle(a: Int1, b: Int1, c: Int1, d: Int1) -> Int1:
 Following Qiskit's little-endian bit ordering convention, the left-most bit (``a``) is the most
 significant bit and the right-most bit (``d``) is the least significant bit.
 
+.. warning::
+
+    The functionality of `qiskit.circuit.classicalfunction` requires `tweedledum`, 
+    which isn't available on all platforms (up to Python version 3.11).
+    See `tweedledum installation guide 
+    <https://github.com/boschmitt/tweedledum/tree/master?tab=readme-ov-file#installation>`_ 
+    for more details.
+
 Supplementary Information
 =========================
 
diff --git a/qiskit/circuit/commutation_checker.py b/qiskit/circuit/commutation_checker.py
index e758674829f8..79f04a65714d 100644
--- a/qiskit/circuit/commutation_checker.py
+++ b/qiskit/circuit/commutation_checker.py
@@ -21,6 +21,7 @@
 from qiskit.circuit.operation import Operation
 from qiskit.circuit.controlflow import CONTROL_FLOW_OP_NAMES
 from qiskit.quantum_info.operators import Operator
+from qiskit._accelerate.circuit import StandardGate
 
 _skipped_op_names = {"measure", "reset", "delay", "initialize"}
 _no_cache_op_names = {"annotated"}
@@ -57,6 +58,23 @@ def __init__(self, standard_gate_commutations: dict = None, cache_max_entries: i
         self._cache_miss = 0
         self._cache_hit = 0
 
+    def commute_nodes(
+        self,
+        op1,
+        op2,
+        max_num_qubits: int = 3,
+    ) -> bool:
+        """Checks if two DAGOpNodes commute."""
+        qargs1 = op1.qargs
+        cargs1 = op2.cargs
+        if not isinstance(op1._raw_op, StandardGate):
+            op1 = op1.op
+        qargs2 = op2.qargs
+        cargs2 = op2.cargs
+        if not isinstance(op2._raw_op, StandardGate):
+            op2 = op2.op
+        return self.commute(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits)
+
     def commute(
         self,
         op1: Operation,
@@ -255,9 +273,15 @@ def is_commutation_skipped(op, qargs, max_num_qubits):
     if getattr(op, "is_parameterized", False) and op.is_parameterized():
         return True
 
+    from qiskit.dagcircuit.dagnode import DAGOpNode
+
     # we can proceed if op has defined: to_operator, to_matrix and __array__, or if its definition can be
     # recursively resolved by operations that have a matrix. We check this by constructing an Operator.
-    if (hasattr(op, "to_matrix") and hasattr(op, "__array__")) or hasattr(op, "to_operator"):
+    if (
+        isinstance(op, DAGOpNode)
+        or (hasattr(op, "to_matrix") and hasattr(op, "__array__"))
+        or hasattr(op, "to_operator")
+    ):
         return False
 
     return False
@@ -409,6 +433,15 @@ def _commute_matmul(
     first_qarg = tuple(qarg[q] for q in first_qargs)
     second_qarg = tuple(qarg[q] for q in second_qargs)
 
+    from qiskit.dagcircuit.dagnode import DAGOpNode
+
+    # If we have a DAGOpNode here we've received a StandardGate definition from
+    # rust and we can manually pull the matrix to use for the Operators
+    if isinstance(first_ops, DAGOpNode):
+        first_ops = first_ops.matrix
+    if isinstance(second_op, DAGOpNode):
+        second_op = second_op.matrix
+
     # try to generate an Operator out of op, if this succeeds we can determine commutativity, otherwise
     # return false
     try:
diff --git a/qiskit/circuit/library/data_preparation/__init__.py b/qiskit/circuit/library/data_preparation/__init__.py
index 38611c911fa4..192308a3a7f5 100644
--- a/qiskit/circuit/library/data_preparation/__init__.py
+++ b/qiskit/circuit/library/data_preparation/__init__.py
@@ -41,7 +41,14 @@
 from .pauli_feature_map import PauliFeatureMap
 from .z_feature_map import ZFeatureMap
 from .zz_feature_map import ZZFeatureMap
-from .state_preparation import StatePreparation
+from .state_preparation import StatePreparation, UniformSuperpositionGate
 from .initializer import Initialize
 
-__all__ = ["PauliFeatureMap", "ZFeatureMap", "ZZFeatureMap", "StatePreparation", "Initialize"]
+__all__ = [
+    "PauliFeatureMap",
+    "ZFeatureMap",
+    "ZZFeatureMap",
+    "StatePreparation",
+    "UniformSuperpositionGate",
+    "Initialize",
+]
diff --git a/qiskit/circuit/library/data_preparation/state_preparation.py b/qiskit/circuit/library/data_preparation/state_preparation.py
index 1d9ad7f7b082..26c37334cfe4 100644
--- a/qiskit/circuit/library/data_preparation/state_preparation.py
+++ b/qiskit/circuit/library/data_preparation/state_preparation.py
@@ -25,7 +25,9 @@
 from qiskit.circuit.library.standard_gates.s import SGate, SdgGate
 from qiskit.circuit.library.generalized_gates import Isometry
 from qiskit.circuit.exceptions import CircuitError
-from qiskit.quantum_info.states.statevector import Statevector  # pylint: disable=cyclic-import
+from qiskit.quantum_info.states.statevector import (
+    Statevector,
+)  # pylint: disable=cyclic-import
 
 _EPS = 1e-10  # global variable used to chop very small numbers to zero
 
@@ -240,3 +242,95 @@ def validate_parameter(self, parameter):
 
     def _return_repeat(self, exponent: float) -> "Gate":
         return Gate(name=f"{self.name}*{exponent}", num_qubits=self.num_qubits, params=[])
+
+
+class UniformSuperpositionGate(Gate):
+    r"""Implements a uniform superposition state.
+
+    This gate is used to create the uniform superposition state
+    :math:`\frac{1}{\sqrt{M}} \sum_{j=0}^{M-1}  |j\rangle` when it acts on an input
+    state :math:`|0...0\rangle`. Note, that `M` is not required to be
+    a power of 2, in which case the uniform superposition could be
+    prepared by a single layer of Hadamard gates.
+
+    .. note::
+
+        This class uses the Shukla-Vedula algorithm [1], which only needs
+        :math:`O(\log_2 (M))` qubits and :math:`O(\log_2 (M))` gates,
+        to prepare the superposition.
+
+    **References:**
+    [1]: A. Shukla and P. Vedula (2024), An efficient quantum algorithm for preparation
+    of uniform quantum superposition states, `Quantum Inf Process 23, 38
+    <https://link.springer.com/article/10.1007/s11128-024-04258-4>`_.
+    """
+
+    def __init__(
+        self,
+        num_superpos_states: int = 2,
+        num_qubits: Optional[int] = None,
+    ):
+        r"""
+        Args:
+            num_superpos_states (int):
+                A positive integer M = num_superpos_states (> 1) representing the number of computational
+                basis states with an amplitude of 1/sqrt(M) in the uniform superposition
+                state (:math:`\frac{1}{\sqrt{M}} \sum_{j=0}^{M-1}  |j\rangle`, where
+                :math:`1< M <= 2^n`). Note that the remaining (:math:`2^n - M`) computational basis
+                states have zero amplitudes. Here M need not be an integer power of 2.
+
+            num_qubits (int):
+                A positive integer representing the number of qubits used.  If num_qubits is None
+                or is not specified, then num_qubits is set to ceil(log2(num_superpos_states)).
+
+        Raises:
+            ValueError: num_qubits must be an integer greater than or equal to log2(num_superpos_states).
+
+        """
+        if num_superpos_states <= 1:
+            raise ValueError("num_superpos_states must be a positive integer greater than 1.")
+        if num_qubits is None:
+            num_qubits = int(math.ceil(math.log2(num_superpos_states)))
+        else:
+            if not (isinstance(num_qubits, int) and (num_qubits >= math.log2(num_superpos_states))):
+                raise ValueError(
+                    "num_qubits must be an integer greater than or equal to log2(num_superpos_states)."
+                )
+        super().__init__("USup", num_qubits, [num_superpos_states])
+
+    def _define(self):
+
+        qc = QuantumCircuit(self._num_qubits)
+
+        num_superpos_states = self.params[0]
+
+        if (
+            num_superpos_states & (num_superpos_states - 1)
+        ) == 0:  # if num_superpos_states is an integer power of 2
+            m = int(math.log2(num_superpos_states))
+            qc.h(range(m))
+            self.definition = qc
+            return
+
+        n_value = [int(x) for x in reversed(np.binary_repr(num_superpos_states))]
+        k = len(n_value)
+        l_value = [index for (index, item) in enumerate(n_value) if item == 1]  # Locations of '1's
+
+        qc.x(l_value[1:k])
+        m_current_value = 2 ** l_value[0]
+        theta = -2 * np.arccos(np.sqrt(m_current_value / num_superpos_states))
+
+        if l_value[0] > 0:  # if num_superpos_states is even
+            qc.h(range(l_value[0]))
+        qc.ry(theta, l_value[1])
+        qc.ch(l_value[1], range(l_value[0], l_value[1]), ctrl_state="0")
+
+        for m in range(1, len(l_value) - 1):
+            theta = -2 * np.arccos(
+                np.sqrt(2 ** l_value[m] / (num_superpos_states - m_current_value))
+            )
+            qc.cry(theta, l_value[m], l_value[m + 1], ctrl_state="0")
+            qc.ch(l_value[m + 1], range(l_value[m], l_value[m + 1]), ctrl_state="0")
+            m_current_value = m_current_value + 2 ** l_value[m]
+
+        self.definition = qc
diff --git a/qiskit/circuit/library/standard_gates/rzz.py b/qiskit/circuit/library/standard_gates/rzz.py
index 119dd370e20c..554ad4954a31 100644
--- a/qiskit/circuit/library/standard_gates/rzz.py
+++ b/qiskit/circuit/library/standard_gates/rzz.py
@@ -72,7 +72,7 @@ class RZZGate(Gate):
 
         .. math::
 
-            R_{ZZ}(\theta = \pi) = - Z \otimes Z
+            R_{ZZ}(\theta = \pi) = - i Z \otimes Z
 
         .. math::
 
diff --git a/qiskit/dagcircuit/dagcircuit.py b/qiskit/dagcircuit/dagcircuit.py
index a0d2c42be17e..b93a90e47f7b 100644
--- a/qiskit/dagcircuit/dagcircuit.py
+++ b/qiskit/dagcircuit/dagcircuit.py
@@ -54,6 +54,7 @@
 from qiskit.dagcircuit.dagnode import DAGNode, DAGOpNode, DAGInNode, DAGOutNode
 from qiskit.circuit.bit import Bit
 from qiskit.pulse import Schedule
+from qiskit._accelerate.euler_one_qubit_decomposer import collect_1q_runs_filter
 
 BitLocations = namedtuple("BitLocations", ("index", "registers"))
 # The allowable arguments to :meth:`DAGCircuit.copy_empty_like`'s ``vars_mode``.
@@ -642,17 +643,17 @@ def _check_wires(self, args: Iterable[Bit | expr.Var], amap: dict[Bit | expr.Var
             if wire not in amap:
                 raise DAGCircuitError(f"wire {wire} not found in {amap}")
 
-    def _increment_op(self, op):
-        if op.name in self._op_names:
-            self._op_names[op.name] += 1
+    def _increment_op(self, op_name):
+        if op_name in self._op_names:
+            self._op_names[op_name] += 1
         else:
-            self._op_names[op.name] = 1
+            self._op_names[op_name] = 1
 
-    def _decrement_op(self, op):
-        if self._op_names[op.name] == 1:
-            del self._op_names[op.name]
+    def _decrement_op(self, op_name):
+        if self._op_names[op_name] == 1:
+            del self._op_names[op_name]
         else:
-            self._op_names[op.name] -= 1
+            self._op_names[op_name] -= 1
 
     def copy_empty_like(self, *, vars_mode: _VarsMode = "alike"):
         """Return a copy of self with the same structure but empty.
@@ -724,7 +725,7 @@ def _apply_op_node_back(self, node: DAGOpNode):
             additional = set(_additional_wires(node)).difference(node.cargs)
 
         node._node_id = self._multi_graph.add_node(node)
-        self._increment_op(node)
+        self._increment_op(node.name)
 
         # Add new in-edges from predecessors of the output nodes to the
         # operation node while deleting the old in-edges of the output nodes
@@ -780,7 +781,7 @@ def apply_operation_back(
 
         node = DAGOpNode(op=op, qargs=qargs, cargs=cargs, dag=self)
         node._node_id = self._multi_graph.add_node(node)
-        self._increment_op(op)
+        self._increment_op(op.name)
 
         # Add new in-edges from predecessors of the output nodes to the
         # operation node while deleting the old in-edges of the output nodes
@@ -832,7 +833,7 @@ def apply_operation_front(
 
         node = DAGOpNode(op=op, qargs=qargs, cargs=cargs, dag=self)
         node._node_id = self._multi_graph.add_node(node)
-        self._increment_op(op)
+        self._increment_op(op.name)
 
         # Add new out-edges to successors of the input nodes from the
         # operation node while deleting the old out-edges of the input nodes
@@ -1379,10 +1380,10 @@ def replace_block_with_op(
                 "Replacing the specified node block would introduce a cycle"
             ) from ex
 
-        self._increment_op(op)
+        self._increment_op(op.name)
 
         for nd in node_block:
-            self._decrement_op(nd.op)
+            self._decrement_op(nd.name)
 
         return new_node
 
@@ -1593,7 +1594,7 @@ def edge_weight_map(wire):
         node_map = self._multi_graph.substitute_node_with_subgraph(
             node._node_id, in_dag._multi_graph, edge_map_fn, filter_fn, edge_weight_map
         )
-        self._decrement_op(node.op)
+        self._decrement_op(node.name)
 
         variable_mapper = _classical_resource_map.VariableMapper(
             self.cregs.values(), wire_map, add_register=self.add_creg
@@ -1624,7 +1625,7 @@ def edge_weight_map(wire):
             new_node = DAGOpNode(m_op, qargs=m_qargs, cargs=m_cargs, dag=self)
             new_node._node_id = new_node_index
             self._multi_graph[new_node_index] = new_node
-            self._increment_op(new_node.op)
+            self._increment_op(new_node.name)
 
         return {k: self._multi_graph[v] for k, v in node_map.items()}
 
@@ -1696,17 +1697,17 @@ def substitute_node(self, node: DAGOpNode, op, inplace: bool = False, propagate_
 
         if inplace:
             if op.name != node.op.name:
-                self._increment_op(op)
-                self._decrement_op(node.op)
+                self._increment_op(op.name)
+                self._decrement_op(node.name)
             node.op = op
             return node
 
         new_node = copy.copy(node)
         new_node.op = op
         self._multi_graph[node._node_id] = new_node
-        if op.name != node.op.name:
-            self._increment_op(op)
-            self._decrement_op(node.op)
+        if op.name != node.name:
+            self._increment_op(op.name)
+            self._decrement_op(node.name)
         return new_node
 
     def separable_circuits(
@@ -1987,7 +1988,7 @@ def remove_op_node(self, node):
         self._multi_graph.remove_node_retain_edges(
             node._node_id, use_outgoing=False, condition=lambda edge1, edge2: edge1 == edge2
         )
-        self._decrement_op(node.op)
+        self._decrement_op(node.name)
 
     def remove_ancestors_of(self, node):
         """Remove all of the ancestor operation nodes of node."""
@@ -2152,19 +2153,7 @@ def filter_fn(node):
 
     def collect_1q_runs(self) -> list[list[DAGOpNode]]:
         """Return a set of non-conditional runs of 1q "op" nodes."""
-
-        def filter_fn(node):
-            return (
-                isinstance(node, DAGOpNode)
-                and len(node.qargs) == 1
-                and len(node.cargs) == 0
-                and isinstance(node.op, Gate)
-                and hasattr(node.op, "__array__")
-                and getattr(node.op, "condition", None) is None
-                and not node.op.is_parameterized()
-            )
-
-        return rx.collect_runs(self._multi_graph, filter_fn)
+        return rx.collect_runs(self._multi_graph, collect_1q_runs_filter)
 
     def collect_2q_runs(self):
         """Return a set of non-conditional runs of 2q "op" nodes."""
diff --git a/qiskit/synthesis/clifford/clifford_decompose_greedy.py b/qiskit/synthesis/clifford/clifford_decompose_greedy.py
index 784e6706d62e..0a679a8a7a6f 100644
--- a/qiskit/synthesis/clifford/clifford_decompose_greedy.py
+++ b/qiskit/synthesis/clifford/clifford_decompose_greedy.py
@@ -12,22 +12,16 @@
 """
 Circuit synthesis for the Clifford class.
 """
-# pylint: disable=invalid-name
 
 # ---------------------------------------------------------------------
 # Synthesis based on Bravyi et. al. greedy clifford compiler
 # ---------------------------------------------------------------------
 
-
-import numpy as np
 from qiskit.circuit import QuantumCircuit
-from qiskit.exceptions import QiskitError
-from qiskit.quantum_info import Clifford, Pauli
-from qiskit.quantum_info.operators.symplectic.clifford_circuits import (
-    _append_cx,
-    _append_h,
-    _append_s,
-    _append_swap,
+from qiskit.quantum_info import Clifford
+
+from qiskit._accelerate.synthesis.clifford import (
+    synth_clifford_greedy as synth_clifford_greedy_inner,
 )
 
 
@@ -56,296 +50,8 @@ def synth_clifford_greedy(clifford: Clifford) -> QuantumCircuit:
            *Clifford Circuit Optimization with Templates and Symbolic Pauli Gates*,
            `arXiv:2105.02291 [quant-ph] <https://arxiv.org/abs/2105.02291>`_
     """
-
-    num_qubits = clifford.num_qubits
-    circ = QuantumCircuit(num_qubits, name=str(clifford))
-    qubit_list = list(range(num_qubits))
-    clifford_cpy = clifford.copy()
-
-    # Reducing the original Clifford to identity
-    # via symplectic Gaussian elimination
-    while len(qubit_list) > 0:
-        # Calculate the adjoint of clifford_cpy without the phase
-        clifford_adj = clifford_cpy.copy()
-        tmp = clifford_adj.destab_x.copy()
-        clifford_adj.destab_x = clifford_adj.stab_z.T
-        clifford_adj.destab_z = clifford_adj.destab_z.T
-        clifford_adj.stab_x = clifford_adj.stab_x.T
-        clifford_adj.stab_z = tmp.T
-
-        list_greedy_cost = []
-        for qubit in qubit_list:
-            pauli_x = Pauli("I" * (num_qubits - qubit - 1) + "X" + "I" * qubit)
-            pauli_x = pauli_x.evolve(clifford_adj, frame="s")
-
-            pauli_z = Pauli("I" * (num_qubits - qubit - 1) + "Z" + "I" * qubit)
-            pauli_z = pauli_z.evolve(clifford_adj, frame="s")
-            list_pairs = []
-            pauli_count = 0
-
-            # Compute the CNOT cost in order to find the qubit with the minimal cost
-            for i in qubit_list:
-                typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, i)
-                list_pairs.append(typeq)
-                pauli_count += 1
-            cost = _compute_greedy_cost(list_pairs)
-            list_greedy_cost.append([cost, qubit])
-
-        _, min_qubit = (sorted(list_greedy_cost))[0]
-
-        # Gaussian elimination step for the qubit with minimal CNOT cost
-        pauli_x = Pauli("I" * (num_qubits - min_qubit - 1) + "X" + "I" * min_qubit)
-        pauli_x = pauli_x.evolve(clifford_adj, frame="s")
-
-        pauli_z = Pauli("I" * (num_qubits - min_qubit - 1) + "Z" + "I" * min_qubit)
-        pauli_z = pauli_z.evolve(clifford_adj, frame="s")
-
-        # Compute the decoupling operator of cliff_ox and cliff_oz
-        decouple_circ, decouple_cliff = _calc_decoupling(
-            pauli_x, pauli_z, qubit_list, min_qubit, num_qubits, clifford_cpy
-        )
-        circ = circ.compose(decouple_circ)
-
-        # Now the clifford acts trivially on min_qubit
-        clifford_cpy = decouple_cliff.adjoint().compose(clifford_cpy)
-        qubit_list.remove(min_qubit)
-
-    # Add the phases (Pauli gates) to the Clifford circuit
-    for qubit in range(num_qubits):
-        stab = clifford_cpy.stab_phase[qubit]
-        destab = clifford_cpy.destab_phase[qubit]
-        if destab and stab:
-            circ.y(qubit)
-        elif not destab and stab:
-            circ.x(qubit)
-        elif destab and not stab:
-            circ.z(qubit)
-
-    return circ
-
-
-# ---------------------------------------------------------------------
-# Helper functions for Bravyi et. al. greedy clifford compiler
-# ---------------------------------------------------------------------
-
-# Global arrays of the 16 pairs of Pauli operators
-# divided into 5 equivalence classes under the action of single-qubit Cliffords
-
-# Class A - canonical representative is 'XZ'
-A_class = [
-    [[False, True], [True, True]],  # 'XY'
-    [[False, True], [True, False]],  # 'XZ'
-    [[True, True], [False, True]],  # 'YX'
-    [[True, True], [True, False]],  # 'YZ'
-    [[True, False], [False, True]],  # 'ZX'
-    [[True, False], [True, True]],
-]  # 'ZY'
-
-# Class B - canonical representative is 'XX'
-B_class = [
-    [[True, False], [True, False]],  # 'ZZ'
-    [[False, True], [False, True]],  # 'XX'
-    [[True, True], [True, True]],
-]  # 'YY'
-
-# Class C - canonical representative is 'XI'
-C_class = [
-    [[True, False], [False, False]],  # 'ZI'
-    [[False, True], [False, False]],  # 'XI'
-    [[True, True], [False, False]],
-]  # 'YI'
-
-# Class D - canonical representative is 'IZ'
-D_class = [
-    [[False, False], [False, True]],  # 'IX'
-    [[False, False], [True, False]],  # 'IZ'
-    [[False, False], [True, True]],
-]  # 'IY'
-
-# Class E - only 'II'
-E_class = [[[False, False], [False, False]]]  # 'II'
-
-
-def _from_pair_paulis_to_type(pauli_x, pauli_z, qubit):
-    """Converts a pair of Paulis pauli_x and pauli_z into a type"""
-
-    type_x = [pauli_x.z[qubit], pauli_x.x[qubit]]
-    type_z = [pauli_z.z[qubit], pauli_z.x[qubit]]
-    return [type_x, type_z]
-
-
-def _compute_greedy_cost(list_pairs):
-    """Compute the CNOT cost of one step of the algorithm"""
-
-    A_num = 0
-    B_num = 0
-    C_num = 0
-    D_num = 0
-
-    for pair in list_pairs:
-        if pair in A_class:
-            A_num += 1
-        elif pair in B_class:
-            B_num += 1
-        elif pair in C_class:
-            C_num += 1
-        elif pair in D_class:
-            D_num += 1
-
-    if (A_num % 2) == 0:
-        raise QiskitError("Symplectic Gaussian elimination fails.")
-
-    # Calculate the CNOT cost
-    cost = 3 * (A_num - 1) / 2 + (B_num + 1) * (B_num > 0) + C_num + D_num
-    if list_pairs[0] not in A_class:  # additional SWAP
-        cost += 3
-
-    return cost
-
-
-def _calc_decoupling(pauli_x, pauli_z, qubit_list, min_qubit, num_qubits, cliff):
-    """Calculate a decoupling operator D:
-    D^{-1} * Ox * D = x1
-    D^{-1} * Oz * D = z1
-    and reduce the clifford such that it will act trivially on min_qubit
-    """
-
-    circ = QuantumCircuit(num_qubits)
-
-    # decouple_cliff is initialized to an identity clifford
-    decouple_cliff = cliff.copy()
-    num_qubits = decouple_cliff.num_qubits
-    decouple_cliff.phase = np.zeros(2 * num_qubits)
-    decouple_cliff.symplectic_matrix = np.eye(2 * num_qubits)
-
-    qubit0 = min_qubit  # The qubit for the symplectic Gaussian elimination
-
-    # Reduce the pair of Paulis to a representative in the equivalence class
-    # ['XZ', 'XX', 'XI', 'IZ', 'II'] by adding single-qubit gates
-    for qubit in qubit_list:
-
-        typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, qubit)
-
-        if typeq in [
-            [[True, True], [False, False]],  # 'YI'
-            [[True, True], [True, True]],  # 'YY'
-            [[True, True], [True, False]],
-        ]:  # 'YZ':
-            circ.s(qubit)
-            _append_s(decouple_cliff, qubit)
-
-        elif typeq in [
-            [[True, False], [False, False]],  # 'ZI'
-            [[True, False], [True, False]],  # 'ZZ'
-            [[True, False], [False, True]],  # 'ZX'
-            [[False, False], [False, True]],
-        ]:  # 'IX'
-            circ.h(qubit)
-            _append_h(decouple_cliff, qubit)
-
-        elif typeq in [
-            [[False, False], [True, True]],  # 'IY'
-            [[True, False], [True, True]],
-        ]:  # 'ZY'
-            circ.s(qubit)
-            circ.h(qubit)
-            _append_s(decouple_cliff, qubit)
-            _append_h(decouple_cliff, qubit)
-
-        elif typeq == [[True, True], [False, True]]:  # 'YX'
-            circ.h(qubit)
-            circ.s(qubit)
-            _append_h(decouple_cliff, qubit)
-            _append_s(decouple_cliff, qubit)
-
-        elif typeq == [[False, True], [True, True]]:  # 'XY'
-            circ.s(qubit)
-            circ.h(qubit)
-            circ.s(qubit)
-            _append_s(decouple_cliff, qubit)
-            _append_h(decouple_cliff, qubit)
-            _append_s(decouple_cliff, qubit)
-
-    # Reducing each pair of Paulis (except of qubit0) to 'II'
-    # by adding two-qubit gates and single-qubit gates
-    A_qubits = []
-    B_qubits = []
-    C_qubits = []
-    D_qubits = []
-
-    for qubit in qubit_list:
-        typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, qubit)
-        if typeq in A_class:
-            A_qubits.append(qubit)
-        elif typeq in B_class:
-            B_qubits.append(qubit)
-        elif typeq in C_class:
-            C_qubits.append(qubit)
-        elif typeq in D_class:
-            D_qubits.append(qubit)
-
-    if len(A_qubits) % 2 != 1:
-        raise QiskitError("Symplectic Gaussian elimination fails.")
-
-    if qubit0 not in A_qubits:  # SWAP qubit0 and qubitA
-        qubitA = A_qubits[0]
-        circ.swap(qubit0, qubitA)
-        _append_swap(decouple_cliff, qubit0, qubitA)
-        if qubit0 in B_qubits:
-            B_qubits.remove(qubit0)
-            B_qubits.append(qubitA)
-            A_qubits.remove(qubitA)
-            A_qubits.append(qubit0)
-        elif qubit0 in C_qubits:
-            C_qubits.remove(qubit0)
-            C_qubits.append(qubitA)
-            A_qubits.remove(qubitA)
-            A_qubits.append(qubit0)
-        elif qubit0 in D_qubits:
-            D_qubits.remove(qubit0)
-            D_qubits.append(qubitA)
-            A_qubits.remove(qubitA)
-            A_qubits.append(qubit0)
-        else:
-            A_qubits.remove(qubitA)
-            A_qubits.append(qubit0)
-
-    # Reduce pairs in Class C to 'II'
-    for qubit in C_qubits:
-        circ.cx(qubit0, qubit)
-        _append_cx(decouple_cliff, qubit0, qubit)
-
-    # Reduce pairs in Class D to 'II'
-    for qubit in D_qubits:
-        circ.cx(qubit, qubit0)
-        _append_cx(decouple_cliff, qubit, qubit0)
-
-    # Reduce pairs in Class B to 'II'
-    if len(B_qubits) > 1:
-        for qubit in B_qubits[1:]:
-            qubitB = B_qubits[0]
-            circ.cx(qubitB, qubit)
-            _append_cx(decouple_cliff, qubitB, qubit)
-
-    if len(B_qubits) > 0:
-        qubitB = B_qubits[0]
-        circ.cx(qubit0, qubitB)
-        circ.h(qubitB)
-        circ.cx(qubitB, qubit0)
-        _append_cx(decouple_cliff, qubit0, qubitB)
-        _append_h(decouple_cliff, qubitB)
-        _append_cx(decouple_cliff, qubitB, qubit0)
-
-    # Reduce pairs in Class A (except of qubit0) to 'II'
-    Alen = int((len(A_qubits) - 1) / 2)
-    if Alen > 0:
-        A_qubits.remove(qubit0)
-    for qubit in range(Alen):
-        circ.cx(A_qubits[2 * qubit + 1], A_qubits[2 * qubit])
-        circ.cx(A_qubits[2 * qubit], qubit0)
-        circ.cx(qubit0, A_qubits[2 * qubit + 1])
-        _append_cx(decouple_cliff, A_qubits[2 * qubit + 1], A_qubits[2 * qubit])
-        _append_cx(decouple_cliff, A_qubits[2 * qubit], qubit0)
-        _append_cx(decouple_cliff, qubit0, A_qubits[2 * qubit + 1])
-
-    return circ, decouple_cliff
+    circuit = QuantumCircuit._from_circuit_data(
+        synth_clifford_greedy_inner(clifford.tableau.astype(bool))
+    )
+    circuit.name = str(clifford)
+    return circuit
diff --git a/qiskit/synthesis/linear/cnot_synth.py b/qiskit/synthesis/linear/cnot_synth.py
index 699523a7e75d..f900ff86d17e 100644
--- a/qiskit/synthesis/linear/cnot_synth.py
+++ b/qiskit/synthesis/linear/cnot_synth.py
@@ -17,33 +17,37 @@
 """
 
 from __future__ import annotations
-import copy
+
 import numpy as np
 from qiskit.circuit import QuantumCircuit
-from qiskit.exceptions import QiskitError
+
+from qiskit._accelerate.synthesis.linear import synth_cnot_count_full_pmh as fast_pmh
 
 
 def synth_cnot_count_full_pmh(
-    state: list[list[bool]] | np.ndarray[bool], section_size: int = 2
+    state: list[list[bool]] | np.ndarray[bool], section_size: int | None = None
 ) -> QuantumCircuit:
-    """
+    r"""
     Synthesize linear reversible circuits for all-to-all architecture
     using Patel, Markov and Hayes method.
 
     This function is an implementation of the Patel, Markov and Hayes algorithm from [1]
     for optimal synthesis of linear reversible circuits for all-to-all architecture,
-    as specified by an :math:`n \\times n` matrix.
+    as specified by an :math:`n \times n` matrix.
 
     Args:
-        state: :math:`n \\times n` boolean invertible matrix, describing
-            the state of the input circuit
+        state: :math:`n \times n` boolean invertible matrix, describing
+            the state of the input circuit.
         section_size: The size of each section in the Patel–Markov–Hayes algorithm [1].
+            If ``None`` it is chosen to be :math:`\max(2, \alpha\log_2(n))` with
+            :math:`\alpha = 0.56`, which approximately minimizes the upper bound on the number
+            of row operations given in [1] Eq. (3).
 
     Returns:
-        QuantumCircuit: a CX-only circuit implementing the linear transformation.
+        A CX-only circuit implementing the linear transformation.
 
     Raises:
-        QiskitError: when variable ``state`` isn't of type ``numpy.ndarray``
+        ValueError: When ``section_size`` is larger than the number of columns.
 
     References:
         1. Patel, Ketan N., Igor L. Markov, and John P. Hayes,
@@ -51,92 +55,15 @@ def synth_cnot_count_full_pmh(
            Quantum Information & Computation 8.3 (2008): 282-294.
            `arXiv:quant-ph/0302002 [quant-ph] <https://arxiv.org/abs/quant-ph/0302002>`_
     """
-    if not isinstance(state, (list, np.ndarray)):
-        raise QiskitError(
-            f"state should be of type list or numpy.ndarray, but was of the type {type(state)}"
+    normalized = np.asarray(state).astype(bool)
+    if section_size is not None and normalized.shape[1] < section_size:
+        raise ValueError(
+            f"The section_size ({section_size}) cannot be larger than the number of columns "
+            f"({normalized.shape[1]})."
         )
-    state = np.array(state)
-    # Synthesize lower triangular part
-    [state, circuit_l] = _lwr_cnot_synth(state, section_size)
-    state = np.transpose(state)
-    # Synthesize upper triangular part
-    [state, circuit_u] = _lwr_cnot_synth(state, section_size)
-    circuit_l.reverse()
-    for i in circuit_u:
-        i.reverse()
-    # Convert the list into a circuit of C-NOT gates
-    circ = QuantumCircuit(state.shape[0])
-    for i in circuit_u + circuit_l:
-        circ.cx(i[0], i[1])
-    return circ
-
-
-def _lwr_cnot_synth(state, section_size):
-    """
-    This function is a helper function of the algorithm for optimal synthesis
-    of linear reversible circuits (the Patel–Markov–Hayes algorithm). It works
-    like gaussian elimination, except that it works a lot faster, and requires
-    fewer steps (and therefore fewer CNOTs). It takes the matrix "state" and
-    splits it into sections of size section_size. Then it eliminates all non-zero
-    sub-rows within each section, which are the same as a non-zero sub-row
-    above. Once this has been done, it continues with normal gaussian elimination.
-    The benefit is that with small section sizes (m), most of the sub-rows will
-    be cleared in the first step, resulting in a factor m fewer row row operations
-    during Gaussian elimination.
 
-    The algorithm is described in detail in the following paper
-    "Optimal synthesis of linear reversible circuits."
-    Patel, Ketan N., Igor L. Markov, and John P. Hayes.
-    Quantum Information & Computation 8.3 (2008): 282-294.
-
-    Note:
-    This implementation tweaks the Patel, Markov, and Hayes algorithm by adding
-    a "back reduce" which adds rows below the pivot row with a high degree of
-    overlap back to it. The intuition is to avoid a high-weight pivot row
-    increasing the weight of lower rows.
-
-    Args:
-        state (ndarray): n x n matrix, describing a linear quantum circuit
-        section_size (int): the section size the matrix columns are divided into
-
-    Returns:
-        numpy.matrix: n by n matrix, describing the state of the output circuit
-        list: a k by 2 list of C-NOT operations that need to be applied
-    """
-    circuit = []
-    num_qubits = state.shape[0]
-    cutoff = 1
+    # call Rust implementation with normalized input
+    circuit_data = fast_pmh(normalized, section_size)
 
-    # Iterate over column sections
-    for sec in range(1, int(np.floor(num_qubits / section_size) + 1)):
-        # Remove duplicate sub-rows in section sec
-        patt = {}
-        for row in range((sec - 1) * section_size, num_qubits):
-            sub_row_patt = copy.deepcopy(state[row, (sec - 1) * section_size : sec * section_size])
-            if np.sum(sub_row_patt) == 0:
-                continue
-            if str(sub_row_patt) not in patt:
-                patt[str(sub_row_patt)] = row
-            else:
-                state[row, :] ^= state[patt[str(sub_row_patt)], :]
-                circuit.append([patt[str(sub_row_patt)], row])
-        # Use gaussian elimination for remaining entries in column section
-        for col in range((sec - 1) * section_size, sec * section_size):
-            # Check if 1 on diagonal
-            diag_one = 1
-            if state[col, col] == 0:
-                diag_one = 0
-            # Remove ones in rows below column col
-            for row in range(col + 1, num_qubits):
-                if state[row, col] == 1:
-                    if diag_one == 0:
-                        state[col, :] ^= state[row, :]
-                        circuit.append([row, col])
-                        diag_one = 1
-                    state[row, :] ^= state[col, :]
-                    circuit.append([col, row])
-                # Back reduce the pivot row using the current row
-                if sum(state[col, :] & state[row, :]) > cutoff:
-                    state[col, :] ^= state[row, :]
-                    circuit.append([row, col])
-    return [state, circuit]
+    # construct circuit from the data
+    return QuantumCircuit._from_circuit_data(circuit_data)
diff --git a/qiskit/synthesis/one_qubit/one_qubit_decompose.py b/qiskit/synthesis/one_qubit/one_qubit_decompose.py
index c84db761b7f0..f60f20f9524e 100644
--- a/qiskit/synthesis/one_qubit/one_qubit_decompose.py
+++ b/qiskit/synthesis/one_qubit/one_qubit_decompose.py
@@ -161,29 +161,19 @@ def build_circuit(self, gates, global_phase) -> QuantumCircuit | DAGCircuit:
         if len(gates) > 0 and isinstance(gates[0], tuple):
             lookup_gate = True
 
-        if self.use_dag:
-            from qiskit.dagcircuit import dagcircuit
-
-            dag = dagcircuit.DAGCircuit()
-            dag.global_phase = global_phase
-            dag.add_qubits(qr)
-            for gate_entry in gates:
-                if lookup_gate:
-                    gate = NAME_MAP[gate_entry[0]](*gate_entry[1])
-                else:
-                    gate = gate_entry
-
-                dag.apply_operation_back(gate, (qr[0],), check=False)
-            return dag
-        else:
-            circuit = QuantumCircuit(qr, global_phase=global_phase)
-            for gate_entry in gates:
-                if lookup_gate:
-                    gate = NAME_MAP[gate_entry[0]](*gate_entry[1])
-                else:
-                    gate = gate_entry
-                circuit._append(gate, [qr[0]], [])
-            return circuit
+        from qiskit.dagcircuit import dagcircuit
+
+        dag = dagcircuit.DAGCircuit()
+        dag.global_phase = global_phase
+        dag.add_qubits(qr)
+        for gate_entry in gates:
+            if lookup_gate:
+                gate = NAME_MAP[gate_entry[0].name](*gate_entry[1])
+            else:
+                gate = gate_entry.name
+
+            dag.apply_operation_back(gate, (qr[0],), check=False)
+        return dag
 
     def __call__(
         self,
@@ -225,11 +215,17 @@ def __call__(
         return self._decompose(unitary, simplify=simplify, atol=atol)
 
     def _decompose(self, unitary, simplify=True, atol=DEFAULT_ATOL):
-        circuit_sequence = euler_one_qubit_decomposer.unitary_to_gate_sequence(
-            unitary, [self.basis], 0, None, simplify, atol
+        if self.use_dag:
+            circuit_sequence = euler_one_qubit_decomposer.unitary_to_gate_sequence(
+                unitary, [self.basis], 0, None, simplify, atol
+            )
+            circuit = self.build_circuit(circuit_sequence, circuit_sequence.global_phase)
+            return circuit
+        return QuantumCircuit._from_circuit_data(
+            euler_one_qubit_decomposer.unitary_to_circuit(
+                unitary, [self.basis], 0, None, simplify, atol
+            )
         )
-        circuit = self.build_circuit(circuit_sequence, circuit_sequence.global_phase)
-        return circuit
 
     @property
     def basis(self):
diff --git a/qiskit/transpiler/layout.py b/qiskit/transpiler/layout.py
index 4117e2987bb6..bece19671794 100644
--- a/qiskit/transpiler/layout.py
+++ b/qiskit/transpiler/layout.py
@@ -454,7 +454,7 @@ class TranspileLayout:
     qubits in the circuit as it fits the circuit to the target backend. For example,
     let the input circuit be:
 
-    .. plot:
+    .. plot::
        :include-source:
 
        from qiskit.circuit import QuantumCircuit, QuantumRegister
@@ -469,7 +469,7 @@ class TranspileLayout:
 
     Suppose that during the layout stage the transpiler reorders the qubits to be:
 
-    .. plot:
+    .. plot::
        :include-source:
 
        from qiskit import QuantumCircuit
@@ -497,7 +497,7 @@ class TranspileLayout:
     the transpiler needs to insert swap gates, and the output circuit
     becomes:
 
-    .. plot:
+    .. plot::
        :include-source:
 
        from qiskit import QuantumCircuit
diff --git a/qiskit/transpiler/passes/optimization/commutation_analysis.py b/qiskit/transpiler/passes/optimization/commutation_analysis.py
index eddb659f0a25..61c77de552b9 100644
--- a/qiskit/transpiler/passes/optimization/commutation_analysis.py
+++ b/qiskit/transpiler/passes/optimization/commutation_analysis.py
@@ -72,14 +72,7 @@ def run(self, dag):
                         does_commute = (
                             isinstance(current_gate, DAGOpNode)
                             and isinstance(prev_gate, DAGOpNode)
-                            and self.comm_checker.commute(
-                                current_gate.op,
-                                current_gate.qargs,
-                                current_gate.cargs,
-                                prev_gate.op,
-                                prev_gate.qargs,
-                                prev_gate.cargs,
-                            )
+                            and self.comm_checker.commute_nodes(current_gate, prev_gate)
                         )
                         if not does_commute:
                             break
diff --git a/qiskit/transpiler/passes/optimization/commutative_cancellation.py b/qiskit/transpiler/passes/optimization/commutative_cancellation.py
index 4c6c487a0ea3..5c0b7317aabd 100644
--- a/qiskit/transpiler/passes/optimization/commutative_cancellation.py
+++ b/qiskit/transpiler/passes/optimization/commutative_cancellation.py
@@ -24,7 +24,7 @@
 from qiskit.circuit.library.standard_gates.rx import RXGate
 from qiskit.circuit.library.standard_gates.p import PhaseGate
 from qiskit.circuit.library.standard_gates.rz import RZGate
-from qiskit.circuit import ControlFlowOp
+from qiskit.circuit.controlflow import CONTROL_FLOW_OP_NAMES
 
 
 _CUTOFF_PRECISION = 1e-5
@@ -138,14 +138,14 @@ def run(self, dag):
                 total_phase = 0.0
                 for current_node in run:
                     if (
-                        getattr(current_node.op, "condition", None) is not None
+                        current_node.condition is not None
                         or len(current_node.qargs) != 1
                         or current_node.qargs[0] != run_qarg
                     ):
                         raise RuntimeError("internal error")
 
                     if current_node.name in ["p", "u1", "rz", "rx"]:
-                        current_angle = float(current_node.op.params[0])
+                        current_angle = float(current_node.params[0])
                     elif current_node.name in ["z", "x"]:
                         current_angle = np.pi
                     elif current_node.name == "t":
@@ -159,8 +159,8 @@ def run(self, dag):
 
                     # Compose gates
                     total_angle = current_angle + total_angle
-                    if current_node.op.definition:
-                        total_phase += current_node.op.definition.global_phase
+                    if current_node.definition:
+                        total_phase += current_node.definition.global_phase
 
                 # Replace the data of the first node in the run
                 if cancel_set_key[0] == "z_rotation":
@@ -200,7 +200,9 @@ def _handle_control_flow_ops(self, dag):
         """
 
         pass_manager = PassManager([CommutationAnalysis(), self])
-        for node in dag.op_nodes(ControlFlowOp):
+        for node in dag.op_nodes():
+            if node.name not in CONTROL_FLOW_OP_NAMES:
+                continue
             mapped_blocks = []
             for block in node.op.blocks:
                 new_circ = pass_manager.run(block)
diff --git a/qiskit/transpiler/passes/optimization/optimize_1q_decomposition.py b/qiskit/transpiler/passes/optimization/optimize_1q_decomposition.py
index 3f8d07839c0f..04d95312aa6d 100644
--- a/qiskit/transpiler/passes/optimization/optimize_1q_decomposition.py
+++ b/qiskit/transpiler/passes/optimization/optimize_1q_decomposition.py
@@ -33,6 +33,7 @@
     XGate,
 )
 from qiskit.circuit import Qubit
+from qiskit.circuit.quantumcircuitdata import CircuitInstruction
 from qiskit.dagcircuit.dagcircuit import DAGCircuit
 from qiskit.dagcircuit.dagnode import DAGOpNode
 
@@ -110,16 +111,7 @@ def _build_error_map(self):
         else:
             return None
 
-    def _resynthesize_run(self, matrix, qubit=None):
-        """
-        Re-synthesizes one 2x2 `matrix`, typically extracted via `dag.collect_1q_runs`.
-
-        Returns the newly synthesized circuit in the indicated basis, or None
-        if no synthesis routine applied.
-
-        When multiple synthesis options are available, it prefers the one with the lowest
-        error when the circuit is applied to `qubit`.
-        """
+    def _get_decomposer(self, qubit=None):
         # include path for when target exists but target.num_qubits is None (BasicSimulator)
         if self._target is not None and self._target.num_qubits is not None:
             if qubit is not None:
@@ -133,6 +125,19 @@ def _resynthesize_run(self, matrix, qubit=None):
                 decomposers = _possible_decomposers(available_1q_basis)
         else:
             decomposers = self._global_decomposers
+        return decomposers
+
+    def _resynthesize_run(self, matrix, qubit=None):
+        """
+        Re-synthesizes one 2x2 `matrix`, typically extracted via `dag.collect_1q_runs`.
+
+        Returns the newly synthesized circuit in the indicated basis, or None
+        if no synthesis routine applied.
+
+        When multiple synthesis options are available, it prefers the one with the lowest
+        error when the circuit is applied to `qubit`.
+        """
+        decomposers = self._get_decomposer(qubit)
 
         best_synth_circuit = euler_one_qubit_decomposer.unitary_to_gate_sequence(
             matrix,
@@ -149,10 +154,13 @@ def _gate_sequence_to_dag(self, best_synth_circuit):
         out_dag.global_phase = best_synth_circuit.global_phase
 
         for gate_name, angles in best_synth_circuit:
-            out_dag.apply_operation_back(NAME_MAP[gate_name](*angles), qubits, check=False)
+            op = CircuitInstruction(gate_name, qubits=qubits, params=angles)
+            out_dag.apply_operation_back(op.operation, qubits, check=False)
         return out_dag
 
-    def _substitution_checks(self, dag, old_run, new_circ, basis, qubit):
+    def _substitution_checks(
+        self, dag, old_run, new_circ, basis, qubit, old_error=None, new_error=None
+    ):
         """
         Returns `True` when it is recommended to replace `old_run` with `new_circ` over `basis`.
         """
@@ -176,11 +184,14 @@ def _substitution_checks(self, dag, old_run, new_circ, basis, qubit):
         # if we're outside of the basis set, we're obligated to logically decompose.
         # if we're outside of the set of gates for which we have physical definitions,
         #    then we _try_ to decompose, using the results if we see improvement.
-        new_error = 0.0
-        old_error = 0.0
         if not uncalibrated_and_not_basis_p:
-            new_error = self._error(new_circ, qubit)
-            old_error = self._error(old_run, qubit)
+            if new_error is None:
+                new_error = self._error(new_circ, qubit)
+            if old_error is None:
+                old_error = self._error(old_run, qubit)
+        else:
+            new_error = 0.0
+            old_error = 0.0
 
         return (
             uncalibrated_and_not_basis_p
@@ -198,32 +209,47 @@ def run(self, dag):
         Returns:
             DAGCircuit: the optimized DAG.
         """
-        runs = dag.collect_1q_runs()
-        for run in runs:
+        runs = []
+        qubits = []
+        bases = []
+        for run in dag.collect_1q_runs():
             qubit = dag.find_bit(run[0].qargs[0]).index
-            operator = run[0].op.to_matrix()
-            for node in run[1:]:
-                operator = node.op.to_matrix().dot(operator)
-            best_circuit_sequence = self._resynthesize_run(operator, qubit)
-
+            runs.append(run)
+            qubits.append(qubit)
+            bases.append(self._get_decomposer(qubit))
+        best_sequences = euler_one_qubit_decomposer.optimize_1q_gates_decomposition(
+            runs, qubits, bases, simplify=True, error_map=self.error_map
+        )
+        for index, best_circuit_sequence in enumerate(best_sequences):
+            run = runs[index]
+            qubit = qubits[index]
             if self._target is None:
                 basis = self._basis_gates
             else:
                 basis = self._target.operation_names_for_qargs((qubit,))
-
-            if best_circuit_sequence is not None and self._substitution_checks(
-                dag, run, best_circuit_sequence, basis, qubit
-            ):
-                for gate_name, angles in best_circuit_sequence:
-                    op = NAME_MAP[gate_name](*angles)
-                    node = DAGOpNode(NAME_MAP[gate_name](*angles), run[0].qargs, dag=dag)
-                    node._node_id = dag._multi_graph.add_node(node)
-                    dag._increment_op(op)
-                    dag._multi_graph.insert_node_on_in_edges(node._node_id, run[0]._node_id)
-                dag.global_phase += best_circuit_sequence.global_phase
-                # Delete the other nodes in the run
-                for current_node in run:
-                    dag.remove_op_node(current_node)
+            if best_circuit_sequence is not None:
+                (old_error, new_error, best_circuit_sequence) = best_circuit_sequence
+                if self._substitution_checks(
+                    dag,
+                    run,
+                    best_circuit_sequence,
+                    basis,
+                    qubit,
+                    old_error=old_error,
+                    new_error=new_error,
+                ):
+                    first_node_id = run[0]._node_id
+                    qubit = run[0].qargs
+                    for gate, angles in best_circuit_sequence:
+                        op = CircuitInstruction(gate, qubits=qubit, params=angles)
+                        node = DAGOpNode.from_instruction(op, dag=dag)
+                        node._node_id = dag._multi_graph.add_node(node)
+                        dag._increment_op(gate.name)
+                        dag._multi_graph.insert_node_on_in_edges(node._node_id, first_node_id)
+                    dag.global_phase += best_circuit_sequence.global_phase
+                    # Delete the other nodes in the run
+                    for current_node in run:
+                        dag.remove_op_node(current_node)
 
         return dag
 
@@ -241,10 +267,7 @@ def _error(self, circuit, qubit):
                 circuit, qubit, self.error_map
             )
         else:
-            circuit_list = [(x.op.name, []) for x in circuit]
-            return euler_one_qubit_decomposer.compute_error_list(
-                circuit_list, qubit, self.error_map
-            )
+            return euler_one_qubit_decomposer.compute_error_list(circuit, qubit, self.error_map)
 
 
 def _possible_decomposers(basis_set):
diff --git a/releasenotes/notes/oxidize-acg-0294a87c0d5974fa.yaml b/releasenotes/notes/oxidize-acg-0294a87c0d5974fa.yaml
index 6bd6761e0355..532d3e8fa55e 100644
--- a/releasenotes/notes/oxidize-acg-0294a87c0d5974fa.yaml
+++ b/releasenotes/notes/oxidize-acg-0294a87c0d5974fa.yaml
@@ -1,5 +1,5 @@
 ---
-upgrade_synthesis:
+features_synthesis:
   - |
     Port :func:`.synth_permutation_acg`, used to synthesize qubit permutations, to Rust.
     This produces an approximate 3x performance improvement on 1000 qubit circuits.
diff --git a/releasenotes/notes/oxidize-permbasic-be27578187ac472f.yaml b/releasenotes/notes/oxidize-permbasic-be27578187ac472f.yaml
index e770aa1ca31b..bd8c969be934 100644
--- a/releasenotes/notes/oxidize-permbasic-be27578187ac472f.yaml
+++ b/releasenotes/notes/oxidize-permbasic-be27578187ac472f.yaml
@@ -1,4 +1,4 @@
 ---
-upgrade_synthesis:
+features_synthesis:
   - |
     Port :func:`.synth_permutation_basic`, used to synthesize qubit permutations, to Rust.
diff --git a/releasenotes/notes/oxidize-pmh-ec74e4002510eaad.yaml b/releasenotes/notes/oxidize-pmh-ec74e4002510eaad.yaml
new file mode 100644
index 000000000000..b25adbec21f9
--- /dev/null
+++ b/releasenotes/notes/oxidize-pmh-ec74e4002510eaad.yaml
@@ -0,0 +1,17 @@
+---
+features_synthesis:
+  - |
+    Port :func:`.synth_cnot_full_pmh`, used to synthesize a linear function
+    into a CX network, to Rust. This produces approximately 44x speedup,
+    as measured on 100 qubit circuits. 
+  - |
+    The function :func:`.synth_cnot_full_pmh` now allows choosing the 
+    (heuristically) optimal ``section_size`` by setting it to ``None``. Then, a value is 
+    chosen which attempts to minimize the upper bound on the number of CX gates, that is
+    :math:`\alpha \log_2(n)` where :math:`n` is the number of qubits and :math:`\alpha \approx 0.56`.
+fixes:
+  - |
+    Fixed a bug in :func:`.synth_cnot_full_pmh` where providing a ``section_size`` that did 
+    not divide the number of qubits without remainder could lead to wrong results. 
+    Now any ``section_size`` (at most equal to the number of qubits) synthesizes the correct
+    circuit. For a (heuristically) optimal value, set ``section_size=None``.
diff --git a/releasenotes/notes/oxidize-synth-clifford-greedy-0739e9688bc4eedd.yaml b/releasenotes/notes/oxidize-synth-clifford-greedy-0739e9688bc4eedd.yaml
new file mode 100644
index 000000000000..ea492e29a7fc
--- /dev/null
+++ b/releasenotes/notes/oxidize-synth-clifford-greedy-0739e9688bc4eedd.yaml
@@ -0,0 +1,6 @@
+---
+features_synthesis:
+  - |
+    The function :func:`.synth_clifford_greedy` that synthesizes :class:`.Clifford` operators
+    was ported to Rust, leading to a significant increase in performance for all numbers of
+    qubits. For Cliffords over 50 qubits, the speedup is on the order of 1000 times.
diff --git a/releasenotes/notes/uniform-superposition-gate-3bd95ffdf05ef18c.yaml b/releasenotes/notes/uniform-superposition-gate-3bd95ffdf05ef18c.yaml
new file mode 100644
index 000000000000..6017979748ea
--- /dev/null
+++ b/releasenotes/notes/uniform-superposition-gate-3bd95ffdf05ef18c.yaml
@@ -0,0 +1,27 @@
+---
+features:
+  - |
+    Implemented :class:`.UniformSuperpositionGate` class, which allows
+    the creation of a uniform superposition state using 
+    the Shukla-Vedula algorithm. This feature facilitates the 
+    creation of quantum circuits that produce a uniform superposition 
+    state :math:`\frac{1}{\sqrt{M}} \sum_{j=0}^{M-1}  |j\rangle`, where 
+    :math:`M` is a positive integer representing the number of 
+    computational basis states with an amplitude of 
+    :math:`\frac{1}{\sqrt{M}}`. This implementation supports the 
+    efficient creation of uniform superposition states, 
+    requiring only :math:`O(\log_2 (M))` qubits and 
+    :math:`O(\log_2 (M))` gates. Usage example:
+
+    .. code-block:: python
+
+       from qiskit import QuantumCircuit 
+       from qiskit.circuit.library.data_preparation import UniformSuperpositionGate
+
+       M = 5
+       num_qubits = 3
+       usp_gate = UniformSuperpositionGate(M, num_qubits)
+       qc = QuantumCircuit(num_qubits)
+       qc.append(usp_gate, list(range(num_qubits)))
+
+       qc.draw()
diff --git a/test/python/circuit/test_gate_definitions.py b/test/python/circuit/test_gate_definitions.py
index 38bf7046cae5..c5df22a0e8a1 100644
--- a/test/python/circuit/test_gate_definitions.py
+++ b/test/python/circuit/test_gate_definitions.py
@@ -283,6 +283,7 @@ class TestGateEquivalenceEqual(QiskitTestCase):
         "ClassicalFunction",
         "ClassicalElement",
         "StatePreparation",
+        "UniformSuperpositionGate",
         "LinearFunction",
         "PermutationGate",
         "Commuting2qBlock",
diff --git a/test/python/circuit/test_uniform_superposition_gate.py b/test/python/circuit/test_uniform_superposition_gate.py
new file mode 100644
index 000000000000..019d6144ac5e
--- /dev/null
+++ b/test/python/circuit/test_uniform_superposition_gate.py
@@ -0,0 +1,98 @@
+# 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.
+
+"""
+Uniform Superposition Gate test.
+"""
+
+import unittest
+import math
+from test import QiskitTestCase
+import numpy as np
+from ddt import ddt, data
+
+from qiskit import QuantumCircuit
+from qiskit.quantum_info import Operator, Statevector
+
+from qiskit.circuit.library.data_preparation import (
+    UniformSuperpositionGate,
+)
+
+
+@ddt
+class TestUniformSuperposition(QiskitTestCase):
+    """Test initialization with UniformSuperpositionGate class"""
+
+    @data(2, 3, 5)
+    def test_uniform_superposition_gate(self, num_superpos_states):
+        """Test Uniform Superposition Gate"""
+        n = int(math.ceil(math.log2(num_superpos_states)))
+        desired_sv = (1 / np.sqrt(num_superpos_states)) * np.array(
+            [1.0] * num_superpos_states + [0.0] * (2**n - num_superpos_states)
+        )
+        gate = UniformSuperpositionGate(num_superpos_states, n)
+        actual_sv = Statevector(gate)
+        np.testing.assert_allclose(desired_sv, actual_sv)
+
+    @data(2, 3, 5, 13)
+    def test_inverse_uniform_superposition_gate(self, num_superpos_states):
+        """Test Inverse Uniform Superposition Gate"""
+        n = int(math.ceil(math.log2(num_superpos_states)))
+        gate = UniformSuperpositionGate(num_superpos_states, n)
+        qc = QuantumCircuit(n)
+        qc.append(gate, list(range(n)))
+        qc.append(gate.inverse(annotated=True), list(range(n)))
+        actual_unitary_matrix = np.array(Operator(qc).data)
+        desired_unitary_matrix = np.eye(2**n)
+        np.testing.assert_allclose(desired_unitary_matrix, actual_unitary_matrix, atol=1e-14)
+
+    @data(-2, -1, 0, 1)
+    def test_incompatible_num_superpos_states(self, num_superpos_states):
+        """Test error raised if num_superpos_states not valid"""
+        n = 1
+        with self.assertRaises(ValueError):
+            UniformSuperpositionGate(num_superpos_states, n)
+
+    @data(1, 2, 3, 4)
+    def test_incompatible_int_num_superpos_states_and_qubit_args(self, n):
+        """Test error raised if number of qubits not compatible with integer
+        state num_superpos_states (n >= log2(num_superpos_states) )"""
+        num_superpos_states = 50
+        with self.assertRaises(ValueError):
+            UniformSuperpositionGate(num_superpos_states, n)
+
+    @data(2, 3, 5)
+    def test_extra_qubits(self, num_superpos_states):
+        """Tests for cases where n >= log2(num_superpos_states)"""
+        num_extra_qubits = 2
+        n = int(math.ceil(math.log2(num_superpos_states))) + num_extra_qubits
+        desired_sv = (1 / np.sqrt(num_superpos_states)) * np.array(
+            [1.0] * num_superpos_states + [0.0] * (2**n - num_superpos_states)
+        )
+        gate = UniformSuperpositionGate(num_superpos_states, n)
+        actual_sv = Statevector(gate)
+        np.testing.assert_allclose(desired_sv, actual_sv)
+
+    @data(2, 3, 5)
+    def test_no_qubit_args(self, num_superpos_states):
+        """Test Uniform Superposition Gate without passing the number of qubits as an argument"""
+        n = int(math.ceil(math.log2(num_superpos_states)))
+        desired_sv = (1 / np.sqrt(num_superpos_states)) * np.array(
+            [1.0] * num_superpos_states + [0.0] * (2**n - num_superpos_states)
+        )
+        gate = UniformSuperpositionGate(num_superpos_states)
+        actual_sv = Statevector(gate)
+        np.testing.assert_allclose(desired_sv, actual_sv)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/python/synthesis/test_clifford_sythesis.py b/test/python/synthesis/test_clifford_sythesis.py
index 887f1af5ad99..8ca11f1ef251 100644
--- a/test/python/synthesis/test_clifford_sythesis.py
+++ b/test/python/synthesis/test_clifford_sythesis.py
@@ -16,6 +16,7 @@
 import numpy as np
 from ddt import ddt
 from qiskit.circuit.random import random_clifford_circuit
+from qiskit.quantum_info import random_clifford
 from qiskit.quantum_info.operators import Clifford
 from qiskit.synthesis.clifford import (
     synth_clifford_full,
@@ -99,8 +100,7 @@ def test_synth_greedy(self, num_qubits):
         rng = np.random.default_rng(1234)
         samples = 50
         for _ in range(samples):
-            circ = random_clifford_circuit(num_qubits, 5 * num_qubits, seed=rng)
-            target = Clifford(circ)
+            target = random_clifford(num_qubits, rng)
             synth_circ = synth_clifford_greedy(target)
             value = Clifford(synth_circ)
             self.assertEqual(value, target)
diff --git a/test/python/synthesis/test_linear_synthesis.py b/test/python/synthesis/test_linear_synthesis.py
index bbfab20a30fc..3aec57647b22 100644
--- a/test/python/synthesis/test_linear_synthesis.py
+++ b/test/python/synthesis/test_linear_synthesis.py
@@ -139,6 +139,42 @@ def test_synth_lnn_kms(self, num_qubits):
                 dist = abs(q0 - q1)
                 self.assertEqual(dist, 1)
 
+    @data(None, 2, 4, 5)
+    def test_synth_full_pmh(self, section_size):
+        """Test the PMH synthesis on pseudo-random matrices."""
+        num_qubits = 5
+        rng = np.random.default_rng(1234)
+        num_trials = 10
+        seeds = rng.integers(100000, size=num_trials, dtype=np.uint64)
+        for seed in seeds:
+            mat = random_invertible_binary_matrix(num_qubits, seed=seed)
+            mat = np.array(mat, dtype=bool)
+
+            qc = synth_cnot_count_full_pmh(mat, section_size)
+            self.assertEqual(LinearFunction(qc), LinearFunction(mat))
+
+    @data(5, 11)
+    def test_pmh_section_sizes(self, num_qubits):
+        """Test the PMH algorithm for different section sizes.
+
+        Regression test of Qiskit/qiskit#12166.
+        """
+        qc = QuantumCircuit(num_qubits)
+        for i in range(num_qubits - 1):
+            qc.cx(i, i + 1)
+        lin = LinearFunction(qc)
+
+        for section_size in range(1, num_qubits + 1):
+            synthesized = synth_cnot_count_full_pmh(lin.linear, section_size=section_size)
+            with self.subTest(section_size=section_size):
+                self.assertEqual(lin, LinearFunction(synthesized))
+
+    def test_pmh_invalid_section_size(self):
+        """Test that a section size larger than the number of qubits raises an error."""
+        mat = [[True, False], [True, False]]
+        with self.assertRaises(ValueError):
+            _ = synth_cnot_count_full_pmh(mat, section_size=3)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/utils/base.py b/test/utils/base.py
index ce9509709bad..e4dbd7f67072 100644
--- a/test/utils/base.py
+++ b/test/utils/base.py
@@ -14,7 +14,7 @@
 
 """Base TestCases for the unit tests.
 
-Implementors of unit tests for Terra are encouraged to subclass
+Implementors of unit tests for Qiskit should subclass
 ``QiskitTestCase`` in order to take advantage of utility functions (for example,
 the environment variables for customizing different options), and the
 decorators in the ``decorators`` package.
@@ -65,11 +65,97 @@ class BaseTestCase(unittest.TestCase):
 
 
 @enforce_subclasses_call(["setUp", "setUpClass", "tearDown", "tearDownClass"])
-class BaseQiskitTestCase(BaseTestCase):
-    """Additions for test cases for all Qiskit-family packages.
+class QiskitTestCase(BaseTestCase):
+    """Additions for Qiskit test cases."""
 
-    The additions here are intended for all packages, not just Terra.  Terra-specific logic should
-    be in the Terra-specific classes."""
+    @classmethod
+    def setUpClass(cls):
+        super().setUpClass()
+        # Set logging to file and stdout if the LOG_LEVEL envar is set.
+        cls.log = logging.getLogger(cls.__name__)
+
+        if log_level := os.getenv("LOG_LEVEL"):
+            log_fmt = f"{cls.log.name}.%(funcName)s:%(levelname)s:%(asctime)s: %(message)s"
+            formatter = logging.Formatter(log_fmt)
+            file_handler = logging.FileHandler(f"{os.path.splitext(inspect.getfile(cls))[0]}.log")
+            file_handler.setFormatter(formatter)
+            cls.log.addHandler(file_handler)
+
+            if os.getenv("STREAM_LOG"):
+                # Set up the stream handler.
+                stream_handler = logging.StreamHandler()
+                stream_handler.setFormatter(formatter)
+                cls.log.addHandler(stream_handler)
+
+            # Set the logging level from the environment variable, defaulting
+            # to INFO if it is not a valid level.
+            level = logging._nameToLevel.get(log_level, logging.INFO)
+            cls.log.setLevel(level)
+
+        warnings.filterwarnings("error", category=DeprecationWarning)
+        warnings.filterwarnings("error", category=QiskitWarning)
+
+        # Numpy 2 made a few new modules private, and have warnings that trigger if you try to
+        # access attributes that _would_ have existed.  Unfortunately, Python's `warnings` module
+        # adds a field called `__warningregistry__` to any module that triggers a warning, and
+        # `unittest.TestCase.assertWarns` then queries said fields on all existing modules.  On
+        # macOS ARM, we see some (we think harmless) warnings come out of `numpy.linalg._linalg` (a
+        # now-private module) during transpilation, which means that subsequent `assertWarns` calls
+        # can spuriously trick Numpy into sending out a nonsense `DeprecationWarning`.
+        # Tracking issue: https://github.com/Qiskit/qiskit/issues/12679
+        warnings.filterwarnings(
+            "ignore",
+            category=DeprecationWarning,
+            message=r".*numpy\.(\w+\.)*__warningregistry__",
+        )
+
+        # We only use pandas transitively through seaborn, so it's their responsibility to mark if
+        # their use of pandas would be a problem.
+        warnings.filterwarnings(
+            "default",
+            category=DeprecationWarning,
+            # The `(?s)` magic is to force use of the `re.DOTALL` flag, because the Pandas message
+            # includes hard-break newlines all over the place.
+            message="(?s).*Pyarrow.*required dependency.*next major release of pandas",
+            module=r"seaborn(\..*)?",
+        )
+
+        # Safe to remove once https://github.com/Qiskit/qiskit-aer/pull/2179 is in a release version
+        # of Aer.
+        warnings.filterwarnings(
+            "default",
+            category=DeprecationWarning,
+            message="Treating CircuitInstruction as an iterable is deprecated",
+            module=r"qiskit_aer(\.[a-zA-Z0-9_]+)*",
+        )
+
+        allow_DeprecationWarning_modules = [
+            "test.python.pulse.test_builder",
+            "test.python.pulse.test_block",
+            "test.python.quantum_info.operators.symplectic.test_legacy_pauli",
+            "qiskit.quantum_info.operators.pauli",
+            "pybobyqa",
+            "numba",
+            "qiskit.utils.measurement_error_mitigation",
+            "qiskit.circuit.library.standard_gates.x",
+            "qiskit.pulse.schedule",
+            "qiskit.pulse.instructions.instruction",
+            "qiskit.pulse.instructions.play",
+            "qiskit.pulse.library.parametric_pulses",
+            "qiskit.quantum_info.operators.symplectic.pauli",
+        ]
+        for mod in allow_DeprecationWarning_modules:
+            warnings.filterwarnings("default", category=DeprecationWarning, module=mod)
+        allow_DeprecationWarning_message = [
+            r"elementwise comparison failed.*",
+            r"The jsonschema validation included in qiskit-terra.*",
+            r"The DerivativeBase.parameter_expression_grad method.*",
+            r"The property ``qiskit\.circuit\.bit\.Bit\.(register|index)`` is deprecated.*",
+            # Caused by internal scikit-learn scipy usage
+            r"The 'sym_pos' keyword is deprecated and should be replaced by using",
+        ]
+        for msg in allow_DeprecationWarning_message:
+            warnings.filterwarnings("default", category=DeprecationWarning, message=msg)
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -99,6 +185,13 @@ def tearDown(self):
             )
         self.__teardown_called = True
 
+        # Reset the default providers, as in practice they acts as a singleton
+        # due to importing the instances from the top-level qiskit namespace.
+        from qiskit.providers.basic_provider import BasicProvider
+
+        with self.assertWarns(DeprecationWarning):
+            BasicProvider()._backends = BasicProvider()._verify_backends()
+
     def assertQuantumCircuitEqual(self, qc1, qc2, msg=None):
         """Extra assertion method to give a better error message when two circuits are unequal."""
         if qc1 == qc2:
@@ -165,112 +258,10 @@ def set_parallel_env(name, value):
         os.environ["QISKIT_PARALLEL"] = parallel_default
 
 
-class QiskitTestCase(BaseQiskitTestCase):
-    """Terra-specific extra functionality for test cases."""
-
-    def tearDown(self):
-        super().tearDown()
-        # Reset the default providers, as in practice they acts as a singleton
-        # due to importing the instances from the top-level qiskit namespace.
-        from qiskit.providers.basic_provider import BasicProvider
-
-        with self.assertWarns(DeprecationWarning):
-            BasicProvider()._backends = BasicProvider()._verify_backends()
-
-    @classmethod
-    def setUpClass(cls):
-        super().setUpClass()
-        # Set logging to file and stdout if the LOG_LEVEL envar is set.
-        cls.log = logging.getLogger(cls.__name__)
-
-        if log_level := os.getenv("LOG_LEVEL"):
-            log_fmt = f"{cls.log.name}.%(funcName)s:%(levelname)s:%(asctime)s: %(message)s"
-            formatter = logging.Formatter(log_fmt)
-            file_handler = logging.FileHandler(f"{os.path.splitext(inspect.getfile(cls))[0]}.log")
-            file_handler.setFormatter(formatter)
-            cls.log.addHandler(file_handler)
-
-            if os.getenv("STREAM_LOG"):
-                # Set up the stream handler.
-                stream_handler = logging.StreamHandler()
-                stream_handler.setFormatter(formatter)
-                cls.log.addHandler(stream_handler)
-
-            # Set the logging level from the environment variable, defaulting
-            # to INFO if it is not a valid level.
-            level = logging._nameToLevel.get(log_level, logging.INFO)
-            cls.log.setLevel(level)
-
-        warnings.filterwarnings("error", category=DeprecationWarning)
-        warnings.filterwarnings("error", category=QiskitWarning)
-
-        # Numpy 2 made a few new modules private, and have warnings that trigger if you try to
-        # access attributes that _would_ have existed.  Unfortunately, Python's `warnings` module
-        # adds a field called `__warningregistry__` to any module that triggers a warning, and
-        # `unittest.TestCase.assertWarns` then queries said fields on all existing modules.  On
-        # macOS ARM, we see some (we think harmless) warnings come out of `numpy.linalg._linalg` (a
-        # now-private module) during transpilation, which means that subsequent `assertWarns` calls
-        # can spuriously trick Numpy into sending out a nonsense `DeprecationWarning`.
-        # Tracking issue: https://github.com/Qiskit/qiskit/issues/12679
-        warnings.filterwarnings(
-            "ignore",
-            category=DeprecationWarning,
-            message=r".*numpy\.(\w+\.)*__warningregistry__",
-        )
-
-        # We only use pandas transitively through seaborn, so it's their responsibility to mark if
-        # their use of pandas would be a problem.
-        warnings.filterwarnings(
-            "default",
-            category=DeprecationWarning,
-            # The `(?s)` magic is to force use of the `re.DOTALL` flag, because the Pandas message
-            # includes hard-break newlines all over the place.
-            message="(?s).*Pyarrow.*required dependency.*next major release of pandas",
-            module=r"seaborn(\..*)?",
-        )
-
-        # Safe to remove once https://github.com/Qiskit/qiskit-aer/pull/2179 is in a release version
-        # of Aer.
-        warnings.filterwarnings(
-            "default",
-            category=DeprecationWarning,
-            message="Treating CircuitInstruction as an iterable is deprecated",
-            module=r"qiskit_aer(\.[a-zA-Z0-9_]+)*",
-        )
-
-        allow_DeprecationWarning_modules = [
-            "test.python.pulse.test_builder",
-            "test.python.pulse.test_block",
-            "test.python.quantum_info.operators.symplectic.test_legacy_pauli",
-            "qiskit.quantum_info.operators.pauli",
-            "pybobyqa",
-            "numba",
-            "qiskit.utils.measurement_error_mitigation",
-            "qiskit.circuit.library.standard_gates.x",
-            "qiskit.pulse.schedule",
-            "qiskit.pulse.instructions.instruction",
-            "qiskit.pulse.instructions.play",
-            "qiskit.pulse.library.parametric_pulses",
-            "qiskit.quantum_info.operators.symplectic.pauli",
-        ]
-        for mod in allow_DeprecationWarning_modules:
-            warnings.filterwarnings("default", category=DeprecationWarning, module=mod)
-        allow_DeprecationWarning_message = [
-            r"elementwise comparison failed.*",
-            r"The jsonschema validation included in qiskit-terra.*",
-            r"The DerivativeBase.parameter_expression_grad method.*",
-            r"The property ``qiskit\.circuit\.bit\.Bit\.(register|index)`` is deprecated.*",
-            # Caused by internal scikit-learn scipy usage
-            r"The 'sym_pos' keyword is deprecated and should be replaced by using",
-        ]
-        for msg in allow_DeprecationWarning_message:
-            warnings.filterwarnings("default", category=DeprecationWarning, message=msg)
-
-
 class FullQiskitTestCase(QiskitTestCase):
-    """Terra-specific further additions for test cases, if ``testtools`` is available.
+    """further additions for Qiskit test cases, if ``testtools`` is available.
 
-    It is not normally safe to derive from this class by name; on import, Terra checks if the
+    It is not normally safe to derive from this class by name; on import, Qiskit checks if the
     necessary packages are available, and binds this class to the name :obj:`~QiskitTestCase` if so.
     If you derive directly from it, you may try and instantiate the class without satisfying its
     dependencies."""