Reminder: HO energies given by $$ E_n = \left(n + \frac{1}{2}\right)\hbar\omega. $$ Working with dimensionless quantities: $$ x^\prime \equiv \frac{x}{a},$$ where $a\equiv \sqrt{\hbar/(m\omega)}$, and $$ E^\prime \equiv \frac{E}{\hbar\omega}, $$ and $$ t^\prime \equiv \frac{t}{T}. $$ This means that $$e^{-iEt/\hbar} \rightarrow e^{-i 2\pi E^\prime t^\prime},$$ effectively setting $h\rightarrow 1$.
Harmonic oscillator eigenfunctions given by
$$ \psi_n(x) = \exp\left(-\frac{m\omega}{2\hbar}x^2\right)\, h_n\left(\sqrt{\frac{m\omega}{\hbar}}x\right) \frac{\sqrt[4]{m\omega/\hbar}}{\sqrt{2^n n! \sqrt{\pi}}}, $$where $h_n(x)$ is the $n^\text{th}$ Hermite polynomial.
Switching to dimensionless length parameter $x^\prime$ the expression for the eigenfunctions simplifies to
$$ \psi_n(x^\prime) = \exp\left(-\frac{1}{2}(x^\prime)^2\right)\, h_n\left(x^\prime\right) \frac{1}{\sqrt{2^n n! \sqrt{\pi}}} $$The so-called coherent states of a harmonic oscillator are special linear combinations of the energy eigenstates that have several nice properties:
import numpy as np
from scipy import special
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Following is an Ipython magic command that puts figures in the notebook.
%matplotlib notebook
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size
def psiHO(n,x):
'''Harmonic oscillator wavefunctions of eigenstates. Uses:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.eval_hermite.html'''
return special.eval_hermite(n, x)*np.exp(-x**2/2.)/ \
np.sqrt(2**n*special.factorial(n)*np.sqrt(np.pi))
def psiTotal(x,t):
'''Total time-dependent wave function as sum of eigenstates with appropriate
time dependences'''
alpha = 2.
m = 40
total = 0
for n in range(m):
total += np.exp(-1j*2*np.pi*(0.5+n)*t)*alpha**n*psiHO(n,x)/ \
np.sqrt(special.factorial(n))
return np.exp(-alpha**2/2)*total
def animate(i): # Update time, title (containing time), and data
dt = 0.01
t0 = 0
t = t0 + i*dt
plt.title("t={0:.2f}".format(t))
P.set_ydata(abs(psiTotal(x,t))**2) # update the data
return P,
fig = plt.figure()
N = 50 # number of intervals
x = np.linspace(-5, 5, N+1, endpoint=True)
plt.grid(True)
plt.ylim(-0.2, 0.7)
plt.xlabel('$x$') # Label for horizontal axis
plt.ylabel("$|\psi|^2$") # Label for vertical axis
plt.axhline(0,color='green') # Makes solid green x-axis
y = psiTotal(x,0)
P, = plt.plot(x, abs(y)**2)
ani = animation.FuncAnimation(fig, animate, interval = 100)
# interval --> time between frames
%version_information
is an IPython magic extension for showing version information for dependency modules in a notebook;
%version_information
is available on Bucknell computers on the linux network. You can easily install it on any computer.
%load_ext version_information
version_information numpy, scipy, matplotlib