The harmonic oscillator is often used as an approximate model for the behaviour of some quantum systems, for example the vibrations of a diatomic molecule. The Schrödinger equation for a particle of mass $m$ moving in one dimension in a potential $V(x) = \frac{1}{2}kx^2$ is
$$
-\frac{\hbar^2}{2m}\frac{\mathrm{d}^2\psi}{\mathrm{d}x^2} + \frac{1}{2}kx^2\psi = E\psi.
$$
With the change of variable, $q = (mk/\hbar^2)^{1/4}x$, this equation becomes
$$
-\frac{1}{2}\frac{\mathrm{d}^2\psi}{\mathrm{d}q^2} + \frac{1}{2}q^2\psi = \frac{E}{\hbar\omega}\psi,
$$
where $\omega = \sqrt{k/m}$. This differential equation has an exact solution in terms of a quantum number $v=0,1,2,\cdots$:
$$
\psi(q) = N_vH_v(q)\exp(-q^2/2),
$$
where $N_v = (\sqrt{\pi}2^vv! )^{-1/2}$ is a normalization constant and $H_v(q)$ is the *Hermite polynomial* of order $v$, defined by:
$$
H_v(q) = (-1)^ve^{q^2}\frac{\mathrm{d}^v}{\mathrm{d}q^v}\left(e^{-q^2}\right).
$$
The Hermite polynomials obey a useful recursion formula:
$$
H_{n+1}(q) = 2qH_n(q) - 2nH_{n-1}(q),
$$
so given the first two: $H_0 = 1$ and $H_1 = 2q$, we can calculate all the others.

The following code plots the harmonic oscillator wavefunctions (`PLOT_PROB = False`

) or probability densities (`PLOT_PROB = True`

) for vibrational levels up to `VMAX`

on the vibrational energy levels depicted within the potential, $V = q^2/2$.

```
import numpy as np
from matplotlib import rc
from scipy.misc import factorial
import matplotlib.pyplot as plt
rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 14})
rc('text', usetex=True)
# PLOT_PROB=False plots the wavefunction, psi; PLOT_PROB=True plots |psi|^2
PLOT_PROB = False
# Maximum vibrational quantum number to calculate wavefunction for
VMAX = 6
# Some appearance settings
# Pad the q-axis on each side of the maximum turning points by this fraction
QPAD_FRAC = 1.3
# Scale the wavefunctions by this much so they don't overlap
SCALING = 0.7
# Colours of the positive and negative parts of the wavefunction
COLOUR1 = (0.6196, 0.0039, 0.2588, 1.0)
COLOUR2 = (0.3686, 0.3098, 0.6353, 1.0)
# Normalization constant and energy for vibrational state v
N = lambda v: 1./np.sqrt(np.sqrt(np.pi)*2**v*factorial(v))
get_E = lambda v: v + 0.5
def make_Hr():
"""Return a list of np.poly1d objects representing Hermite polynomials."""
# Define the Hermite polynomials up to order VMAX by recursion:
# H_[v] = 2qH_[v-1] - 2(v-1)H_[v-2]
Hr = [None] * (VMAX + 1)
Hr[0] = np.poly1d([1.,])
Hr[1] = np.poly1d([2., 0.])
for v in range(2, VMAX+1):
Hr[v] = Hr[1]*Hr[v-1] - 2*(v-1)*Hr[v-2]
return Hr
Hr = make_Hr()
def get_psi(v, q):
"""Return the harmonic oscillator wavefunction for level v on grid q."""
return N(v)*Hr[v](q)*np.exp(-q*q/2.)
def get_turning_points(v):
"""Return the classical turning points for state v."""
qmax = np.sqrt(2. * get_E(v + 0.5))
return -qmax, qmax
def get_potential(q):
"""Return potential energy on scaled oscillator displacement grid q."""
return q**2 / 2
fig, ax = plt.subplots()
qmin, qmax = get_turning_points(VMAX)
xmin, xmax = QPAD_FRAC * qmin, QPAD_FRAC * qmax
q = np.linspace(qmin, qmax, 500)
V = get_potential(q)
def plot_func(ax, f, scaling=1, yoffset=0):
"""Plot f*scaling with offset yoffset.
The curve above the offset is filled with COLOUR1; the curve below is
filled with COLOUR2.
"""
ax.plot(q, f*scaling + yoffset, color=COLOUR1)
ax.fill_between(q, f*scaling + yoffset, yoffset, f > 0.,
color=COLOUR1, alpha=0.5)
ax.fill_between(q, f*scaling + yoffset, yoffset, f < 0.,
color=COLOUR2, alpha=0.5)
# Plot the potential, V(q).
ax.plot(q, V, color='k', linewidth=1.5)
# Plot each of the wavefunctions (or probability distributions) up to VMAX.
for v in range(VMAX+1):
psi_v = get_psi(v, q)
E_v = get_E(v)
if PLOT_PROB:
plot_func(ax, psi_v**2, scaling=SCALING*1.5, yoffset=E_v)
else:
plot_func(ax, psi_v, scaling=SCALING, yoffset=E_v)
# The energy, E = (v+0.5).hbar.omega.
ax.text(s=r'$\frac{{{}}}{{2}}\hbar\omega$'.format(2*v+1), x=qmax+0.2,
y=E_v, va='center')
# Label the vibrational levels.
ax.text(s=r'$v={}$'.format(v), x=qmin-0.2, y=E_v, va='center', ha='right')
# The top of the plot, plus a bit.
ymax = E_v+0.5
if PLOT_PROB:
ylabel = r'$|\psi(q)|^2$'
else:
ylabel = r'$\psi(q)$'
ax.text(s=ylabel, x=0, y=ymax, va='bottom', ha='center')
ax.set_xlabel('$q$')
ax.set_xlim(xmin, xmax)
ax.set_ylim(0, ymax)
ax.spines['left'].set_position('center')
ax.set_yticks([])
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.savefig('sho-psi{}-{}.png'.format(PLOT_PROB+1, VMAX))
plt.show()
```

## Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

## Max Wegerhoff 2 years, 4 months ago

Hello Christian,

Link | ReplyI found this plot to be the nicest looking all over the web - may I use an adapted version on my own website? It would be much appreciated :)

Cheers,

Max Wegerhoff

## Moritz 1 year, 2 months ago

Hello, I discovered your website and I think it's amazing, the description is logic and I love that you provided the code .

Link | ReplyGreetings from Germany

Moritz

## christian 1 year, 2 months ago

Thanks, Moritz – that's very kind of you. I'm glad you found the code useful!

Link | ReplyCheers, Christian

## Doug Brown 11 months, 2 weeks ago

I get this error when running. I tried to instal "latex" but it didn't help

Link | ReplyFile "C:\Users\opabr\venv\Py10Exercises\lib\site-packages\matplotlib\texmanager.py", line 237, in _run_checked_subprocess

raise RuntimeError(

RuntimeError: Failed to process string with tex because latex could not be found

## christian 11 months, 1 week ago

I'm guessing you're using Windows? This might help: https://stackoverflow.com/questions/58121461/runtimeerror-failed-to-process-string-with-tex-because-latex-could-not-be-found – you apparently need to ensure that the Latex executable and additional applications are on your system PATH.

Link | Reply## New Comment