Solving the radial Dirac equations: a numerical odyssey

We discuss, in a pedagogical way, how to solve for relativistic wave functions from the radial Dirac equations. After an brief introduction, in Section II we solve the equations for a linear Lorentz scalar potential, V_s(r), that provides for confinement of a quark. The case of massless u and d quarks is treated first, as these are necessarily quite relativistic. We use an iterative procedure to find the eigenenergies and the upper and lower component wave functions for the ground state and then, later, some excited states. Solutions for the massive quarks (s, c, and b) are also presented. In Section III we solve for the case of a Coulomb potential, which is a time-like component of a Lorentz vector potential, V_v(r). We re-derive, numerically, the (analytically well-known) relativistic hydrogen atom eigenenergies and wave functions, and later extend that to the cases of heavier one-electron atoms and muonic atoms. Finally, Section IV finds solutions for a combination of the V_s and V_v potentials. We treat two cases. The first is one in which V_s is the linear potential used in Sec. II and V_v is Coulombic, as in Sec. III. The other is when both V_s and V_v are linearly confining, and we establish when these potentials give a vanishing spin-orbit interaction (as has been shown to be the case in quark models of the hadronic spectrum).

💡 Research Summary

The paper presents a step‑by‑step numerical methodology for solving the radial Dirac equation under several physically relevant potentials, and demonstrates its applicability to both quark confinement models and relativistic atomic systems. After a concise introduction that reviews the Dirac equation in the presence of Lorentz‑scalar (Vₛ) and Lorentz‑vector (Vᵥ) potentials, the authors devote Section II to a purely scalar linear potential Vₛ(r)=k r, which is the standard phenomenological model for quark confinement. They first treat massless u and d quarks, emphasizing that the absence of a rest‑mass term forces the upper (g) and lower (f) radial components to be strongly coupled. An initial guess for the wave‑function near the origin is obtained from a Taylor expansion, and a fourth‑order Runge‑Kutta integrator propagates the solution outward. The eigenenergy is found by a shooting method: the energy is varied until both components decay smoothly to zero at large r, with convergence accelerated by a bisection search combined with Newton‑Raphson corrections. The procedure yields the ground‑state energy and radial profiles, which show a linear rise near the origin and exponential damping beyond a few femtometers.

The same algorithm is then applied to massive quarks (s, c, b). The inclusion of a finite mass term reduces the amplitude of the lower component and shifts the spectrum upward in a non‑linear fashion. By adjusting the string tension k, the authors reproduce the observed ordering of heavy‑meson masses, demonstrating that a value around 0.9 GeV fm⁻¹ provides a reasonable fit. Excited radial states (n = 2, 3) are also obtained with comparable numerical stability, confirming that the method scales to higher quantum numbers.

Section III switches to a time‑like vector potential Vᵥ(r)=−α/r, i.e., the Coulomb interaction. Here the problem reduces to the relativistic hydrogen atom, whose analytic eigenvalues are known:
Eₙ = mc²