Quantum Mechanics with Python Coding Solutions

Perform the following calculations Python or another programming language. You may use a linear algebra library for your work. Compare your results to those calculated on the conceptual component of the homework.

Part A

Consider the following Hamiltonian (where \(\hbar\) and \(\omega\) are positive, real constants):

\[\hat{H} = \hbar\omega\begin{bmatrix} \frac{1}{2} & 0 & 0 \\ 0 & \frac{3}{2} & 0 \\ 0 & 0 & \frac{5}{2}\end{bmatrix},\]

which describes the energy of the following wavefunction:

\[|\psi\rangle = i|1\rangle + 5|2\rangle - 2i|3\rangle,\]

where \(|1\rangle\) corresponds to the eigenvector of \(\hat{H}\) with the lowest eigenvalue, \(|2\rangle\) to the eigenvector with the middle eigenvalue, and \(|3\rangle\) to the eigenvector with the highest eigenvalue.

Is \(|\psi\rangle\) normalized? If so, prove it. If not, normalize it.
If you measure the energy of the system, what are the possible results and with what probabilities will they occur?
Find \(\langle \psi | \hat{H} | \psi \rangle = \langle E \rangle\).
Show that \(|1\rangle\), \(|2\rangle\), and \(|3\rangle\) have the properties of a computational basis.
Define an operator which could be applied to \(|\psi\rangle\). It must have the properties of an operator but it does not need to correspond to an actual operator. What are the possible observables with this new operator and with what probabilities will they occur?

# Import numpy for vectors, matrices, and linear algebra
import numpy as np

# Set the constants hbar and omega equal to 1
hbar = omega = 1

# Define the Hamiltonian
H = (hbar*omega)/2*np.array([[1/2,0,0],[0,3/2,0],[0,0,5/2]])


# Define the state vectors
state1 = np.array([1,0,0])
state2 = np.array([0,1,0])
state3 = np.array([0,0,1])

# Define psi
psi = 1.0j*state1 + 5*state2 - 2.0j*state3

##############################
##         PART A.1         ##
##############################

# Find the magnitude of psi
psi_mag = np.linalg.norm(psi)

# Check to see if psi is normalized (if psi_mag is 1)
# If it is normalized then say so, if not normalized then
# normalize it
if psi_mag == 1.0:
    print("Psi is normalized!")
    print()
else:
    psi /= psi_mag
    print("Psi is now normalized")
    print()

Psi is now normalized

##############################
##         PART A.2         ##
##############################

# Find the eigenvalues and eigenvectors of the Hamiltonian
eigenvalues, eigenvectors = np.linalg.eig(H)

# Print the eigenvalues as they are the possible observables
print("When measuring energy, the possible observables values are the eigenvalues of the Hamiltonian.")
print("These values are:", eigenvalues)
print()

When measuring energy, the possible observables values are the eigenvalues of the Hamiltonian.
These values are: [0.25 0.75 1.25]

##############################
##         PART A.3         ##
##############################

# Need to start by finding the bra vector of psi
psi_bra = psi.T.conj()

# Compute the expectation value
expectation_H = psi_bra@H@psi

# Print the result
print("The expected (or average) value for the energy is:", np.round(expectation_H,3))
print()

The expected (or average) value for the energy is: (0.8+0j)

##############################
##         PART A.4         ##
##############################

# Do the dot products of each state vector and itself
one_one = np.dot(state1, state1)
two_two = np.dot(state2, state2)
three_three = np.dot(state3, state3)

# Do the dot products for the opposite vector pairs
# Remember that <i|j> = <j|i>^*
one_two = np.dot(state1, state2)
one_three = np.dot(state1, state3)
two_one = one_two.conj()
two_three = np.dot(state2, state3)
three_one = one_three.conj()
three_two = two_three.conj()

if one_one == two_two == three_three == 1 and one_two == one_three == two_one ==\
    two_three == three_one == three_two == 0:
    
    print("The state form a basis!")
    print()
else:
    print("The states do not form a basis.")
    print()

The state form a basis!

##############################
##         PART A.5         ##
##############################

# An operator needs to be Hermitain
O = np.array([[4,0,0], [0,5,0], [0,0,6]])

# Check to see if it is Hermitian
print("The operator is Hermitian:\n", O == O.T.conj())
print()

# Find the eigenvalues and eigenvectors of the operator
eigenvalues, eigenvectors = np.linalg.eig(O)

# Iterate through the eigenvalues and eigenvectors to find
# the observables and the probabilities
# Create a variable to store the total probability
tot_prob = 0
for i in range(len(eigenvalues)):
    # Print the observable/eigenvalue
    print ("The observable is", eigenvalues[i])

    # Extract the eigenvector corresponding to the eigenvalue
    vec = eigenvectors[:,i]

    # Compute the coefficient with the current eigenvector and the state psi
    coef = np.dot(vec.T.conj(), psi)

    # Calculate the probablity and since it must be real, only keep the real part
    prob = np.real(coef*coef.conj())

    # Add the current probability to the total probability
    tot_prob += prob

    # Print the probability
    print("The probability is", np.round(prob, 3))
    print()

print("The total probability is:", np.round(tot_prob, 3))

The operator is Hermitian:
 [[ True  True  True]
 [ True  True  True]
 [ True  True  True]]

The observable is 4.0
The probability is 0.033

The observable is 5.0
The probability is 0.833

The observable is 6.0
The probability is 0.133

The total probability is: 1.0

Part B

Consider the operators which measure spin in the x-direction, y-direction, and z-direction:

\[S_x = \frac{\hbar}{2}\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \quad S_y = \frac{\hbar}{2}\begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix} \quad S_z = \frac{\hbar}{2}\begin{bmatrix}1 & 0 \\ 0 & -1 \end{bmatrix}.\]

