Analytic model of a multi-electron atom

TL;DR

This paper introduces a fully analytical approximation for multi-electron atoms using a hydrogen-like basis with a single-effective charge, enabling accurate, convergent calculations of atomic properties including correlation effects.

Contribution

It develops a novel analytical perturbation method employing a hydrogen-like basis that captures correlation effects and allows closed-form summations over states.

Findings

Accurately describes the entire spectrum of multi-electron atoms.

Provides results comparable to multi-configuration Hartree-Fock at second order.

Achieves convergence and accuracy independent of the number of electrons.

Abstract

A fully analytical approximation for the observable characteristics of many-electron atoms is developed via a complete and orthonormal hydrogen-like basis with a single-effective charge parameter for all electrons of a given atom. The basis completeness allows us to employ the secondary-quantized representation for the construction of regular perturbation theory, which includes in a natural way correlation effects, converges fast and enables an effective calculation of the subsequent corrections. The hydrogen-like basis set provides a possibility to perform all summations over intermediate states in closed form, including both the discrete and continuous spectra. This is achieved with the help of the decomposition of the multi-particle Green function in a convolution of single-electronic Coulomb Green functions. We demonstrate that our fully analytical zeroth-order approximation describes the whole spectrum of the system, provides accuracy, which is independent of the number of electrons and is important for applications where the Thomas-Fermi model is still utilized. In addition already in second-order perturbation theory our results become comparable with those via a multi-configuration Hartree-Fock approach.

Tables

Table 1. Table 1: The comparison of the energy in a.u. of the second-order perturbation theory in the hydrogen-like basis Eq. ( 19 ) with the values via HF Saito ( 2009 ); da Silva and Trsic ( 1996 ); Froese ( 1967 ); *cohen10.1139.65.178 ; *doi:10.1139/p66-263 ; *TrefftzA1957oscillator , MCHF Froese ( 1967 ) , variational estimation and 1 / Z 1 𝑍 1/Z expansion ( H − superscript H \mathrm{H}^{-} , He, para- and ortho-He Zamastil  et al. ( 2010 ); Frolov ( 2015 ); Pekeris ( 1962 ) and Li Puchalski and Pachucki ( 2006 ) ).

	$E^{\left(2\right)}$	$E_{\mathrm{var}}$	$E_{\mathrm{MCHF}}$	$E_{\mathrm{HF}}$
${\text{H}}^{-}$	-0.532	-0.528	-0.528	-0.488
He	-2.907	-2.904	-2.903	-2.861
Li	-7.467	-7.478	-7.477	-7.433
He $2^{3}S$	-2.172	-2.175	-2.175	-2.174
He $2^{1}S$	-2.154	-2.146	-2.146	-2.143

Table 2. Table 2: An effective charge Z ∗ superscript 𝑍 Z^{*} and the comparison of the energy in a.u. of the zeroth-order approximation and the second-order perturbation theory in the hydrogen-like basis (single-particle excitation) Eq. ( 19 ) with the values via HF Saito ( 2009 ) . The line after Z = 20 𝑍 20 Z=20 indicates that for larger Z 𝑍 Z the inclusion of relativistic corrections is important Huang  et al. ( 1976 ) .

