Repository: Qiskit/qiskit-terra (python)

Issue:
`PauliList.evolve` raises exception when evolving certain Clifford circuits.
### Environment

- **Qiskit Terra version**:  0.21.1
- **Python version**: 3.9
- **Operating system**: MacOS 11.6.7


### What is happening?

There is a bug in the `BasePauli/PauliList/Pauli` classes where the phase vector is sometimes set as dtype `uint8` instead of `int64` leading to errors when calling the `evolve` function for circuits containing gates that modify the pauli phase vector.

Whether this happens or not seems to depend on how the Pauli is initialized (sometimes it will be initialized with uint8 phase vector, other times with int64 phase vector).

### How can we reproduce the issue?

Minimal example

```python
from qiskit import QuantumCircuit
import qiskit.quantum_info as qi

# Create circuit with gate that modifies Pauli phase
qc = QuantumCircuit(2)
qc.cz(0, 1)

# Generate random Pauli
pauli = qi.random_pauli(2, seed=123)

# Try and evolve
pauli.evolve(qc)
```

Give the exception:
```
---------------------------------------------------------------------------
UFuncTypeError                            Traceback (most recent call last)
/var/folders/lp/ftt3j6qs7s90j6h9wszq3jcw0000gn/T/ipykernel_71613/3628800509.py in <module>
      7 pauli = qi.random_pauli(2, seed=123)
      8 #qc = qi.random_clifford(4, seed=123).to_circuit()
----> 9 pauli.evolve(qc)

~/anaconda3/envs/qiskit/lib/python3.9/site-packages/qiskit/quantum_info/operators/symplectic/pauli.py in evolve(self, other, qargs, frame)
    597             other = Pauli(other)
    598 
--> 599         return Pauli(super().evolve(other, qargs=qargs, frame=frame))
    600 
    601     # ---------------------------------------------------------------------

~/anaconda3/envs/qiskit/lib/python3.9/site-packages/qiskit/quantum_info/operators/symplectic/base_pauli.py in evolve(self, other, qargs, frame)
    277         if frame == "s":
    278             return self.copy()._append_circuit(other, qargs=qargs)
--> 279         return self.copy()._append_circuit(other.inverse(), qargs=qargs)
    280 
    281     def _evolve_clifford(self, other, qargs=None, frame="h"):

~/anaconda3/envs/qiskit/lib/python3.9/site-packages/qiskit/quantum_info/operators/symplectic/base_pauli.py in _append_circuit(self, circuit, qargs)
    585             # Get the integer position of the flat register
    586             new_qubits = [qargs[bit_indices[tup]] for tup in instruction.qubits]
--> 587             self._append_circuit(instruction.operation, new_qubits)
    588 
    589         # Since the individual gate evolution functions don't take mod

~/anaconda3/envs/qiskit/lib/python3.9/site-packages/qiskit/quantum_info/operators/symplectic/base_pauli.py in _append_circuit(self, circuit, qargs)
    558             if len(qargs) != 2:
    559                 raise QiskitError("Invalid qubits for 2-qubit gate.")
--> 560             return basis_2q[name](self, qargs[0], qargs[1])
    561 
    562         # If not a Clifford basis gate we try to unroll the gate and

~/anaconda3/envs/qiskit/lib/python3.9/site-packages/qiskit/quantum_info/operators/symplectic/base_pauli.py in _evolve_cz(base_pauli, q1, q2)
    661     base_pauli._z[:, q1] ^= x2
    662     base_pauli._z[:, q2] ^= x1
--> 663     base_pauli._phase += 2 * np.logical_and(x1, x2).T
    664     return base_pauli
    665 

UFuncTypeError: Cannot cast ufunc 'add' output from dtype('int64') to dtype('uint8') with casting rule 'same_kind'
```

### What should happen?

The dtype of the phase vector should always be the same, so that evolution and other methods that attempt to inplace modify the vector don't cause exceptions for mismatched dtypes.

### Any suggestions?

One fix would be to explicitly cast the phase vector when initializing base pauli

Pull request description:
Fix pauli.evolve dtype casting error
<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add a release note file using reno if the change needs to be
  documented in the release notes.
