Relativistic Two-Center Methods

• Relativistic two-center approaches are formulations that solve the Dirac equation in systems with two fixed Coulomb centers, capturing key molecular and collision dynamics.
• They utilize diverse numerical methods such as DKB finite-basis expansions, finite-element techniques, and cubic Hermite splines to ensure stability and prevent variational collapse.
• Time-dependent implementations and fully relativistic two-body formulations extend the approach to heavy-ion collisions, QED corrections, and accurate center-of-mass reductions.

In the literature considered here, the relativistic two-center approach denotes a class of fully relativistic formulations for systems governed by two nuclei, two fixed Coulomb centers, or two interacting constituents. In atomic and molecular applications, it is centered on the stationary or time-dependent Dirac equation in a two-center potential for quasimolecules, heavy-ion collisions, and molecular ions such as ${\rm H}_2^+$ and its isotopologues; in complementary two-body mechanics, it concerns covariant center-of-mass and relative-motion reductions, interaction-dependent boosts, and the fate of integrability under relativity (Kotov et al., 2023, McConnell et al., 2012, Kullie et al., 2022, Shields et al., 2013).

1. Conceptual definition and mathematical setting

The canonical one-electron formulation is the relativistic two-center Dirac equation

$$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$

with a two-center nuclear potential of the form

$$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$

For fixed internuclear distance in the Born-Oppenheimer approximation, this problem is axially symmetric around the internuclear axis. For general geometries, the potential is often expanded as

$V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$

so that the monopole term $V_0(r)$ defines a spherically symmetric auxiliary problem and the higher multipoles reconstruct the full two-center interaction (Kotov et al., 2023, McConnell et al., 2012, Solovyev et al., 2023).

This framework is used for heteronuclear quasimolecules such as Bi–Au, U–Pb, and Cf–U, for light one-electron quasi-molecular ions H-H$^+$, He$^+$-He$^{2+}$, and He$^+$-H$^+$, and for the two-center Coulomb problem relevant to $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$0 and isotopologues (Kotov et al., 2023, Solovyev et al., 2023, Kullie et al., 2022). In time-dependent settings, the same structural decomposition appears in electron dynamics in the field of moving ions, in slow ion-ion collisions, and in ionization during nuclear $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$1-decay (McConnell et al., 2012, Deyneka et al., 2013).

A recurring technical theme is that the full two-center problem is not spherically symmetric, whereas the monopole approximation is. The monopole problem therefore serves both as a computational approximation and as a basis-generation device for solving the full axially symmetric or fully time-dependent problem (McConnell et al., 2012, Kotov et al., 2023).

2. Stationary formulations and basis constructions

The stationary relativistic two-center problem has been attacked with several numerically distinct, but structurally related, basis constructions. A common requirement is control of the Dirac spectrum so that spurious states or variational collapse are avoided while retaining accurate access to bound and continuum sectors.

The dual-kinetically balanced finite-basis-set expansion is a central example. In the heteronuclear Dirac problem it is used in a rigorous two-center treatment and in monopole calculations, with the DKB constraint ensuring a proper relation between large and small components and suppressing spurious states. In the axially symmetric implementation termed A-DKB, the bispinor is expanded in B-splines in $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$2 and Legendre polynomials in $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$3, with the DKB constraint imposed through a $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$4 matrix. The reported advantages include applicability for arbitrary internuclear distance, direct solution of the Dirac equation, numerical stability, and the absence of partial-wave decomposition; the reported disadvantages include high computational cost and a rapidly increasing basis size when higher accuracy is sought (Kotov et al., 2023, Solovyev et al., 2023).

A different algebraic route uses Slater-type spinor orbitals. In that approach, STSOs are employed in an LCAO-MO expansion, the generalized matrix equation

$$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$5