$Z$	$Z^{\ast }$	$E^{\left(0\right)}$	$E_{\mathrm{single}}^{\left(2\right)}$	$E_{\mathrm{HF}}$	$Z$	$Z^{\ast }$	$E^{\left(0\right)}$	$E_{\mathrm{single}}^{\left(2\right)}$	$E_{\mathrm{HF}}$
1	1.	-0.5	-0.5	-0.5	51	40.3872	-5974.	-6274.4	-6313.49
2	1.6875	-2.8477	-2.8610	-2.86168	52	41.2295	-6259.8	-6571.53	-6611.8
3	2.5454	-7.2891	-7.4114	-7.43273	53	42.0706	-6553.2	-6876.53	-6917.98
4	3.3716	-14.2096	-14.5212	-14.573	54	42.9104	-6854.3	-7189.47	-7232.1
5	4.1511	-23.6936	-24.4115	-24.5291	55	43.7925	-7165.6	-7510.2	-7553.9
6	4.9127	-36.2016	-37.4927	-37.6886	56	44.6732	-7484.4	-7838.7	-7883.5
7	5.6605	-52.0662	-54.1107	-54.4009	57	45.4977	-7804.6	-8174.35	-8221.1
8	6.3823	-71.2844	-74.3812	-74.8094	58	46.2332	-8125.8	-8516.62	-8566.9
9	7.0975	-94.4525	-98.8188	-99.4093	59	46.8783	-8447.6	-8865.53	-8921.2
10	7.8073	-121.908	-127.769	-128.547	60	47.6094	-8783.9	-9224.48	-9283.9
11	8.6561	-154.020	-160.894	-161.859	61	48.3384	-9127.99	-9591.84	-9655.1
12	9.4972	-190.415	-198.448	-199.615	62	49.0657	-9479.96	-9967.76	-10035.0
13	10.3161	-230.579	-240.453	-241.877	63	49.7914	-9839.95	-10352.3	-10423.5
14	11.1294	-275.254	-287.171	-288.854	64	50.6075	-10216.4	-10747.6	-10820.7
15	11.9377	-324.603	-338.769	-340.719	65	51.2340	-10582.4	-11146.7	-11226.6
16	12.7366	-378.517	-395.236	-397.505	66	51.9530	-10965.9	-11557.2	-11641.5
17	13.5314	-437.400	-456.884	-459.482	67	52.6702	-11357.4	-11976.5	-12065.3
18	14.3222	-501.418	-523.879	-526.818	68	53.3856	-11757.1	-12404.7	-12498.2
19	15.1910	-571.305	-595.918	-599.165	69	54.0996	-12165.2	-12841.9	-12940.2
20	16.0556	-646.244	-673.183	-676.758	70	54.8124	-12581.8	-13288.3	-13391.5
21	16.8063	-723.779	-755.341	-759.736	71	55.6210	-13017.6	-13746.6	-13851.8
22	17.5526	-806.609	-843.167	-848.406	72	56.4286	-13462.0	-14213.9	-14321.2
23	18.2939	-894.773	-936.754	-942.884	73	57.2350	-13915.1	-14690.2	-14799.8
24	18.9135	-984.973	-1035.65	-1043.36	74	58.0403	-14376.8	-15175.6	-15287.5
25	19.7636	-1087.71	-1141.81	-1149.87	75	58.8447	-14847.3	-15670.3	-15784.5
26	20.4882	-1192.25	-1253.26	-1262.44	76	59.6472	-15326.2	-16173.9	-16290.6
27	21.2099	-1302.72	-1371.06	-1381.41	77	60.4487	-15813.9	-16686.8	-16806.1
28	21.9279	-1419.13	-1495.28	-1506.87	78	61.1879	-16300.8	-17208.2	-17331.1
29	22.5146	-1536.57	-1625.17	-1638.96	79	61.9874	-16806.4	-17739.9	-17865.4
30	23.3548	-1670.43	-1763.60	-1777.85	80	62.8473	-17330.8	-18281.7	-18409.0
31	24.1826	-1809.22	-1908.47	-1923.26	81	63.7020	-17861.7	-18832.7	-18961.8
32	25.0083	-1954.42	-2060.00	-2075.36	82	64.5560	-18401.7	-19393.0	-19524.0
33	25.8319	-2106.13	-2218.27	-2234.24	83	65.4092	-18950.8	-19962.7	-20095.6
34	26.6516	-2264.11	-2383.19	-2399.87	84	66.2610	-19508.5	-20541.5	-20676.5
35	27.4694	-2428.77	-2555.03	-2572.44	85	67.1120	-20075.4	-21129.8	-21266.9
36	28.2853	-2600.19	-2733.88	-2752.05	86	67.9623	-20651.5	-21727.6	-21866.8
37	29.1585	-2780.21	-2919.52	-2938.36	87	68.8470	-21241.1	-22334.9	-22475.9
38	30.0296	-2966.85	-3111.99	-3131.55	88	69.7309	-21839.6	-22951.5	-23094.3
39	30.8213	-3155.02	-3310.59	-3331.68	89	70.5707	-22437.9	-23576.6	-23722.2
40	31.6110	-3350.00	-3516.37	-3539.00	90	71.4096	-23045.4	-24211.2	-24359.6
41	32.3199	-3546.33	-3728.64	-3753.6	91	72.1161	-23639.5	-24850.7	-25007.1
42	33.1052	-3755.01	-3949.05	-3975.55	92	72.8875	-24254.1	-25502.5	-25664.3
43	33.9674	-3976.23	-4177.47	-4204.79	93	73.6578	-24878.0	-26164.2	-26331.5
44	34.6664	-4192.63	-4411.72	-4441.54	94	74.3604	-25499.2	-26833.3	-27008.7
45	35.4442	-4422.14	-4654.37	-4685.88	95	75.1286	-26141.7	-27515.0	-27695.9
46	36.1379	-4652.44	-4903.88	-4937.92	96	75.9626	-26805.4	-28209.0	-28392.8
47	36.9945	-4902.97	-5162.74	-5197.7	97	76.7280	-27466.1	-28910.3	-29099.8
48	37.8493	-5160.83	-5429.17	-5465.13	98	77.4251	-28124.0	-29619.4	-29817.4
49	38.6966	-5424.43	-5703.19	-5740.17	99	78.1888	-28803.8	-30341.2	-30545.0
50	39.5426	-5695.47	-5984.91	-6022.93	100	78.9514	-29493.1	-31073.2	-31282.8

Equations

$$\mathsf{H}$$

$${\mathsf{H}}_0$$

$$\mathsf{W}$$

$$+\frac{1}{2}\sum _{\nu {\nu }_1\mu {\mu }_1}⟨\nu |⟨{\nu }_1|\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}|{\mu }_1⟩|\mu ⟩{\mathsf{a}}_{\nu }^{†}{\mathsf{a}}_{{\nu }_1}^{†}{\mathsf{a}}_{\mu }{\mathsf{a}}_{{\mu }_1}.$$

$$E(Z^{\ast })$$

$$=-Z^{\ast }(2Z-Z^{\ast })A+Z^{\ast }B,$$

$$Z^{\ast }$$

$$E^{(0)}$$

$$\Delta E^{(1)}$$

$$=-Z^{\ast }(Z-Z^{\ast })2A+Z^{\ast }B.$$

$$\Delta E^{(2)}=-\frac{1}{2}\sum _{k<l}\text{sideset}{}^{\prime }{\sum }_{{\sigma }_k{\sigma }_l}\frac{|W_{{\lambda }_k{\lambda }_l,{\sigma }_k{\sigma }_l}{|}^2}{E_{{\sigma }_k}+E_{{\sigma }_l}-E_{{\lambda }_k}-E_{{\lambda }_l}},$$

$$⟨{\lambda }_1\dots {\lambda }_N|$$

$$=⟨{\lambda }_k|U_{{\lambda }_k}|{\sigma }_k⟩-\sum _{l\ne k}⟨{\lambda }_l|V_{{\lambda }_k{\lambda }_l}^{(2)}|{\sigma }_k⟩,$$

