\neq Z$ significantly increases the accuracy of the zeroth-order approximation, while rendering the complexity of all calculations low.

However, the fully relativistic description becomes much more complicated due to the structure of Dirac equation. First, the Dirac equation possesses negative energy states and, therefore, there exist matrix elements between electronic states and negative energy ones. This significantly increases the number of required matrix elements to be taken into account. Second, due to the mass term of the Dirac equation the energy of the system does not have a universal behavior of the quadratic function of the effective charge as in the nonrelativistic case. Therefore, it is impossible to separate explicitly the dependence of the matrix elements on the effective charge. For this reason, before investigating the second-order correction to the energy of the system we focused on the zeroth-order approximation first.

At the same time, we point out that since the Dirac-Coulomb Green’s function is known in analytical closed-form [@wong_dirac_1985; @swainson_unified_1991], it is still possible to express all sums over intermediate states through it and therefore to perform calculations of higher-order perturbative corrections.

ODS is grateful to A. Leonau, S. Bragin and U. Sinha for helpful discussions. This article comprises 