is solved, and one- and two-center integrals are evaluated in ellipsoidal coordinates by a numerical global-adaptive method with Gauss-Kronrod extension. The screening constants are treated as variational parameters for each nuclear separation, and the construction is designed to enforce kinetic balance and to avoid variational collapse (Bagci et al., 2014).

For ultra-high precision in the Coulomb two-center problem, the finite-element method has been combined with the 2-spinor minmax principle. There the Dirac equation is recast into a nonlinear 2-spinor problem, solved iteratively after elimination of the small component. Prolate spheroidal coordinates and additional singular coordinate transformations are used to regularize the nuclear singularities, triangular elements discretize the domain, and quadruple precision arithmetic is employed (Kullie et al., 2022).

A separate time-dependent lattice-based strategy uses cubic Hermite splines on a three-dimensional grid. The local support of the splines yields sparse overlap and Hamiltonian matrices, and the Crank-Nicolson propagator preserves the norm during time evolution. That method has been used in monopole, axially symmetric, and full 3D variants for U$$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$6–U$$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$7 collisions (Deyneka et al., 2013).

Framework	Problem class	Distinctive feature
DKB finite-basis-set expansion (Kotov et al., 2023)	Heteronuclear quasimolecules	Rigorous TC and monopole calculations
A-DKB (Solovyev et al., 2023)	Light one-electron quasi-molecular ions	B-splines and Legendre polynomials with DKB
Monopole-basis coupled-channel (McConnell et al., 2012)	Time-dependent moving-ion fields	Full multipole expansion on monopole eigenfunctions
Cubic Hermite splines on a lattice (Deyneka et al., 2013)	Time-dependent heavy-ion collisions	Sparse 1D, 2D, and 3D propagators
2-spinor minmax FEM (Kullie et al., 2022)	Coulomb two-center precision spectroscopy	Variationally stable high-precision Dirac solver
STSO algebraic method (Bagci et al., 2014)	One-electron atoms and diatomics	Ellipsoidal-coordinate integral evaluation

3. Monopole reduction, coordinate origin, and precision hierarchy

The monopole approximation is not a unique object in heteronuclear systems. In the Bi–Au, U–Pb, and Cf–U calculations, three coordinate-system origins were considered: halfway between the nuclei, at the heavy nucleus, and at the light nucleus. A substantial difference between the resulting monopole energies was found. The midpoint choice, MA(1), gives results closest to the full two-center energies, but a steadily growing deviation appears as the internuclear distance grows; MA(2) and MA(3) yield progressively larger deviations and differ substantially from one another; all three monopole schemes agree with the full two-center result as $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$8 (Kotov et al., 2023).

For U–Pb, two dominant uncertainty sources were identified: basis-set error and nuclear model error. The basis-set error was found to be minor, with the quoted value $$\Big[ \vec{\alpha} \cdot \vec{p} + \beta + V(\vec{r}) \Big] \Psi_n(\vec{r}) = E_n \Psi_n(\vec{r}),$$9, whereas uncertainty in the RMS nuclear radii and nuclear modeling dominates at small internuclear distance (Kotov et al., 2023). The same work evaluates the leading one-electron QED effects—self-energy and vacuum polarization—within MA(1) only, and reports that these QED corrections are generally an order of magnitude smaller than the difference between the rigorous TC and MA(1) Dirac energies (Kotov et al., 2023).

At the high-precision end of the two-center Coulomb problem, the finite-element minmax approach yields total energies with estimated fractional uncertainties of a few times $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$0 for unit charges and bond length $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$1 atomic units, while the purely relativistic contribution carries a fractional uncertainty of $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$2. The result is described as relevant for future precision experiments, although at present the dominant uncertainties in rovibrational transition frequencies arise from the QED treatment (Kullie et al., 2022).

QED precision requirements also motivate the calculation of the relativistic Bethe logarithm in the Coulomb two-center problem. A variational approach has achieved 10 significant digits for the nonrelativistic Bethe logarithm and 3–4 significant digits for the relativistic corrections starting from $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$3 bohr. The stated implication is that transition frequencies in systems such as $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$4, $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$5, and antiprotonic helium can reach $$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$6 relative uncertainty once these terms are incorporated (Korobov et al., 2013).

