Ionization potentials in the limit of large atomic number

A central problem of electronic structure is the calculation of the ground-state energy of the electrons of any atom, molecule, or solid, within the non-relativistic Born-Oppenheimer limit. Density functional theory (DFT) is a popular choice, balancing computational efficiency with useful accuracy. The original DFT was that of Thomas Thomas and Fermi Fermi , TF theory, in which a local density approximation is made for the the kinetic energy and the electron-electron repulsion is approximated by the simple Coulomb energy of the charge density. In the 1970’s, Lieb and co-workers Lieb1 showed that the TF energy becomes relatively exact for neutral matter as $Z\rightarrow\infty$ in a specific way. The energy, $E$, grows in magnitude as $Z^{7/3}$, where $Z$ is the total charge. For atoms and their ions, the leading corrections in powers of $Z^{-1/3}$ were found by Scott S52 , Dirac Dir , Schwinger and others Sa80 ; ESa84 , as summarized by Englert En88 . These corrections are given exactly by extended Thomas-Fermi (ETF) theory, which includes both the gradient correction for the kinetic energy (one-ninth the von Weisacker functional vW ), and the local density approximation for exchange (LSDX KS ).

However, TF theory and its extensions are insufficiently accurate to predict chemical propertiesTeller . Modern DFT uses the Kohn-Sham (KS) scheme, in which only a very small fraction of the total energy, the exchange-correlation (XC), needs be approximated. But the idea of asymptotic correctness was recently extended to KS, relating the success of exchange GGA’s such as PBE PBE for total energies to their recovery of the $(Z^{1})$ term in the expansion of the exchange energy PCSB ; EB09 . The relation between semiclassical and local density approximations ELCB contributed to the creation of PBEsol PRCVSCZB , a PBE-like functional that is nearly optimum for solids near equilibrium but not for atoms and molecules, and to revTPSS PRCCS , a nearly optimal semi-local functional for all three kinds of systems.

But total electronic energies are irrelevant to chemistry. Only differences matter, such as the ionization potential of an atom ($I$ is the energy difference between the positive ion and the neutral) or the dissociation energy of a chemical bond. How relevant are asymptotic expansions for these quantities? The asymptotic expansion for $E$ is in powers of $Z^{-1/3}$, so if $I$ remains finite as $Z\to\infty$, the neutral and ion energies must agree for the first seven powers in such an expansion, a truly remarkable balancing act between quantum effects, the Pauli principle, and the Coulomb forces of nuclear attraction and inter-electron repulsion. In this letter, we demonstrate by both calculation and analysis that (i) $I$ has no limit as $Z\to\infty$, but remains column-dependent (ii) that each column has a finite limit; (iii) the local (spin) density approximation (LSD KS ) of KS theory becomes very accurate (if not exact) for $I$ for certain cases; (iv) ETF theory becomes very accurate (if not exact) for the average of $I$ over an entire shell; (v) the shell-averaged difference in density between the neutral and its ion approches that of TF. We demonstrate these statements in the limit of weak interelectron repulsion.

Our most important results are shown in Fig. 1. We plot $I$ from various calculations, extrapolated to infinite row number, versus the column number for main-group elements ($s$ and $p$ valence shells). We calculate exchange exactly KLI ; TS , using the optimized effective potential (OEP, which here should be indistinguishable from Hartree-Fock primer ), extrapolating all values to $Z\to\infty$. At the exchange level, LSD and PBE are almost exact for $p$-valence elements, and are highly accurate but inexact for the $s$-valence cases. Furthermore, ETF yields En88 a single number (3.15 eV), very close to the $s$- and $p$-average (3.02 eV). When correlation is included, gradient effects are slight, and it is in the regime of large electron number that approximate density functionals work best, sometimes exactly PCSB ; ELCB ; x . We speculate that LSDX on accurate densities becomes almost exact in this limit for $p$-shell cases, that ETF is exact for some shell-average, and that our XC results are extremely accurate and practically impossible to calculate with any other method.

To understand why local functionals become accurate in this limit, begin with total energies of neutral atoms, whose large-$Z$ expansion is

where $E_{Q}(Z)$ is the energy of an atom with atomic number $Z$ and charge $Q$, and the $c^{(j)}$ are coefficients depending on the degree of ionization, $q=Q/Z$. We use atomic units throughout. The neutral coefficients were derived via semiclassical analysis by Schwinger and Englert Sa80 ; ESa84 . The TF energy is exactly $-c^{(0)}_{q}\;Z^{7/3}$. The second term S52 comes from the $s$-orbitals at the nucleus and must be treated quantum mechanically. The third term is derivable in ETF theory Sa80 , of which 2/11 arises from the gradient correction to the kinetic energy, and 9/11 from LSDX. When treated in terms of the potential, the divergence at the nucleus is avoided En88 .

