Chapter 38: Capstone — Hydrogen Atom from First Principles

"The hydrogen atom has already been solved. Nobody knows the solution to the helium atom. But hydrogen — a single proton and a single electron — is where quantum mechanics proves itself, down to the twelfth decimal place." — Victor Weisskopf, in conversation

"The spectrum of hydrogen has proved to be the Rosetta stone of modern physics." — Edward Condon, The Theory of Atomic Spectra (1935)

You have reached the summit.

Throughout this textbook, we have returned to the hydrogen atom again and again — as a warmup (Chapter 2), as an exactly solvable system (Chapter 5), as a testbed for perturbation theory (Chapters 17-18), as a benchmark for variational methods (Chapter 19), and as a bridge to multi-electron atoms (Chapter 16) and relativistic quantum mechanics (Chapter 29). Each visit revealed new structure, new corrections, new physics. But we have never assembled the complete picture in one place.

This capstone chapter changes that. We will build the hydrogen atom from the ground up, layer by layer, using every major technique you have learned:

1. Exact solution — the Coulomb problem in all its analytical glory
2. Perturbation corrections — fine structure, hyperfine structure, and radiative corrections
3. Variational bounds — independent verification using trial wavefunctions
4. Numerical computation — direct solution of the radial Schrödinger equation on a grid

We will compare all four approaches, quantify their agreements and disagreements, and connect everything to real spectroscopic data. By the end, you will have a complete hydrogen simulation — a piece of software that predicts the hydrogen spectrum to extraordinary accuracy, built entirely from the principles developed in this book.

This is not a review chapter. It is a synthesis chapter. The whole is more than the sum of its parts.

🏃 Fast Track: This chapter is designed to be worked through sequentially; each section builds on the previous one. The code in Section 38.8 integrates everything. If you are short on time, read Sections 38.1-38.2 carefully, skim the derivations in 38.3-38.5, and focus on the comparison table in Section 38.6 and the code in Section 38.8.

38.1 The Hydrogen Atom: What We Know Now

Why Hydrogen Matters

Hydrogen is the simplest atom: one proton, one electron, one Coulomb interaction. No electron-electron repulsion. No exchange symmetry complications. No screening. Just:

$$\hat{H} = -\frac{\hbar^2}{2m_e}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 r}$$

This simplicity is precisely what makes hydrogen the most important atom in physics. It is:

• The only atom with an exact analytical solution (ignoring QED corrections)
• The benchmark against which all approximation methods are tested
• The proving ground for quantum electrodynamics — the Lamb shift discrepancy between Dirac theory and experiment drove the development of QED
• The most abundant element in the universe — its spectral lines are the primary tool of observational astronomy
• The basis for understanding all other atoms — from helium onward, atomic physics is perturbation theory around hydrogen-like systems

Let us take stock of everything we know about hydrogen, organized by the precision layer at which each effect operates.

The Energy Hierarchy

The energy levels of hydrogen form a hierarchy of corrections, each smaller than the last by roughly a factor of $\alpha^2 \approx 1/18,769$, where $\alpha = e^2/(4\pi\epsilon_0\hbar c) \approx 1/137.036$ is the fine-structure constant:

Each layer splits degeneracies that were present at the previous level. The gross structure gives $n^2$-fold degeneracy. Fine structure partially lifts this, leaving states with the same $j$ degenerate. The Lamb shift splits the $2S_{1/2}$ and $2P_{1/2}$ levels that fine structure leaves degenerate. Hyperfine structure splits each $j$ level into $F$ sublevels.

💡 Key Insight: The fine-structure constant $\alpha$ is the organizing parameter of atomic physics. It governs the ratio of each correction layer to the one above it. That $\alpha \approx 1/137$ is small is what makes perturbation theory work for atoms. If $\alpha$ were of order unity, atoms would be relativistic and the perturbative hierarchy would collapse.

Quantum Numbers: The Complete Set

A hydrogen eigenstate is fully specified by the quantum numbers:

$$|n, l, j, m_j, F, m_F\rangle$$

where:

In spectroscopic notation, we write $n^{2S+1}L_J$ — for example, the ground state is $1^2S_{1/2}$, and the states involved in the Balmer-$\alpha$ transition are $3^2S_{1/2}$, $3^2P_{1/2}$, $3^2P_{3/2}$, $3^2D_{3/2}$, $3^2D_{5/2}$ → $2^2S_{1/2}$, $2^2P_{1/2}$, $2^2P_{3/2}$.

🔗 Connection: The quantum numbers $n, l, m_l$ were introduced in Chapter 5 via separation of variables. The spin quantum number $m_s$ appeared in Chapter 13. The coupled representation $(j, m_j)$ was constructed in Chapter 14 via addition of angular momentum. The hyperfine quantum numbers $(F, m_F)$ were introduced in Chapter 18. Here we use them all simultaneously for the first time.

38.2 Exact Solution: The Coulomb Problem Revisited

Setting Up the Problem

We begin from the time-independent Schrödinger equation for hydrogen. Using the reduced mass $\mu = m_e m_p/(m_e + m_p) \approx m_e(1 - m_e/m_p)$ to account for the finite proton mass:

$$\left[-\frac{\hbar^2}{2\mu}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 r}\right]\psi(\mathbf{r}) = E\psi(\mathbf{r})$$

In spherical coordinates, exploiting the central-force nature of the Coulomb potential, we separate:

$$\psi_{nlm}(r, \theta, \phi) = R_{nl}(r) \, Y_l^m(\theta, \phi)$$

where the $Y_l^m$ are the spherical harmonics (Chapter 12) that satisfy:

$$\hat{L}^2 Y_l^m = \hbar^2 l(l+1) Y_l^m, \quad \hat{L}_z Y_l^m = \hbar m Y_l^m$$

The Radial Equation

Substituting the separated form and dividing by $Y_l^m$, we obtain the radial equation (Chapter 5):

$$-\frac{\hbar^2}{2\mu}\left[\frac{d^2 u}{dr^2} - \frac{l(l+1)}{r^2}u\right] - \frac{e^2}{4\pi\epsilon_0 r} u = Eu$$

where $u(r) = rR(r)$ is the reduced radial wavefunction, with boundary conditions $u(0) = 0$ and $u(r) \to 0$ as $r \to \infty$ for bound states.

Solving with Dimensionless Variables

Introducing the dimensionless variable $\rho = r/a_0$ where $a_0 = 4\pi\epsilon_0\hbar^2/(\mu e^2)$ is the Bohr radius (corrected for reduced mass), and defining $\epsilon = -E/E_1$ where $E_1 = \mu e^4/(2(4\pi\epsilon_0)^2\hbar^2) = 13.6$ eV, the radial equation becomes:

$$\frac{d^2 u}{d\rho^2} = \left[\frac{l(l+1)}{\rho^2} - \frac{2}{\rho} + \frac{1}{n^2}\right]u$$

after anticipating that $\epsilon = 1/n^2$.

The Solution: Laguerre Polynomials

The complete solution (derived in full in Chapter 5, Section 5.4) is:

$$R_{nl}(r) = -\sqrt{\left(\frac{2}{na_0}\right)^3 \frac{(n-l-1)!}{2n[(n+l)!]^3}} \; e^{-r/(na_0)} \left(\frac{2r}{na_0}\right)^l L_{n-l-1}^{2l+1}\left(\frac{2r}{na_0}\right)$$

where $L_{n-l-1}^{2l+1}$ is an associated Laguerre polynomial.

⚠️ Common Misconception: There are multiple conventions for associated Laguerre polynomials in the literature. Griffiths, Sakurai, and Shankar use different normalizations. When comparing with a reference, always verify which convention is being used. Our convention follows Griffiths 3rd edition: $L_q^p(x)$ where $q$ is the degree of the polynomial (non-negative integer) and $p$ is the order.

Energy Eigenvalues

The bound-state energies depend only on $n$:

$$E_n = -\frac{\mu e^4}{2(4\pi\epsilon_0)^2\hbar^2} \cdot \frac{1}{n^2} = -\frac{13.6 \text{ eV}}{n^2}$$

This is the famous $1/n^2$ spectrum, with the $n^2$-fold degeneracy in $l$ and $m$ (and an additional factor of 2 from spin):

