charge of the nucleus, we can easily obtain expressions for electronic densities of any atom or ion. Using the density operator: $\hat{\rho}({\boldsymbol}r) = \hat{\Psi}^\dag({\boldsymbol}r)\hat{\Psi}({\boldsymbol}r)$, we get the zeroth-order density of a given multi-electron state as: $$\begin{aligned}
  \rho^{(0)}({\boldsymbol}r)
  &= \langle \lambda_1, \ldots, \lambda_N | \hat{\rho}({\boldsymbol}{r})
  |\lambda_1, \ldots, \lambda_N\rangle \nonumber
  \\
  &= \sum_{k=1}^{N}
    |\Psi_{\lambda_k}^{\epsilon}(Z^{*}, {\boldsymbol}{r})|^2, \label{eq:14}\end{aligned}$$ with the explicit expressions for hydrogen-like wave functions $\Psi_{\lambda_k}^{\epsilon}({\boldsymbol}r)$ given in Appendix \[sec:appendix-a.-details\]. In fact, expression (\[eq:14\]) represents very simple products of exponentials and polynomials, thus allowing its use in numerical plasma [@chung_flychk:_2005] or description of ionization in particle-in-cell (PIC) computer codes that describe laser-matter interactions [@arber_contemporary_2015].

Accuracy of the zeroth-order approximation {#sec:zeroth-order-appr}
==========================================

Ground state energies and electronic densities of neutral atoms {#sec:ground-state-energ}
---------------------------------------------------------------

In order to calculate the ground-state energies of multi-electron atoms or ions, one needs to specify a set of occupation numbers, characterizing the state within our Dirac hydrogen-like basis. However, in the zeroth-order approximation the particular choice may not have the correct ordering due to relatively low ionization energies of heavy atoms.

The comparison of our analytical expressions obtained via Eq. (\[eq:9\]) for the values of ground state energies with the results obtained via solutions of DHF equations demonstrates that similarly 