The extension of these ideas to $I$ has proven more difficult. Terms of higher order than those shown in Eq. (1) oscillate ESa84 with $Z$, as a precursor to the periodic variation of chemical properties that is missed by ETF, but well-described in KS DFT. The oscillations in $I$ dominate over trends with $Z^{-1/3}$. While numerous studies exist in the literature CGDPFF for fixed electron number $N$ with $Z\to\infty$, we are interested in $I(Z)=E_{1}(Z)-E_{0}(Z)$ as $Z\to\infty$. Within TF theory, Lieb proved Lieb2 that $I$ does not grow with $Z$, and by considering $c^{(0)}$ as $q\to 0$, Englert showed $I^{\rm TF}\to 3\Lambda^{-2/3}/7a\approx 1.29$ eV, where $\Lambda=32.729416$ is a known constant En88 , and $a=(9\pi^{2}/128)^{1/3}$. Even this simple result requires explanation, because $\mu$, the chemical potential, is zero for the neutral atom in TF theory, suggesting $I$ should be too. But the TF energy is the smooth envelope of $E_{Q}(Z)$ as a function of $q$, whereas the true energy consists of line segments between integer values PPLB . Thus $\mu=\partial E/\partial q=-I$ for the exact system, but the TF energy behaves as $q^{7/3}$ for small $q$. So $\mu_{TF}=0$, but the better value of $I^{TF}$ is the energy difference En88 with $Q=1$.

Next we discuss KS DFT, in which the (non-interacting) kinetic energy is not approximated, but is found exactly from the KS orbitals. We perform KS self-consistent calculations for atoms and ions up to 2938 electrons using LSD and PBE XC functional approximations, as well as the exact OEP exchange. These were done using the Engel code primer , but with tightened convergence criteria and maximum numbers of orbitals, and a logarithmic radial grid with 800 points. In Fig. 2 we show $I$ versus $Z^{-1/3}$ for each main-group column of the periodic table. In all cases, the behavior is almost linear as a function of $Z^{-1/3}$ for all $Z\gtrsim 169$, so we extrapolated these curves using a parabolic fit in $Z^{-1/3}$ and found the ionization energy for $Z\rightarrow\infty$ as shown in Fig. 1. The spherical approximations of the density (LSD, PBE) and of the potential (OEP), used in the Engel code (see Ref. TS ), give errors less that 0.1 eV for $I$. We use electronic configurations based on the aufbau principle and Madelung rule Made . For the noble gases, $Z={n}(n^{2}+6n+14)/6-\Delta(n)(n/2+1)$, where $n$ is the row number and $\Delta(n)=0$ for even and 1 for odd rows.

To understand in detail the results shown in Fig. 1, which are also tabulated in Table I, we begin at the exchange level. Both PBE and LSD exchange are almost identical to the OEP values for the $p$-group elements, with a maximum difference between them of 0.02 eV, and of either from OEP of 0.08 eV. This is not so for the alkalis and alkali earths, presumably because they have only one or two electrons outside a closed shell, with accompanying self-interaction error of approximate functionals. The ionization of $p$-elements involves removing electrons from a full (or almost full) shell with $\sim(n+2)^{2}/2$ electrons, where $n$ is the row number for even $n$, and rows $n$ and $n+1$ have the same structure.

In fact, Englert also showed that the TF result is not correct as $Z\to\infty$. The terms of $O(Z^{5/3})$ in Eq. (1) also yield a finite contribution, which is included in ETF, yielding an $I$ of 3.15 eV, very close to the average over both $s$ and $p$-shell values (3.02 eV).

To check this is no accident, consider the simpler system of atoms with an infinitesimal electron-electron repulsion, $\lambda$, sometimes called Bohr atoms. The orbitals are hydrogenic, requiring no self-consistency and simplifying the integrals HL . One finds that $I^{\rm TF}$ is exact for large $Z$ at $\lambda=0$. In Fig. 3, we show the exchange correction (divided by $\lambda$) to $I$ for LSDX applied to the TF density (yielding $8(2/3)^{1/3}/(3\pi^{2})\approx 0.2360$), to the exact densities (each averaged over entire shells), and exactly. All three match as $Z\to\infty$, but a small error remains if, e.g., just the $s$-shell is used. Thus we speculate that, for real atoms, LSDX (in a KS calculation) matches the average over the entire shell as $Z\to\infty$.