$$G_{ZE}(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$G_{ZEl}(r,r^{\prime })$$

$${\tilde{G}}_{E_{{\lambda }_k}-\mathrm{i}\delta }^{{\lambda }_k}=G_{E_{{\lambda }_k}-\mathrm{i}\delta }-\sum _l{}^{\ast }\frac{|{\lambda }_l⟩⟨{\lambda }_l|}{E_{{\lambda }_l}-(E_{{\lambda }_k}-\mathrm{i}\delta )},$$

$$E_{\text{single}}^{(2)}$$

$$\Delta E_{\text{single},{\lambda }_k}^{(2)}=$$

$$+2\mathrm{Re}\sum _{r\ne k}⟨{\lambda }_k|U_{{\lambda }_k}{\tilde{G}}_{E_{{\lambda }_k}-\mathrm{i}\delta }^{{\lambda }_k}{V}_{{\lambda }_k{\lambda }_r}^{(2)\ast }|{\lambda }_r⟩$$

$$-\sum _{\begin{array}{c}l\ne k\\ r\ne k\end{array}}⟨{\lambda }_l|V_{{\lambda }_k{\lambda }_l}^{(2)}{\tilde{G}}_{E_{{\lambda }_k}-\mathrm{i}\delta }^{{\lambda }_k}{V}_{{\lambda }_k{\lambda }_r}^{(2)\ast }|{\lambda }_r⟩.$$

$$⟨{\lambda }_1\dots {\lambda }_N|$$

$$=⟨{\lambda }_k|⟨{\lambda }_l|V^{(2)}|{\sigma }_k⟩|{\sigma }_l⟩-⟨{\lambda }_l|⟨{\lambda }_k|V^{(2)}|{\sigma }_k⟩|{\sigma }_l⟩.$$

$${\int }_{-\infty }^{\infty }\frac{dt}{(t+a-\mathrm{i}\delta /2)(t-b+\mathrm{i}\delta /2)}$$

$$G_{E-\mathrm{i}\delta }(1,2)=-{\int }_{-\infty }^{\infty }\frac{dt}{2\pi \mathrm{i}}{G}_{t+\frac{E-\mathrm{i}\delta }{2}}\otimes G_{-t+\frac{E-\mathrm{i}\delta }{2}}.$$

$$\Delta E_{\mathrm{multi}}^{(2)}$$

$$\displaystyle\mspace{150.0mu}-\operatorname{Re}\int_{-\infty}^{\infty}\frac{dt}{2\pi\mathrm{i}}\langle\lambda_{k}|\langle\lambda_{l}|V^{(2)}\tilde{G}_{t+(E_{\lambda_{k}}+E_{\lambda_{l}}-\mathrm{i}\delta)/2}^{\lambda_{k}}\otimes\tilde{G}_{-t+(E_{\lambda_{k}}+E_{\lambda_{l}}-\mathrm{i}\delta)/2}^{\lambda_{l}}V^{(2)}|\lambda_{l}\rangle|\lambda_{k}\rangle\delta_{m_{s\lambda_{k}}m_{s\lambda_{l}}}\Bigg{)},$$

$$E^{(2)}=E^{(0)}+\Delta E_{\mathrm{single}}^{(2)}+\Delta E_{\mathrm{multi}}^{(2)}.$$

$$|{\psi }^{(1)}⟩$$

$$⟨{\psi }^{(1)}|\text{uprho}(Z^{\ast }\mathbf{r})|{\psi }^{(1)}⟩$$

$$\text{mspace}-90.0mu-2Z^{\ast 2}\mathrm{Re}\sum _{i=1}^N{\phi }_{{\lambda }_i}^{†}(Z^{\ast }\mathbf{r})\left[\int d{\mathbf{r}}^{\prime }{\tilde{G}}_{E_{{\lambda }_i}-\mathrm{i}\delta }^{{\lambda }_i}(Z^{\ast }\mathbf{r},r^{\prime })U_{{\lambda }_i}({\mathbf{r}}^{\prime }){\psi }_{{\lambda }_i}({\mathbf{r}}^{\prime })-\sum _{\begin{array}{c}l=1\\ l\ne i\end{array}}^N\int d{\mathbf{r}}^{\prime }{\tilde{G}}_{E_{{\lambda }_i}-\mathrm{i}\delta }^{{\lambda }_i}(Z^{\ast }\mathbf{r},r^{\prime })V_{{\lambda }_l{\lambda }_i}^{(2)}({\mathbf{r}}^{\prime }){\psi }_{{\lambda }_l}({\mathbf{r}}^{\prime })\right].$$

$${\rho }_{\mathrm{Ne}}^{(0)}=\frac{Z^{\ast 3}}{8\pi }{e}^{-2Z^{\ast }r}(16+e^{Z^{\ast }r}(2+Z^{\ast }r(Z^{\ast }r-2))),$$

$$f_{\text{sph.sim}}^{(0)}(q)$$

$$F_{nl}^{\text{sph.sim}}$$

$$J$$

$$K$$

$$\rho (\mathbf{r})$$

$${\rho }_{{\lambda }_k}(\mathbf{r})$$

$${\rho }_{{\lambda }_k}(\mathbf{r},{\mathbf{r}}^{\prime })$$

$${\psi }_{{\lambda }_k}(\mathbf{r})=R_{nl}(r)Y_{lm}(\Omega ).$$

Full text

basis set for chemistry, atomic perturbation theory, correlation effects, Coulomb Green function

