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/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..79e665e1b38c 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
@@ -561,13 +562,19 @@ 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
+            # Non-integer powers of operators with an eigenvalue of -1 have a branch cut in the
+            # complex plane, which makes the calculation of the principal root around this cut
+            # finnicky and subject to precision / differences in BLAS implementation.  For example,
+            # the square root of Pauli Y can return the pi/2 or -pi/2 Y rotation depending on
+            # whether the -1 eigenvalue is found as `complex(-1, tiny)` or `complex(-1, -tiny)`.
+            # Such -1 eigenvalues are really common, so we first phase the matrix slightly to move
+            # common eigenvalues away from the branch-cut point of the power function.  The exact
+            # value of the epsilon doesn't matter much, but shouldn't coincide with any "nice"
+            # eigenvalues we expect to encounter. This isn't perfect, just pragmatic.
+            eps = cmath.pi * 1e-3
+            ret._data = cmath.rect(1, eps * n) * scipy.linalg.fractional_matrix_power(
+                cmath.rect(1, -eps) * self.data, n
+            )
         return ret
 
     def tensor(self, other: Operator) -> Operator:
diff --git a/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml
new file mode 100644
index 000000000000..939b8fe8dee6
--- /dev/null
+++ b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml
@@ -0,0 +1,18 @@
+---
+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.
+
+    See `#13305 <https://github.com/Qiskit/qiskit/issues/13305>`__.
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..f206da551169 100644
--- a/test/python/quantum_info/operators/test_operator.py
+++ b/test/python/quantum_info/operators/test_operator.py
@@ -20,9 +20,10 @@
 
 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
@@ -97,7 +98,7 @@ def simple_circuit_with_measure(self):
         return circ
 
 
-@ddt
+@ddt.ddt
 class TestOperator(OperatorTestCase):
     """Tests for Operator linear operator class."""
 
@@ -290,7 +291,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 +496,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 +514,35 @@ 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_expand(self):
         """Test expand method."""
         mat1 = self.UX