$$g_n = 2n^2$$

The degeneracy in $l$ is specific to the $1/r$ potential — it reflects a hidden $SO(4)$ symmetry (the Runge-Lenz vector is a conserved quantity). This "accidental" degeneracy is broken by every correction we will consider.

The Hidden Symmetry: The Runge-Lenz Vector

The $n^2$-fold degeneracy of the hydrogen spectrum — the fact that states with different $l$ but the same $n$ have the same energy — is not explained by the obvious $SO(3)$ rotational symmetry of the Coulomb potential. Rotational symmetry explains the $m$-degeneracy (states with different $m_l$ have the same energy), but not the $l$-degeneracy.

The deeper symmetry is generated by the Runge-Lenz vector:

$$\hat{\mathbf{A}} = \frac{1}{2m_e}(\hat{\mathbf{p}} \times \hat{\mathbf{L}} - \hat{\mathbf{L}} \times \hat{\mathbf{p}}) - \frac{e^2}{4\pi\epsilon_0}\frac{\hat{\mathbf{r}}}{r}$$

This vector operator commutes with the hydrogen Hamiltonian: $[\hat{H}, \hat{\mathbf{A}}] = 0$. Together, $\hat{\mathbf{L}}$ and $\hat{\mathbf{A}}$ generate an $SO(4)$ symmetry group (for bound states), which is isomorphic to $SU(2) \times SU(2)$. The irreducible representations of this group are labeled by a single quantum number $n$, and each representation has dimension $n^2$ — exactly the observed degeneracy.