The choice of the initial approximation for the single-electron wave functions (SEWF), plays an important role in modern quantum chemistry, both in the framework of the density functional theory or the solution of the Schrödinger equation. It is well known that the numerical solutions of the Hartree-Fock (HF) equations Hartree (1928a); *HartreeA1928wave-1; *FockA1930naeherungs provide the best possible zeroth-order approximation for SEWF. The inclusion of many HF configurations (MCHF) or application of various post-HF methods Dolg and Cao (2012); Bartlett and Musiał (2007); Bartlett and Stanton (1994); Puchalski et al. (2015); Filippin et al. (2016); Møller and Plesset (1934); Hirschfelder et al. (1964) allows one to include corrections due to correlation effects.

However, despite the great efficiency of modern numerical algorithms Grant (2007); Fischer (1977), simple analytical approximations Petersson et al. (1988); *doi:10.1063/1.460447; Jensen (2007); Szabo and Ostlund (1989); Roothaan (1951) still play an important role for many applications, where there is no need for extremely high accuracy, but a simple algorithm of repeated calculations of atomic characteristics is required. For example, the models based on, e.g., the Thomas-Fermi Thomas (1927); *FermiA1926medoto or multi-parametric screening hydrogen Hau-Riege (2012) approximations are widely used in computational plasma Ciricosta et al. (2012); Scott (2001); Chung et al. (2005); Lee (1987); Stewart and Kedar D. (1966) and X-ray physics Hau-Riege (2012); Toraya (2016), crystallography Feranchuk et al. (2002); Waasmaier and Kirfel (1995); Toraya (2016) or semiconductors physics Smith et al. (2014); Schulze-Halberg et al. (2013); Gaggero-Sager et al. (1998). In addition, the simplest possible inclusion of screening corrections in various cross sections like bremsstrahlung Seltzer and Berger (1985) or pair production Bethe and Maximon (1954); Davies et al. (1954) is required for later usage in particle-in-cell computer codes for simulation of strong laser-matter interaction Di Piazza et al. (2012), where computational efficiency is crucial.

In the present work we suggest a new basis set of fully analytical SEWF, which on the one hand provides a sufficiently accurate analytical zeroth-order approximation and on the other hand allows one to construct regular perturbation theory (RPT) for the inclusion of higher-order corrections. Our basis set includes the hydrogen-like wave functions with a single-variational parameter, namely the effective charge $Z^{*}$, which is identical for all SEWF of a given atom. The fact that the effective charge is identical for all SEWF is the principal difference of our approach in comparison with the inclusion of the multi-parametric screening corrections Slater (1930); Vainshtein and Safronova (1985); Feranchuk et al. (2002) or the quantum defect method Seaton (1958).

The identical effective charge for all wave functions automatically provides the complete and orthonormal basis and, consequently, renders the transition into the secondary-quantized representation natural. We have demonstrated that the analytical zeroth-order approximation contains the whole spectrum of a multi-electron atom and constructed a perturbation theory series, which converges fast with the rate $\sim 1/10$. In addition, we stress here that the accuracy of our results does not depend on the number of electrons in an atom, i.e., our approximation is uniformly available for all atoms or ions. Moreover, the results via second-order perturbation theory, are comparable with those via MCHF.

In addition, our approach is distinct from the one based on the application of the Coulomb-Sturmian basis set, which was successfully employed for the approximation of the SEWF in a variety of nonrelativistic Rotenberg (1962); Shull and Löwdin (1955) and relativistic scattering problems Fursa and Bray (2008); Bostock et al. (2014, 2011, 2013, 2012a, 2012b, 2012c); Bostock (2011). However, as highlighted in Ref. Rotenberg (1962) only the ground state wave function has a direct physical meaning. Consequently, it is problematic to interpret the occupation numbers for other non-ground states and to calculate the observable characteristics of multi-electron atoms such as their densities or form factors. Furthermore, the significant advantage of the hydrogen-like basis set is the knowledge of the closed-form expression through Whittaker functions for the Coulomb Green function, which in the Coulomb-Sturmian case is represented as a sum over Sturmian wave functions Shakeshaft (1985); Hostler (1970); Veselov and Labzowsky (1986). This analytical expression for the Coulomb Green function allowed us to perform all summations via intermediate states in perturbation theory including both the discrete and continuous spectra in closed form.

First of all, let us demonstrate the effectiveness of our basis for the calculation of the atomic ground-state energies of nonrelativistic atoms. For this purpose, we write down the Hamiltonian of an atomic system with a nucleus charge $Z$ and $N$ electrons in atomic units in the secondary-quantized representation Landau and Lifshitz (1977)

[TABLE]

Here the greek letters represent the collective quantum number $\nu=nlmm_{s}$ (or $\bm{k}lmm_{s}$ for the continuous spectrum) for the hydrogen-like wave function $\varphi_{\binom{nlm}{\bm{k}lm}}(Z^{*}\bm{r})\chi_{m_{s}}(s)=\langle\bm{r}|\binom{nlm}{\bm{k}lm}\rangle\langle s|m_{s}\rangle=\langle\bm{r}|\langle s|\nu\rangle$ with the effective charge $Z^{*}$ and $\langle\nu|\nu_{1}\rangle=\delta_{\nu\nu_{1}}$. The fermionic operators anticommute $\{\mathsf{a}_{\nu},\mathsf{a}^{{\dagger}}_{\nu^{\prime}}\}=\delta_{\nu\nu^{\prime}}$ and by acting on the $N$ particle state create the $N+1$ particle state $|\nu\lambda_{1}\ldots\lambda_{N}\rangle=\mathsf{a}^{{\dagger}}_{\nu}|\lambda_{1}\ldots\lambda_{N}\rangle$ 111We use the notations of Ref. Feynman (1998) for the secondary-quantized representation.. The Hamiltonian (1) is the exact expression written in the hydrogen-like basis, since we have only added and subtracted the term $Z^{*}/r$.

