Two-center Dirac equation in a Gaussian basis set

A method for solving the two-center Dirac equation using a Gaussian-product basis set tailored to the axial symmetry of diatomic systems is presented. For the ground state of H${}_2^{+}$ at an internuclear distance $R=2.0a_0$, we obtain a relativistic energy with relative precision at the ${10}^{-18}$ level. The properties of the basis set allow for an analytic Fourier transform of the wave function. As a first application we evaluate, for the first time, the zero-potential contribution to the one-loop self-energy for H${}_2^{+}$ in Feynman and Coulomb gauges. Although this contribution is not separately gauge invariant, its calculation provides a proof of principle for the use of the present framework in prospective QED calculations.