The classical counterpart of this symmetry is the fact that Keplerian orbits do not precess. In a general central potential, the elliptical orbit precesses slowly (the perihelion advances, as Mercury's does in general relativity). But for the exact $1/r$ potential, the orbit closes perfectly after each revolution. The Runge-Lenz vector points from the focus to the perihelion and is constant — this is the classical $SO(4)$ symmetry.

💡 Key Insight: The $l$-degeneracy of hydrogen is not an accident of the mathematics — it reflects a deep symmetry of the Coulomb potential. When this symmetry is broken (by fine structure, by screening in multi-electron atoms, by external fields), the $l$-degeneracy lifts, and states with different $l$ acquire different energies. Understanding why the degeneracy existed in the first place makes the pattern of its breaking more transparent.

Explicit Wavefunctions

Let us write out the first few radial wavefunctions explicitly, as we will need them for perturbation theory:

$n = 1, l = 0$: $$R_{10}(r) = 2\left(\frac{1}{a_0}\right)^{3/2} e^{-r/a_0}$$

$n = 2, l = 0$: $$R_{20}(r) = \frac{1}{2\sqrt{2}}\left(\frac{1}{a_0}\right)^{3/2}\left(2 - \frac{r}{a_0}\right)e^{-r/(2a_0)}$$

$n = 2, l = 1$: $$R_{21}(r) = \frac{1}{2\sqrt{6}}\left(\frac{1}{a_0}\right)^{3/2}\frac{r}{a_0}e^{-r/(2a_0)}$$

$n = 3, l = 0$: $$R_{30}(r) = \frac{2}{3\sqrt{3}}\left(\frac{1}{a_0}\right)^{3/2}\left(1 - \frac{2r}{3a_0} + \frac{2r^2}{27a_0^2}\right)e^{-r/(3a_0)}$$

These wavefunctions encode all the physics of the Coulomb problem: the exponential decay (set by $n$), the polynomial structure (encoding the number of radial nodes $n - l - 1$), and the centrifugal barrier ($r^l$ behavior near the origin).

✅ Checkpoint: Verify that $R_{10}(r)$ is normalized: $\int_0^\infty |R_{10}(r)|^2 r^2 \, dr = 1$. You should find that the integral evaluates to $4/a_0^3 \times a_0^3/4 = 1$. If you cannot do this integral by hand, you need to review integration by parts with exponentials before proceeding.

Key Expectation Values

For the exact Coulomb wavefunctions, several important expectation values have closed-form expressions:

$$\langle r \rangle_{nl} = \frac{a_0}{2}\left[3n^2 - l(l+1)\right]$$

$$\langle r^2 \rangle_{nl} = \frac{a_0^2 n^2}{2}\left[5n^2 + 1 - 3l(l+1)\right]$$

$$\left\langle \frac{1}{r} \right\rangle_{nl} = \frac{1}{n^2 a_0}$$

$$\left\langle \frac{1}{r^2} \right\rangle_{nl} = \frac{1}{n^3 a_0^2(l + 1/2)}$$

$$\left\langle \frac{1}{r^3} \right\rangle_{nl} = \frac{1}{n^3 a_0^3 l(l+1/2)(l+1)} \quad (l \neq 0)$$

The Kramers relation, $\langle r^{s+1}\rangle = \frac{(2s+1)a_0}{s+1}\langle r^s\rangle - \frac{s}{4(s+1)}[(2l+1)^2 - s^2]a_0^2\langle r^{s-1}\rangle$, provides a powerful recursion that generates all moments $\langle r^s \rangle$ from the known values above.

📊 By the Numbers: For the ground state ($n = 1, l = 0$): $\langle r \rangle = \frac{3}{2}a_0 = 0.794$ Å, $\Delta r = \sqrt{\langle r^2\rangle - \langle r\rangle^2} = \frac{\sqrt{3}}{2}a_0 = 0.459$ Å. The electron "lives" at about $1.5$ Bohr radii from the proton, with a spread comparable to the mean distance — the hydrogen atom is intrinsically fuzzy.

38.3 Perturbation Corrections: Fine Structure

The exact Coulomb solution treats the electron non-relativistically and ignores spin. Real hydrogen has corrections that arise from special relativity and the electron's magnetic moment. These fine-structure corrections were the subject of Chapter 18; here we assemble them into a complete picture.

The Three Fine-Structure Terms

The fine-structure Hamiltonian consists of three terms, each of order $\alpha^2 E_n$:

$$\hat{H}_{\text{fs}} = \hat{H}_{\text{rel}} + \hat{H}_{\text{SO}} + \hat{H}_{\text{Darwin}}$$

1. Relativistic kinetic energy correction:

The non-relativistic kinetic energy $p^2/(2m_e)$ is the first term in the expansion of the relativistic kinetic energy $\sqrt{p^2c^2 + m_e^2c^4} - m_ec^2$. The next term gives:

$$\hat{H}_{\text{rel}} = -\frac{\hat{p}^4}{8m_e^3 c^2}$$

Using first-order perturbation theory (Chapter 17):

$$E_{\text{rel}}^{(1)} = \langle nlm | \hat{H}_{\text{rel}} | nlm \rangle = -\frac{1}{2m_ec^2}\left[\left(\frac{E_n}{1}\right)^2 + 2E_n \frac{e^2}{4\pi\epsilon_0}\left\langle\frac{1}{r}\right\rangle + \left(\frac{e^2}{4\pi\epsilon_0}\right)^2\left\langle\frac{1}{r^2}\right\rangle\right]$$

After substituting the expectation values:

$$E_{\text{rel}}^{(1)} = -\frac{E_n^2}{2m_ec^2}\left(\frac{4n}{l+1/2} - 3\right)$$

2. Spin-orbit coupling:

The electron, moving through the proton's electric field, experiences a magnetic field in its rest frame. This field couples to the electron's magnetic moment:

$$\hat{H}_{\text{SO}} = \frac{e^2}{8\pi\epsilon_0 m_e^2 c^2 r^3}\hat{\mathbf{L}} \cdot \hat{\mathbf{S}}$$

The factor here includes the Thomas precession correction (a factor of 1/2). In the coupled basis $|n, l, j, m_j\rangle$ (Chapter 14):

$$\langle \hat{\mathbf{L}} \cdot \hat{\mathbf{S}} \rangle = \frac{\hbar^2}{2}[j(j+1) - l(l+1) - 3/4]$$

giving:

$$E_{\text{SO}}^{(1)} = \frac{E_n^2}{m_ec^2} \cdot \frac{n[j(j+1) - l(l+1) - 3/4]}{l(l+1/2)(l+1)} \quad (l \neq 0)$$

3. Darwin term:

This is a relativistic correction to the potential energy, arising from the Zitterbewegung — the electron's "trembling motion" that smears out the point-like Coulomb interaction:

$$\hat{H}_{\text{Darwin}} = \frac{\pi\hbar^2 e^2}{2m_e^2 c^2 (4\pi\epsilon_0)} \delta^3(\mathbf{r})$$

This term is nonzero only for $l = 0$ states (which have $|\psi(0)|^2 \neq 0$):

$$E_{\text{Darwin}}^{(1)} = \frac{2E_n^2}{m_ec^2} \cdot n \quad (l = 0 \text{ only})$$

The Complete Fine-Structure Formula

Remarkably, when all three terms are combined, the result depends only on $n$ and $j$ — not on $l$ separately:

$$E_{\text{fs}}^{(1)} = \frac{\alpha^2 E_n}{n^2}\left(\frac{n}{j + 1/2} - \frac{3}{4}\right)$$

or equivalently:

$$E_{nj} = E_n\left[1 + \frac{\alpha^2}{n^2}\left(\frac{n}{j+1/2} - \frac{3}{4}\right)\right]$$

💡 Key Insight: The remarkable fact that the fine-structure correction depends only on $j$ (not on $l$) is not a coincidence — it reflects the hidden symmetry of the Dirac equation. The Dirac equation (Chapter 29) gives this same result exactly, not as a perturbative correction. The $l$-independence of the fine-structure formula is a consequence of the exact degeneracy in $j$ that the Dirac equation predicts. This degeneracy is broken by the Lamb shift (QED corrections), which we discuss in Section 38.3.3.

Numerical Fine-Structure Corrections

Let us compute the fine-structure splittings for the $n = 2$ level:

For $n = 2$, $E_2 = -3.40$ eV:

The fine-structure splitting between $2P_{3/2}$ and $2P_{1/2}$ (equivalently $2S_{1/2}$) is:

$$\Delta E_{\text{fs}}(n=2) = \frac{\alpha^2 |E_2|}{16} = 4.53 \times 10^{-5} \text{ eV} = 0.365 \text{ cm}^{-1} = 10.9 \text{ GHz}$$

This matches experiment to within the Lamb shift correction.

The Lamb Shift

The Dirac equation predicts that the $2S_{1/2}$ and $2P_{1/2}$ states are exactly degenerate (same $n$ and $j$). In 1947, Willis Lamb and Robert Retherford discovered experimentally that the $2S_{1/2}$ level lies above the $2P_{1/2}$ level by about 1057 MHz (0.035 cm$^{-1}$, or $4.4 \times 10^{-6}$ eV).

This tiny discrepancy — about 1 part in $10^6$ of the $n = 2$ energy — was the clue that drove the development of quantum electrodynamics. The Lamb shift arises from:

1. Electron self-energy: The electron interacts with its own virtual photon cloud, effectively smearing its charge over a distance $\sim \alpha a_0 \ln(1/\alpha)$. This smearing modifies the Coulomb potential experienced by $s$-states (which have nonzero probability at the origin) more than $p$-states.

2. Vacuum polarization: Virtual electron-positron pairs in the vacuum are polarized by the proton's field, partially screening the charge. This effect shifts $s$-states down and partially counteracts the self-energy shift.

The dominant contribution to the Lamb shift for the $2S_{1/2}$ state is:

$$\Delta E_{\text{Lamb}}(2S_{1/2}) \approx \frac{\alpha^5 m_e c^2}{6\pi}\left[\ln\frac{1}{\alpha^2} - \ln k_0(2,0) + \frac{19}{30}\right]$$

where $k_0(2,0)$ is a Bethe logarithm that must be computed numerically. The full calculation gives:

$$\Delta E_{\text{Lamb}}(2S_{1/2} - 2P_{1/2}) = 1057.845(9) \text{ MHz}$$

in agreement with the experimental value $1057.845(3)$ MHz.

🔵 Historical Note: The measurement of the Lamb shift was announced at the famous Shelter Island Conference in June 1947. Hans Bethe, on the train ride from Shelter Island to Schenectady, performed the first theoretical calculation of the Lamb shift using non-relativistic QED, obtaining approximately 1040 MHz — within 4% of the experimental value. Bethe later called this calculation "the most important thing I ever did." Julian Schwinger, Richard Feynman, and Sin-Itiro Tomonaga independently developed the full relativistic QED treatment over the next two years, for which they shared the 1965 Nobel Prize.

Hyperfine Structure

The proton has a magnetic moment $\boldsymbol{\mu}_p = g_p \mu_N \hat{\mathbf{I}}/\hbar$ where $g_p = 5.586$, $\mu_N = e\hbar/(2m_p) = 5.051 \times 10^{-27}$ J/T is the nuclear magneton, and $\hat{\mathbf{I}}$ is the proton spin operator ($I = 1/2$).

The interaction of the electron's magnetic moment with the proton's magnetic dipole field gives the hyperfine Hamiltonian. For $s$-states:

$$E_{\text{hf}}^{(1)} = \frac{4}{3} g_p \alpha^2 \frac{m_e}{m_p} E_n \cdot \frac{1}{n^3} \times \begin{cases} +1/4 & F = 1 \\ -3/4 & F = 0 \end{cases}$$

The hyperfine splitting of the ground state ($n = 1, l = 0, j = 1/2$) between $F = 1$ and $F = 0$ is:

$$\Delta E_{\text{hf}} = \frac{4}{3}g_p \alpha^2 \frac{m_e}{m_p} \frac{|E_1|}{1} = 5.88 \times 10^{-6} \text{ eV}$$

This corresponds to a frequency of $\nu = 1420.405751768(1)$ MHz, or a wavelength of 21.1 cm.

🧪 Experiment: The hydrogen 21-cm line is one of the most important spectral lines in all of astronomy. It is used to map the distribution of neutral hydrogen throughout the universe, including the spiral structure of our own Milky Way galaxy. The frequency is so precisely known that it has been proposed as a universal standard for interstellar communication — the famous "hydrogen maser" on the Pioneer plaque and Voyager golden record.

38.4 Variational Bounds

The variational principle (Chapter 19) provides a completely independent way to bound the hydrogen ground-state energy. This serves as a powerful check on the exact solution and demonstrates the method for systems where exact solutions do not exist.

The Variational Principle: Quick Review

For any normalized trial wavefunction $|\psi_{\text{trial}}\rangle$:

$$E_{\text{ground}} \leq \langle \psi_{\text{trial}} | \hat{H} | \psi_{\text{trial}} \rangle$$

The closer $|\psi_{\text{trial}}\rangle$ is to the true ground state, the tighter the bound. Equality holds if and only if $|\psi_{\text{trial}}\rangle$ is the exact ground state.

Trial Wavefunction 1: Simple Exponential

The simplest reasonable trial wavefunction for the hydrogen ground state is:

$$\psi_{\text{trial}}(r) = A e^{-\beta r}$$

where $\beta$ is the variational parameter and $A = (\beta^3/\pi)^{1/2}$ ensures normalization. This captures the essential physics: an exponential decay from the nucleus, with a length scale to be optimized.

Computing the energy expectation value:

$$\langle E \rangle = \frac{\hbar^2\beta^2}{2m_e} - \frac{e^2\beta}{4\pi\epsilon_0}$$

The integrals are:

$$\langle T \rangle = -\frac{\hbar^2}{2m_e}\frac{\beta^3}{\pi}\int_0^\infty e^{-\beta r}\nabla^2(e^{-\beta r}) 4\pi r^2\,dr = \frac{\hbar^2\beta^2}{2m_e}$$

$$\langle V \rangle = -\frac{e^2}{4\pi\epsilon_0}\frac{\beta^3}{\pi}\int_0^\infty e^{-2\beta r} \cdot 4\pi r \,dr = -\frac{e^2\beta}{4\pi\epsilon_0}$$

Minimizing with respect to $\beta$:

$$\frac{d\langle E\rangle}{d\beta} = \frac{\hbar^2\beta}{m_e} - \frac{e^2}{4\pi\epsilon_0} = 0$$

$$\beta_{\text{opt}} = \frac{m_e e^2}{4\pi\epsilon_0 \hbar^2} = \frac{1}{a_0}$$

Substituting back:

$$E_{\text{var}} = -\frac{m_e e^4}{2(4\pi\epsilon_0)^2\hbar^2} = -13.6 \text{ eV}$$

The variational method recovers the exact answer! This is because the trial wavefunction $e^{-r/a_0}$ is actually the exact ground-state wavefunction. The variational principle "finds" the right answer because it has enough flexibility to reach it.

💡 Key Insight: The variational principle is not always this clean. For hydrogen, the exponential trial function happens to be the exact answer. For helium (Chapter 19), the best single-exponential gives only about 98% of the exact ground-state energy. The power of the variational method is that it always gives an upper bound, and you can systematically improve it by adding more parameters.

Trial Wavefunction 2: Gaussian

What if we use a Gaussian instead?

$$\psi_{\text{trial}}(r) = A e^{-\alpha r^2}$$

with $A = (2\alpha/\pi)^{3/4}$. This is the wrong functional form (exponential in $r^2$ rather than $r$), so it cannot give the exact answer. But the variational principle guarantees that it gives an upper bound.

Computing:

$$\langle T \rangle = \frac{3\hbar^2 \alpha}{2m_e}$$

$$\langle V \rangle = -\frac{e^2}{4\pi\epsilon_0}\sqrt{\frac{2\alpha}{\pi}} \cdot 2 = -\frac{e^2}{4\pi\epsilon_0}\sqrt{\frac{8\alpha}{\pi}}$$

Minimizing:

$$\alpha_{\text{opt}} = \frac{8m_e^2 e^4}{9\pi(4\pi\epsilon_0)^2\hbar^4}$$

$$E_{\text{var}}^{\text{Gauss}} = -\frac{4}{3\pi}\frac{m_e e^4}{(4\pi\epsilon_0)^2\hbar^2} \cdot 2 = -\frac{8}{3\pi}E_1 = -11.5 \text{ eV}$$

The Gaussian trial gives $E_{\text{var}} = -11.5$ eV compared to the exact $-13.6$ eV. This is a bound that is 15% too high — certainly an upper bound, as guaranteed, but not a very tight one. The Gaussian decays too fast at large $r$ (missing the exponential tail) and too slowly at small $r$ (missing the cusp at the origin).

Trial Wavefunction 3: Two-Parameter Exponential

To demonstrate the power of adding parameters, try:

$$\psi_{\text{trial}}(r) = A(1 + cr)e^{-\beta r}$$

with both $c$ and $\beta$ as variational parameters. This has enough flexibility to capture the radial node structure. For the ground state (which has no nodes), the optimization finds $c \to 0$, recovering the exact answer. But for the $2s$ state (with appropriate orthogonality constraints), this form can capture the single radial node.

Variational Results Summary

The lesson: the exponential $e^{-\beta r}$ is the natural functional form for the Coulomb problem. Any trial function containing it can reach the exact answer. Trial functions that lack the correct functional form (Gaussians, polynomials) give rigorous upper bounds but converge more slowly.

⚠️ Common Misconception: Students sometimes think the variational method is "approximate" and the exact solution is "exact," so the variational method is inferior. This is misleading. For hydrogen, both methods give the same answer. For every other atom and molecule, the exact solution is unknown and the variational method is the primary tool. Density functional theory, Hartree-Fock, configuration interaction, coupled cluster — all modern computational chemistry methods are variational or closely related to variational techniques.

38.5 Numerical Computation

Now let us solve the hydrogen atom numerically, directly on a computer, without using any analytical results. This serves three purposes:

1. Verification — independent check of the analytical and perturbative results
2. Generalization — the same numerical methods work for any central potential, not just Coulomb
3. Preparation — you will need numerical methods for every real-world quantum problem beyond hydrogen

Why Numerical Methods?

You might ask: if we have exact analytical solutions for hydrogen, why bother with numerical methods? Three reasons:

1. Verification. Independent numerical computation provides a cross-check on analytical results. This is not pedantry — sign errors in perturbation theory are common, and numerical verification has caught real mistakes in published calculations.

2. Generalization. The same numerical methods that solve hydrogen also solve helium, lithium, and every other atom — systems where no exact solution exists. Mastering numerical methods on a problem with a known answer prepares you for problems without known answers.

3. Realistic potentials. Real atoms in real environments experience modified potentials — screened Coulomb potentials in plasmas, atoms in external fields, atoms near surfaces, atoms inside quantum dots. Numerical methods handle all of these without modification.

4. Computational thinking. The ability to translate a physical problem into code, solve it numerically, and interpret the results is an essential skill for the modern physicist. Hydrogen is the ideal training ground.

Method 1: Finite Difference (Matrix Diagonalization)

We discretize the reduced radial equation on a grid $r_i = i \Delta r$ for $i = 0, 1, \ldots, N$:

$$-\frac{\hbar^2}{2\mu}\frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta r)^2} + V_{\text{eff}}(r_i) u_i = E u_i$$

