---
title: Two-Center Dirac Equation
url: https://www.emergentmind.com/topics/two-center-dirac-equation
type: topic
---

# Two-Center Dirac Equation

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$ [2310.04057][2411.12427]. 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^+$ [2204.07087][1208.4731].

## 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 $1\sigma_g$ [2310.04057][2204.07087].

A basic structural feature is the four-component spinor decomposition into large and small Pauli-spinor components, $\psi=(\phi_+,\phi_-)^T$. In block form, one commonly writes
\[
\begin{pmatrix}
V & L\\
L & V-2c^2
\end{pmatrix}
\begin{pmatrix}
\phi_+\\
\phi_-
\end{pmatrix}
=
\epsilon
\begin{pmatrix}
\phi_+\\
\phi_-
\end{pmatrix},
\qquad
L=-ic\,\boldsymbol{\sigma}\!\cdot\!\nabla,
\]
with $\epsilon=E-mc^2$ in the convention used by the minmax finite-element formulations [2411.12427][2204.07087].

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 [2411.12427][1507.07398].

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 $r_l^{-1+\gamma_{l,\kappa}}$, where
\[
\gamma_{l,\kappa}=\sqrt{\kappa^2-(\alpha Z_l)^2}
\]
for the corresponding angular channel [2310.04057][2411.12427].

## 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
\[
\psi(r,\theta,\varphi)=\frac1r
\begin{pmatrix}
G_1(r,\theta)e^{i(m_J-\frac12)\varphi}\\
G_2(r,\theta)e^{i(m_J+\frac12)\varphi}\\
iF_1(r,\theta)e^{i(m_J-\frac12)\varphi}\\
iF_2(r,\theta)e^{i(m_J+\frac12)\varphi}
\end{pmatrix},
\]
where $G_{1,2}$ are large-component functions and $F_{1,2}$ are small-component functions. Substitution into the stationary Dirac equation yields four coupled first-order PDEs in $(r,\theta)$ [2310.04057].

A distinct, widely used representation adopts prolate spheroidal coordinates $(\xi,\eta,\phi)$ adapted to the two foci. Their standard definition is
\[
x=\frac R2\,u(\xi,\eta)\cos\phi,\qquad
y=\frac R2\,u(\xi,\eta)\sin\phi,\qquad
z=\frac R2\,\xi\eta,
\]
with
\[
u(\xi,\eta)=\sqrt{(\xi^2-1)(1-\eta^2)},\qquad
r_1=\frac R2(\xi+\eta),\qquad
r_2=\frac R2(\xi-\eta).
\]
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 [2411.12427][2204.07087].

A third representation uses Cassini coordinates $(w,\delta,\varphi)$, where
\[
w=\frac{\sqrt{r_1r_2}}{a},\qquad \delta=\frac{\theta_1+\theta_2}{2},\qquad a=\frac12|\mathbf R_1-\mathbf R_2|.
\]
In that representation, the reduced Hamiltonian is written as
\[
\mathcal H_{\mathrm{Cassini}}^{(\mu)}
=
-i\frac{D^{1/4}}{aw}
\left(
\alpha_3\partial_w+\alpha_1\frac1w\partial_\delta
\right)
+\alpha_2\frac{\mu}{\rho}
+V(w,\delta;a)+\beta,
\]
with geometry-dependent functions $D(w,\delta)$ and $\rho(w,\delta)$. An important feature of this formulation is that the position-dependent spinor rotation contains a discontinuity across the internuclear line, so $\Psi/\sqrt{\rho}$ develops step-like behavior; this directly motivated mixed bases containing both B-splines and step-like functions [1610.05263].

The time-dependent collision literature also employs spherical coordinates with a multipole expansion of the two-center potential,
\[
V_{\mathrm{TC}}(r,\theta;R)=\sum_{\ell=0}^\infty V_\ell(r,R)P_\ell(\cos\theta),
\]
where the monopole term $V_0$ defines a central reference Hamiltonian and the higher multipoles generate channel couplings [1208.4731]. 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
\[
\chi \approx \frac{\boldsymbol{\sigma}\!\cdot\!\mathbf p}{2mc}\,\phi,
\]
or, in units with $\hbar=c=m_e=1$,
\[
\chi \approx \frac{\boldsymbol{\sigma}\!\cdot\!\mathbf p}{2}\,\phi.
\]
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 [2310.04057].

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,
\[
\phi(r,\theta)\approx
\sum_{u=1}^{4}\sum_{i_r=1}^{N_r}\sum_{i_\theta=1}^{N_\theta}
\Lambda\,C^u_{i_ri_\theta}\,B_{i_r}(r)\,Q_{i_\theta}(\theta)\,e_u,
\]
with $Q_{i_\theta}(\theta)=P_{i_\theta-1}(2\theta/\pi-1)$. In the cited calculations, second-order B-splines and Legendre polynomials are used, with typical grids $N_r=550$ and $N_\theta=80$, yielding about 7 significant digits for energies [2310.04057].

