energies satisfy the standard anticommutation relations $\{a_\nu, a^\dag_{\nu'}\}=\{b_\nu, b^\dag_{\nu'}\}=\delta_{\nu,\nu'}$, which ensures that multi-particle states $|\nu_1...\nu_N\rangle$ correspond to fully antisymmetric Slater determinants. In the following, we will denote electron occupation numbers with $\lambda_{i}$ and negative energy occupation numbers as $\mu_{i}$.

In order to compute the energy of the system we specify a set of $N$ occupation numbers $\lambda_{1}, \ldots, \lambda_{N}$, which characterizes the state of $N$ electrons $|\lambda_{1},\ldots,\lambda_{N}\rangle$ and construct a perturbation series by considering the operator $\hat{W} = \hat{W}_{1} + \hat{W}_{2}$ as a perturbation $$\begin{aligned}
  E
  &= E^{(0)} (Z^*) + \Delta E^{(1)} (Z^*) + ... , \label{eq:6}
  \\
  \Delta E^{(1)} (Z^*)
  &= \langle\lambda_{1},\ldots,\lambda_{N}| \hat{W}
  |\lambda_{1},\ldots,\lambda_{N}\rangle. \label{eq:7}\end{aligned}$$

Since the state $|\lambda_{1},\ldots,\lambda_{N}\rangle$ is the eigenstate of $\hat{H}_{0}$ the zeroth-order energy $E^{(0)}(Z^{*})$ is a sum of hydrogen-like Dirac energies over the occupied states. The expression for the Dirac energy of the hydrogen-like atom is well-known [@flugge_practical_1971] and is given by $$E^D_{n_{r}jm_j} (Z^*) = \frac{c^{2}}{\sqrt{1 +
      \left(\frac{\alpha Z^{*}}{n_{r} + \sqrt{(j+1/2)^{2} - (\alpha
            Z^{*})^{2}}}\right)^{2}}}. \label{eq:8}$$ For this reason the energy $E^{(0)}(Z^{*})$ reads $$\begin{aligned}
  E^{(0)}(Z^{*}) = \sum_{{\lambda_{i}}}
  E_{n_{r}jm_{j}}^{D}(Z^{*}). \label{eq:9} \end{aligned}$$

Evaluation of the first-order correction to the energy of the system is straightforward: $$\begin{aligned}
  \Delta E^{(1)} (Z^*)
  &= \langle \lambda_{1},...,\lambda_{N} | \hat{W}
    |\lambda_{1},...,\lambda_{N}\rangle \nonumber
  \\
  &= 