diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml
index 0169e3ed4cb4..a9c1839487a3 100644
--- a/.github/workflows/coverage.yml
+++ b/.github/workflows/coverage.yml
@@ -71,7 +71,7 @@ jobs:
           lcov --add-tracefile python.info --add-tracefile rust.info --output-file coveralls.info
 
       - name: Coveralls
-        uses: coverallsapp/github-action@v2.3.0
+        uses: coverallsapp/github-action@v2.3.4
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           path-to-lcov: coveralls.info
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 20c40dc38321..39439f7dd059 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -20,7 +20,7 @@ jobs:
       matrix:
         # Normally we test min and max version but we can't run python
         # 3.9 on arm64 until actions/setup-python#808 is resolved
-        python-version: ["3.10", "3.12"]
+        python-version: ["3.10", "3.13"]
     steps:
       - uses: actions/checkout@v4
       - name: Install Rust toolchain
@@ -44,7 +44,7 @@ jobs:
           python -m pip install -U -r requirements.txt -c constraints.txt
           python -m pip install -U -r requirements-dev.txt -c constraints.txt
           python -m pip install -c constraints.txt -e .
-        if: matrix.python-version == '3.12'
+        if: matrix.python-version == '3.13'
       - name: 'Install optionals'
         run: |
           python -m pip install -r requirements-optional.txt -c constraints.txt
diff --git a/.github/workflows/wheels-build.yml b/.github/workflows/wheels-build.yml
index 30fdcd84bbb5..a7453eb05568 100644
--- a/.github/workflows/wheels-build.yml
+++ b/.github/workflows/wheels-build.yml
@@ -27,7 +27,7 @@ on:
         default: "default"
         required: false
 
-      wheels-32bit: 
+      wheels-32bit:
         description: >-
           The action to take for Tier 1 wheels.
           Choose from 'default', 'build' or 'skip'.
@@ -36,7 +36,7 @@ on:
         default: "default"
         required: false
 
-      wheels-linux-s390x: 
+      wheels-linux-s390x:
         description: >-
           The action to take for Linux s390x wheels.
           Choose from 'default', 'build' or 'skip'.
@@ -44,7 +44,7 @@ on:
         default: "default"
         required: false
 
-      wheels-linux-ppc64le: 
+      wheels-linux-ppc64le:
         description: >-
           The action to take for Linux ppc64le wheels.
           Choose from 'default', 'build' or 'skip'.
@@ -52,7 +52,7 @@ on:
         default: "default"
         required: false
 
-      wheels-linux-aarch64: 
+      wheels-linux-aarch64:
         description: >-
           The action to take for Linux AArch64 wheels.
           Choose from 'default', 'build' or 'skip'.
@@ -77,7 +77,7 @@ on:
         type: boolean
         default: true
         required: false
-      
+
 
 jobs:
   wheels-tier-1:
@@ -126,7 +126,7 @@ jobs:
         env:
           PGO_WORK_DIR: ${{ github.workspace }}/pgo-data
           PGO_OUT_PATH: ${{ github.workspace }}/merged.profdata
