Dirac–Coulomb System: Theory and Applications

• The Dirac–Coulomb system is a family of models coupling Dirac dynamics with a 1/r potential across various formulations including external fields, many-body no-pair Hamiltonians, and self-consistent nonlinear approaches.
• Its spectral theory employs orthogonal polynomial methods and perturbative re-expansions to derive bound-state energies, scattering phase shifts, and non-relativistic limits across dimensions.
• The framework provides deep insights into supercritical fields, vacuum polarization, and high-precision computational techniques in both atomic physics and emergent Dirac materials.

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 (Alhaidari et al., 2012, Davydov et al., 2017, Jeszenszki et al., 2021, Baskin et al., 21 Jul 2025, Comech et al., 2012). 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 (Okuma et al., 2014).

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 $$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),$$0 (Baskin et al., 21 Jul 2025). In the massless case, Baskin–Booth–Gell-Redman consider

$$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),$$1

or equivalently

$$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),$$2

for $$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),$$3, which guarantees essential self-adjointness of the Dirac–Coulomb Hamiltonian (Baskin et al., 2021).

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 $$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),$$4, with a two-component Dirac equation

$$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),$$5

or, in polar variables with half-integer $$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),$$6, as a radial first-order system for $$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),$$7 (Davydov et al., 2017). A related 2+1-dimensional construction uses

$$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),$$8

with a short-distance cutoff at $$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),$$9, again leading to partial-wave radial Dirac equations and a Jost-function formulation (Davydov et al., 2017). In 1+1 dimensions, the Hamiltonian

$$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$0

provides a supercritical Dirac–Coulomb model with an explicit short-distance smearing parameter $$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$1 (Davydov et al., 2017).

A self-consistent nonlinear version replaces the external Coulomb field by the Coulomb potential generated by the spinor density. In units $$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$2, the Dirac–Coulomb system may be written as

$$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$3

or, after eliminating $$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$4,

$$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$5

With the solitary-wave ansatz $$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$6, this becomes the nonlinear eigenvalue problem

$$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$7

(Comech et al., 2012).

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 $$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$8-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

$$A_0(t,x)=\frac{\kappa}{|x|}+A_{0,\infty}(t,x),\qquad |\kappa|<1,$$9

and a kinetic-balance relation for $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$0, the Hamiltonian becomes a real symmetric tridiagonal infinite matrix, and the resulting three-term recursion is recognized as the defining recurrence for Pollaczek polynomials (Alhaidari et al., 2012). Darboux analysis of the generating function yields the familiar bound-state spectrum

$\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$1

while the scattering regime gives the phase shift

$\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$2

with $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$3 and $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$4 (Alhaidari et al., 2012).

A complementary approach is the non-relativistic expansion of the Dirac–Coulomb energy in powers of $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$5. Zhou, Yang and Qiao expand the large and small components and the eigenvalue as

$\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$6

derive iterative equations for $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$7, $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$8, and $\supp A_j,\supp A_{0,\infty}\subset\{t\in(t_-,t_+)\}$9, and obtain explicit operator formulas through $t$0 (Zhou et al., 2024). For one electron, the iterative equations are numerically carried to order $t$1 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 $t$2 (Zhou et al., 2024).

The same expansion framework also incorporates the non-retarded Breit interaction. The $t$3 order correction to the Dirac Coulomb energy and non-retarded Breit interaction corresponds precisely to the $t$4 order relativistic correction, while higher even orders contribute at $t$5 and represent the contributions from all Coulomb photons and single transverse photons under the non-retarded approximation (Zhou et al., 2024).

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

$t$6

where $t$7 subtracts the linear-in-$t$8 first Born part of the Green function. This subtraction removes all ultraviolet divergences and enforces

$t$9

Whenever a bound level dives into the lower continuum at $H(t)$0, the additional nonperturbative vacuum-shell density

$H(t)$1

appears and the total induced charge jumps by $H(t)$2 (Davydov et al., 2017). The vacuum energy admits a phase-shift representation plus a $H(t)$3 counterterm,

$H(t)$4

and its large-$H(t)$5 behavior depends on the sign of $H(t)$6: for $H(t)$7, $H(t)$8 eventually becomes large and negative, asymptotically $H(t)$9 (Davydov et al., 2017).