If the effective charge $Z^{*}$ is known, then the single-particle Hamiltonian of the zeroth-order approximation $\mathsf{H}_{0}$ is well defined. Consequently, the first question that we need to answer is how to calculate $Z^{*}$. For this we firstly performed the variational calculation, i.e., we choose the trial state vector $|\lambda_{1}\ldots\lambda_{N}\rangle$, which is characterized with a set of occupation numbers $g_{\lambda_{1}},\ldots,g_{\lambda_{N}}$, $g_{\lambda_{k}}=0$ or $1$ and calculated the energy of the system

[TABLE]

where $A=\sum_{k=1}^{N}g_{\lambda_{k}}/(2n_{k}^{2})$, $B=J+K$ is the sum of the Coulomb integral $J=1/2\int\rho(\bm{r})\rho(\bm{r}^{\prime})/|\bm{r}-\bm{r}^{\prime}|d\bm{r}d\bm{r}^{\prime}$ and the exchange integral $K=-1/2\int\sum_{k,l=1}^{N}g_{\lambda_{k}}g_{\lambda_{l}}\rho_{\lambda_{k}}(\bm{r},\bm{r}^{\prime})\rho_{\lambda_{l}}^{*}(\bm{r},\bm{r}^{\prime})/|\bm{r}-\bm{r}^{\prime}|d\bm{r}d\bm{r}^{\prime}$. In addition, $\rho(\bm{r})=\sum_{k=1}^{N}g_{\lambda_{k}}\rho_{\lambda_{k}}(\bm{r})$, $\rho_{\lambda_{k}}(\bm{r})=|\psi_{\lambda_{k}}(\bm{r})|^{2}$, $\rho_{\lambda_{k}}(\bm{r},\bm{r}^{\prime})=\psi_{\lambda_{k}}(\bm{r})\psi_{\lambda_{k}}^{*}(\bm{r}^{\prime})$ and $\psi_{\lambda_{k}}(\bm{r})$ is the hydrogen wave function, i.e., $Z=1$. The quantity $B$ can be calculated analytically for an arbitrary atom (see Appendix A).

Here we stress the extremely important fact that due to the property of the Coulomb potential the variable change $\tilde{\bm{r}}=Z^{*}\bm{r}$ leads to the simple energy dependence on the effective charge, which is given explicitly. Consequently, the analytical minimization of $E(Z^{*})$ with respect to the $Z^{*}$ yields the desired values of the effective charge and the energy of the system

[TABLE]

Proceeding further, let us construct the perturbation theory due to the operator $\mathsf{W}$. For this we utilize the eigenstates of the zeroth-order Hamiltonian $\mathsf{H}_{0}$ and calculate the first correction to the energy of the system. Since $\mathsf{H}_{0}$ is exactly diagonalizable, this is trivial to perform

[TABLE]

By observing Eq. (7) we can conclude that the above value of the effective charge $Z^{*}$ (5), found from the variational estimation leads to the vanishing first-order correction to the energy of the system $\Delta E^{(1)}=0$.

Before discussing the second-order correction let us briefly consider the accuracy of the zeroth-order approximation. Since the effective charge $Z^{*}$ is defined via Eq. (5), the energy of the system $E^{(0)}(Z,\{g_{\lambda}\})$ in Eq. (6) depends only on the charge of the nucleus and the set of occupation numbers $\{g_{\lambda}\}$, which should determine the minimal energy and satisfy the normalization condition $\sum_{k}g_{\lambda_{k}}=N$. We stress here that all our results, both for the zeroth-order approximation and for the second-order perturbation theory, are valid not only for atoms, but also for ions. Therefore, we write in the normalization condition $N$ and not $Z$, since in the general case $N\neq Z$.

Simple analytical calculations based on Eqs. (3) - (6) demonstrate that the optimal choice of the occupation numbers is given according to the “Aufbau” or Madelung-Janet-Klechkovskii rule Klechkovskii (1950); Madelung (1936); Wong (1979). For example, from the two sets of the occupation numbers $[\mathrm{Ar}]4\mathrm{s}^{1}$ and $[\mathrm{Ar}]3\mathrm{d}^{1}$ for potassium ($Z=19$), the former set possesses minimal energy (compare $-571.305$ versus $-568.473$). In Fig. 1 and Appendix C (Table LABEL:tab:2) the energy $E^{(0)}$ is compared with the results obtained via nonrelativistic HF equations Saito (2009). As can be concluded from Fig. 1 the chosen SEWF basis leads to a uniform approximation, i.e., it provides a relative accuracy of $\sim 5\%$ with respect to HF for all elements of the periodic table, which is considerably better than the Thomas-Fermi approximation Landau and Lifshitz (1977); Thomas (1927); *FermiA1926medoto.

Let us proceed with the calculation 222The details of all calculations are presented in supplementary information. of the correction to the energy of the system in second-order perturbation theory. Due to the two-particle structure of the perturbation potential $\mathsf{W}$ the only non-vanishing intermediate states in second-order are described by the state vectors $|\lambda_{1}\ldots\sigma_{k}\ldots\sigma_{l}\ldots\lambda_{N}\rangle$, in which $\lambda_{k}$ and $\lambda_{l}$ are replaced via intermediate states $\sigma_{k}$ and $\sigma_{l}$, respectively. Consequently, the second-order correction to the energy of the system reads