-      - uses: pypa/cibuildwheel@v2.19.2
+      - uses: pypa/cibuildwheel@v2.21.3
       - uses: actions/upload-artifact@v4
         with:
           path: ./wheelhouse/*.whl
@@ -151,7 +151,7 @@ jobs:
         with:
           components: llvm-tools-preview
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.19.2
+        uses: pypa/cibuildwheel@v2.21.3
         env:
           CIBW_SKIP: 'pp* cp36-* cp37-* cp38-* *musllinux* *amd64 *x86_64'
       - uses: actions/upload-artifact@v4
@@ -173,7 +173,7 @@ jobs:
       - uses: docker/setup-qemu-action@v3
         with:
           platforms: all
-      - uses: pypa/cibuildwheel@v2.19.2
+      - uses: pypa/cibuildwheel@v2.21.3
         env:
           CIBW_ARCHS_LINUX: s390x
           CIBW_TEST_SKIP: "cp*"
@@ -196,7 +196,7 @@ jobs:
       - uses: docker/setup-qemu-action@v3
         with:
           platforms: all
-      - uses: pypa/cibuildwheel@v2.19.2
+      - uses: pypa/cibuildwheel@v2.21.3
         env:
           CIBW_ARCHS_LINUX: ppc64le
           CIBW_TEST_SKIP: "cp*"
@@ -218,7 +218,7 @@ jobs:
       - uses: docker/setup-qemu-action@v3
         with:
           platforms: all
-      - uses: pypa/cibuildwheel@v2.19.2
+      - uses: pypa/cibuildwheel@v2.21.3
         env:
           CIBW_ARCHS_LINUX: aarch64
           CIBW_TEST_COMMAND: cp -r {project}/test . && QISKIT_PARALLEL=FALSE stestr --test-path test/python run --abbreviate -n test.python.compiler.test_transpiler
diff --git a/asv.conf.json b/asv.conf.json
index a75bc0c59c3a..379a01cd7edd 100644
--- a/asv.conf.json
+++ b/asv.conf.json
@@ -17,7 +17,7 @@
     "dvcs": "git",
     "environment_type": "virtualenv",
     "show_commit_url": "http://github.com/Qiskit/qiskit/commit/",
-    "pythons": ["3.9", "3.10", "3.11", "3.12"],
+    "pythons": ["3.9", "3.10", "3.11", "3.12", "3.13"],
     "benchmark_dir": "test/benchmarks",
     "env_dir": ".asv/env",
     "results_dir": ".asv/results"
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 49df19808498..d182ba32f1b3 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -37,7 +37,7 @@ parameters:
   - name: "supportedPythonVersions"
     displayName: "All supported versions of Python"
     type: object
-    default: ["3.9", "3.10", "3.11", "3.12"]
+    default: ["3.9", "3.10", "3.11", "3.12", "3.13"]
 
   - name: "minimumPythonVersion"
     displayName: "Minimum supported version of Python"
@@ -47,7 +47,7 @@ parameters:
   - name: "maximumPythonVersion"
     displayName: "Maximum supported version of Python"
     type: string
-    default: "3.12"
+    default: "3.13"
 
   # These two versions of Python can be chosen somewhat arbitrarily, but we get
   # slightly better coverage per PR if they're neither the maximum nor minimum
diff --git a/constraints.txt b/constraints.txt
index cef0c15114b7..42e11f9e1460 100644
--- a/constraints.txt
+++ b/constraints.txt
@@ -3,10 +3,6 @@
 # https://github.com/Qiskit/qiskit-terra/issues/10345 for current details.
 scipy<1.11; python_version<'3.12'
 
-# Temporary pin to avoid CI issues caused by scipy 1.14.0
-# See https://github.com/Qiskit/qiskit/issues/12655 for current details.
-scipy==1.13.1; python_version=='3.12'
-
 # z3-solver from 4.12.3 onwards upped the minimum macOS API version for its
 # wheels to 11.7. The Azure VM images contain pre-built CPythons, of which at
 # least CPython 3.8 was compiled for an older macOS, so does not match a
diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs
index ed3b75d309d6..e04f6d63a936 100644
--- a/crates/accelerate/src/lib.rs
+++ b/crates/accelerate/src/lib.rs
@@ -37,6 +37,7 @@ pub mod nlayout;
 pub mod optimize_1q_gates;
 pub mod pauli_exp_val;
 pub mod remove_diagonal_gates_before_measure;
+pub mod remove_identity_equiv;
 pub mod results;
 pub mod sabre;
 pub mod sampled_exp_val;
diff --git a/crates/accelerate/src/remove_identity_equiv.rs b/crates/accelerate/src/remove_identity_equiv.rs
new file mode 100644
index 000000000000..a3eb921628e2
--- /dev/null
+++ b/crates/accelerate/src/remove_identity_equiv.rs
@@ -0,0 +1,149 @@
+// 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 num_complex::Complex64;
+use num_complex::ComplexFloat;
+use pyo3::prelude::*;
+use rustworkx_core::petgraph::stable_graph::NodeIndex;
+
+use crate::nlayout::PhysicalQubit;
+use crate::target_transpiler::Target;
+use qiskit_circuit::dag_circuit::DAGCircuit;
+use qiskit_circuit::operations::Operation;
+use qiskit_circuit::operations::OperationRef;
+use qiskit_circuit::operations::Param;
+use qiskit_circuit::operations::StandardGate;
+use qiskit_circuit::packed_instruction::PackedInstruction;
+
+#[pyfunction]
+#[pyo3(signature=(dag, approx_degree=Some(1.0), target=None))]
+fn remove_identity_equiv(
+    dag: &mut DAGCircuit,
+    approx_degree: Option<f64>,
+    target: Option<&Target>,
+) {
+    let mut remove_list: Vec<NodeIndex> = Vec::new();
+
+    let get_error_cutoff = |inst: &PackedInstruction| -> f64 {
+        match approx_degree {
+            Some(degree) => {
+                if degree == 1.0 {
+                    f64::EPSILON
+                } else {
+                    match target {
+                        Some(target) => {
+                            let qargs: Vec<PhysicalQubit> = dag
+                                .get_qargs(inst.qubits)
+                                .iter()
+                                .map(|x| PhysicalQubit::new(x.0))
+                                .collect();
+                            let error_rate = target.get_error(inst.op.name(), qargs.as_slice());
+                            match error_rate {
+                                Some(err) => err * degree,
+                                None => f64::EPSILON.max(1. - degree),
+                            }
+                        }
+                        None => f64::EPSILON.max(1. - degree),
+                    }
+                }
+            }
+            None => match target {
+                Some(target) => {
+                    let qargs: Vec<PhysicalQubit> = dag
+                        .get_qargs(inst.qubits)
+                        .iter()
+                        .map(|x| PhysicalQubit::new(x.0))
+                        .collect();
+                    let error_rate = target.get_error(inst.op.name(), qargs.as_slice());
+                    match error_rate {
+                        Some(err) => err,
+                        None => f64::EPSILON,
+                    }
+                }
+                None => f64::EPSILON,
+            },
+        }
+    };
+
+    for op_node in dag.op_nodes(false) {
+        let inst = dag.dag()[op_node].unwrap_operation();
+        match inst.op.view() {
+            OperationRef::Standard(gate) => {
+                let (dim, trace) = match gate {
+                    StandardGate::RXGate | StandardGate::RYGate | StandardGate::RZGate => {
+                        if let Param::Float(theta) = inst.params_view()[0] {
+                            let trace = (theta / 2.).cos() * 2.;
+                            (2., trace)
+                        } else {
+                            continue;
+                        }
+                    }
+                    StandardGate::RXXGate
+                    | StandardGate::RYYGate
+                    | StandardGate::RZZGate
+                    | StandardGate::RZXGate => {
+                        if let Param::Float(theta) = inst.params_view()[0] {
+                            let trace = (theta / 2.).cos() * 4.;
+                            (4., trace)
+                        } else {
+                            continue;
+                        }
+                    }
+                    _ => {
+                        // Skip global phase gate
+                        if gate.num_qubits() < 1 {
+                            continue;
+                        }
+                        if let Some(matrix) = gate.matrix(inst.params_view()) {
+                            let dim = matrix.shape()[0] as f64;
+                            let trace = matrix.diag().iter().sum::<Complex64>().abs();
+                            (dim, trace)
+                        } else {
+                            continue;
+                        }
+                    }
+                };
+                let error = get_error_cutoff(inst);
+                let f_pro = (trace / dim).powi(2);
+                let gate_fidelity = (dim * f_pro + 1.) / (dim + 1.);
+                if (1. - gate_fidelity).abs() < error {
+                    remove_list.push(op_node)
+                }
+            }
+            OperationRef::Gate(gate) => {
+                // Skip global phase like gate
+                if gate.num_qubits() < 1 {
+                    continue;
+                }
+                if let Some(matrix) = gate.matrix(inst.params_view()) {
+                    let error = get_error_cutoff(inst);
+                    let dim = matrix.shape()[0] as f64;
+                    let trace: Complex64 = matrix.diag().iter().sum();
+                    let f_pro = (trace / dim).abs().powi(2);
+                    let gate_fidelity = (dim * f_pro + 1.) / (dim + 1.);
+                    if (1. - gate_fidelity).abs() < error {
+                        remove_list.push(op_node)
+                    }
+                }
+            }
+            _ => continue,
+        }
+    }
+    for node in remove_list {
+        dag.remove_op_node(node);
+    }
+}
+
+pub fn remove_identity_equiv_mod(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(remove_identity_equiv))?;
+    Ok(())
+}
diff --git a/crates/accelerate/src/target_transpiler/mod.rs b/crates/accelerate/src/target_transpiler/mod.rs
index 54ab34341a4c..9e6c220818a3 100644
--- a/crates/accelerate/src/target_transpiler/mod.rs
+++ b/crates/accelerate/src/target_transpiler/mod.rs
@@ -21,7 +21,7 @@ use std::ops::Index;
 use ahash::RandomState;
 
 use hashbrown::HashSet;
-use indexmap::{IndexMap, IndexSet};
+use indexmap::IndexMap;
 use itertools::Itertools;
 use nullable_index_map::NullableIndexMap;
 use pyo3::{
@@ -57,7 +57,7 @@ type GateMapState = Vec<(String, Vec<(Option<Qargs>, Option<InstructionPropertie
 
 /// Represents a Qiskit `Gate` object or a Variadic instruction.
 /// Keeps a reference to its Python instance for caching purposes.
-#[derive(Debug, Clone, FromPyObject)]
+#[derive(FromPyObject, Debug, Clone)]
 pub(crate) enum TargetOperation {
     Normal(NormalOperation),
     Variadic(PyObject),
@@ -82,19 +82,23 @@ impl ToPyObject for TargetOperation {
 }
 
 impl TargetOperation {
-    fn num_qubits(&self) -> u32 {
+    /// Gets the number of qubits of a [TargetOperation], will panic if the operation is [TargetOperation::Variadic].
+    pub fn num_qubits(&self) -> u32 {
         match &self {
-            Self::Normal(normal) => normal.operation.view().num_qubits(),
+            Self::Normal(normal) => normal.operation.num_qubits(),
             Self::Variadic(_) => {
-                unreachable!("'num_qubits' property is reserved for normal operations only.")
+                panic!("'num_qubits' property doesn't exist for Variadic operations")
             }
         }
     }
 
-    fn params(&self) -> &[Param] {
+    /// Gets the parameters of a [TargetOperation], will panic if the operation is [TargetOperation::Variadic].
+    pub fn params(&self) -> &[Param] {
         match &self {
             TargetOperation::Normal(normal) => normal.params.as_slice(),
-            TargetOperation::Variadic(_) => &[],
+            TargetOperation::Variadic(_) => {
+                panic!("'parameters' property doesn't exist for Variadic operations")
+            }
         }
     }
 }
@@ -173,7 +177,6 @@ pub(crate) struct Target {
     #[pyo3(get)]
     _gate_name_map: IndexMap<String, TargetOperation, RandomState>,
     global_operations: IndexMap<u32, HashSet<String>, RandomState>,
-    variable_class_operations: IndexSet<String, RandomState>,
     qarg_gate_map: NullableIndexMap<Qargs, Option<HashSet<String>>>,
     non_global_strict_basis: Option<Vec<String>>,
     non_global_basis: Option<Vec<String>>,
@@ -269,7 +272,6 @@ impl Target {
             concurrent_measurements,
             gate_map: GateMap::default(),
             _gate_name_map: IndexMap::default(),
-            variable_class_operations: IndexSet::default(),
             global_operations: IndexMap::default(),
             qarg_gate_map: NullableIndexMap::default(),
             non_global_basis: None,
@@ -302,16 +304,15 @@ impl Target {
             )));
         }
         let mut qargs_val: PropsMap;
-        match instruction {
+        match &instruction {
             TargetOperation::Variadic(_) => {
                 qargs_val = PropsMap::with_capacity(1);
                 qargs_val.extend([(None, None)]);
-                self.variable_class_operations.insert(name.to_string());
             }
-            TargetOperation::Normal(_) => {
+            TargetOperation::Normal(normal) => {
                 if let Some(mut properties) = properties {
                     qargs_val = PropsMap::with_capacity(properties.len());
-                    let inst_num_qubits = instruction.num_qubits();
+                    let inst_num_qubits = normal.operation.view().num_qubits();
                     if properties.contains_key(None) {
                         self.global_operations
                             .entry(inst_num_qubits)
@@ -619,7 +620,7 @@ impl Target {
         } else if let Some(operation_name) = operation_name {
             if let Some(parameters) = parameters {
                 if let Some(obj) = self._gate_name_map.get(&operation_name) {
-                    if self.variable_class_operations.contains(&operation_name) {
+                    if matches!(obj, TargetOperation::Variadic(_)) {
                         if let Some(_qargs) = qargs {
                             let qarg_set: HashSet<PhysicalQubit> = _qargs.iter().cloned().collect();
                             return Ok(_qargs
@@ -1053,8 +1054,8 @@ impl Target {
         if let Some(Some(qarg_gate_map_arg)) = self.qarg_gate_map.get(qargs).as_ref() {
             res.extend(qarg_gate_map_arg.iter().map(|key| key.as_str()));
         }
-        for name in self._gate_name_map.keys() {
-            if self.variable_class_operations.contains(name) {
+        for (name, obj) in self._gate_name_map.iter() {
+            if matches!(obj, TargetOperation::Variadic(_)) {
                 res.insert(name);
             }
         }
@@ -1160,34 +1161,40 @@ impl Target {
                     }
                     if gate_prop_name.contains_key(None) {
                         let obj = &self._gate_name_map[operation_name];
-                        if self.variable_class_operations.contains(operation_name) {
+                        match obj {
+                            TargetOperation::Variadic(_) => {
+                                return qargs.is_none()
+                                    || _qargs.iter().all(|qarg| {
+                                        qarg.index() <= self.num_qubits.unwrap_or_default()
+                                    }) && qarg_set.len() == _qargs.len();
+                            }
+                            TargetOperation::Normal(obj) => {
+                                let qubit_comparison = obj.operation.num_qubits();
+                                return qubit_comparison == _qargs.len() as u32
+                                    && _qargs.iter().all(|qarg| {
+                                        qarg.index() < self.num_qubits.unwrap_or_default()
+                                    });
+                            }
+                        }
+                    }
+                } else {
+                    // Duplicate case is if it contains none
+                    let obj = &self._gate_name_map[operation_name];
+                    match obj {
+                        TargetOperation::Variadic(_) => {
                             return qargs.is_none()
                                 || _qargs.iter().all(|qarg| {
                                     qarg.index() <= self.num_qubits.unwrap_or_default()
                                 }) && qarg_set.len() == _qargs.len();
-                        } else {
-                            let qubit_comparison = obj.num_qubits();
+                        }
+                        TargetOperation::Normal(obj) => {
+                            let qubit_comparison = obj.operation.num_qubits();
                             return qubit_comparison == _qargs.len() as u32
                                 && _qargs.iter().all(|qarg| {
                                     qarg.index() < self.num_qubits.unwrap_or_default()
                                 });
                         }
                     }
-                } else {
-                    // Duplicate case is if it contains none
-                    if self.variable_class_operations.contains(operation_name) {
-                        return qargs.is_none()
-                            || _qargs
-                                .iter()
-                                .all(|qarg| qarg.index() <= self.num_qubits.unwrap_or_default())
-                                && qarg_set.len() == _qargs.len();
-                    } else {
-                        let qubit_comparison = self._gate_name_map[operation_name].num_qubits();
-                        return qubit_comparison == _qargs.len() as u32
-                            && _qargs
-                                .iter()
-                                .all(|qarg| qarg.index() < self.num_qubits.unwrap_or_default());
-                    }
                 }
             } else {
                 return true;
diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs
index 560541455138..49070e85db70 100644
--- a/crates/pyext/src/lib.rs
+++ b/crates/pyext/src/lib.rs
@@ -50,6 +50,7 @@ fn _accelerate(m: &Bound<PyModule>) -> PyResult<()> {
     add_submodule(m, ::qiskit_accelerate::optimize_1q_gates::optimize_1q_gates, "optimize_1q_gates")?;
     add_submodule(m, ::qiskit_accelerate::pauli_exp_val::pauli_expval, "pauli_expval")?;
     add_submodule(m, ::qiskit_accelerate::remove_diagonal_gates_before_measure::remove_diagonal_gates_before_measure, "remove_diagonal_gates_before_measure")?;
+    add_submodule(m, ::qiskit_accelerate::remove_identity_equiv::remove_identity_equiv_mod, "remove_identity_equiv")?;
     add_submodule(m, ::qiskit_accelerate::results::results, "results")?;
     add_submodule(m, ::qiskit_accelerate::sabre::sabre, "sabre")?;
     add_submodule(m, ::qiskit_accelerate::sampled_exp_val::sampled_exp_val, "sampled_exp_val")?;
diff --git a/pyproject.toml b/pyproject.toml
index 70a6e7c9d0cf..326d28f2a273 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,6 +31,7 @@ classifiers = [
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Topic :: Scientific/Engineering",
 ]
 # These are configured in the `tool.setuptools.dynamic` table.
diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 63e5a2a48a47..202ebc32be85 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -106,6 +106,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.remove_identity_equiv"] = _accelerate.remove_identity_equiv
 
 from qiskit.exceptions import QiskitError, MissingOptionalLibraryError
 
diff --git a/qiskit/quantum_info/operators/channel/transformations.py b/qiskit/quantum_info/operators/channel/transformations.py
index 8f429cad8cea..657ee62703ec 100644
--- a/qiskit/quantum_info/operators/channel/transformations.py
+++ b/qiskit/quantum_info/operators/channel/transformations.py
@@ -220,32 +220,39 @@ def _kraus_to_choi(data):
 
 def _choi_to_kraus(data, input_dim, output_dim, atol=ATOL_DEFAULT):
     """Transform Choi representation to Kraus representation."""
-    from scipy import linalg as la
+    import scipy.linalg
 
     # Check if hermitian matrix
     if is_hermitian_matrix(data, atol=atol):
-        # Get eigen-decomposition of Choi-matrix
-        # This should be a call to la.eigh, but there is an OpenBlas
-        # threading issue that is causing segfaults.
-        # Need schur here since la.eig does not
-        # guarantee orthogonality in degenerate subspaces
-        w, v = la.schur(data, output="complex")
-        w = w.diagonal().real
-        # Check eigenvalues are non-negative
-        if len(w[w < -atol]) == 0:
-            # CP-map Kraus representation
-            kraus = []
-            for val, vec in zip(w, v.T):
-                if abs(val) > atol:
-                    k = np.sqrt(val) * vec.reshape((output_dim, input_dim), order="F")
-                    kraus.append(k)
-            # If we are converting a zero matrix, we need to return a Kraus set
-            # with a single zero-element Kraus matrix
+        # Ideally we'd use `eigh`, but `scipy.linalg.eigh` has stability problems on macOS (at a
+        # minimum from SciPy 1.1 to 1.13 with the bundled OpenBLAS, or ~0.3.6 before they started
+        # bundling one in).  The Schur form of a Hermitian matrix is guaranteed diagonal:
+        #
+        #       H = U T U+   for upper-triangular T.
+        #   => H+ = U T+ U+
+        #   =>  T = T+       because H = H+, and thus T cannot have super-diagonal elements.
+        #
+        # So the eigenvalues are on the diagonal, therefore the basis-transformation matrix must be
+        # a spanning set of the eigenspace.
+        triangular, vecs = scipy.linalg.schur(data)
+        values = triangular.diagonal().real
+        # If we're not a CP map, fall-through back to the generalization handling.  Since we needed
+        # to get the eigenvalues anyway, we can do the CP check manually rather than deferring to a
+        # separate re-calculation.
+        if all(values >= -atol):
+            kraus = [
+                math.sqrt(value) * vec.reshape((output_dim, input_dim), order="F")
+                for value, vec in zip(values, vecs.T)
+                if abs(value) > atol
+            ]
+            # If we are converting a zero matrix, we need to return a Kraus set with a single
+            # zero-element Kraus matrix
             if not kraus:
-                kraus.append(np.zeros((output_dim, input_dim), dtype=complex))
+                kraus = [np.zeros((output_dim, input_dim), dtype=complex)]
             return kraus, None
-    # Non-CP-map generalized Kraus representation
-    mat_u, svals, mat_vh = la.svd(data)
+        # Fall through.
+    # Non-CP-map generalized Kraus representation.
+    mat_u, svals, mat_vh = scipy.linalg.svd(data)
     kraus_l = []
     kraus_r = []
     for val, vec_l, vec_r in zip(svals, mat_u.T, mat_vh.conj()):
diff --git a/qiskit/quantum_info/operators/operator.py b/qiskit/quantum_info/operators/operator.py
index 42593626f2cc..fc3dd9ee187a 100644
--- a/qiskit/quantum_info/operators/operator.py
+++ b/qiskit/quantum_info/operators/operator.py
@@ -16,6 +16,7 @@
 
 from __future__ import annotations
 
+import cmath
 import copy as _copy
 import re
 from numbers import Number
@@ -540,11 +541,37 @@ def compose(self, other: Operator, qargs: list | None = None, front: bool = Fals
         ret._op_shape = new_shape
         return ret
 
-    def power(self, n: float) -> Operator:
+    def power(self, n: float, branch_cut_rotation=cmath.pi * 1e-12) -> Operator:
         """Return the matrix power of the operator.
 
+        Non-integer powers of operators with an eigenvalue whose complex phase is :math:`\\pi` have
+        a branch cut in the complex plane, which makes the calculation of the principal root around
+        this cut subject to precision / differences in BLAS implementation.  For example, the square
+        root of Pauli Y can return the :math:`\\pi/2` or :math:`\\-pi/2` Y rotation depending on
+        whether the -1 eigenvalue is found as ``complex(-1, tiny)`` or ``complex(-1, -tiny)``. Such
+        eigenvalues are really common in quantum information, so this function first phase-rotates
+        the input matrix to shift the branch cut to a far less common point.  The underlying
+        numerical precision issues around the branch-cut point remain, if an operator has an
+        eigenvalue close to this phase.  The magnitude of this rotation can be controlled with the
+        ``branch_cut_rotation`` parameter.
+
+        The choice of ``branch_cut_rotation`` affects the principal root that is found.  For
+        example, the square root of :class:`.ZGate` will be calculated as either :class:`.SGate` or
+        :class:`.SdgGate` depending on which way the rotation is done::
+
+            from qiskit.circuit import library
+            from qiskit.quantum_info import Operator
+
+            z_op = Operator(library.ZGate())
+            assert z_op.power(0.5, branch_cut_rotation=1e-3) == Operator(library.SGate())
+            assert z_op.power(0.5, branch_cut_rotation=-1e-3) == Operator(library.SdgGate())
+
         Args:
             n (float): the power to raise the matrix to.
+            branch_cut_rotation (float): The rotation angle to apply to the branch cut in the
+                complex plane.  This shifts the branch cut away from the common point of :math:`-1`,
+                but can cause a different root to be selected as the principal root.  The rotation
+                is anticlockwise, following the standard convention for complex phase.
 
         Returns:
             Operator: the resulting operator ``O ** n``.
@@ -561,13 +588,11 @@ def power(self, n: float) -> Operator:
         else:
             import scipy.linalg
 
-            # Experimentally, for fractional powers this seems to be 3x faster than
-            # calling scipy.linalg.fractional_matrix_power(self.data, n)
-            decomposition, unitary = scipy.linalg.schur(self.data, output="complex")
-            decomposition_diagonal = decomposition.diagonal()
-            decomposition_power = [pow(element, n) for element in decomposition_diagonal]
-            unitary_power = unitary @ np.diag(decomposition_power) @ unitary.conj().T
-            ret._data = unitary_power
+            ret._data = cmath.rect(
+                1, branch_cut_rotation * n
+            ) * scipy.linalg.fractional_matrix_power(
+                cmath.rect(1, -branch_cut_rotation) * self.data, n
+            )
         return ret
 
     def tensor(self, other: Operator) -> Operator:
diff --git a/qiskit/transpiler/passes/__init__.py b/qiskit/transpiler/passes/__init__.py
index ca4e3545a98b..5bc1ae555a5e 100644
--- a/qiskit/transpiler/passes/__init__.py
+++ b/qiskit/transpiler/passes/__init__.py
@@ -93,6 +93,7 @@
    NormalizeRXAngle
    OptimizeAnnotated
    Split2QUnitaries
+   RemoveIdentityEquivalent
 
 Calibration
 =============
@@ -247,6 +248,7 @@
 from .optimization import ElidePermutations
 from .optimization import NormalizeRXAngle
 from .optimization import OptimizeAnnotated
+from .optimization import RemoveIdentityEquivalent
 from .optimization import Split2QUnitaries
 
 # circuit analysis
diff --git a/qiskit/transpiler/passes/optimization/__init__.py b/qiskit/transpiler/passes/optimization/__init__.py
index 8e2883b27781..c0e455b2065b 100644
--- a/qiskit/transpiler/passes/optimization/__init__.py
+++ b/qiskit/transpiler/passes/optimization/__init__.py
@@ -38,5 +38,6 @@
 from .elide_permutations import ElidePermutations
 from .normalize_rx_angle import NormalizeRXAngle
 from .optimize_annotated import OptimizeAnnotated
+from .remove_identity_equiv import RemoveIdentityEquivalent
 from .split_2q_unitaries import Split2QUnitaries
 from .collect_and_collapse import CollectAndCollapse
diff --git a/qiskit/transpiler/passes/optimization/remove_identity_equiv.py b/qiskit/transpiler/passes/optimization/remove_identity_equiv.py
new file mode 100644
index 000000000000..4f8551f12442
--- /dev/null
+++ b/qiskit/transpiler/passes/optimization/remove_identity_equiv.py
@@ -0,0 +1,69 @@
+# 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.
+
+"""Transpiler pass to drop gates with negligible effects."""
+
+from __future__ import annotations
+
+from qiskit.dagcircuit import DAGCircuit
+from qiskit.transpiler.target import Target
+from qiskit.transpiler.basepasses import TransformationPass
+from qiskit._accelerate.remove_identity_equiv import remove_identity_equiv
+
+
+class RemoveIdentityEquivalent(TransformationPass):
+    r"""Remove gates with negligible effects.
+
+    Removes gates whose effect is close to an identity operation, up to the specified
+    tolerance. Zero qubit gates such as :class:`.GlobalPhaseGate` are not considered
+    by this pass.
+
+    For a cutoff fidelity :math:`f`, this pass removes gates whose average
+    gate fidelity with respect to the identity is below :math:`f`. Concretely,
+    a gate :math:`G` is removed if :math:`\bar F < f` where
+
+    .. math::
+
+        \bar{F} = \frac{1 + F_{\text{process}}{1 + d}
+
+        F_{\text{process}} = \frac{|\mathrm{Tr}(G)|^2}{d^2}
+
+    where :math:`d = 2^n` is the dimension of the gate for :math:`n` qubits.
+    """
+
+    def __init__(
+        self, *, approximation_degree: float | None = 1.0, target: None | Target = None
+    ) -> None:
+        """Initialize the transpiler pass.
+
+        Args:
+            approximation_degree: The degree to approximate for the equivalence check. This can be a
+                floating point value between 0 and 1, or ``None``. If the value is 1 this does not
+                approximate above floating point precision. For a value < 1 this is used as a scaling
+                factor for the cutoff fidelity. If the value is ``None`` this approximates up to the
+                fidelity for the gate specified in ``target``.
+
+            target: If ``approximation_degree`` is set to ``None`` and a :class:`.Target` is provided
+                for this field the tolerance for determining whether an operation is equivalent to
+                identity will be set to the reported error rate in the target. If
+                ``approximation_degree`` (the default) this has no effect, if
+                ``approximation_degree=None`` it uses the error rate specified in the ``Target`` for
+                the gate being evaluated, and a numeric value other than 1 with ``target`` set is
+                used as a scaling factor of the target's error rate.
+        """
+        super().__init__()
+        self._approximation_degree = approximation_degree
+        self._target = target
+
+    def run(self, dag: DAGCircuit) -> DAGCircuit:
+        remove_identity_equiv(dag, self._approximation_degree, self._target)
+        return dag
diff --git a/qiskit/visualization/dag_visualization.py b/qiskit/visualization/dag_visualization.py
index a6a111fe75e9..8f8b8fc89097 100644
--- a/qiskit/visualization/dag_visualization.py
+++ b/qiskit/visualization/dag_visualization.py
@@ -81,7 +81,7 @@ def dag_drawer(dag, scale=0.7, filename=None, style="color"):
     ``rustworkx`` package to draw the DAG.
 
     Args:
-        dag (DAGCircuit): The dag to draw.
+        dag (DAGCircuit or DAGDependency): The dag to draw.
         scale (float): scaling factor
         filename (str): file path to save image to (format inferred from name)
         style (str): 'plain': B&W graph
diff --git a/releasenotes/notes/remove_identity_equiv-9c627c8c35b2298a.yaml b/releasenotes/notes/remove_identity_equiv-9c627c8c35b2298a.yaml
new file mode 100644
index 000000000000..90d016803d62
--- /dev/null
+++ b/releasenotes/notes/remove_identity_equiv-9c627c8c35b2298a.yaml
@@ -0,0 +1,29 @@
+---
+features_transpiler:
+  - |
+    Added a new transpiler pass, :class:`.RemoveIdentityEquivalent` that is used
+    to remove gates that are equivalent to an identity up to some tolerance.
+    For example if you had a circuit like:
+
+        .. plot::
+
+           from qiskit.circuit import QuantumCircuit
+
+           qc = QuantumCircuit(2)
+           qc.cp(1e-20, 0, 1)
+           qc.draw("mpl")
+
+    running the pass would eliminate the :class:`.CPhaseGate`:
+
+        .. plot::
+           :include-source:
+
+           from qiskit.circuit import QuantumCircuit
+           from qiskit.transpiler.passes import RemoveIdentityEquivalent
+
+           qc = QuantumCircuit(2)
+           qc.cp(1e-20, 0, 1)
+
+           removal_pass = RemoveIdentityEquivalent()
+           result = removal_pass(qc)
+           result.draw("mpl")
diff --git a/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml
new file mode 100644
index 000000000000..41f4b8790286
--- /dev/null
+++ b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml
@@ -0,0 +1,25 @@
+---
+fixes:
+  - |
+    Fixed :meth:`.Operator.power` when called with non-integer powers on a matrix whose Schur form
+    is not diagonal (for example, most non-unitary matrices).
+  - |
+    :meth:`.Operator.power` will now more reliably return the expected principal value from a
+    fractional matrix power of a unitary matrix with a :math:`-1` eigenvalue.  This is tricky in
+    general, because floating-point rounding effects can cause a matrix to _truly_ have an eigenvalue
+    on the negative side of the branch cut (even if its exact mathematical relation would not), and
+    imprecision in various BLAS calls can falsely find the wrong side of the branch cut.
+
+    :meth:`.Operator.power` now shifts the branch-cut location for matrix powers to be a small
+    complex rotation away from :math:`-1`.  This does not solve the problem, it just shifts it to a
+    place where it is far less likely to be noticeable for the types of operators that usually
+    appear.  Use the new ``branch_cut_rotation`` parameter to have more control over this.
+
+    See `#13305 <https://github.com/Qiskit/qiskit/issues/13305>`__.
+features_quantum_info:
+  - |
+    The method :meth:`.Operator.power` has a new parameter ``branch_cut_rotation``.  This can be
+    used to shift the branch-cut point of the root around, which can affect which matrix is chosen
+    as the principal root.  By default, it is set to a small positive rotation to make roots of
+    operators with a real-negative eigenvalue (like Pauli operators) more stable against numerical
+    precision differences.
diff --git a/test/python/quantum_info/operators/channel/test_kraus.py b/test/python/quantum_info/operators/channel/test_kraus.py
index 5d50ee9b4759..3b75b2dd614b 100644
--- a/test/python/quantum_info/operators/channel/test_kraus.py
+++ b/test/python/quantum_info/operators/channel/test_kraus.py
@@ -19,7 +19,7 @@
 
 from qiskit import QiskitError
 from qiskit.quantum_info.states import DensityMatrix
-from qiskit.quantum_info import Kraus
+from qiskit.quantum_info import Kraus, Operator
 from .channel_test_case import ChannelTestCase
 
 
@@ -68,7 +68,14 @@ def test_circuit_init(self):
         circuit, target = self.simple_circuit_no_measure()
         op = Kraus(circuit)
         target = Kraus(target)
-        self.assertEqual(op, target)
+        # The given separable circuit should only have a single Kraus operator.
+        self.assertEqual(len(op.data), 1)
+        self.assertEqual(len(target.data), 1)
+        kraus_op = Operator(op.data[0])
+        kraus_target = Operator(target.data[0])
+        # THe Kraus representation is not unique, but for a single operator, the only gauge freedom
+        # is the global phase.
+        self.assertTrue(kraus_op.equiv(kraus_target))
 
     def test_circuit_init_except(self):
         """Test initialization from circuit with measure raises exception."""
diff --git a/test/python/quantum_info/operators/channel/test_stinespring.py b/test/python/quantum_info/operators/channel/test_stinespring.py
index 9bcc886a026c..693e85d7c1cc 100644
--- a/test/python/quantum_info/operators/channel/test_stinespring.py
+++ b/test/python/quantum_info/operators/channel/test_stinespring.py
@@ -19,7 +19,7 @@
 
 from qiskit import QiskitError
 from qiskit.quantum_info.states import DensityMatrix
-from qiskit.quantum_info import Stinespring
+from qiskit.quantum_info import Stinespring, Operator
 from .channel_test_case import ChannelTestCase
 
 
@@ -61,7 +61,10 @@ def test_circuit_init(self):
         circuit, target = self.simple_circuit_no_measure()
         op = Stinespring(circuit)
         target = Stinespring(target)
-        self.assertEqual(op, target)
+        # If the Stinespring is CPTP (and it should be), it's defined in terms of a single
+        # rectangular operator, which has global-phase gauge freedom.
+        self.assertTrue(op.is_cptp())
+        self.assertTrue(Operator(op.data).equiv(Operator(target.data)))
 
     def test_circuit_init_except(self):
         """Test initialization from circuit with measure raises exception."""
diff --git a/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py b/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py
index dedd84279a8d..65f19eb8e44c 100644
--- a/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py
+++ b/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py
@@ -271,7 +271,7 @@ def test_to_matrix_zero(self):
 
         zero_sparse = zero.to_matrix(sparse=True)
         self.assertIsInstance(zero_sparse, scipy.sparse.csr_matrix)
-        np.testing.assert_array_equal(zero_sparse.A, zero_numpy)
+        np.testing.assert_array_equal(zero_sparse.todense(), zero_numpy)
 
     def test_to_matrix_parallel_vs_serial(self):
         """Parallel execution should produce the same results as serial execution up to
diff --git a/test/python/quantum_info/operators/test_operator.py b/test/python/quantum_info/operators/test_operator.py
index 725c46576a9d..d9423d0ec141 100644
--- a/test/python/quantum_info/operators/test_operator.py
+++ b/test/python/quantum_info/operators/test_operator.py
@@ -20,12 +20,14 @@
 
 from test import combine
 import numpy as np
-from ddt import ddt
+import ddt
 from numpy.testing import assert_allclose
-import scipy.linalg as la
+import scipy.stats
+import scipy.linalg
 
 from qiskit import QiskitError
 from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
+from qiskit.circuit import library
 from qiskit.circuit.library import HGate, CHGate, CXGate, QFT
 from qiskit.transpiler import CouplingMap
 from qiskit.transpiler.layout import Layout, TranspileLayout
@@ -97,7 +99,7 @@ def simple_circuit_with_measure(self):
         return circ
 
 
-@ddt
+@ddt.ddt
 class TestOperator(OperatorTestCase):
     """Tests for Operator linear operator class."""
 
@@ -290,7 +292,7 @@ def test_copy(self):
     def test_is_unitary(self):
         """Test is_unitary method."""
         # X-90 rotation
-        X90 = la.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
+        X90 = scipy.linalg.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
         self.assertTrue(Operator(X90).is_unitary())
         # Non-unitary should return false
         self.assertFalse(Operator([[1, 0], [0, 0]]).is_unitary())
@@ -495,7 +497,7 @@ def test_compose_front_subsystem(self):
 
     def test_power(self):
         """Test power method."""
-        X90 = la.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
+        X90 = scipy.linalg.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
         op = Operator(X90)
         self.assertEqual(op.power(2), Operator([[0, -1j], [-1j, 0]]))
         self.assertEqual(op.power(4), Operator(-1 * np.eye(2)))
@@ -513,6 +515,43 @@ def test_floating_point_power(self):
 
         self.assertEqual(op.power(0.25), expected_op)
 
+    def test_power_of_nonunitary(self):
+        """Test power method for a non-unitary matrix."""
+        data = [[1, 1], [0, -1]]
+        powered = Operator(data).power(0.5)
+        expected = Operator([[1.0 + 0.0j, 0.5 - 0.5j], [0.0 + 0.0j, 0.0 + 1.0j]])
+        assert_allclose(powered.data, expected.data)
+
+    @ddt.data(0.5, 1.0 / 3.0, 0.25)
+    def test_root_stability(self, root):
+        """Test that the root of operators that have eigenvalues that are -1 up to floating-point
+        imprecision stably choose the positive side of the principal-root branch cut."""
+        rng = np.random.default_rng(2024_10_22)
+
+        eigenvalues = np.array([1.0, -1.0], dtype=complex)
+        # We have the eigenvalues exactly, so this will safely find the principal root.
+        root_eigenvalues = eigenvalues**root
+
+        # Now, we can construct a bunch of Haar-random unitaries with our chosen eigenvalues.  Since
+        # we already know their eigenvalue decompositions exactly (up to floating-point precision in
+        # the matrix multiplications), we can also compute the principal values of the roots of the
+        # complete matrices.
+        bases = scipy.stats.unitary_group.rvs(2, size=16, random_state=rng)
+        matrices = [basis.conj().T @ np.diag(eigenvalues) @ basis for basis in bases]
+        expected = [basis.conj().T @ np.diag(root_eigenvalues) @ basis for basis in bases]
+        self.assertEqual(
+            [Operator(matrix).power(root) for matrix in matrices],
+            [Operator(single) for single in expected],
+        )
+
+    def test_root_branch_cut(self):
+        """Test that we can choose where the branch cut appears in the root."""
+        z_op = Operator(library.ZGate())
+        # Depending on the direction we move the branch cut, we should be able to select the root to
+        # be either of the two valid options.
+        self.assertEqual(z_op.power(0.5, branch_cut_rotation=1e-3), Operator(library.SGate()))
+        self.assertEqual(z_op.power(0.5, branch_cut_rotation=-1e-3), Operator(library.SdgGate()))
+
     def test_expand(self):
         """Test expand method."""
         mat1 = self.UX
diff --git a/test/python/transpiler/test_remove_identity_equivalent.py b/test/python/transpiler/test_remove_identity_equivalent.py
new file mode 100644
index 000000000000..1db392d3654b
--- /dev/null
+++ b/test/python/transpiler/test_remove_identity_equivalent.py
@@ -0,0 +1,185 @@
+# 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.
+
+"""Tests for the DropNegligible transpiler pass."""
+
+import numpy as np
+
+from qiskit.circuit import Parameter, QuantumCircuit, QuantumRegister, Gate
+from qiskit.circuit.library import (
+    CPhaseGate,
+    RXGate,
+    RXXGate,
+    RYGate,
+    RYYGate,
+    RZGate,
+    RZZGate,
+    XXMinusYYGate,
+    XXPlusYYGate,
+    GlobalPhaseGate,
+)
+from qiskit.quantum_info import Operator
+from qiskit.transpiler.passes import RemoveIdentityEquivalent
+from qiskit.transpiler.target import Target, InstructionProperties
+
+from test import QiskitTestCase  # pylint: disable=wrong-import-order
+
+
+class TestDropNegligible(QiskitTestCase):
+    """Test the DropNegligible pass."""
+
+    def test_drops_negligible_gates(self):
+        """Test that negligible gates are dropped."""
+        qubits = QuantumRegister(2)
+        circuit = QuantumCircuit(qubits)
+        a, b = qubits
+        circuit.append(CPhaseGate(1e-5), [a, b])
+        circuit.append(CPhaseGate(1e-8), [a, b])
+        circuit.append(RXGate(1e-5), [a])
+        circuit.append(RXGate(1e-8), [a])
+        circuit.append(RYGate(1e-5), [a])
+        circuit.append(RYGate(1e-8), [a])
+        circuit.append(RZGate(1e-5), [a])
+        circuit.append(RZGate(1e-8), [a])
+        circuit.append(RXXGate(1e-5), [a, b])
+        circuit.append(RXXGate(1e-8), [a, b])
+        circuit.append(RYYGate(1e-5), [a, b])
+        circuit.append(RYYGate(1e-8), [a, b])
+        circuit.append(RZZGate(1e-5), [a, b])
+        circuit.append(RZZGate(1e-8), [a, b])
+        circuit.append(XXPlusYYGate(1e-5, 1e-8), [a, b])
+        circuit.append(XXPlusYYGate(1e-8, 1e-8), [a, b])
+        circuit.append(XXMinusYYGate(1e-5, 1e-8), [a, b])
+        circuit.append(XXMinusYYGate(1e-8, 1e-8), [a, b])
+        transpiled = RemoveIdentityEquivalent()(circuit)
+        self.assertEqual(circuit.count_ops()["cp"], 2)
+        self.assertEqual(transpiled.count_ops()["cp"], 1)
+        self.assertEqual(circuit.count_ops()["rx"], 2)
+        self.assertEqual(transpiled.count_ops()["rx"], 1)
+        self.assertEqual(circuit.count_ops()["ry"], 2)
+        self.assertEqual(transpiled.count_ops()["ry"], 1)
+        self.assertEqual(circuit.count_ops()["rz"], 2)
+        self.assertEqual(transpiled.count_ops()["rz"], 1)
+        self.assertEqual(circuit.count_ops()["rxx"], 2)
+        self.assertEqual(transpiled.count_ops()["rxx"], 1)
+        self.assertEqual(circuit.count_ops()["ryy"], 2)
+        self.assertEqual(transpiled.count_ops()["ryy"], 1)
+        self.assertEqual(circuit.count_ops()["rzz"], 2)
+        self.assertEqual(transpiled.count_ops()["rzz"], 1)
+        self.assertEqual(circuit.count_ops()["xx_plus_yy"], 2)
+        self.assertEqual(transpiled.count_ops()["xx_plus_yy"], 1)
+        self.assertEqual(circuit.count_ops()["xx_minus_yy"], 2)
+        self.assertEqual(transpiled.count_ops()["xx_minus_yy"], 1)
+        np.testing.assert_allclose(
+            np.array(Operator(circuit)), np.array(Operator(transpiled)), atol=1e-7
+        )
+
+    def test_handles_parameters(self):
+        """Test that gates with parameters are ignored gracefully."""
+        qubits = QuantumRegister(2)
+        circuit = QuantumCircuit(qubits)
+        a, b = qubits
+        theta = Parameter("theta")
+        circuit.append(CPhaseGate(theta), [a, b])
+        circuit.append(CPhaseGate(1e-5), [a, b])
+        circuit.append(CPhaseGate(1e-8), [a, b])
+        transpiled = RemoveIdentityEquivalent()(circuit)
+        self.assertEqual(circuit.count_ops()["cp"], 3)
+        self.assertEqual(transpiled.count_ops()["cp"], 2)
+
+    def test_approximation_degree(self):
+        """Test that approximation degree handled correctly."""
+        qubits = QuantumRegister(2)
+        circuit = QuantumCircuit(qubits)
+        a, b = qubits
+        circuit.append(CPhaseGate(1e-4), [a, b])
+        # fidelity 0.9999850001249996 which is above the threshold and not excluded
+        # so 1e-2 is the only gate remaining
+        circuit.append(CPhaseGate(1e-2), [a, b])
+        circuit.append(CPhaseGate(1e-20), [a, b])
+        transpiled = RemoveIdentityEquivalent(approximation_degree=0.9999999)(circuit)
+        self.assertEqual(circuit.count_ops()["cp"], 3)
+        self.assertEqual(transpiled.count_ops()["cp"], 1)
+        self.assertEqual(transpiled.data[0].operation.params[0], 1e-2)
+
+    def test_target_approx_none(self):
+        """Test error rate with target."""
+
+        target = Target()
+        props = {(0, 1): InstructionProperties(error=1e-10)}
+        target.add_instruction(CPhaseGate(Parameter("theta")), props)
+        circuit = QuantumCircuit(2)
+        circuit.append(CPhaseGate(1e-4), [0, 1])
+        circuit.append(CPhaseGate(1e-2), [0, 1])
+        circuit.append(CPhaseGate(1e-20), [0, 1])
+        transpiled = RemoveIdentityEquivalent(approximation_degree=None, target=target)(circuit)
+        self.assertEqual(circuit.count_ops()["cp"], 3)
+        self.assertEqual(transpiled.count_ops()["cp"], 2)
+
+    def test_target_approx_approx_degree(self):
+        """Test error rate with target."""
+
+        target = Target()
+        props = {(0, 1): InstructionProperties(error=1e-10)}
+        target.add_instruction(CPhaseGate(Parameter("theta")), props)
+        circuit = QuantumCircuit(2)
+        circuit.append(CPhaseGate(1e-4), [0, 1])
+        circuit.append(CPhaseGate(1e-2), [0, 1])
+        circuit.append(CPhaseGate(1e-20), [0, 1])
+        transpiled = RemoveIdentityEquivalent(approximation_degree=0.9999999, target=target)(
+            circuit
+        )
+        self.assertEqual(circuit.count_ops()["cp"], 3)
+        self.assertEqual(transpiled.count_ops()["cp"], 2)
+
+    def test_custom_gate_no_matrix(self):
+        """Test that opaque gates are ignored."""
+
+        class CustomOpaqueGate(Gate):
+            """Custom opaque gate."""
+
+            def __init__(self):
+                super().__init__("opaque", 2, [])
+
+        qc = QuantumCircuit(3)
+        qc.append(CustomOpaqueGate(), [0, 1])
+        transpiled = RemoveIdentityEquivalent()(qc)
+        self.assertEqual(qc, transpiled)
+
+    def test_custom_gate_identity_matrix(self):
+        """Test that custom gates with matrix are evaluated."""
+
+        class CustomGate(Gate):
+            """Custom gate."""
+
+            def __init__(self):
+                super().__init__("custom", 3, [])
+
+            def to_matrix(self):
+                return np.eye(8, dtype=complex)
+
+        qc = QuantumCircuit(3)
+        qc.append(CustomGate(), [0, 1, 2])
+        transpiled = RemoveIdentityEquivalent()(qc)
+        expected = QuantumCircuit(3)
+        self.assertEqual(expected, transpiled)
+
+    def test_global_phase_ignored(self):
+        """Test that global phase gate isn't considered."""
+
+        qc = QuantumCircuit(1)
+        qc.id(0)
+        qc.append(GlobalPhaseGate(0))
+        transpiled = RemoveIdentityEquivalent()(qc)
+        expected = QuantumCircuit(1)
+        expected.append(GlobalPhaseGate(0))
+        self.assertEqual(transpiled, expected)
diff --git a/test/python/transpiler/test_target.py b/test/python/transpiler/test_target.py
index fc9eb2a6923c..980924224d15 100644
--- a/test/python/transpiler/test_target.py
+++ b/test/python/transpiler/test_target.py
@@ -11,6 +11,7 @@
 # that they have been altered from the originals.
 
 # pylint: disable=missing-docstring
+from pickle import loads, dumps
 
 import math
 import numpy as np
@@ -30,6 +31,7 @@
     CCXGate,
     RZXGate,
     CZGate,
+    UnitaryGate,
 )
 from qiskit.circuit import IfElseOp, ForLoopOp, WhileLoopOp, SwitchCaseOp
 from qiskit.circuit.measure import Measure
@@ -1166,6 +1168,23 @@ def test_instruction_supported_no_args(self):
     def test_instruction_supported_no_operation(self):
         self.assertFalse(self.ibm_target.instruction_supported(qargs=(0,), parameters=[math.pi]))
 
+    def test_target_serialization_preserve_variadic(self):
+        """Checks that variadics are still seen as variadic after serialization"""
+
+        target = Target("test", 2)
+        # Add variadic example gate with no properties.
+        target.add_instruction(UnitaryGate, None, "u_var")
+
+        # Check that this this instruction is compatible with qargs (0,). Should be
+        # true since variadic operation can be used with any valid qargs.
+        self.assertTrue(target.instruction_supported("u_var", (0, 1)))
+
+        # Rebuild the target using serialization
+        deserialized_target = loads(dumps(target))
+
+        # Perform check again, should not throw exception
+        self.assertTrue(deserialized_target.instruction_supported("u_var", (0, 1)))
+
 
 class TestPulseTarget(QiskitTestCase):
     def setUp(self):
diff --git a/test/python/utils/test_deprecation.py b/test/python/utils/test_deprecation.py
index 5d3b82e75aa2..08872b9e7a41 100644
--- a/test/python/utils/test_deprecation.py
+++ b/test/python/utils/test_deprecation.py
@@ -14,6 +14,7 @@
 
 from __future__ import annotations
 
+import sys
 from textwrap import dedent
 
 from qiskit.utils.deprecation import (
@@ -409,17 +410,21 @@ def func3():
             to a new line."""
 
         add_deprecation_to_docstring(func3, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func3.__doc__,
-            (
-                f"""Docstring extending
+        expected_doc = f"""Docstring extending
             to a new line.
 {indent}
             .. deprecated:: 9.99
               Deprecated!
 {indent}"""
-            ),
-        )
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """Docstring extending
+to a new line.
+
+.. deprecated:: 9.99
+  Deprecated!
+"""
+
+        self.assertEqual(func3.__doc__, expected_doc)
 
         def func4():
             """
@@ -427,10 +432,7 @@ def func4():
             """
 
         add_deprecation_to_docstring(func4, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func4.__doc__,
-            (
-                f"""\
+        expected_doc = f"""\
 
             Docstring starting on a new line.
 {indent}
@@ -438,7 +440,18 @@ def func4():
             .. deprecated:: 9.99
               Deprecated!
 {indent}"""
-            ),
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """\
+
+Docstring starting on a new line.
+
+.. deprecated:: 9.99
+  Deprecated!
+"""
+
+        self.assertEqual(
+            func4.__doc__,
+            expected_doc,
         )
 
         def func5():
@@ -450,11 +463,7 @@ def func5():
 
             """
 
-        add_deprecation_to_docstring(func5, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func5.__doc__,
-            (
-                f"""\
+        expected_doc = f"""\
 
             Paragraph 1, line 1.
             Line 2.
@@ -466,7 +475,24 @@ def func5():
             .. deprecated:: 9.99
               Deprecated!
 {indent}"""
-            ),
+
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """\
+
+Paragraph 1, line 1.
+Line 2.
+
+Paragraph 2.
+
+
+.. deprecated:: 9.99
+  Deprecated!
+"""
+
+        add_deprecation_to_docstring(func5, msg="Deprecated!", since="9.99", pending=False)
+        self.assertEqual(
+            func5.__doc__,
+            expected_doc,
         )
 
         def func6():
@@ -478,11 +504,7 @@ def func6():
                 continued
             """
 
-        add_deprecation_to_docstring(func6, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func6.__doc__,
-            (
-                f"""Blah.
+        expected_doc = f"""Blah.
 
             A list.
               * element 1
@@ -493,7 +515,23 @@ def func6():
             .. deprecated:: 9.99
               Deprecated!
 {indent}"""
-            ),
+
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """Blah.
+
+A list.
+  * element 1
+  * element 2
+    continued
+
+.. deprecated:: 9.99
+  Deprecated!
+"""
+
+        add_deprecation_to_docstring(func6, msg="Deprecated!", since="9.99", pending=False)
+        self.assertEqual(
+            func6.__doc__,
+            expected_doc,
         )
 
     def test_add_deprecation_docstring_meta_lines(self) -> None:
@@ -511,10 +549,7 @@ def func1():
             """
 
         add_deprecation_to_docstring(func1, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func1.__doc__,
-            (
-                f"""\
+        expected_doc = f"""\
 {indent}
             .. deprecated:: 9.99
               Deprecated!
@@ -526,7 +561,22 @@ def func1():
             Raises:
                 SomeError
 {indent}"""
-            ),
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """\
+
+.. deprecated:: 9.99
+  Deprecated!
+
+
+Returns:
+    Content.
+
+Raises:
+    SomeError"""
+
+        self.assertEqual(
+            func1.__doc__,
+            expected_doc,
         )
 
         def func2():
@@ -537,10 +587,7 @@ def func2():
             """
 
         add_deprecation_to_docstring(func2, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func2.__doc__,
-            (
-                f"""Docstring.
+        expected_doc = f"""Docstring.
 {indent}
             .. deprecated:: 9.99
               Deprecated!
@@ -549,8 +596,17 @@ def func2():
             Returns:
                 Content.
 {indent}"""
-            ),
-        )
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """Docstring.
+
+.. deprecated:: 9.99
+  Deprecated!
+
+
+Returns:
+    Content."""
+
+        self.assertEqual(func2.__doc__, expected_doc)
 
         def func3():
             """
@@ -562,11 +618,7 @@ def func3():
                 Content.
             """
 
-        add_deprecation_to_docstring(func3, msg="Deprecated!", since="9.99", pending=False)
-        self.assertEqual(
-            func3.__doc__,
-            (
-                f"""\
+        expected_doc = f"""\
 
             Docstring starting on a new line.
 
@@ -579,7 +631,24 @@ def func3():
             Examples:
                 Content.
 {indent}"""
-            ),
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """\
+
+Docstring starting on a new line.
+
+Paragraph 2.
+
+.. deprecated:: 9.99
+  Deprecated!
+
+
+Examples:
+    Content."""
+
+        add_deprecation_to_docstring(func3, msg="Deprecated!", since="9.99", pending=False)
+        self.assertEqual(
+            func3.__doc__,
+            expected_doc,
         )
 
     def test_add_deprecation_docstring_multiple_entries(self) -> None:
@@ -613,10 +682,7 @@ def func2():
 
         add_deprecation_to_docstring(func2, msg="Deprecated #1!", since="9.99", pending=False)
         add_deprecation_to_docstring(func2, msg="Deprecated #2!", since="9.99", pending=False)
-        self.assertEqual(
-            func2.__doc__,
-            (
-                f"""\
+        expected_doc = f"""\
 
             Docstring starting on a new line.
 {indent}
@@ -628,7 +694,21 @@ def func2():
             .. deprecated:: 9.99
               Deprecated #2!
 {indent}"""
-            ),
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """\
+
+Docstring starting on a new line.
+
+.. deprecated:: 9.99
+  Deprecated #1!
+
+.. deprecated:: 9.99
+  Deprecated #2!
+"""
+
+        self.assertEqual(
+            func2.__doc__,
+            expected_doc,
         )
 
         def func3():
@@ -638,12 +718,7 @@ def func3():
                 Content.
             """
 
-        add_deprecation_to_docstring(func3, msg="Deprecated #1!", since="9.99", pending=False)
-        add_deprecation_to_docstring(func3, msg="Deprecated #2!", since="9.99", pending=False)
-        self.assertEqual(
-            func3.__doc__,
-            (
-                f"""Docstring.
+        expected_doc = f"""Docstring.
 {indent}
             .. deprecated:: 9.99
               Deprecated #1!
@@ -656,7 +731,26 @@ def func3():
             Yields:
                 Content.
 {indent}"""
-            ),
+
+        if sys.version_info >= (3, 13, 0):
+            expected_doc = """Docstring.
+
+.. deprecated:: 9.99
+  Deprecated #1!
+
+
+.. deprecated:: 9.99
+  Deprecated #2!
+
+
+Yields:
+    Content."""
+
+        add_deprecation_to_docstring(func3, msg="Deprecated #1!", since="9.99", pending=False)
+        add_deprecation_to_docstring(func3, msg="Deprecated #2!", since="9.99", pending=False)
+        self.assertEqual(
+            func3.__doc__,
+            expected_doc,
         )
 
     def test_add_deprecation_docstring_pending(self) -> None:
diff --git a/tox.ini b/tox.ini
index 5456e7731920..d4d822766675 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,6 +1,6 @@
 [tox]
 minversion = 4.0
-envlist = py39, py310, py311, py312, lint-incr
+envlist = py39, py310, py311, py312, py313, lint-incr
 isolated_build = true
 
 [testenv]