In 2+1 dimensions, vacuum polarization is formulated through partial waves labeled by $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$0. The exact vacuum-polarization density can be written as a contour integral of the full Green function and then as

$C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$1

with each $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$2 expressed as an integral over $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$3 (Davydov et al., 2017). The renormalized density is

$C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$4

where subtraction of $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$5 removes the only divergent piece. By construction, $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$6, all $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$7 are finite, the total induced charge vanishes for $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$8, and when a pair of bound levels dives at $C_c^\infty(\mathbb R^3\setminus\{0\};\mathbb C^4)$9, each contributes $H^1(\mathbb R^3;\mathbb C^4)$0 to the total induced charge in exact agreement with Furry’s theorem (Davydov et al., 2017).

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 $H^1(\mathbb R^3;\mathbb C^4)$1 and bound-state energies $H^1(\mathbb R^3;\mathbb C^4)$2. Each partial contribution $H^1(\mathbb R^3;\mathbb C^4)$3 is finite, but the series over $H^1(\mathbb R^3;\mathbb C^4)$4 diverges linearly because for large $H^1(\mathbb R^3;\mathbb C^4)$5,

$H^1(\mathbb R^3;\mathbb C^4)$6

The unique divergent piece is proportional to $H^1(\mathbb R^3;\mathbb C^4)$7, and renormalization is imposed channel by channel by matching the $H^1(\mathbb R^3;\mathbb C^4)$8 limit to the first Born approximation $H^1(\mathbb R^3;\mathbb C^4)$9 (Davydov et al., 2017). In the overcritical regime $U(t,s)$0, only a finite set of channels contributes shells, with $U(t,s)$1, and the renormalized energy scales as

$U(t,s)$2

with $U(t,s)$3 (Davydov et al., 2017).

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 $U(t,s)$4, the number of vacuum shells grows roughly as $U(t,s)$5, and correspondingly

$U(t,s)$6

Since the classical self-energy of a uniformly charged sphere is $U(t,s)$7, qualitative arguments are presented in favor of the possibility for complete screening for $U(t,s)$8 in 3+1 D (Davydov et al., 2017).

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 $U(t,s)$9 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 $(i\partial_t+H(t))U(t,s)=0$0, the causal propagator $(i\partial_t+H(t))U(t,s)=0$1 and the two-point functions

$(i\partial_t+H(t))U(t,s)=0$2

satisfy $(i\partial_t+H(t))U(t,s)=0$3, $(i\partial_t+H(t))U(t,s)=0$4, and $(i\partial_t+H(t))U(t,s)=0$5 (Baskin et al., 21 Jul 2025). On the blown-up space $(i\partial_t+H(t))U(t,s)=0$6, the relevant regularity is measured in Melrose’s $(i\partial_t+H(t))U(t,s)=0$7-Sobolev spaces. The main propagation theorem states that if $(i\partial_t+H(t))U(t,s)=0$8 and $(i\partial_t+H(t))U(t,s)=0$9, then

$$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),$$00

and absence of singularity propagates along diffractive bicharacteristics, including broken null rays meeting $$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),$$01 (Baskin et al., 21 Jul 2025). One consequence is that the in- and out-Dirac–Coulomb vacua are Hadamard on all of $$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),$$02. Another is that for any two Hadamard states, the relative charge density

$$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),$$03

is well-defined as a locally integrable function including near $$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),$$04, with

$$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),$$05

(Baskin et al., 21 Jul 2025).

For the massless Dirac–Coulomb equation with $$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),$$06, the long-time behavior is described by a polyhomogeneous expansion on a suitable blow-up of future null infinity. After the rescaling $$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),$$07, the forward Friedlander radiation field

$$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),$$08

is $$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),$$09 and admits an expansion

$$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),$$10

while the slowest decaying off-light-cone term is $$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),$$11 (Baskin et al., 2021). 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” (Baskin et al., 2021).

Time dependence can also enter through moving Coulomb centers. Cacciafesta–de Suzzoni–Noja study a Dirac electron field coupled to classically moving point nuclei,

$$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),$$12

together with Newton equations for the nuclei. Under $$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),$$13, well-separated centers, and bounded $$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),$$14-norm of the trajectories, they construct a unique two-parameter unitary flow