where $V_{\text{eff}}(r) = -e^2/(4\pi\epsilon_0 r) + \hbar^2 l(l+1)/(2\mu r^2)$ is the effective potential including the centrifugal barrier.

This converts the differential equation into an $N \times N$ matrix eigenvalue problem:

$$\mathbf{H}\mathbf{u} = E\mathbf{u}$$

where the Hamiltonian matrix has the tridiagonal form:

$$H_{ii} = \frac{\hbar^2}{\mu(\Delta r)^2} + V_{\text{eff}}(r_i), \quad H_{i,i\pm 1} = -\frac{\hbar^2}{2\mu(\Delta r)^2}$$

Diagonalizing this matrix yields the eigenvalues (energies) and eigenvectors (wavefunctions on the grid).

🔴 Warning: The naive uniform grid has a problem: the Coulomb potential diverges as $r \to 0$, and many grid points are wasted in the large-$r$ region where the wavefunction has decayed to zero. A better approach uses a logarithmic grid: $r_i = r_{\min} e^{i\Delta s}$, which concentrates grid points near the nucleus. This is what our code in Section 38.8 implements.

Method 2: Shooting Method

An alternative numerical approach is the shooting method:

1. Pick a trial energy $E$.
2. Integrate the radial equation outward from $r = 0$ (where $u \sim r^{l+1}$) to some matching point $r_m$.
3. Integrate inward from $r = r_{\max}$ (where $u \sim e^{-\kappa r}$ with $\kappa = \sqrt{-2\mu E}/\hbar$) to the same matching point.
4. Adjust $E$ until the outward and inward solutions match smoothly at $r_m$ (continuous value and continuous derivative).

The shooting method is conceptually simple and gives excellent results. It works for any potential and can find excited states by counting the number of nodes.

Node counting rule: The radial wavefunction for quantum numbers $(n, l)$ has exactly $n - l - 1$ nodes (zeros, excluding the boundary points at $r = 0$ and $r = \infty$). This provides a robust way to identify which state you have found.

Method 3: Basis Set Expansion (Spectral Method)

Instead of a grid, expand the wavefunction in a basis of known functions:

$$u(r) = \sum_{k=1}^{N_{\text{basis}}} c_k \phi_k(r)$$

Common basis choices for the Coulomb problem include:

• Slater-type orbitals (STOs): $\phi_k(r) = r^{n_k} e^{-\zeta_k r}$ — the natural choice, matching the true functional form
• Gaussian-type orbitals (GTOs): $\phi_k(r) = r^{n_k} e^{-\alpha_k r^2}$ — wrong functional form but much easier to compute integrals (all integrals are Gaussian and have closed forms)

Substituting into the Schrödinger equation and projecting onto each basis function gives the generalized eigenvalue problem:

$$\mathbf{H}\mathbf{c} = E\mathbf{S}\mathbf{c}$$

where $H_{jk} = \langle \phi_j | \hat{H} | \phi_k \rangle$ and $S_{jk} = \langle \phi_j | \phi_k \rangle$ is the overlap matrix.

📊 By the Numbers: Convergence comparison for the hydrogen ground-state energy with different methods:

Method	$N$ (basis/grid size)	$E_0$ (eV)	Relative Error
Uniform grid	100	$-13.55$	$3.7 \times 10^{-3}$
Uniform grid	1000	$-13.603$	$2.2 \times 10^{-4}$
Log grid	100	$-13.604$	$1.5 \times 10^{-4}$
Log grid	1000	$-13.6057$	$2.2 \times 10^{-6}$
Shooting	(tolerance $10^{-10}$)	$-13.60570$	$< 10^{-6}$
5 STOs	5	$-13.60570$	$< 10^{-6}$
5 GTOs	5	$-13.601$	$3.4 \times 10^{-4}$
20 GTOs	20	$-13.60569$	$< 10^{-5}$

STOs converge fastest because they match the true functional form. GTOs require more basis functions but the integrals are analytically tractable, which is why GTOs dominate computational chemistry.

Basis Set Expansion: From Hydrogen to Chemistry

The basis-set method deserves special attention because it is the foundation of all modern computational quantum chemistry. The idea is to expand the unknown wavefunction in a set of known functions:

$$\psi(r) = \sum_{k=1}^{N} c_k \phi_k(r)$$

and determine the coefficients $c_k$ by solving the generalized eigenvalue problem $\mathbf{H}\mathbf{c} = E\mathbf{S}\mathbf{c}$, where $H_{jk} = \langle \phi_j|\hat{H}|\phi_k\rangle$ and $S_{jk} = \langle \phi_j|\phi_k\rangle$.

For hydrogen, two basis choices are natural:

Slater-type orbitals (STOs): $\phi_k(r) = N_k r^{n_k-1} e^{-\zeta_k r}$, where $\zeta_k$ are optimized exponents. STOs have the correct functional form (exponential decay, cusp at the origin) and converge extremely fast — just 3-5 STOs can reproduce the hydrogen ground-state energy to 6 significant figures. However, the two-electron integrals needed for atoms beyond hydrogen have no closed-form expressions in the STO basis, making them computationally expensive.

Gaussian-type orbitals (GTOs): $\phi_k(r) = N_k r^{n_k-1} e^{-\alpha_k r^2}$. GTOs have the wrong functional form (Gaussian decay instead of exponential, no cusp at the origin) and converge more slowly. But they have a crucial computational advantage: all multi-center integrals can be evaluated analytically using the Gaussian product theorem. This is why GTOs dominate computational chemistry — the modest loss in convergence rate is more than compensated by the dramatic speedup in integral evaluation.

The progression from hydrogen (one electron, exact solution) to helium (two electrons, need basis sets and correlation methods) to lithium and beyond (need Hartree-Fock, configuration interaction, coupled cluster) is a journey from textbook physics to research-grade computational chemistry. Hydrogen is where that journey begins.

🔗 Connection: The basis-set expansion method used here for hydrogen is mathematically identical to the Ritz variational method of Chapter 19. The variational principle guarantees that the lowest eigenvalue of $\mathbf{H}\mathbf{c} = E\mathbf{S}\mathbf{c}$ is an upper bound on the true ground-state energy, and that the bound improves monotonically as the basis is enlarged. This is why basis-set calculations always converge from above.

Numerical Fine-Structure Corrections

We can also compute fine-structure corrections numerically by evaluating the expectation values $\langle 1/r \rangle$, $\langle 1/r^2 \rangle$, $\langle 1/r^3 \rangle$ on the numerical wavefunctions and plugging into the perturbation-theory formulas.

For the $n = 2$ states with a 1000-point logarithmic grid:

The numerical results agree with the analytical values to 4-5 significant figures, providing a strong consistency check.

38.6 Comparison: Analytical vs. Perturbative vs. Variational vs. Numerical

Now we arrive at the heart of this capstone: a systematic comparison of all four methods.

Ground-State Energy Comparison

💡 Key Insight: The discrepancy between the Coulomb-model value ($-13.6057$ eV) and the experimental ionization energy ($-13.5984$ eV) is about 0.05% and is due to: (1) reduced mass correction (already included if using $\mu$ instead of $m_e$), (2) fine structure, (3) Lamb shift, and (4) recoil corrections. The fact that a model as simple as "one electron in a Coulomb potential" gets within 0.05% of experiment is a triumph of quantum mechanics.

Method Comparison: Strengths and Weaknesses

⚖️ Interpretation: No single method is "the best." Exact solutions provide maximum insight but exist only for the simplest systems. Perturbation theory reveals the physics of corrections (which effects are important and why) but may diverge at high orders. Variational methods give rigorous bounds and scale to complex systems, but the quality depends entirely on the choice of trial function. Numerical methods handle any system but can be black-box-like. The practicing physicist uses all four, choosing the tool that best fits the question.

The Fine-Structure Splitting: All Methods Agree

For the $n = 2$ fine-structure splitting ($2P_{3/2}$ vs. $2S_{1/2}/2P_{1/2}$):

The agreement to 7 significant figures between theory (with QED) and experiment is one of the greatest achievements of theoretical physics.

38.7 Spectroscopy Connections

The hydrogen spectrum was the first quantum system to be understood, and it remains the most precisely measured. Let us connect our theoretical results to actual spectroscopic data.

The Electromagnetic Spectrum of Hydrogen

Hydrogen's spectrum is the Rosetta Stone of atomic physics. Every spectral line is a transition between two energy levels, and the collection of all transitions encodes the complete energy-level structure. Let us develop this connection carefully.

The energy of a photon emitted in a transition from level $n$ to level $n'$ ($n > n'$) is:

$$E_{\text{photon}} = E_n - E_{n'} = 13.6\text{ eV}\left(\frac{1}{n'^2} - \frac{1}{n^2}\right)$$

The corresponding photon wavelength is:

