Wavepacket lab

The equation

A particle of mass m moving along a line is described by a complex wavefunction \psi(x, t) that obeys the time-dependent Schrödinger equation

i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \frac{\partial^2 \psi}{\partial x^2} + V(x)\,\psi .

Schrödinger wrote down this time-dependent form in the fourth communication of his 1926 series [1]. The probability of finding the particle between x and x + dx is |\psi|^2 dx, and because the right-hand side is a Hermitian operator H acting on \psi, the total probability \int |\psi|^2 dx never changes [2, Ch. 1].

The lab works in units where \hbar = m = 1. Lengths are in an arbitrary unit L, energies in \hbar^2 / m L^2 and times in m L^2 / \hbar. With these choices a plane wave e^{ikx} carries momentum p = k, kinetic energy k^2/2 and moves at velocity k. You lose nothing physical: dimensional analysis puts the constants back.

Wavepackets and the uncertainty relation

A plane wave has a perfectly sharp momentum, but it's spread over all of space, so you can't normalise it. To pin a particle down you add plane waves together into a wavepacket. The lab launches the Gaussian

\psi(x, 0) = (2\pi\sigma^2)^{-1/4} \exp\!\left[-\frac{(x - x_0)^2}{4\sigma^2} + i k_0 x\right],

whose probability density is a normal distribution centred on x_0 with standard deviation \sigma. Its Fourier transform is again Gaussian, centred on k_0 with standard deviation \sigma_k = 1/(2\sigma), so

\sigma_x \sigma_p = \frac{\hbar}{2}.

Heisenberg argued in 1927 that position and momentum cannot both be sharp [3]. Kennard proved the same year that \sigma_x \sigma_p \ge \hbar/2 for any state, with equality only for Gaussians, and worked out how Gaussian packets move for the simplest potentials [4]. Squeeze a packet in space and it spreads in momentum, so its fast and slow parts pull apart as it travels. For a free particle the width grows as

\sigma(t) = \sigma \sqrt{1 + \left(\frac{\hbar t}{2 m \sigma^2}\right)^2},

which doubles after t = 2\sqrt{3}\, m\sigma^2/\hbar [2, Ch. 2; 4]. The mean energy of the packet is \langle E \rangle = k_0^2/2 + 1/(8\sigma^2) + \langle V \rangle: the extra 1/(8\sigma^2) is the kinetic energy of the momentum spread.

In a harmonic well V = \tfrac{1}{2} m \omega^2 x^2 something special happens. A Gaussian with exactly the ground-state width \sigma = \sqrt{\hbar / 2m\omega}, displaced from the centre, oscillates like a classical particle and keeps its shape forever. Schrödinger found this in 1926 while looking for the bridge between wave mechanics and classical motion [5]; today these packets are called coherent states.

Tunnelling

Send a particle with energy E at a rectangular barrier of height V_0 and width a. A classical particle gets over if E > V_0 and bounces back if not. A quantum one doesn't stop dead at the wall. Inside the barrier its wavefunction decays as e^{-\kappa x} with \kappa = \sqrt{2m(V_0 - E)}/\hbar instead of vanishing, and whatever amplitude survives to the far side carries on as a transmitted wave. Matching \psi and \psi' at both edges gives the transmission coefficient [2, Ch. 2; 6, Ch. 6]

T = \left[1 + \frac{V_0^2 \sinh^2(\kappa a)}{4E(V_0 - E)}\right]^{-1} \quad (E < V_0),