⚠️ If your pull request fixes an open issue, please link to the issue.

- [x] I have added the tests to cover my changes.
- [x] I have updated the documentation accordingly.
- [x] I have read the CONTRIBUTING document.
-->

### Summary

Fixes #8438 

Fixes issue where `Pauli.evolve` and `PauliList.evolve` would raise a dtype error when evolving by certain Clifford gates which modified the Pauli's phase.

### Details and comments

This fix make it so that rather than cast the dtype of the phase vector to a fixed type, it will accept any signed or unsigned integer dtype and work with that.

Patch:
diff --git a/qiskit/quantum_info/operators/symplectic/base_pauli.py b/qiskit/quantum_info/operators/symplectic/base_pauli.py
index 3e03838cd..1f9e74036 100644
--- a/qiskit/quantum_info/operators/symplectic/base_pauli.py
+++ b/qiskit/quantum_info/operators/symplectic/base_pauli.py
@@ -131,9 +131,9 @@ class BasePauli(BaseOperator, AdjointMixin, MultiplyMixin):
         # Get phase shift
         phase = self._phase + other._phase
         if front:
-            phase += 2 * _count_y(x1, z2)
+            phase += 2 * _count_y(x1, z2, dtype=phase.dtype)
         else:
-            phase += 2 * _count_y(x2, z1)
+            phase += 2 * _count_y(x2, z1, dtype=phase.dtype)
 
         # Update Pauli
         x = np.logical_xor(x1, x2)
@@ -185,7 +185,7 @@ class BasePauli(BaseOperator, AdjointMixin, MultiplyMixin):
     def transpose(self):
         """Return the transpose of each Pauli in the list."""
         # Transpose sets Y -> -Y. This has effect on changing the phase
-        parity_y = self._count_y() % 2
+        parity_y = self._count_y(dtype=self._phase.dtype) % 2
         if np.all(parity_y == 0):
             return self
         return BasePauli(self._z, self._x, np.mod(self._phase + 2 * parity_y, 4))
@@ -388,7 +388,8 @@ class BasePauli(BaseOperator, AdjointMixin, MultiplyMixin):
             raise QiskitError("z and x vectors are different size.")
 
         # Convert group phase convention to internal ZX-phase conversion.
-        base_phase = np.mod(_count_y(base_x, base_z) + phase, 4)
+        dtype = getattr(phase, "dtype", None)
+        base_phase = np.mod(_count_y(base_x, base_z, dtype=dtype) + phase, 4)
         return base_z, base_x, base_phase
 
     @staticmethod
@@ -601,7 +602,7 @@ def _evolve_h(base_pauli, qubit):
     z = base_pauli._z[:, qubit].copy()
     base_pauli._x[:, qubit] = z
     base_pauli._z[:, qubit] = x
-    base_pauli._phase += 2 * np.logical_and(x, z).T
+    base_pauli._phase += 2 * np.logical_and(x, z).T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
@@ -609,7 +610,7 @@ def _evolve_s(base_pauli, qubit):
     """Update P -> S.P.Sdg"""
     x = base_pauli._x[:, qubit]
     base_pauli._z[:, qubit] ^= x
-    base_pauli._phase += x.T
+    base_pauli._phase += x.T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
@@ -617,7 +618,7 @@ def _evolve_sdg(base_pauli, qubit):
     """Update P -> Sdg.P.S"""
     x = base_pauli._x[:, qubit]
     base_pauli._z[:, qubit] ^= x
-    base_pauli._phase -= x.T
+    base_pauli._phase -= x.T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
@@ -629,19 +630,21 @@ def _evolve_i(base_pauli, qubit):
 
 def _evolve_x(base_pauli, qubit):
     """Update P -> X.P.X"""
-    base_pauli._phase += 2 * base_pauli._z[:, qubit].T
+    base_pauli._phase += 2 * base_pauli._z[:, qubit].T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
 def _evolve_y(base_pauli, qubit):
     """Update P -> Y.P.Y"""
