---
title: 'Dirac–Coulomb System: Theory and Applications'
url: https://www.emergentmind.com/topics/dirac-coulomb-system
type: topic
---

# Dirac–Coulomb System: Theory and Applications

In the literature summarized here, the Dirac–Coulomb system appears as relativistic Dirac dynamics coupled to a Coulomb potential in several distinct but closely related senses: the standard external-field bound-state and scattering problem; supercritical external-field QED in 1+1 and 2+1 dimensions; no-pair projected many-electron and two-body Hamiltonians for atoms and molecules; time-dependent Dirac propagation in the presence of Coulomb singularities; and self-consistent nonlinear models in which the Coulomb field is generated by the Dirac density itself [1204.6493] [1712.02704] [2110.06638] [2507.16085] [1207.2870]. The same label also covers effective condensed-matter realizations, such as a two-dimensional massless Dirac surface state with unscreened long-range Coulomb interaction on a three-dimensional topological insulator surface [1412.0797].

## 1. Operator formulations and dimensional variants

A canonical 3+1-dimensional formulation is the Dirac Hamiltonian on Minkowski space \(M=\mathbb R_t\times\mathbb R^3_x\),
\[
H(t)=\sum_{j=1}^3\alpha^j\bigl(i^{-1}\partial_{x^j}+A_j(t,x)\bigr)+m\,\beta-A_0(t,x),
\]
with
\[
A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,
\]
and \(\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}\). For each fixed \(t\), \(H(t)\) is essentially self-adjoint on \(C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)\), its closure has domain \(H^1(\mathbb R^3;\mathbb C^4)\), and one may define the unitary propagator \(U(t,s)\) solving \((i\partial_t+H(t))U(t,s)=0\) with \(U(s,s)=\mathrm Id\) [2507.16085]. In the massless case, Baskin–Booth–Gell-Redman consider
\[
D=i\gamma^\mu\partial_\mu-\frac{Z}{r},
\]
or equivalently
\[
{}_{Z/r}=\gamma^0\Bigl(\partial_t+i\,\frac{Z}{r}\Bigr)+\sum_{j=1}^3\gamma^j\partial_{x_j},
\]
for \(|Z|<\tfrac12\), which guarantees essential self-adjointness of the Dirac–Coulomb Hamiltonian [2112.06111].

Lower-dimensional realizations are structurally similar but technically distinct. In 2+1 dimensions, the external Coulomb source can be taken as the projection onto the plane of a uniformly charged sphere of radius \(R\), with a two-component Dirac equation
\[
(-i\,\gamma^i\partial_i+V(r)+\gamma^0)\psi(\vec r)=\epsilon\,\psi(\vec r),
\]
or, in polar variables with half-integer \(m_j\), as a radial first-order system for \(f_1,f_2\) [1712.02703]. A related 2+1-dimensional construction uses
\[
A_0^{ext}(r)=V(r)=-Z\alpha/r
\]
with a short-distance cutoff at \(r=R\), again leading to partial-wave radial Dirac equations and a Jost-function formulation [1712.02704]. In 1+1 dimensions, the Hamiltonian
\[
H=-\,i\,\alpha\,\frac{d}{dx}+\beta\,m+V(x),\qquad V(x)=-\frac{Z\alpha}{|x|+a},
\]
provides a supercritical Dirac–Coulomb model with an explicit short-distance smearing parameter \(a>0\) [1709.04239].

A self-consistent nonlinear version replaces the external Coulomb field by the Coulomb potential generated by the spinor density. In units \(\hbar=c=1\), the Dirac–Coulomb system may be written as
\[
i\,\partial_t\psi=D_m\psi+e\,\phi\,\psi,\qquad \Delta\phi=e\,\psi^*\psi,
\]
or, after eliminating \(\phi\),
\[
i\,\partial_t\psi=D_m\psi-e^2\,(|x|^{-1}\ast\psi^*\psi)\psi.
\]
With the solitary-wave ansatz \(\psi(x,t)=e^{-i\omega t}\phi(x)\), this becomes the nonlinear eigenvalue problem
\[
(D_m-\omega)\phi-e^2\bigl(|x|^{-1}\ast|\phi|^2\bigr)\phi=0
\]
[1207.2870].

