# Wave Packets



In [1]:

%pylab inline




Populating the interactive namespace from numpy and matplotlib




In [2]:

import matplotlib.animation




In [3]:

from IPython.display import HTML



A particle with total energy $E$ in a region of constant potential $V_0$ has a wave number $$k = \pm \frac{2m}{\hbar^2}(E - V_0)$$ and dispersion relation $$\omega(k) = \frac{E(k)}{\hbar} = \frac{\hbar k^2}{2m} + \frac{V_0}{\hbar}$$ leading to phase and group velocities $$v_p(k) = \frac{\omega(k)}{k} = \frac{\hbar k}{2 m} + \frac{V_0}{\hbar k} \quad, \quad v_g(k) = \frac{d\omega}{dk}(k) = \frac{\hbar k}{m} \; .$$

Consider an initial state at $t=0$ that is a normalized Gaussian wavepacket $$\Psi(x,0) = \pi^{-1/4} \sigma_x^{-1/2} e^{i k_0 x}\, \exp\left(-\frac{1}{2} \frac{x^2}{\sigma_x^2}\right)$$ with a central group velocity $v_g = \hbar k_0 / m$.

We can expand this state's time evolution in plane waves as $$\Psi(x,t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} c(k) \exp\left[ i (k x - \omega(k) t)\right]\, dk$$ with coefficients $$c(k) = \frac{1}{\sqrt{2\pi}}\,\int_{-\infty}^{+\infty} \Psi(x,0)\, e^{-i k x} dx = \pi^{-1/4} \sigma_x^{1/2}\, \exp\left( -\frac{1}{2} (k - k_0)^2 \sigma_x^2\right) \; .$$

Approximate the integral over $k$ with a discrete sum of values $k_i$ centered on $k_0$ then the (un-normalized) real part of the wave function is $$\text{Re}\Psi(x,t) \simeq \sum_{i=-N}^{+N} c(k_i) \cos\left[ k x_i - \left(\frac{\hbar k_i^2}{2m} + \frac{V_0}{\hbar}\right) t \right] \; .$$



In [4]:

def solve(k0=10., sigmax=0.25, V0=0., mass=1., tmax=0.25, nwave=15, nx=500, nt=10):
"""
Solve for the evolution of a 1D Gaussian wave packet.

Parameters
----------
k0 : float
Central wavenumber, which determines the group velocity of the wave
packet and can be negative, zero, or positive.
sigmax : float
V0 : float
Constant potential in units of hbar. Changes the (unphysical)
phase velocities but not the group velocity.
mass : float
Particle mass in units of hbar.
tmax : float
Amount of time to simulate on a uniform grid in arbitrary units.
nwave : int
Wave packet is approximated by 2*nwave+1 plane waves centered on k0.
nx : int
Number of grid points to use in x.
nt : int
Number of grid points to use in t.
"""
t = np.linspace(0, tmax, nt).reshape(-1, 1)

# Calculate the group velocity at k0.
vgroup0 = k0 / mass

# Calculate the distance traveled by the wave packet.
dist = np.abs(vgroup0) * tmax

# Calculate the spreading of the wave packet.
spread = np.sqrt((sigmax ** 4 + (t / mass) ** 2) / sigmax ** 2)

# Calculate an x-range that keeps the packet visible during tmax.
nsigmas = 1.5
xrange = max(tails + dist, 2 * tails)
x0 = nsigmas * spread[0] + 0.5 * (xrange - tails) - 0.5 * vgroup0 * tmax - 0.5 * xrange
x = np.linspace(-0.5 * xrange, 0.5 * xrange, nx) - x0

# Build grid of k values to use, centered on k0.
nsigmas = 2.0
sigmak = 1. / sigmax
k = k0 + sigmak * np.linspace(-nsigmas, +nsigmas, 2 * nwave + 1).reshape(-1, 1, 1)

# Calculate coefficients c(k).
ck = np.exp(-0.5 * (k - k0) ** 2 * sigmax ** 2)

# Calculate omega(k)
omega = k ** 2 / (2 * mass) + V0

# Calculate the (un-normalized) evolution of each wavenumber.
psi = ck * np.cos(k * x - omega * t)

# Calculate the (x,y) coordinates of a tracer for each wavenumber.
xtrace = np.zeros((nt, 2 * nwave + 1))
nz = k != 0
xtrace[:, nz.ravel()] = ((k[nz] / (2 * mass) + V0 / k[nz]) * t)
ytrace = ck.reshape(-1)

