Dirac equation in a Coulomb field

Intuition Beginner

The non-relativistic hydrogen atom has a beautifully simple spectrum: energy levels labelled by a single integer $n$, with $E_n=-13.6/n^2$ eV. Every state with the same $n$ has the same energy. The $2s$ and $2p$ orbitals sit at exactly the same level. The electron's spin and its orbital motion are independent bookkeeping.

A spectrometer with enough resolution tells a different story. Each "single" line of the Balmer series splits into several closely spaced lines — the fine structure of hydrogen, first observed by Michelson and Morley in 1887. The splittings are small: of order $13.6\times {\alpha }^2$ eV, where $\alpha \approx 1/137$ is the fine-structure constant. That puts the gap between $2p_{3/2}$ and $2p_{1/2}$ at about $10^{-4}$ eV — small, but unmistakably there.

Sommerfeld guessed the splittings in 1916 from a relativistic correction to Bohr's quantised circular orbits, getting a formula that matched experiment to surprising precision. The old quantum theory was wrong about almost everything, but Sommerfeld's formula stuck. When Dirac wrote down his relativistic equation for the electron in 1928, the very first application — done independently by Darwin and Gordon within weeks — was to solve it in the Coulomb field of a proton. Out came the Sommerfeld formula, no fitting required, with the spin and the relativistic kinematics fused into one calculation.

The picture the Dirac calculation paints is this. The electron now has four components — two for spin, two for the particle-antiparticle structure imposed by special relativity. The Coulomb attraction $-Z{e}^2/r$ couples to all four. Total angular momentum $\mathbf{J}=\mathbf{L}+\mathbf{S}$ is conserved, but orbital $\mathbf{L}$ and spin $\mathbf{S}$ separately are not. The energy depends on the total angular momentum quantum number $j$, not on $\ell$ alone — so the $2s_{1/2}$ (which has $\ell =0$, $j=1/2$) and the $2p_{1/2}$ (which has $\ell =1$, $j=1/2$) come out degenerate in the Dirac equation, with $2p_{3/2}$ a little higher.

This $j$-degeneracy is exactly what experiment shows — to a part in $10^4$. The remaining splitting between $2s_{1/2}$ and $2p_{1/2}$, the famous Lamb shift of about 1058 MHz measured by Lamb and Retherford in 1947, comes from quantum-electrodynamic effects beyond the single-particle Dirac equation. The Dirac-Coulomb calculation is the foundation; quantum electrodynamics is the next storey.

Visual Beginner

The non-relativistic spectrum has the $n=2$ shell as one flat line: $2s$ and $2p$ share an energy. The Dirac equation breaks this into two: $2p_{3/2}$ at the top, $2s_{1/2}$ and $2p_{1/2}$ jointly at the bottom. The gap between them is the fine-structure splitting, of order ${\alpha }^{4}m{c}^2\approx 4.5\times 10^{-5}$ eV in hydrogen. The remaining $2s_{1/2}$–$2p_{1/2}$ degeneracy is a prediction of the Dirac equation that experiment violates by about 1 part in $10^7$ — a tiny but real discrepancy that quantum electrodynamics accounts for via the Lamb shift.

Worked example Beginner

Compute the leading relativistic correction to the hydrogen ground-state energy.

The Sommerfeld formula, expanded to leading order in $Z\alpha$, reads

$$E_{n,j}=m{c}^2\left[1-\frac{(Z\alpha {)}^2}{2n^2}-\frac{(Z\alpha {)}^4}{2n^3}\left(\frac{1}{j+1/2}-\frac{3}{4n}\right)+\cdots \right].$$

The first term is the rest energy $m{c}^2$. The second is the familiar Bohr energy $E_n^{\text{Bohr}}=-m{c}^2(Z\alpha {)}^2/(2n^2)$. The third is the fine-structure correction.

For hydrogen ($Z=1$) in the ground state ($n=1$, $j=1/2$): the Bohr energy is $-m{c}^2{\alpha }^2/2=-13.606$ eV (using $m{c}^2=511000$ eV and $\alpha =1/137.036$). The fine-structure correction is

$$\Delta E_{1,1/2}^{\text{fs}}=-\frac{m{c}^2{\alpha }^4}{2}\left(\frac{1}{1}-\frac{3}{4}\right)=-\frac{m{c}^2{\alpha }^4}{8}\approx -1.81\times 10^{-4}\text{eV}.$$

So the relativistic ground state sits about 0.18 meV below the Bohr ground state. The corresponding wavelength shift is small but observable. The general pattern: the fine-structure correction lowers all levels (the bracket is positive for the cases that matter), and the magnitude scales as ${\alpha }^2$ times the non-relativistic energy, with extra $1/n$ and $j$-dependence inside.

What this tells us: Dirac's equation does not just add a tiny correction. It changes the labelling. The state that used to be called "$2s$" and the state that used to be called "$2p_{1/2}$" are now degenerate because both have $j=1/2$. The "$2p_{3/2}$" — same $\ell =1$ as $2p_{1/2}$ but with $j=3/2$ — sits above them. Angular momentum quantum numbers reorganise themselves into the relativistically correct combination.

Check your understanding Beginner

Formal definition Intermediate+

Work in natural units $\hslash =c=1$ with the mostly-minus metric ${\eta }^{\mu \nu }=\mathrm{diag}(+1,-1,-1,-1)$. The Dirac equation for an electron of mass $m$ and charge $-e$ ($e>0$) in an external electromagnetic four-potential $A_{\mu }$ is obtained by minimal coupling 12.11.01:

$$\boxed{\left[{\gamma }^{\mu }(i{\partial }_{\mu }+e{A}_{\mu })-m\right]\psi (x)=0.}$$

For a static central Coulomb potential created by a point nucleus of charge $+Ze$ at the origin,

$$A_0=-\frac{Ze}{r},\mathbf{A}=0,e{A}_0=-\frac{Z\alpha }{r},$$

where $\alpha =e^2/(4\pi )$ in natural Heaviside-Lorentz units (numerically $\alpha \approx 1/137.036$). For a stationary state $\psi (x)=e^{-iEt}\psi (\mathbf{r})$ with energy $E$, the equation reduces to

$$\left[\mathbf{\alpha }\cdot \mathbf{p}+\beta m-\frac{Z\alpha }{r}\right]\psi (\mathbf{r})=E\psi (\mathbf{r}),$$

where $\mathbf{\alpha }={\gamma }^0\mathbf{\gamma }$ and $\beta ={\gamma }^0$ are the Dirac matrices in standard ("Dirac") representation. This is a stationary first-order matrix differential equation on ${\mathbb{R}}^3$ for the four-component spinor $\psi (\mathbf{r})\in {\mathbb{C}}^4$.

Conserved operators and the κ-operator