This range of formulations suggests that the phrase “Dirac–Coulomb system” functions less as a single equation than as a family of relativistic Coulomb-coupled models whose common structure is the interaction between Dirac kinematics and a \(1/r\)-type potential.

## 2. Spectral theory, tridiagonalization, and non-relativistic expansion

For the standard one-electron problem, Alhaidari–Bahlouli–Ismail obtain a symmetric tridiagonal matrix representation of the Dirac–Coulomb operator in a square-integrable Laguerre spinor basis. With
\[
\phi_n^+(r)=A_n(\omega r)^{\gamma+1}e^{-\omega r/2}L_n^{2\gamma+1}(\omega r),\qquad \gamma=\sqrt{\kappa^2-(Z\alpha)^2},
\]
and a kinetic-balance relation for \(\phi_n^-(r)\), the Hamiltonian becomes a real symmetric tridiagonal infinite matrix, and the resulting three-term recursion is recognized as the defining recurrence for Pollaczek polynomials [1204.6493]. Darboux analysis of the generating function yields the familiar bound-state spectrum
\[
E_{n,\kappa}
=
m c^2
\Bigl[1+\frac{(Z\alpha)^2}{(n+\gamma+1)^2}\Bigr]^{-1/2},
\]
while the scattering regime gives the phase shift
\[
\delta_\kappa(\epsilon)
=
\arg\Gamma\!\bigl(\gamma+1+i\eta\bigr)
-\eta\ln(2kR_0)
+\Bigl(\gamma+\tfrac12\Bigr)\Bigl(\theta-\tfrac\pi2\Bigr),
\]
with \(k=\sqrt{\epsilon^2-1}\) and \(\eta=\alpha/k\) [1204.6493].

A complementary approach is the non-relativistic expansion of the Dirac–Coulomb energy in powers of \(\alpha^2\). Zhou, Yang and Qiao expand the large and small components and the eigenvalue as
\[
E=mc^2+\sum_{k=1}^\infty E^{(2k)}\alpha^{2k},\qquad
\varphi=\sum_{k=1}^\infty\varphi^{(2k)},\qquad
\chi=\sum_{k=1}^\infty\chi^{(2k+1)},
\]
derive iterative equations for \(E^{(2n)}\), \(\varphi^{(2n)}\), and \(\chi^{(2n+1)}\), and obtain explicit operator formulas through \(E^{(8)}\) [2409.06195]. For one electron, the iterative equations are numerically carried to order \(\alpha^{20}\) for hydrogen and are reported to converge rapidly to the analytical results of the hydrogen atom. For the two-electron Dirac–Coulomb system, the same paper presents iterative equations for high-order energy corrections and ground-state energy corrections up to order \(\alpha^8\) [2409.06195].

The same expansion framework also incorporates the non-retarded Breit interaction. The \(\alpha^4\) order correction to the Dirac Coulomb energy and non-retarded Breit interaction corresponds precisely to the \(\alpha^4\) order relativistic correction, while higher even orders contribute at \(\alpha^6,\alpha^8,\dots\) and represent the contributions from all Coulomb photons and single transverse photons under the non-retarded approximation [2409.06195].

Taken together, these results show that the spectral theory of the Dirac–Coulomb system admits both an exact algebraic formulation in terms of orthogonal polynomials and a systematic perturbative re-expansion that interfaces directly with Breit–Pauli and nrQED-type structures.

## 3. Supercritical fields, vacuum polarization, and shell effects

