fix PauliEvo for all identities

crates/accelerate/src/circuit_library/pauli_evolution.rs

@@ -192,6 +192,7 @@ fn multi_qubit_evolution(
 /// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
 /// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
 /// Pauli string XIYZ we have the circuit
+///
 ///                 ┌───┐┌───────┐┌───┐
 /// 0: ─────────────┤ X ├┤ Rz(2) ├┤ X ├───────────
 ///    ┌──────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌────┐
@@ -205,9 +206,9 @@ fn multi_qubit_evolution(
 /// Args:
 ///     num_qubits: The number of qubits in the Hamiltonian.
 ///     sparse_paulis: The Paulis to implement. Given in a sparse-list format with elements
-///         ``(pauli_string, qubit_indices, coefficient)``. An element of the form
-///         ``("IXYZ", [0,1,2,3], 0.2)``, for example, is interpreted in terms of qubit indices as
-///          I_q0 X_q1 Y_q2 Z_q3 and will use a RZ rotation angle of 0.4.
+///         ``(pauli_string, qubit_indices, rz_rotation_angle)``. An element of the form
+///         ``("XIYZ", [0,1,2,3], 2)``, for example, is interpreted in terms of qubit indices as
+///          I_q0 X_q1 Y_q2 Z_q3 and will use a RZ rotation angle of 2.
 ///     insert_barriers: If ``true``, insert a barrier in between the evolution of individual
 ///         Pauli terms.
 ///     do_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
@@ -244,7 +245,7 @@ pub fn py_pauli_evolution(
         }
 
         paulis.push(pauli);
-        times.push(time); // note we do not multiply by 2 here, this is done Python side!
+        times.push(time); // note we do not multiply by 2 here, this is already done Python side!
         indices.push(tuple.get_item(1)?.extract::<Vec<u32>>()?)
     }
 
@@ -266,12 +267,12 @@ pub fn py_pauli_evolution(
         },
     );
 
-    // When handling all-identity Paulis above, we added the time as global phase.
-    // However, the all-identity Paulis should add a negative phase, as they implement
-    // exp(-i t I). We apply the negative sign here, to only do a single (-1) multiplication,
-    // instead of doing it every time we find an all-identity Pauli.
+    // When handling all-identity Paulis above, we added the RZ rotation angle as global phase,
+    // meaning that we have implemented of exp(i 2t I). However, what we want it to implement
+    // exp(-i t I). To only use a single multiplication, we apply a factor of -0.5 here.
+    // This is faster, in particular as long as the parameter expressions are in Python.
     if modified_phase {
-        global_phase = multiply_param(&global_phase, -1.0, py);
+        global_phase = multiply_param(&global_phase, -0.5, py);
     }
 
     CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)

crates/circuit/src/operations.rs

@@ -2340,8 +2340,14 @@ pub fn add_param(param: &Param, summand: f64, py: Python) -> Param {
 }
 
 pub fn radd_param(param1: Param, param2: Param, py: Python) -> Param {
-    match [param1, param2] {
+    match [&param1, &param2] {
         [Param::Float(theta), Param::Float(lambda)] => Param::Float(theta + lambda),
+        [Param::Float(theta), Param::ParameterExpression(_lambda)] => {
+            add_param(&param2, *theta, py)
+        }
+        [Param::ParameterExpression(_theta), Param::Float(lambda)] => {
+            add_param(&param1, *lambda, py)
+        }
         [Param::ParameterExpression(theta), Param::ParameterExpression(lambda)] => {
             Param::ParameterExpression(
                 theta

qiskit/synthesis/evolution/suzuki_trotter.py

@@ -138,7 +138,7 @@ def expand(
 
         .. code-block:: text
 
-            ("X", [0], t), ("ZZ", [0, 1], 2t), ("X", [0], 2)
+            ("X", [0], t), ("ZZ", [0, 1], 2t), ("X", [0], t)
 
         Note that the rotation angle contains a factor of 2, such that that evolution
         of a Pauli :math:`P` over time :math:`t`, which is :math:`e^{itP}`, is represented

releasenotes/notes/fix-pauli-evo-all-identity-b129acd854d8c391.yaml

@@ -0,0 +1,8 @@
+---
+fixes:
+  - |
+    The :class:`.PauliEvolutionGate`, if used with a product formula synthesis (this is the default),
+    did not correctly handle all-identity terms in the operator. The all-identity term
+    should introduce a global phase equal to ``-evolution_time``, but was off by a factor of 2
+    and could break for parameterized times. This behavior is now fixed.
+    Fixed `#13625 <https://github.com/Qiskit/qiskit/issues/13625>`__.

test/python/circuit/library/test_evolution_gate.py

@@ -479,6 +479,36 @@ def atomic_evolution(pauli, time):
         decomposed = evo_gate.definition.decompose()
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)
 
+    def test_all_identity(self):
+        """Test circuit with all identity Paulis works correctly."""
+        evo = PauliEvolutionGate(I ^ I, time=1).definition
+        expected = QuantumCircuit(2, global_phase=-1)
+        self.assertEqual(expected, evo)
+
+    def test_global_phase(self):
+        """Test a circuit with parameterized global phase terms.
+
+        Regression test of #13625.
+        """
+        pauli = (X ^ X) + (I ^ I) + (I ^ X)
+        time = Parameter("t")
+        evo = PauliEvolutionGate(pauli, time=time)
+
+        expected = QuantumCircuit(2, global_phase=-time)
+        expected.rxx(2 * time, 0, 1)
+        expected.rx(2 * time, 0)
+
+        with self.subTest(msg="check circuit"):
+            self.assertEqual(expected, evo.definition)
+
+        # since all terms in the Pauli operator commute, we can compare to an
+        # exact matrix exponential
+        time_value = 1.76123
+        bound = evo.definition.assign_parameters([time_value])
+        exact = scipy.linalg.expm(-1j * time_value * pauli.to_matrix())
+        with self.subTest(msg="check correctness"):
+            self.assertEqual(Operator(exact), Operator(bound))
+
 
 def exact_atomic_evolution(circuit, pauli, time):
     """An exact atomic evolution for Suzuki-Trotter.