---
abstract: 'A relativistic effective charge model has been developed for computation of observable characteristics of multi-electron atoms and ions. A complete and orthogonal Dirac hydrogen basis set, depending on one parameter — effective nuclear charge $Z^{*}$ — identical for all single-electron wave functions of a given atom or ion, is employed for the construction of the secondary-quantized representation. The effective charge is uniquely determined by the charge of the nucleus and a set of electron occupation numbers for a given state. We thoroughly study the accuracy of the leading-order approximation for the total binding energy and demonstrate that it is independent of the number of electrons of a multi-electron atom. In addition, it is shown that the fully analytical leading-order approximation is especially suited for the description of highly charged ions since our wave functions are almost coincident with the Dirac-Hartree-Fock ones for the complete spectrum. Finally, we envision that our effective charge model is more accurate and thus can replace the Thomas-Fermi-Dirac model for all applications where it is still utilized.'
author:
- 'K. Dzikowski'
- 'O. D. Skoromnik'
- 'I. D. Feranchuk'
- 'N. S. Oreshkina'
- 'C. H. Keitel'
bibliography:
- 'relativistic\_zeroth.bib'
title: 'Relativistic effective charge model of a multi-electron atom'
---

Introduction {#sec:introduction}
============

Despite great advances in the accuracy and efficiency of modern numerical methods [@filippin_multiconfiguration_2016; @puchalski_relativistic_2017; @peach_fractional_2015; @jensen2007introduction], simple analytical models for describing multi-electron atoms and ions are still actively developed [@gomez_simple_2019; @tatewaki_relativistic_2018; @karpov_atomic-number_2017; @hey_forms_2017]. The most popular of them include the multi-parametric screening hydrogen orbitals [@hau2012high; @seaton_quantum_1983], Slater orbitals [@slater_atomic_1930] and the semi-classical Thomas-Fermi (TF) model [@thomas_1927; @landau1981quantum]. These models are widely used, whenever certain level of precision has to be sacrificed for improved computation time, e.g. in computational plasma physics [@ciricosta_direct_2012; @starrett_thomas-fermi_2017; @lee_model_1987; @scott_cretinradiative_2001; @chung_flychk:_2005; @dyachkov_region_2016], semiconductors [@schulze-halberg_position-dependent_2013; @smith_electronic_2014], screening effects in the cross sections of bremsstrahlung and pair production [@tsai_pair_1974; @olsen_photon_1959; @davies_theory_1954; @seltzer_bremsstrahlung_1985; @poskus_brems:_2018], X-ray scattering and diffraction [@hau2012high; @feranchuk_new_2002] and crystallography [@toraya_new_2016].

In addition, the wave functions obtained from these simple analytical models are used as the initial approximation for the solution of Hartree-Fock (HF) [@FockA1930naeherungs; @fischer1977hartree; @grant2007relativistic; @jonsson_new_2013] and post HF equations [@bostock_fully_2011], convergence of which is strongly dependent on the choice of the trial functions. Therefore, improving the choice of the initial wave functions can significantly reduce the number of iterations required for obtaining the desired solution or enable the convergence at all.

At the same time a lot of effort has gone into improving the accuracy of these models, while keeping their complexity low. For example, the inclusion of relativistic corrections [@dolg_relativistic_2012; @nakajima_douglaskrollhess_2012; @marini_relativistic_1981; @waber_relativistic_1975] and inhomogeneity corrections [@chau_systematic_2018] into the TF model was performed. However, solving the TF equation to produce numerically stable results of electronic densities is in general nontrivial problem and is an active area of research [@akgul_constructing_2017; @parand_accurate_2017; @liu_laguerre_2015].

Consequently, finding a set of analytical wave functions, which provide better accuracy than the models described above would be advantageous. As was recently demonstrated for nonrelativistic atoms or ions, this can be achieved within the effective charge model [@skoromnik_analytic_2017], where the hydrogen-like wave functions can be used to analytically describe observable characteristics with high accuracy. In this approach, one specifies a complete and orthonormal hydrogen-like basis set with a single parameter — effective nuclear charge $Z^{*}$ — identical for all single-electron wave functions of a given atom. The basis completeness allows one to perform a transition into the secondary quantized representation, which is especially suited for the description of many-body problems.

Then in order to compute observable characteristics one specifies the set of occupation numbers for the state in question and the charge of the nucleus. This allows one to determine the effective charge from the vanishing first-order correction to the energy of the system analytically. Then observable characteristics are expressed through the expectation values, computed with the hydrogen-like wave functions of this charge.

We point out here that perturbation theory in $Z^{-1}$ was developed in many works (see [@vainshtein_energy_1985; @fischer1977hartree] and citations therein), however, exactly the introduction of the effective charge $Z^{*}$, instead of the usage of the nuclear charge $Z$, significantly increases the accuracy of the leading zeroth-order approximation, while keeping the complexity of all calculations low as it contains the first-order correction implicitly.

As a result, it was demonstrated that the analytical leading-order approximation describes the whole spectrum of a nonrelativistic multi-electron atom or ion and the associated perturbation theory series was constructed, which converges fast with the rate of $\sim 1/10$. The accuracy of these results does not depend on the number of electrons in an atom, i.e., the effective charge model is uniformly available for all states, both ground and excited ones, of all atoms or ions. Moreover, the results via second-order perturbation theory are comparable with those via multi-configuration HF (MCHF) [@fischer1977hartree].

We emphasize that the effective charge model is a completely *ab-initio* theory and thus crucially distinct from semi-empirical models, such as those based on quantum defects of Rydberg atoms [@foot2005atomic], screened hydrogen [@hau2012high] or empirical interpolation of X-ray measurements [@genoni_can_2017]. In addition, contrary to numerical methods such as solving the Hatree-Fock equations [@hartree_1928], the effective charge model provides analytical expressions for electron densities and scattering factors.

So far, the effective charge description has only been performed for the nonrelativistic Schrödinger equation. However, it is well-known that when a product of a fine structure constant $\alpha$ and the charge of the nucleus $Z$, i.e., $\alpha Z$ is not much smaller than unity the relativistic corrections become important. Since the Dirac equation for the hydrogen-like atom can be solved exactly and the Dirac-Coulomb Green’s function is known in closed form, it is possible to develop a relativistic version of the effective charge model. Consequently, this generalization is the main achievement of the present work.

Similarly to the nonrelativistic case, we employ Dirac hydrogen basis set with a single parameter — the effective charge and require it to be identical for all wave functions of a given atom. The basis completeness allows us to work in the secondary-quantized representation. We characterize states via electron occupation numbers and determine the value of the effective charge, by requiring the first-order perturbative correction to the energy of the system to vanish. With this we obtain a description of energy spectra, both ground and excited states together with wave functions. We have shown that the analytical description within the effective charge model agrees with the high accuracy with the numerical solution of Dirac-Hartree-Fock (DHF) equations for all ions and atoms of the periodic table.

The article is organized in the following way: in Sec. \[sec:relat-effect-charge\] we introduce the effective charge description and construct a perturbation theory series. In Sec. \[sec:zeroth-order-appr\] we analyze the accuracy of the zeroth- and first-order approximations. Consequently, we start from computing numerical values of effective charges by solving Eq. (\[eq:13\]). After this in Sec. \[sec:ground-state-energ\] we determine the ground state energies for 118 atoms of the periodic table and the electronic densities for some selected atoms. In Sec. \[sec:atom-scatt-fact\] we analyze atomic scattering factors. In Sec. \[sec:highly-charged-ions\] we demonstrate that the leading-order approximation is especially suited for the description of highly charged ions, where our analytical functions are almost coincident with numerical DHF wave functions. In Sec. \[sec:conclusions\] we provide conclusions and outlook. Finally, the details of the calculations are summarized in Appendices \[sec:appendix-a.-details\]-\[sec:solut-relat-tf\] and values of the energies are given in Appendix \[sec:values-ground-state\] and the atomic units with $\hbar = m = e
= 1$ are employed through the paper.

Relativistic effective charge model {#sec:relat-effect-charge}
===================================

As was discussed in the introduction, we can write down the Hamiltonian of the relativistic multi-electron atom in the secondary-quantized representation, using the complete and orthonormal Dirac hydrogen-like basis set, with a single parameter — effective charge $Z^{*}$ — identical for all wave functions of the basis set. Consequently, the Hamiltonian of the system reads $$\begin{aligned}
  \hat{H}
  &= \hat{H}_0 + \hat{W} = \hat{H}_{0} + \hat{W}_{1} +
    \hat{W}_{2},\label{eq:1}
  \\
  \hat{H}_0
  &=\int \hat{\Psi}^\dag({\boldsymbol}{r}) \left(c{\boldsymbol}{\alpha} \cdot
    \hat{{\boldsymbol}{p}} + \beta c^{2} - \frac{Z^*}{|{\boldsymbol}{r}|}\right)
    \hat{\Psi}({\boldsymbol}{r})d{\boldsymbol}{r}, \label{eq:2}
  \\
  \hat{W}_1
  &=\int  \hat{\Psi}^\dag({\boldsymbol}{r}) \frac{Z^*-Z}{|{\boldsymbol}{r}|}
    \hat{\Psi}({\boldsymbol}{r})d{\boldsymbol}{r}, \label{eq:3}
  \\
  \hat{W}_2
  &= \frac{1}{2}\int \frac{\hat{\Psi}^\dag({\boldsymbol}{r})
    \hat{\Psi}^\dag({\boldsymbol}{r'}) \hat{\Psi}({\boldsymbol}{r})
    \hat{\Psi}({\boldsymbol}{r'})}{|{\boldsymbol}{r}-{\boldsymbol}{r}'|}
    d{\boldsymbol}{r}d{\boldsymbol}{r'}, \label{eq:4}\end{aligned}$$ where $c$ is the speed of light ($c = 1/\alpha = 137.035999084$), $Z$ the charge of the nucleus, $Z^{*}$ the effective charge, which will be determined later, $\hat{{\boldsymbol}p} = -{\mathrm{i}}\boldsymbol{\nabla}$ the momentum operator, ${\boldsymbol}\alpha$, $\beta$ the Dirac matrices and $\hat\Psi({\boldsymbol}r)$ the secondary-quantized wave function [@FeynmanB1998statistical]. Furthermore, we split the perturbation operator $\hat{W}$ into the sum of single-electron $\hat{W}_1$ and double-electron $\hat{W}_2$ operator components. Moreover, we do not take into account the Breit part in the potential of the electron-electron interaction and radiative quantum electrodynamics corrections, thus limiting ourselves to the description of the multi-electron atom within DHF approximation. However, we mention here that for large $Z$ these corrections become relevant and would have to be included for all models. We stress here that Eq. (\[eq:1\]) is the exact transformation of the Hamiltonian of the multi-electron atom, since we just added and subtracted the term proportional to $Z^{*}$. The eigenvalues of $\hat{H}$ define the values of the energy $E$ of a multi-electron atom.

The secondary-quantized wave function $\hat{\Psi}({\boldsymbol}r)$ is expanded in the complete Dirac hydrogen-like basis, which is dependent on $Z^{*}$ $$\hat{\Psi}({\boldsymbol}r) = \sum_\nu\left( \hat{a}_\nu \Psi_\nu^\epsilon(Z^{*},
    {\boldsymbol}r) + \hat{b}^\dag_\nu  \Psi_\nu^{-\epsilon}(Z^{*}, {\boldsymbol}r)\right), \label{eq:5}$$ with $\Psi^\epsilon_{\nu}(Z^{*}, {\boldsymbol}r)$ and $\Psi^{-\epsilon}_{\nu}(Z^{*}, {\boldsymbol}r)$ being the positive and negative energy eigenstates of the single-particle Dirac Hamiltonian respectively, identified by the set of collective quantum numbers $\nu=n_r l j m_{j}$ and $\nu = {\boldsymbol}p l j m_{j}$ for the discrete and continuous spectra correspondingly (actual expressions for the functions are given in Appendix \[sec:appendix-a.-details\]). The fermionic operators of positive and negative 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
  \\
  &= (Z^*-Z)\sum_{k = 1}^{N} A_{\lambda_k}(Z^*) \label{eq:10}
  \\
  &+ \sum_{k<l = 1}^{N}
    \left(B_{\lambda_k,\lambda_l}^{\lambda_k,\lambda_l}(Z^*) -
    B_{\lambda_k,\lambda_l}^{\lambda_l,\lambda_k}(Z^*)\right),
    \nonumber\end{aligned}$$ where we have defined: $$\begin{aligned}
  A_{\lambda_{k}}
  &= \int
    \frac{|\Psi_{\lambda_{k}}^{\epsilon}({\boldsymbol}{r})|^2}{|{\boldsymbol}{r}|}
    d{\boldsymbol}{r}, \label{eq:11}
  \\
  B_{\lambda_3,\lambda_4}^{\lambda_1,\lambda_2}
  &=
    \int\frac{\Psi_{\lambda_1}^{\epsilon*}({\boldsymbol}{r})
    \Psi^{\epsilon*}_{\lambda_3}({\boldsymbol}{r}')
    \Psi_{\lambda_2}^{\epsilon}({\boldsymbol}{r})
    \Psi_{\lambda_4}^{\epsilon}({\boldsymbol}{r}')}{|{\boldsymbol}{r}-{\boldsymbol}{r}'|}
    d{\boldsymbol}{r} d{\boldsymbol}{r}', \label{eq:12}\end{aligned}$$ with the implied dependence on the effective charge ommited for simplicity. The $A_{\lambda_{k}}$ describes the contribution from the single-particle operator $\hat{W}_{1}$, while $B_{\lambda_3,\lambda_4}^{\lambda_1,\lambda_2}$ from the double-particle one $\hat{W}_{2}$. The last term in Eq. (\[eq:10\]) is the difference of the Coulomb and exchange integrals. Both $A_{\lambda_{k}}$ and $B_{\lambda_3,\lambda_4}^{\lambda_1,\lambda_2}$ can be evaluated analytically for a given set of occupation numbers (See Appendix \[sec:appendix-a.-details\]).

![(Color online) The relative difference between the values of ground state energies of the first 118 neutral atoms obtained via the effective charge model ($E^{(0)}$) and by solving the Dirac-Hartree-Fock equations ($E_{\mathrm{DHF}}$)[@DESCLAUX1973311].[]{data-label="accuracy"}](error_comparison.pdf){width="80mm"}

In addition, we mention here that fermionic operators of negative energies do not contribute to the energy of the system in the zeroth- and first-order as we are considering the corrections only to the electronic states, that is the states that do not contain $\mu_{i}$ occupation numbers. However, starting from second-order perturbation theory, the negative energy states will contribute to the observable characteristics, since there exist nonvanishing matrix elements due to the 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 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 to the nonrelativistic case, the optimal choice of the occupation numbers is given according to the “Aufbau” or Madelung–Janet–Klechkovskii rule [@Madelung1936; @KlechkovskiiA1962justification; @doi:10.1021/ed056p714], as it provides a simple and consistent choice of occupation numbers for any number of electrons, while resulting in the lowest zeroth-order ground state energies in almost all cases. For example, from the three sets of occupation numbers $[\mathrm{Xe}]6\mathrm{s}^{1}$, $[\mathrm{Xe}]4\mathrm{f}^{1}$ and $[\mathrm{Xe}]5\mathrm{d}^{1}$ for cesium ($Z=55$), the first state possesses the lowest energy, while the remaining two are excited states. Compare the values for the energy $E_{[\mathrm{Xe}]6\mathrm{s}^{1}} = -7361.18$ a.u versus $E_{[\mathrm{Xe}]4\mathrm{f}^{1}} = -7346.19$ a.u and $E_{[\mathrm{Xe}]5\mathrm{d}^{1}} = -7354.40$ a.u, which corresponds to the “Aufbau” rule.

The accuracy of the zeroth-order calculation is presented in Fig. \[accuracy\] and numerical values of effective charges and energies are given in Appendix \[sec:values-ground-state\] in Tab. \[tab:ground\_state\_energies\]. We compare our results with solutions of DHF equations from Ref. [@DESCLAUX1973311]. As can be concluded from Fig. \[accuracy\] the effective charge description leads to a uniform approximation, i.e., the relative accuracy is independent of the number of electrons of an atom and $\sim 6.5\%$ with respect to DHF for all elements of the periodic table. This is considerably better than the TF approximation [@landau1981quantum].

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 that the zeroth-order approximation correctly reproduces all qualitative features of the corresponding DHF result and contrary to the relativistic TFD model provides density oscillations and a correct asymptotic behavior for large values of $r$. In addition, we point out that the relativistic TFD model [@dolg_relativistic_2012; @nakajima_douglaskrollhess_2012; @marini_relativistic_1981; @waber_relativistic_1975], unlike its nonrelativistic counterpart [@landau1981quantum], does not have a universal dependence on the charge of the nucleus. For this reason, in order to obtain the electronic density, the relativistic TFD equation needs to be repeatedly solved numerically, which is a nontrivial procedure. (See Appendix \[sec:solut-relat-tf\]).

Atomic scattering factors {#sec:atom-scatt-fact}
-------------------------

Another observable characteristics that can be extracted from the effective charge model are the atomic scattering factors. According to their definition [@hau2012high], they are expressed by the Fourier transforms of electronic density [@prince_international_2006]: $$\begin{aligned}
  f({\boldsymbol}{q}) = \int \rho({\boldsymbol}{r})e^{i {\boldsymbol}{q} \cdot {\boldsymbol}{r}}
  d{\boldsymbol}{r}, \label{eq:15}\end{aligned}$$ which in our approach can be calculated analytically (see Appendix \[sec:appendix-a.-details\]).

The atomic scattering factors are very important for crystallography and X-ray physics, since the crystal polarizability $\chi$ as the function of X-ray frequency $\omega_{\mathrm{r}}$, can be evaluated by employing the relation [^1] [@ahmadi_parametric_2013; @baryshevsky2005parametric] $$\begin{aligned}
  \chi({\boldsymbol}g, \omega_{\mathrm{r}}) = \frac{4 \pi S({\boldsymbol}g)}{\Omega_{0}
  \omega_{\mathrm{r}}^{2}}
  \left(-\frac{e^{2}}{m}f({\boldsymbol}g)\right), \label{eq:17}\end{aligned}$$ where $\Omega_0$ is the volume of a crystal cell, ${\boldsymbol}g$ the reciprocal lattice vector, $S({\boldsymbol}g)$ the structure factor of the crystal, $m$ the electron mass and $e$ the electron charge.

In Fig. \[fig:ASF\] we present the results for gold, lead, uranium and xenon. Our analytical expressions for atomic scattering factors are comparable to the DHF calculation to within $25\%$.

As a realistic example, we also evaluated the relative difference for the emission intensity of parametric X-ray radiation (PXR) [@feranchuk_theoretical_1985; @baryshevsky2005parametric; @skoromnik_radical_2017] from relativistic electrons in a tungsten crystal. The intensity of PXR is proportional to the square of the absolute value of crystal polarizability. Therefore, the relative difference for the intensity of radiation with crystal polarizabilities from DHF calculation and from effective charge model is given by $(|\chi^{(0)}|^{2} -
|\chi_{\mathrm{DHF}}|^{2})/|\chi_{\mathrm{DHF}}|^{2}$. Consequently, for the relevant range of frequencies and Bragg’s angles the relative accuracy was $\sim 20\%$.

Highly charged ions {#sec:highly-charged-ions}
===================

The effective charge model is also suitable for the description of ions, since it does not require the charge of the nucleus to be equal to the number of electrons of an atom. Consequently, we fix the nuclear charge as $Z$ and use a state vector $|\lambda_{1},\ldots,\lambda_{N}\rangle$, where $N \neq Z$. As an example we evaluated energies of ground and excited states of He-, Li- and B-like uranium ($Z = 92$) ions.

We want to stress here again, that the effective charge $Z^{*}$ is identical for all single particle states from a set that specifies the state. However, two different states will have different charges. For example, the He-like uranium ground state $1s_{\uparrow}1s_{\downarrow}$ and the He-like uranium excited state $1s_{\uparrow}2s_{\downarrow}$ possess unequal effective charges.

The results of the calculation of electronic densities are presented in Fig. \[fig:HCI\]. The electronic densities, despite being analytical expressions in the zeroth-order approximation, coincide remarkably well with the ones obtained from numerical solutions of the DHF equations.

Furthermore, the energy eigenvalues are presented in Tab. \[tab:Ions\]. It is clear that the accuracy is significantly better (below $0.02\%$) than for neutral atoms and in fact sufficient to obtain correct ordering, even for very closely spaced excited states. It is worth pointing out that the accuracy can be further improved by forming linear combinations of all sub-configurations arising from the single $JLS$ configuration and subsequently diagonalizing the Hamiltonian exactly in such finite basis — a technique that has been successfully employed in the nonrelativistic case [@skoromnik_analytic_2017].

Finally, we note that $Z^{-1}$ expansion gives reasonable results for highly charged ions. However, it is less accurate than the effective charge model by at least $50\%$.

Conclusions and Outlook {#sec:conclusions}
=======================

We have demonstrated that the relativistic effective charge model describes multi-electron atoms and ions and provides analytic expressions for energies, electronic densities and scattering factors with accuracy comparable to the numerical solutions of Dirac-Hartree-Fock equations. We have also shown that our approach is suitable for the description of arbitrary excited states and provides accuracy independent of the number of electrons already in the zeroth-order approximation. Furthermore, the zeroth-order approximation can be easily modified to include interactions with external fields.

[|c|c|c|c|c|c|c|]{} Configuration, $J^{\pm}$ & $Z^{*}$ & $E^{(0)}$ & $E_{\mathrm{DHF}}$ & $\Delta E/E_{\mathrm{DHF}} \cdot 100\%$ & CI\
$1s^2,J=0^+$ & 91.7130 & -9651.35 & -9651.39 & 0.0004$\%$&-9651.45\
$1s^1 2s^1,J=1^+$ & 91.8623 &-6096.96 & -6097.01 & 0.0008$\%$&-6097.01\
$1s^1 2s^1,J=0^+$ & 91.8404 &-6093.51 & -6093.41 & -0.0016$\%$&-6097.01\
$1s^1 2p^1_{1/2},J=0^-$ & 91.8148 & -6089.50 & -6090.16 &0.0108$\%$ & -6090.17\
$1s^1 2p^1_{1/2},J=1^-$ & 91.8113 & -6088.95 & -6089.10 & 0.0025$\%$ & -6089.11\
$1s^1 2p^1_{3/2},J=2^-$ & 91.8487 & -5928.36 & -5928.44 & 0.0013$\%$ & -5928.44\
$1s^1 2p^1_{3/2},J=1^-$ & 91.8395 & -5927.00 & -5926.68 &-0.0054$\%$ & -5928.44\
$1s^1 3s^1,J=1^+$ & 91.9086 & -5389.83 & -5389.88 & 0.0009$\%$ & -5389.88\
$1s^1 3s^1,J=0^+$ & 91.9177 & -5388.95 & -5388.73 & -0.0041$\%$ & -5389.88\
$1s^2 2s^1,J=\frac{1}{2}^+$ & 91.5805 &-10862.2 & -10862.4 & 0.0018$\%$&-10862.5\
$1s^2 2p^1_{1/2},J=\frac{1}{2}^-$ & 91.5378 & -10850.3 & -10849.1 & -0.0104$\%$&-10850.8\
$1s^2 2p^1_{3/2},J=\frac{3}{2}^-$ & 91.5685 & -10694.8 & -10693.9 & -0.0084$\%$&-10695.2\
$1s^2 3s^1,J=\frac{1}{2}^+$ & 91.6447 & -10168.7 & -10168.9 & 0.0020$\%$&-10169.0\
$1s^2 3p^1_{1/2},J=\frac{1}{2}^-$ & 91.6323 & -10165.5 & -10165.4 & -0.0015$\%$ &-10165.8\
$1s^2 3p^1_{3/2},J=\frac{3}{2}^-$ & 91.6429 & -10119.2 & -10119.0 & 0.0004$\%$&-10119.4\
$1s^2 3d^1_{3/2},J=\frac{3}{2}^+$ & 91.6399 & -10118.4 & -10118.6 & -0.0020$\%$&-10118.7\
$1s^2 3d^1_{5/2},J=\frac{5}{2}^+$ & 91.6435 & -10106.7 & -10106.9 & 0.0022$\%$ &-10107.0\
$1s^2 2s^2 2p^1_{1/2},J=\frac{1}{2}^-$ & 91.2080 & -13221.1 & -13222.5 & 0.0104$\%$&-13222.7\
$1s^2 2s^1 2p^2_{1/2},J=\frac{1}{2}^+$ & 91.1604 & -13204.8 & -13206.7 & 0.0140$\%$&-13206.9\
$1s^2 2s^2 2p^1_{3/2},J=\frac{3}{2}^-$ & 91.2339 & -13068.9 & -13070.2 & 0.0099$\%$&-13070.6\
$1s^2 2s^1 2p^1_{1/2} 2p^1_{3/2},J=\frac{3}{2}$ & 91.1972 & -13056.6 & -13058.8 & 0.0168$\%$&-13058.9\
$1s^2 2s^1 2p^1_{1/2} 2p^1_{3/2},J=\frac{5}{2}$ & 91.1914 & -13054.7 & -13056.3 & 0.0128$\%$&-13056.5\
$1s^2 2p^2_{1/2} 2p^1_{3/2},J=\frac{3}{2}^-$ & 91.1327 & -13035.1 & -13037.2 & 0.0159$\%$ &-13037.1\
$1s^2 2s^2 3s^1,J=\frac{1}{2}^+$ & 91.3255 & -12556.2 & -12557.2 & 0.0080$\%$&-12557.7\

Exactly the introduction of the effective charge $Z^{*}$, instead of the usage of the nuclear charge $Z$, i.e., $Z^{*} \neq Z$ significantly increases the accuracy of the zeroth-order approximation, while rendering the complexity of all calculations low.

However, the fully relativistic description becomes much more complicated due to the structure of Dirac equation. First, the Dirac equation possesses negative energy states and, therefore, there exist matrix elements between electronic states and negative energy ones. This significantly increases the number of required matrix elements to be taken into account. Second, due to the mass term of the Dirac equation the energy of the system does not have a universal behavior of the quadratic function of the effective charge as in the nonrelativistic case. Therefore, it is impossible to separate explicitly the dependence of the matrix elements on the effective charge. For this reason, before investigating the second-order correction to the energy of the system we focused on the zeroth-order approximation first.

At the same time, we point out that since the Dirac-Coulomb Green’s function is known in analytical closed-form [@wong_dirac_1985; @swainson_unified_1991], it is still possible to express all sums over intermediate states through it and therefore to perform calculations of higher-order perturbative corrections.

ODS is grateful to A. Leonau, S. Bragin and U. Sinha for helpful discussions. This article comprises parts of the PhD thesis work of Kamil Dzikowski to be submitted to the Heidelberg University, Germany.

Details of the calculation {#sec:appendix-a.-details}
==========================

Here, we present all formulas needed for the fully analytical evaluation of all integrals encountered in the calculation of the first-order correction to the energy of the system, as well as the computation of atomic scattering factors and electronic densities.

The key idea, allowing for the analytic calculation is the observation that both the basis wave functions and all appearing potentials can be split into their radial and angular components by simple expansion in spherical harmonics. Thus we get the hydrogen-like Dirac wave functions as [@flugge_practical_1971] $$\begin{aligned}
  \Psi_{n_{r}ljm}
  &({\boldsymbol}{r},Z^*) \nonumber
  \\
  &=
  \frac{1}{r} \left(
  \begin{array}{c}
    G_{n_{r},\kappa}(r,Z^*)\Omega_{\kappa,m}(\theta,\phi)
    \\
    F_{n_{r},\kappa}(r,Z^*)\Omega_{-\kappa,m}(\theta,\phi)
  \end{array}
  \right), \label{eq:28}\end{aligned}$$ with the angular part given by spherical harmonics $$\begin{aligned}
  \Omega_{\kappa,m}
  &(\theta,\phi) \nonumber
  \\
  &= \left(
  \begin{array}{c}
    \sqrt{\frac{1}{2} - \frac{m}{2\kappa+1}}
    Y_\kappa^{m-1/2}(\theta,\phi)
    \\
    -\frac{\kappa}{|\kappa|}\sqrt{\frac{1}{2} + \frac{m}{2\kappa+1}}
    Y_\kappa^{m+1/2}(\theta,\phi)
  \end{array}
  \right), \label{eq:29}\end{aligned}$$ and the radial part by the Whittaker functions $M_{\kappa,\mu}(z)$ of the first kind [@whittaker_expression_1903]

$$\begin{aligned}
    \begin{pmatrix}
      G_{n_r,\kappa}(Z^{*}, r)
      \\
      F_{n_r,\kappa}(Z^{*}, r)
    \end{pmatrix}
     = N\left\{\frac{n_r}{2\gamma}\left(
      \begin{array}{c}
        Z^*\alpha
        \\
        i (\gamma-\kappa)
      \end{array}\right) M_{\gamma+n_r,\gamma+\frac{1}{2}}(2 \chi  Z^*r)
    -\rho(1+2\gamma)\left(
      \begin{array}{c}
        \frac{\gamma-\kappa}{Z^* \alpha}
        \\
        i\end{array}
    \right) M_{\gamma+n_r,\gamma-\frac{1}{2}}(2 \chi
    Z^*r)\right\}, \label{eq:30}
  \end{aligned}$$

where $$\begin{aligned}
  \rho
  &=\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)
  Y_{l}^{m}(\Omega')\end{aligned}$$ and note that the integration of three spherical harmonics yields $3j$ symbols [@Brink:104381]:

$$\begin{aligned}
    \int Y_{l_{1}}^{m_{1}}(\Omega) Y_{l_{2}}^{m_{2}}(\Omega)
    Y_{l_{3}}^{m_{3}}(\Omega) d\Omega =
    \sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4 \pi}}
    \begin{pmatrix}
      l_1 & l_2 & l_3
      \\
      0 & 0 & 0
    \end{pmatrix}
              \begin{pmatrix}
                l_1 & l_2 & l_3
                \\
                m_1 & m_2 & m_3
              \end{pmatrix}. \label{eq:34}
  \end{aligned}$$

This allows us to split the result into the radial and angular parts as $$\begin{aligned}
    B_{\nu_2,\nu_4}^{\nu_1,\nu_3} =
    -\delta_{m_{j2}-m_{j1}}^{m_{j3}-m_{j4}} (-1)^{2
    m_{j1}-m_{j2}-m_{j3}} \sum_p \left( \Phi_{\nu_1,\nu_2}^p \otimes
    \Phi_{\nu_3,\nu_4}^p \right) \cdot
    \sigma_{\nu_1,\nu_2,\nu_3,\nu_4}^p, \label{eq:35}
  \end{aligned}$$ with the angular part given by $$\begin{aligned}
    \Phi_{\nu_1,\nu_2}^p = \left(
    \begin{array}{c}
      \begin{pmatrix}
        \kappa_1 & \kappa_2 & p
        \\
        0 & 0 & 0
      \end{pmatrix}
                \left(
                \phi_{\kappa_1,m_{j1} -
                \frac{1}{2}}^{\kappa_2,m_{j2}-\frac{1}{2}}
                (p) -
                \phi_{\kappa_1,\frac{1}{2} -
                m_{j1}}^{\kappa_2,\frac{1}{2}-m_{j2}}
                (p)\right)
      \\
      \begin{pmatrix}
        -\kappa_1 & -\kappa_2 & p
        \\ 0 & 0 & 0
      \end{pmatrix}
      \left(\phi_{-\kappa_1,m_{j1} -
                   \frac{1}{2}}^{-\kappa_2,m_{j2}-\frac{1}{2}}(p) -
                   \phi_{-\kappa_1,\frac{1}{2} -
                   m_{j1}}^{-\kappa_2,\frac{1}{2}-m_{j2}}(p)\right)
    \end{array}\right), \label{eq:36}
  \end{aligned}$$ where $$\begin{aligned}
    \phi_{k_1,m_1}^{k_2,m_2}(p) = \sqrt{(k_1-m_1)(k_2-m_2)}
    \begin{pmatrix}
      \kappa_1 & \kappa_2 & p
      \\
      -m_{1} & m_{2} & m_{1}-m_{2}
    \end{pmatrix}. \label{eq:37}
  \end{aligned}$$

Furthermore, the radial part is given by $$\begin{aligned}
    \sigma_{\nu_1,\nu_2,\nu_3,\nu_4}^p = \int
    {{G_{\nu_1}^{*}(r)G_{\nu_2}(r)}\choose{F_{\nu_1}^{*}(r) F_{\nu_2}(r)}}
    \otimes
    {{G_{\nu_3}^{*}(r')G_{\nu_4}(r')}\choose{F_{\nu_3}^{*}(r') F_{\nu_4}(r')}}
    \frac{r_<^p}{r_>^{p+1}}drdr', \label{eq:38}
  \end{aligned}$$ where the dependence on $Z^*$ has been omitted for clarity. Eq. (\[eq:38\]) can be calculated analytically by noting that the integral of the four Whittaker functions reads [@gradstejn_table_2009] $$\begin{aligned}
    \int
    &M_{a_1+b_1,b_1-1/2}(q_1 r) M_{a_2+b_2,b_2-1/2}(q_2 r)
      M_{a_3+b_3,b_3-1/2}(q_3 r') M_{a_4+b_4,b_4-1/2}(q_4 r')
      \frac{r_<^l}{r_>^{l+1}} dr dr' \nonumber
    \\
    &=\sum_{i_1=0}^{a_1} \sum_{i_2=0}^{a_2} \sum_{i_3=0}^{a_3}
      \sum_{i_4=0}^{a_4} T_{{\boldsymbol}{a},{\boldsymbol}{b},{\boldsymbol}{q}}({\boldsymbol}{i})
      \left(f_{i_1+i_2+b_1+b_2+l+1}^{i_3+i_4+b_3+b_4-l}
      \left(\frac{q_3+q_4}{2},\frac{q_1+q_2}{2}\right)
      + f_{i_3+i_4+b_3+b_4+l+1}^{i_1+i_2+b_1+b_2-l}
      \left(\frac{q_1+q_2}{2},\frac{q_3+q_4}{2}\right)\right),
  \end{aligned}$$ where $$\begin{aligned}
    T_{{\boldsymbol}{a},{\boldsymbol}{b},{\boldsymbol}{q}}({\boldsymbol}{i}) =
    \prod_{k=1}^4{{a_k}\choose{i_k}}
    \frac{\Gamma(2b_k)}{\Gamma(2b_k+i_k)}
    (-1)^{i_k}q_k^{b_k+i_k}
  \end{aligned}$$ and we have made use of $$\begin{aligned}
    f_{a}^{b}(x,y) = \int_0^\infty \int_r^\infty e^{-\lambda r -
    \lambda' r'} r^{a-1}
    {r'}^{b-1} dr' dr=\frac{\Gamma(a+b)}{a{\lambda'}^{a+b}}
    {_2F_1}\left(a, a+b, a + 1, -\frac{\lambda}{\lambda'}\right).
  \end{aligned}$$ The properties of $f_{a}^{b}(x,y)$ functions are described in some detail in [@dzikowski_generating_2019]. The bold ${\boldsymbol}a$, ${\boldsymbol}b$, ${\boldsymbol}q$ and ${\boldsymbol}i$ are lists of four values, i.e., ${\boldsymbol}a = \{a_{1}, a_{2}, a_{3}, a_{4}\}$ with similar expressions for ${\boldsymbol}b$, ${\boldsymbol}q$ and ${\boldsymbol}i$.

Now we would like to compute the scattering factors. For this we consider the electronic density, which is averaged over the angular variables. By using the orthogonality of spherical spinors the angular dependence integrates out trivially. Therefore, we get the radial density as $$\begin{aligned}
    r^2 \rho_{n_r,\kappa} (r,Z^*) =
    |G_{n_r,\kappa}(r,Z^*)|^2 + |F_{n_r,\kappa}(r,Z^*)|^2, \label{eq:39}
  \end{aligned}$$ which gives the scattering factors as $$\begin{aligned}
    f_{n_r,\kappa}(q,Z^*) = \int \left(|G_{n_r,\kappa}(r,Z^*)|^2 +
    |F_{n_r,\kappa}(r,Z^*)|^2\right) e^{i {{\boldsymbol}q\cdot{\boldsymbol}r}}
    d{\boldsymbol}{r}. \label{eq:40}
  \end{aligned}$$ Expanding Whittaker functions in a finite series in Eq. (\[eq:30\]) and using [@gradstejn_table_2009] $$\begin{aligned}
    \int e^{-\alpha r}r^{n-2} e^{i {\boldsymbol}{q} \cdot {\boldsymbol}{r}} d{\boldsymbol}{r} =
    4\pi \Gamma(n) \frac{\sin(n
    \tan^{-1}(\frac{q}{\alpha}))}{\sqrt{(\alpha^2+q^2)^n}}.
    \label{eq:42}
  \end{aligned}$$ we get $$\begin{aligned}
    f_{n_r,\kappa}(q,Z^*) = (N(2\gamma+1)\Gamma(2\gamma))^2 \left(
    2\kappa (\kappa-\gamma) n_r^2\sigma_1 + 4(\kappa-\gamma) \rho n_r
    \sigma_2 +
    \frac{2\kappa}{\kappa+\gamma}\rho^2\sigma_3\right), \label{eq:41}
  \end{aligned}$$ where $$\begin{aligned}
    \sigma_{1}
    &= \sum_{\substack{i=1\\j=1}}^{n_r} {{n_r-1}\choose{i-1}}
    {{n_r-1}\choose{j-1}}
    \frac{\Gamma(i+j+2\gamma)}{\Gamma(2\gamma+i+1)!
    \Gamma(2\gamma+j+1)!} \xi_{i,j}(q,Z^*),
    \\
    \sigma_{2}
    &= \sum_{\substack{i=1\\j=0}}^{n_r} {{n_r-1}\choose{i-1}}
    {{n_r}\choose{j}} \frac{\Gamma(i+j+2\gamma)}{\Gamma(2\gamma+i+1)!
    \Gamma(2\gamma+j)!} \xi_{i,j}(q,Z^*),
    \\
    \sigma_{3}
    &= \sum_{\substack{i=0\\j=0}}^{n_r} {{n_r}\choose{i}}
    {{n_r}\choose{j}}
    \frac{\Gamma(i+j+2\gamma)}{\Gamma(2\gamma+i)!\Gamma(2\gamma+j)!}
    \xi_{i,j}(q,Z^*),
  \end{aligned}$$ and $$\begin{aligned}
    \xi_{i,j}(q,Z^*)
    &= \frac{(-1)^{i+j}}{q} \sin\left((i+j+2\gamma)
    \tan^{-1}\left(\frac{q}{2\chi Z^*}\right)\right) \left(\frac{2 \chi  
    Z^*}{\sqrt{(2\chi Z^*)^2+q^2}}\right)^{i+j+2\gamma}.
  \end{aligned}$$

Solution of the relativistic Thomas-Fermi equation {#sec:solut-relat-tf}
==================================================

In this appendix we describe the solution of the relativistic TF equation [@gilvarry_relativistic_1954]. The equation written in atomic units reads [@waber_relativistic_1975] $$\begin{aligned}
    x^{1/2}\chi''(x) = \chi^{3/2}(x)\left(1 +
    \left(\frac{128}{9\pi^{2}}\right)^{1/3} \frac{Z^{4/3}}{c^{2}}
    \chi'(x) \left(1 - \frac{x
    \chi'(x)}{2\chi(x)}\right)\right)^{3/2}, \label{eq:18}
  \end{aligned}$$

where $x = r / (b Z^{-1/3})$, $b = (9\pi^{2} / 128)^{1/3}$ and the dimensionless self-consistent potential $\chi(x)$ is related to the self-consistent potential of the TF model as $\phi(r) = Z \chi(rZ^{1/3}/b) - \phi_{0}$, with the constant $\phi_{0}$ defined from the normalization. For neutral atoms $\phi_{0}$ equals zero. For ions it is chosen such that the self-consistent potential vanishes not at infinity, but rather at some finite value $x_{c}$. In the nonrelativistic limit, i.e., when the speed of light tends to infinity the relativistic TF equation coincides with its nonrelativistic counterpart.

The TF equation must be complemented with boundary conditions, which for neutral atoms are given by [@waber_relativistic_1975] $$\begin{aligned}
  \chi(0) = 1, \quad \chi(\infty) = 0, \label{eq:19}\end{aligned}$$ and for ions [@marini_relativistic_1981] $$\begin{aligned}
  \chi(0) = 1, \quad -x_{c}\chi'(x_{c}) = 1 - N/Z, \label{eq:20}\end{aligned}$$ respectively. Here we also employed the fact that $\chi(x_{c}) = 0$.

As was mentioned in the introduction solution of the TF equation is a nontrivial mathematical problem since it represents a boundary value problem on a semi-infinite interval. In order to solve the equation, we used the shooting method. For neutral atoms we reformulated the boundary value problem as an initial value one $$\begin{aligned}
  \chi(0) = 1, \quad \chi'(0) = \mu, \label{eq:21}\end{aligned}$$ where $\mu$ represents a parameter. Consequently, we were seeking for the root of the equation $\chi(X, \mu) = 0$, where $X$ we changed from some small value to the very large one. For every $X$ we were solving Eq. (\[eq:21\]) by varying $\mu$. With this we obtained the following solutions $$\begin{aligned}
  \mu_{\mathrm{Xe}}
  &= -1.50965873266, \quad \chi(80,
  \mu_{\mathrm{Xe}}) \approx 10^{-6}, \label{eq:22}
  \\
  \mu_{\mathrm{U}}
  &= -1.49103044294, \quad \chi(80,
  \mu_{\mathrm{U}}) \approx 10^{-6} \label{eq:23}\end{aligned}$$ for atoms.

For ions we used a similar strategy, however, we “shot” from infinity. In this case the boundary value problem is already written as the initial value one $$\begin{aligned}
  \chi(x_{c}) = 0, \quad \chi'(x_{c}) =
  -\frac{1-N/Z}{x_{c}}. \label{eq:24}\end{aligned}$$

For this reason we simply varied the value of $x_{c}$ till the value of $\chi$ at zero becomes one. With this we got $$\begin{aligned}
  x_{c}
  &= 0.34635, \quad \chi(10^{-6}) \approx 1, \label{eq:25}
  \\
  x_{c}
  &= 0.47890, \quad \chi(10^{-6}) \approx 1. \label{eq:26}\end{aligned}$$

Finally, the density of the atom or ion is expressed through the self-consistent potential as

$$\begin{aligned}
    \rho(r) = \frac{8\sqrt{2}}{3\pi} \left(\frac{Z
    \chi(x)}{r} - \phi_{0}\right)^{3/2} \left(1 +
    \left(\frac{128}{9\pi^{2}}\right)^{1/3} \frac{Z^{4/3}}{c^{2}}
    \chi'(x) \left(1 - \frac{x
    \chi'(x)}{2\chi(x)}\right)\right)^{3/2}. \label{eq:27}
  \end{aligned}$$

Values of ground state energies {#sec:values-ground-state}
===============================

[|ccccc|ccccc|]{} $Z$ &$Z^{*}_{\mathrm{R}}$ & $E^{(0)}_{\mathrm{NR}}$ & $E^{(0)}_{\mathrm{R}}$ & $E_{\mathrm{DHF}}$ & $Z$ & $Z^{*}_{\mathrm{R}}$ & $E^{(0)}_{\mathrm{NR}}$ & $E^{(0)}_{\mathrm{R}}$ & $E_{\mathrm{DHF}}$\
1 & 1.00000 & -0.50000 & -0.50000 & -0.50000 & 51 & 40.7062 & -5974.00 & -6160.94 & -6475.24\
2 & 1.68749 & -2.84766 & -2.84772 & -2.86175 & 52 & 41.5615 & -6259.79 & -6461.59 & -6788.06\
3 & 2.54539 & -7.28906 & -7.28951 & -7.43327 & 53 & 42.4156 & -6553.19 & -6770.65 & -7109.76\
4 & 3.37163 & -14.2096 & -14.2121 & -14.5752 & 54 & 43.2653 & -6854.26 & -7087.18 & -7440.46\
5 & 4.15118 & -23.6936 & -23.7003 & -24.5350 & 55 & 44.1573 & -7165.57 & -7415.35 & -7779.91\
6 & 4.90693 & -36.2016 & -36.1296 & -37.6732 & 56 & 45.0484 & -7484.40 & -7752.08 & -8128.34\
7 & 5.64987 & -52.0662 & -51.8941 & -54.3229 & 57 & 45.7984 & -7804.64 & -8083.37 & -8485.87\
8 & 6.42240 & -71.2844 & -72.2209 & -74.8172 & 58 & 46.5481 & -8125.81 & -8423.60 & -8852.82\
9 & 7.17595 & -94.4525 & -96.6125 & -99.4897 & 59 & 47.2966 & -8447.55 & -8772.63 & -9229.40\
10 & 7.88116 & -121.908 & -124.316 & -128.674 & 60 & 48.0432 & -8783.92 & -9130.29 & -9615.86\
11 & 8.72835 & -154.020 & -156.740 & -162.053 & 61 & 48.7880 & -9127.99 & -9496.65 & -10012.3\
12 & 9.56796 & -190.415 & -193.471 & -199.901 & 62 & 49.5310 & -9479.96 & -9871.77 & -10418.8\
13 & 10.3870 & -230.579 & -234.059 & -242.286 & 63 & 50.2713 & -9839.95 & -10255.2 & -10835.5\
14 & 11.1991 & -275.254 & -279.124 & -289.403 & 64 & 51.0151 & -10216.4 & -10649.8 & -11262.6\
15 & 12.0048 & -324.603 & -328.816 & -341.420 & 65 & 51.7609 & -10582.4 & -11055.1 & -11700.3\
16 & 12.8193 & -378.517 & -384.172 & -398.503 & 66 & 52.5067 & -10965.9 & -11470.4 & -12148.7\
17 & 13.6272 & -437.400 & -444.551 & -460.821 & 67 & 53.2510 & -11357.4 & -11895.2 & -12607.8\
18 & 14.4170 & -501.418 & -509.263 & -528.540 & 68 & 53.9927 & -11757.1 & -12329.1 & -13078.0\
19 & 15.2858 & -571.305 & -579.971 & -601.352 & 69 & 54.7306 & -12165.2 & -12771.4 & -13559.3\
20 & 16.1505 & -646.244 & -655.816 & -679.502 & 70 & 55.4635 & -12581.8 & -13221.8 & -14051.9\
21 & 16.9029 & -723.779 & -734.443 & -763.133 & 71 & 56.2925 & -13017.6 & -13695.5 & -14555.9\
22 & 17.6518 & -806.609 & -818.525 & -852.531 & 72 & 57.1215 & -13462.0 & -14179.9 & -15071.3\
23 & 18.3949 & -894.773 & -907.982 & -947.852 & 73 & 57.9500 & -13915.1 & -14674.8 & -15598.3\
24 & 19.1329 & -984.973 & -1002.97 & -1049.21 & 74 & 58.7777 & -14376.8 & -15180.3 & -16136.9\
25 & 19.8643 & -1087.71 & -1103.38 & -1156.87 & 75 & 59.6048 & -14847.3 & -15696.0 & -16687.4\
26 & 20.6041 & -1192.25 & -1211.07 & -1270.88 & 76 & 60.4345 & -15326.2 & -16224.1 & -17249.9\
27 & 21.3454 & -1302.72 & -1325.53 & -1391.42 & 77 & 61.2652 & -15813.9 & -16764.2 & -17824.6\
28 & 22.0817 & -1419.13 & -1446.14 & -1518.64 & 78 & 62.0954 & -16300.8 & -17315.4 & -18400.7\
29 & 22.8086 & -1536.57 & -1572.35 & -1652.71 & 79 & 62.9241 & -16806.4 & -17877.3 & -19011.3\
30 & 23.5219 & -1670.43 & -1703.53 & -1793.78 & 80 & 63.7501 & -17330.8 & -18449.2 & -19623.5\
31 & 24.3524 & -1809.22 & -1845.33 & -1941.63 & 81 & 64.6315 & -17861.7 & -19042.6 & -20248.3\
32 & 25.1804 & -1954.42 & -1993.69 & -2096.42 & 82 & 65.5128 & -18401.7 & -19647.9 & -20886.0\
33 & 26.0059 & -2106.13 & -2148.64 & -2258.28 & 83 & 66.3952 & -18950.8 & -20264.7 & -21536.7\
34 & 26.8353 & -2264.11 & -2311.42 & -2427.30 & 84 & 67.2796 & -19508.5 & -20894.9 & -22200.7\
35 & 27.6621 & -2428.77 & -2481.13 & -2603.59 & 85 & 68.1641 & -20075.4 & -21537.4 & -22878.2\
36 & 28.4813 & -2600.19 & -2656.87 & -2787.28 & 86 & 69.0467 & -20651.5 & -22191.1 & -23561.1\
37 & 29.3581 & -2780.21 & -2841.74 & -2978.07 & 87 & 69.9636 & -21241.1 & -22863.5 & -24237.8\
38 & 30.2331 & -2966.85 & -3033.62 & -3176.18 & 88 & 70.8808 & -21839.6 & -23548.8 & -24992.3\
39 & 31.0303 & -3155.02 & -3227.45 & -3381.68 & 89 & 71.6925 & -22437.9 & -24221.6 & -25724.9\
40 & 31.8263 & -3350.00 & -3428.61 & -3594.81 & 90 & 72.5050 & -23045.4 & -24907.7 & -26471.9\
41 & 32.6201 & -3546.33 & -3636.95 & -3815.67 & 91 & 73.3179 & -23639.5 & -25606.9 & -27233.7\
42 & 33.4115 & -3755.01 & -3852.55 & -4044.45 & 92 & 74.1309 & -24254.1 & -26319.2 & -28010.5\
43 & 34.2000 & -3976.23 & -4075.25 & -4281.19 & 93 & 74.9439 & -24878.0 & -27044.6 & -28802.9\
44 & 34.9923 & -4192.63 & -4306.84 & -4526.11 & 94 & 75.7568 & -25499.2 & -27783.2 & -29610.8\
45 & 35.7855 & -4422.14 & -4546.83 & -4779.23 & 95 & 76.5699 & -26141.7 & -28534.8 & -30434.9\
46 & 36.5766 & -4652.44 & -4794.58 & -5040.71 & 96 & 77.3858 & -26805.4 & -29302.3 & -31275.1\
47 & 37.3634 & -4902.97 & -5049.59 & -5310.66 & 97 & 78.2039 & -27466.1 & -30085.3 & -32132.1\
48 & 38.1440 & -5160.83 & -5311.30 & -5589.05 & 98 & 79.0232 & -28124.0 & -30883.5 & -33006.0\
49 & 38.9993 & -5424.43 & -5586.42 & -5875.84 & 99 & 79.8430 & -28803.8 & -31696.2 & -33897.2\
50 & 39.8532 & -5695.47 & -5869.68 & -6171.21 & 100 & 80.6627 & -29493.1 & -32523.3 & -34806.3\

[^1]: In this relation $\omega_{\mathrm{r}}$ and $m$ are measured in $\mathrm{cm}^{-1}$ and $\Omega_{0}$ in $\mathrm{cm}^{3}$.
