The Hydrogen Atom

Pass a current through hydrogen gas and it glows pink. Send that glow through a prism and the pink resolves into something far stranger than a rainbow: four razor-sharp lines — red, blue-green, and two violets — and blackness in between. Every element does this, each with its own barcode, and through the whole nineteenth century nobody had the faintest idea why.

The apparatus that reveals the barcode is simple, and it is worth picturing before we read anything off it. Excite the gas with a current so it glows; send the glow through a narrow slit to make a clean pencil of light; pass that through a prism (or a diffraction grating), which bends each wavelength by a different amount; and catch the result on a detector. A continuous source would smear into a full rainbow — but hydrogen paints only a handful of sharp, separated lines, because it emits only a handful of sharp, separated wavelengths.

In 1885 Johann Balmer — a Swiss schoolteacher, sixty years old, fitting numbers as a hobby — noticed that hydrogen's four visible wavelengths all fit one absurdly simple rule: each is proportional to $n^2/(n^2-4)$ for $n=3,4,5,6$. Not roughly. To the precision of the measurements. Integers — counting numbers — hiding in the light of a gas. It took physics forty more years to understand what the integers were counting, and you already know the answer from Chapter 5: spectral lines are beat notes between energy levels, so Balmer had unknowingly measured the level ladder itself. Here is the whole structure his formula was pointing at:

The task of this chapter is to earn that ladder — to start from one proton, one electron, and the Coulomb attraction, and have $-13.6\text{eV}/n^2$ forced on us the way the box's $n^2$ and the oscillator's $n+\frac{1}{2}$ were: by a wave that has to fit.

Two bodies, not one: the electron does not orbit a nailed-down proton; both swing about their common center of mass, like a waltzing couple. Exactly as in classical mechanics, the two-body problem splits into free center-of-mass motion plus a one-body problem for the relative coordinate, with the reduced mass

A twentieth of a percent — file it away; it is about to buy us two decimal places of glory. From here on we use atomic units ($\hslash =m_e=e=1$: lengths in Bohr radii, energies in Hartrees) and the problem is a single particle in the central potential $V(r)=-1/r$. “Central” is the magic word: the Hamiltonian commutes with every rotation, so by Chapter 9 we may hunt for states of definite $l$ and $m$, and the angular dependence is already solved:

where the substitution $u=rR$ is the little gift that keeps the chapter honest: it turns the messy three-dimensional radial operator into a plain one-dimensional Schrödinger equation — Chapter 6 on the half-line:

Look at the effective potential. The Coulomb well pulls in; the centrifugal barrier $l(l+1)/2r^2$ — the quantum cost of carrying angular momentum — pushes out, and wins at small r. An electron with $l>0$ is walled away from the nucleus: its wavefunction must vanish there like $\psi \sim r^l$, so it is only the combination $u=rR$ that the barrier forces to go as $r^{l+1}$ — which is why the seed in §10.5's solver carries that extra power. An s-electron, with no barrier, presses right up against the proton (its $\psi$ is flatly finite at $r=0$). Chemistry will spend that difference for the rest of the book.

Now the squeeze, the same one that quantized the box and the well. Far away, (10.3) says $u^{\prime \prime }\approx -2Eu$; for a bound state ($E<0$) the solutions are $e^{\pm \kappa r}$ with $\kappa =\sqrt{-2E}$, and physics keeps only the decay. Near the origin the centrifugal term rules and forces $u\sim r^{l+1}$. Write the solution as decay × regularity × a polynomial correction, feed it back into (10.3), and the recursion for the polynomial has a fatal flaw: unless it terminates, its tail rebuilds $e^{+2\kappa r}$ and the wave explodes (exercise 4 walks every step). Termination happens only when

There is Balmer's integer, caught red-handed: it counts the allowed decay rates of a wave pinned between “vanish at the nucleus” and “vanish at infinity.” The natural length that falls out, $a={\hslash }^2/m_e{e}^2=0.0529\text{nm}$ — the Bohr radius — sets the size of every atom in the universe, and the ground state is the simplest wavefunction in this book:

No orbit, no plane, no ellipse — a fuzzy exponential ball. And one oddity to savor before the lab confirms it numerically: equation (10.4) contains $n$ but not $l$. The 2s state (no angular momentum, one radial node) and the 2p state (one unit of angular momentum, no radial node) are completely different waves with identical energy. For a general central potential that never happens — each $l$ gets its own ladder. The Coulomb potential alone hides an extra conserved quantity (the Runge–Lenz vector; exercise 5 comes at the same symmetry from the Hellmann–Feynman side) that forces the ladders into registry, giving level $n$ its full ${\sum }_{l<n}(2l+1)=n^2$ degeneracy. Physicists call it “accidental.” It is nothing of the sort — it is a symmetry wearing a trench coat, and the periodic table is what happens when multi-electron screening tears the coat off.

Time to look at the things. First the radial story alone — where, at what distance, the electron lives. The honest density is not $|R{|}^2$ but $P(r)=r^2|R{|}^2$ (a thin shell at radius r has area $4\pi r^2$ — more room farther out), and that geometric factor is what makes even the ground state peak away from the nucleus:

Now marry the radial factor to Chapter 9's spherical harmonics and meet the complete object — the orbital, chemistry's alphabet:

One bookkeeping note with large consequences: at each n there are $n^2$ of these states, and (once Part IV gives the electron its spin and its exclusion rule) each holds two electrons, so shell n has capacity $2n^2=2,8,18,32$. Be careful with the next step, because it is where the folklore goes wrong: those are shell capacities, not the lengths of the periodic table's rows. The rows run 2, 8, 8, 18, 18, 32, 32 — each of the last three repeated — because real atoms do not fill a shell before starting the next. Orbitals fill in order of increasing $n+l$ (and, within a tie, increasing $n$), which puts 4s below 3d and so closes a row at argon after 3p rather than after 3d. That reordering is exactly what happens when screening tears the trench coat off, and exercise 7 has you break the degeneracy yourself and watch the ordering appear. The shelving is built here; the filling rule is not ours yet.

The lab Rust-QP/ch10-hydrogen never lets equation (10.4) answer a question. It treats (10.3) as Chapter 6 treated its wells: a boundary-value problem, to be shot. One honest caveat, since the lab admits it in a comment: the search bracket for each level is placed using $-1/2n^2$, because you have to start the hunt somewhere. But the bracket is wide — it runs to the midpoints of the neighbouring levels — and 64 halvings inside it put the answer wherever the differential equation wants it, to about one part in $10^{19}$ of the bracket width. If the integrator were wrong the root would move and the first referee would say so. Integrate outward from the origin with the Numerov method — a classic scheme for second-order equations of the form $u^{\prime \prime }=f(r)u$ that, by leaning on the second derivative it is already computing, reaches $O(h^4)$ accuracy for barely more work than a crude Euler step. A wrong energy explodes at large r, positive on one side of an eigenvalue and negative on the other; bisect 64 times on that sign. The Bohr formula is not an input — it is the verdict:

1/// Numerov integration of u'' = f(r) u outward. Returns (u_grid, r_grid).
2fn shoot(l: u32, e: f64, r_max: f64) -> Vec<f64> {
3    let n_steps = (r_max / H) as usize;
4    let f = |r: f64| -> f64 {
5        (l * (l + 1)) as f64 / (r * r) - 2.0 / r - 2.0 * e
6    };
7    let mut u = vec![0.0f64; n_steps + 1];
8    // seed near the origin with TWO terms of the exact series
9    // u ~ r^{l+1} (1 - r/(l+1) + ...): the second term matters — with the
10    // bare power law the eigenvalues came out 3x worse than tolerance
11    let seed = |r: f64| r.powi(l as i32 + 1) * (1.0 - r / (l + 1) as f64);
12    u[1] = seed(H);
13    u[2] = seed(2.0 * H);
14    let h2_12 = H * H / 12.0;
15    for i in 2..n_steps {
16        let r_m = i as f64 * H - H;
17        let r_0 = i as f64 * H;
18        let r_p = i as f64 * H + H;
19        let num = 2.0 * u[i] * (1.0 + 5.0 * h2_12 * f(r_0)) - u[i - 1] * (1.0 - h2_12 * f(r_m));
20        u[i + 1] = num / (1.0 - h2_12 * f(r_p));
21        // renormalize mid-flight if the trial solution explodes, so we
22        // never overflow — only the SIGN at r_max matters for bisection
23        if u[i + 1].abs() > 1e250 {
24            let s = 1e-100;
25            for v in u[..=i + 1].iter_mut() {
26                *v *= s;
27            }
28        }
29    }
30    u
31}