In supercritical external fields, the Dirac–Coulomb system ceases to be adequately described by perturbative vacuum polarization alone. In 1+1 dimensions, the renormalized vacuum density is written as
\[
\rho_{\rm vac}^R(x)=\rho_{\rm vac}^{(1)}(x)+\rho_{\rm vac}^{(3+)}(x),
\]
where \(\rho_{\rm vac}^{(3+)}\) subtracts the linear-in-\(Z\) first Born part of the Green function. This subtraction removes all ultraviolet divergences and enforces
\[
\int\rho_{\rm vac}^R(x)\,dx=0\qquad (Z<Z_{cr,1}).
\]
Whenever a bound level dives into the lower continuum at \(Z=Z_{cr,k}\), the additional nonperturbative vacuum-shell density
\[
\Delta\rho_{\rm vac}(x)=-\,|e|\,|\psi_n(x)|^2
\]
appears and the total induced charge jumps by \(-|e|\) [1709.04239]. The vacuum energy admits a phase-shift representation plus a \(Z^2\) counterterm,
\[
E_{\rm vac}^R(Z)=E_{\rm vac}(Z)+\lambda(a)Z^2,
\]
and its large-\(Z\) behavior depends on the sign of \(\lambda(a)\): for \(\lambda(a)<0\), \(E_{\rm vac}^R\) eventually becomes large and negative, asymptotically \(-|\lambda|Z^2\) [1709.04239].

In 2+1 dimensions, vacuum polarization is formulated through partial waves labeled by \(m_j=\pm \tfrac12,\pm\tfrac32,\dots\). The exact vacuum-polarization density can be written as a contour integral of the full Green function and then as
\[
\rho_{VP}(r)=2\sum_{m_j=1/2,3/2,\dots}\rho_{VP,|m_j|}(r),
\]
with each \(\rho_{VP,|m_j|}\) expressed as an integral over \(\operatorname{Re}\operatorname{Tr}G_{|m_j|}(r,r;iy)\) [1712.02704]. The renormalized density is
\[
\rho_{VP}^{ren}(r)=2\Bigl[\rho_{VP}^{(1)}(r)+\sum_{m_j=1/2,3/2,\dots}\rho_{VP,|m_j|}^{(3+)}(r)\Bigr],
\]
where subtraction of \(2\,\mathrm{Tr}\,G^{(1)}_{m_j}\) removes the only divergent piece. By construction, \(\int d^2r\,\rho_{VP}^{(1)}_{ }=0\), all \(\rho^{(3+)}\) are finite, the total induced charge vanishes for \(Z<Z_{cr,1}\), and when a pair of bound levels dives at \(Z=Z_{cr,k}\), each contributes \(-2|e|\) to the total induced charge in exact agreement with Furry’s theorem [1712.02704].

The vacuum energy in the 2+1-dimensional supercritical system is defined by a mode sum with free-vacuum subtraction and rewritten through scattering phase shifts \(\delta_{tot,m_j}(k)\) and bound-state energies \(\epsilon_{n,\pm m_j}\). Each partial contribution \(E_{VP,m_j}\) is finite, but the series over \(m_j\) diverges linearly because for large \(|m_j|\),
\[
E_{VP,m_j}\sim \frac1\pi\int_0^\infty V^2(r)\,dr + O(1/|m_j|^3).
\]
The unique divergent piece is proportional to \(Z^2\), and renormalization is imposed channel by channel by matching the \(Z\to0\) limit to the first Born approximation \(E_{VP}^{(1)}\) [1712.02703]. In the overcritical regime \(Z>Z_{cr,1}\), only a finite set of channels contributes shells, with \(m_j^{max}(Z)\sim Z\), and the renormalized energy scales as
\[
E_{VP}^{ren}(Z)\sim -\eta_{eff}\,\frac{Z^3}{R},\qquad Z\gg Z_{cr,1},
\]
with \(\eta_{eff}>0\) [1712.02703].

The same analysis is used to motivate a 3+1-dimensional screening picture. The paper states that in 3+1 D the analogous expansion carries an extra degeneracy factor \(2j+1\), the number of vacuum shells grows roughly as \(N(Z)\sim Z^{3.2\ldots 3.3}\), and correspondingly
\[
E_{VP,3+1}^{ren}(Z)\sim -\eta_{3+1}\,\frac{Z^4}{R}.
\]
Since the classical self-energy of a uniformly charged sphere is \(E_{class}=Z^2/(2R)\), qualitative arguments are presented in favor of the possibility for complete screening for \(Z\gg Z_{cr,1}\) in 3+1 D [1712.02703].

A plausible implication is that “vacuum shells” are the central nonperturbative organizing principle of the supercritical Dirac–Coulomb problem: they govern both induced charge and the transition from perturbative \(Z^2\) behavior to strongly negative renormalized vacuum energy.

## 4. Propagators, wave-front structure, and time-dependent singularities

