Expectation value of operator - python

Question

I am suposed to solve following problem:

Calculate the probability $P_n(x > X)$ that a particle in then n-th eigenstate is found at a position with an x-value larger than X. Here it is convenient to formulate it as an expectation value of an operator. We make an operator $P(x > X)$ which is a matrix with unity on the diagonal if $x > X$.

To do this I used the command > which returns True if the number on the left is larger than the one on the right and False otherwise. To turn the True or False into 0 or 1 we use the method astype(int) on the NumPy array, which we then transform to an operator.T ake a few values of X, e.g. X= 3, 5, and 7, and plot $P_n(x > X)$ as a function of $n$ (make sure that the box size $L$ is somewhat larger than $X$).

The excercise continues as follow

Explain what you see: Why is it low when $n$ is small? At which value ofnwould you expect $P_n(x > X)$ to grow from a small to a large value? Hint :which energy would it require classically?

My question is: How do I find the expectation value of the operator?

Answer 1

The expectation of an observable $A$ with respect to the state $|\psi \rangle$ can be calculated as: $$\langle \psi |A| \psi \rangle = Tr( \langle \psi | A| \psi \rangle) = Tr( A|\psi \rangle \langle \psi|) = Tr(A \rho) = Tr(\rho A)$$

$\rho = | \psi \rangle \langle \psi |$ is the density matrix formulation of $|\psi\rangle$.

Now, if giving a $|\psi\rangle$ and $A$, and you want to calculate $\langle \psi | A| \psi \rangle$ in Python, you can do it as follows:

import numpy as np 
norm_psi = [1., 1., 1., 1.]/np.linalg.norm(psi)
A  = np.matrix( '3,0,0,1; 0,-1,0,0; 0,0,-1,0; 1,0,0,1' )
print('Operator A=\n', A)
expectation  = np.inner(np.conj(norm_psi).T, np.matmul(A,norm_psi) ) #Calculate <norm_psi|A|norm_psi> 
print('\n expectation value calculated by <norm_psi|A|norm_psi> :\n ' , expectation)
#-------
rho =  np.outer(norm_psi, np.conj(norm_psi) )
expectation = np.trace(rho*A)
print('\n expectation value calculated by Tr(rho*A):\n ' , expectation)

The output would be something like:

Operator A=
 [[ 3  0  0  1]
 [ 0 -1  0  0]
 [ 0  0 -1  0]
 [ 1  0  0  1]]

 expectation value calculated by <norm_psi|A|norm_psi> :
  [[1.]]

 expectation value calculated by Tr(rho*A):
  1.0

---

You can use Qiskit to do this as well. And it is quite simply if your operator $A$ is straight forward decomposition of Pauli strings, like $A = Z\otimes Z$ or $A = X \otimes Z$ etc. For instance, if $$A = X \otimes Z = \begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} \otimes \begin{pmatrix} 1 & 0\\ 0 & -1 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 & 0\\ 0 & 0& 0 & -1 \\ 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \end{pmatrix}$$

and $\psi$ is the same as previously: $|\psi \rangle = \dfrac{1}{2}\begin{pmatrix} 1 \\ 1 \\ 1\\ 1 \end{pmatrix}$. Note that this state can be prepared on a quantum circuit as:

then you can calculate the expectation as:

from qiskit import QuantumCircuit
from qiskit.aqua.operators import X, Y, Z, I
from qiskit.aqua.operators import StateFn

operator = X ^ Z  #Note that ^ represents tensor product 
qc = QuantumCircuit(2)
qc.h(0)
qc.h(1)
psi = StateFn(qc) 
expectation = (~psi @ operator @ psi).eval()
print('expectation value is:', expectation)

output: expectation value is: 0j

Which is what you would expect by looking at the matrix form of $A$.