Two-Center Dirac Equation

• The two-center Dirac equation is the relativistic eigenvalue problem for a spin-½ electron in the Coulomb field of two fixed nuclei within the Born–Oppenheimer approximation.
• It employs multiple coordinate representations and symmetry reductions to address nonseparability and overcome spurious-state contamination in numerical methods.
• Advanced discretization approaches such as A-DKB and minmax-FEM achieve high-precision benchmarks, resolving Coulomb singularities and enabling QED corrections in heavy-ion studies.

Searching arXiv for recent and foundational papers on the two-center Dirac equation to ground the article in the current literature. The two-center Dirac equation is the relativistic eigenvalue or evolution equation for a spin-$\tfrac12$ electron in the field of two fixed nuclei, usually within the Born–Oppenheimer approximation. In its stationary form, it is written as

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$

with a two-center Coulomb potential

$$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$

In relativistic atomic or natural units, this becomes

$$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$

or, equivalently in atomic units, $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$ (Solovyev et al., 2023, Kullie, 2024). Unlike the one-center Dirac–Coulomb problem, the two-center problem is generically nonseparable in spherical coordinates because the potential is not central; axial symmetry reduces the dimensionality, but does not restore a single radial equation. The equation is central to relativistic molecular structure, quasi-molecular collision theory, supercritical heavy-ion physics, and precision spectroscopy of one-electron molecular ions such as ${\rm H}_2^+$ (Kullie et al., 2022, McConnell et al., 2012).

1. Operator structure and spectral setting

The standard geometry places the nuclei on the internuclear axis,

$$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$

so the system is axially symmetric about the $z$-axis. The conserved quantum number is the projection of total electronic angular momentum on the molecular axis,

$J_z\psi=m_J\psi,$

or, in prolate-spheroidal formulations, $j_z=\Omega$ for the molecular state labels such as $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$0 (Solovyev et al., 2023, Kullie et al., 2022).

A basic structural feature is the four-component spinor decomposition into large and small Pauli-spinor components, $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$1. In block form, one commonly writes

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$2

with $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$3 in the convention used by the minmax finite-element formulations (Kullie, 2024, Kullie et al., 2022).

The spectral difficulty is intrinsic to the Dirac operator: its spectrum is unbounded above and below, with a negative-energy continuum. As a result, naive Rayleigh–Ritz minimization is unstable and can generate variational collapse or spurious states. The literature represented here treats this as a defining numerical issue of the two-center Dirac equation rather than a secondary implementation detail (Kullie, 2024, Fillion-Gourdeau et al., 2015).

Boundary conditions depend on the nuclear model. For bound states one imposes square-integrability and decay at infinity. For extended-charge nuclei, the solutions remain finite at the nuclear centers. For point nuclei, the local behavior is singular but controlled, with near-nuclear asymptotics of the form $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$4, where

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$5

for the corresponding angular channel (Solovyev et al., 2023, Kullie, 2024).

2. Symmetry reduction and coordinate representations

Axial symmetry permits analytic separation of the azimuthal dependence. In spherical coordinates, a standard ansatz for the four-component spinor is

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$6

where $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$7 are large-component functions and $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$8 are small-component functions. Substitution into the stationary Dirac equation yields four coupled first-order PDEs in $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$9 (Solovyev et al., 2023).

A distinct, widely used representation adopts prolate spheroidal coordinates $$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$0 adapted to the two foci. Their standard definition is

$$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$1

with

$$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$2

This coordinate system is central in high-precision finite-element treatments because it follows the two-center geometry and converts the three-dimensional problem, after separation of the azimuthal phase, into a two-dimensional problem in transformed coordinates (Kullie, 2024, Kullie et al., 2022).

A third representation uses Cassini coordinates $$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$3, where

$$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$4

In that representation, the reduced Hamiltonian is written as

$$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$5

with geometry-dependent functions $$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$6 and $$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$7. An important feature of this formulation is that the position-dependent spinor rotation contains a discontinuity across the internuclear line, so $$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$8 develops step-like behavior; this directly motivated mixed bases containing both B-splines and step-like functions (Hahn et al., 2016).

The time-dependent collision literature also employs spherical coordinates with a multipole expansion of the two-center potential,

$$V(\mathbf r)=-\frac{Z_1 e^2}{|\mathbf r-\mathbf R_1|}-\frac{Z_2 e^2}{|\mathbf r-\mathbf R_2|}.$$9