A major recent development is the microlocal analysis of Dirac–Coulomb propagators and states in the presence of the Coulomb singularity. For the time-dependent Hamiltonian with \(A_0(t,x)=\kappa/|x|+A_{0,\infty}(t,x)\), the causal propagator \(S\) and the two-point functions
\[
\Lambda^-(t_1,t_2)=U(t_1,t_0)\gamma U(t_0,t_2),\qquad
\Lambda^+(t_1,t_2)=U(t_1,t_0)(1-\gamma)U(t_0,t_2)
\]
satisfy \(D\Lambda^\pm=0\), \(\Lambda^\pm D=0\), and \(\Lambda^++\Lambda^-=S\) [2507.16085]. On the blown-up space \(M=[\mathbb R_t\times\mathbb R^3;\{r=0\}]\), the relevant regularity is measured in Melrose’s \(b\)-Sobolev spaces. The main propagation theorem states that if \(Du=0\) and \(u\in H_b^{1,-\infty}\), then
\[
\mathrm{WF}_b^{1,\infty}(u)\subset\Sigma,
\]
and absence of singularity propagates along diffractive bicharacteristics, including broken null rays meeting \(r=0\) [2507.16085]. One consequence is that the in- and out-Dirac–Coulomb vacua are Hadamard on all of \(M\setminus\{r=0\}\). Another is that for any two Hadamard states, the relative charge density
\[
\rho(x)=\operatorname{Tr}_{\mathbb C^4}\bigl[(\gamma-\gamma_{ref})(x,x)\bigr]
\]
is well-defined as a locally integrable function including near \(r=0\), with
\[
\rho\in L^1_{loc}(\mathbb R^3)\cap C^\infty(\mathbb R^3\setminus\{0\}),\qquad
\rho(r)=O(r^{-2-\varepsilon})
\]
[2507.16085].

For the massless Dirac–Coulomb equation with \(|Z|<\tfrac12\), the long-time behavior is described by a polyhomogeneous expansion on a suitable blow-up of future null infinity. After the rescaling \(u=\rho^{-1-iZ}\psi\), the forward Friedlander radiation field
\[
\mathcal R_+[\psi](s,\theta)=\lim_{r\to\infty}r^{1+iZ}\psi(s+r,r\theta)
\]
is \(C^\infty(\mathbb R_s\times S^2_\theta)\) and admits an expansion
\[
\mathcal R_+[\psi](s,\theta)\sim
s^{iZ}\sum_{j=1}^\infty\sum_{k=1}^\infty s^{-j-\sqrt{k^2-Z^2}}\,a_{j,k}(\theta),
\]
while the slowest decaying off-light-cone term is \((t-r)^{-(1+\sqrt{1-Z^2})}\) [2112.06111]. The same paper proves a diffractive propagation statement formulated as “there is no diffraction-generated singularity at the origin: no ‘new wave-front’ is created when an incoming singularity transits the Coulomb singularity” [2112.06111].

Time dependence can also enter through moving Coulomb centers. Cacciafesta–de Suzzoni–Noja study a Dirac electron field coupled to classically moving point nuclei,
\[
i\partial_tu=(D_0+\beta)u-\sum_{k=1}^N\frac{Z_k}{|x-q_k(t)|}u+\bigl(|x|^{-1}\ast|u|^2\bigr)u,
\]
together with Newton equations for the nuclei. Under \(|Z_k|<\sqrt3/2\), well-separated centers, and bounded \(W^{2,1}\)-norm of the trajectories, they construct a unique two-parameter unitary flow
\[
U_q(t,s):H^\sigma\to H^\sigma,\qquad 0\le s\le t\le T,\quad \sigma\in[0,2),
\]
and prove local existence of the coupled system for \(u(0)\in H^\sigma\), \(1\le \sigma<\tfrac32\) [1709.05317].

Taken together, these results rigorously justify the usual physical picture in external-field Dirac–Coulomb theory: the Coulomb singularity does not destroy Hadamard structure away from \(r=0\), propagation through the singularity can be controlled microlocally, and time-dependent singular configurations still admit a well-posed propagator theory.

## 5. No-pair many-body Hamiltonians and high-precision computation