$$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),$$15

and prove local existence of the coupled system for $$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),$$16, $$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),$$17 (Cacciafesta et al., 2017).

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 $$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),$$18, 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

$$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),$$19

Because this operator couples the electronic, Brown–Ravenhall, and positronic continua, Ferenc, Jeszenszki, and Mátyus introduce the no-pair projector

$$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),$$20

and solve the projected problem variationally in an explicitly correlated Gaussian basis with restricted kinetic balance (Jeszenszki et al., 2021). Three projection algorithms are implemented: Complex-Coordinate Rotation, Energy-Cutting, and “Punching.” For low $$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),$$21 systems, the Hermitian energy-cutting procedure reproduces CCR results to $$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),$$22 (Jeszenszki et al., 2021).

The working generalized eigenvalue problem has dimension $$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),$$23. In the QUANTEN implementation, basis parameters are first optimized at the non-relativistic level and then refined for the no-pair Dirac–Coulomb energy (Jeszenszki et al., 2021). Reported benchmarks include helium $$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),$$24,

$$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),$$25

with CCR- and cutting-projected energies agreeing to $$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),$$26; $$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),$$27 at $$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),$$28 bohr,

$$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),$$29

and $$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),$$30 at $$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),$$31 bohr,

$$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),$$32

(Jeszenszki et al., 2021). For isoelectronic atoms with $$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),$$33, the no-pair Dirac–Coulomb energies lie systematically above the $$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),$$34-only values by 10–200 nE, and the $$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),$$35 two-photon term recovers most of this difference for $$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),$$36, while for $$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),$$37 significant higher-order terms are required (Jeszenszki et al., 2021).

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(t)=\sum_{j=1}^3\alpha^j\bigl(i^{-1}\partial_{x^j}+A_j(t,x)\bigr)+m\,\beta-A_0(t,x),$$38

acting on a sixteen-component spinor (Ferenc et al., 2023). In the zero-total-momentum frame, the spatial basis is taken as spherical Gaussians $$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),$$39, kinetic balance is imposed by a block-diagonal transformation $$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),$$40, and the generalized eigenvalue problem is solved with basis size up to $$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),$$41, giving sub-ppb convergence of the lowest eigenvalue (Ferenc et al., 2023). By fitting the $$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),$$42-dependence,

$$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),$$43

the fitted coefficients for positronium, hydrogen, muonium, and muonic hydrogen are stated to agree with known analytic nrQED coefficients within the variational-convergence uncertainty (Ferenc et al., 2023).

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

$$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),$$44

arises as the Euler–Lagrange equation of the energy-minus-frequency functional

$$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),$$45

under the charge constraint $$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),$$46 (Comech et al., 2012). In the small-charge, nonrelativistic regime with $$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),$$47, the no-node branch bifurcates from the unique positive, radial ground state of the Choquard equation, with scaling

$$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),$$48

Linearization leads to a real-linear spectral problem $$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),$$49, and for $$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),$$50 sufficiently close to $$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),$$51 the point spectrum has no eigenvalues with $$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),$$52; the argument uses a limiting-absorption principle, exclusion of bifurcation at $$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),$$53, and the Vakhitov–Kolokolov condition $$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),$$54 inherited from the Choquard limit (Comech et al., 2012).

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 $$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),$$55,

$$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),$$56

with

$$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),$$57

and study the long-range Coulomb interaction by Wilsonian renormalization group (Okuma et al., 2014). At one loop, the self-energy produces

$$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),$$58

After adding the Zeeman term $$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),$$59, the vertex correction gives

$$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),$$60

(Okuma et al., 2014). 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 $$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),$$61-factor enhancement, which is described there as a unique property of the surface Dirac system. Numerically, the $$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),$$62-factor enhancement dominates over the $$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),$$63-driven suppression in the spin susceptibility, so $$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),$$64 is overall increased by the long-range Coulomb interaction at low $$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),$$65 (Okuma et al., 2014).

The contrast drawn in that paper is explicit: in graphene the Pauli matrices in $$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),$$66 act on sublattice pseudo-spin and therefore commute with the real-spin Zeeman operator, so no vertex correction to $$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),$$67 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 $$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),$$68 (Okuma et al., 2014).

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.