$$\lambda = \frac{hc}{E_{\text{photon}}} = \frac{hc}{13.6\text{ eV}}\left(\frac{1}{n'^2} - \frac{1}{n^2}\right)^{-1}$$

This is the Rydberg formula, first discovered empirically by Balmer (1885) for the case $n' = 2$ and generalized by Rydberg. That an empirical formula for spectral lines, discovered decades before quantum mechanics, turns out to be an exact consequence of the Schrödinger equation is one of the most satisfying results in all of physics.

The Spectral Series

Transitions between hydrogen energy levels produce photons with frequencies:

$$\nu = \frac{|E_n - E_{n'}|}{h} = cR_H\left(\frac{1}{n'^2} - \frac{1}{n^2}\right) \quad (n > n')$$

where $R_H = 1.0967758 \times 10^7$ m$^{-1}$ is the Rydberg constant for hydrogen. The spectral series are:

Selection Rules

Not all transitions are equally probable. Electric dipole transitions — the dominant mechanism for photon emission and absorption — obey selection rules derived from the matrix elements $\langle n', l', m' | \hat{\mathbf{r}} | n, l, m \rangle$ (Chapter 21):

$$\Delta l = \pm 1, \quad \Delta m_l = 0, \pm 1, \quad \Delta j = 0, \pm 1 \; (j = 0 \to j' = 0 \text{ forbidden})$$

There is no restriction on $\Delta n$. The selection rule $\Delta l = \pm 1$ is the reason why the $2S_{1/2}$ state is metastable — it cannot decay to $1S_{1/2}$ by single-photon emission because $\Delta l = 0$. It decays instead by two-photon emission, with a lifetime of $\sim 0.14$ seconds — enormous by atomic standards (compare the $2P_{1/2}$ lifetime of $1.6 \times 10^{-9}$ seconds).

The Rydberg Constant: From Hydrogen to Fundamental Physics

The Rydberg constant for hydrogen is:

$$R_H = R_\infty \frac{\mu}{m_e} = R_\infty \left(1 - \frac{m_e}{m_p}\right) = R_\infty \times 0.999456$$

where $R_\infty = m_e e^4/(8\epsilon_0^2 h^3 c) = 10,973,731.568160(21)$ m$^{-1}$ is the Rydberg constant for infinite nuclear mass.

The Rydberg constant is one of the most precisely known fundamental constants. It connects to other constants through:

$$R_\infty = \frac{\alpha^2 m_e c}{2h}$$

This means that a precise measurement of hydrogen spectral lines, combined with the theoretical expression above, determines the fine-structure constant $\alpha$. In practice, the best value of $\alpha$ comes from the electron anomalous magnetic moment, but hydrogen provides an independent check.

The isotope dependence of $R$ is striking: deuterium ($m_D = 2m_p$) has $R_D = R_\infty \times 0.999728$, differing from hydrogen by $0.027\%$. This difference shifts the Balmer-$\alpha$ line by 1.79 Angstroms — easily measurable with a modest spectrograph. Harold Urey discovered deuterium in 1931 by observing this shift, winning the 1934 Nobel Prize in Chemistry.

🔵 Historical Note: The hydrogen spectrum played a decisive role in the development of quantum mechanics at almost every stage: Balmer's formula (1885) demanded explanation, Bohr's model (1913) provided the first quantum derivation, Sommerfeld's extension (1916) introduced the fine structure, and the Lamb shift (1947) launched quantum electrodynamics. No other physical system has so consistently driven theoretical progress.

Transition Rates

The spontaneous emission rate (Einstein $A$ coefficient) for a transition from state $|n, l, j\rangle$ to $|n', l', j'\rangle$ is (Chapter 21):

$$A_{n'l'j' \leftarrow nlj} = \frac{4\alpha\omega^3}{3c^2}\left|\langle n', l', j' \| \hat{\mathbf{r}} \| n, l, j \rangle\right|^2$$

For the Lyman-$\alpha$ transition ($2P \to 1S$):

$$A = 6.27 \times 10^8 \text{ s}^{-1}, \quad \tau = 1.60 \text{ ns}$$

This is one of the strongest atomic transitions in nature. A hydrogen atom in the $2P$ state survives less than 2 nanoseconds before emitting a 121.6-nm ultraviolet photon.

For the Balmer-$\alpha$ transition ($3D_{5/2} \to 2P_{3/2}$, the dominant component):

$$A = 6.47 \times 10^7 \text{ s}^{-1}, \quad \tau = 15.5 \text{ ns}$$

The natural linewidth of a transition is related to the excited-state lifetime by the energy-time uncertainty relation:

$$\Delta\nu = \frac{1}{2\pi\tau}$$

For Lyman-$\alpha$: $\Delta\nu = 99.7$ MHz. For Balmer-$\alpha$: $\Delta\nu = 10.3$ MHz. These linewidths set the fundamental limit on the resolution of hydrogen spectroscopy using single-photon transitions. Sub-Doppler techniques (saturation spectroscopy, two-photon spectroscopy) can eliminate Doppler broadening but not the natural linewidth.

High-Resolution Spectroscopy: What Fine Structure Looks Like

In high-resolution spectroscopy, the Balmer-$\alpha$ line (usually quoted as "656.28 nm") actually consists of seven components, resolved with modern instruments:

1. $3D_{5/2} \to 2P_{3/2}$ — strongest component
2. $3D_{3/2} \to 2P_{3/2}$
3. $3D_{3/2} \to 2P_{1/2}$
4. $3S_{1/2} \to 2P_{3/2}$
5. $3S_{1/2} \to 2P_{1/2}$
6. $3P_{3/2} \to 2S_{1/2}$
7. $3P_{1/2} \to 2S_{1/2}$

The total spread of these components is about 0.36 cm$^{-1}$ — roughly one part in 43,000 of the line's frequency. Resolving these components requires a spectrometer with resolving power $R > 40,000$.

The 1S-2S Transition: The Most Precise Measurement in Physics

The $1S \to 2S$ two-photon transition in hydrogen has been measured with extraordinary precision:

$$\nu_{1S-2S} = 2,466,061,413,187,035(10) \text{ Hz}$$

That is 15 significant figures — one part in $10^{14}$. This measurement, performed by Theodor Hänsch's group in Munich using an optical frequency comb (for which Hänsch shared the 2005 Nobel Prize), is used to:

1. Determine the Rydberg constant to 12 digits
2. Test QED to the highest precision achievable
3. Search for time-variation of fundamental constants
4. Constrain the proton charge radius (leading to the "proton radius puzzle" — see Case Study 2)

📊 By the Numbers: The theoretical prediction for the $1S-2S$ transition frequency, including all known QED corrections up to order $\alpha^7$, agrees with the experimental measurement to 12 significant figures. This is the single most precisely verified prediction in all of science.

38.8 The Complete Hydrogen Simulation

We now describe the capstone code: a comprehensive Python simulation that implements all four methods and compares them against each other and against experimental data. The complete code is in code/example-01-hydrogen-complete.py.

Design Philosophy

The capstone code is designed to be both educational and practical. Every calculation is implemented twice: once using analytical formulas (which are fast and exact) and once using numerical methods (which are general and extensible). This dual implementation serves as a continuous cross-check — if the analytical and numerical results disagree, there is a bug.

The code is organized as a set of independent, composable modules, each corresponding to a section of this chapter. This modularity reflects good software engineering practice, but it also reflects the structure of the physics: the exact solution, perturbation theory, variational methods, and numerical methods are genuinely independent approaches that can be developed and tested separately.

Architecture Overview

The simulation is organized into modules:

class HydrogenAtom:
    """Complete hydrogen atom simulation.

    Methods:
        exact_energy(n)         -> Coulomb energy E_n
        exact_wavefunction(n,l) -> R_nl(r) on grid
        fine_structure(n,l,j)   -> fine-structure correction
        lamb_shift(n,l,j)       -> approximate Lamb shift
        hyperfine(n,l,j,F)      -> hyperfine splitting
        variational(trial_func) -> optimized energy bound
        numerical_solve(l, N)   -> grid-based eigenvalues
        shooting_solve(l, E)    -> shooting method solution
        transition_rate(n1,l1,j1, n2,l2,j2) -> Einstein A
        spectrum(n_max)         -> all transition frequencies
    """

Key Implementation Details

Exact solution module:

The exact wavefunctions are implemented using scipy.special.genlaguerre for the associated Laguerre polynomials and scipy.special.sph_harm for the spherical harmonics. The radial normalization factor is computed using logarithms of factorials to avoid overflow for large $n$.

Perturbation theory module:

Fine-structure corrections are computed from the analytical formula. The Lamb shift is approximated using Bethe's formula with tabulated Bethe logarithms for low $n$. Hyperfine structure uses the exact contact-interaction formula.

Variational module:

The variational solver accepts an arbitrary callable trial_func(r, *params) and uses scipy.optimize.minimize to find the parameters that minimize the energy expectation value. Integrals are computed with scipy.integrate.quad.

Numerical module:

The finite-difference solver constructs the Hamiltonian matrix on a logarithmic grid and diagonalizes it with numpy.linalg.eigh. The shooting solver uses scipy.integrate.solve_ivp with the Dormand-Prince (RK45) method.

Running the Simulation

The complete simulation produces:

1. Energy level diagram — showing gross, fine, Lamb, and hyperfine structure for $n = 1, 2, 3$
2. Wavefunction comparison — analytical vs. numerical vs. variational for the ground state and first few excited states
3. Convergence plot — numerical accuracy vs. grid size/basis size
4. Spectral line positions — predicted transition frequencies compared to NIST data
5. Fine-structure splitting — computed from perturbation theory, Dirac formula, and numerically

Here is a summary of what the code outputs:

=== HYDROGEN ATOM: COMPLETE SIMULATION ===

--- Gross Structure ---
E(1) = -13.6057 eV (exact: -13.6057 eV)
E(2) =  -3.4014 eV (exact:  -3.4014 eV)
E(3) =  -1.5117 eV (exact:  -1.5117 eV)

--- Fine Structure (n=2) ---
2S_1/2: dE = -5.665e-05 eV  (analytical: -5.665e-05 eV)
2P_1/2: dE = -5.665e-05 eV  (analytical: -5.665e-05 eV)
2P_3/2: dE = -1.133e-05 eV  (analytical: -1.133e-05 eV)
Splitting: 0.365 cm^-1 (exp: 0.365 cm^-1)

--- Lamb Shift ---
2S_1/2 - 2P_1/2: 1057.8 MHz (exp: 1057.845 MHz)

--- Hyperfine (1S) ---
F=1 - F=0: 1420.4 MHz (exp: 1420.406 MHz)

--- Variational ---
Trial: exp(-beta*r)  -> E = -13.606 eV (exact)
Trial: exp(-alpha*r^2) -> E = -11.49 eV (bound)

--- Numerical (1000-point log grid) ---
E(1,0) = -13.6054 eV (error: 2.2e-04 rel)
E(2,0) =  -3.4012 eV (error: 5.9e-05 rel)
E(2,1) =  -3.4013 eV (error: 2.9e-05 rel)

--- Spectral Lines ---
Lyman-alpha:  121.567 nm (exp: 121.567 nm)
Balmer-alpha: 656.279 nm (exp: 656.281 nm)
Balmer-beta:  486.133 nm (exp: 486.135 nm)

✅ Checkpoint: Run the code yourself and verify that your output matches the above. If any values disagree by more than 0.01%, check your grid parameters and physical constants. The most common errors are: (1) using $m_e$ instead of the reduced mass $\mu$, (2) getting the Laguerre polynomial convention wrong, (3) insufficient grid resolution near $r = 0$.

Visualization Gallery

The simulation generates several key plots:

Plot 1: Energy Level Diagram. A Grotrian diagram showing the hydrogen energy levels with gross structure, fine structure, and Lamb shift on progressively expanded scales. This is the classic "hydrogen atom poster" that every physics department has on its wall.

Plot 2: Radial Wavefunctions. The first few $R_{nl}(r)$ and $|R_{nl}(r)|^2 r^2$ (radial probability density) plotted on the same axes. The exponential decay, the nodes ($n - l - 1$), and the centrifugal push (the $r^l$ behavior near the origin for $l > 0$) are all visible.

Plot 3: Probability Density Maps. Color-coded cross-sections of $|\psi_{nlm}(r, \theta, \phi)|^2$ in the $xz$-plane for the first 16 states ($n = 1, 2, 3, 4$ and all allowed $l, m$). These are the famous "orbital shapes" that appear in every chemistry textbook.

Plot 4: Method Comparison. A convergence plot showing how the ground-state energy from each method (variational with increasing parameters, numerical with increasing grid size) converges to the exact value. The variational energies approach from above; the numerical values oscillate.

Plot 5: Spectral Lines. A stick spectrum showing the positions and relative intensities of all allowed transitions between $n = 1, 2, 3, 4$, color-coded by series (Lyman = UV/purple, Balmer = visible, Paschen = red/IR).

The Progressive Project: Toolkit Integration

The capstone code integrates modules from throughout the toolkit:

• From Ch 2: Wavefunction class (normalization, probability density)
• From Ch 5: hydrogen_Rnl(), solve_radial(), plot_orbital_3d()
• From Ch 8: Ket, Bra, inner_product()
• From Ch 9: eigensolve()
• From Ch 12: AngularMomentum, wigner_d()
• From Ch 13: SpinState, pauli_matrices()
• From Ch 14: clebsch_gordan()
• From Ch 17: perturbation_1st(), perturbation_2nd()
• From Ch 18: fine_structure(), spin_orbit()
• From Ch 19: variational_solve(), optimize_trial()
• From Ch 21: transition_prob(), dipole_matrix()

The capstone demonstrates that the toolkit you have been building chapter by chapter is now a powerful quantum simulation library.

🔗 Connection: This complete hydrogen simulation will be contrasted with the Bell test simulator (Chapter 39) and the quantum circuit simulator (Chapter 40). Together, the three capstones demonstrate the breadth of quantum mechanics: Chapter 38 is about structure (the static energy levels of an atom), Chapter 39 is about entanglement (the nonlocal correlations that distinguish quantum from classical), and Chapter 40 is about computation (using quantum mechanics to process information).

38.9 The Zeeman Effect: Hydrogen in a Magnetic Field

No treatment of the hydrogen atom is complete without considering its behavior in an external magnetic field — the Zeeman effect. This effect connects atomic structure to one of the most powerful experimental techniques in atomic physics: spectroscopy in controlled magnetic fields.

The Zeeman Hamiltonian

An external magnetic field $\mathbf{B} = B\hat{z}$ couples to both the orbital and spin angular momentum of the electron:

$$\hat{H}_Z = -\boldsymbol{\mu} \cdot \mathbf{B} = \frac{eB}{2m_e}(\hat{L}_z + 2\hat{S}_z) = \frac{\mu_B B}{\hbar}(\hat{L}_z + 2\hat{S}_z)$$

where $\mu_B = e\hbar/(2m_e) = 9.274 \times 10^{-24}$ J/T is the Bohr magneton and the factor of 2 in front of $\hat{S}_z$ is the electron spin g-factor ($g_s \approx 2$).

The behavior depends critically on the relative magnitudes of $\hat{H}_Z$ and $\hat{H}_{\text{fs}}$.

The Weak-Field Regime (Anomalous Zeeman Effect)

When $\mu_B B \ll E_{\text{fs}}$ (which requires $B \ll 0.5$ T for the $n = 2$ states), the Zeeman effect is a perturbation on the fine-structure levels. The good quantum numbers are $n, l, j, m_j$, and the first-order energy correction is:

$$E_Z^{(1)} = g_j m_j \mu_B B$$

where $g_j$ is the Lande g-factor:

$$g_j = 1 + \frac{j(j+1) + 3/4 - l(l+1)}{2j(j+1)}$$

For the ground state ($l = 0, j = 1/2$): $g_j = 2$. For $2P_{1/2}$ ($l = 1, j = 1/2$): $g_j = 2/3$. For $2P_{3/2}$ ($l = 1, j = 3/2$): $g_j = 4/3$.

Each fine-structure level with total angular momentum $j$ splits into $2j + 1$ equally spaced sublevels, but the spacing (determined by $g_j$) differs between levels. This is the "anomalous" Zeeman effect — anomalous because the classical prediction (equal splitting for all levels) fails due to the spin contribution.

The Strong-Field Regime (Paschen-Back Effect)

When $\mu_B B \gg E_{\text{fs}}$ (which requires $B \gg 0.5$ T for $n = 2$), the magnetic field dominates over spin-orbit coupling. The good quantum numbers become $m_l$ and $m_s$ separately, and the energy correction is:

$$E_Z^{(1)} = (m_l + 2m_s)\mu_B B$$

The spin-orbit coupling is then treated as a perturbation on the Zeeman levels, giving:

$$E_{\text{SO}}^{(1)} = \frac{A}{2}m_l m_s$$

where $A$ depends on $n$ and $l$. The pattern simplifies dramatically: levels cluster into groups with the same $m_l + 2m_s$.

The Intermediate Regime

For intermediate fields, neither the weak-field nor the strong-field approximation applies. The full Hamiltonian $\hat{H}_0 + \hat{H}_{\text{fs}} + \hat{H}_Z$ must be diagonalized numerically. The resulting energy levels as a function of $B$ form the Breit-Rabi diagram — a beautiful plot showing the smooth transition from weak-field (anomalous Zeeman) to strong-field (Paschen-Back) behavior. The code in example-01-hydrogen-complete.py generates this diagram for $n = 2$.

🧪 Experiment: The Zeeman effect is not just a textbook exercise — it is one of the most important tools in atomic physics. Zeeman spectroscopy is used to measure magnetic fields in distant stars, to trap and cool atoms with laser light (magneto-optical traps, essential for Bose-Einstein condensation), and to separate nuclear spin states in NMR and MRI imaging. The anomalous Zeeman effect was also one of the key pieces of evidence that led to the discovery of electron spin.

38.10 Pushing Further: Beyond What We Have Computed

Second-Order Perturbation Theory and Higher Corrections

Our fine-structure calculation used only first-order perturbation theory. What happens at higher orders? The second-order correction to the energy involves a sum over all intermediate states:

$$E_n^{(2)} = \sum_{k \neq n} \frac{|\langle k | \hat{H}' | n \rangle|^2}{E_n^{(0)} - E_k^{(0)}}$$

For hydrogen, the second-order fine-structure correction gives terms of order $\alpha^4 E_n$ — the same order as the Lamb shift. At this level, perturbation theory alone is insufficient; one needs the full machinery of quantum electrodynamics.

The lesson is important: perturbation theory is an asymptotic series, not a convergent one. For the fine-structure constant $\alpha \approx 1/137$, the series converges rapidly at low orders, but at very high orders ($\sim 1/\alpha \approx 137$ terms) the terms begin to grow and the series diverges. This divergence is a deep feature of quantum field theory, connected to the renormalization program and the existence of non-perturbative effects (instantons).

For hydrogen, the practical situation is excellent: the first few orders of perturbation theory give extraordinary accuracy, and the higher-order QED corrections have been computed (with heroic effort) to sufficient precision to match experiment.

The hydrogen atom, despite being "solved," continues to surprise. Let us briefly survey the frontiers.

The Proton Radius Puzzle

The Rydberg constant and the proton charge radius are entangled in the theoretical prediction for hydrogen energy levels. In 2010, Randolf Pohl and collaborators measured the Lamb shift in muonic hydrogen (where the electron is replaced by a muon, 207 times heavier), obtaining a proton charge radius 4% smaller than the accepted value from electronic hydrogen spectroscopy and electron scattering. This discrepancy, if real, would imply either new physics or a systematic error in decades of measurements. The puzzle has been partially resolved by newer electronic hydrogen measurements that now agree with the muonic value, but the story is not over. (See Case Study 2.)

Two-Photon Spectroscopy and Optical Frequency Combs

The extraordinary precision of the $1S-2S$ measurement relies on two key innovations: two-photon spectroscopy (which eliminates first-order Doppler broadening) and the optical frequency comb (which directly links optical frequencies to the microwave cesium clock standard). These tools have made hydrogen the testing ground for fundamental physics at the precision frontier.

The Stark Effect: Hydrogen in an Electric Field

Just as the Zeeman effect describes hydrogen in a magnetic field, the Stark effect describes hydrogen in an electric field $\mathbf{E} = \mathcal{E}\hat{z}$:

$$\hat{H}_{\text{Stark}} = e\mathcal{E}z = e\mathcal{E}r\cos\theta$$

The linear Stark effect in hydrogen (first-order energy shift proportional to $\mathcal{E}$) is unique to hydrogen because of the $l$-degeneracy. In atoms without this degeneracy (all atoms other than hydrogen), the Stark effect is quadratic in $\mathcal{E}$ at leading order.

For the $n = 2$ level, the four-fold degenerate states ($2S_{1/2}$ and $2P_{1/2,3/2}$) mix under the electric field. The non-zero matrix element is:

$$\langle 200 | e\mathcal{E}z | 210 \rangle = -3e\mathcal{E}a_0$$

This produces a linear splitting of $\pm 3e\mathcal{E}a_0$ — the linear Stark effect, observable even for modest electric fields.

The Stark effect is important in plasma physics (Stark broadening of hydrogen lines in high-density plasmas is used to measure electron density), in Rydberg atom physics (highly excited hydrogen atoms are extremely sensitive to electric fields), and in astrophysics (interstellar electric fields produce observable Stark shifts in hydrogen masers).

Hydrogen in Extreme Environments

Hydrogen atoms in strong magnetic fields (on neutron star surfaces, $B \sim 10^8$ T) exhibit qualitatively different behavior: the Coulomb and magnetic energies become comparable, the spherical symmetry is broken, and the atom becomes elongated along the field direction. This is the "quadratic Zeeman effect" regime, where perturbation theory fails and numerical methods become essential.

Hydrogen in intense laser fields exhibits above-threshold ionization, high-harmonic generation, and attosecond electron dynamics. These phenomena are at the forefront of ultrafast physics.

Antihydrogen

The ALPHA, ATRAP, and AEgIS experiments at CERN have produced and trapped antihydrogen (a positron bound to an antiproton). Comparing the spectrum of antihydrogen to hydrogen tests CPT symmetry — the most fundamental symmetry in quantum field theory. Any difference would be revolutionary. As of 2024, the $1S-2S$ transition in antihydrogen has been measured to agree with hydrogen to 2 parts in $10^{12}$.

🧪 Experiment: In 2018, the ALPHA collaboration at CERN measured the $1S-2S$ transition in antihydrogen at $2,466,061,103,079.8(5.4)$ kHz, consistent with the hydrogen value to 2 parts in $10^{12}$. This is the most precise test of CPT symmetry in the baryon sector. Improving this measurement is an active area of research.

38.11 Synthesis: What Hydrogen Teaches Us About Quantum Mechanics

Let us step back and consider what the hydrogen atom, taken as a whole, teaches us about quantum mechanics as a theory.

1. Quantum mechanics makes precise, quantitative predictions. The hydrogen spectrum is not vaguely "consistent" with quantum mechanics — it agrees to 12 decimal places. No other scientific theory has been tested to this precision. When people say "quantum mechanics works," this is what they mean.

2. The hierarchy of approximations is the key to understanding. Real physics is never about finding "the" solution. It is about identifying the dominant effect (Coulomb), computing it exactly, and then systematically adding corrections (fine structure, Lamb shift, hyperfine) ordered by their importance. The fine-structure constant $\alpha$ provides the natural organizing parameter.

3. Different methods illuminate different aspects. The exact solution reveals the symmetry structure (the $SO(4)$ symmetry, the degeneracy, the quantum numbers). Perturbation theory reveals the physics of each correction (what causes it, how large it is, what it depends on). Variational methods provide rigorous bounds. Numerical methods handle anything. A complete understanding requires all four.

4. Hydrogen is the foundation of all atomic physics. Every atom beyond hydrogen is, at some level, a perturbation of a hydrogen-like system. The shell structure of the periodic table, the rules governing chemical bonding, the spectra of stars — all trace back to the physics of this single atom.

5. "Solved" does not mean "understood." Despite 100 years of study, hydrogen continues to yield surprises (the proton radius puzzle, tests of CPT symmetry, behavior in extreme fields). A "solved" system is not a closed book — it is a launching pad.

💡 Key Insight: The hydrogen atom is the Drosophila of physics — a model system that is simple enough to study in extraordinary detail, yet rich enough to teach us general principles that apply across all of quantum mechanics. Every new technique in this book was first tested on hydrogen. Every student of quantum mechanics must, at some point, sit down and compute hydrogen from scratch. You have now done so.

Chapter Summary

In this capstone chapter, we have assembled a complete picture of the hydrogen atom by bringing together:

1. The exact Coulomb solution (Chapter 5): energy eigenvalues $E_n = -13.6/n^2$ eV, wavefunctions $\psi_{nlm} = R_{nl}Y_l^m$, quantum numbers $n, l, m$, and the $n^2$-fold degeneracy arising from the hidden $SO(4)$ symmetry.

2. Perturbation corrections (Chapters 17-18): fine structure (relativistic + spin-orbit + Darwin terms, splitting levels by $j$ at order $\alpha^2 E_n$), the Lamb shift (QED corrections splitting $2S_{1/2}$ from $2P_{1/2}$ by 1058 MHz), and hyperfine structure (proton spin interaction splitting the $1S$ ground state into $F = 0, 1$ with the famous 21-cm line at 1420 MHz).

3. Variational bounds (Chapter 19): independent verification of the ground-state energy using trial wavefunctions, demonstrating that the exponential is the natural functional form and that Gaussians converge more slowly.

4. Numerical methods (Chapters 3, 9): finite-difference and shooting-method solutions on a grid, basis-set expansion methods, and convergence analysis — the tools that generalize to any quantum system.

5. Spectroscopy (Chapter 21): transition rates, selection rules, spectral series, and the extraordinary precision of modern hydrogen spectroscopy ($1S-2S$ measured to 15 significant figures).

The hydrogen atom is not just a textbook problem. It is the testing ground for quantum mechanics, quantum electrodynamics, and fundamental physics. The agreement between theory and experiment — across 12 orders of magnitude in precision — is the strongest evidence we have that quantum mechanics correctly describes the physical world.

You have now seen the hydrogen atom from every angle available in this book. The progressive project toolkit you have built contains all the pieces needed for this simulation. Run the code, make the plots, and verify the numbers for yourself. That is the final exam.

"We have always found that each group of phenomena that has been explained points to new phenomena not yet explained." — Richard Feynman, QED: The Strange Theory of Light and Matter (1985)