[TABLE]

where $W_{\lambda_{k}\lambda_{l},\sigma_{k}\sigma_{l}}=\langle\lambda_{1}\ldots\lambda_{N}|\mathsf{W}|\lambda_{1}\ldots\sigma_{k}\ldots\sigma_{l}\ldots\lambda_{N}\rangle$.

In this equation the sum is performed over all substitutions $\lambda_{k},\lambda_{l}$ with $\sigma_{k},\sigma_{l}$, $k,l=\{1\ldots N\}$ and the primed sum over $\sigma_{k},\sigma_{l}$ represents the sum over all possible quantum numbers excluding the ground state. It is convenient to split the total second-order correction to the energy as the sum of $\Delta E^{(2)}_{\text{single}}$, when only one electron goes into an intermediate state, and $\Delta E^{(2)}_{\text{multi}}$, when two electrons undergo the transition into intermediate states. In the first case the intermediate-state vector is $|\lambda_{1}\ldots\sigma_{k}\ldots\lambda_{N}\rangle$ with $\sigma_{k}\neq\lambda_{k}$ and in the second case $|\lambda_{1}\ldots\sigma_{k}\ldots\sigma_{l}\ldots\lambda_{N}\rangle$ with $\sigma_{k}\neq\lambda_{k}\lambda_{l}$ and $\sigma_{l}\neq\lambda_{k}\lambda_{l}$.

We continue the calculation of the single-electron excitation. The required matrix elements can be easily evaluated and represented in compact form

[TABLE]

where $U_{\lambda_{k}}=V^{(1)}+\sum_{l\neq k}V^{(2)}_{\lambda_{l}\lambda_{l}}$, $V^{(1)}=-Z^{*}(Z-Z^{*})/r$, $V^{(2)}=Z^{*}/|\bm{r}-\bm{r}^{\prime}|$ and $V^{(2)}_{\lambda_{k}\lambda_{l}}=\langle\lambda_{k}|V^{(2)}|\lambda_{l}\rangle$. Here we also carried out the variable change $\bm{r}\to Z^{*}\bm{r}$ in order to separate out the explicit dependence on $Z^{*}$.

Since the required matrix element has been evaluated, we have to perform the summation over intermediate states. Here we face one of the largest advantages of using a hydrogen-like basis set as the required summation can be performed in closed form Dalgarno and Lewis (1955); Hostler (1969); Dalgarno (1959), since the single-particle Coulomb Green function is well known Johnson and Hirschfelder (1979); Hostler and Pratt (1963); Hostler (1964); Veselov and Labzowsky (1986) and is expressed through the decomposition over spherical harmonics $Y_{lm}(\Omega)$ and Whittaker functions Whittaker and Watson (1996) $W_{\kappa,\mu}(r,r^{\prime}),M_{\kappa,\mu}(r,r^{\prime})$

[TABLE]

where $\nu=\frac{Z}{\sqrt{-2E}}$, $r_{>}=\max(r,r^{\prime})$ and $r_{<}=\min(r,r^{\prime})$.

In order to calculate $\Delta E^{(2)}_{\text{single}}$ one needs to take into account that in Eq. (8) not all states are present. Moreover, since we are dealing with a multi-electron system we have to take into account the Pauli exclusion principle. It reveals itself here in the subtractions of the occupied states from the Green function of the electrons with the same spin as the electron that undergoes the transition into intermediate states. Consequently, we introduce the reduced Coulomb Green function

[TABLE]

where $\lambda_{k}$ in $\tilde{G}^{\lambda_{k}}_{E_{\lambda_{k}}-\mathrm{i}\delta}$ and the starred sum describe the subtractions with the same spin as the electron $\lambda_{k}$. With the help of Eqs. (8)–(12) we obtain the energy of the system including the single-electron excitation correction

[TABLE]

where

[TABLE]

Here we also pay attention to the fact that both the Green function and the states $|\lambda_{k,l}\rangle$ in Eq. (14) are referred to the hydrogen wave functions, since due to the variable change $\bm{r}\to Z^{*}\bm{r}^{\prime}$ the dependence on $Z^{*}$ is given explicitly and the $Z^{*}$ in the denominator of Eq. (10) is cancelled with the one from the matrix element.

The results obtained via Eqs. (13), (14) are presented in Fig. 1 and in Appendix C (Table LABEL:tab:2) for nuclear charges $Z=1..100$, however, for $Z$ larger than $20$ the incorporation of relativistic corrections would be required for precise values. One can observe that the inclusion of the single-electron excitation reduces the relative error by one order of magnitude. The inclusion of $\Delta E^{(2)}_{\text{single}}$ does not take into account correlation effects and consequently the corrected wave functions still remain in the class of SEWF. For this reason the condition $|E^{(2)}_{\text{single}}|<|E_{\mathrm{HF}}|$ holds as the solution of HF equations realizes a minimum of the functional. In addition, single-electron correction in third-order perturbation theory should also not be taken into account, since its value is smaller than the corresponding second-order correlation correction.

Let us proceed with the calculation of the correlation correction $\Delta E^{(2)}_{\mathrm{multi}}$. The required matrix elements in this case are represented as

[TABLE]

Here only the two-particle part of the operator $\mathsf{W}$ contributes, since both $\sigma_{k}$ and $\sigma_{l}$ can not be equal to both $\lambda_{k}$ and $\lambda_{l}$, respectively.

In order to perform a summation over intermediate states in this case we need to obtain a two-particle Coulomb Green function, which is not known. However, with the help of the identity

[TABLE]

the two-particle Green function can be represented as a convolution of two single-particle ones

[TABLE]