Now consider the following states, which could describe qubits of a quantum computer:

\[|\psi\rangle = \begin{bmatrix}3i\\5\end{bmatrix} \quad |\phi\rangle = \begin{bmatrix}4i\\0\end{bmatrix}\]

Are \(|\psi\rangle\) and \(|\phi\rangle\) normalized? If yes, prove it. If not, normalize them.
Calculate the following expectation values: \(\langle\psi|S_z|\psi\rangle\) and \(\langle\phi|S_x|\phi\rangle\).
Calculate \(\sigma_{S_z}\) (the uncertainity on the measurement of \(S_z\)) for both \(|\psi\rangle\) and \(|\phi\rangle\). Hint: \(\langle S_z^2\rangle\) means to square \(S_z\) before calculating the expectation value.
What is the probability of measuring \(S_y\) for \(|\psi\rangle\) and \(|\phi\rangle\) and obtaining a result of \(-\frac{\hbar}{2}\)?
What is the probability of measuring \(S_z\) for \(|\psi\rangle\) and \(|\phi\rangle\) and obtaining a result of \(0\)?

# Set the constant hbar equal to 1
hbar = 1

# Define the spin matrices
Sx = (hbar/2)*np.array([[0,1],[1,0]])
Sy = (hbar/2)*np.array([[0,-1.0j],[1.0j,0]])
Sz = (hbar/2)*np.array([[1,0],[0,-1]])

# Define the two states
psi = np.array([3.0j,5])
phi = np.array([4.0j,0])

##############################
##         PART B.1         ##
##############################

# PSI FIRST
# Find the magnitude of psi
psi_mag = np.linalg.norm(psi)

# Check to see if psi is normalized (if psi_mag is 1)
# If it is normalized then say so, if not normalized then
# normalize it
if psi_mag == 1.0:
    print("Psi is normalized!")
    print()
else:
    psi /= psi_mag
    print("Psi is now normalized")
    print()

## NOW PHI
# Find the magnitude of psi
phi_mag = np.linalg.norm(phi)

# Check to see if psi is normalized (if psi_mag is 1)
# If it is normalized then say so, if not normalized then
# normalize it
if phi_mag == 1.0:
    print("Phi is normalized!")
    print()
else:
    phi /= phi_mag
    print("Phi is now normalized")
    print()

Psi is now normalized

Phi is now normalized

##############################
##         PART B.2         ##
##############################

# First need to calculate the bras of both states
psi_bra = psi.T.conj()
phi_bra = phi.T.conj()

# Calculate the expecation values
expectation_Sz_psi = psi_bra@Sz@psi
expectation_Sx_phi = phi_bra@Sx@phi

# Print the expectation values
print("The expectation value of S_z and Psi is:", expectation_Sz_psi)
print("The expectation value of S_x and Phi is:", expectation_Sx_phi)

The expectation value of S_z and Psi is: (-0.2352941176470589+0j)
The expectation value of S_x and Phi is: 0j

##############################
##         PART B.3         ##
##############################

# Calculate Sz^2
Sz_squared = Sz@Sz

# PSI FIRST
# Compute the expectation values of Sz and Sz^2
expectation_Sz_psi = psi_bra@Sz@psi
expectation_Sz_squared_psi = psi_bra@Sz_squared@psi

# Calculate the uncertainity and then print it
uncertainity_psi = np.sqrt(expectation_Sz_squared_psi-expectation_Sz_psi**2)
print("The uncertainity when measuring Sz using psi is:", uncertainity_psi)
print()

# NOW PHI
# Compute the expectation values of Sz and Sz^2
expectation_Sz_phi = phi_bra@Sz@phi
expectation_Sz_squared_phi = phi_bra@Sz_squared@phi

# Calculate the uncertainity and then print it
uncertainity_phi = np.sqrt(expectation_Sz_squared_phi-expectation_Sz_phi**2)
print("The uncertainity when measuring Sz using phi is:", uncertainity_phi)
print()

The uncertainity when measuring Sz using psi is: (0.4411764705882353+0j)

The uncertainity when measuring Sz using phi is: 0j

##############################
##         PART B.4         ##
##############################

# Find the eigenvalues and eigenvectors of Sy
eigenvalues, eigenvectors = np.linalg.eig(Sy)

# Determine the index of the negative eigenvalue
print("Eigenvalues of Sy:", eigenvalues)

# Index of the negative eigenvalue
index = 1

# Get the corresponding eigenvector
state_neg = eigenvectors[:,index]

# FIRST PSI
# Find the coefficient
coeff = np.dot(state_neg, psi)

# Find the probability
prob = coeff*coeff.conj()

# Print the results
print("The probability of measuring Sy and getting -hbar/2 for psi is", prob)
print()

# NOW PHI
# Find the coefficient
coeff = np.dot(state_neg, phi)

# Find the probability
prob = coeff*coeff.conj()

# Print the results
print("The probability of measuring Sy and getting -hbar/2 for phi is", prob)
print()

Eigenvalues of Sy: [ 0.5+0.j -0.5+0.j]
The probability of measuring Sy and getting -hbar/2 for psi is (0.058823529411764705+0j)

The probability of measuring Sy and getting -hbar/2 for phi is (0.5000000000000001+0j)

##############################
##         PART B.4         ##
##############################

# Find the eigenvalues and eigenvectors of Sz
eigenvalues, eigenvectors = np.linalg.eig(Sz)

# Print the eigenvalues
print("Eigenvalues of Sz:", eigenvalues)

# 0 is not an eigenvalue of Sz
# Therefore it is not a possible observable
# So both probabilities are zero

Eigenvalues of Sz: [ 0.5 -0.5]