We point out here that if we take the energy of the system in first-order perturbation theory and consider the effective charge to be equal to the full nuclear charge, one obtains the perturbation theory over $Z^{-1}$. For this reason, we have compared the values of energies from the effective charge model with the ones from $Z^{-1}$ expansion. For example, for He one gets $-2.75$ a.u vs $-2.86$ a.u using $Z^{-1}$ expansion and effective charge model respectively, while for larger atoms the accuracy drops significantly.

In order to compute the electronic density, one needs to use Eq. (\[eq:14\]), i.e., to compute the sum of squares of absolute values of Dirac hydrogen-like wave functions. After simplification, this sum is just a product of exponentials and polynomials. Consequently, our method, unlike the DHF and TF calculations, provides fully analytical expressions for electronic densities and is therefore particularly useful for applications requiring repeated calculations. Relatively simple expressions resulting from our model can, therefore, be incorporated into existing software, used in the description of X-ray scattering on Mössbauer crystals [@Sturhahn2000].

In Fig. \[fig:densities\_neutral\_atoms\] we plot the resulting dependence of electronic densities on the radial coordinate $r$ for selected neutral atoms. It is immediately clear 