We'll use $\hbar=1$ for numerical results, this is standard practice, but it is important to remember.
from numpy import sqrt, pi, exp, angle
from qutip import *
Our usual spin operators (spin-1/2)
pz = basis(2,0)
mz = basis(2,1)
px = 1/sqrt(2)*(pz + mz)
mx = 1/sqrt(2)*(pz - mz)
py = 1/sqrt(2)*(pz + 1j*mz)
my = 1/sqrt(2)*(pz - 1j*mz)
Sx = 1/2.0*sigmax()
Sy = 1/2.0*sigmay()
Sz = 1/2.0*sigmaz()
Splus = Sx + 1j*Sy
Sminus = Sx - 1j*Sy
Define the spin-1 operators. We use J just to keep them apart. Could be S instead.
Jx = 1/sqrt(2)*Qobj([[0,1,0],[1,0,1],[0,1,0]])
Jy = 1/sqrt(2)*Qobj([[0,-1j,0],[1j,0,-1j],[0,1j,0]])
Jz = Qobj([[1,0,0],[0,0,0],[0,0,-1]])
Jplus = Jx + 1j*Jy # raising operator
Jminus = Jx - 1j*Jy # lowering operator
Jz
To define an exponentiated operator, we apply the .expm()
method to the argument. So to write the rotation operator:
we do the following. It looks a little odd to read since you might expect the exp
to come first, but looks like this:
Rz90 = (-1j*pi/2*Sz).expm()
Rz90
Rz90*px
Doesn't look like example 7.4 because it's not simplified. How to check:
exp(-1j*pi/4)*py
Yes, this agrees.
Can also use inner product:
py.dag()*Rz90*px
(py.dag()*Rz90*px).norm()
0.9999999999999999
angle(0.707 - 0.707j) == -pi/4
True
According to the postulates, the allowed values of an observable are the eigenvalues with corresponding eigenstates.
We know the matrix representation of Jx, Jy, Jz in the Z-basis so we can find all 9 states (in the Z basis). Why nine? There are three possible values $\hbar$, 0, $-\hbar$ for each of three directions.
yevals, (yd,y0,yu) = Jy.eigenstates()
zevals, (zd,z0,zu) = Jz.eigenstates()
xevals, (xd,x0,xu) = Jx.eigenstates()
# Fix the signs to match book:
# (if first value of eigenstate is negative, flip all signs)
states = [xd, x0, xu, yd, y0, yu, zd, z0, zu]
states = [-s if s[0][0][0]<0 else s for s in states]
# Python doesn't change the values in the list so we have
# to replace the old values with the new ones:
[xd, x0, xu, yd, y0, yu, zd, z0, zu] = states
xd
Rz90 = (-1j*pi/2*Jz).expm()
Rz90*xu == -1j*yu
True
Rz90*xu
yd
(y0.dag()*Rz90*x0).norm()
1.0
(z0.dag()*yu).norm()**2
0.4999999999999999
(zu.dag()*yu).norm()**2
0.2500000000000002
(zd.dag()*yu).norm()**2
0.2500000000000001
(y0.dag()*xd).norm()**2
0.4999999999999998
(xd.dag()*y0).norm()**2
0.4999999999999998
psi = (1/3)*(2*zu -1j*z0 +2*zd)
psi
Jsquared = Jx*Jx + Jy*Jy + Jz*Jz
Jsquared
psi.dag()*Jsquared*psi
psi.dag()*Jy*psi
psi.dag()*Jz*psi