Repository: Qiskit/qiskit
Issue category: Bug Report

Issue:
Statevector._expectation_value_pauli returns incorrect result in certain cases
The calculation for the expectation value of the __identity operator__ seems incorrect when the statevector is __not normalized__ to 1.

I think the reason is that one should use the norm squared rather than the norm in the following lines.
https://github.com/Qiskit/qiskit/blob/1be2c5aa103f1991f2fdc75a77a09f82635b8431/qiskit/quantum_info/states/statevector.py#L478-L479

---

The whole function is attached below to provide context.

https://github.com/Qiskit/qiskit/blob/1be2c5aa103f1991f2fdc75a77a09f82635b8431/qiskit/quantum_info/states/statevector.py#L458-L491


Hints:
Please can you show a reproducing code block following the bug report template, and the output you expected?
Sure. Here it is

```python
from qiskit.quantum_info import Statevector, SparsePauliOp, Pauli

state_normed = Statevector.from_label("0")
state_not_normed = 2 * state_normed

identity_op = SparsePauliOp.from_list([("I", 1)])

# using the public interface (expectation_value)
print(state_not_normed.expectation_value(identity_op))

# using _expectation_value_pauli (the backend)
print(state_not_normed._expectation_value_pauli(Pauli('I')))
```

## Environment
Python 3.11.6
qiskit 0.45.3

Judging from the source code on github, the issue should be platform-independent and exists across different versions.

## What is the expected result
The two print statement should give `4+0j` and `4.0` respectively.

## What is the actual result
The two print statement gave `2+0j` and `2.0` respectively.
Nice catch, that seems to be unique only to the Pauli-I operator. Would you like to open a PR to fix the issue? 🙂 
Sure, I have created one PR.

Patch:
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index 901ce95af424..1265d677abe4 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -476,7 +476,7 @@ def _expectation_value_pauli(self, pauli, qargs=None):
         pauli_phase = (-1j) ** pauli.phase if pauli.phase else 1
 
         if x_mask + z_mask == 0:
-            return pauli_phase * np.linalg.norm(self.data)
+            return pauli_phase * np.linalg.norm(self.data) ** 2
 
         if x_mask == 0:
             return pauli_phase * expval_pauli_no_x(self.data, self.num_qubits, z_mask)


Test patch:
diff --git a/releasenotes/notes/fix_identity_operator_9e2ec9770ac046a6.yaml b/releasenotes/notes/fix_identity_operator_9e2ec9770ac046a6.yaml
new file mode 100644
index 000000000000..69779a793cc3
--- /dev/null
+++ b/releasenotes/notes/fix_identity_operator_9e2ec9770ac046a6.yaml
@@ -0,0 +1,5 @@
+---
+fixes:
+  - |
+    Fixed a bug that caused :meth:`.Statevector.expectation_value` to yield incorrect results
+    for the identity operator when the statevector was not normalized.
diff --git a/test/python/quantum_info/states/test_statevector.py b/test/python/quantum_info/states/test_statevector.py
index 29d5d42f3783..d16b3b3453ec 100644
--- a/test/python/quantum_info/states/test_statevector.py
+++ b/test/python/quantum_info/states/test_statevector.py
@@ -1152,6 +1152,31 @@ def test_expval_pauli_qargs(self, qubits):
         expval = state.expectation_value(op, qubits)
         self.assertAlmostEqual(expval, target)
 
+    def test_expval_identity(self):
+        """Test whether the calculation for identity operator has been fixed"""
+
+        # 1 qubit case test
+        state_1 = Statevector.from_label("0")
+        state_1_n1 = 2 * state_1  # test the same state with different norms
+        state_1_n2 = (1 + 2j) * state_1
+        identity_op_1 = SparsePauliOp.from_list([("I", 1)])
+        expval_state_1 = state_1.expectation_value(identity_op_1)
+        expval_state_1_n1 = state_1_n1.expectation_value(identity_op_1)
+        expval_state_1_n2 = state_1_n2.expectation_value(identity_op_1)
+        self.assertAlmostEqual(expval_state_1, 1.0 + 0j)
+        self.assertAlmostEqual(expval_state_1_n1, 4 + 0j)
+        self.assertAlmostEqual(expval_state_1_n2, 5 + 0j)
+
+        # Let's try a multi-qubit case
+        n_qubits = 3
+        state_coeff = 3 - 4j
+        op_coeff = 2 - 2j
+        state_test = state_coeff * Statevector.from_label("0" * n_qubits)
+        op_test = SparsePauliOp.from_list([("I" * n_qubits, op_coeff)])
+        expval = state_test.expectation_value(op_test)
+        target = op_coeff * np.abs(state_coeff) ** 2
+        self.assertAlmostEqual(expval, target)
+
     @data(*(qargs for i in range(4) for qargs in permutations(range(4), r=i + 1)))
     def test_probabilities_qargs(self, qargs):
         """Test probabilities method with qargs"""