Two war stories this time, both instructive. The origin seed $u\approx r^{l+1}$ alone left the energies hundreds of times worse than the run now demands; adding the second series term (one line) bought a factor of about 300. And the first bracketing scheme (“±35% around the guess”) worked until n = 5, where hydrogen's crowding levels put two eigenvalues inside one bracket and the bisection lost its sign change — the fix is to bracket by midpoints toward the neighboring levels, which contains exactly one eigenvalue by construction:

1/// Bisect on the sign of u(r_max) inside a bracket around the expected
2/// energy window. The window seeds the SEARCH; the ODE decides the answer.
3/// Bracket = midpoints toward the neighboring levels, so it contains
4/// exactly one eigenvalue (the levels crowd together as n grows — a
5/// naive +/-35% window would swallow two at n = 5).
6fn solve_state(n: u32, l: u32) -> Solved {
7    let e_of = |k: u32| -0.5 / (k * k) as f64;
8    let e_below = if n > 1 { e_of(n - 1) } else { 1.5 * e_of(1) };
9    let e_above = e_of(n + 1);
10    let (mut e_lo, mut e_hi) = (
11        0.5 * (e_of(n) + e_below),
12        0.5 * (e_of(n) + e_above),
13    );
14    let r_max = 28.0 * n as f64;
15    let tail_sign = |e: f64| -> f64 {
16        let u = shoot(l, e, r_max);
17        u[u.len() - 1].signum()
18    };
19    let (s_lo, s_hi) = (tail_sign(e_lo), tail_sign(e_hi));
20    assert!(
21        s_lo != s_hi,
22        "no sign change in bracket for n={n}, l={l} — bad window"
23    );
24    for _ in 0..64 {
25        let e_mid = 0.5 * (e_lo + e_hi);
26        if tail_sign(e_mid) == s_lo {
27            e_lo = e_mid;
28        } else {
29            e_hi = e_mid;
30        }
31    }
32    let e = 0.5 * (e_lo + e_hi);
33    let u = shoot(l, e, r_max);
34    Solved { e, u, r_max }
35}

Then the scoreboard. Six referees, each one a comparison between two things that could genuinely disagree: the shot eigenvalues against $-1/2n^2$; the counted nodes of $u$ against n−l−1; the two different differential equations for 2s and 2p landing on one root; the virial theorem, which compares an integral over the numeric wavefunction to the bisected energy; $⟨r⟩$ against its closed form; and the telescope test. Every tolerance is the achieved value rounded up by about a quarter, and none of them is decorative — perturb the Coulomb strength by one percent and five of the six fail at once. Note the finiteness pass that runs first: a NaN compares false against every tolerance, so it would sneak past a carelessly written check rather than announce itself.

1// -- the referees ------------------------------------------------------
2    // Meta-referee: FINITE before small. NaN compares false against every
3    // tolerance, so an unnoticed 0/0 would read as a silent pass in a
4    // `!(x > tol)` world rather than the loud failure it deserves.
5    for (label, v) in [
6        ("worst_energy_err", referee.worst_energy_err),
7        // …
8    ] {
9        assert!(v.is_finite(), "{label} is not finite: {v}");
10    }
11    // …
12    let referees = vec![
13        Referee {
14            name: format!("worst |E_shot − (−1/2n²)| over {n_states} states, Hartree"),
15            value: referee.worst_energy_err,
16            tol: 6e-8,
17            pass: referee.worst_energy_err < 6e-8,
18        },
19        // …
20        Referee {
21            name: format!(
22                "telescope test: worst |λ − λ_measured| over {} Balmer lines, nm",
23                lines.len()
24            ),
25            value: referee.worst_lambda_err_nm,
26            tol: 0.012,
27            pass: referee.worst_lambda_err_nm < 0.012,
28        },
29    ];
30    // …
31    let mut all_passed = true;
32    for r in &data.referees {
33        // …
34        all_passed &= r.pass;
35    }
36    assert!(all_passed, "a referee failed");