The only nontrivial operation remained is to perform the required subtractions, taking into account the Pauli exclusion principle. Let us illustrate this in the lithium case. For example, for the matrix element $\langle\lambda_{1}\lambda_{2}\lambda_{3}|\mathsf{W}|\lambda_{1}\sigma_{1}\sigma_{2}\rangle$, $\sigma_{1}$ can not be equal to $\lambda_{1}\lambda_{2}\lambda_{3}$. The same applies to $\sigma_{2}$ 333The situation when $\sigma_{2}=\lambda_{3}$ and $\sigma_{1}\neq\lambda_{1}\lambda_{2}\lambda_{3}$ or equivalently $\sigma_{1}=\lambda_{2}$ and $\sigma_{2}\neq\lambda_{1}\lambda_{2}\lambda_{3}$ is taken into account in the single-electron excitation.. Consequently, taking into account the spin orthogonality, we need to subtract $\frac{|\lambda_{2}\rangle\langle\lambda_{2}|}{(t+E_{0}/2-\mathrm{i}\delta/2)-E_{\lambda_{2}}}$ from the Green function with the index $t$, with $E_{0}=E_{\lambda_{2}}+E_{\lambda_{3}}$. The Green function with the index $-t$ undergoes two subtractions, namely $\frac{|\lambda_{1}\rangle\langle\lambda_{1}|}{(-t+E_{0}/2-\mathrm{i}\delta/2)-E_{\lambda_{1}}}$ and $\frac{|\lambda_{3}\rangle\langle\lambda_{3}|}{(-t+E_{0}/2-\mathrm{i}\delta/2)-E_{\lambda_{3}}}$.

With the help of the above notation for the reduced Green function Eq. (12), the correlation correction is written as

[TABLE]

which is valid for an arbitrary atom or ion due to the pairwise character of the correlation contribution.

Finally, combining all together we obtain the total energy of the system in second-order perturbation theory

[TABLE]

In order to demonstrate the effectiveness of our basis set we have evaluated the energy Eq. (19) in second-order perturbation theory for the ground states of $\mathrm{H}^{-}$, He and Li and the excited states for He, namely ortho- and para-helium (see Table. 1). We note here, that due to the degeneracy of ortho- and para-helium the perturbation theory should be modified, i.e., the zeroth-order state is defined as Feranchuk and Triguk (2011) $|\psi^{2\,^{3}S,2\,^{1}S}\rangle=\frac{1}{\sqrt{2}}(|\lambda_{1}\uparrow,\lambda_{2}\downarrow\rangle\pm|\lambda_{1}\downarrow,\lambda_{2}\uparrow\rangle)$, the so called perturbation theory for the doubly degenerate energy levels Landau and Lifshitz (1977) (See also supplementary material). In addition, for some energy levels in Table 1 ($\mathrm{H}^{-}$, He, He $2\,^{1}\mathrm{S}$) their energy values within a complete second-order perturbation theory are smaller than the corresponding exact results via the variational method. This is related to the fact that in our calculations we employ perturbation theory series, which convergence to the exact value can be oscillatory in some problems Feranchuk et al. (2015, 1995); An et al. (1986), i.e., the absolute value of the difference between the exact and the approximate results is decreasing in each order of perturbation theory, however, in the second order the approximate value is smaller than the exact result, while in the third order it is larger respectively. The mathematical proof of this convergence property requires additional investigations.

At last we want to demonstrate that our basis set provides a good approximation not only for the integral characteristics of the system but also for the local ones. For this we have evaluated the radial electron density $4\pi r^{2}\rho(r)$, with the wave function

[TABLE]

which includes the first-order single-electron excitation correction over $\mathsf{W}$. Consequently, one finds the expectation value of the density operator $\langle\psi^{(1)}|\uprho(\bm{r})|\psi^{(1)}\rangle=\langle\psi^{(1)}|\sum_{\nu\nu^{\prime}}\psi^{{\dagger}}_{\nu}(\bm{r})\psi_{\nu^{\prime}}(\bm{r})\mathsf{a}^{{\dagger}}_{\nu}\mathsf{a}_{\nu^{\prime}}|\psi^{(1)}\rangle$ up to first-order in $\mathsf{W}$:

[TABLE]

The dependence of the density on the radial variable $r$ for Ne and Ar, which possess spherically symmetric radial density, is presented in Fig. 2. We immediately observe that already in the fully analytical zeroth-order approximation the error in the density does not exceed $\sim 20\%$ in comparison with the corresponding HF value Wang and Parr (1977). Moreover, our fully analytical result provides much better agreement than the quasi-classical Thomas-Fermi model Landau and Lifshitz (1977); Thomas (1927); *FermiA1926medoto. At the same time, the inclusion of the single-electron excitation correction improves the agreement with HF significantly.

As was mentioned in the introduction, when a large number of repeated calculations need to be performed, the simplest possible expressions for the electronic density and the spectrum of the system are required. This happens in, e.g., computer codes for plasma simulations like CRETIN Scott (2001), FLYCHK Chung et al. (2005), LASNEX-DCA Lee (1987) where Thomas-Fermi model is used for the determination of the electrostatic potential Stewart and Kedar D. (1966) or for the calculation of X-ray scattering factors Hau-Riege (2012) in crystallography, where fits of an electronic density are employed. Since our zeroth-order approximation is fully analytical it can be perfectly used for these purposes. For example, for the density of Ne, the atom with the largest number of electrons in the first row of the periodic table, one obtains in the zeroth-order approximation

[TABLE]

