structure of the interaction operator $\hat{W}$.

a\) b) c) d)

To find the effective charge we proceed in analogy with Ref. [@skoromnik_analytic_2017] and choose it from the condition that the first-order correction to the energy of the system for a given state is vanishing, i.e., $$\begin{aligned}
  \Delta E^{(1)}(Z^*) = 0. \label{eq:13}\end{aligned}$$ For this reason the expression for the energy of the system in first-order perturbation theory is given via a sum of hydrogen-like energies, Eq. (\[eq:9\]) with the effective charge $Z^{*}$, defined as a solution of Eq. (\[eq:13\]).

It is worth noting here, that the nontrivial dependence of Dirac hydrogen wave functions on the nuclear charge makes it impossible to separate the effective charge from the above integrals. This means that contrary to the nonrelativistic case, $A_{\lambda_{k}}$ and $B^{\lambda_{1},\lambda_{2}}_{\lambda_{3},\lambda_{4}}$ are implicitly dependent on $Z^*$. This is related to the fact that the Dirac equation, unlike the Schrödinger equation, is not scale invariant. In fact, rescaling the radial variable $r$ in a Schrödinger hydrogen atom, effectively changes its charge, while in the Dirac hydrogen atom, it effectively changes its mass. Therefore, in order to solve Eq. (\[eq:13\]) a numerical method has to be used.

Following the procedure of calculating the effective 