In atomic and molecular applications, the many-electron Dirac–Coulomb operator is
\[
\mathcal H_{DC}
=
\sum_{i=1}^N
\Bigl[c\,\boldsymbol\alpha_i\!\cdot\!\mathbf p_i+\beta_i c^2+U(\mathbf r_i)\Bigr]
+
\sum_{i<j}\frac1{r_{ij}},
\qquad
U(\mathbf r)= -\sum_{A=1}^{nuc}\frac{Z_A}{|\mathbf r-\mathbf R_A|}.
\]
Because this operator couples the electronic, Brown–Ravenhall, and positronic continua, Ferenc, Jeszenszki, and Mátyus introduce the no-pair projector
\[
\Lambda^+=\sum_{\mu\in E^+}|\Psi_{0,\mu}\rangle\langle\tilde\Psi_{0,\mu}|,
\qquad
\bar{\mathcal H}_{DC}=\Lambda^+\mathcal H_{DC}\Lambda^+,
\]
and solve the projected problem variationally in an explicitly correlated Gaussian basis with restricted kinetic balance [2110.06638]. Three projection algorithms are implemented: Complex-Coordinate Rotation, Energy-Cutting, and “Punching.” For low \(Z\le4\) systems, the Hermitian energy-cutting procedure reproduces CCR results to \(\sim10^{-11}E\) [2110.06638].

The working generalized eigenvalue problem has dimension \(16N_b\). In the QUANTEN implementation, basis parameters are first optimized at the non-relativistic level and then refined for the no-pair Dirac–Coulomb energy [2110.06638]. Reported benchmarks include helium \(1\,{}^1S\),
\[
E_{nr}=-2.903\,724\,377\,E,\qquad
E_{DC}^{proj}=-2.903\,856\,631\,E,
\]
with CCR- and cutting-projected energies agreeing to \(10^{-11}E\); \(\mathrm H_2\) at \(R_{pp}=1.4\) bohr,
\[
E_{nr}=-1.174\,475\,714\,E,\qquad
E_{DC}^{proj}=-1.174\,489\,754\,E;
\]
and \(\mathrm H_3^+\) at \(R_{pp}=1.65\) bohr,
\[
E_{DC}^{proj}=-1.343\,850\,527\,E
\]
[2110.06638]. For isoelectronic atoms with \(Z=1\ldots26\), the no-pair Dirac–Coulomb energies lie systematically above the \(\alpha^2\)-only values by 10–200 nE, and the \(\alpha^3\) two-photon term recovers most of this difference for \(Z\le4\), while for \(Z>4\) significant higher-order terms are required [2110.06638].

A two-body, pre-Born–Oppenheimer extension starts from the Salpeter–Sucher form of the Bethe–Salpeter equation and the no-pair Dirac–Coulomb–Breit Hamiltonian
\[
H=\Lambda_{++}\bigl[h_D(1)+h_D(2)+V_C(r_{12})+V_B(r_{12})\bigr]\Lambda_{++},
\]
acting on a sixteen-component spinor [2301.13477]. In the zero-total-momentum frame, the spatial basis is taken as spherical Gaussians \(f_i(r)=e^{-\zeta_i r^2}\), kinetic balance is imposed by a block-diagonal transformation \(X\), and the generalized eigenvalue problem is solved with basis size up to \(50\), giving sub-ppb convergence of the lowest eigenvalue [2301.13477]. By fitting the \(\alpha\)-dependence,
\[
F(\alpha)=E_0+\alpha^2\varepsilon_2+\alpha^3\varepsilon_3+\alpha^4\ln(\alpha)\varepsilon_4'+\alpha^4\varepsilon_4,
\]
the fitted coefficients for positronium, hydrogen, muonium, and muonic hydrogen are stated to agree with known analytic nrQED coefficients within the variational-convergence uncertainty [2301.13477].

This computational literature suggests a distinct modern meaning of the Dirac–Coulomb system: not only a solvable one-center equation, but also a reference Hamiltonian for precision spectroscopy, variational relativistic quantum chemistry, and non-perturbative baselines for perturbative QED corrections.

## 6. Self-consistent Coulomb fields and condensed-matter realizations