-    base_pauli._phase += 2 * base_pauli._x[:, qubit].T + 2 * base_pauli._z[:, qubit].T
+    xp = base_pauli._x[:, qubit].T.astype(base_pauli._phase.dtype)
+    zp = base_pauli._z[:, qubit].T.astype(base_pauli._phase.dtype)
+    base_pauli._phase += 2 * (xp + zp)
     return base_pauli
 
 
 def _evolve_z(base_pauli, qubit):
     """Update P -> Z.P.Z"""
-    base_pauli._phase += 2 * base_pauli._x[:, qubit].T
+    base_pauli._phase += 2 * base_pauli._x[:, qubit].T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
@@ -658,7 +661,7 @@ def _evolve_cz(base_pauli, q1, q2):
     x2 = base_pauli._x[:, q2].copy()
     base_pauli._z[:, q1] ^= x2
     base_pauli._z[:, q2] ^= x1
-    base_pauli._phase += 2 * np.logical_and(x1, x2).T
+    base_pauli._phase += 2 * np.logical_and(x1, x2).T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
@@ -670,7 +673,7 @@ def _evolve_cy(base_pauli, qctrl, qtrgt):
     base_pauli._x[:, qtrgt] ^= x1
     base_pauli._z[:, qtrgt] ^= x1
     base_pauli._z[:, qctrl] ^= np.logical_xor(x2, z2)
-    base_pauli._phase += x1 + 2 * np.logical_and(x1, x2).T
+    base_pauli._phase += x1 + 2 * np.logical_and(x1, x2).T.astype(base_pauli._phase.dtype)
     return base_pauli
 
 
@@ -687,7 +690,4 @@ def _evolve_swap(base_pauli, q1, q2):
 
 def _count_y(x, z, dtype=None):
     """Count the number of I Pauli's"""
-    axis = 1
-    if dtype is None:
-        dtype = np.min_scalar_type(x.shape[axis])
-    return (x & z).sum(axis=axis, dtype=dtype)
+    return (x & z).sum(axis=1, dtype=dtype)
diff --git a/qiskit/quantum_info/operators/symplectic/pauli.py b/qiskit/quantum_info/operators/symplectic/pauli.py
index 941ae4d52..b2663222d 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli.py
@@ -286,12 +286,12 @@ class Pauli(BasePauli):
     def phase(self):
         """Return the group phase exponent for the Pauli."""
         # Convert internal ZX-phase convention of BasePauli to group phase
-        return np.mod(self._phase - self._count_y(), 4)[0]
+        return np.mod(self._phase - self._count_y(dtype=self._phase.dtype), 4)[0]
 
     @phase.setter
     def phase(self, value):
         # Convert group phase convention to internal ZX-phase convention
-        self._phase[:] = np.mod(value + self._count_y(), 4)
+        self._phase[:] = np.mod(value + self._count_y(dtype=self._phase.dtype), 4)
 
     @property
     def x(self):
@@ -639,7 +639,9 @@ class Pauli(BasePauli):
             raise QiskitError(f"{op} is not an N-qubit identity")
         base_z = np.zeros((1, op.num_qubits), dtype=bool)
         base_x = np.zeros((1, op.num_qubits), dtype=bool)
-        base_phase = np.mod(cls._phase_from_complex(op.coeff) + _count_y(base_x, base_z), 4)
+        base_phase = np.mod(
+            cls._phase_from_complex(op.coeff) + _count_y(base_x, base_z), 4, dtype=int
+        )
         return base_z, base_x, base_phase
 
     @classmethod
diff --git a/qiskit/quantum_info/operators/symplectic/pauli_list.py b/qiskit/quantum_info/operators/symplectic/pauli_list.py
index 188711e4b..293e0a641 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli_list.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli_list.py
@@ -237,12 +237,12 @@ class PauliList(BasePauli, LinearMixin, GroupMixin):
     def phase(self):
         """Return the phase exponent of the PauliList."""
         # Convert internal ZX-phase convention to group phase convention
-        return np.mod(self._phase - self._count_y(), 4)
+        return np.mod(self._phase - self._count_y(dtype=self._phase.dtype), 4)
 
     @phase.setter
     def phase(self, value):
         # Convert group phase convetion to internal ZX-phase convention
-        self._phase[:] = np.mod(value + self._count_y(), 4)
+        self._phase[:] = np.mod(value + self._count_y(dtype=self._phase.dtype), 4)
 
     @property
     def x(self):