A second route is the 2-spinor minmax formulation. Eliminating the small component,
\[
\phi_-=(\epsilon+2c^2-V)^{-1}L\phi_+,
\]
gives the energy-dependent second-order equation
\[
L\!\left[(\epsilon+2c^2-V)^{-1}L\phi_+\right]
=
(\epsilon-V)\phi_+,
\]
or, in weak form,
\[
\int \frac{c^{-2}|L\phi_+|^2}{2+(\epsilon-V)/c^2}\,d^3r
=
\int (\epsilon-V)|\phi_+|^2\,d^3r.
\]
This weak form underlies the high-precision minmax-FEM calculations and converges from above to physical electronic eigenvalues while excluding spurious states [2411.12427][2204.07087].

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,
\[
G_1^k(s,t)=[(\xi^2-1)(1-\eta^2)]^{m_k/2},\qquad
G_2(s,t)=r_1^{-1+\gamma_{1,\kappa}}r_2^{-1+\gamma_{2,\kappa}},
\]
and singular coordinate transformations concentrate mesh points near the nuclei. In the 2024 minmax benchmark, $p=10$ elements are used throughout, with typical mapping exponents $\nu=8$ for ${\rm H}_2^+$ and $\nu=10$ for ${\rm Th}_2^{179+}$ [2411.12427].

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
\[
C\,\mathbf c = E\,S\,\mathbf c
\]
or the semi-discrete time-dependent system
\[
i\,S\,\dot{\mathbf c}(t)=\big(C+D(t)\big)\mathbf c(t),
\]
with PETSc/SLEPc for sparse algebra [1507.07398]. 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 [1407.6720].

## 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 ${\rm H}$–${\rm H}^+$ (${\rm H}_2^+$), ${\rm He}^+$–${\rm He}^{2+}$ (${\rm He}_2^{3+}$), and ${\rm He}^+$–${\rm H}^+$ across internuclear distances from atomic scales down to the femtometer regime [2310.04057].

For ${\rm H}_2^+$, A-DKB yields ground-state energies that reproduce known relativistic corrections relative to high-precision nonrelativistic data. Representative values are
\[
R=0.1:\quad E_{\mathrm{A-DKB}}=-1.9783363,\qquad
R=2.0:\quad E_{\mathrm{A-DKB}}=-1.10264160,
\]
with the relativistic shift increasing in magnitude as $R$ decreases. At $R=2$ a.u., the A-DKB ground-state energy agrees with fully relativistic calculations to a relative deviation of $2.08\times10^{-8}$ [2310.04057].

The minmax-FEM benchmarks sharpen this substantially. For ${\rm H}_2^+$ at $R=2$ a.u., the reported values are
\[
E_{\mathrm{rel}}=-1.10264158103257716411812499995\ \mathrm{au},
\]
\[
E_{\mathrm{nrel}}=-1.10263421449494646150896894154\ \mathrm{au},
\]
\[
\Delta E=-7.3665376307026091560584\times10^{-6}\ \mathrm{au},
\]
with estimated fractional uncertainties of $\sim10^{-23}$ for ${\rm H}_2^+$ in the 2024 minmax-FEM study, and a few times $10^{-20}$ for the total energy with a relativistic-shift fractional uncertainty of $1\times10^{-17}$ in the 2022 FEM study [2411.12427][2204.07087].

The physically expected asymptotic limits are explicit in the data. As $R\to0$, homonuclear levels tend to the spectrum of a united atom with total nuclear charge $Z_1+Z_2$; as $R\to\infty$, they approach separated hydrogenic levels. The adiabatic potential curves exhibit gerade/ungerade splitting and relativistic fine-structure resolution of the $n=2$ manifold into $2s_{1/2}$, $2p_{1/2}$, and $2p_{3/2}$. For ${\rm H}_2^+$, an explicit crossing is found around $R\approx500$ fm where the $2p_{3/2}$ term meets the degenerate $2s_{1/2}/2p_{1/2}$ pair [2310.04057].

The same framework extends to ${\rm He}_2^{3+}$ and ${\rm HeH}^{2+}$. For ${\rm He}_2^{3+}$, the ground-state energies include
\[
R=0.1:\ E=-7.715707556,\qquad
R=2.0:\ E=-3.184424360,\qquad
R=10.0:\ E=-2.200146300.
\]
For ${\rm He}^+$–${\rm H}^+$,
\[
R=0.1:\ E=-4.411760188,\qquad
R=2.0:\ E=-2.512296432,\qquad
R=10.0:\ E=-2.100100978.
\]
The reported interpretation is that the adiabatic curves are intermediate between ${\rm H}_2^+$ and ${\rm He}_2^{3+}$, consistent with $Z_1+Z_2=3$ [2310.04057].

## 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 ${\rm Th}_2^{179+}$ 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 [2411.12427][2305.04233].