where the monopole term $$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$0 defines a central reference Hamiltonian and the higher multipoles generate channel couplings (McConnell et al., 2012). This makes explicit that “two-center Dirac equation” denotes not a single numerical representation, but a family of equivalent formulations distinguished by symmetry exploitation and basis construction.

3. Variational principles, balance conditions, and discretization strategies

A central numerical issue is the relation between large and small components. In the nonrelativistic limit, kinetic balance imposes

$$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$1

or, in units with $$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$2,

$$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$3

Dual-kinetic balance generalizes this to treat positive- and negative-energy sectors on equal footing and is used to eliminate spurious states in finite-basis Dirac calculations (Solovyev et al., 2023).

One major route is the A-DKB method for axially symmetric systems. There the bispinor components are expanded in a basis of radial B-splines and angular polynomials,

$$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$4

with $$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$5. In the cited calculations, second-order B-splines and Legendre polynomials are used, with typical grids $$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$6 and $$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$7, yielding about 7 significant digits for energies (Solovyev et al., 2023).

A second route is the 2-spinor minmax formulation. Eliminating the small component,

$$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$8

gives the energy-dependent second-order equation

$$\hat H_D=\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta+V_{\mathrm{nucl}}(\mathbf r),\qquad V_{\mathrm{nucl}}(\mathbf r)=-\frac{\alpha Z_1}{|\mathbf r-\mathbf R_1|}-\frac{\alpha Z_2}{|\mathbf r-\mathbf R_2|},$$9

or, in weak form,

$$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$0

This weak form underlies the high-precision minmax-FEM calculations and converges from above to physical electronic eigenvalues while excluding spurious states (Kullie, 2024, Kullie et al., 2022).

The finite-element implementations use high-order complete polynomials on triangular elements in transformed two-dimensional domains. Global prefactors encode both axial-angular behavior and Coulomb cusps,

$$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$1

and singular coordinate transformations concentrate mesh points near the nuclei. In the 2024 minmax benchmark, $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$2 elements are used throughout, with typical mapping exponents $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$3 for $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$4 and $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$5 for $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$6 (Kullie, 2024).

Other discretization strategies remain relevant. The unsplit 3-D Galerkin method in prolate spheroidal coordinates uses atomically or kinetically balanced B-spline bases and solves either the stationary generalized eigenproblem

$$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$7

or the semi-discrete time-dependent system

$$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$8

with PETSc/SLEPc for sparse algebra (Fillion-Gourdeau et al., 2015). An algebraic STSO approach instead builds molecular orbitals from Slater-type spinor orbitals and solves a generalized eigenvalue problem with overlap, kinetic, and nuclear-attraction matrix blocks evaluated in ellipsoidal coordinates (Bagci et al., 2014).

4. Light one-electron quasi-molecules and precision benchmarks

Light one-electron systems provide the cleanest testing ground for the stationary two-center Dirac equation. The literature summarized here considers $$H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+c^2\beta+V$$9–${\rm H}_2^+$0 (${\rm H}_2^+$1), ${\rm H}_2^+$2–${\rm H}_2^+$3 (${\rm H}_2^+$4), and ${\rm H}_2^+$5–${\rm H}_2^+$6 across internuclear distances from atomic scales down to the femtometer regime (Solovyev et al., 2023).

For ${\rm H}_2^+$7, A-DKB yields ground-state energies that reproduce known relativistic corrections relative to high-precision nonrelativistic data. Representative values are

${\rm H}_2^+$8

with the relativistic shift increasing in magnitude as ${\rm H}_2^+$9 decreases. At $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$0 a.u., the A-DKB ground-state energy agrees with fully relativistic calculations to a relative deviation of $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$1 (Solovyev et al., 2023).

The minmax-FEM benchmarks sharpen this substantially. For $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$2 at $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$3 a.u., the reported values are

$$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$4

$$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$5

$$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$6

with estimated fractional uncertainties of $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$7 for $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$8 in the 2024 minmax-FEM study, and a few times $$\mathbf R_1=(0,0,-R/2),\qquad \mathbf R_2=(0,0,+R/2),$$9 for the total energy with a relativistic-shift fractional uncertainty of $z$0 in the 2022 FEM study (Kullie, 2024, Kullie et al., 2022).

