diff --git a/Cargo.lock b/Cargo.lock
index 67bf99d2a083..f973ba5f2a2e 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -2528,6 +2528,7 @@ dependencies = [
  "rustworkx-core",
  "smallvec",
  "thiserror 2.0.18",
+ "thread_local",
  "uuid",
 ]
 
diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs
index 0b70641b36cc..c8e8bcfa71e8 100644
--- a/crates/pyext/src/lib.rs
+++ b/crates/pyext/src/lib.rs
@@ -76,6 +76,11 @@ fn _accelerate(m: &Bound<PyModule>) -> PyResult<()> {
     add_submodule(m, ::qiskit_transpiler::passes::high_level_synthesis_mod, "high_level_synthesis")?;
     add_submodule(m, ::qiskit_transpiler::passes::remove_diagonal_gates_before_measure_mod, "remove_diagonal_gates_before_measure")?;
     add_submodule(m, ::qiskit_transpiler::passes::remove_identity_equiv_mod, "remove_identity_equiv")?;
+    add_submodule(
+        m,
+        ::qiskit_transpiler::passes::two_qubit_peephole_mod,
+        "two_qubit_peephole",
+    )?;
     add_submodule(m, ::qiskit_accelerate::results::results, "results")?;
     add_submodule(m, ::qiskit_transpiler::passes::sabre::sabre, "sabre")?;
     add_submodule(m, ::qiskit_accelerate::sampled_exp_val::sampled_exp_val, "sampled_exp_val")?;
diff --git a/crates/synthesis/src/qsd.rs b/crates/synthesis/src/qsd.rs
index ce6a3ace9a95..9d2fb688ff35 100644
--- a/crates/synthesis/src/qsd.rs
+++ b/crates/synthesis/src/qsd.rs
@@ -109,7 +109,7 @@ pub fn quantum_shannon_decomposition(
         smallvec![],
         aview2(&qiskit_circuit::gate_matrix::CX_GATE),
         1.0,
-        "U",
+        EulerBasis::U,
         None,
     )?;
     let one_qubit_decomposer = one_qubit_decomposer_basis_set.unwrap_or(&default_1q_basis);
diff --git a/crates/synthesis/src/two_qubit_decompose/basis_decomposer.rs b/crates/synthesis/src/two_qubit_decompose/basis_decomposer.rs
index 42db5797fe9a..86220b3955dd 100644
--- a/crates/synthesis/src/two_qubit_decompose/basis_decomposer.rs
+++ b/crates/synthesis/src/two_qubit_decompose/basis_decomposer.rs
@@ -10,6 +10,8 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
+use std::str::FromStr;
+
 use approx::{abs_diff_eq, relative_eq};
 use num_complex::{Complex64, ComplexFloat};
 use num_traits::Zero;
@@ -139,6 +141,12 @@ impl TwoQubitBasisDecomposer {
         __num_basis_gates(self.basis_decomposer.b, self.basis_fidelity, u)
     }
 
+    /// Is the gate super controlled
+    #[inline]
+    pub fn super_controlled(&self) -> bool {
+        self.super_controlled
+    }
+
     fn decomp1_inner(
         &self,
         target: &TwoQubitWeylDecomposition,
@@ -524,7 +532,7 @@ impl TwoQubitBasisDecomposer {
         gate_params: SmallVec<[f64; 3]>,
         gate_matrix: ArrayView2<Complex64>,
         basis_fidelity: f64,
-        euler_basis: &str,
+        euler_basis: EulerBasis,
         pulse_optimize: Option<bool>,
     ) -> PyResult<Self> {
         let basis_decomposer =
@@ -610,7 +618,7 @@ impl TwoQubitBasisDecomposer {
             gate,
             gate_params,
             basis_fidelity,
-            euler_basis: EulerBasis::__new__(euler_basis)?,
+            euler_basis,
             pulse_optimize,
             basis_decomposer,
             super_controlled,
@@ -787,7 +795,7 @@ impl TwoQubitBasisDecomposer {
             gate_params?,
             gate_matrix.as_array(),
             basis_fidelity,
-            euler_basis,
+            EulerBasis::from_str(euler_basis).map_err(PyValueError::new_err)?,
             pulse_optimize,
         )
     }
@@ -1020,7 +1028,7 @@ pub fn two_qubit_decompose_up_to_diagonal(
         smallvec![],
         aview2(&CX_GATE),
         1.0,
-        "U",
+        EulerBasis::U,
         None,
     )?;
 
diff --git a/crates/synthesis/src/two_qubit_decompose/gate_sequence.rs b/crates/synthesis/src/two_qubit_decompose/gate_sequence.rs
index cc170d4386a4..8d8d79e4b49b 100644
--- a/crates/synthesis/src/two_qubit_decompose/gate_sequence.rs
+++ b/crates/synthesis/src/two_qubit_decompose/gate_sequence.rs
@@ -17,8 +17,8 @@ pub(crate) type TwoQubitSequenceVec = Vec<(PackedOperation, SmallVec<[f64; 3]>,
 
 #[derive(Clone, Debug)]
 pub struct TwoQubitGateSequence {
-    pub(crate) gates: TwoQubitSequenceVec,
-    pub(crate) global_phase: f64,
+    pub gates: TwoQubitSequenceVec,
+    pub global_phase: f64,
 }
 
 impl TwoQubitGateSequence {
diff --git a/crates/transpiler/Cargo.toml b/crates/transpiler/Cargo.toml
index 6221abb5f06b..f244782059fd 100644
--- a/crates/transpiler/Cargo.toml
+++ b/crates/transpiler/Cargo.toml
@@ -31,6 +31,7 @@ thiserror.workspace = true
 bytemuck.workspace = true
 fixedbitset.workspace = true
 anyhow.workspace = true
+thread_local = "1.1"
 
 [dependencies.uuid]
 workspace = true
diff --git a/crates/transpiler/src/passes/consolidate_blocks.rs b/crates/transpiler/src/passes/consolidate_blocks.rs
index 50c5b6178a9f..67b9e461ef26 100644
--- a/crates/transpiler/src/passes/consolidate_blocks.rs
+++ b/crates/transpiler/src/passes/consolidate_blocks.rs
@@ -31,6 +31,7 @@ use qiskit_circuit::interner::Interned;
 use qiskit_circuit::operations::StandardGate;
 use qiskit_circuit::operations::{ArrayType, Operation, Param, UnitaryGate};
 use qiskit_circuit::packed_instruction::PackedOperation;
+use qiskit_synthesis::euler_one_qubit_decomposer::EulerBasis;
 use qiskit_synthesis::linalg::nalgebra_array_view;
 use qiskit_synthesis::matrix::two_qubit::{
     blocks_to_matrix, get_1q_matrix_from_inst, get_2q_matrix_from_inst, get_matrix_from_inst,
@@ -145,7 +146,7 @@ fn get_decomposer_and_basis_gate(
                         SmallVec::default(),
                         get_matrix(&gate),
                         approximation_degree,
-                        "U",
+                        EulerBasis::U,
                         None,
                     )
                     .unwrap_or_else(|_| {
@@ -167,7 +168,7 @@ fn get_decomposer_and_basis_gate(
                 SmallVec::default(),
                 aview2(&CX_GATE),
                 1.0,
-                "U",
+                EulerBasis::U,
                 None,
             )
             .expect("Error while creating TwoQubitBasisDecomposer using a 'cx' gate."),
diff --git a/crates/transpiler/src/passes/high_level_synthesis.rs b/crates/transpiler/src/passes/high_level_synthesis.rs
index 750d8f890a6d..bab207ada29f 100644
--- a/crates/transpiler/src/passes/high_level_synthesis.rs
+++ b/crates/transpiler/src/passes/high_level_synthesis.rs
@@ -769,7 +769,7 @@ fn extract_definition(op: &PackedOperation, params: &[Param]) -> PyResult<Option
                         SmallVec::new(),
                         aview2(&CX_GATE),
                         1.0,
-                        "U",
+                        EulerBasis::U,
                         None,
                     )?;
                     let two_qubit_sequence =
diff --git a/crates/transpiler/src/passes/mod.rs b/crates/transpiler/src/passes/mod.rs
index 64c71b75fbe8..d1432e773b17 100644
--- a/crates/transpiler/src/passes/mod.rs
+++ b/crates/transpiler/src/passes/mod.rs
@@ -50,6 +50,7 @@ mod schedule_analysis;
 mod split_2q_unitaries;
 mod substitute_pi4_rotations;
 mod synthesize_rz_rotations;
+mod two_qubit_peephole;
 pub mod unitary_synthesis;
 mod unroll_3q_or_more;
 pub mod vf2;
@@ -104,6 +105,10 @@ pub use schedule_analysis::scheduling_mod;
 pub use split_2q_unitaries::{run_split_2q_unitaries, split_2q_unitaries_mod};
 pub use substitute_pi4_rotations::{py_run_substitute_pi4_rotations, substitute_pi4_rotations_mod};
 pub use synthesize_rz_rotations::{py_run_synthesize_rz_rotations, synthesize_rz_rotations_mod};
+pub use two_qubit_peephole::{
+    py_two_qubit_unitary_peephole_optimize, two_qubit_peephole_mod,
+    two_qubit_unitary_peephole_optimize,
+};
 pub use unitary_synthesis::{
     UnitarySynthesisConfig, UnitarySynthesisState, run_unitary_synthesis, unitary_synthesis_mod,
 };
diff --git a/crates/transpiler/src/passes/two_qubit_peephole.rs b/crates/transpiler/src/passes/two_qubit_peephole.rs
new file mode 100644
index 000000000000..1049db006561
--- /dev/null
+++ b/crates/transpiler/src/passes/two_qubit_peephole.rs
@@ -0,0 +1,371 @@
+// 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 https://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 std::cell::RefCell;
+use std::sync::Mutex;
+#[cfg(feature = "cache_pygates")]
+use std::sync::OnceLock;
+use std::sync::atomic::AtomicBool;
+
+use nalgebra::U4;
+use num_complex::Complex64;
+use pyo3::Python;
+use pyo3::intern;
+use pyo3::prelude::*;
+use rayon::prelude::*;
+use rustworkx_core::petgraph::algo::toposort;
+use rustworkx_core::petgraph::stable_graph::NodeIndex;
+use rustworkx_core::petgraph::visit::NodeIndexable;
+
+use qiskit_circuit::dag_circuit::{DAGCircuit, NodeType};
+use qiskit_circuit::instruction::Parameters;
+use qiskit_circuit::operations::{
+    Operation, OperationRef, Param, PyOperationTypes, PythonOperation,
+};
+use qiskit_circuit::packed_instruction::{PackedInstruction, PackedOperation};
+use qiskit_circuit::{BlocksMode, Qubit, VarsMode};
+
+use super::unitary_synthesis::Direction2q;
+use crate::passes::unitary_synthesis::{
+    Approximation, QpuConstraint, TwoQSynthesisResult, fidelity_2q_sequence, synthesize_2q_matrix,
+};
+use crate::passes::{UnitarySynthesisConfig, UnitarySynthesisState};
+use crate::target::Target;
+use qiskit_circuit::PhysicalQubit;
+use qiskit_synthesis::linalg::nalgebra_array_view;
+use qiskit_synthesis::matrix::two_qubit::blocks_to_matrix;
+use qiskit_synthesis::two_qubit_decompose::TwoQubitGateSequence;
+use qiskit_util::getenv_use_multiple_threads;
+use thread_local::ThreadLocal;
+
+type MappingIterItem = Option<(TwoQSynthesisResult<f64>, [Qubit; 2])>;
+
+/// A python entry-point to the pass function
+///
+/// This function explicitly releases the GIL prior to entering the parallel portion of the pass.
+/// This is necessary because if there are any Python defined and owned objects in the circuit
+/// the pass will need GIL access to interact with that object in parallel.
+#[pyfunction(name = "two_qubit_unitary_peephole_optimize")]
+pub fn py_two_qubit_unitary_peephole_optimize(
+    py: Python,
+    dag: &DAGCircuit,
+    target: &Target,
+    approximation_degree: Option<f64>,
+) -> PyResult<Option<DAGCircuit>> {
+    let result = py.detach(move || {
+        two_qubit_unitary_peephole_optimize_analysis(dag, target, approximation_degree)
+    })?;
+    let Some(result) = result else {
+        return Ok(None);
+    };
+    two_qubit_unitary_peephole_optimize_apply(dag, result)
+}
+
+/// A non-python entry-point to the pass function.
+///
+/// This function is not safe in the context where Python owned objects are in the circuit.
+/// It will hang/deadlock on the GIL if called in these contexts and should not be used if
+/// there are Python owned objects in the circuit. If you're using this from python you should call
+/// `py_two_qubit_unitary_peephole_optimize` instead.
+pub fn two_qubit_unitary_peephole_optimize(
+    dag: &DAGCircuit,
+    target: &Target,
+    approximation_degree: Option<f64>,
+) -> PyResult<Option<DAGCircuit>> {
+    let result = two_qubit_unitary_peephole_optimize_analysis(dag, target, approximation_degree)?;
+    let Some(result) = result else {
+        return Ok(None);
+    };
+    two_qubit_unitary_peephole_optimize_apply(dag, result)
+}
+
+fn score_sequence(
+    dir: &Direction2q,
+    sequence: &TwoQubitGateSequence,
+    constraint: &QpuConstraint,
+    qargs: [PhysicalQubit; 2],
+) -> (i64, f64, i64) {
+    let fidelity = fidelity_2q_sequence(dir, sequence, constraint, qargs);
+    // Make the gate counts negative because synthesize_2q_matrix picks the largest value
+    // we want to minimize the gate counts.
+    let gate_count = -(sequence.gates.len() as i64);
+    let twoq_gate_count = -(sequence
+        .gates
+        .iter()
+        .filter(|x| x.0.num_qubits() == 2)
+        .count() as i64);
+    (twoq_gate_count, fidelity, gate_count)
+}
+
+struct PeepholeResult {
+    run_mapping: Vec<MappingIterItem>,
+    node_mapping: Vec<usize>,
+}
+
+fn two_qubit_unitary_peephole_optimize_analysis(
+    dag: &DAGCircuit,
+    target: &Target,
+    approximation_degree: Option<f64>,
+) -> PyResult<Option<PeepholeResult>> {
+    let runs: Vec<Vec<NodeIndex>> = dag.collect_2q_runs().unwrap();
+    if runs.is_empty() {
+        return Ok(None);
+    }
+    let node_mapping: Vec<usize> = vec![usize::MAX; dag.dag().node_bound()];
+    let locked_node_mapping = Mutex::new(node_mapping);
+    let substitution_made = AtomicBool::new(false);
+    let physical_qubits = (0..dag.num_qubits() as u32)
+        .map(PhysicalQubit::new)
+        .collect::<Vec<_>>();
+    let approximation = Approximation::from_py_approximation_degree(approximation_degree);
+    let unitary_synthesis_config = UnitarySynthesisConfig {
+        approximation,
+        ..Default::default()
+    };
+    let thread_local_states = ThreadLocal::new();
+    let find_best_sequence =
+        |run_index: usize, node_indices: &[NodeIndex]| -> PyResult<MappingIterItem> {
+            let q_virt = node_indices
+                .iter()
+                .find_map(|node_index| {
+                    let inst = dag.dag()[*node_index].unwrap_operation();
+                    let qubits = dag.get_qargs(inst.qubits);
+                    if qubits.len() == 2 {
+                        Some([qubits[0], qubits[1]])
+                    } else {
+                        None
+                    }
+                })
+                .unwrap();
+            let q_phys = q_virt.map(|q| physical_qubits[q.index()]);
+            let matrix = blocks_to_matrix(dag, node_indices, q_virt)?;
+            let synthesis_state: &RefCell<UnitarySynthesisState> = thread_local_states
+                .get_or(|| RefCell::new(UnitarySynthesisState::new(unitary_synthesis_config)));
+
+            let result = synthesize_2q_matrix(
+                nalgebra_array_view::<Complex64, U4, U4>(matrix.as_view()).into(),
+                q_phys,
+                &mut synthesis_state.borrow_mut(),
+                QpuConstraint::Target(target),
+                score_sequence,
+            )?;
+            if result.is_none() {
+                return Ok(None);
+            }
+            let result = result.unwrap();
+            let mut original_fidelity: f64 = 1.;
+            let mut original_2q_count: usize = 0;
+            let original_total_count: usize = node_indices.len();
+            let mut outside_target = false;
+            for node_index in node_indices {
+                let NodeType::Operation(ref inst) = dag.dag()[*node_index] else {
+                    unreachable!("All run nodes will be ops")
+                };
+                let qubits: &[_] = match dag.get_qargs(inst.qubits) {
+                    [q] => &[PhysicalQubit(q.0)],
+                    [q0, q1] => &[PhysicalQubit(q0.0), PhysicalQubit(q1.0)],
+                    _ => panic!("Runs should only contain 1q and 2q gates"),
+                };
+                if qubits.len() == 2 {
+                    original_2q_count += 1;
+                }
+                let name = inst.op.name();
+                let gate_fidelity = match target.get_error(name, qubits) {
+                    Some(err) => 1. - err,
+                    None => {
+                        // If error rate is None this can mean either the gate is not supported
+                        // in the target or the gate is ideal. We need to do a second lookup
+                        // to determine if the gate is supported, and if it isn't we don't need
+                        // to finish scoring because we know we'll use the synthesis output
+                        if !target.instruction_supported(name, qubits, inst.params_view(), true) {
+                            outside_target = true;
+                            break;
+                        }
+                        1.
+                    }
+                };
+                original_fidelity *= gate_fidelity;
+            }
+            let original_score = (
+                original_2q_count,
+                1. - original_fidelity,
+                original_total_count,
+            );
+            let new_2q_count = result
+                .score
+                .map(|score| -score.0 as usize)
+                .unwrap_or_else(|| {
+                    result
+                        .sequence
+                        .gates
+                        .iter()
+                        .filter(|x| x.0.num_qubits() == 2)
+                        .count()
+                });
+            let new_gate_count = result.sequence.gates.len();
+            let new_score = (
+                new_2q_count,
+                1. - result.score.map(|score| score.1).unwrap_or_else(|| {
+                    fidelity_2q_sequence(
+                        &result.dir,
+                        &result.sequence,
+                        &QpuConstraint::Target(target),
+                        q_phys,
+                    )
+                }),
+                new_gate_count,
+            );
+            // If we are not outside the target and the new score isn't any better just use the
+            // original (this includes a tie).
+            if !outside_target && new_score >= original_score {
+                return Ok(None);
+            }
+            // This is done at the end of the map in some attempt to minimize
+            // lock contention. If this were serial code it'd make more sense
+            // to do this as part of the iteration building the
+            let mut node_mapping = locked_node_mapping.lock().unwrap();
+            for node in node_indices {
+                node_mapping[node.index()] = run_index;
+            }
+            substitution_made.store(true, std::sync::atomic::Ordering::Relaxed);
+            let result = TwoQSynthesisResult {
+                sequence: result.sequence,
+                dir: result.dir,
+                score: result.score.map(|score| score.1),
+            };
+            Ok(Some((result, q_virt)))
+        };
+
+    let run_mapping: PyResult<Vec<MappingIterItem>> = if getenv_use_multiple_threads() {
+        // Build a vec of all the best synthesized two qubit gate sequences from the collected runs.
+        // This is done in parallel
+        runs.par_iter()
+            .enumerate()
+            .map(|(index, sequence)| find_best_sequence(index, sequence.as_slice()))
+            .collect()
+    } else {
+        runs.iter()
+            .enumerate()
+            .map(|(index, sequence)| find_best_sequence(index, sequence.as_slice()))
+            .collect()
+    };
+    let run_mapping = run_mapping?;
+    if !substitution_made.into_inner() {
+        return Ok(None);
+    }
+    Ok(Some(PeepholeResult {
+        node_mapping: locked_node_mapping.into_inner().unwrap(),
+        run_mapping,
+    }))
+}
+
+/// This function runs the two qubit unitary peephole optimization pass
+///
+/// It returns None if there is no modifications/optimiations made to the input dag and the pass
+/// function calling this should just return the input dag from the pass.
+fn two_qubit_unitary_peephole_optimize_apply(
+    dag: &DAGCircuit,
+    result: PeepholeResult,
+) -> PyResult<Option<DAGCircuit>> {
+    // After we've computed all the sequences to execute now serially build up a new dag.
+    let mut processed_runs: Vec<bool> = vec![false; result.run_mapping.len()];
+    let out_dag = dag.copy_empty_like_with_same_capacity(VarsMode::Alike, BlocksMode::Keep)?;
+    let mut out_dag_builder = out_dag.into_builder();
+    for node in toposort(dag.dag(), None).unwrap() {
+        if !matches!(dag.dag()[node], NodeType::Operation(_)) {
+            continue;
+        }
+        let run_index = result.node_mapping[node.index()];
+        if run_index != usize::MAX {
+            if processed_runs[run_index] {
+                continue;
+            }
+            // If this is not a two qubit gate then there is a chance this will cause the
+            // insertion to happen too early. We skip the nodes in a run until we encounter
+            // a 2q gate which ensure the block is inserted into the correct location in the
+            // circuit.
+            if dag.dag()[node].unwrap_operation().op.num_qubits() != 2 {
+                continue;
+            }
+            // A None is inserted into the run_mapping as the value for a run that we don't
+            // substitute but was identified so we added an explicit None to preserve the
+            // indexing with the vec. This shouldn't be possible to hit the else condition
+            // since node mapping will never contain a value for a run_mapping index that
+            // is set to None.
+            let Some((result, qargs_virt)) = result.run_mapping[run_index].as_ref() else {
+                unreachable!(
+                    "node_mapping can't contain a value pointing to an unpoluated run in run_mapping"
+                );
+            };
+            let order = result.dir.as_indices();
+            let out_qargs = [qargs_virt[order[0] as usize], qargs_virt[order[1] as usize]];
+            let qubit_keys = [
+                out_dag_builder.insert_qargs(&[out_qargs[0]]),
+                out_dag_builder.insert_qargs(&[out_qargs[1]]),
+                out_dag_builder.insert_qargs(&[out_qargs[0], out_qargs[1]]),
+                out_dag_builder.insert_qargs(&[out_qargs[1], out_qargs[0]]),
+            ];
+            for (gate, params, local_qubits) in &result.sequence.gates {
+                let qubits = match local_qubits.as_slice() {
+                    [0] => qubit_keys[0],
+                    [1] => qubit_keys[1],
+                    [0, 1] => qubit_keys[2],
+                    [1, 0] => qubit_keys[3],
+                    _ => panic!(
+                        "internal logic error: decomposed sequence contained unexpected qargs"
+                    ),
+                };
+                let op = match gate.view() {
+                    OperationRef::StandardGate(gate) => PackedOperation::from(gate),
+                    OperationRef::Gate(py_gate) => Python::attach(|py| -> PyResult<_> {
+                        let gate = py_gate.py_copy(py)?;
+                        gate.instruction
+                            .setattr(py, intern!(py, "params"), params)?;
+                        Ok(PackedOperation::from(Box::new(PyOperationTypes::Gate(
+                            gate,
+                        ))))
+                    })?,
+                    _ => {
+                        panic!("internal logic error: decomposed sequence contains a non-gate")
+                    }
+                };
+                let params = (!params.is_empty()).then(|| {
+                    Box::new(Parameters::Params(
+                        params.iter().copied().map(Param::Float).collect(),
+                    ))
+                });
+                out_dag_builder.push_back(PackedInstruction {
+                  op,
+                  qubits,
+                  clbits: Default::default(),
+                  params,
+                  label: None,
+                  #[cfg(feature = "cache_pygates")] // W: code is inactive due to #[cfg] directives: feature …
+                  py_op: OnceLock::new(),
+                })?;
+            }
+            out_dag_builder.add_global_phase(&Param::Float(result.sequence.global_phase()))?;
+            processed_runs[run_index] = true;
+        } else {
+            let NodeType::Operation(ref instr) = dag.dag()[node] else {
+                unreachable!("Must be an op node")
+            };
+            out_dag_builder.push_back(instr.clone())?;
+        }
+    }
+    Ok(Some(out_dag_builder.build()))
+}
+
+pub fn two_qubit_peephole_mod(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(py_two_qubit_unitary_peephole_optimize))?;
+    Ok(())
+}
diff --git a/crates/transpiler/src/passes/unitary_synthesis/decomposers.rs b/crates/transpiler/src/passes/unitary_synthesis/decomposers.rs
index b16d527dcd13..e7535bf6ecae 100644
--- a/crates/transpiler/src/passes/unitary_synthesis/decomposers.rs
+++ b/crates/transpiler/src/passes/unitary_synthesis/decomposers.rs
@@ -94,7 +94,18 @@ impl hash::Hash for ContinuousPauliConstructor {
 }
 impl ContinuousPauliConstructor {
     fn try_build(&self) -> PyResult<TwoQubitControlledUDecomposer> {
-        TwoQubitControlledUDecomposer::new_inner(self.source.clone(), self.euler.as_str())
+        match self.source {
+            RXXEquivalent::Standard(gate) => TwoQubitControlledUDecomposer::new_inner(
+                RXXEquivalent::Standard(gate),
+                self.euler.as_str(),
+            ),
+            RXXEquivalent::CustomPython(ref gate) => Python::attach(|py| {
+                TwoQubitControlledUDecomposer::new_inner(
+                    RXXEquivalent::CustomPython(gate.clone_ref(py)),
+                    self.euler.as_str(),
+                )
+            }),
+        }
     }
 }
 
@@ -203,7 +214,7 @@ impl StaticKakConstructor {
             self.source.params.clone(),
             matrix.view(),
             self.fidelity.get(),
-            self.euler.as_str(),
+            self.euler,
             self.use_pulse_optimizer.to_py_pulse_optimize(),
         )
     }
@@ -650,9 +661,15 @@ fn get_2q_decomposers(
                 // _all_ of the 1q bases simultaneously without further decompositions, but don't
                 // expose that functionality.  This wastes huge amounts of time and needs a fix.
                 for euler in euler_bases.get_bases() {
+                    let rxx_copy = match rxx_equivalent {
+                        RXXEquivalent::Standard(gate) => RXXEquivalent::Standard(gate),
+                        RXXEquivalent::CustomPython(ref gate) => {
+                            RXXEquivalent::CustomPython(Python::attach(|py| gate.clone_ref(py)))
+                        }
+                    };
                     let constructor =
                         Decomposer2qConstructor::ContinuousPauli(ContinuousPauliConstructor {
-                            source: rxx_equivalent.clone(),
+                            source: rxx_copy,
                             euler,
                         });
                     let flip = choose_flip(candidate.direction, &constructor);
diff --git a/crates/transpiler/src/passes/unitary_synthesis/mod.rs b/crates/transpiler/src/passes/unitary_synthesis/mod.rs
index 5021009b9d24..c84a9deb4e21 100644
--- a/crates/transpiler/src/passes/unitary_synthesis/mod.rs
+++ b/crates/transpiler/src/passes/unitary_synthesis/mod.rs
@@ -23,7 +23,9 @@ use numpy::PyReadonlyArray2;
 use pyo3::prelude::*;
 use pyo3::{intern, wrap_pyfunction};
 
-use self::decomposers::{Decomposer2q, DecomposerCache, Direction2q, FlipDirection};
+pub(crate) use self::decomposers::Direction2q;
+
+use self::decomposers::{Decomposer2q, DecomposerCache, FlipDirection};
 use crate::QiskitError;
 use crate::target::Target;
 use qiskit_circuit::bit::QuantumRegister;
@@ -486,14 +488,79 @@ fn synthesize_1q_matrix_onto(
     Ok(true)
 }
 
-fn synthesize_2q_matrix_onto(
-    out: &mut DAGCircuitBuilder,
+#[derive(Debug)]
+pub struct TwoQSynthesisResult<S> {
+    pub sequence: TwoQubitGateSequence,
+    pub dir: Direction2q,
+    pub score: Option<S>,
+}
+
+/// Estimate the fidelity of a synthesized two qubit unitary synthesis output
+///
+/// The estimated fidelity is computed by taking the error rate from the target
+/// and computing the product of the fidelities of each gate in the circuit. In
+/// the absence of error rates (either QpuConstraint doesn't have a target or
+/// the target contains no error rates) then a fidelity of 1 is returned.
+///
+/// # Args
+/// `dir` - The direction of the synthesis (forward or reverse)
+/// `sequence` - The synthesis sequence
+/// `constraint` - The qpu constraints used for the synthesis, typically just a Target
+/// `qargs_phys` - The qpu qargs the unitary is run on
+#[inline]
+pub(crate) fn fidelity_2q_sequence(
+    dir: &Direction2q,
+    sequence: &TwoQubitGateSequence,
+    constraint: &QpuConstraint,
+    qargs_phys: [PhysicalQubit; 2],
+) -> f64 {
+    let QpuConstraint::Target(target) = &constraint else {
+        return 1.;
+    };
+    let order = dir.as_indices();
+    let phys = [qargs_phys[order[0] as usize], qargs_phys[order[1] as usize]];
+    sequence
+        .gates()
+        .iter()
+        .map(|(op, _, qubits)| {
+            let qargs: &[_] = match *qubits.as_slice() {
+                [q] => &[phys[q as usize]],
+                [q1, q2] => &[phys[q1 as usize], phys[q2 as usize]],
+                _ => panic!("sequences should only contain 1q and 2q gates"),
+            };
+            // TODO: this does not handle the possibility of a 2q decomposer (like the
+            // XXDecomposer) using specialised instructions whose operation names do not match
+            // their target key.
+            1. - target.get_error(op.name(), qargs).unwrap_or(0.)
+        })
+        .product()
+}
+
+/// Synthesize a given two qubit unitary matrix into a gate sequence.
+///
+/// For overcomplete targets this will run all compatible decomposers and pick the output synthesis
+/// with the highest estimated fidelity.
+///
+/// # Args
+/// `unitary` - The unitary to synthesize
+/// `qargs_phys` - The physical qubits the unitary is being applied to
+/// `state` - The internal state of the pass, this includes the configured synthesizers
+/// `constraint` - The qpu constraints used for the synthesis, typically just a Target
+/// `fidelity_calculation` - A callable that is called with a two qubit synthesis output sequence and
+///   the context around it. It is expected to return a fidelity score to compare that synthesis output
+///   against other decomposers. The synthesis output with the maximum score is selected and is
+///   what is returned by the `synthesize_2q_matrix`.
+pub(crate) fn synthesize_2q_matrix<F, S>(
     mut unitary: CowArray<Complex64, Ix2>,
     qargs_phys: [PhysicalQubit; 2],
-    qargs_virt: [Qubit; 2],
     state: &mut UnitarySynthesisState,
     constraint: QpuConstraint,
-) -> PyResult<bool> {
+    mut fidelity_calculation: F,
+) -> PyResult<Option<TwoQSynthesisResult<S>>>
+where
+    F: FnMut(&Direction2q, &TwoQubitGateSequence, &QpuConstraint, [PhysicalQubit; 2]) -> S,
+    S: PartialOrd,
+{
     let decomposer_cache = &mut state.cache;
     let config = &state.config;
 
@@ -546,31 +613,7 @@ fn synthesize_2q_matrix_onto(
         // inconsistent; either it should be an error in _all_ circumstances if synthesis fails or
         // in _none_.  It's tricky to recreate the pre-Qiskit-2.4 behaviour bug-for-bug in the new
         // refactor because of how the split between decomposer construction and use works now.
-        return Ok(false);
-    };
-
-    let fidelity = |pair: &(Direction2q, TwoQubitGateSequence)| -> f64 {
-        let QpuConstraint::Target(target) = &constraint else {
-            return 1.;
-        };
-        let (dir, sequence) = pair;
-        let order = dir.as_indices();
-        let phys = [qargs_phys[order[0] as usize], qargs_phys[order[1] as usize]];
-        sequence
-            .gates()
-            .iter()
-            .map(|(op, _, qubits)| {
-                let qargs: &[_] = match *qubits.as_slice() {
-                    [q] => &[phys[q as usize]],
-                    [q1, q2] => &[phys[q1 as usize], phys[q2 as usize]],
-                    _ => panic!("sequences should only contain 1q and 2q gates"),
-                };
-                // TODO: this does not handle the possibility of a 2q decomposer (like the
-                // XXDecomposer) using specialised instructions whose operation names do not match
-                // their target key.
-                1. - target.get_error(op.name(), qargs).unwrap_or(0.)
-            })
-            .product()
+        return Ok(None);
     };
 
     // We only need to calculate the best score if there's more than one sequence.
@@ -578,8 +621,10 @@ fn synthesize_2q_matrix_onto(
     let mut best_pair = first;
     for sequence in sequences {
         let sequence = sequence?;
-        let prev_fidelity = best_fidelity.unwrap_or_else(|| fidelity(&best_pair));
-        let this_fidelity = fidelity(&sequence);
+        let prev_fidelity = best_fidelity.unwrap_or_else(|| {
+            fidelity_calculation(&best_pair.0, &best_pair.1, &constraint, qargs_phys)
+        });
+        let this_fidelity = fidelity_calculation(&sequence.0, &sequence.1, &constraint, qargs_phys);
         if this_fidelity > prev_fidelity {
             best_fidelity = Some(this_fidelity);
             best_pair = sequence;
@@ -587,10 +632,28 @@ fn synthesize_2q_matrix_onto(
             best_fidelity = Some(prev_fidelity);
         }
     }
+    Ok(Some(TwoQSynthesisResult {
+        sequence: best_pair.1,
+        dir: best_pair.0,
+        score: best_fidelity,
+    }))
+}
 
+fn synthesize_2q_matrix_onto(
+    out: &mut DAGCircuitBuilder,
+    unitary: CowArray<Complex64, Ix2>,
+    qargs_phys: [PhysicalQubit; 2],
+    qargs_virt: [Qubit; 2],
+    state: &mut UnitarySynthesisState,
+    constraint: QpuConstraint,
+) -> PyResult<bool> {
+    let Some(result) =
+        synthesize_2q_matrix(unitary, qargs_phys, state, constraint, fidelity_2q_sequence)?
+    else {
+        return Ok(false);
+    };
     // ... now apply the best sequence.
-    let (dir, sequence) = best_pair;
-    let order = dir.as_indices();
+    let order = result.dir.as_indices();
     let out_qargs = [qargs_virt[order[0] as usize], qargs_virt[order[1] as usize]];
     let qubit_keys = [
         out.insert_qargs(&[out_qargs[0]]),
@@ -598,8 +661,8 @@ fn synthesize_2q_matrix_onto(
         out.insert_qargs(&[out_qargs[0], out_qargs[1]]),
         out.insert_qargs(&[out_qargs[1], out_qargs[0]]),
     ];
-    out.add_global_phase(&Param::Float(sequence.global_phase()))?;
-    for (gate, params, qubits) in sequence.gates() {
+    out.add_global_phase(&Param::Float(result.sequence.global_phase()))?;
+    for (gate, params, qubits) in result.sequence.gates() {
         let qubits = match qubits.as_slice() {
             [0] => qubit_keys[0],
             [1] => qubit_keys[1],
diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 79e09acb8369..37242055de44 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -128,6 +128,7 @@
 sys.modules["qiskit._accelerate.inverse_cancellation"] = _accelerate.inverse_cancellation
 sys.modules["qiskit._accelerate.check_map"] = _accelerate.check_map
 sys.modules["qiskit._accelerate.filter_op_nodes"] = _accelerate.filter_op_nodes
+sys.modules["qiskit._accelerate.two_qubit_peephole"] = _accelerate.two_qubit_peephole
 sys.modules["qiskit._accelerate.twirling"] = _accelerate.twirling
 sys.modules["qiskit._accelerate.high_level_synthesis"] = _accelerate.high_level_synthesis
 sys.modules["qiskit._accelerate.remove_identity_equiv"] = _accelerate.remove_identity_equiv
diff --git a/qiskit/transpiler/passes/__init__.py b/qiskit/transpiler/passes/__init__.py
index 6787f2c95158..fe08fe9ab521 100644
--- a/qiskit/transpiler/passes/__init__.py
+++ b/qiskit/transpiler/passes/__init__.py
@@ -99,6 +99,7 @@
    Split2QUnitaries
    SubstitutePi4Rotations
    TemplateOptimization
+   TwoQubitPeepholeOptimization
 
 Scheduling
 =============
@@ -252,6 +253,7 @@
 from .optimization import Split2QUnitaries
 from .optimization import SubstitutePi4Rotations
 from .optimization import TemplateOptimization
+from .optimization import TwoQubitPeepholeOptimization
 
 # circuit analysis
 from .analysis import CountOps
diff --git a/qiskit/transpiler/passes/optimization/__init__.py b/qiskit/transpiler/passes/optimization/__init__.py
index 7c08b5eb0c1b..0dc12dd27698 100644
--- a/qiskit/transpiler/passes/optimization/__init__.py
+++ b/qiskit/transpiler/passes/optimization/__init__.py
@@ -39,6 +39,7 @@
 from .remove_identity_equiv import RemoveIdentityEquivalent
 from .split_2q_unitaries import Split2QUnitaries
 from .collect_and_collapse import CollectAndCollapse
+from .two_qubit_peephole import TwoQubitPeepholeOptimization
 from .contract_idle_wires_in_control_flow import ContractIdleWiresInControlFlow
 from .optimize_clifford_t import OptimizeCliffordT
 from .litinski_transformation import LitinskiTransformation
diff --git a/qiskit/transpiler/passes/optimization/two_qubit_peephole.py b/qiskit/transpiler/passes/optimization/two_qubit_peephole.py
new file mode 100644
index 000000000000..9d5b381d0663
--- /dev/null
+++ b/qiskit/transpiler/passes/optimization/two_qubit_peephole.py
@@ -0,0 +1,124 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 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 https://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.
+
+"""A transpiler pass to optimize 2q blocks in a circuit."""
+
+from __future__ import annotations
+
+from qiskit.transpiler.basepasses import TransformationPass
+from qiskit.transpiler.target import Target
+from qiskit.dagcircuit.dagcircuit import DAGCircuit
+from qiskit._accelerate.two_qubit_peephole import two_qubit_unitary_peephole_optimize
+
+
+class TwoQubitPeepholeOptimization(TransformationPass):
+    """Unified two qubit unitary peephole optimization
+
+    This transpiler pass is designed to perform two qubit unitary peephole
+    optimization. This pass finds all the 2 qubit blocks in the circuit,
+    computes the unitary of that block, and then synthesizes that unitary.
+    If the synthesized two qubit unitary is "better" than the original
+    subcircuit that subcircuit is used to replace the original. The heuristic
+    used to determine if it's better first looks at the two qubit gate count
+    in the circuit, and prefers the synthesis with fewer two qubit gates, if
+    the two qubit gate counts are the same then it looks at the estimated
+    fidelity of the circuit and picks the subcircuit with higher estimated
+    fidelity, and finally if needed it picks the subcircuit with the fewest
+    total gates.
+
+    In case the target is overcomplete the pass will try all the
+    decomposers supported for all the gates supported on a given qubit.
+    The decomposition that has the best expected performance using the above
+    heuristic will be selected and used to replace the block.
+
+    This pass is designed to be run on a physical circuit and the details of
+    operations on a given qubit is assumed to be the hardware qubit from the
+    target. However, the output of the pass might not use hardware operations,
+    specifically single qubit gates might be emitted outside the target's supported
+    operations, typically only if a parameterized gate supported by the
+    :class:`.TwoQubitControlledUDecomposer` is used for synthesis. As such if running
+    this pass in a physical optimization stage (such as :ref:`transpiler-preset-stage-optimization`)
+    this should be paired with passes such as :class:`.BasisTranslator` and/or
+    :class:`.Optimize1qGatesDecomposition` to ensure that these errant single qubit
+    gates are replaced with hardware supported operations prior to exiting the stage.
+
+    This pass is multithreaded, and will perform the analysis in parallel
+    and use all the cores available on your local system. You can refer to
+    the `configuration guide <https://docs.quantum.ibm.com/guides/configure-qiskit-local>`__
+    for details on how to control the threading behavior for Qiskit more broadly
+    which will also control this pass
+
+    This pass is similar in functionality to running :class:`.ConsolidateBlocks`
+    and :class:`.UnitarySynthesis` sequentially in your pass manager. However,
+    this pass offers improved runtime performance by performing the synthesis in
+    parallel. It also has improved heuristics enabled by doing the optimization in
+    a single step which can result in better quality output, especially in cases
+    of overcomplete and/or heterogeneous targets. However, these heuristics and this
+    pass as a whole are only valid for physical circuits. Additionally, unlike
+    :class:`.UnitarySynthesis` pass this does not use :ref:`unitary-synth-plugin`.
+    This is a tradeoff for performance and it forgoes the pluggability exposed
+    via that interface. Internally it currently only uses the :class:`.TwoQubitBasisDecomposer`
+    and :class:`.TwoQubitControlledUDecomposer` for synthesizing the two qubit unitaries.
+    You should not use this pass if you need to use the pluggable interface and the ability
+    to use different synthesis algorithms, instead you should use a combination of
+    :class:`.ConsolidateBlocks` and :class:`.UnitarySynthesis` to leverage the plugin
+    mechanism in :class:`.UnitarySynthesis`.
+
+    .. plot::
+      :include-source:
+      :alt: Optimized circuit
+
+      from qiskit.circuit import QuantumCircuit
+      from qiskit.transpiler.passes import TwoQubitPeepholeOptimization
+      from qiskit.providers.fake_provider import GenericBackendV2
+
+      # Build an unoptimized 2 qubit circuit
+      unoptimized = QuantumCircuit(2)
+      for i in range(10):
+        if i % 2:
+          unoptimized.cx(0, 1)
+        else:
+          unoptimized.cx(1, 0)
+      # Generate a target with random error rates
+      backend = GenericBackendV2(2, ["u", "cx"], coupling_map=[[0, 1], [1, 0]])
+      # Instantiate pass
+      peephole_pass = TwoQubitPeepholeOptimization(backend.target)
+      # Run pass and visualize output
+      optimized = peephole_pass(unoptimized)
+      optimized.draw("mpl")
+    """
+
+    def __init__(
+        self,
+        target: Target,
+        approximation_degree: float | None = 1.0,
+    ):
+        """Initialize the pass
+
+        Args:
+            target: The target to run the pass for
+            approximation_degree: heuristic dial used for circuit approximation (1.0=no
+                approximation, 0.0=maximal approximation). Approximation can decrease the number
+                of gates used in the synthesized unitaries smaller at the cost of straying from the
+                original unitary. If ``None``, approximation is done based on gate fidelities
+                specified in the ``target``.
+        """
+
+        super().__init__()
+        self._target = target
+        self._approximation_degree = approximation_degree
+
+    def run(self, dag: DAGCircuit) -> DAGCircuit:
+        result = two_qubit_unitary_peephole_optimize(dag, self._target, self._approximation_degree)
+        if result is None:
+            return dag
+        return result
diff --git a/qiskit/transpiler/passes/synthesis/plugin.py b/qiskit/transpiler/passes/synthesis/plugin.py
index 11e6c180036a..1f41837e7a3b 100644
--- a/qiskit/transpiler/passes/synthesis/plugin.py
+++ b/qiskit/transpiler/passes/synthesis/plugin.py
@@ -254,6 +254,8 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **
 Using Plugins
 =============
 
+.. _unitary-synth-plugin:
+
 Unitary Synthesis Plugins
 -------------------------
 
diff --git a/qiskit/transpiler/preset_passmanagers/builtin_plugins.py b/qiskit/transpiler/preset_passmanagers/builtin_plugins.py
index b53a6ff39e22..fc60ef6ba0f8 100644
--- a/qiskit/transpiler/preset_passmanagers/builtin_plugins.py
+++ b/qiskit/transpiler/preset_passmanagers/builtin_plugins.py
@@ -34,6 +34,7 @@
 from qiskit.transpiler.passes import ElidePermutations
 from qiskit.transpiler.passes import RemoveDiagonalGatesBeforeMeasure
 from qiskit.transpiler.passes import CommutativeOptimization
+from qiskit.transpiler.passes import TwoQubitPeepholeOptimization
 from qiskit.transpiler.passes import BasisTranslator
 from qiskit.transpiler.passes import SynthesizeRZRotations
 from qiskit.transpiler.passes import OptimizeCliffordT
@@ -504,18 +505,9 @@ def pass_manager(self, pass_manager_config, optimization_level=None):
                 loop_check, continue_loop = _optimization_check_fixed_point()
             case 2:
                 pre_loop = [
-                    ConsolidateBlocks(
-                        basis_gates=pass_manager_config.basis_gates,
-                        target=pass_manager_config.target,
-                        approximation_degree=pass_manager_config.approximation_degree,
-                    ),
-                    UnitarySynthesis(
-                        pass_manager_config.basis_gates,
+                    TwoQubitPeepholeOptimization(
+                        pass_manager_config.target,
                         approximation_degree=pass_manager_config.approximation_degree,
-                        coupling_map=pass_manager_config.coupling_map,
-                        method=pass_manager_config.unitary_synthesis_method,
-                        plugin_config=pass_manager_config.unitary_synthesis_plugin_config,
-                        target=pass_manager_config.target,
                     ),
                 ]
                 loop = [
@@ -534,18 +526,9 @@ def pass_manager(self, pass_manager_config, optimization_level=None):
             case 3:
                 pre_loop = []
                 loop = [
-                    ConsolidateBlocks(
-                        basis_gates=pass_manager_config.basis_gates,
-                        target=pass_manager_config.target,
-                        approximation_degree=pass_manager_config.approximation_degree,
-                    ),
-                    UnitarySynthesis(
-                        pass_manager_config.basis_gates,
+                    TwoQubitPeepholeOptimization(
+                        pass_manager_config.target,
                         approximation_degree=pass_manager_config.approximation_degree,
-                        coupling_map=pass_manager_config.coupling_map,
-                        method=pass_manager_config.unitary_synthesis_method,
-                        plugin_config=pass_manager_config.unitary_synthesis_plugin_config,
-                        target=pass_manager_config.target,
                     ),
                     RemoveIdentityEquivalent(
                         approximation_degree=pass_manager_config.approximation_degree,
diff --git a/releasenotes/notes/two_qubit_peephole-de6d3438ed7df6a9.yaml b/releasenotes/notes/two_qubit_peephole-de6d3438ed7df6a9.yaml
new file mode 100644
index 000000000000..27891fbed377
--- /dev/null
+++ b/releasenotes/notes/two_qubit_peephole-de6d3438ed7df6a9.yaml
@@ -0,0 +1,38 @@
+---
+features_transpiler:
+  - |
+    Added a new transpiler pass :class:`.TwoQubitPeepholeOptimization` which
+    is intended to perform two qubit unitary peephole optimization. The pass
+    analyzes the circuit to find two qubit blocks in the circuit, it then
+    computes the unitary of that subcircuit, and will replace the original
+    block with the synthesized unitary which uses fewer operations. For example:
+
+    .. plot::
+
+        from qiskit.circuit import QuantumCircuit
+        from qiskit.transpiler.passes import TwoQubitPeepholeOptimization
+        from qiskit.providers.fake_provider import GenericBackendV2
+
+        # Build an unoptimized 2 qubit circuit
+        unoptimized = QuantumCircuit(2)
+        for i in range(10):
+          if i % 2:
+            unoptimized.cx(0, 1)
+          else:
+            unoptimized.cx(1, 0)
+
+        # Generate a target with random error rates
+        target = GenericBackendV2(2, ["u", "cx"], coupling_map=[[0, 1]]).target
+        # Instantiate pass
+        peephole_pass = TwoQubitPeepholeOptimization(target)
+        # Run pass and visualize output
+        optimized = peephole_pass(unoptimized)
+        optimized.draw("mpl")
+
+    This functionality could previously have been perfomed by running
+    :class:`.ConsolidateBlocks`, and:class:`.UnitarySynthesis` sequentially
+    in your pass manager. However, this new pass offers improved runtime
+    performance by performing the synthesis in parallel. It also has improved
+    heuristics enabled by doing the optimization in a single step which can
+    result in better quality output, especially in cases of overcomplete and/or
+    hetergeneous targets.
diff --git a/test/python/circuit/test_scheduled_circuit.py b/test/python/circuit/test_scheduled_circuit.py
index 3972f0a681f0..3a6d1a794002 100644
--- a/test/python/circuit/test_scheduled_circuit.py
+++ b/test/python/circuit/test_scheduled_circuit.py
@@ -157,8 +157,13 @@ def test_transpile_delay_circuit_with_backend(self):
             qc, backend=self.backend_with_dt, scheduling_method="alap", layout_method="trivial"
         )
         target_durations = self.backend_with_dt.target.durations()
+        # The longest path output is 3 sx gates 1 cx(1,) and two delays for
+        # alignment
+        cx_duration = target_durations.get("cx", (1, 0))
+        sx_duration = target_durations.get("sx", (0,))
+        expected_duration = sx_duration * 3 + cx_duration + 270 + 15
         with self.assertWarns(DeprecationWarning):
-            self.assertEqual(scheduled.duration, target_durations.get("cx", (0, 1)) + 450)
+            self.assertEqual(scheduled.duration, expected_duration)
 
     def test_transpile_circuit_with_custom_instruction(self):
         """See: https://github.com/Qiskit/qiskit-terra/issues/5154"""
@@ -239,8 +244,11 @@ def test_unit_seconds_when_using_backend_durations(self):
             qc, backend=self.backend_with_dt, scheduling_method="alap", layout_method="trivial"
         )
         target_durations = self.backend_with_dt.target.durations()
+        cx_duration = target_durations.get("cx", (1, 0))
+        sx_duration = target_durations.get("sx", (0,))
+        expected_duration = cx_duration + 3 * sx_duration + 320 + 15
         with self.assertWarns(DeprecationWarning):
-            self.assertEqual(scheduled.duration, target_durations.get("cx", (0, 1)) + 500)
+            self.assertEqual(scheduled.duration, expected_duration)
 
     def test_per_qubit_durations(self):
         """Test target with custom instruction_durations"""
diff --git a/test/python/transpiler/test_preset_passmanagers.py b/test/python/transpiler/test_preset_passmanagers.py
index fcf0df5dd4c2..789e1200f51b 100644
--- a/test/python/transpiler/test_preset_passmanagers.py
+++ b/test/python/transpiler/test_preset_passmanagers.py
@@ -48,10 +48,10 @@
 from qiskit.providers.fake_provider import GenericBackendV2
 from qiskit.converters import circuit_to_dag
 from qiskit.circuit.library import GraphStateGate, UnitaryGate
-from qiskit.quantum_info import random_unitary
+from qiskit.quantum_info import random_unitary, Operator
 from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
 from qiskit.transpiler.preset_passmanagers import level0, level1, level2, level3
-from qiskit.transpiler.passes import ConsolidateBlocks, GatesInBasis
+from qiskit.transpiler.passes import ConsolidateBlocks, GatesInBasis, TwoQubitPeepholeOptimization
 from qiskit.transpiler.passes.scheduling.alignments.check_durations import InstructionDurationCheck
 from qiskit.transpiler.preset_passmanagers.builtin_plugins import OptimizationPassManager
 from qiskit.transpiler.timing_constraints import TimingConstraints
@@ -291,15 +291,15 @@ def test_unroll_only_if_not_gates_in_basis(self):
         )
         qv_circuit = quantum_volume(3)
         gates_in_basis_true_count = 0
-        consolidate_blocks_count = 0
+        peephole_count = 0
 
         def counting_callback_func(pass_, dag, time, property_set, count):
             nonlocal gates_in_basis_true_count
-            nonlocal consolidate_blocks_count
+            nonlocal peephole_count
             if isinstance(pass_, GatesInBasis) and property_set["all_gates_in_basis"]:
                 gates_in_basis_true_count += 1
-            if isinstance(pass_, ConsolidateBlocks):
-                consolidate_blocks_count += 1
+            if isinstance(pass_, TwoQubitPeepholeOptimization):
+                peephole_count += 1
 
         transpile(
             qv_circuit,
@@ -308,7 +308,7 @@ def counting_callback_func(pass_, dag, time, property_set, count):
             callback=counting_callback_func,
             translation_method="synthesis",
         )
-        self.assertEqual(gates_in_basis_true_count + 2, consolidate_blocks_count)
+        self.assertEqual(gates_in_basis_true_count, peephole_count)
 
     @data(0, 1, 2, 3)
     def test_layout_registers_preserved(self, optimization_level):
@@ -1059,7 +1059,7 @@ def test_layout_tokyo_fully_connected_cx(self, level):
             [0, 1, 2, 3, 4],
             [5, 6, 10, 0, 11],
             [5, 6, 10, 0, 11],
-            [5, 11, 6, 0, 10],
+            [6, 7, 1, 5, 2],
         ]
         backend = GenericBackendV2(num_qubits=20, coupling_map=TOKYO_CMAP, seed=42)
         result = transpile(qc, backend, optimization_level=level, seed_transpiler=42)
diff --git a/test/python/transpiler/test_two_qubit_peephole.py b/test/python/transpiler/test_two_qubit_peephole.py
new file mode 100644
index 000000000000..b66d1dccd9e3
--- /dev/null
+++ b/test/python/transpiler/test_two_qubit_peephole.py
@@ -0,0 +1,781 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2025
+#
+# 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 https://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.
+
+# pylint: disable=missing-function-docstring
+
+"""
+Tests for the TwoQubitPeepholeOptimization transpiler pass.
+"""
+
+import math
+import numpy as np
+import ddt
+
+from qiskit import generate_preset_pass_manager, transpile
+from qiskit.providers.fake_provider import GenericBackendV2
+from qiskit.circuit import QuantumCircuit, QuantumRegister
+from qiskit.circuit.random import random_circuit
+from qiskit.circuit.parameterexpression import ParameterValueType
+from qiskit.converters import circuit_to_dag, dag_to_circuit
+from qiskit.quantum_info.operators import Operator
+from qiskit.quantum_info.random import random_unitary
+from qiskit.transpiler import PassManager, CouplingMap, Target, InstructionProperties
+from qiskit.transpiler.passes import (
+    TwoQubitPeepholeOptimization,
+    TrivialLayout,
+    ConsolidateBlocks,
+    UnitarySynthesis,
+)
+from qiskit.circuit.library import (
+    IGate,
+    CXGate,
+    RZGate,
+    RXGate,
+    SXGate,
+    XGate,
+    iSwapGate,
+    ECRGate,
+    UGate,
+    ZGate,
+    RYYGate,
+    CZGate,
+    RXXGate,
+    RZZGate,
+    RZXGate,
+    CPhaseGate,
+    CRZGate,
+    CRXGate,
+    CRYGate,
+    CUGate,
+    XXPlusYYGate,
+    XXMinusYYGate,
+)
+from qiskit.circuit import Measure
+from qiskit.circuit.controlflow import IfElseOp
+from qiskit.circuit import Parameter, Gate
+
+from test import combine  # pylint: disable=wrong-import-order
+from test import QiskitTestCase  # pylint: disable=wrong-import-order
+from test.python.providers.fake_mumbai_v2 import (  # pylint: disable=wrong-import-order
+    FakeMumbaiFractionalCX,
+)
+from ..legacy_cmaps import YORKTOWN_CMAP
+
+
+class FakeBackend2QV2(GenericBackendV2):
+    """A 2-qubit fake backend"""
+
+    def __init__(self):
+        super().__init__(num_qubits=2, basis_gates=["rx", "u"], seed=42)
+        cx_props = {
+            (0, 1): InstructionProperties(duration=5.23e-7, error=0.00098115),
+        }
+        self._target.add_instruction(CXGate(), cx_props)
+        ecr_props = {
+            (1, 0): InstructionProperties(duration=4.52e-9, error=0.0000132115),
+        }
+        self._target.add_instruction(ECRGate(), ecr_props)
+
+
+class FakeBackend5QV2(GenericBackendV2):
+    """A 5-qubit fake backend"""
+
+    def __init__(self, bidirectional=True):
+        super().__init__(num_qubits=5, basis_gates=["u"], seed=42)
+        cx_props = {
+            (0, 1): InstructionProperties(duration=5.23e-7, error=0.00098115),
+            (3, 4): InstructionProperties(duration=5.23e-7, error=0.00098115),
+        }
+        if bidirectional:
+            cx_props[(1, 0)] = InstructionProperties(duration=6.23e-7, error=0.00099115)
+            cx_props[(4, 3)] = InstructionProperties(duration=7.23e-7, error=0.00099115)
+        self._target.add_instruction(CXGate(), cx_props)
+        ecr_props = {
+            (1, 2): InstructionProperties(duration=4.52e-9, error=0.0000132115),
+            (2, 3): InstructionProperties(duration=4.52e-9, error=0.0000132115),
+        }
+        if bidirectional:
+            ecr_props[(2, 1)] = InstructionProperties(duration=5.52e-9, error=0.0000232115)
+            ecr_props[(3, 2)] = InstructionProperties(duration=5.52e-9, error=0.0000232115)
+
+
+@ddt.ddt
+class TestTwoQubitPeepholeOptimization(QiskitTestCase):
+    """Test TwoQubitPeepholeOptimization."""
+
+    @combine(
+        bidirectional=[True, False],
+        dsc=(
+            "test natural_direction works with transpile using a"
+            "target with multiple 2q gates with bidirectional={bidirectional}"
+        ),
+        name="bidirectional_{bidirectional}",
+    )
+    def test_coupling_map_transpile_with_backendv2(self, bidirectional):
+        backend = FakeBackend5QV2(bidirectional)
+        qr = QuantumRegister(2)
+        circ = QuantumCircuit(qr)
+        circ.append(random_unitary(4, seed=1), [0, 1])
+        circ_01 = TwoQubitPeepholeOptimization(backend.target)(circ)
+        circ_01_index = {qubit: index for index, qubit in enumerate(circ_01.qubits)}
+        self.assertGreaterEqual(len(circ_01.get_instructions("cx")), 1)
+        for instr in circ_01.get_instructions("cx"):
+            self.assertEqual(
+                (0, 1), (circ_01_index[instr.qubits[0]], circ_01_index[instr.qubits[1]])
+            )
+
+    @combine(
+        bidirectional=[True, False],
+        dsc=(
+            "Test direction with transpile using a "
+            "target with multiple 2q gates with bidirectional={bidirectional}"
+            "direction [0, 1] is lower error and should be picked."
+        ),
+        name="bidirectional_{bidirectional}",
+    )
+    def test_coupling_unequal_duration_with_backendv2(self, bidirectional):
+        qr = QuantumRegister(2)
+        circ = QuantumCircuit(qr)
+        circ.append(random_unitary(4, seed=1), [1, 0])
+        backend = FakeBackend5QV2(bidirectional)
+        tqc = TwoQubitPeepholeOptimization(backend.target)(circ)
+        tqc_index = {qubit: index for index, qubit in enumerate(tqc.qubits)}
+        self.assertGreaterEqual(len(tqc.get_instructions("cx")), 1)
+        for instr in tqc.get_instructions("cx"):
+            self.assertEqual((0, 1), (tqc_index[instr.qubits[0]], tqc_index[instr.qubits[1]]))
+
+    def test_non_overlapping_kak_gates_with_backendv2(self):
+        qr = QuantumRegister(2)
+        circ = QuantumCircuit(qr)
+        circ.append(random_unitary(4, seed=1), [1, 0])
+        backend = FakeBackend2QV2()
+        tqc = TwoQubitPeepholeOptimization(backend.target)(circ)
+
+        tqc_index = {qubit: index for index, qubit in enumerate(tqc.qubits)}
+        self.assertGreaterEqual(len(tqc.get_instructions("ecr")), 1)
+        for instr in tqc.get_instructions("ecr"):
+            self.assertEqual((1, 0), (tqc_index[instr.qubits[0]], tqc_index[instr.qubits[1]]))
+
+    def test_fractional_cx_with_backendv2(self):
+        """Test fractional CX gets used if present in target."""
+        qr = QuantumRegister(2)
+        circ = QuantumCircuit(qr)
+        circ.append(random_unitary(4, seed=1), [0, 1])
+        backend = FakeMumbaiFractionalCX()
+        synth_pass = TwoQubitPeepholeOptimization(target=backend.target)
+        tqc = synth_pass(circ)
+        np.testing.assert_allclose(Operator.from_circuit(tqc), Operator.from_circuit(circ))
+
+    def test_reverse_direction(self):
+        target = Target(2)
+        target.add_instruction(CXGate(), {(0, 1): InstructionProperties(error=1.2e-6)})
+        target.add_instruction(ECRGate(), {(0, 1): InstructionProperties(error=1.2e-7)})
+        target.add_instruction(
+            UGate(Parameter("theta"), Parameter("phi"), Parameter("lam")), {(0,): None, (1,): None}
+        )
+        qr = QuantumRegister(2)
+        circ = QuantumCircuit(qr)
+        circ.append(random_unitary(4, seed=1), [1, 0])
+        tqc = TwoQubitPeepholeOptimization(target)(circ)
+        tqc_index = {qubit: index for index, qubit in enumerate(tqc.qubits)}
+        self.assertGreaterEqual(len(tqc.get_instructions("ecr")), 1)
+        for instr in tqc.get_instructions("ecr"):
+            self.assertEqual((0, 1), (tqc_index[instr.qubits[0]], tqc_index[instr.qubits[1]]))
+
+    def test_controlled_basis(self):
+        target = Target(2)
+        target.add_instruction(RYYGate(np.pi / 8), {(0, 1): InstructionProperties(error=1.2e-6)})
+        target.add_instruction(
+            UGate(Parameter("theta"), Parameter("phi"), Parameter("lam")), {(0,): None, (1,): None}
+        )
+        qr = QuantumRegister(2)
+        circ = QuantumCircuit(qr)
+        circ.append(random_unitary(4, seed=1), [1, 0])
+        tqc = TwoQubitPeepholeOptimization(target)(circ)
+        # Until XX decomposer is ported we can't synthesize using a RYY(pi/8) as the only 2q
+        # gate in the target
+        self.assertGreaterEqual(len(tqc.get_instructions("unitary")), 1)
+
+    def test_mapping_control_flow(self):
+        """Test that inner dags use proper qubit mapping."""
+        qr = QuantumRegister(3, "q")
+        qc = QuantumCircuit(qr)
+
+        # Create target that supports CX only between 0 and 2.
+        fake_target = Target()
+        fake_target.add_instruction(CXGate(), {(0, 2): None})
+        fake_target.add_instruction(
+            UGate(Parameter("t"), Parameter("p"), Parameter("l")),
+            {
+                (0,): None,
+                (1,): None,
+                (2,): None,
+            },
+        )
+
+        qc_uni1 = QuantumCircuit(2)
+        qc_uni1.swap(0, 1)
+        qc_uni1_mat = Operator(qc_uni1)
+
+        loop_body = QuantumCircuit(2)
+        loop_body.unitary(qc_uni1_mat, [0, 1])
+
+        # Loop body uses qubits 0 and 2, mapped to 0 and 1 in the block.
+        # If synthesis doesn't handle recursive mapping, it'll incorrectly
+        # look for a CX on (0, 1) instead of on (0, 2).
+        qc.for_loop((0,), None, loop_body, [0, 2], [])
+
+        result = TwoQubitPeepholeOptimization(fake_target)(qc)
+        self.assertIsInstance(result, QuantumCircuit)
+
+    def test_single_qubit_with_target(self):
+        """Test input circuit with only 1q works with target."""
+        qc = QuantumCircuit(1)
+        qc.append(ZGate(), [qc.qubits[0]])
+        dag = circuit_to_dag(qc)
+        backend = GenericBackendV2(num_qubits=5, seed=42)
+        unitary_synth_pass = TwoQubitPeepholeOptimization(target=backend.target)
+        result_dag = unitary_synth_pass.run(dag)
+        result_qc = dag_to_circuit(result_dag)
+        self.assertEqual(qc, result_qc)
+
+    def test_two_qubit_identity_with_target(self):
+        """Test input single qubit identity works with target."""
+        qc = QuantumCircuit(2)
+        qc.unitary(np.eye(4, dtype=complex), [0, 1])
+        dag = circuit_to_dag(qc)
+        backend = GenericBackendV2(num_qubits=5)
+        unitary_synth_pass = TwoQubitPeepholeOptimization(target=backend.target)
+        result_dag = unitary_synth_pass.run(dag)
+        result_qc = dag_to_circuit(result_dag)
+        self.assertEqual(result_qc, QuantumCircuit(2))
+
+    def test_unitary_synthesis_with_ideal_and_variable_width_ops(self):
+        """Test unitary synthesis works with a target that contains ideal and variadic ops."""
+        qc = QuantumCircuit(2)
+        qc.unitary(np.eye(4), [0, 1])
+        dag = circuit_to_dag(qc)
+        target = GenericBackendV2(num_qubits=5).target
+        target.add_instruction(IfElseOp, name="if_else")
+        target.add_instruction(ZGate())
+        target.add_instruction(ECRGate())
+        unitary_synth_pass = TwoQubitPeepholeOptimization(target=target)
+        result_dag = unitary_synth_pass.run(dag)
+        result_qc = dag_to_circuit(result_dag)
+        self.assertEqual(result_qc, QuantumCircuit(2))
+
+    def test_unitary_synthesis_custom_gate_target(self):
+        qc = QuantumCircuit(2)
+        qc.unitary(np.eye(4), [0, 1])
+        dag = circuit_to_dag(qc)
+
+        class CustomGate(Gate):
+            """Custom Opaque Gate"""
+
+            def __init__(self):
+                super().__init__("custom", 2, [])
+
+        target = Target(num_qubits=2)
+        target.add_instruction(
+            UGate(Parameter("t"), Parameter("p"), Parameter("l")), {(0,): None, (1,): None}
+        )
+        target.add_instruction(CustomGate(), {(0, 1): None, (1, 0): None})
+        unitary_synth_pass = TwoQubitPeepholeOptimization(target=target)
+        result_dag = unitary_synth_pass.run(dag)
+        result_qc = dag_to_circuit(result_dag)
+        self.assertEqual(result_qc, qc)
+
+    def test_iswap_no_cx_synthesis_succeeds(self):
+        """Test basis set with iswap but no cx can synthesize a circuit"""
+        target = Target()
+        theta = Parameter("theta")
+
+        i_props = {
+            (0,): InstructionProperties(duration=35.5e-9, error=0.000413),
+            (1,): InstructionProperties(duration=35.5e-9, error=0.000502),
+        }
+        target.add_instruction(IGate(), i_props)
+        rz_props = {
+            (0,): InstructionProperties(duration=0, error=0),
+            (1,): InstructionProperties(duration=0, error=0),
+        }
+        target.add_instruction(RZGate(theta), rz_props)
+        sx_props = {
+            (0,): InstructionProperties(duration=35.5e-9, error=0.000413),
+            (1,): InstructionProperties(duration=35.5e-9, error=0.000502),
+        }
+        target.add_instruction(SXGate(), sx_props)
+        x_props = {
+            (0,): InstructionProperties(duration=35.5e-9, error=0.000413),
+            (1,): InstructionProperties(duration=35.5e-9, error=0.000502),
+        }
+        target.add_instruction(XGate(), x_props)
+        iswap_props = {
+            (0, 1): InstructionProperties(duration=519.11e-9, error=0.01201),
+            (1, 0): InstructionProperties(duration=554.66e-9, error=0.01201),
+        }
+        target.add_instruction(iSwapGate(), iswap_props)
+        measure_props = {
+            (0,): InstructionProperties(duration=5.813e-6, error=0.0751),
+            (1,): InstructionProperties(duration=5.813e-6, error=0.0225),
+        }
+        target.add_instruction(Measure(), measure_props)
+
+        qc = QuantumCircuit(2)
+        cxmat = Operator(CXGate()).to_matrix()
+        qc.unitary(cxmat, [0, 1])
+        unitary_synth_pass = TwoQubitPeepholeOptimization(target=target)
+        dag = circuit_to_dag(qc)
+        result_dag = unitary_synth_pass.run(dag)
+        result_qc = dag_to_circuit(result_dag)
+        self.assertTrue(np.allclose(Operator(result_qc.to_gate()).to_matrix(), cxmat))
+
+    def test_rxx_gate_in_target(self):
+        """Test synthesis with custom parameterized gate in target."""
+
+        theta = Parameter("θ")
+        lam = Parameter("λ")
+        phi = Parameter("ϕ")
+
+        target = Target(num_qubits=2)
+        target.add_instruction(RZGate(lam))
+        target.add_instruction(RXGate(phi))
+        target.add_instruction(RXXGate(theta))
+
+        qc = QuantumCircuit(2)
+        qc.unitary(random_unitary(4, seed=1234), [0, 1])
+        qc_transpiled = TwoQubitPeepholeOptimization(target=target)(qc)
+        opcount = qc_transpiled.count_ops()
+        self.assertTrue(set(opcount).issubset({"rz", "rx", "rxx"}))
+
+        self.assertTrue(np.allclose(Operator(qc_transpiled), Operator(qc)))
+
+    def test_custom_parameterized_gate_in_target_skips(self):
+        """Test that synthesis is skipped with custom parameterized
+        gate in target that is not RXX equivalent."""
+
+        class CustomXYGate(Gate):
+            """Custom Gate subclass that's not a standard gate and not RXX equivalent"""
+
+            _standard_gate = None
+
+            def __init__(self, theta: ParameterValueType, label=None):
+                """Create new custom rotstion XY gate."""
+                super().__init__("MyCustomXYGate", 2, [theta])
+
+            def __array__(self, dtype=None):
+                """Return a Numpy.array for the custom gate."""
+                theta = self.params[0]
+                cos = math.cos(theta)
+                isin = 1j * math.sin(theta)
+                return np.array(
+                    [[1, 0, 0, 0], [0, cos, -isin, 0], [0, -isin, cos, 0], [0, 0, 0, 1]],
+                    dtype=dtype,
+                )
+
+            def inverse(self, annotated: bool = False):
+                return CustomXYGate(-self.params[0])
+
+        theta = Parameter("θ")
+        lam = Parameter("λ")
+        phi = Parameter("ϕ")
+
+        target = Target(num_qubits=2)
+        target.add_instruction(RZGate(lam))
+        target.add_instruction(RXGate(phi))
+        target.add_instruction(CustomXYGate(theta))
+
+        qc = QuantumCircuit(2)
+        qc.unitary(random_unitary(4, seed=1234), [0, 1])
+        qc_transpiled = TwoQubitPeepholeOptimization(target=target)(qc)
+        opcount = qc_transpiled.count_ops()
+        self.assertTrue(set(opcount).issubset({"unitary"}))
+        self.assertTrue(np.allclose(Operator(qc_transpiled), Operator(qc)))
+
+    def test_determinism(self):
+        """Test that the decomposition is deterministic."""
+        gate_counts = {"rx": 6, "rz": 12, "iswap": 2}
+        basis_gates = ["rx", "rz", "iswap"]
+        target = Target.from_configuration(basis_gates=basis_gates)
+        pm = generate_preset_pass_manager(target=target, optimization_level=2, seed_transpiler=42)
+
+        qc = QuantumCircuit(2)
+        qc.h(0)
+        qc.cx(0, 1)
+
+        for _ in range(10):
+            out = pm.run(qc)
+            self.assertTrue(Operator(out).equiv(qc))
+            self.assertTrue(set(out.count_ops()).issubset(basis_gates))
+            for basis_gate in basis_gates:
+                self.assertLessEqual(out.count_ops()[basis_gate], gate_counts[basis_gate])
+
+    @combine(gate=["unitary", "swap"])
+    def test_two_qubit_synthesis_to_directional_cx_target(self, gate):
+        """Verify two qubit unitaries are synthesized to match basis gates."""
+        # TODO: should make check more explicit e.g. explicitly set gate
+        # direction in test instead of using specific fake backend
+        backend = GenericBackendV2(
+            num_qubits=5,
+            basis_gates=["id", "rz", "sx", "x", "cx", "reset"],
+            coupling_map=YORKTOWN_CMAP,
+            seed=1,
+        )
+        coupling_map = CouplingMap(backend.coupling_map)
+        triv_layout_pass = TrivialLayout(coupling_map)
+
+        qr = QuantumRegister(2)
+        qc = QuantumCircuit(qr)
+        if gate == "unitary":
+            qc.unitary(random_unitary(4, seed=12), [0, 1])
+        elif gate == "swap":
+            qc.swap(qr[0], qr[1])
+
+        unisynth_pass = TwoQubitPeepholeOptimization(
+            target=backend.target,
+        )
+        pm = PassManager([triv_layout_pass, unisynth_pass])
+        qc_out = pm.run(qc)
+        self.assertEqual(Operator(qc), Operator(qc_out))
+
+    def test_swap_on_cz_target(self):
+        """Test that a cz target synthesizes simple circuits correctly."""
+        target = Target(num_qubits=2)
+        target.add_instruction(
+            CZGate(),
+            {
+                (0, 1): InstructionProperties(error=3.058e-3, duration=6.8e-8),
+                (1, 0): InstructionProperties(error=3.058e-3, duration=6.8e-8),
+            },
+        )
+        target.add_instruction(
+            RZGate(Parameter("theta")),
+            {
+                (0,): InstructionProperties(error=0, duration=0),
+                (1,): InstructionProperties(error=0, duration=0),
+            },
+        )
+        target.add_instruction(
+            SXGate(),
+            {
+                (0,): InstructionProperties(error=2.421e-4, duration=2.4e-8),
+                (1,): InstructionProperties(error=2.229e-4, duration=2.41e-8),
+            },
+        )
+        target.add_instruction(
+            XGate(),
+            {
+                (0,): InstructionProperties(error=2.41e-4, duration=3.4e-8),
+                (1,): InstructionProperties(error=2.29e-4, duration=3.41e-8),
+            },
+        )
+        peephole = TwoQubitPeepholeOptimization(target)
+        qc = QuantumCircuit(2)
+        qc.swap(0, 1)
+        qc = transpile(qc, target=target, seed_transpiler=1234, optimization_level=0)
+        res = peephole(qc)
+        np.testing.assert_allclose(Operator(qc), Operator(res), atol=1e-12, rtol=0)
+        # Check run of swaps
+        qc_duplicated = QuantumCircuit(2)
+        for _ in range(100):
+            qc_duplicated.swap(0, 1)
+        qc_duplicated = transpile(
+            qc_duplicated, target=target, seed_transpiler=1234, optimization_level=0
+        )
+        res = peephole(qc_duplicated)
+        self.assertEqual(res, QuantumCircuit(2))
+
+        qc_duplicated = QuantumCircuit(2)
+        for _ in range(101):
+            qc_duplicated.swap(0, 1)
+        qc_duplicated = transpile(
+            qc_duplicated, target=target, seed_transpiler=1234, optimization_level=0
+        )
+        res = peephole(qc_duplicated)
+        self.assertEqual(Operator(res), Operator(qc))
+
+    def test_pass_respects_directionality(self):
+        """Test that a cz target synthesizes simple circuits correctly."""
+        target = Target(num_qubits=2)
+        target.add_instruction(
+            CXGate(),
+            {
+                (1, 0): InstructionProperties(error=3.058e-3, duration=6.8e-8),
+                (0, 1): InstructionProperties(error=0.99999),
+            },
+        )
+        target.add_instruction(
+            RZGate(Parameter("theta")),
+            {
+                (0,): InstructionProperties(error=0, duration=0),
+                (1,): InstructionProperties(error=0, duration=0),
+            },
+        )
+        target.add_instruction(
+            SXGate(),
+            {
+                (0,): InstructionProperties(error=2.421e-4, duration=2.4e-8),
+                (1,): InstructionProperties(error=2.229e-4, duration=2.41e-8),
+            },
+        )
+        target.add_instruction(
+            XGate(),
+            {
+                (0,): InstructionProperties(error=2.41e-4, duration=3.4e-8),
+                (1,): InstructionProperties(error=2.29e-4, duration=3.41e-8),
+            },
+        )
+        peephole = TwoQubitPeepholeOptimization(target)
+        qc = QuantumCircuit(2)
+        qc.swap(0, 1)
+        res = peephole(qc)
+        self.assertTrue(self.all_inst_in_target(res, target))
+        self.assertEqual(Operator(res), Operator(qc))
+        for inst in res.data:
+            if len(inst.qubits) == 2:
+                self.assertEqual(inst.qubits, (res.qubits[1], res.qubits[0]))
+        # Check run of swaps
+        qc_duplicated = QuantumCircuit(2)
+        for _ in range(100):
+            qc_duplicated.swap(0, 1)
+        qc_duplicated = transpile(
+            qc_duplicated, target=target, seed_transpiler=1234, optimization_level=0
+        )
+        res = peephole(qc_duplicated)
+        self.assertEqual(Operator(res), Operator(qc_duplicated))
+        self.assertEqual(Operator(res), Operator(np.eye(4, dtype=complex)))
+        self.assertTrue(self.all_inst_in_target(res, target))
+
+        qc_duplicated = QuantumCircuit(2)
+        for _ in range(101):
+            qc_duplicated.swap(0, 1)
+        qc_duplicated = transpile(
+            qc_duplicated, target=target, seed_transpiler=1234, optimization_level=0
+        )
+        res = peephole(qc_duplicated)
+        self.assertTrue(self.all_inst_in_target(res, target))
+        self.assertEqual(Operator(res), Operator(qc_duplicated))
+        for inst in res.data:
+            if len(inst.qubits) == 2:
+                self.assertEqual(inst.qubits, (res.qubits[1], res.qubits[0]))
+
+    def all_inst_in_target(self, circuit: QuantumCircuit, target: Target, allow_inverse=False):
+        for inst in circuit.data:
+            if not target.instruction_supported(
+                inst.name, tuple(circuit.find_bit(x).index for x in inst.qubits)
+            ):
+                if allow_inverse:
+                    if not target.instruction_supported(
+                        inst.operation.inverse().name,
+                        tuple(circuit.find_bit(x).index for x in inst.qubits),
+                    ):
+                        raise self.fail(
+                            f"{inst.name} {tuple(circuit.find_bit(x).index for x in inst.qubits)} not supported"
+                        )
+
+                else:
+                    raise self.fail(
+                        f"{inst.name} {tuple(circuit.find_bit(x).index for x in inst.qubits)} not supported"
+                    )
+        return True
+
+    @combine(
+        gate=[
+            RXXGate(0.1),
+            RYYGate(0.1),
+            RZZGate(0.1),
+            RZXGate(0.1),
+            CPhaseGate(0.1),
+            CRZGate(0.1),
+            CRXGate(0.1),
+            CRYGate(0.1),
+            CUGate(0.1, 0.2, 0.3, 0.4),
+            XXPlusYYGate(0.1),
+            XXMinusYYGate(0.1),
+        ],
+        target_gate=[
+            CXGate(),
+            CZGate(),
+            ECRGate(),
+        ],
+        add_noise=[True, False],
+        name="{gate}_{target_gate}_noise={add_noise}",
+    )
+    def test_two_qubit_parametrized_gates_basis_decomp_target(self, gate, target_gate, add_noise):
+        """Test the synthesis of a circuit containing a 2-qubit parametrized gate
+        on a target with a CX gate"""
+        theta = Parameter("θ")
+        target = Target(num_qubits=2)
+        if add_noise:
+            target.add_instruction(
+                target_gate, {(i, i + 1): InstructionProperties(error=0.001) for i in [0]}
+            )
+            target.add_instruction(RZGate(theta))
+            target.add_instruction(
+                SXGate(), {(i,): InstructionProperties(error=0.0001) for i in [0, 1]}
+            )
+        else:
+            target.add_instruction(target_gate)
+            target.add_instruction(RZGate(theta))
+            target.add_instruction(SXGate())
+
+        qc = QuantumCircuit(2)
+        qc.append(gate, [0, 1])
+
+        peephole = TwoQubitPeepholeOptimization(target)
+        transpiled_circuit = peephole(qc)
+
+        legacy_path = PassManager(
+            [
+                ConsolidateBlocks(target=target),
+                UnitarySynthesis(
+                    target=target,
+                ),
+            ]
+        )
+        legacy = legacy_path.run(qc)
+        self.all_inst_in_target(transpiled_circuit, target)
+        np.testing.assert_allclose(Operator(transpiled_circuit), Operator(qc), atol=1e-12, rtol=0)
+        self.assertDictEqual(
+            dict(sorted(transpiled_circuit.count_ops().items())),
+            dict(sorted(legacy.count_ops().items())),
+        )
+
+    @combine(
+        gate=[
+            RXXGate(0.1),
+            RYYGate(0.1),
+            RZZGate(0.1),
+            RZXGate(0.1),
+            CPhaseGate(0.1),
+            CRZGate(0.1),
+            CRXGate(0.1),
+            CRYGate(0.1),
+            CUGate(0.1, 0.2, 0.3, 0.4),
+        ],
+        target_gate_cls=[
+            RZZGate,
+            RXXGate,
+            RZXGate,
+            RYYGate,
+        ],
+        add_noise=[True, False],
+        name="{gate}_{target_gate_cls}_noise={add_noise}",
+    )
+    def test_two_qubit_parametrized_gates_controlled_u_target(
+        self, gate, target_gate_cls, add_noise
+    ):
+        """Test the synthesis of a circuit containing a 2-qubit parametrized gate
+        on a target with a RZZ gate"""
+        theta = Parameter("θ")
+        lam = Parameter("λ")
+        phi = Parameter("ϕ")
+        target = Target(num_qubits=2)
+        target_gate = target_gate_cls(phi)
+        if add_noise:
+            target.add_instruction(
+                RXGate(lam), {(i,): InstructionProperties(error=0.0001) for i in [0, 1]}
+            )
+            target.add_instruction(RZGate(theta))
+            target.add_instruction(
+                target_gate, {(i, i + 1): InstructionProperties(error=0.001) for i in [0]}
+            )
+        else:
+            target.add_instruction(RXGate(lam))
+            target.add_instruction(RZGate(theta))
+            target.add_instruction(target_gate)
+        qc = QuantumCircuit(2)
+        qc.append(gate, [0, 1])
+        peephole = TwoQubitPeepholeOptimization(target)
+        transpiled_circuit = peephole(qc)
+        legacy_path = PassManager(
+            [
+                ConsolidateBlocks(target=target),
+                UnitarySynthesis(
+                    target=target,
+                ),
+            ]
+        )
+
+        legacy = legacy_path.run(qc)
+        self.all_inst_in_target(transpiled_circuit, target)
+        self.assertEqual(Operator(transpiled_circuit), Operator(qc))
+        self.assertDictEqual(
+            dict(sorted(transpiled_circuit.count_ops().items())),
+            dict(sorted(legacy.count_ops().items())),
+        )
+        self.assertEqual(transpiled_circuit.size(lambda x: x.operation.num_qubits == 2), 1)
+
+    def test_two_qubit_rzz_cz_gates_rzz_target(self):
+        """Test the synthesis of a circuit containing a RZZ and CZ gates
+        on a target with RZZ and CZ gates"""
+        theta = Parameter("θ")
+        lam = Parameter("λ")
+        phi = Parameter("ϕ")
+        target = Target(num_qubits=2)
+        target.add_instruction(RXGate(lam))
+        target.add_instruction(RZGate(theta))
+        target.add_instruction(RZZGate(phi))
+        target.add_instruction(CZGate())
+
+        qc = QuantumCircuit(2)
+        qc.rzz(0.2, 0, 1)
+        qc.cz(0, 1)
+
+        peephole = TwoQubitPeepholeOptimization(target)
+        transpiled_circuit = peephole(qc)
+        self.all_inst_in_target(transpiled_circuit, target)
+        self.assertEqual(Operator(transpiled_circuit), Operator(qc))
+        self.assertTrue(set(transpiled_circuit.count_ops()).issubset({"rz", "rx", "rzz"}))
+        self.assertEqual(transpiled_circuit.count_ops()["rzz"], 1)
+
+    @combine(
+        target_2q_gates=[
+            ["cx"],
+            ["cz"],
+            ["ecr"],
+            ["cx", "cz", "ecr"],
+            ["rxx", "cx"],
+            ["rzz", "cz"],
+            ["rzx", "ecr"],
+            ["rzz", "rzx", "rxx", "ecr", "cx", "cz"],
+            ["cz", "swap"],
+            ["cx", "swap"],
+        ],
+        target_1q_gates=[
+            ["rz", "sx"],
+            ["rz", "sx", "x"],
+            ["rz", "rx", "sx", "x"],
+            ["ry", "rx"],
+            ["ry", "rx", "rz"],
+            ["rz", "rx"],
+            ["u"],
+            ["p", "sx"],
+            ["u", "u1", "u2", "u3", "rz", "sx", "x", "ry", "rx", "r", "p"],
+        ],
+        name="{target_2q_gates}_{target_1q_gates}",
+    )
+    def test_random_circuit(self, target_2q_gates, target_1q_gates):
+        cmap = CouplingMap.from_grid(2, 3)
+        basis_gates = target_1q_gates + target_2q_gates
+        backend = GenericBackendV2(
+            cmap.size(), basis_gates=basis_gates, seed=2024, coupling_map=cmap
+        )
+        peephole = TwoQubitPeepholeOptimization(backend.target)
+        qc = random_circuit(cmap.size(), 79, max_operands=2, seed=12345_42)
+        qc = transpile(qc, target=backend.target, optimization_level=0, seed_transpiler=2025)
+        result = peephole(qc)
+        if "rzx" in target_2q_gates or "rzz" in target_2q_gates or "rxx" in target_2q_gates:
+            self.all_inst_in_target(result, backend.target, allow_inverse=True)
+        else:
+            self.all_inst_in_target(result, backend.target)
+        np.testing.assert_allclose(Operator(result), Operator(qc), atol=1e-12, rtol=0)
diff --git a/test/python/transpiler/test_vf2_post_layout.py b/test/python/transpiler/test_vf2_post_layout.py
index ea5c5f15bf2a..1e6bb82af604 100644
--- a/test/python/transpiler/test_vf2_post_layout.py
+++ b/test/python/transpiler/test_vf2_post_layout.py
@@ -567,11 +567,10 @@ def test_2q_circuit_5q_backend_v2(self):
             coupling_map=YORKTOWN_CMAP,
             seed=self.seed,
         )
-
         qr = QuantumRegister(2, "qr")
         circuit = QuantumCircuit(qr)
         circuit.cx(qr[1], qr[0])  # qr1 -> qr0
-        tqc = transpile(circuit, backend, layout_method="dense")
+        tqc = transpile(circuit, backend, initial_layout=[1, 0])
         initial_layout = tqc._layout
         dag = circuit_to_dag(tqc)
         pass_ = VF2PostLayout(target=backend.target, seed=self.seed, strict_direction=False)