These results establish a clear precision hierarchy: the full two-center geometry can dominate over monopole-based simplifications, while high-order radiative terms become decisive only after the underlying Dirac problem is controlled to comparable accuracy.

4. Time-dependent dynamics in strong fields and collisions

For moving nuclei, the relativistic two-center approach becomes explicitly time dependent. One formulation writes

$$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$7

with

$$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$8

In the spherical-coordinate method, eigenfunctions of the full two-center Hamiltonian are expanded in a quasi-complete set of monopole eigenfunctions,

$$V(\vec{r}) = V^A_{\text{nucl}}(|\vec{r}-\vec{R}_A|) + V^B_{\text{nucl}}(|\vec{r}-\vec{R}_B|).$$9

and the time-dependent wavefunction is propagated in a coupled-channel expansion,

$V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$0

The nuclei move on Rutherford trajectories, and the channel amplitudes satisfy a coupled system driven by the time dependence of the two-center basis states (McConnell et al., 2012).

This non-perturbative scheme was applied to K- and L-shell ionization of hydrogen-like ions during nuclear $V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$1-decay and to slow ion-ion collisions. The full multipole expansion was found to be essential, especially at large internuclear distance. In the $V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$2-decay application, the monopole approximation underestimates ionization probabilities by up to 30% for large $V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$3, while the non-perturbative results reproduce first-order perturbative predictions within about 5% for K and L shells at large $V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$4 (McConnell et al., 2012).

The cubic-Hermite-spline lattice approach provides a second time-dependent implementation. There the wavefunction is expanded in a time-independent basis on a 3D grid, the matrix equation

$V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$5

is propagated with the Crank-Nicolson scheme, and separate 1D, 2D, and 3D realizations correspond respectively to the monopole, axially symmetric, and full three-dimensional problems (Deyneka et al., 2013). In U$V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$6–U$V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$7 collisions at $V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$8MeV/u, the 2D and 3D charge-transfer probabilities were reported to agree closely with one another and with previous Dirac-Sturm calculations. The same study found the negative-energy continuum contribution to be small and interpreted the high ground-state survival probabilities as evidence of adiabatic dynamics near the Coulomb barrier (Deyneka et al., 2013).

A central methodological conclusion is that the monopole problem can serve as a useful starting point for time propagation, but quantitative observables in dynamical strong-field regimes depend on higher multipoles, explicit geometry, and, in charge-transfer problems, full multi-center structure.

5. Relativistic two-body formulations, boosts, and clock behavior

Outside the electronic Dirac problem, the two-center idea also appears in fully relativistic two-body mechanics. One classical construction uses the Bakamjian-Thomas framework, where the generators $V(r,\theta)=\sum_l V_l(r) P_l(\cos\theta),$9, $V_0(r)$0, and $V_0(r)$1 satisfy the Poincaré Lie algebra and the interacting Hamiltonian is written as

$V_0(r)$2

For a regularized Coulomb potential $V_0(r)$3, numerical simulations compared two dynamics: one with the center of mass at rest and one obtained by actively boosting the initial state with the full, interaction-dependent boost generator. The reported result is that standard kinematic Lorentz formulas fail to transform particle observables between the two frames. The explanation is the inevitable interaction dependence of the boost generator in the instant form of relativistic dynamics, presented as a direct illustration of the Currie-Jordan-Sudarshan No-Interaction Theorem. At the same time, the orbital periods are related by Einstein’s time-dilation factor,

$V_0(r)$4

for all interaction strengths (Shields et al., 2013).

In predictive relativistic mechanics, an isolated Hamiltonian two-body system is formulated on a covariant phase space with an equal-time condition. The natural center is a center of energy rather than a Newtonian center of mass. At equal times in the rest frame it reduces to