The physically expected asymptotic limits are explicit in the data. As $z$1, homonuclear levels tend to the spectrum of a united atom with total nuclear charge $z$2; as $z$3, they approach separated hydrogenic levels. The adiabatic potential curves exhibit gerade/ungerade splitting and relativistic fine-structure resolution of the $z$4 manifold into $z$5, $z$6, and $z$7. For $z$8, an explicit crossing is found around $z$9 fm where the $J_z\psi=m_J\psi,$0 term meets the degenerate $J_z\psi=m_J\psi,$1 pair (Solovyev et al., 2023).

The same framework extends to $J_z\psi=m_J\psi,$2 and $J_z\psi=m_J\psi,$3. For $J_z\psi=m_J\psi,$4, the ground-state energies include

$J_z\psi=m_J\psi,$5

For $J_z\psi=m_J\psi,$6–$J_z\psi=m_J\psi,$7,

$J_z\psi=m_J\psi,$8

The reported interpretation is that the adiabatic curves are intermediate between $J_z\psi=m_J\psi,$9 and $j_z=\Omega$0, consistent with $j_z=\Omega$1 (Solovyev et al., 2023).

5. Heavy quasi-molecules, supercriticality, and QED

For heavy systems, the two-center Dirac equation is used to study quasi-molecular ground states, continuum approach, and the onset of supercritical phenomena. Benchmark calculations have been reported for $j_z=\Omega$2 in a point-nucleus minmax-FEM treatment and for heteronuclear systems such as Bi–Au, U–Pb, and Cf–U in a DKB finite-basis treatment with finite nuclear size (Kullie, 2024, Kotov et al., 2023).

In the heavy homonuclear case $j_z=\Omega$3, the 2024 minmax-FEM work reports at $j_z=\Omega$4 a.u.

$j_z=\Omega$5

$j_z=\Omega$6

$j_z=\Omega$7

with estimated fractional uncertainty of $j_z=\Omega$8 for the relativistic shift (Kullie, 2024). The combined nuclear charge is $j_z=\Omega$9, which exceeds the nominal single-center supercritical threshold $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$00, but at the studied distance the bound-state energy remains well above $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$01, so no diving occurs in the reported data (Kullie, 2024).

The heteronuclear DKB study emphasizes a different issue: the monopole approximation is not unique for heteronuclear systems because the spherically averaged potential depends on the placement of the coordinate-system origin. Three choices are analyzed: the midpoint, the heavy nucleus, and the light nucleus. For U–Pb at $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$02 fm, the $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$03 energies are

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$04

and for Cf–U at $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$05 fm,

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$06

This documents multi-keV origin dependence in the monopole approximation and establishes that the full two-center calculation is required for a stable non-QED baseline in heteronuclear systems (Kotov et al., 2023).

The same heteronuclear study evaluates leading one-electron QED corrections within the midpoint monopole approximation. The vacuum-polarization contribution is represented by the Uehling potential, while the self-energy is computed through the renormalized bound-state self-energy operator. Representative totals are

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$07

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$08

The reported comparison is that these QED shifts are typically an order of magnitude smaller than the TC–MA(1) discrepancy of the underlying Dirac energy (Kotov et al., 2023).

A recurring theme is the distinction between supercriticality and merely large relativistic shifts. The light-ion study states explicitly that no supercritical behavior occurs for $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$09, $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$10, or $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$11; that phenomenon requires much larger effective charges (Solovyev et al., 2023). The heavy-ion literature, by contrast, treats near-continuum approach and eventual resonance formation as a principal physical motivation, while noting that resonance widths would require techniques such as complex scaling or exterior complex scaling beyond the present bound-state setups (Kullie, 2024).

6. Time dependence, coupled channels, and reduced models

The time-dependent two-center Dirac equation is used when the nuclei move, as in slow ion–ion collisions or $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$12 decay. In one line of work, the full time-dependent equation

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$13

is solved in prolate spheroidal coordinates by an unsplit 4-component Galerkin method with atomically or kinetically balanced B-spline bases. Time propagation is performed through a norm-conserving Crank–Nicolson scheme,

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$14

and the framework has been applied to spectral computation and to driven evolution in external electromagnetic fields (Fillion-Gourdeau et al., 2015).

A complementary approach expands the stationary two-center Hamiltonian in a monopole basis and retains the full multipole expansion of the electron–nuclei interaction in spherical coordinates. There the exact stationary two-center eigenfunctions are expanded in eigenfunctions of the monopole Hamiltonian, producing a generalized eigenvalue problem for each internuclear separation. The time-dependent amplitudes then satisfy coupled-channel equations containing radial nonadiabatic couplings and, for nonzero impact parameter, rotational couplings through $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$15 (McConnell et al., 2012).

