&=\frac{1}{n_r+2\gamma} \left(\kappa(n_r+\gamma) -
    \frac{\gamma}{\chi}\right),
  \\
  \chi
  &=(n_r^2+\kappa^2+2n_r\gamma)^{-\frac{1}{2}},
  \\
  \gamma
  &= \sqrt{\kappa^2-(Z^* \alpha)^2},\end{aligned}$$ and $N$ is the corresponding normalization constant $$\begin{aligned}
  N=\frac{\chi}{\Gamma(2\gamma+2)} \sqrt{\Gamma(2\gamma+n_r+1)
  \frac{(\kappa+\gamma)}{n_r!} \frac{\chi Z^*}{\rho}}.\end{aligned}$$

The quantity $\kappa$ is defined as $$\begin{aligned}
  \kappa = \left\{
  \begin{aligned}
    l, \quad \mathrm{if}\ j = l - \frac{1}{2},
    \\
    -(l+1), \quad  \mathrm{if}\ j = l + \frac{1}{2},
  \end{aligned}\right. \label{eq:31}\end{aligned}$$ and $\Gamma(z)$ is the Gamma function. In addition, we employ the spherical harmonics with a negative first index, which is used in Mathematica [@Mathematica] and is defined by the following identity $Y_{l}^{m} = Y_{-(l+1)}^{m}$ for $l \leq -1$ and $(l+1) \leq m \leq -(l+1)$ for $l \leq -1$.

We also mention here that for the convenience we expressed the Dirac wave functions through the Whittaker functions and not through more commonly used hypergeometric ones [@flugge_practical_1971]. However, our definition in Eq. (\[eq:30\]) is in full agreement with the commonly used wave functions from [@flugge_practical_1971].

Now, by noting that $$\int M_{a+\gamma,\gamma-1/2}(r) M_{a+\gamma,\gamma-1/2}(r)
  \frac{dr}{r} = \frac{\Gamma(2\gamma)^2 a!}{\Gamma(a-2\gamma)},$$ we get (after lengthy but straightforward simplifications) $$\begin{aligned}
  A_{n_r \kappa}(Z^*) =
  Z^*\chi ^3\left(n_r+\frac{\kappa^2}{\gamma}\right). \label{eq:32}\end{aligned}$$

Furthermore, in order to evaluate $B_{\nu_2,\nu_4}^{\nu_1,\nu_3}$, we employ the expansion of the electron-electron interaction potential in spherical harmonics [@landau1981quantum] $$\begin{aligned}
  \frac{1}{|{\boldsymbol}{r}-{\boldsymbol}{r}'|} = \sum_{l=0}^\infty \sum_{m=-l}^{l}
  \frac{4\pi}{2l+1} \frac{r_<^l}{r_>^{l+1}} Y_{l}^{m*}(\Omega)
  