# Calculate the motion of the center of the wave packet.
xcenter = vgroup0 * t

return x, psi, xtrace, ytrace, xcenter



Build an animation using the solver above. Each frame shows:

• The real part of the component plane waves in blue, with increasing $k$ (decreasing $\lambda$) moving up the plot. Each component is vertically offset and has an amplitude of $c(k)$.
• The sum of each plane is shown in red and represents the real part of the wave packet wave function.
• Blue dots trace the motion of each plane from the initial wave packet center, each traveling at their phase velocity $v_p(k)$.
• A red vertical line moves with the central group velocity $v_g(k_0)$.

The main points to note are:

• In the intial frame, the blue plane waves all interfere constructively at the center of the wave packet, but become progressively more incoherent moving away from the peak.
• Each blue plane wave propagates at a different phase velocity $v_p(k)$, causing the blue tracer dots to separate horizontally over time.
• Changes in the constant potential $V_0$ lead to different phase velocities but an identical combined red wave function and group velocity.
• The center of the wave packet travels at the central group velocity $v_g(k_0)$ indicated by the vertical red line.
• The red wave packet spreads as it propagates.


In [5]:

def animate(k0=10., sigmax=0.25, V0=0., mass=1., nt=30, save=None, height=480, width=720):

x, psi, xt, yt, xc = solve(k0, sigmax, V0, mass, nt=nt)
nw, nt, nx = psi.shape
nwave = (nw - 1) // 2
psi_sum = np.sum(psi, axis=0)

ymax = 1.02 * np.max(np.abs(psi_sum))
dy = 0.95 * ymax * np.arange(-nwave, nwave + 1) / nwave
assert len(dy) == nw

artists = []
dpi = 100.
dot = 0.006 * height
figure = plt.figure(figsize=(width / dpi, height / dpi), dpi=dpi, frameon=False)
ax = figure.add_axes([0, 0, 1, 1])
plt.axis('off')
for i in range(2 * nwave + 1):
artists += ax.plot(x, psi[i, 0] + dy[i], lw=2 * dot, c='b', alpha=0.2)
artists.append(ax.axvline(xc[0], c='r', ls='-', lw=dot, alpha=0.4))
artists += ax.plot(xt[0], yt + dy, 'b.', ms=2.5 * dot, alpha=0.4, lw=0)
artists += ax.plot(x, psi_sum[0], 'r-', lw=2.5 * dot, alpha=0.5)
ax.set_xlim(x[0], x[-1])
ax.set_ylim(-ymax, +ymax)

def init():
return artists

def update(j):
for i in range(2 * nwave + 1):
artists[i].set_ydata(psi[i, j] + dy[i])
artists[-3].set_xdata([xc[j], xc[j]])
artists[-2].set_xdata(xt[j])
artists[-1].set_ydata(psi_sum[j])
return artists

animation = matplotlib.animation.FuncAnimation(
figure, update, init_func=init, frames=nt, repeat=True, blit=True)
if save:
meta = dict(
title='Gaussian quantum wavepacket superposition in 1D',
artist='David Kirkby <dkirkby@uci.edu>',
comment='https://dkirkby.github.io/quantum-demo/',
engine = 'imagemagick' if save.endswith('.gif') else 'ffmpeg'
animation.save(save, writer=writer)

return animation



Simulate with different values of the constant potential: \begin{aligned} V_0 &= -50 \quad &\Rightarrow \quad v_p(k_0) &= 0 \\ V_0 &= 0 \quad &\Rightarrow \quad v_p(k_0) &= v_g(k_0) / 2 \\ V_0 &= +50 \quad &\Rightarrow \quad v_p(k_0) &= v_g(k_0) \\ \end{aligned}



In [6]:

animate(V0=-50, nt=300, save='wavepacket0.mp4');







In [7]:

animate(V0=0, nt=300, height=480, width=720, save='wavepacket1.mp4');







In [8]:

animate(V0=50, nt=300, height=480, width=720, save='wavepacket2.mp4');






Uncomment and run the line below to display an animation inline:



In [9]:

#HTML(animate().to_html5_video())



Convert video to the open-source Theora format using, e.g.

ffmpeg -i wavepacket0.mp4 -codec:v libtheora -qscale:v 7 wavepacket0.ogv

Produce a smaller GIF animation for wikimedia. Note that this file is slightly larger (~1Mb) than the MP4 files, despite the lower quality.



In [10]:

animate(V0=0, nt=100, height=150, width=200, save='wavepacket1.gif');