In the nonlinear self-consistent setting, solitary waves of the Dirac–Coulomb system are analyzed as polarons. The stationary problem
\[
(D_m-\omega)\phi-e^2\bigl(|x|^{-1}\ast|\phi|^2\bigr)\phi=0
\]
arises as the Euler–Lagrange equation of the energy-minus-frequency functional
\[
E[\phi]
=
\int \phi^*(D_m-\omega)\phi\,dx
-\frac{e^2}{2}\iint
\frac{\phi^*(x)\phi(x)\phi^*(y)\phi(y)}{|x-y|}\,dx\,dy
\]
under the charge constraint \(\int |\phi|^2dx=Q\) [1207.2870]. In the small-charge, nonrelativistic regime with \(\epsilon^2=m^2-\omega^2\ll1\), the no-node branch bifurcates from the unique positive, radial ground state of the Choquard equation, with scaling
\[
\phi_e(x)\approx \epsilon^2u(\epsilon x),\qquad
\phi_p(x)=O(\epsilon^3),\qquad
\omega=\sqrt{m^2-\epsilon^2}.
\]
Linearization leads to a real-linear spectral problem \(J L(\omega)\xi=\lambda\xi\), and for \(\omega\) sufficiently close to \(m\) the point spectrum has no eigenvalues with \(\operatorname{Re}\lambda>0\); the argument uses a limiting-absorption principle, exclusion of bifurcation at \(\lambda=\pm2mi\), and the Vakhitov–Kolokolov condition \(dQ/d\omega<0\) inherited from the Choquard limit [1207.2870].

A different effective realization occurs on the surface of a three-dimensional topological insulator. Okuma and Ogata consider the two-dimensional massless Dirac surface state at chemical potential \(\mu=0\),
\[
H_0=\int d^2k\,\Psi^\dagger(k)\,[\,v_F(k_x\sigma_y-k_y\sigma_x)\,]\Psi(k),
\qquad
H_{int}=\frac12\sum_qV(q)\,n(-q)n(q),
\]
with
\[
V(q)=\frac{2\pi e^2}{\epsilon |q|},
\]
and study the long-range Coulomb interaction by Wilsonian renormalization group [1412.0797]. At one loop, the self-energy produces
\[
\frac{dv_F}{dl}=\frac{e^2}{4\epsilon},
\qquad
v_F(l)=v_F(0)+\frac{e^2}{4\epsilon}l.
\]
After adding the Zeeman term \(H_Z=(g\mu_B/2)B_z\sigma_z\), the vertex correction gives
\[
\frac{dg}{dl}=g\,\frac{e^2}{2\epsilon}\,\frac1{v_F(l)},
\qquad
g(l)=g(0)\Bigl[\frac{v_F(l)}{v_F(0)}\Bigr]^2
=
g(0)\Bigl[1+\frac{e^2}{4\epsilon v_F(0)}\,l\Bigr]^2
\]
[1412.0797]. The remarkable feature is that the Dirac Hamiltonian contains not pseudo spin but real spin Pauli matrices. Because of this feature, the paper finds the \(g\)-factor enhancement, which is described there as a unique property of the surface Dirac system. Numerically, the \(g\)-factor enhancement dominates over the \(v_F\)-driven suppression in the spin susceptibility, so \(\chi_s(T)\) is overall increased by the long-range Coulomb interaction at low \(T\) [1412.0797].

The contrast drawn in that paper is explicit: in graphene the Pauli matrices in \(H_0\) act on sublattice pseudo-spin and therefore commute with the real-spin Zeeman operator, so no vertex correction to \(g\) appears at one loop, whereas on a topological-insulator surface the Pauli matrices are the true electron spins and the Coulomb-induced self-energy and vertex diagrams “talk” to \(\sigma_z\) [1412.0797].

This suggests that the Dirac–Coulomb label extends naturally from relativistic atomic Hamiltonians to emergent Dirac media whenever Coulomb interactions remain unscreened or self-consistently generated, and that the same formal ingredients—Dirac kinematics, Coulomb kernels, spectral thresholds, and renormalization—recur across high-energy, mathematical, chemical, and condensed-matter settings.

Source: https://www.emergentmind.com/topics/dirac-coulomb-system