Pre-Born-Oppenheimer Dirac-Coulomb-Breit computations for two-body systems

The sixteen-component, no-pair Dirac–Coulomb–Breit equation, derived from the Bethe–Salpeter equation, is solved in a variational procedure using Gaussian-type basis functions for the example of positronium, muonium, hydrogen atom, and muonic hydrogen. The $\alpha$ fine-structure-constant dependence of the variational energies, through fitting a function of ${\alpha }^n$ and ${\alpha }^n\text{ln}\alpha$ terms, shows excellent agreement with the relevant energy expressions of the (perturbative) non-relativistic QED framework, and thereby, establishes a solid reference for the development of a computational relativistic QED approach.