The Hamiltonian $H=\mathbf{\alpha }\cdot \mathbf{p}+\beta m-Z\alpha /r$ commutes with neither orbital angular momentum $\mathbf{L}=\mathbf{r}\times \mathbf{p}$ nor with the $4\times 4$ spin $\mathbf{S}=\mathbf{\Sigma }/2$ where $\mathbf{\Sigma }=\mathrm{diag}(\mathbf{\sigma },\mathbf{\sigma })$ in the standard basis. The conserved angular momentum is total $\mathbf{J}=\mathbf{L}+\mathbf{S}$, with $[H,\mathbf{J}]=0$, $[H,J^2]=0$, and $[H,J_z]=0$. Parity $P\psi (\mathbf{r}):=\beta \psi (-\mathbf{r})$ also commutes with $H$.

A fifth conserved operator, due to Dirac, is the κ-operator

$$K:=\beta (\mathbf{\Sigma }\cdot \mathbf{L}+1).$$

Direct computation using $\{\mathbf{\alpha },\beta \}=0$ and $[\mathbf{\alpha }\cdot \mathbf{p},\mathbf{\Sigma }\cdot \mathbf{L}]=i\mathbf{\alpha }\cdot (\mathbf{L}\times \mathbf{p}-\mathbf{p}\times \mathbf{L})$ verifies $[H,K]=0$. The eigenvalues of $K$ on a state of orbital angular momentum $\ell$ and total angular momentum $j=\ell \pm 1/2$ are