Sit with what just happened, because it reaches far past hydrogen. A boundary-value problem on a laptop reproduced light measured in a discharge tube to five figures — which means the logic runs backwards too. Point a spectrograph at a star, a nebula, or a quasar, read the pattern of Balmer lines out of its light, and you are reading the atoms that made it: their identity from which lines appear, their velocity from how far the whole barcode is Doppler-shifted, the expansion of the universe from the redshift of hydrogen in the most distant galaxies. Every composition, temperature, and recession speed in the observable cosmos is inferred from equation (10.4) run in reverse. This is the founding instrument of astrophysics.

And look once more at the residual — the last digits where $E_n=-13.606\text{eV}/n^2$ and the real atom quietly part ways. Run the energy referee's bound through $\lambda =hc/\Delta E$ and the solver can only be blamed for a small fraction of it; the rest of the gap sits at the scale where the next layer of physics we have not put in yet lives. The electron's spin and its relativistic speed split each level into fine structure; the proton's own magnetism adds the hyperfine splitting that gives the 21 cm line radio astronomers map the galaxy with; even the vacuum itself nudges the levels (the Lamb shift). Each is a thread the next chapters pull. The clean $1/n^2$ ladder was never the end of the story — it is the baseline precise enough to make the next layer of the atom visible.

1. (F) Do the u = rR substitution honestly: start from the full 3D Laplacian's radial part $\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)$ and show every first-derivative term cancels, leaving (10.3). Where exactly does the boundary condition $u(0)=0$ come from, and why would $R(0)=\text{finite},u(0)\ne 0$ be physically fatal?
2. (F) Prove the virial theorem the lab checks: for $V\propto r^k$, show $2⟨T⟩=k⟨V⟩$ using equation (5.4) on the operator $\hat{\vec{r}}\cdot \hat{\vec{p}}$ in a stationary state. For Coulomb (k = −1) conclude $⟨V⟩=2E$ — so a bound electron speeds up as it binds deeper, paying with double the potential.
3. (C) In the radial explorer, verify by eye that P(r) for 1s peaks at exactly r = 1 (prove it: maximize $r^2{e}^{-2r}$). Then find the 2p peak and show it sits at r = 4 — while ⟨r⟩ = 5. Why is the mean beyond the peak? (Look at the tail.)
4. (F, hard) The full termination argument: substitute $u=r^{l+1}{e}^{-\kappa r}v(r)$ into (10.3), derive the recursion $\frac{c_{k+1}}{c_k}=\frac{2\kappa (k+l+1)-2}{(k+1)(k+2l+2)}$ for the series $v=\sum c_k{r}^k$, show that a non-terminating series behaves like $e^{2\kappa r}$ at large r (compare the ratio to that exponential's), and conclude κ = 1/n. You have just done, in full, what the Numerov solver does blind.
5. (F, hard) The trench coat, via Hellmann–Feynman: for a Hamiltonian depending on a parameter, $\partial E/\partial \lambda =⟨\partial \hat{H}/\partial \lambda ⟩$ (prove it — two lines with normalization). Apply it twice to hydrogen: differentiate with respect to the nuclear charge Z to get $⟨1/r⟩=Z/n^2$, then treat $l$ as the parameter (it enters V_eff continuously!) to extract $⟨1/r^2⟩$. Check both against the lab's wavefunctions numerically.
6. (P) Discover deuterium. Urey found it (1931) as a faint Balmer line displaced from hydrogen's — the nucleus is twice as heavy, so the reduced mass shifts. Change one constant in the lab (deuteron mass ≈ 3670.48 mₑ), recompute Hα, and predict the H–D splitting. (You should land near 0.18 nm — resolvable on a good grating, and worth a Nobel prize if you'd done it first.)
7. (P, hard) Break the trench coat yourself: add screening to the lab's potential, $V=-(1+b{e}^{-r})/r$ with b = 1 (a crude sodium: one valence electron outside a screening core). Re-shoot the “3s” and “3p” states and watch the degeneracy shatter — the s state, which penetrates the core, drops far below p. Quantify it as a quantum defect $E=-1/2(n-{\delta }_l{)}^2$ and compare the ordering you find with the actual periodic table's 4s-before-3d anomaly.