T = \left[1 + \frac{V_0^2 \sin^2(k' a)}{4E(E - V_0)}\right]^{-1}, \quad k' = \frac{\sqrt{2m(E - V_0)}}{\hbar} \quad (E > V_0).

For a thick barrier, \kappa a \gg 1, the first form becomes T \approx 16 \tfrac{E}{V_0}(1 - \tfrac{E}{V_0}) e^{-2\kappa a}: tunnelling falls off exponentially with width and with the square root of the energy deficit. The second form holds two surprises. Above the barrier T is generally less than one, so particles that classically always pass are partly reflected. And T = 1 exactly when k' a = n\pi, when a whole number of half wavelengths fits inside the barrier. A step of height V_0 reflects too: above it, T = 4kk'/(k + k')^2, where the factor k'/k in the probability current matters as much as the amplitudes [6, Ch. 6].

Within two years of the Schrödinger equation, tunnelling cracked a real puzzle. Alpha particles leave heavy nuclei with energies of a few MeV, well below the height of the Coulomb barrier that holds them in, and their lifetimes range over many orders of magnitude while their energies change only by a factor of about two. In 1928 George Gamow [7], and independently Ronald Gurney and Edward Condon [8, 9], treated the alpha particle as a wave leaking through the barrier. The exponential factor \exp(-2\int \kappa\, dx) turns a small change in energy into a huge change in the escape rate, which accounted for the empirical relation between decay constant and range that Geiger and Nuttall had found. Friedrich Hund had analysed the related problem of a symmetric double well in his 1927 work on molecular spectra [10].

In a double well the two lowest states are the symmetric and antisymmetric combinations of the states localised in each well, split by a small energy \Delta E. A packet placed in one well is a superposition of the two, so it moves to the other well and back with period 2\pi\hbar/\Delta E. The deeper or wider the central hump, the smaller \Delta E and the slower the exchange.

Resonant tunnelling

Put two barriers in a row. Between them a particle can bounce back and forth, and at energies close to the quasi-bound levels of the gap the multiply reflected waves add up in phase. For a symmetric pair the transmission at resonance is exactly one, however opaque each barrier is on its own [11]; the lab's transfer-matrix test finds a peak above 0.999 for the preset barriers. Between resonances the pair transmits far less than either barrier alone.

Raphael Tsu and Leo Esaki computed the transmission of electrons through a finite semiconductor superlattice in 1973 and predicted resonant peaks, with negative differential resistance in the current-voltage characteristic [11]. Chang, Esaki and Tsu observed resonant tunnelling through GaAs and GaAlAs double barriers the following year [12]. The resonant tunnelling diode grew out of that line of work.

The resonance width tells you how long the particle lingers between the barriers. A narrow resonance means a long-lived quasi-bound state, and a wavepacket aimed at it leaks out of the gap slowly, in both directions.

How the simulation works

Split-operator Fourier method. The formal solution over one time step is \psi(t + \Delta t) = e^{-iH\Delta t/\hbar}\psi(t) with H = T + V. The kinetic operator T = p^2/2m is diagonal in momentum space and V is diagonal in position space, but they do not commute. Feit, Fleck and Steiger used the symmetric splitting [13]

e^{-iH\Delta t} \approx e^{-iV\Delta t/2}\, e^{-iT\Delta t}\, e^{-iV\Delta t/2},

applying each factor where it is a simple multiplication and jumping between the two spaces with a fast Fourier transform. Each factor is a phase, so the scheme is unitary and conserves the norm to rounding error. The symmetric arrangement cancels the first-order commutator term, leaving a local error of order \Delta t^3 and a global error of order \Delta t^2, as Strang showed for splittings of this kind [14]. One of the lab's tests checks that halving \Delta t cuts the phase error of a stationary state by four. Kosloff's review places the method among other time-dependent techniques [19].

The grid. The lab uses 4096 points on a periodic box of length 240, so \Delta x \approx 0.059 and the largest representable momentum is \pi/\Delta x \approx 54. The FFT is a radix-2 Cooley and Tukey transform [15] written in plain JavaScript and tested against a naive discrete Fourier transform. Each time step costs two FFTs of 4096 points.

Absorbing boundaries. A periodic box would bring a transmitted packet back in from the other side. In the outer 20 units at each end the lab adds a negative imaginary potential -iW(x) with W = 3\,s^2, where s runs from 0 to 1 across the absorbing layer. This damps the wave gradually instead of reflecting it, the approach of Kosloff and Kosloff [16] and of Neuhauser and Baer [17]. Applied as a factor e^{-W\Delta t} each step, it returns less than 10^{-3} of the probability for momenta between 0.7 and 5, which the tests check. Slower packets have wavelengths comparable to the absorber and reflect a few per cent. The lab counts every bit of probability the absorbers remove, on each side, so nothing goes missing.

Softened edges. A perfectly sharp barrier puts weight at the shortest wavelengths the grid can hold, where the splitting error is largest, and the measured \langle E \rangle then wanders by several per cent during a collision. The lab rounds each edge with a \tanh profile 0.05 units wide. That keeps \langle E \rangle within 0.5% at \Delta t = 0.005 and moves T by less than 0.01 from the sharp-edged formula.

Phase colouring. In the default view each column of the band has height |\psi|^2 and a colour set by the phase \arg\psi. The colours trace an ellipse in CIELAB that varies only lightness and the blue-yellow axis. People with protanopia or deuteranopia see colour along roughly those two directions [20], so the wheel stays readable for them, and cyclic maps should have no seam where -\pi meets \pi [21]. A packet with momentum k_0 shows bands of colour repeating every 2\pi/k_0.

The first animated wavepacket collisions were computed by Goldberg, Schey and Schwartz in 1967, who filmed Gaussian packets hitting barriers and wells [18]. This lab descends from those films, except now you get to steer.

Things to try

1. Tunnel through a thin barrier. V_0 = 2, a = 1, E \approx 1.28. About 29% of the packet appears on the far side even though every component of it has less energy than the barrier. Double the width and T drops to a few per cent. That's the exponential at work.
2. Bounce off a low barrier. With E \approx 2.9 above a barrier of 2, about a quarter reflects. Slide k_0 to find the energy where k'a = \pi and reflection vanishes.
3. Resonant tunnelling. Two barriers each transmitting about a third pass about two thirds together at the resonance near E = 0.81. Watch probability build up between the barriers and leak out on both sides. Move k_0 to 1.0 and almost everything reflects. Then run the energy sweep: the exact curve shows needle-like peaks that the packet, with its finite momentum spread, only partly resolves.
4. Spread freely. A packet with \sigma = 1 doubles its width by t \approx 3.5. The width plot lays the measured \Delta x over the exact \sigma(t); they agree to six digits, because free motion is the case where the split-operator method is exact.
5. Coherent state. In a harmonic well with \omega = 0.25, a packet of width 1/\sqrt{2\omega} \approx 1.41 released at x_0 = -20 swings with period 2\pi/\omega \approx 25 and keeps its width. Change \sigma and the width starts to breathe at twice the oscillation frequency.
6. Double well splitting. A packet resting in the left well tunnels to the right and back. T here counts the probability right of the centre, so the live plot shows the exchange directly. Raise the hump to slow it down.
7. Step up. With E \approx 2.4 above a step of 2, about a fifth reflects. Lower k_0 below 2 and the packet dips briefly into the step, an evanescent tail, before returning whole.
8. Draw your own. Drag across the canvas to paint V(x). Draw a comb of barriers and sweep T(E) to see the beginnings of band gaps, the effect Tsu and Esaki computed [11].

Reading the lab's instruments

The canvas shows the whole box from x = -120 to 120. The shaded strips at the ends are the absorbers. The heavy solid line is V(x) and the thin horizontal line is V = 0, which is also the baseline of the probability band. The dashed line marks \langle E \rangle on the same energy axis, so you can see at a glance whether the packet is above or below a barrier. The band is scaled so the initial peak of |\psi|^2 fills about a third of the space above zero, leaving room for interference peaks up to four times higher where incoming and reflected waves overlap; as the packet spreads, the band flattens. In the Re and Im view the solid line is \mathrm{Re}\,\psi, the dashed line \mathrm{Im}\,\psi and the shaded area |\psi|^2.

Norm on grid is \int |\psi|^2 dx over the box. It starts at 1 and falls only as the absorbers remove probability.

Transmitted T is the probability to the right of the potential, plus everything absorbed at the right edge. Reflected R is the same on the left. For a barrier, well or double barrier, "right" and "left" are measured from the outer edges of the potential; for the step, from the step; for the oscillator, double well and free particle, from x = 0. Once the packet has gone, T + R = 1 and T is the measured transmission.

Mean position, mean momentum and mean energy are the expectation values \langle x \rangle, \langle p \rangle and \langle E \rangle over the probability still on the grid, renormalised to it. Momentum and kinetic energy come from an extra FFT of \psi each frame. \langle E \rangle is conserved in exact quantum mechanics, so drift in it measures the time-step error: raise \Delta t to 0.05 during a collision and watch it.

Width is \Delta x, the standard deviation of |\psi|^2. Free width, exact is the analytic \sigma(t) free-particle width for the launched \sigma, meaningful only for the free particle and shown for comparison.

T and R against time follow the collision. For a single barrier both rise and settle. For a resonance T keeps creeping up while the trapped part leaks out. In the double well T oscillates.

The energy sweep fires 30 packets at energies E = k_0^2/2 from 0.25 to twice the barrier height, each with the current \sigma, and plots the final T as dots. The line is the exact plane-wave transmission, from the formula above for a barrier and from a transfer matrix for other piecewise-constant shapes, including a drawn one [6, Ch. 6]. The dots average the line over the packet's momentum spread \sigma_k = 1/(2\sigma), so sharp features come out smoothed. Make \sigma larger for sharper energy resolution, at the price of a longer packet.

References

[1] E. Schrödinger, "Quantisierung als Eigenwertproblem (Vierte Mitteilung)," Annalen der Physik 386, 109-139 (1926). DOI: 10.1002/andp.19263861802

[2] D. J. Griffiths and D. F. Schroeter, Introduction to Quantum Mechanics, 3rd ed. (Cambridge University Press, 2018). Ch. 1 "The Wave Function" and Ch. 2 "Time-Independent Schrödinger Equation", including its end-of-chapter problems on the rectangular barrier.

[3] W. Heisenberg, "Über den anschaulichen Inhalt der quantentheoretischen Kinematik und Mechanik," Zeitschrift für Physik 43, 172-198 (1927). DOI: 10.1007/BF01397280

[4] E. H. Kennard, "Zur Quantenmechanik einfacher Bewegungstypen," Zeitschrift für Physik 44, 326-352 (1927). DOI: 10.1007/BF01391200

[5] E. Schrödinger, "Der stetige Übergang von der Mikro- zur Makromechanik," Die Naturwissenschaften 14, 664-666 (1926). DOI: 10.1007/BF01507634

[6] E. Merzbacher, Quantum Mechanics, 3rd ed. (Wiley, 1998). Ch. 6 "Sectionally Constant Potentials in One Dimension."

[7] G. Gamow, "Zur Quantentheorie des Atomkernes," Zeitschrift für Physik 51, 204-212 (1928). DOI: 10.1007/BF01343196

[8] R. W. Gurney and E. U. Condon, "Wave Mechanics and Radioactive Disintegration," Nature 122, 439 (1928). DOI: 10.1038/122439a0

[9] R. W. Gurney and E. U. Condon, "Quantum Mechanics and Radioactive Disintegration," Physical Review 33, 127-140 (1929). DOI: 10.1103/PhysRev.33.127

[10] F. Hund, "Zur Deutung der Molekelspektren. I," Zeitschrift für Physik 40, 742-764 (1927). DOI: 10.1007/BF01400234

[11] R. Tsu and L. Esaki, "Tunneling in a finite superlattice," Applied Physics Letters 22, 562-564 (1973). DOI: 10.1063/1.1654509

[12] L. L. Chang, L. Esaki and R. Tsu, "Resonant tunneling in semiconductor double barriers," Applied Physics Letters 24, 593-595 (1974). DOI: 10.1063/1.1655067

[13] M. D. Feit, J. A. Fleck Jr. and A. Steiger, "Solution of the Schrödinger equation by a spectral method," Journal of Computational Physics 47, 412-433 (1982). DOI: 10.1016/0021-9991(82)90091-2

[14] G. Strang, "On the Construction and Comparison of Difference Schemes," SIAM Journal on Numerical Analysis 5, 506-517 (1968). DOI: 10.1137/0705041

[15] J. W. Cooley and J. W. Tukey, "An algorithm for the machine calculation of complex Fourier series," Mathematics of Computation 19, 297-301 (1965). DOI: 10.1090/S0025-5718-1965-0178586-1

[16] R. Kosloff and D. Kosloff, "Absorbing boundaries for wave propagation problems," Journal of Computational Physics 63, 363-376 (1986). DOI: 10.1016/0021-9991(86)90199-3

[17] D. Neuhauser and M. Baer, "The time-dependent Schrödinger equation: Application of absorbing boundary conditions," Journal of Chemical Physics 90, 4351-4355 (1989). DOI: 10.1063/1.456646

[18] A. Goldberg, H. M. Schey and J. L. Schwartz, "Computer-Generated Motion Pictures of One-Dimensional Quantum-Mechanical Transmission and Reflection Phenomena," American Journal of Physics 35, 177-186 (1967). DOI: 10.1119/1.1973991

[19] R. Kosloff, "Time-dependent quantum-mechanical methods for molecular dynamics," Journal of Physical Chemistry 92, 2087-2100 (1988). DOI: 10.1021/j100319a003

[20] H. Brettel, F. Viénot and J. D. Mollon, "Computerized simulation of color appearance for dichromats," Journal of the Optical Society of America A 14, 2647-2655 (1997). DOI: 10.1364/JOSAA.14.002647

[21] P. Kovesi, "Good Colour Maps: How to Design Them," arXiv: 1509.03700 (2015)