Expose branch_cut_rotation parameter in Operator.power

The rotation used to stabilise matrix roots has an impact on which
matrix is selected as the principal root.  Exposing it to users to allow
control makes the most sense.

qiskit/quantum_info/operators/operator.py

@@ -541,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``.
@@ -562,18 +588,10 @@ def power(self, n: float) -> Operator:
         else:
             import scipy.linalg
 
-            # 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
+            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
 

releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml

@@ -13,6 +13,13 @@ fixes:
     :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.
+    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.

test/python/quantum_info/operators/test_operator.py

@@ -27,6 +27,7 @@
 
 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
@@ -543,6 +544,14 @@ def test_root_stability(self, root):
             [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