which is extremely simple and provides better qualitative agreement with HF than the Thomas-Fermi model, see Fig. 2. For larger atoms, expressions for the density look similar and have a structure of an exponential function which multiplies a polynomial, as can be found in Appendix B.

Moreover, the Fourier transform $f^{(0)}(\bm{q})=\int\rho^{(0)}(\bm{r})e^{\mathrm{i}\bm{q}\cdot\bm{r}}d\bm{r}$ of an electronic density for an arbitrary atom or ion can also be calculated analytically in the zeroth-order approximation (see Appendix. B). However, specifically for the atoms with spherically symmetric density the following closed form expression can be employed Benediktovitch et al. (2013)

[TABLE]

Usually for numerical applications the fit of the form $f(s)=\sum_{i}a_{i}\exp{(-b_{i}s^{2})}$ is used for the atomic scattering factors Waasmaier and Kirfel (1995); Rez et al. (1994); Prince (2006). Here $s=\sin\theta/\lambda\,[\text{\AA}^{-1}]$, $\theta$ is the scattering angle and $\lambda$ the wavelength of the X-ray radiation. The parameter $s$ is related to $q$ as $4\pi s\cdot 0.529177=q$. Consequently, in Fig. 3 we compare highly accurate Gaussian fits Waasmaier and Kirfel (1995) of results from relativistic HF calculations with our simple analytical expressions. For light elements we find good agreement; however, for larger elements the discrepancy increases, indicating the importance of corrections.

Lastly, we conclude that the proposed procedure can be straightforwardly generalized for the relativistic hydrogen-like basis set, since the analytical form of the relativistic Coulomb Green function is known Hylton (1984); Swainson and Drake (1991a); *0305-4470-24-8-022. The analytical zeroth-order wave functions can substitute the less accurate Thomas-Fermi approximation. We also remind that our calculation scheme is valid for ions, since the number of electron is not necessary equal to $Z$, which can be useful for the statistical theory of plasma. Our approach can be easily modified to include interactions with external fields, employed for the classification of excited states and applied in time dependent problems, where the knowledge of the system’s spectrum is required. In addition our results provide a uniform approximation for the observable characteristics, i.e., independent of the number electrons in an atom. Furthermore, since our zeroth-order approximation yields algebraic expressions for electronic densities and scattering factors, our results can be useful for repeated calculations in other particle-in-cell and plasma computer codes for simulation laser-matter interactions.

Finally, the authors are working toward the release of the computer software, which will incorporate relativistic corrections, will allow an automated calculation of the energies of the excited states, transition matrix elements and oscillator strengths.

Appendix A Explicit expressions for the zeroth-order energy

In this Appendix we present an analytical calculation of the zeroth-order approximation for the energy of the system.

According to the discussion before Eq. (4), in order to calculate an effective charge $Z^{*}$ one needs to evaluate the expectation value of the Hamiltonian with a trial state vector $|\lambda_{1}\ldots\lambda_{N}\rangle$. The evaluation of the matrix elements is presented in Ref. Szabo and Ostlund (1989). As a result, the calculation of the energy is reduced to the calculation of the Coulomb and exchange integrals

[TABLE]

where we have introduced the diagonal and of-diagonal elements of the density matrix

[TABLE]

The hydrogen wave function $\psi_{\lambda_{k}}(\bm{r})$ in the expressions (26)–(28) is the product of the spherical harmonic $Y_{lm}(\Omega)$ and the radial wave function

[TABLE]

The actual expression for $R_{nl}(r)$ is given, for example, in Landau and Lifshitz (1977).

In the following we will use the expansion of the Coulomb interaction between electrons over spherical harmonics Landau and Lifshitz (1977)

[TABLE]

with $r_{>}=\max(r,r^{\prime})$ and $r_{<}=\min(r,r^{\prime})$. In addition, the integration of the product of three spherical harmonics yields 3j symbols Cowan (1981)

[TABLE]

Consequently, the calculation of the quantities $J$ and $K$ reduces to the computation of the two integrals and summation over occupation numbers.

We proceed with the evaluation of the Coulomb part. For this, we firstly calculate

[TABLE]

By plugging Eq. (29) into Eq. (32), using expansion of Eq. (30) and integrating out the angular variables one obtains

[TABLE]

where

[TABLE]

and

[TABLE]

In addition the quantum numbers $\lambda_{k}=nlmm_{s}$ and $\lambda_{l}=n_{1}l_{1}m_{1}m_{s1}$. Therefore, the Coulomb integral is equal to

[TABLE]

which can be fast evaluated using computer algebra software of a choice, e.g., MATHEMATICA. We would like to mention here, that the same set of integrals appears constantly and in order to speed up the evaluation it makes sense to precalculate the integrals and store the values in an array as a function of quantum numbers. Consequently, once this is performed, the evaluation of the Coulomb integral for any set of quantum numbers can be done almost instantly.

In a full analogy one can calculate the exchange integral, yielding

[TABLE]

where

[TABLE]

and

[TABLE]

In addition, the above discussion about calculation efficiency of the Coulomb integrals is fully applicable for the exchange integral.

Appendix B Explicit expressions for the zeroth-order electronic density

In this appendix we present the explicit expressions for the electronic density in the analytical zeroth-order approximation together with its Fourier transforms for a number of selected atoms. In the following we use the notation $Z^{*}r=u$.

He

[TABLE]

C

[TABLE]

O

[TABLE]

Ne

[TABLE]

Xe

[TABLE]

while the Fourier transform can be obtained by applying the expression (23).

For all other atoms or ions similar expressions can be obtained, yielding the product of an exponential by a polynomial functions.

Appendix C Calculated values of the effective charges, ground state energies and their comparison with HF