Next, we discuss the DFT calculations with correlation, which remains finite as $Z\to\infty$ and varies across a row. The differences between PBE and LSD are relatively small, giving greater confidence in both. The maximum deviation between them for $p$-elements is 0.14 eV, comparable to the deviations of these functionals at the exchange-only level from OEP for the alkali and alkali earths. Thus the gradient corrections are not vanishing, suggesting that while both calculations are accurate, neither is exact. The PBE average, 3.61 eV, is our best estimate of a universal ionization potential, defined as the limit of $I$ averaged over the $n$-th shell, as $n\to\infty$.

The other major descriptor of chemistry is the electron affinity $A(Z)=E_{0}(Z)-E_{-1}(Z)$. Within LSD or PBE, the first negative atomic ion of energy $E_{-1}(Z)$ has no stable solution, but $A(Z)$ can still be estimated P88 via a charged conductor model, in which $I-A=1/\langle r\rangle$, and $\langle r\rangle$ is the centroid of the added charge. Define the radial ionization density as

which integrates to $1$. Then choose $\langle r\rangle=\int^{\infty}_{0}\;dr\;r\Delta n_{R}(Z,r)$. Table I shows PBE $Z\rightarrow\infty$ limits for $I$, $\langle r\rangle$, and $A$. Averaging over $s$ and $p$, our best estimate for a universal value of $A$ is 0.78 eV.

We next extrapolate the ionization density via

which correctly integrates over $r$ to $1$. Here $Z_{0}=2935$, $\beta=1+b(Z^{-1/3}-Z^{-1/3}_{0})$, and $\gamma=c(Z^{-1/3}-Z^{-1/3}_{0})$, with fit parameters $b=5$, and $c=-2$. Finally, we also averaged over the 6 $p$-shell curves, to find the results shown in Fig. 4. The TF solution for the infinitesimally charged ion has a finite sizeEn88 ; Lieb2 , i.e.,

Beyond this radius, $\Delta n_{R}^{\rm TF}(r)$ is just the radial density of the neutral, which has reached its asymptotic form, decaying as $1/r^{4}$. The maximum of this curve is about 0.1830 at $r=8.855\;\rm{bohr}$. The agreement between the extrapolated $p$-shell densities and the TF theory is remarkably good, but not exact, while the extrapolated alkali ionization density is very different. We speculate that averaging over an entire shell would yield perfect agreement, as we find numerically for the Bohr atom.

Finally, we justify why such large atomic numbers (larger by a factor of 10 than those of Ref. CMA ) are needed to get these results. Because of the scaling with $Z^{-1/3}$, even $Z=125$ only makes $Z^{-1/3}=0.2$, while $Z>1000$ brings $Z^{-1/3}$ below 0.1, making the extrapolation much more reliable. In Fig. 5, we show accurate ionization densities for the 8-th column of the extended table at finite $Z$. The scaling of the TF ionization density is quite different from that of the exact solutions: Before extrapolation, even at $Z=2935$, the TF ionization density agrees much better with that of the alkalis, not the $p$-shell average. For the same reasons, having HF energies for only $Z\lesssim 100$, Englert erroneously concluded that $I^{\rm ETF}$ was the limit of the alkalis, not the shell-average (see Fig. 4-8 and its discussion of Ref. En88 ).

Thomas-Fermi theory produces the first term of Eq. (1) and extended TF yields an average $Z\to\infty$ limit for the ionization energy, but no periodic variation of chemical properties and no binding Teller of atoms to form molecules or solids. Within non-relativistic KS theory, any reasonable approximation to the XC energy with the correct uniform-density limit for exchange will produce the total-energy expansion of Eq. (1) and a finite column-dependent $Z\to\infty$ limit for the ionization energy. It appears that LSD is extremely accurate and possibly exact in certain cases for $I$. We have shown this for shell-averages in the limit of weak inter-electron repulsion. But in that case, the last shell is spread throughout the entire atom, and average gradients contributing to ionization vanish as $Z\to\infty$, which is not true for real atoms.

Thus we have established that, in the large-$Z$ limit, the periodic table becomes perfectly periodic. Moreover, local approximations appear to become exact, even for energy differences that are relatively vanishingly small in this limit. These are new, numerically relevant, exact conditions that approximate functionals should satisfy.

All conclusions are based upon numerical calculations and extrapolation. Proving them rigorously is a challenge to mathematical physics. This work was supported by DOE grant DE-FG02-08ER46496 at Irvine, and NSF (Grants DMR-0501588 and DMR-0854769) at Tulane.