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.