The Harmonic Oscillator

Here is a piece of luck so large it looks like a conspiracy. Take any potential with a stable resting point — a chemical bond, an atom in a crystal, a pendulum, the charge in a circuit — and zoom in on the bottom of its valley. Taylor says:

The constant shifts nothing; the slope vanishes because the point is a minimum; and what remains — the first term with any content — is a parabola. Every gentle vibration in nature, near its resting point, is the same problem:

Solve this Hamiltonian once and you have pre-solved the vibrations of every molecule, the sound waves in every solid (phonons, Chapter 15's guests), the sloshing of currents in circuits — and, when Part IV quantizes the electromagnetic field, each mode of light will turn out to be exactly this oscillator, with photons as its rungs. No other page of algebra in physics pays a dividend like this one. So let us do it properly — and by “properly” we mean: without solving any differential equation at all.

Look at (7.2) with an algebraist's eye: it is a sum of two squares, and sums of squares ache to be factored, $u^2+v^2=(u-iv)(u+iv)$. Try it. Define the natural dimensionless combination (one for each square) and its product:

The factorization almost works — the leftover ½ is the toll charged by $[\hat{x},\hat{p}]=i\hslash$, because $u$ and $v$ here do not commute. Work out the same toll between the factors themselves and you get the one equation this whole chapter runs on:

Now watch the card trick. Compute how $\hat{H}$ talks to $\hat{a}$: $[\hat{H},\hat{a}]=-\hslash \omega \hat{a}$. Read that as a sentence: if $|E⟩$ is an energy eigenstate, then

so $\hat{a}|E⟩$ is another eigenstate, one rung down. The operator $\hat{a}$ descends the spectrum in steps of $\hslash \omega$; ${\hat{a}}^{†}$ climbs it. But the descent cannot go on forever — energies are bounded below, because $⟨{\hat{a}}^{†}\hat{a}⟩=\Vert \hat{a}|\psi ⟩{\Vert }^2\ge 0$ (a squared length, the same non-negativity that proved the uncertainty principle). The only escape from an infinite descent is a state the ladder annihilates:

That is the entire solution. An evenly spaced ladder — compare the box's widening $n^2$ staircase — with a ground rung lifted $\hslash \omega /2$ off the floor: the zero-point energy, now as an algebraic inevitability rather than Chapter 6's wave argument. Keeping the bookkeeping of norms, the ladder's matrix elements are $\hat{a}|n⟩=\sqrt{n}|n-1⟩$ and ${\hat{a}}^{†}|n⟩=\sqrt{n+1}|n+1⟩$. Notice what never happened: we never solved for a single wavefunction. The spectrum came from one commutator and one inequality. When Part IV needs photons, this algebra — verbatim — will be the theory of light; the operator ${\hat{a}}^{†}$ will be the thing a glowing atom does.

The wavefunctions are now an afterthought — which is the luxury the algebra bought us. The bottom rung obeys $\hat{a}|0⟩=0$, which in position language is a first-order equation, $\left(x+\frac{\hslash }{m\omega }\frac{d}{dx}\right){\psi }_0=0$, and first-order equations surrender in one line:

— a Gaussian, and not just any Gaussian: the minimum-uncertainty state of Chapter 4, sitting exactly at $\Delta x\Delta p=\hslash /2$. The ground state of nature's favorite potential is the most certain state the uncertainty principle permits. Every higher rung is a stroke of the raising operator — each application of ${\hat{a}}^{†}$ multiplies by another power of x (minus a derivative), building the Hermite polynomials — the specific family of polynomials $H_n(x)$ that multiply the Gaussian in the $n$-th eigenstate, ${\psi }_n\propto H_n(x){\psi }_0$, and whose $n$ zeros are exactly the wavefunction's $n$ nodes — one wiggle at a time. Before you climb, here is the whole picture in one frame:

Climb high enough and the quantum oscillator must shake hands with the classical one — but the handshake is stranger than “they look alike.” Ask a classical pendulum: where do you spend your time? Mostly at the edges, where it moves slowest — a stroboscope photo piles up at the turning points. Now put the n = 20 eigenstate (from the lab's matrix method) on top of that prediction:

Do not skate past what this picture means. The quantum comb and the classical tent are not approximately the same curve — the wiggles never go away — yet every measurable average built from them agrees, and agrees better the higher you climb. That is the correspondence principle with a face: classical mechanics is not overthrown by the quantum world, it is the quantum world seen coarsely, at large quantum numbers where the interference fringes are finer than any ruler. And the trick of reading a classical trajectory off the smooth envelope of a high-$n$ wavefunction is not just a sanity check — sharpened into the WKB approximation, it becomes a working tool: a way to get energy levels, tunneling rates, and scattering phases for potentials far too messy to solve exactly, by threading a semiclassical wave through the classically allowed region. Squinting, done carefully, is a method.

One puzzle remains, and it is the best one. Every energy eigenstate is a statue — yet real springs swing. Chapter 5 told us where motion lives: in superpositions, beating at level differences. For the oscillator every level difference is the same $\omega$ — all the beats agree! — so the oscillator can support motion of unprecedented discipline. Which superposition swings most faithfully? Here is a candidate with impeccable manners: ask for an eigenstate of the lowering operator,

A strange request — $\hat{a}$ is not Hermitian, so $\alpha$ may be complex, and no measurement “reads” it directly. Its rung populations, though, are perfectly ordinary: squaring the coefficients in (7.8) gives $P(n)=e^{-|\alpha {|}^2}|\alpha {|}^{2n}/n!$, the Poisson distribution — the law obeyed by counts of independent rare events (raindrops on a tile, clicks in a counter), here spreading the coherent state over the ladder's rungs with mean $|\alpha {|}^2$. But the payoff is immediate. Under time evolution each $|n⟩$ turns at $e^{-in\omega t}$ (dropping the common zero-point phase), and the sum in (7.8) re-assembles itself exactly:

A coherent state stays a coherent state forever — its label just circles the complex plane at the classical frequency. And that complex plane is nothing exotic: the real and imaginary parts of $\alpha$ are the (scaled) mean position and mean momentum, so the plane is exactly phase space — the position–momentum plane in which classical mechanics draws its orbits. The coherent state traces the same circle Newton would. Concretely: a rigid Gaussian of ground-state width whose center runs Newton's trajectory, with $⟨x⟩(t)=\sqrt{2\hslash /m\omega }|\alpha |\cos (\omega t-\phi )$. No spreading, ever — the parabola's linear force refocuses precisely as fast as dispersion defocuses. Fly it:

And the making of one is almost comically easy: displace the ground state and let go — which is what plucking a guitar string, kicking a pendulum, or driving a microwave cavity does. That is why the classical world is full of coherent states: they are what happens when anything pushes on an oscillator's vacuum. When a laser does the pushing to a mode of light, the Poisson bars above become the photon counting statistics of the beam — measured in every quantum optics lab, exactly Poissonian, mean $|\alpha {|}^2$.

In 1925 Heisenberg solved this system with infinite matrices; months later Schrödinger solved it with a wave equation; the answers agreed and nobody at first understood why (Chapter 3 knew: one arrow, two shadows).

The lab Rust-QP/ch07-oscillator restages the showdown with modern referees. Method one is Heisenberg's, updated: the Hamiltonian as a 1600×1600 symmetric matrix, handed to an eigensolver —

1// ---------------------------------------------------------------------------
2// Method 1: matrix mechanics — diagonalize the discretized Hamiltonian
3// ---------------------------------------------------------------------------
4// …
5fn matrix_mechanics(n: usize, l: f64, keep_states: &[usize]) -> MatrixResult {
6    let dx = l / n as f64;
7    let x: Vec<f64> = (0..n).map(|i| -l / 2.0 + (i as f64 + 0.5) * dx).collect();
8
9    // Kinetic term -psi''/2 via a 4TH-ORDER 5-POINT Laplacian stencil.
10    // The second derivative is approximated by
11    //     psi''(x_i) ~ (-psi_{i-2} + 16 psi_{i-1} - 30 psi_i
12    //                    + 16 psi_{i+1} - psi_{i+2}) / (12 h^2),
13    // whose leading error is O(h^4) — versus O(h^2) for the plain
14    // 3-point [1, -2, 1]/h^2 stencil. On the same 1600-point grid this
15    // buys ~2 extra digits of accuracy in the eigenvalues for free.
16    // Multiplying by -1/2 gives the pentadiagonal kinetic operator:
17    //     diagonal      : -1/2 * (-30/12) / h^2 =  5/4  / h^2
18    //     first  off-diag: -1/2 * ( 16/12) / h^2 = -2/3  / h^2
19    //     second off-diag: -1/2 * ( -1/12) / h^2 =  1/24 / h^2
20    // Boundary rows: the bound states decay to machine-zero long before
21    // the wall (grid half-width 15 >> classical turning point ~6.4 for
22    // n=20), so we impose Dirichlet walls — psi = 0 outside the grid —
23    // by simply omitting any stencil neighbor that falls off the ends.
24    let inv_h2 = 1.0 / (dx * dx);
25    let d0 = 5.0 / 4.0 * inv_h2;
26    let d1 = -2.0 / 3.0 * inv_h2;
27    let d2 = 1.0 / 24.0 * inv_h2;
28    let mut h = DMatrix::<f64>::zeros(n, n);
29    for i in 0..n {
30        h[(i, i)] = d0 + 0.5 * x[i] * x[i];
31        if i + 1 < n {
32            h[(i, i + 1)] = d1;
33            h[(i + 1, i)] = d1;
34        }
35        if i + 2 < n {
36            h[(i, i + 2)] = d2;
37            h[(i + 2, i)] = d2;
38        }
39    }
40    let eig = SymmetricEigen::new(h);
41    // …
42}

Method two is Schrödinger's: Chapter 6's split-operator engine descending in imaginary time. Method three tests section 7.5's boldest claim in real time — that a displaced Gaussian obeys Newton without spreading:

1fn coherent(x0: f64) -> (Vec<CoherentPoint>, CoherentStats) {
2    let eng = WaveEngine::new(1024, 40.0);
3    // displaced ground state: the definition of a coherent state
4    let mut psi: Vec<C> = eng
5        .x
6        .iter()
7        .map(|&x| c((-(x - x0) * (x - x0) / 2.0).exp() / PI.powf(0.25), 0.0))
8        .collect();
9    let dt = 0.002;
10    let t_end = 4.0 * PI; // two classical periods
11    let steps = (t_end / dt) as usize;
12    let sigma_exact: f64 = 0.5f64.sqrt();
13    // <H> for a coherent state displaced to x0: the classical energy of the
14    // turning point plus the zero-point half-quantum. Analytic, so the
15    // integrator is being marked against something it never computed.
16    let energy_exact: f64 = 0.5 * x0 * x0 + 0.5;
17    let mut out = Vec::new();
18    let mut st = CoherentStats { worst_x: 0.0, worst_sigma: 0.0, worst_norm: 0.0, worst_energy: 0.0 };
19    for i in 0..=steps {
20        let t = i as f64 * dt;
21        if i % (steps / 80) == 0 {
22            let (mean, sigma) = eng.moments(&psi);
23            let classical = x0 * t.cos();
24            st.worst_x = st.worst_x.max((mean - classical).abs());
25            st.worst_sigma = st.worst_sigma.max((sigma - sigma_exact).abs());
26            st.worst_norm = st.worst_norm.max((eng.norm2(&psi) - 1.0).abs());
27            st.worst_energy = st.worst_energy.max((eng.energy(&psi) - energy_exact).abs());
28            out.push(CoherentPoint { t, x_mean: mean, classical, sigma });
29        }
30        eng.step(&mut psi, c(dt, 0.0));
31    }
32    (out, st)
33}

Notice what that loop watches besides $⟨x⟩$ and σ. It also re-measures the norm and the energy at every sample, inside the potential — Chapter 6's norm check ran where V ≡ 0, which only exercises the FFT round trip, whereas here the potential half-steps are in the loop too. And it marks $⟨H⟩$ against $x_0^2/2+\frac{1}{2}$, a number the integrator never computes and could not have tuned itself to. Then eight referees, finiteness before smallness, tolerances set to what the run actually achieves:

1// -- the referees ------------------------------------------------------
2    // Finiteness FIRST — a NaN compares false against every bound and would
3    // otherwise sail through as a pass.
4    for (label, v) in [
5        ("worst_matrix", worst_matrix),
6        ("worst_wave", worst_wave),
7        ("worst_cross", worst_cross),
8        ("worst_x", cs.worst_x),
9        ("worst_sigma", cs.worst_sigma),
10        ("worst_norm", cs.worst_norm),
11        ("worst_energy", cs.worst_energy),
12        ("worst_corr", worst_corr),
13    ] {
14        assert!(v.is_finite(), "{label} is not finite: {v}");
15    }
16
17    // Tolerances are the honestly-achieved values with a small factor of
18    // headroom for platform rounding — never a round number chosen to pass.
19    let mut referees: Vec<Referee> = Vec::new();
20    let mut push = |name: &str, value: f64, tol: f64| {
21        referees.push(Referee { name: name.into(), value, tol, pass: value < tol });
22    };
23    // 4th-order 5-point Laplacian: worst error is ~1.5e-6 on this grid, down
24    // from ~2e-3 with the old 3-point stencil.
25    push("Matrix mechanics vs the algebra: worst |E_n − (n+½)|, n < 10", worst_matrix, 2e-6);
26    push("Wave mechanics vs the algebra: worst |E_n − (n+½)|, n < 4", worst_wave, 5e-11);
27    push("Matrix vs wave, head to head: worst |E_n^mat − E_n^wave|, n < 4", worst_cross, 2e-7);
28    // limited by the O(dt^2) Trotter error at dt = 0.002, not by physics
29    push("Coherent ⟨x⟩(t) vs Newton's x₀ cos t, two periods", cs.worst_x, 1e-5);
30    push("Coherent width σ(t) vs 1/√2: the packet is rigid", cs.worst_sigma, 1e-6);
31    push("Real-time unitarity ‖ψ‖² = 1 inside the parabola", cs.worst_norm, 3e-12);
32    push("Real-time energy ⟨H⟩(t) vs x₀²/2 + ½", cs.worst_energy, 1e-5);
33    push("Correspondence: P(|x| < f·A) at n = 20 vs the arcsine law", worst_corr, 3e-2);

1. (F) Verify $[\hat{a},{\hat{a}}^{†}]=1$ from $[\hat{x},\hat{p}]=i\hslash$, and compute $[\hat{H},{\hat{a}}^{†}]$ to confirm the raising direction.
2. (F) Using only the ladder rules, compute $⟨n|\hat{x}|m⟩$ and show it vanishes unless $m=n\pm 1$. This selection rule is why a vibrating molecule absorbs light at ω and (to first order) nowhere else — infrared spectroscopy in one matrix element.
3. (C) In the ladder widget, superpose n = 0 with n = 1 and time the beat against the n = 4 with n = 5 superposition. Same period? Explain why in one sentence — and say what property of the spectrum would make them differ (look ahead: every real molecule's anharmonicity).
4. (F, hard) Pure ladder-algebra workout: write ${\hat{x}}^4$ in terms of $\hat{a},{\hat{a}}^{†}$ and show $⟨n|{\hat{x}}^4|n⟩=\frac{3}{4}{\left(\frac{\hslash }{m\omega }\right)}^2(2n^2+2n+1)$. (Bookkeeping discipline: only terms with equal numbers of raisings and lowerings survive.) You have just computed the raw ingredient of the anharmonic correction that Chapter 12 will turn into real molecular spectra.
5. (F, hard) Prove the three claims of (7.8) and (7.9): that the displaced-vacuum expansion satisfies $\hat{a}|\alpha ⟩=\alpha |\alpha ⟩$; that $P(n)=|⟨n|\alpha ⟩{|}^2$ is Poisson with mean $|\alpha {|}^2$; and that time evolution maps $|\alpha ⟩\to e^{-i\omega t/2}|\alpha e^{-i\omega t}⟩$. Then show two coherent states are never orthogonal: $|⟨\beta |\alpha ⟩{|}^2=e^{-|\alpha -\beta {|}^2}$ — they overfill the Hilbert space (an “overcomplete basis,” the working language of quantum optics).
6. (P) Squeezed states: start the lab's real-time run with a Gaussian of the wrong width (σ = 2 × ground). Watch σ(t) breathe at frequency 2ω — derive that factor of two, then verify the breathing minimum and maximum against your formula.
7. (P, hard) Build Schrödinger's cat: the superposition $(|\alpha ⟩+|-\alpha ⟩)/\mathcal{N}$ with α = 3, as two displaced Gaussians. Evolve it. Twice per period the halves pass through each other at the origin — plot |ψ|² at that instant and find the interference fringes (spacing π/α√2... derive it!). Then add Chapter 4's γ-dial as a crude decoherence model and watch the fringes — the state's “catness” — die while the two humps survive. You have simulated why big superpositions are hard and big mixtures are everywhere.