$V_0(r)$5

with $V_0(r)$6 the individual energies in the center-of-mass frame. Relative motion has a structure similar to that of a nonrelativistic one-body motion in a stationary external potential, but the evolution parameter is generally not a linear function of the center-of-mass time unless the relative motion is circular; in the circular case, the motion is periodic in the center-of-mass time (Droz-Vincent, 2010).

Several covariant wave-equation programs pursue analogous reductions. One derives exact covariant operator equations in center-mass and relative variables beyond the Breit frame and identifies a class of external vector fields, linear in the coordinates, for which the problem can be reformulated in the $V_0(r)$7 phase space of relative motion variables uncoupled from the center-of-mass motion. For fermions, a new $V_0(r)$8-matrix basis is introduced so that the 16-component Dirac-like equation decouples into four independent equations for four-component spinors (Koshelkin, 7 Jan 2025). A related extension of the Feshbach-Villars formalism to two spin-0 particles separates the center-of-mass motion and yields a relative-coordinate equation of Feshbach-Villars type with the correct total-mass rest term $V_0(r)$9 (Papp, 4 May 2026). A broader review of two-body relativistic wave equations presents covariant kinematics, reduced phase space, and explicit scalar-scalar, fermion-scalar, and fermion-fermion equations with central scalar and vector potentials, recovering the one-particle and nonrelativistic limits (Giachetti et al., 2018).

Taken together, these results show that relativistic two-center reasoning is not restricted to fixed nuclei. It also encompasses the problem of how a composite relativistic system is boosted, separated into intrinsic and collective motion, and represented in a covariant Hamiltonian or wave-equation form.

6. Symmetry, non-integrability, and one-center versus two-center approximations

The relativistic two-center problem does not generally inherit the hidden symmetries of its Newtonian analogue. In Newtonian gravity, the Euler two-center problem admits a third Carter-like constant of motion when the mass multipole moments satisfy the Kerr-like pattern $^+$0, $^+$1. Its relativistic analogue, the Bach-Weyl solution for two point masses held fixed apart, does not admit such a constant: no non-trivial second-rank Killing tensor exists, and the existence of a Carter-like constant already fails at first post-Newtonian order (Mirshekari et al., 2010). This places a structural limit on attempts to transfer Newtonian two-center integrability into general relativity.

A different, but related, limitation appears in relativistic electronic-structure approximations. Within exact two-component theory, one-center approximations to relativistic two-electron picture-change corrections can be highly effective for localized spin-dependent effects, but not for more delocalized scalar-relativistic effects. The X2CAMF scheme, which keeps one-center spin-dependent Coulomb and Breit contributions, is reported to recover nearly all relevant two-electron spin-orbit effects for standard basis sets, whereas one-center scalar-relativistic approximations can become unstable or inaccurate, especially with diffuse basis sets. The recommendation given there is that scalar-relativistic two-electron effects require at least two-center treatments when high accuracy is needed, and that MP/AMF corrections must be treated as one-electron operators in the energy expression (Wang et al., 28 Apr 2025).

This suggests that “two-center” is not merely a descriptor of geometry. It also marks an approximation threshold. In some problems, one-center reductions are numerically efficient and physically adequate; in others, especially when nonlocal couplings, large internuclear separations, or scalar-relativistic corrections dominate, the full two-center structure is indispensable.

The relativistic two-center approach is therefore best understood as a family of methods tied together by a common demand: relativity must be enforced simultaneously with the presence of two distinguished centers or constituents. The resulting formulations differ in representation—Dirac, Feshbach-Villars, Hamiltonian, coupled-channel, finite-element, spline, or algebraic basis—but they converge on the same technical lesson: boosts, basis balance, multipole content, center-of-mass reduction, and approximation level are all system dependent, and the accuracy of any reduced description must be validated against the full relativistic two-center problem.