In the heavy homonuclear case ${\rm Th}_2^{179+}$, the 2024 minmax-FEM work reports at $R=2/90$ a.u.
\[
E_{\mathrm{rel}}=-9504.756648434009500737\ \mathrm{au},
\]
\[
E_{\mathrm{nrel}}=-8931.3371374090663382226\ \mathrm{au},
\]
\[
\Delta E=-573.4195110249431625138\ \mathrm{au},
\]
with estimated fractional uncertainty of $\sim10^{-21}$ for the relativistic shift [2411.12427]. The combined nuclear charge is $180$, which exceeds the nominal single-center supercritical threshold $Z_{\mathrm{crit}}\approx173$, but at the studied distance the bound-state energy remains well above $-mc^2$, so no diving occurs in the reported data [2411.12427].

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 $D=25$ fm, the $1\sigma$ energies are
\[
\mathrm{TC}=-399{,}528\ \mathrm{eV},\qquad
\mathrm{MA(1)}=-386{,}390\ \mathrm{eV},\qquad
\mathrm{MA(2)}=-325{,}462\ \mathrm{eV},\qquad
\mathrm{MA(3)}=-307{,}855\ \mathrm{eV},
\]
and for Cf–U at $D=50$ fm,
\[
\mathrm{TC}=-491{,}640\ \mathrm{eV},\qquad
\mathrm{MA(1)}=-459{,}027\ \mathrm{eV},\qquad
\mathrm{MA(2)}=-383{,}320\ \mathrm{eV},\qquad
\mathrm{MA(3)}=-366{,}680\ \mathrm{eV}.
\]
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 [2305.04233].

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
\[
\Delta E_{\mathrm{QED}}=+1148\ \mathrm{eV}\quad \text{for U--Pb at }D=25\ \mathrm{fm},
\]
\[
\Delta E_{\mathrm{QED}}=+1814\ \mathrm{eV}\quad \text{for Cf--U at }D=50\ \mathrm{fm}.
\]
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 [2305.04233].

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 ${\rm H}_2^+$, ${\rm He}_2^{3+}$, or ${\rm HeH}^{2+}$; that phenomenon requires much larger effective charges [2310.04057]. 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 [2411.12427].

## 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 $\alpha$ decay. In one line of work, the full time-dependent equation
\[
i\,\partial_t\Psi(\mathbf r,t)=H(t)\Psi(\mathbf r,t),
\qquad
H(t)=c\,\boldsymbol{\alpha}\!\cdot\!\big(\mathbf p-e\mathbf A(\mathbf r,t)\big)+\beta mc^2+V(\mathbf r,t),
\]
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,
\[
S\,\mathbf c^{n+1}
=
S\,\mathbf c^n
-\frac{i\Delta t}{2}
\Big[
\big(C+D^n\big)\mathbf c^n+\big(C+D^{n+1}\big)\mathbf c^{n+1}
\Big],
\]
and the framework has been applied to spectral computation and to driven evolution in external electromagnetic fields [1507.07398].

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 $j_y$ [1208.4731].

This multipole coupled-channel method has been used to calculate K- and L-shell ionization probabilities in $\alpha$ decay and in slow ${\rm U}^{91+}$–${\rm U}^{92+}$ collisions. For example, in $\alpha$ decay of hydrogen-like $^{210}{\rm Po}^{+83}$, the asymptotic K-shell ionization probability is reported as $2.1\times10^{-6}$ in the full multipole treatment versus $1.4\times10^{-6}$ in the monopole-only approximation, while the $2p_{3/2}$ L-shell probability is $0.61\times10^{-5}$ in the full calculation and $0.04\times10^{-5}$ in the monopole-only approximation [1208.4731]. This establishes that higher multipoles are quantitatively important at large $R$ even when the monopole basis is a useful starting point.

A reduced but conceptually illuminating variant is the one-dimensional Dirac equation with two $\delta$-function centers. For symmetric configurations with centers at $\pm R$, 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
\[
\frac{\sqrt{m^2-E^2}}{E}
=
\tan g
\mp
\frac{m}{E}e^{-2R\sqrt{m^2-E^2}},
\]
while the dipole-like antisymmetric case yields
\[
E=\pm m\sqrt{\cos^2 g+e^{-4z}\sin^2 g},
\qquad z=R\sqrt{m^2-E^2}.
\]
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 [2309.00856].

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 [2310.04057]. 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 [1507.07398]. The minmax literature states that convergence is from above to physical electronic eigenvalues and that the negative-energy continuum is projected out [2411.12427][2204.07087].

The second issue is Coulomb-singularity resolution. The best-performing methods build cusp information directly into the discretization through global factors such as
\[
r_1^{-1+\gamma_{1,\kappa}}r_2^{-1+\gamma_{2,\kappa}},
\]
through singular coordinate mappings, or through repeated knots and radial clustering near the nuclei [2411.12427][2204.07087][1507.07398]. 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 [1610.05263].

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 [1208.4731]. In heteronuclear stationary problems, the monopole approximation becomes origin-dependent and can differ from the full two-center result by tens of keV [2305.04233]. 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 ${\rm H}_2^+$ [2204.07087], benchmark fractional uncertainties of $\sim10^{-23}$ for light systems and $\sim10^{-21}$ for heavy ones in the minmax-FEM framework [2411.12427], and fully relativistic adiabatic spectra for light quasi-molecules over broad $R$ ranges in A-DKB calculations [2310.04057]. 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 [2411.12427][2305.04233].

Source: https://www.emergentmind.com/topics/two-center-dirac-equation