Almost-Solvable Worlds

In 1887, Michelson — the man whose interferometers ended the luminiferous aether — pointed one at the red Balmer line of hydrogen and found that it is not a line. It is a doublet: two lines about 0.016 nm apart. Set that against the pattern Balmer had fit — Hα to Hβ is 170 nm — and the doublet is ten thousand times finer than the structure the formula was built to explain. The exact solution of Chapter 10, the one that matched the telescope to five figures, says nothing about this. And here is the delicious part: remember that our own telescope test carried a worst-case miss of a hundredth of a nanometer? That residual was not numerical noise. It was physics knocking — the fine structure, sitting exactly at the scale where our Bohr-formula precision ran out.

This is how real physics almost always presents itself: an exactly solvable skeleton — box, oscillator, hydrogen — wearing corrections like layers of increasingly fine clothing. Relativity tweaks the electron's kinetic energy. The electron's spin feels the atom's own internal magnetism. The proton's tiny magnet whispers to the electron's. Each layer is small, and that smallness is not an obstacle — it is a gift, because small things can be computed order by order. The machinery for doing so, perturbation theory, is the single most-used tool in all of physics; most working physicists spend most of their lives inside it. Time to build it.

The setup: a Hamiltonian $\hat{H}={\hat{H}}_0+\lambda {\hat{H}}^{\prime }$ whose first piece we have solved completely — ${\hat{H}}_0|n⟩=E_n|n⟩$ — and a dial λ we imagine turning up from zero. Expand everything in powers of λ,

Before the algebra, hold the whole story as a picture. Figure 12.1 is what turning the dial λ up from zero actually does to a spectrum: every level of the solved ${\hat{H}}_0$ slides a little (a shift), and any level that packed several states at one energy fans them out (a split). The concrete case drawn here — hydrogen's degenerate 2p level splitting into $2p_{3/2}$ and $2p_{1/2}$ — is the fine structure of §12.4, and it is exactly Michelson's doublet from §12.1. Everything below is just how to compute the sizes of those shifts and splits.

feed the expansion into $\hat{H}|\psi ⟩=E|\psi ⟩$, and sort by powers of λ like sorting coins. At order λ, project onto $⟨n|$ and the unknown $|{\psi }^{(1)}⟩$ drops out entirely (exercise 1 — two lines), leaving the most useful formula in quantum mechanics:

averaged over the unperturbed state. No new problem to solve; just an integral over a wavefunction you already own. Project instead onto the other states and you get the first-order wavefunction — the state leans toward its neighbors, each admixture weighted by coupling over distance — and feeding that lean back in gives the second-order energy:

Read (12.3) the way you would read a sociology of energy levels. Every other state pushes on level n, with strength (coupling)² over (gap): levels repel. A neighbor above pushes you down; a neighbor below pushes you up; the ground state, with everyone above it, is pushed only down — second order never raises a ground state, a theorem you get for free from the sign of the denominator. And the formula announces its own limits: it is a good approximation exactly when couplings are small compared to gaps, $|⟨m|{\hat{H}}^{\prime }|n⟩|\ll |E_n-E_m|$. Which raises an awkward question — what if two levels share the same energy, and the gap is zero?

Then (12.3) divides by zero, and the failure is trying to tell you something true: when several states share an energy, the question “which state does the correction apply to?” is ill-posed. Any blend of degenerate states is as good an eigenstate as any other — until the perturbation arrives and prefers some blends to others. The cure is to ask the perturbation first: diagonalize ${\hat{H}}^{\prime }$ within the degenerate subspace. Its eigenvectors are the blends that survive contact with the nuisance; its eigenvalues are their first-order splittings; and with the right basis chosen, the divide-by-zero terms simply never arise.

You have already done this without the name. Chapter 11's total-j states are exactly the blends that diagonalize the spin-orbit coupling $\hat{\vec{L}}\cdot \hat{\vec{S}}$ inside a degenerate hydrogen level — the CG machinery was degenerate perturbation theory's basis-picker, built one chapter early. Watch both regimes live:

Now collect winnings. Ride along with the electron: from its seat, the proton circles it, and a circulating charge is a current loop — the electron sits in a magnetic field of the atom's own making, and its spin magnet pays an energy ${\hat{H}}_{SO}=\xi (r)\hat{\vec{L}}\cdot \hat{\vec{S}}$ (the spin-orbit coupling, ξ ∝ α²/r³ in atomic units). Degenerate PT says: diagonalize $\hat{\vec{L}}\cdot \hat{\vec{S}}$ in the 2p level — and Chapter 11 already did, handing us the j-basis with eigenvalues from the operator identity $\hat{\vec{L}}\cdot \hat{\vec{S}}=\frac{1}{2}({\hat{J}}^2-{\hat{L}}^2-{\hat{S}}^2)$:

The lab shoots the 2p wavefunction (Chapter 10's Numerov, unchanged), integrates $⟨r^{-3}⟩$ numerically, and lands the interval within the 0.2% its referee enforces of the measured 10.969 GHz — the panel below prints the number the run actually produced. What is left over has a name: the electron's anomalous g-factor plus QED, physics we gladly leave outstanding. Translated to wavelength, this splitting smears Hα by about 0.016 nm: Michelson's doublet, and our Chapter 10 residual, both accounted for.

One piece of bookkeeping before the winnings are banked, because it is where level diagrams usually lie. $\hat{\vec{L}}\cdot \hat{\vec{S}}$ is traceless — its shifts, weighted by degeneracy, cancel — so spin-orbit coupling alone would leave the doublet straddling the Bohr level, one member up and one down. It does not, because spin-orbit is not travelling alone. The electron's relativistic kinetic correction $-{\hat{p}}^4/8m^3{c}^2$ is the same order in α² and depends on n and l but not on j, so it slides the whole doublet bodily downward. Add the two and both members come out below $E_2$, at $-\frac{1}{4}\Delta E$ and $-\frac{5}{4}\Delta E$. The 1s level, with no orbital angular momentum, gets no spin-orbit term at all — but it does get the relativistic one plus a contact (Darwin) term, and their sum is $4\Delta E$: the level that cannot split moves four times as far as the whole splitting we just computed. The interval is what a spectrometer sees, which is why the lab computes that and the figures draw the rest to scale.

Then the deeper whisper. The proton is also a magnet, 658 times feebler; its interaction with the electron's spin at the nucleus (the Fermi contact term, ∝ $|\psi (0){|}^2\hat{\vec{S}}\cdot \hat{\vec{I}}$) splits Chapter 11's singlet from its triplet. Only s-states feel it — only they touch the origin — and the splitting is the 21 cm quantum:

The lab prices it from the shot 1s wavefunction and five physical constants, and its referee holds the answer to 0.01% of the 1420.406 MHz that radio astronomy is built on — four significant figures from a wavefunction nobody solved on paper. Perturbation theory, first order, one contact integral — and the map of the galaxy's spiral arms rests on it.

Perturbation theory needs a solvable skeleton nearby. What if there is none — two electrons, say, shoving each other around a helium nucleus with full force? There is a second tool, and its guarantee is almost insolent. Take any normalized guess $|\phi ⟩$ whatsoever and expand it over the true (unknown!) eigenstates:

— a weighted average of energies can never undercut the smallest one. Every guess bounds the truth from above. So guessing becomes a competitive sport: build a trial state with knobs, turn the knobs to minimize $⟨\hat{H}⟩$, and the lowest value you reach is a certified ceiling on the ground state. For helium, the physics suggests the knob: each electron partly screens the nucleus from the other, so try two hydrogen-like 1s states with an effective charge:

Now be precise about what is cheap here and what is not, because the textbook version of this calculation quietly hides its only hard step. $⟨\hat{T}⟩$ and the nuclear attraction $⟨1/r⟩$ are one-dimensional radial integrals over a function we already know. The electron-electron repulsion $⟨1/r_{12}⟩$ is not: it is a six-dimensional integral over both electrons at once, and $5Z/8$ is its answer — a result every text quotes and almost none derives on the page. The trick that makes it tractable is that both one-electron densities are spherical, so of the whole multipole expansion of $1/r_{12}$ only the monopole survives the angular integrals and $1/r_{12}$ collapses to $1/\max (r_1,r_2)$, leaving a nested pair of radial integrals. The lab does exactly that, numerically, and refereeing the result against $5Z/8$ is the check that the collapse was done right. Minimizing the closed form instead would have refereed nothing but the minimizer's arithmetic.

And note which referee is doing the real work. That the minimum lands on $Z=27/16$ is a statement about the numerics. That the resulting energy sits above the measured −2.9037 Hartree is a statement about the physics, and it is the one claim in this section that nothing can make true by construction: get any of the three integrals wrong in the direction that lowers the energy and the bound simply breaks. We deliberately broke the repulsion kernel once to check, and it did.

That insolent guarantee — every guess a certified ceiling, better knobs a better ceiling — is the engine room of computational chemistry. Let the knobs be the shapes of molecular orbitals and minimizing $⟨\hat{H}⟩$ is the Hartree–Fock method; let the knob be the electron density itself and it becomes density-functional theory, the workhorse that prices the energy of essentially every molecule and material simulated today. The same inequality now runs on quantum hardware, where a variational quantum eigensolver turns the knobs on a real qubit register. One theorem, and an entire industry of approximation stands on it — the reader who wants to feel it should add a second knob to the lab's helium (let the two electrons see slightly different charges) and watch two percent become one.

One tool remains, for nuisances that oscillate — a light wave washing over an atom. First-order time-dependent perturbation theory asks: starting in $|i⟩$, what amplitude leaks into $|f⟩$ under ${\hat{H}}^{\prime }(t)=\hat{V}{e}^{-i\omega t}$? Integrating the Schrödinger equation once (the same move as 12.2, with time playing λ's role) gives an amplitude concentrated where the drive matches the gap — and you already know this answer exactly: it is the small-time, weak-drive limit of Chapter 5's Rabi formula, $P_{i\to f}\approx \frac{|V_{fi}{|}^2}{{\hslash }^2}{t}^2$ on resonance. The new physics appears when the final state is not one level but a continuum (an ionized electron, an emitted photon's many modes): summing the sharpening resonance over final states turns $t^2$ growth into steady $t$ growth — a constant rate:

Matrix element squared, times how many places there are to go. Every decay rate, cross-section, and transition intensity in atomic physics is this formula wearing different matrix elements — and Chapter 11's selection rules are precisely the statement of which of those matrix elements are zero. The metastable 2s state, the ten-million-year 21 cm lifetime, the brightness ratios within the Balmer lines: all golden-rule arithmetic. We will meet Γ again in Part IV, where the “continuum” is the quantized light field itself and spontaneous emission finally gets a first-principles price.

The lab Rust-QP/ch12-perturb runs four experiments. The first is the referee's showcase: perturbation theory for the anharmonic oscillator $\hat{H}={\hat{p}}^2/2+{\hat{x}}^2/2+\lambda {\hat{x}}^4$, with the PT coefficients computed from ladder matrix elements and the exact answer from an 80-state diagonalization — including a war story where the series' own fourth-order term ambushed the referee:

1// perturbation theory straight from the matrix elements
2    let coeff1 = x4[(0, 0)];
3    let coeff2: f64 = (1..N).map(|m| -x4[(m, 0)] * x4[(m, 0)] / m as f64).sum();
4
5    let exact_e0 = |lam: f64| -> f64 {
6// …
7        let eig = SymmetricEigen::new(h);
8        eig.eigenvalues.iter().cloned().fold(f64::INFINITY, f64::min)
9    };
10// …
11    // Referee: the residual E_exact - PT2 must be ~ c3 lambda^3 with the
12    // known c3 = 333/16 = 20.8125. The naive estimate resid/lambda^3 is
13    // contaminated at first order in lambda by the FOURTH-order term
14    // (c4 = -30885/128 ~ -241, so the naive number is short by roughly
15    // 241 * lambda; our first run tripped on exactly this). Richardson
16    // extrapolation in lambda kills that linear contamination. The run
17    // prints both, and the contamination is refereed as a FLOOR so the
18    // war story cannot quietly stop being true.
19    let c3_at = |lam: f64| -> f64 {
20        let resid = exact_e0(lam) - (0.5 + coeff1 * lam + coeff2 * lam * lam);
21        resid / lam.powi(3)
22    };
23    let c3_naive = c3_at(0.01);
24    let c3_est = 2.0 * c3_at(0.005) - c3_at(0.01);
25    let c3_exact = 333.0 / 16.0;

c3_est is Richardson extrapolation: accelerating a numerical estimate by combining results computed at several step sizes (here two values of λ) so that the leading error term cancels exactly, leaving a far cleaner number. It is the standard trick for squeezing a contaminated estimate — and here it is the thing that lets our referee separate the genuine third-order coefficient from the fourth-order term leaking into it. Note that c3_naive — the estimate before the fix — is kept and refereed too, as a floor: the run asserts that the naive number really is contaminated by more than 5%. A war story you stop testing is just an anecdote, and the panel below prints how badly the naive route actually misses.

The hydrogen experiments then reuse Chapter 10's Numerov solver verbatim — same seeds, same brackets — and feed shot wavefunctions into this chapter's formulas. Nothing analytic is assumed about the states; the 21 cm number rests on an extrapolated $|\psi (0){|}^2$ and five constants of nature:

1let (_, u_raw, r_max) = solve_state(1, 0);
2    let u = normalized(&u_raw, r_max);
3    // R(0) by linear extrapolation of u/r from the first few grid points
4    let g1 = u[1] / H_STEP;
5    let g2 = u[2] / (2.0 * H_STEP);
6    let r0 = 2.0 * g1 - g2; // linear in r: R(0) = 2 g(h) - g(2h)
7    let psi0_sq = r0 * r0 / (4.0 * std::f64::consts::PI);
8    let factor = psi0_sq * std::f64::consts::PI; // 1.0 for the exact 1s
9
10    // Fermi contact: nu = (4/3) g_e g_p (m_e/m_p) alpha^2 cR_inf
11    //                     * (1 + m_e/m_p)^{-3} * [numeric |psi(0)|^2 pi a^3]
12    let nu_hz = (4.0 / 3.0)
13        * G_E
14        * G_P
15        * ME_OVER_MP
16        * ALPHA
17        * ALPHA
18        * C_RYD_HZ
19        * (1.0 + ME_OVER_MP).powi(-3)
20        * factor;
21    let measured = 1420.405_751_768;
22    let nu_mhz = nu_hz / 1e6;
23// …
24        under("π|ψ_1s(0)|²a³ vs 1, extrapolated from the shot u", psi0_err, 5e-5),
25        under(
26            format!(
27                "21 cm frequency vs measured {:.4} MHz (rel)",
28                hf.nu_measured_mhz
29            ),
30            hf_rel,
31            1e-4,
32        ),

The fourth experiment is the helium bound of §12.5, and it is here that the lab does the work a textbook delegates to a footnote. All three expectation values are quadratures on the trial orbital; the $5Z/8$ and $Z^2/2$ closed forms appear only inside the referees, which is the one honest place for them:

1fn he_integrals(z: f64) -> HeIntegrals {
2    let n = (HE_R_MAX / HE_H) as usize;
3// …
4    // composite trapezoid over the whole grid
5    let trap = |f: &dyn Fn(usize) -> f64| -> f64 {
6        let mut s = 0.5 * (f(0) + f(n));
7        for i in 1..n {
8            s += f(i);
9        }
10        s * HE_H
11    };
12    let norm = trap(&|i| u[i] * u[i]);
13    // <T> = (1/2) integral u'^2 dr   (l = 0, u(0) = 0, so no surface term)
14    let t = 0.5 * trap(&|i| du[i] * du[i]) / norm;
15    // <1/r> = integral u^2 / r dr    (u^2/r ~ r, finite at the origin)
16    let v_ne = trap(&|i| if i == 0 { 0.0 } else { u[i] * u[i] / (i as f64 * HE_H) }) / norm;
17// …
18    // <1/r_12>. Both one-electron densities are spherical, so of the whole
19    // multipole expansion of 1/r_12 only the l = 0 term survives and the
20    // angular integrals collapse to 1/max(r1,r2):
21    //   <1/r12> = int P(r1) [ (1/r1) int_0^{r1} P dr2 + int_{r1}^inf P/r2 dr2 ] dr1
22    // with P = u^2 / norm. Cumulative trapezoids make the double integral O(n).
23    let p: Vec<f64> = (0..=n).map(|i| u[i] * u[i] / norm).collect();
24    let g: Vec<f64> = (0..=n)
25        .map(|i| if i == 0 { 0.0 } else { p[i] / (i as f64 * HE_H) })
26        .collect();
27    let mut lo = vec![0.0f64; n + 1]; // int_0^{r_i} P dr
28    for i in 1..=n {
29        lo[i] = lo[i - 1] + 0.5 * HE_H * (p[i - 1] + p[i]);
30    }
31    let mut hi = vec![0.0f64; n + 1]; // int_{r_i}^{R} P/r dr
32    for i in (0..n).rev() {
33        hi[i] = hi[i + 1] + 0.5 * HE_H * (g[i] + g[i + 1]);
34    }
35    let bracket = |i: usize| -> f64 {
36        if i == 0 {
37            hi[0]
38        } else {
39            lo[i] / (i as f64 * HE_H) + hi[i]
40        }
41    };
42    let v_ee = trap(&|i| p[i] * bracket(i));
43// …
44        under(
45            "⟨1/r₁₂⟩ by quadrature vs 5Z/8 at Z_min (rel)",
46            vee_rel,
47            2e-8,
48        ),
49        under(
50            "numeric E(Z) vs the closed form, worst over the Z grid (Ha)",
51            he_form_worst,
52            3e-6,
53        ),
54        under(
55            "golden-section Z_min of the NUMERIC E(Z) vs 27/16",
56            z_err,
57            2e-6,
58        ),
59        over(
60            "variational bound E_trial − E_measured (Ha) — FLOOR, must EXCEED tol",
61            he_margin,
62            0.0,
63        ),
64        under("helium E_min vs measured −2.9037 Ha (rel)", he_rel, 2e-2),

Most of the lab's referees are ceilings — something must be small. Two are floors, and their names say so: the naive λ³ contamination above, and the variational margin $E_{\text{trial}}-E_{\text{measured}}$, which must come out positive. That second one is the chapter's load-bearing check. Every other helium referee compares a number against a closed form and so, at worst, catches a coding slip; the bound catches wrong physics, because a trial energy below the true ground state is not a small error, it is an impossibility. Breaking the repulsion kernel on purpose — dropping the $r_2>r_1$ branch of the collapsed integral — drives the margin negative by nearly half a Hartree and the run aborts. That is what a referee with teeth looks like.

1. (F) Derive (12.2) and (12.3) by the coin-sorting sketched in 12.2: expand, collect powers of λ, and project the order-λ equation onto $⟨n|$ and onto $⟨m\ne n|$ separately. Where exactly does the freedom to choose $⟨n|{\psi }^{(1)}⟩=0$ come from?
2. (F) Level repulsion, exactly: for the 2×2 Hamiltonian with energies ±Δ and coupling v, expand the exact eigenvalues $\pm \sqrt{{\Delta }^2+v^2}$ for small v and check against (12.3). Then take Δ → 0 and check against the degenerate recipe. One matrix, both regimes.
3. (C) In the explorer, follow level 1 and find the λ at which PT2's error first exceeds 10% of the gap to level 2. Compare with the rough criterion λ|H'₁₂| ≈ gap. Then make levels 1–2 degenerate and explain why the two gold lines are straight while every exact curve bends.
4. (F) Compute $⟨\hat{\vec{L}}\cdot \hat{\vec{S}}⟩$ for general $j=l\pm \frac{1}{2}$ and show the fine-structure splitting of any level scales as ${\alpha }^2|E_n|$ — a millionth of the Bohr energy. Why does the same α² govern both the relativistic and spin-orbit corrections? (Both are v²/c² effects, and v/c ~ α in hydrogen — check it.)
5. (F, hard) The linear Stark effect. Apply a field εz to hydrogen's n = 2 level: four degenerate states {2s, 2p₀, 2p₊₁, 2p₋₁}. Use Chapter 11's selection rules (parity and m) to show the only nonzero matrix element is $⟨2s|z|2p_0⟩=-3a$ (compute it with Chapter 10's wavefunctions), diagonalize the block, and find states shifting linearly — ±3aε — with field. Every other atom shifts quadratically. What is special about hydrogen? (The trench coat again: only exact l-degeneracy lets opposite-parity states blend for free.)
6. (P) Divergence, discovered: extend the lab's Richardson trick to extract c₄ from the exact E₀(λ) data (expect ≈ −241). Then look up Dyson's 1952 argument — for λ < 0 the potential is unbounded below, so E₀(λ) cannot be analytic at λ = 0 — and reconcile: how can a series with zero radius of convergence predict energies to five digits? (Asymptotic series: stop at the smallest term.)
7. (P, hard) The deuterium radio line. Deuterium's nucleus has spin I = 1 and g_d = 0.857438: redo the contact calculation (the $\hat{\vec{S}}\cdot \hat{\vec{I}}$ jump between F = 3/2 and F = 1/2 is now 3ħ²/2 — Chapter 11's machinery for 1 ⊗ ½) and predict the DI line near 327 MHz. Radio astronomers found it in 2007, fifty years after the search began; your number should land within a fraction of a percent of 327.384 MHz.