This multipole coupled-channel method has been used to calculate K- and L-shell ionization probabilities in $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$16 decay and in slow $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$17–$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$18 collisions. For example, in $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$19 decay of hydrogen-like $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$20, the asymptotic K-shell ionization probability is reported as $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$21 in the full multipole treatment versus $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$22 in the monopole-only approximation, while the $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$23 L-shell probability is $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$24 in the full calculation and $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$25 in the monopole-only approximation (McConnell et al., 2012). This establishes that higher multipoles are quantitatively important at large $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$26 even when the monopole basis is a useful starting point.

A reduced but conceptually illuminating variant is the one-dimensional Dirac equation with two $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$27-function centers. For symmetric configurations with centers at $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$28, one can derive closed transcendental equations for the bound-state energies either from a Green’s-function determinant or from transfer matrices. For the trigonometric self-adjoint extension, the symmetric double-well spectrum obeys

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$29

while the dipole-like antisymmetric case yields

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$30

The reduced model makes explicit that the spectral problem depends on the chosen self-adjoint extension: different connection matrices, both satisfying the self-adjointness condition, lead to different merging-limit behavior, including whether strength additivity holds when the two centers coalesce (Gusynin et al., 2023).

In the full three-dimensional Coulomb problem, this does not imply analogous ambiguity in the physical Dirac operator itself; rather, it illustrates a broader principle that relativistic two-center problems are sensitive to boundary conditions, cusp structure, and the correct treatment of singular interactions. A plausible implication is that many numerical controversies in the three-dimensional literature—spurious states, origin dependence of monopole approximations, and slow convergence near cusps—are different manifestations of this same structural sensitivity.

7. Numerical pathologies, approximations, and methodological significance

Three recurring numerical issues dominate the modern literature. The first is spurious-state contamination. DKB, atomic balance, and minmax formulations all address this directly. The A-DKB literature attributes spurious-state removal to the balanced construction of the four-component basis (Solovyev et al., 2023). The atomically balanced Galerkin method states that naive Rayleigh–Ritz without balance can produce nonconvergent levels, whereas atomic balance eliminates spectral pollution for Coulomb potentials under the cited conditions (Fillion-Gourdeau et al., 2015). The minmax literature states that convergence is from above to physical electronic eigenvalues and that the negative-energy continuum is projected out (Kullie, 2024, Kullie et al., 2022).

The second issue is Coulomb-singularity resolution. The best-performing methods build cusp information directly into the discretization through global factors such as

$$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$31

through singular coordinate mappings, or through repeated knots and radial clustering near the nuclei (Kullie, 2024, Kullie et al., 2022, Fillion-Gourdeau et al., 2015). The Cassini-coordinate work shows that even geometric spinor discontinuities induced by coordinate representation can materially affect convergence and may require non-smooth basis enrichment (Hahn et al., 2016).

The third issue is approximation hierarchy. Monopole approximations are useful, but the collision and heteronuclear studies show their limitations in different ways. In the time-dependent collision problem, the monopole approximation underestimates ionization probabilities, especially at large internuclear distance (McConnell et al., 2012). In heteronuclear stationary problems, the monopole approximation becomes origin-dependent and can differ from the full two-center result by tens of keV (Kotov et al., 2023). These results make clear that the monopole approximation is a computational device rather than an invariant physical reduction.

Across the reported calculations, the two-center Dirac equation emerges as a benchmark problem for relativistic numerical analysis. It has supported energies with 20+ significant digits for $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$32 (Kullie et al., 2022), benchmark fractional uncertainties of $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$33 for light systems and $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$34 for heavy ones in the minmax-FEM framework (Kullie, 2024), and fully relativistic adiabatic spectra for light quasi-molecules over broad $$\hat H_D \psi(\mathbf r)=E\psi(\mathbf r),\qquad \hat H_D=c\,\boldsymbol{\alpha}\!\cdot\!\mathbf p+\beta m c^2+V(\mathbf r),$$35 ranges in A-DKB calculations (Solovyev et al., 2023). At the same time, the literature emphasizes that extensions to continuum resonances, rigorous two-center QED, finite nuclear size in high-precision heavy-ion work, and multi-electron generalizations remain active directions rather than closed topics (Kullie, 2024, Kotov et al., 2023).