$$\kappa =\mp (j+1/2)=\{\begin{array}{ll}-(\ell +1)& \text{for }j=\ell +1/2,\\ \ell & \text{for }j=\ell -1/2.\end{array}$$

The integer $\kappa$ is the Dirac quantum number. It packages the pair $(\ell ,j)$ into a single signed integer and turns out to be the natural quantum number for the radial problem.

Spinor spherical harmonics

The simultaneous eigenfunctions of $J^2$, $J_z$, $L^2$, and parity on two-component spinors are the spinor spherical harmonics ${\Omega }_{j\ell m}(\theta ,\varphi )$, defined by Clebsch-Gordan combination of orbital spherical harmonics $Y_{\ell ,m_{\ell }}$ with the Pauli spin eigenstates ${\chi }_{\pm }$:

$${\Omega }_{j\ell m}(\theta ,\varphi )=\sum _{m_s=\pm 1/2}⟨\ell ,m-m_s;\frac{1}{2},m_s|j,m⟩Y_{\ell ,m-m_s}(\theta ,\varphi ){\chi }_{m_s}.$$

For $j=\ell +1/2$,

$${\Omega }_{j,\ell ,m}=\left(\begin{array}{c}\sqrt{\frac{\ell +m+1/2}{2\ell +1}}{Y}_{\ell ,m-1/2}\\ \sqrt{\frac{\ell -m+1/2}{2\ell +1}}{Y}_{\ell ,m+1/2}\end{array}\right);$$

for $j=\ell -1/2$,

$${\Omega }_{j,\ell ,m}=\left(\begin{array}{c}-\sqrt{\frac{\ell -m+1/2}{2\ell +1}}{Y}_{\ell ,m-1/2}\\ \sqrt{\frac{\ell +m+1/2}{2\ell +1}}{Y}_{\ell ,m+1/2}\end{array}\right).$$

The key identity that drives the radial reduction is $(\mathbf{\sigma }\cdot \hat{\mathbf{r}}){\Omega }_{j\ell m}=-{\Omega }_{j{\ell }^{\prime }m}$ where ${\ell }^{\prime }=2j-\ell$ is the other orbital value compatible with the same $j$: if $\ell =j-1/2$ then ${\ell }^{\prime }=j+1/2$, and vice versa. The operator $\mathbf{\sigma }\cdot \hat{\mathbf{r}}$ is a parity-odd, $J$-conserving operator that flips $\ell \leftrightarrow {\ell }^{\prime }$ within a fixed $(j,m)$ multiplet.

Angular separation ansatz

Write the four-component spinor in block form,

$$\psi (\mathbf{r})=\left(\begin{array}{c}\varphi (\mathbf{r})\\ \chi (\mathbf{r})\end{array}\right)=\frac{1}{r}\left(\begin{array}{c}g(r){\Omega }_{j\ell m}(\theta ,\varphi )\\ if(r){\Omega }_{j{\ell }^{\prime }m}(\theta ,\varphi )\end{array}\right),$$

where $\varphi$ is the upper ("large") two-component spinor, $\chi$ is the lower ("small") two-component spinor, and ${\ell }^{\prime }=2j-\ell$. The factor $1/r$ in front converts the three-dimensional radial measure $r^{2}dr$ into the one-dimensional measure $dr$, so the radial functions $g(r),f(r)$ are normalised by ${\int }_0^{\infty }(|g{|}^2+|f{|}^2)dr=1$. The factor $i$ in the lower component is a phase convention that makes the radial equations real.

The opposite orbital angular momenta $\ell$ and ${\ell }^{\prime }=2j-\ell$ in the two components are forced by parity: the standard $\beta =\mathrm{diag}(1,-1)$ assigns opposite intrinsic parities to upper and lower spinors, so for $\psi$ to be a parity eigenstate the spatial parts must carry orbital parities that differ by one, i.e., $(-1{)}^{\ell }=-(-1{)}^{{\ell }^{\prime }}$, which gives ${\ell }^{\prime }=\ell \pm 1$. The choice ${\ell }^{\prime }=2j-\ell$ selects the orbital partner with the same total angular momentum $j$.

Radial Dirac system

Inserting the ansatz into the stationary Dirac equation and using the identity $(\mathbf{\sigma }\cdot \mathbf{p})[f(r){\Omega }_{j\ell m}/r]=i[f^{\prime }+(1+\kappa )f/r]{\Omega }_{j{\ell }^{\prime }m}/r$ (proved via $\mathbf{\sigma }\cdot \mathbf{p}=-i(\mathbf{\sigma }\cdot \hat{\mathbf{r}})({\partial }_r+(1-\kappa )/r\cdot \hat{r}\text{-acting})$ on $\Omega$) collapses the angular dependence and yields the coupled radial Dirac system for $g(r)$ and $f(r)$:

$$\boxed{\frac{dg}{dr}+\frac{\kappa }{r}g-\left(m+E+\frac{Z\alpha }{r}\right)f=0,}$$ $$\boxed{\frac{df}{dr}-\frac{\kappa }{r}f+\left(m-E-\frac{Z\alpha }{r}\right)g=0.}$$

These two first-order ODEs determine the spectrum and the radial wavefunctions completely once boundary conditions are imposed: regularity at the origin (the integrals $\int |g{|}^2+|f{|}^2$ must converge near $r=0$) and exponential decay at infinity (for bound states, $E<m$). The Dirac quantum number $\kappa$ enters as the only angular input; the rest of the calculation is a one-dimensional eigenvalue problem.

Key derivation Intermediate+

Theorem (Sommerfeld fine-structure formula). The discrete spectrum of the radial Dirac-Coulomb system for $0<Z\alpha <1$ is

$$\boxed{E_{n,j}=m{\left[1+\frac{(Z\alpha {)}^2}{(n-j-\frac{1}{2}+\sqrt{(j+\frac{1}{2}{)}^2-(Z\alpha {)}^2}{)}^2}\right]}^{-1/2},}$$

where $n=1,2,3,\dots$ is the principal quantum number and $j=1/2,3/2,\dots ,n-1/2$ runs over the allowed total angular momenta. The eigenvalue depends on $n$ and $j$ alone, not on $\ell$, so for each $j\in \{1/2,\dots ,n-1/2\}$ the two states with $\ell =j-1/2$ and $\ell =j+1/2$ (when both exist) are degenerate.

Proof. Eliminate one of the two radial functions to reduce the coupled first-order system to a single second-order ODE, then solve via the standard confluent-hypergeometric ansatz.

Introduce the scaled radial coordinate $\rho :=\lambda r$ where $\lambda :=\sqrt{m^2-E^2}$ (real and positive for bound states $|E|<m$), and the dimensionless parameter $\nu :=Z\alpha E/\lambda$. The radial system becomes

$$\rho g^{\prime }+\kappa g=\left[\frac{m+E}{\lambda }\rho +Z\alpha \right]f,\rho f^{\prime }-\kappa f=-\left[\frac{m-E}{\lambda }\rho +Z\alpha \right]g.$$

(Here primes denote $d/d\rho$.) Multiply the first equation by $\sqrt{(m-E)/(m+E)}$ and the second by $\sqrt{(m+E)/(m-E)}$ and add/subtract to introduce the linear combinations

$${\Phi }_{\pm }(\rho ):=\sqrt{m\pm E}[g(\rho )\pm f(\rho )].$$

These combinations diagonalise the mass term and put the system into the form

$$\rho {\Phi }_{\pm }^{\prime }+(\kappa \mp \nu ){\Phi }_{\pm }=\mp (\rho \mp Z\alpha m/\lambda ){\Phi }_{\mp }.$$

Look for asymptotic behaviour. At large $\rho$, the system reduces to $\rho {\Phi }_{\pm }^{\prime }\sim \mp \rho {\Phi }_{\mp }$, with solutions ${\Phi }_{\pm }\sim e^{\mp \rho }$. The physically acceptable boundary condition at infinity selects ${\Phi }_{\pm }\sim e^{-\rho }$ for bound states. Near $\rho =0$, the system is $\rho {\Phi }_{\pm }^{\prime }+(\kappa \mp \nu ){\Phi }_{\pm }\sim \pm Z\alpha m{\Phi }_{\mp }/\lambda$, giving power-law behaviour ${\Phi }_{\pm }\sim {\rho }^{\gamma }$ with index $\gamma$ determined by the secular condition

$$\det \left(\begin{array}{cc}\gamma +\kappa & -Z\alpha m/\lambda \\ -Z\alpha m/\lambda & \gamma -\kappa \end{array}\right)\stackrel{!}{=}0⟹{\gamma }^2={\kappa }^2-(Z\alpha {)}^2(m/\lambda {)}^2\cdot 0\dots$$

A more careful matching using the full leading-order system around $r=0$ (where the Coulomb $1/r$ term dominates) gives the index

$$\gamma =\sqrt{{\kappa }^2-(Z\alpha {)}^2}.$$

The two roots $\pm \gamma$ correspond to regular and singular behaviour; only the regular root $+\gamma$ keeps ${\int }_0^{ϵ}|g{|}^{2}dr$ finite for $Z\alpha <|\kappa |\le \sqrt{3}/2\cdot 2=\sqrt{3}$ (and more generally for $Z\alpha <1$ with the distinguished self-adjoint extension).

Write ${\Phi }_{\pm }(\rho )={\rho }^{\gamma }{e}^{-\rho }{F}_{\pm }(\rho )$ and substitute. The functions $F_{\pm }$ satisfy a hypergeometric-type system whose solutions are confluent hypergeometric functions ${}_1{F}_1$. The bound-state quantisation condition arises from the requirement that the power series for $F_{\pm }$ terminate (otherwise ${\Phi }_{\pm }\sim e^{+\rho }$ at infinity, violating normalisability). Termination requires

$$\gamma -\nu =-n_r,n_r=0,1,2,\dots ,$$

where $n_r$ is the radial quantum number (the highest power of $\rho$ in $F_{+}$). Substituting $\nu =Z\alpha E/\lambda$ and $\gamma =\sqrt{{\kappa }^2-(Z\alpha {)}^2}$:

$$Z\alpha \frac{E}{\lambda }=n_r+\sqrt{{\kappa }^2-(Z\alpha {)}^2}.$$

Solve for $E$. Using ${\lambda }^2=m^2-E^2$,

$$(Z\alpha {)}^2{E}^2=(m^2-E^2)[n_r+\sqrt{{\kappa }^2-(Z\alpha {)}^2}{]}^2,$$ $$E^2[(Z\alpha {)}^2+(n_r+\sqrt{{\kappa }^2-(Z\alpha {)}^2}{)}^2]=m^2(n_r+\sqrt{{\kappa }^2-(Z\alpha {)}^2}{)}^2,$$ $$E=m{\left[1+\frac{(Z\alpha {)}^2}{(n_r+\sqrt{{\kappa }^2-(Z\alpha {)}^2}{)}^2}\right]}^{-1/2}.$$

Finally, identify the principal quantum number $n:=n_r+|\kappa |=n_r+j+1/2$. Substituting $|\kappa |=j+1/2$,

$$n_r+\sqrt{{\kappa }^2-(Z\alpha {)}^2}=n-(j+1/2)+\sqrt{(j+1/2{)}^2-(Z\alpha {)}^2},$$

which reproduces the boxed Sommerfeld formula. $\square$

Corollary 1 (Bohr limit). Expanding to leading order in $Z\alpha$, $E_{n,j}\approx m-m(Z\alpha {)}^2/(2n^2)$, the non-relativistic Schrödinger-Coulomb spectrum 12.06.01 depending on $n$ alone.

Corollary 2 (Fine-structure expansion). To next order, $$ E_{n,j} = m\left[1 - \frac{(Z\alpha)^2}{2n^2} - \frac{(Z\alpha)^4}{2n^3}\left(\frac{1}{j + 1/2} - \frac{3}{4n}\right) + O((Z\alpha)^6)\right]. $$

The fine-structure correction depends on $j$ but not on $\ell$ separately. The four contributions encoded in the $(Z\alpha {)}^4$ term — relativistic kinetic correction, spin-orbit coupling, Darwin term, and the residual rearrangement — are exactly what the Foldy-Wouthuysen non-relativistic expansion of the Dirac equation produces 12.11.01, and exactly what the Sommerfeld 1916 quantised-orbit calculation guessed empirically.

Corollary 3 (Degeneracy). For each $j\in \{1/2,3/2,\dots ,n-1/2\}$ and each $n\ge 2$, the two states with $\ell =j-1/2$ and $\ell =j+1/2$ (the latter requires $\ell +1/2\le n-1/2$, i.e., $\ell \le n-1$) share the same Sommerfeld energy. The $2s_{1/2}$ ($\ell =0$, $j=1/2$) and $2p_{1/2}$ ($\ell =1$, $j=1/2$) are an example.

Worked example: hydrogen $n=2$ levels

For hydrogen ($Z=1$) at $n=2$, the allowed $(\ell ,j)$ pairs are $(0,1/2)\to 2s_{1/2}$, $(1,1/2)\to 2p_{1/2}$, and $(1,3/2)\to 2p_{3/2}$. Apply the Sommerfeld formula with $\alpha \approx 1/137.036$:

$$E_{2,1/2}=m{\left[1+\frac{{\alpha }^2}{(2-1+\sqrt{1-{\alpha }^2}{)}^2}\right]}^{-1/2}=m{\left[1+\frac{{\alpha }^2}{(1+\sqrt{1-{\alpha }^2}{)}^2}\right]}^{-1/2};$$ $$E_{2,3/2}=m{\left[1+\frac{{\alpha }^2}{(2-2+\sqrt{4-{\alpha }^2}{)}^2}\right]}^{-1/2}=m{\left[1+\frac{{\alpha }^2}{4-{\alpha }^2}\right]}^{-1/2}.$$

Expanding to fourth order in $\alpha$:

$$E_{2,1/2}-m\approx -\frac{m{\alpha }^2}{8}-\frac{5m{\alpha }^4}{128}\approx -3.4014\text{eV}-4.53\times 10^{-5}\text{eV},$$ $$E_{2,3/2}-m\approx -\frac{m{\alpha }^2}{8}-\frac{m{\alpha }^4}{128}\approx -3.4014\text{eV}-0.91\times 10^{-5}\text{eV}.$$

The fine-structure splitting is $E_{2,3/2}-E_{2,1/2}\approx 3.62\times 10^{-5}\text{eV}\approx 10.95\text{GHz}\times h$, in excellent agreement with the spectroscopic value $\sim 10.969$ GHz. The Lamb shift adds an additional $\sim 1058$ MHz splitting between $2s_{1/2}$ and $2p_{1/2}$ that the Sommerfeld formula does not capture.

Exercises Intermediate+

Lean formalization Intermediate+

Mathlib's relevant infrastructure covers parts of the foundation but does not assemble the Dirac-Coulomb result. The components are:

• Mathlib.LinearAlgebra.CliffordAlgebra.Basic provides the abstract Clifford algebra over a commutative ring with quadratic form, the universal property, and basic functoriality. The Cl(1,3) specialisation needed for the gamma matrices is a direct instantiation.
• Mathlib.Analysis.SpecialFunctions.Sphere.SphericalHarmonic (planned; current coverage via Laplacian eigenfunctions) gives orbital spherical harmonics $Y_{\ell m}$, but the spinor spherical harmonics ${\Omega }_{j\ell m}$ require Clebsch-Gordan combination with Pauli spinors, an unimplemented physics-layer construction.
• The radial Dirac system's confluent-hypergeometric solution would build on Mathlib.Analysis.SpecialFunctions.Gamma and a yet-unwritten KummerFunction library.
• The self-adjoint extension theory of the Dirac-Coulomb operator on $C_0^{\infty }({\mathbb{R}}^3\setminus \{0\};{\mathbb{C}}^4)$ is missing entirely. The Kato-Rellich and von Neumann deficiency-index machinery exists in Mathlib.Analysis.NormedSpace.Spectrum but has not been specialised to matrix-valued first-order differential operators with $1/r$ singularities.

lean_status: none. The unit ships without a Lean module path; the Mathlib gap analysis frontmatter field records the specific gap. Aggregated with the other relativistic-QM none units, this feeds the upstream Mathlib physics-formalisation roadmap.

Advanced results Master

The κ-operator and the relativistic Johnson-Lippmann symmetry

The conservation of $K=\beta (\mathbf{\Sigma }\cdot \mathbf{L}+1)$ is straightforward once the algebra is unwound: $[H_0,K]=0$ for $H_0=\mathbf{\alpha }\cdot \mathbf{p}+\beta m$ follows from $\{\mathbf{\alpha },\beta \}=0$ and the angular-momentum commutator $[\mathbf{\alpha }\cdot \mathbf{p},\mathbf{\Sigma }\cdot \mathbf{L}]=i\mathbf{\alpha }\cdot (\mathbf{L}\times \mathbf{p}+\mathbf{p}\times \mathbf{L})=i\mathbf{\alpha }\cdot (2\mathbf{p}\times \mathbf{p}-i\mathbf{p})=-\mathbf{\alpha }\cdot \mathbf{p}$ (using $\mathbf{L}\times \mathbf{p}+\mathbf{p}\times \mathbf{L}=i\mathbf{p}$ from the canonical commutators). Adding the $\beta$ factor and the $+1$ flips the sign appropriately to cancel against $[\mathbf{\alpha }\cdot \mathbf{p},\beta ]=-2\mathbf{\alpha }\cdot \mathbf{p}\beta /\beta \beta =\dots$ The Coulomb potential $-Z\alpha /r$ depends only on $r$, so it commutes term by term with $\mathbf{\Sigma }\cdot \mathbf{L}$ (which acts on angular variables and on the spinor index, neither of which the scalar potential touches), and $[H,K]=0$ extends from the free to the Coulomb case.

The κ-operator is essentially $K=\beta (J^2-L^2-S^2+1)=\beta (J^2-L^2+1/4)$. Its eigenvalues $\kappa =\mp (j+1/2)$ for $j=\ell \pm 1/2$ package $(\ell ,j)$ into a single signed integer that takes the role of "angular momentum" in the radial reduction.

Beyond $K$, the Dirac-Coulomb problem admits a relativistic analogue of the Runge-Lenz vector, the Johnson-Lippmann operator. For the non-relativistic Schrödinger-Coulomb problem, the Runge-Lenz vector $\mathbf{A}=(1/2m)(\mathbf{p}\times \mathbf{L}-\mathbf{L}\times \mathbf{p})-Z\alpha \hat{\mathbf{r}}$ commutes with the Hamiltonian; together with $\mathbf{L}$ it closes into the Lie algebra of $\mathrm{SO}(4)$ and explains the $n^2$-fold accidental degeneracy. The relativistic generalisation, found by Johnson and Lippmann in 1950, is

$${\hat{A}}_{\text{JL}}=\mathbf{\Sigma }\cdot \hat{\mathbf{r}}K-Z\alpha \beta \mathbf{\Sigma }\cdot \hat{\mathbf{r}}.$$

This operator commutes with the Dirac-Coulomb Hamiltonian and connects the $\ell =j-1/2$ and $\ell =j+1/2$ states at fixed $(n,j)$, providing the algebraic origin of the residual $\ell$-degeneracy in the Sommerfeld spectrum. The associated symmetry group is not $\mathrm{SO}(4)$ but a relativistic deformation whose representation theory recovers the Sommerfeld energies as the eigenvalues of a Casimir-like operator. Vacuum polarisation and QED self-energy break this symmetry beyond tree level, which is why the Lamb shift lifts the residual degeneracy.

Radial wavefunctions

Once the Sommerfeld energy is fixed, the radial functions $g(r),f(r)$ are determined up to normalisation. They take the form

$$g(r)=\mathcal{N}(2\lambda r{)}^{\gamma }{e}^{-\lambda r}[(n_r+\gamma ){}_1{F}_1(-n_r,2\gamma +1;2\lambda r)-(\kappa -Z\alpha m/\lambda ){}_1{F}_1(1-n_r,2\gamma +1;2\lambda r)],$$

with a parallel expression for $f(r)$ differing in coefficients and signs (Bethe-Salpeter §14 eq. (14.22) lists the explicit form). The normalisation constant $\mathcal{N}$ involves the Gamma function:

$$\mathcal{N}=\frac{\sqrt{2{\lambda }^{2\gamma +1}\Gamma (2\gamma +n_r+1)/n_r!}}{\Gamma (2\gamma +1)}\sqrt{\frac{m\pm E}{4m(n_r+\gamma )(n_r+\gamma -Z\alpha m/\lambda )}}.$$

For $1s_{1/2}$ ($n=1$, $\kappa =-1$, $n_r=0$, $\gamma =\sqrt{1-{\alpha }^2}$, $\lambda =m\alpha \gamma /(n_r+\gamma )=m\alpha$ exactly), the radial functions simplify dramatically:

$$g_{1s}(r)=N_g(2m\alpha r{)}^{\gamma }{e}^{-m\alpha r},f_{1s}(r)=-\frac{Z\alpha m+\lambda \kappa }{Z\alpha m-\lambda \kappa }{g}_{1s}(r)\cdot (\text{sign factor}).$$

The upper component $g_{1s}/r\cdot {\Omega }_{1/2,0,m}$ resembles the Schrödinger $1s$ wavefunction $(m\alpha {)}^{3/2}{e}^{-m\alpha r}/\sqrt{\pi }$ but with the additional $r^{\gamma -1}$ factor, where $\gamma -1\approx -{\alpha }^2/2\approx -2.6\times 10^{-5}$. The mild singularity at the origin is the relativistic deviation from Schrödinger; the lower-component-to-upper-component ratio $|f/g{|}_{r\sim a_0}\sim Z\alpha$ controls how "small" the small component is in the non-relativistic limit.

For heavy hydrogenic atoms ($Z\gg 1$), the small-component admixture becomes substantial. At $Z=80$ (mercurial $K$-shell), $Z\alpha \approx 0.58$, and $|f/g{|}^2$ near the nucleus reaches $\sim 0.3$ — the electron is genuinely relativistic, and its non-relativistic wavefunction is a poor approximation. This is the structural origin of the chemistry of heavy elements: gold's colour, mercury's liquidity at room temperature, and the inert-pair effect in lead all trace to the contraction of the $6s$ orbital under relativistic kinematics.

Expansion of the Sommerfeld formula

The systematic $Z\alpha$ expansion of the Sommerfeld formula reads

$$E_{n,j}=m[1-\frac{(Z\alpha {)}^2}{2n^2}-\frac{(Z\alpha {)}^4}{2n^3}\left(\frac{1}{j+1/2}-\frac{3}{4n}\right)-\frac{(Z\alpha {)}^6}{8n^4}\left(\frac{1}{(j+1/2{)}^3}-\frac{6}{n(j+1/2{)}^2}\cdot \dots \right)+\cdots ].$$

The terms have the following physical interpretation, identified by the Foldy-Wouthuysen reduction of the Dirac equation 12.11.01:

• The $(Z\alpha {)}^2$ term is the non-relativistic Bohr energy $-m(Z\alpha {)}^2/(2n^2)$.
• The $(Z\alpha {)}^4$ term groups together (i) the relativistic kinetic-energy correction $-⟨p^4⟩/(8m^3)$, (ii) the spin-orbit coupling $⟨(1/2m^2)(1/r)(dV/dr)\mathbf{L}\cdot \mathbf{S}⟩$ (present only for $\ell >0$), and (iii) the Darwin term $⟨(1/8m^2){\nabla }^{2}V⟩=-(Z\alpha )\pi {\delta }^3(\mathbf{r})/(2m^2)$ (present only for $\ell =0$). The remarkable cancellation: the kinetic, spin-orbit, and Darwin contributions sum to a $j$-only result, with the $\ell$-dependence cancelling between spin-orbit (active for $p$ states) and Darwin (active for $s$ states). This is the Dirac-equation prediction of $2s_{1/2}$-$2p_{1/2}$ degeneracy.
• The $(Z\alpha {)}^6$ term is the next correction; at hydrogen scales it is of order $10^{-9}$ eV and well below the experimental precision needed to detect it directly.

The Pomeranchuk-Smorodinsky catastrophe

At $Z\alpha =|\kappa |$, the power-law exponent $\gamma =\sqrt{{\kappa }^2-(Z\alpha {)}^2}$ vanishes; at $Z\alpha >|\kappa |$, $\gamma$ becomes imaginary, $\gamma =i|\gamma |$, and the wavefunction near the origin oscillates as $(2\lambda r{)}^{i|\gamma |}=\cos (|\gamma |\ln (2\lambda r))+i\sin (\dots )$. The bound-state ansatz breaks down in the strict point-Coulomb model.

For $\kappa =\pm 1$ (the $s_{1/2}$ and $p_{1/2}$ states), the critical $Z$ is $1/\alpha \approx 137.036$. Pomeranchuk and Smorodinsky 1945 first identified the catastrophe and proposed that finite nuclear extent is the physical regularisation. Modern analyses (Greiner-Reinhardt) show that for $137<Z<173$, the spectrum extends smoothly when the Coulomb potential is replaced by the actual finite-nucleus charge distribution; for $Z>173$, the $1s_{1/2}$ binding energy reaches $2m{c}^2$ and the level dives into the negative-energy continuum, producing the predicted but as-yet-unobserved supercritical-field decay of the QED vacuum: spontaneous emission of positrons accompanied by the trapping of $1s$ electrons in the supercritical-field bound state.

The catastrophe is also addressed in the language of self-adjoint extensions: for $Z\alpha <\sqrt{3}/2$ ($Z\le 118$), the Dirac-Coulomb operator on $C_0^{\infty }({\mathbb{R}}^3\setminus \{0\};{\mathbb{C}}^4)$ is essentially self-adjoint and the bound-state problem has a unique answer; for $\sqrt{3}/2<Z\alpha <1$, multiple self-adjoint extensions exist (each parameterised by a $U(1)$ phase at $r=0$), and the physical extension is the distinguished extension that matches finite-nucleus regularisation as $R_N\to 0$.

Hyperfine structure and QED corrections

The Sommerfeld formula is the bare-Dirac-Coulomb result. Real hydrogen has three additional level structures:

1. Hyperfine structure from the proton magnetic moment: each Dirac level splits by $\Delta E_{\text{hf}}\sim (m_e/m_p){\alpha }^4{m}_e{c}^2/n^3$. For the $1s_{1/2}$ ground state, this is the famous 21 cm line ($\nu =1420$ MHz), used by radio astronomy to map neutral hydrogen across the galaxy.

2. Lamb shift from QED radiative corrections: vacuum polarisation, self-energy, and vertex contributions lift the residual $2s_{1/2}$-$2p_{1/2}$ Dirac-Coulomb degeneracy by $\sim 1058$ MHz. Bethe's 1947 non-relativistic calculation reproduced the leading-log term $\sim {\alpha }^5{m}_e{c}^2\log (m_e{c}^2/⟨E_{\text{atomic}}⟩)/(6\pi )$; the full one-loop result agrees with experiment at the four-decimal-place level. The Lamb shift is the natural successor unit 12.16.04.

3. Recoil and finite-nucleus corrections: the proton is not infinitely heavy, and the centre-of-mass correction shifts levels by $m_e/m_p\sim 5\times 10^{-4}$ of the Bohr energy. The finite proton radius $r_p\sim 0.84$ fm shifts levels by $\sim (r_p/a_0{)}^2$, comparable to the Lamb shift for $s$-states. The proton radius puzzle of 2010-2019 was a measured discrepancy in $r_p$ extracted from electronic hydrogen vs. muonic hydrogen spectroscopy; the discrepancy is now believed resolved in favour of the muonic value $r_p=0.842$ fm.

The Dirac-Coulomb spectrum is the foundation on which all of hydrogen-atom precision spectroscopy is built. Modern measurements of the $1s$-$2s$ transition frequency in hydrogen reach a precision of $\sim 10^{-15}$, providing one of the tightest tests of QED and of physics beyond the Standard Model.

Full proof set Master

The proof of the Sommerfeld formula given in the Key derivation section is complete in outline. The technical step that warrants further justification is the index calculation $\gamma =\sqrt{{\kappa }^2-(Z\alpha {)}^2}$ at the origin. Here is the detailed argument.

Near $r=0$, the radial Dirac system to leading order in $r$ has the form

$$\rho g^{\prime }+\kappa g=Z\alpha f+O(\rho ),\rho f^{\prime }-\kappa f=-Z\alpha g+O(\rho ).$$

(The bulk of the mass and energy terms $(m\pm E)\rho$ are subdominant compared to the Coulomb $1/r=\lambda /\rho$ as $\rho \to 0$.) Try $g,f\propto {\rho }^{\gamma }$. Substituting,

$$(\gamma +\kappa )g_0=Z\alpha f_0,(\gamma -\kappa )f_0=-Z\alpha g_0.$$

This is a $2\times 2$ linear system; nonzero solutions require the determinant to vanish:

$$(\gamma +\kappa )(\gamma -\kappa )+(Z\alpha {)}^2=0⟹{\gamma }^2={\kappa }^2-(Z\alpha {)}^2.$$

The regular root is $+\sqrt{{\kappa }^2-(Z\alpha {)}^2}$. For the wavefunction $\psi =g(r)\Omega /r$ to satisfy ${\int }_0^{ϵ}|\psi {|}^2{r}^{2}dr<\infty$, we need $g^2$ integrable at the origin, i.e., $g\sim r^{\gamma }$ with $\gamma >-1/2$, which requires $\gamma >-1/2$, comfortably satisfied for $|\kappa |\ge 1$ as long as $Z\alpha <|\kappa |$. The singular root $-\sqrt{{\kappa }^2-(Z\alpha {)}^2}$ gives $g\sim r^{-\gamma }$ which fails integrability for $\gamma >1/2$ (i.e., for ${\kappa }^2-(Z\alpha {)}^2>1/4$), and is discarded.

The asymptotic analysis at $r\to \infty$ proceeds similarly: the Coulomb term is subdominant, the equations reduce to free Dirac with effective mass $m$ and energy $E$, and the bound-state condition forces ${\Phi }_{\pm }\sim e^{-\lambda r}$. Substituting the ansatz ${\Phi }_{\pm }={\rho }^{\gamma }{e}^{-\rho }{F}_{\pm }(\rho )$ produces a confluent-hypergeometric system for $F_{\pm }$. The termination condition for polynomial $F_{\pm }$ (necessary for ${\int }_0^{\infty }|{\Phi }_{\pm }{|}^{2}d\rho <\infty$) is precisely $\nu =n_r+\gamma$ for non-negative integer $n_r$. Solving for $E$ as in the Key derivation produces the Sommerfeld formula.

The proof of $K$-conservation: from the canonical commutators $[L_i,p_j]=i{ϵ}_{ijk}{p}_k$ and $[\mathbf{\Sigma },\mathbf{\alpha }]=2i\mathbf{\alpha }\times \hat{\mathbf{n}}$ (component-wise),

$$[\mathbf{\alpha }\cdot \mathbf{p},\mathbf{\Sigma }\cdot \mathbf{L}]={\alpha }_i{p}_i{\Sigma }_j{L}_j-{\Sigma }_j{L}_j{\alpha }_i{p}_i={\alpha }_i{\Sigma }_j[p_i,L_j]+[{\alpha }_i,{\Sigma }_j]p_i{L}_j.$$

Using $[p_i,L_j]=i{ϵ}_{ijk}{p}_k$ and $[{\alpha }_i,{\Sigma }_j]=2i{ϵ}_{ijk}{\alpha }_k$:

$$=i{ϵ}_{ijk}{\alpha }_i{\Sigma }_j{p}_k+2i{ϵ}_{ijk}{\alpha }_k{p}_i{L}_j=i\mathbf{\alpha }\cdot (\mathbf{\Sigma }\times \mathbf{p})+2i\mathbf{\alpha }\cdot (\mathbf{p}\times \mathbf{L}).$$

Using $\mathbf{\Sigma }\times \mathbf{p}=\mathbf{p}\mathbf{\Sigma }\cdot \widehat{?}\dots$ and $\mathbf{p}\times \mathbf{L}-\mathbf{L}\times \mathbf{p}=2i\mathbf{p}$, plus the identity $\mathbf{\alpha }\times (\mathbf{p}\times \mathbf{L})=\mathbf{\Sigma }\cdot (\dots )$, one arrives at $[\mathbf{\alpha }\cdot \mathbf{p},\mathbf{\Sigma }\cdot \mathbf{L}+1]=-2\mathbf{\alpha }\cdot \mathbf{p}$. Multiplying by $\beta$ on the right and using $\{\beta ,\mathbf{\alpha }\}=0$ gives $[\mathbf{\alpha }\cdot \mathbf{p},\beta (\mathbf{\Sigma }\cdot \mathbf{L}+1)]=0$. Adding $[\beta m,K]=0$ (since $\beta$ commutes with itself and with $\mathbf{\Sigma }\cdot \mathbf{L}$) and $[-Z\alpha /r,K]=0$ (the central potential commutes with the angular operator $K$) gives $[H,K]=0$. $K$ is conserved.

The Johnson-Lippmann symmetry: ${\hat{A}}_{\text{JL}}=\mathbf{\Sigma }\cdot \hat{\mathbf{r}}K-Z\alpha \beta \mathbf{\Sigma }\cdot \hat{\mathbf{r}}$ commutes with $H$ as follows. The first term: $[\mathbf{\Sigma }\cdot \hat{\mathbf{r}},H_0]$ involves the radial gradient and produces an offset $\sim \mathbf{\alpha }\cdot \widehat{?}$ that combines with $\beta \mathbf{\Sigma }\cdot \hat{\mathbf{r}}$ to cancel. The second term cancels against the Coulomb contribution. The detailed calculation (see Greiner 2000 §9.4) verifies $[{\hat{A}}_{\text{JL}},H]=0$. The action on Dirac spinors is ${\hat{A}}_{\text{JL}}{\psi }_{n,j,\ell =j-1/2,m}\propto {\psi }_{n,j,\ell =j+1/2,m}$, which establishes the residual $\ell$-degeneracy at fixed $(n,j)$.

Connections Master

• Dirac equation and relativistic spin 12.11.01 is the prerequisite framework: the gamma matrices, the Clifford algebra, the four-component spinor, the spin operator, and the magnetic-moment readout $g=2$. The present unit specialises the abstract Dirac equation to the central Coulomb potential and extracts the exact bound-state spectrum.

• Hydrogen atom bound states 12.06.01 is the non-relativistic limit: $E_n=-m(Z\alpha {)}^2/(2n^2)$ is the $O((Z\alpha {)}^2)$ truncation of the Sommerfeld formula. The Schrödinger-Coulomb $\mathrm{SO}(4)$ accidental degeneracy is partially lifted by relativistic kinematics, leaving the residual $(n,j)$-degeneracy of the Dirac-Coulomb spectrum.

• Angular momentum operators and $\mathrm{SU}(2)$ 12.05.01 supply the orbital angular momentum $\mathbf{L}$, the Pauli matrices, the Clebsch-Gordan coefficients needed for spinor spherical harmonics, and the total angular momentum $\mathbf{J}=\mathbf{L}+\mathbf{S}$. The $\kappa$-operator's eigenvalues are computed via the standard ${\mathbf{J}}^2-{\mathbf{L}}^2-{\mathbf{S}}^2$ algebra.

• Lamb shift from one-loop QED 12.16.04 is the natural successor: the Dirac-Coulomb spectrum predicts $2s_{1/2}$-$2p_{1/2}$ degeneracy; QED radiative corrections (vacuum polarisation, electron self-energy, vertex) lift this by $\sim 1058$ MHz. The Lamb-Retherford 1947 measurement of this splitting was the first precision test of QED and the empirical anchor for the 1948 Schwinger / Tomonaga / Feynman / Dyson reformulation.

• Stark and Zeeman effects 12.07.05 use the Dirac-Coulomb wavefunctions as the unperturbed basis for external-field perturbation calculations. The fine-structure splittings derived here are the zeroth-order intervals against which Stark and Zeeman terms are compared.

• Klein-Gordon equation in a Coulomb field 12.11.02 is the spin-0 analogue. The Klein-Gordon Sommerfeld formula is obtained by the substitution $j+1/2\to \ell +1/2$ in the Dirac result. The structural similarity is a manifestation of the underlying integrability of both problems; the difference (spin-0 vs spin-1/2) shows up only in the angular sector.

• Heavy-element relativistic effects 14.04.01 in chemistry trace back to the $(Z\alpha {)}^2$ scaling of the Sommerfeld formula. For $Z=80$ (mercury) and above, the inner $1s$ orbital is genuinely relativistic, and its contraction propagates outward through the orbital basis, explaining gold's colour, mercury's liquidity, and the heavy-element periodic-table irregularities.

• Spin-orbit coupling and fine structure in atomic, molecular, and condensed-matter physics is the universal phenomenology of the $j$-dependence revealed here. The $\mathbf{L}\cdot \mathbf{S}$ Hamiltonian in non-relativistic Pauli theory is the leading $1/m^2$ truncation of the exact Dirac structure; for heavy atoms and for the $L$-shell of transition metals, the full Dirac treatment is needed.

Historical & philosophical context Master

Sommerfeld's 1916 derivation of the fine-structure formula was a tour de force of the old quantum theory ^{[Sommerfeld 1916]}. Starting from the Bohr model, he allowed elliptical orbits (with two quantum numbers $n$ and $k$, the latter the azimuthal quantum number controlling the orbit eccentricity) and added relativistic kinematic corrections. The resulting energy levels matched experiment to four significant figures — an extraordinary success for a theory that, with hindsight, has no rigorous foundation. The formula

$$E_{n,k}^{\text{Som 1916}}=m{c}^2{\left[1+\frac{(Z\alpha {)}^2}{(n-k+\sqrt{k^2-(Z\alpha {)}^2}{)}^2}\right]}^{-1/2}$$

was the spectroscopic gold standard from 1916 to 1928. The puzzle: it worked, but no one knew why. Schrödinger's 1926 non-relativistic wave equation could not reproduce the fine structure; first attempts to relativise it produced the Klein-Gordon equation, which gave the wrong fine-structure pattern (with the $j+1/2$ in the Sommerfeld formula replaced by $\ell +1/2$, a different set of splittings).

Dirac's 1928 equation cracked the puzzle ^{[Dirac 1928]}. Within weeks of the Dirac paper's appearance, Darwin in Cambridge and Gordon in Berlin independently solved the Dirac equation in the Coulomb field ^{[Darwin 1928] [Gordon 1928]} and recovered exactly the Sommerfeld formula, with the azimuthal $k$ identified with the Dirac $|\kappa |=j+1/2$. The half-integer shift that distinguishes the Dirac result from the Klein-Gordon result is the fingerprint of the electron's spin-1/2 nature, derived from the equation rather than assumed.

The conceptual significance ran deeper than the numerical agreement. The Dirac equation showed that the spin-orbit coupling and the relativistic kinetic correction — two contributions Pauli had added by hand to the Schrödinger-Pauli equation to fit the fine structure — were actually two faces of the same underlying object: the Foldy-Wouthuysen-reduced Dirac Hamiltonian. The cancellation between spin-orbit and Darwin terms that produces the $j$-only dependence was not "accidental" but a structural feature of the relativistic theory.

The Sommerfeld-Dirac-Coulomb result was the first instance in physics of an exact, closed-form bound-state spectrum derived from a relativistic quantum equation. The integrability of the problem (the existence of the κ-operator and the Johnson-Lippmann symmetry) is the algebraic reason for the closed form, and a relativistic analogue of the integrability of the non-relativistic Kepler problem (whose closed-form spectrum is explained by the $\mathrm{SO}(4)$ Runge-Lenz symmetry). Few relativistic bound-state problems are integrable in this sense; the Dirac-Coulomb is the principal example, with the Dirac equation in a uniform magnetic field (relativistic Landau levels) as the second.

Pomeranchuk and Smorodinsky's 1945 analysis of the supercritical regime ^{[Pomeranchuk-Smorodinsky 1945]} anticipated by three decades the modern understanding of quantum-electrodynamic vacuum decay in supercritical Coulomb fields. The phenomenon — spontaneous positron emission from a $Z>173$ point-like nucleus — has been the target of heavy-ion-collision experiments at GSI Darmstadt since the 1980s; transient supercritical fields are produced when two heavy ions (Uranium-Uranium, total $Z=184$) collide quasi-elastically, and the predicted positron line at $\sim 250$ keV remains a target of present-day experimental searches.

Lamb and Retherford's 1947 microwave measurement of the $2s_{1/2}$-$2p_{1/2}$ splitting in hydrogen ^{[Lamb-Retherford 1947]} revealed the limits of the Dirac-Coulomb description. The Dirac equation predicts strict degeneracy; the experiment measured a $\sim 1058$ MHz splitting. The discrepancy was the catalyst for the 1947 Shelter Island conference and the rapid development of renormalised QED by Schwinger, Tomonaga, Feynman, and Dyson in 1948-1949. The Lamb shift, alongside the anomalous magnetic moment $a_e$, became the empirical anchor for QED's claim to be the most accurate theory in physics. The 21st-century legacy lives in hydrogen $1s$-$2s$ frequency-comb spectroscopy at $10^{-15}$ relative precision, in muonic-hydrogen spectroscopy that revived the proton-radius puzzle in 2010, and in the perpetual tension between Standard-Model predictions and the muon $g-2$ measurement.

Bibliography Master

Primary literature:

• Sommerfeld, A., "Zur Quantentheorie der Spektrallinien", Ann. Phys. (Leipzig) 51 (1916), 1–94. [Originator: the fine-structure formula from quantised elliptical Bohr orbits.]
• Dirac, P. A. M., "The Quantum Theory of the Electron", Proc. Roy. Soc. A 117 (1928), 610–624; Part II, A 118 (1928), 351–361. [The Dirac equation.]
• Darwin, C. G., "The Wave Equations of the Electron", Proc. Roy. Soc. A 118 (1928), 654–680. [Independent Dirac-Coulomb solution.]
• Gordon, W., "Die Energieniveaus des Wasserstoffatoms nach der Diracschen Quantentheorie des Elektrons", Z. Phys. 48 (1928), 11–14. [Independent Dirac-Coulomb solution.]
• Johnson, M. H. & Lippmann, B. A., "Relativistic Kepler Problem", Phys. Rev. 78 (1950), 329. [Relativistic Runge-Lenz operator; residual $\ell$-degeneracy.]
• Pomeranchuk, I. & Smorodinsky, J., "On the energy levels of systems with $Z>137$", J. Phys. USSR 9 (1945), 97. [Point-nucleus catastrophe; supercritical Coulomb fields.]
• Lamb, W. E. & Retherford, R. C., "Fine Structure of the Hydrogen Atom by a Microwave Method", Phys. Rev. 72 (1947), 241–243. [Discovery of the $2s_{1/2}$-$2p_{1/2}$ splitting beyond the Dirac prediction.]
• Bethe, H. A., "The Electromagnetic Shift of Energy Levels", Phys. Rev. 72 (1947), 339. [Non-relativistic estimate of the Lamb shift.]

Textbooks and monographs:

• Bethe, H. A. & Salpeter, E. E., Quantum Mechanics of One- and Two-Electron Atoms (Springer, 1957), Ch. III. [Canonical reference: Dirac equation in a central field, Sommerfeld formula derivation, radial wavefunctions, fine-structure expansion.]
• Greiner, W., Relativistic Quantum Mechanics: Wave Equations, 3rd ed. (Springer, 2000), Ch. 9. [Detailed textbook treatment with worked examples; supercritical Coulomb regime.]
• Berestetskii, V. B., Lifshitz, E. M., Pitaevskii, L. P., Quantum Electrodynamics, 2nd ed. (Pergamon, 1982), §36. [Landau-Lifshitz Vol. 4 treatment; concise and rigorous.]
• Landau, L. D. & Lifshitz, E. M., Quantum Mechanics: Non-Relativistic Theory, 3rd ed. (Pergamon, 1977), §72. [Fine-structure formula derivation via the non-relativistic expansion.]
• Itzykson, C. & Zuber, J.-B., Quantum Field Theory (McGraw-Hill, 1980), §2-3-3. [Spinor spherical harmonics and the Dirac radial problem.]
• Bjorken, J. D. & Drell, S. D., Relativistic Quantum Mechanics (McGraw-Hill, 1964), Ch. 4. [Standard reference for the single-particle Dirac theory in external fields.]
• Sakurai, J. J., Advanced Quantum Mechanics (Addison-Wesley, 1967), Ch. 3. [Accessible derivation of the radial Dirac system and the Sommerfeld formula.]
• Greiner, W. & Reinhardt, J., Quantum Electrodynamics, 4th ed. (Springer, 2009), Ch. 9. [Supercritical-field decay of the QED vacuum.]
• Eides, M. I., Grotch, H. & Shelyuto, V. A., Theory of Light Hydrogenic Bound States (Springer, 2007). [Modern reference on hydrogen QED corrections beyond the Dirac-Coulomb result.]