diff --git a/releasenotes/notes/fix_8438-159e67ecb6765d08.yaml b/releasenotes/notes/fix_8438-159e67ecb6765d08.yaml
new file mode 100644
index 000000000..189827703
--- /dev/null
+++ b/releasenotes/notes/fix_8438-159e67ecb6765d08.yaml
@@ -0,0 +1,7 @@
+---
+fixes:
+  - |
+    Fixes issue where :meth:`~.Pauli.evolve` and `~.PauliList.evolve` would
+    raise a dtype error when evolving by certain Clifford gates which
+    modified the Pauli's phase.
+    Fixes `issue #8438 <https://github.com/Qiskit/qiskit-terra/issues/8438>`_


Test patch:
diff --git a/test/python/quantum_info/operators/symplectic/test_pauli.py b/test/python/quantum_info/operators/symplectic/test_pauli.py
index c952a3a8f..3f529924b 100644
--- a/test/python/quantum_info/operators/symplectic/test_pauli.py
+++ b/test/python/quantum_info/operators/symplectic/test_pauli.py
@@ -404,6 +404,34 @@ class TestPauli(QiskitTestCase):
         self.assertEqual(value, value_h)
         self.assertEqual(value_inv, value_s)
 
+    @data(
+        *it.product(
+            (
+                IGate(),
+                XGate(),
+                YGate(),
+                ZGate(),
+                HGate(),
+                SGate(),
+                SdgGate(),
+                CXGate(),
+                CYGate(),
+                CZGate(),
+                SwapGate(),
+            ),
+            [int, np.int8, np.uint8, np.int16, np.uint16, np.int32, np.uint32, np.int64, np.uint64],
+        )
+    )
+    @unpack
+    def test_phase_dtype_evolve_clifford(self, gate, dtype):
+        """Test phase dtype for evolve method for Clifford gates."""
+        z = np.ones(gate.num_qubits, dtype=bool)
+        x = np.ones(gate.num_qubits, dtype=bool)
+        phase = (np.sum(z & x) % 4).astype(dtype)
+        paulis = Pauli((z, x, phase))
+        evo = paulis.evolve(gate)
+        self.assertEqual(evo.phase.dtype, dtype)
+
     def test_evolve_clifford_qargs(self):
         """Test evolve method for random Clifford"""
         cliff = random_clifford(3, seed=10)
diff --git a/test/python/quantum_info/operators/symplectic/test_pauli_list.py b/test/python/quantum_info/operators/symplectic/test_pauli_list.py
index 70b4bc267..db38d428a 100644
--- a/test/python/quantum_info/operators/symplectic/test_pauli_list.py
+++ b/test/python/quantum_info/operators/symplectic/test_pauli_list.py
@@ -2034,6 +2034,40 @@ class TestPauliListMethods(QiskitTestCase):
         self.assertListEqual(value, value_h)
         self.assertListEqual(value_inv, value_s)
 
+    def test_phase_dtype_evolve_clifford(self):
+        """Test phase dtype during evolve method for Clifford gates."""
+        gates = (
+            IGate(),
+            XGate(),
+            YGate(),
+            ZGate(),
+            HGate(),
+            SGate(),
+            SdgGate(),
+            CXGate(),
+            CYGate(),
+            CZGate(),
+            SwapGate(),
+        )
+        dtypes = [
+            int,
+            np.int8,
+            np.uint8,
+            np.int16,
+            np.uint16,
+            np.int32,
+            np.uint32,
+            np.int64,
+            np.uint64,
+        ]
+        for gate, dtype in itertools.product(gates, dtypes):
+            z = np.ones(gate.num_qubits, dtype=bool)
+            x = np.ones(gate.num_qubits, dtype=bool)
+            phase = (np.sum(z & x) % 4).astype(dtype)
+            paulis = Pauli((z, x, phase))
+            evo = paulis.evolve(gate)
+            self.assertEqual(evo.phase.dtype, dtype)
+
     @combine(phase=(True, False))
     def test_evolve_clifford_qargs(self, phase):
         """Test evolve method for random Clifford"""