Exact density functionals for the exchange and correlation energies are approximated in practical calculations for the ground-state electronic structure of a many-electron system. An important exact constraint for the construction of approximations is to recover the correct non-relativistic large-Z expansions for the corresponding energies of neutral atoms with atomic number Z and electron number N = Z, which are correct to the leading order (−0.221Z5/3 and −0.021Z ln Z, respectively) even in the lowest-rung or local density approximation. We find that hydrogenic densities lead to Ex(N, Z) ≈ −0.354N2/3Z (as known before only for ZN ≫ 1) and Ec ≈ −0.02N ln N. These asymptotic estimates are most correct for atomic ions with large N and ZN, but we find that they are qualitatively and semi-quantitatively correct even for small N and NZ. The large-N asymptotic behavior of the energy is pre-figured in small-N atoms and atomic ions, supporting the argument that widely predictive approximate density functionals should be designed to recover the correct asymptotics. It is shown that the exact Kohn–Sham correlation energy, when calculated from the pure ground-state wavefunction, should have no contribution proportional to Z in the Z limit for any fixed N.

In this work, we will find closed-form formulas for the exchange energy and correlation energy of an atom or atomic ion with electron number N and proton number Z. We will paint with a broad brush, seeking not the most accurate formulas but the simplest and most understandable ones, from which we can draw conclusions relevant to the construction of density functional approximations for these energies.

In exact non-relativistic quantum chemistry,1 the Hartree–Fock ground-state wavefunction is that single Slater determinant that minimizes the expectation value of the Hamiltonian. The quantum chemical correlation energy is the difference between the true total energy and the Hartree–Fock total energy. In exact Hohenberg–Kohn–Sham density functional theory,2–6 the Kohn–Sham ground-state wavefunction is that wavefunction that yields the true ground-state electron density and minimizes the expectation value of the kinetic energy,7 i.e., it is the ground eigenstate of the Kohn–Sham effective Hamiltonian. When that ground-state is degenerate, it can be a linear combination of a few Slater determinants, chosen to connect adiabatically7 to a given interacting ground state. When it is a single Slater determinant, the exact Kohn–Sham exchange and correlation energies of atoms and atomic ions are numerically close to those defined in quantum chemistry, and the quantum chemical correlation energy is an upper bound to the exact Kohn–Sham correlation energy.8 Only for one-electron densities do these two exact theories have exactly the same Slater determinant, the same exchange energy (to cancel the Hartree energy), and the same (zero) correlation energy.

In Kohn–Sham theory, the exchange and correlation energies are functionals of the electron density. Approximations to these functionals are made for the sake of practical computation for real atoms, molecules, and solids. The simplest such approximation is the local density approximation (LDA),3 

EνLDA[n]=n(r)ενunif(n(r))dr,
(1)

where n(r) is the electron density, ν = x (exchange) or c (correlation), and ενunif(n) is the corresponding energy per electron in an electron gas of uniform density n(r). [A spin-polarized system requires ενunif(n,n).] Higher-rung functionals (e.g., Refs. 9 and 10) retain the correct uniform-gas limit while satisfying other exact constraints.

Although the LDA is exact for a density that varies slowly over space, its relevance to real atoms and molecules is not obvious. Dirac11 added LDA exchange to the Thomas–Fermi model, and Schwinger12 may have been the first to realize that LDA exchange becomes relatively exact for neutral atoms in the limit of large atomic numbers. In this limit, the bulk of the density becomes Thomas–Fermi-like, with a locally slow spatial variation,13,14 and the exact energies have large-Z asymptotic expansions,15–18 

Ex(Z,Z)=AxZ5/3+BxZ+(N=Z),
(2)
Ec(Z,Z)=AcZlnZ+BcZ+(N=Z).
(3)

The leading coefficients (Ax = 0.221 hartree and Ac = 0.021 hartree) are known17–19 to be those from the LDA evaluated on the self-consistent Kohn–Sham (or Thomas–Fermi) density for large N = Z, and corrections to the LDA (e.g., Refs. 9 and 10) can give higher-order coefficients. Recent work20–22 has shown that errors of even a few percent in the uniform-density limit can seriously undermine the accuracy of approximate density functionals for the equilibrium properties of atoms and small molecules.

There is evidence that satisfaction of Eqs. (2) and (3) can produce functionals that are notably accurate for the atomization energies of molecules, without being fitted to molecules. The Becke 1988 (B8823) generalized gradient approximation (GGA) for exchange, still widely used in chemistry, was constructed to recover the LDA in the slowly varying limit [recovering Ax of Eq. (2)] and was fitted only to the exchange energies of rare-gas atoms but was noted to be consistent with the leading asymptotic correction to LDA exchange. Elliott and Burke16 showed that this functional recovers a nearly exact (within 1.1%) value for Bx of Eq. (2). The Perdew–Burke–Ernzerhof (PBE9) GGA, while it is more accurate than B88 for solids, is less accurate for molecular atomization energies and for Bx. However, the acGGA18 revision of PBE and the SCAN meta-GGA,10 both constructed to satisfy Eqs. (2) and (3) (as well as other exact constraints) without fitting to molecules, are more accurate for atomization energies than is PBE (and SCAN is remarkably more accurate).

Somewhat related to the large-Z expansion of Eqs. (2) and (3) is the 1/Z perturbation expansion24,25 of the total energy of an N-electron atomic ion,

E(N,Z)=Z2ε1(N)+ε2(N)Z+,
(4)

where the unperturbed problem is N non-interacting electrons in the potential −Z/r and the perturbation is the Coulomb repulsion among the electrons. Only a few of the coefficients in Eq. (4) are known and only for a few small values of N. In the limit NZ, the leading terms of Eq. (4) should be relatively accurate. Moreover, the first term exactly describes N = 1. In the regime NZ, Z is large enough that the electron–nucleus attraction, −Z/r, dominates electron–electron Coulomb repulsion, and the latter is treated as a perturbation. The unperturbed problem (i.e., the non-interacting system) has a hydrogenic density that occupies the N scaled [Z3/2ψi,σ(Zr)] hydrogen-atom spin orbitals ψi,σ(r) of lowest energy. Clearly, for ZN − 1, the −Z/r interaction of the electron with the nucleus will dominate all other terms in the Kohn–Sham one-electron potential. For Z = N, however, the exact Kohn–Sham potential will vary from −Z/r (plus a positive constant from the Hartree and exchange–correlation potentials) to −1/r as r increases from 0 to . Moreover, by the Hellmann–Feynman theorem, Z∂E(N, Z)/∂Z = 2Z2ε1(N) + 2(N) + ⋯ is the electron-nuclear attraction energy.

The various components of the total energy also have 1/Z expansions. For example,

Ex(N,Z)=α1(N)Z+α2(N)+,
(5)
Ec(N,Z)=β1(N)Z+β2(N)+.
(6)

The Kohn–Sham wavefunction was discussed in the second paragraph of this article. Here, the full Hamiltonian Ĥ is the sum of a hydrogenic part Ĥ0 and a weak electron–electron perturbation V^ee. In a quantum chemical calculation, β1(N) = 0 for those electron numbers N (e.g., N = 1, 2, 3, 7, 8, 9, 10, 11) for which the ground state of Ĥ0 is either non-degenerate (N = 2, 10) or is without degenerate configurations of the same multiplet symmetry so that the Kohn–Sham wavefunction is a single Slater determinant. These must also be electron numbers for which β1(N) = 0 in exact Kohn–Sham theory, since there is no qualitative difference between the Kohn–Sham and Hartree–Fock wavefunctions. For N = 4, where the degenerate configurations (1s)2(2s)2 and (1s)2(2p)2 can both belong to the multiplet 1S the Kohn–Sham wavefunction is a linear combination of Slater determinants, differing significantly from the Hartree–Fock wavefunction. A different choice of Kohn–Sham representation (using an ensemble ground-state density rather than a pure ground-state wavefunction as in  Appendix A) for N = 4 and ZN in which the 2s and 2p orbitals are degenerate for Z > 23 is discussed in Ref. 26. The limiting constant values for the correlation energy can arise27 because the hydrogenic density Z3fN(Zr) is uniformly scaled to the high-density limit when Z at fixed N, and this uniform-scaling behavior is built into the PBE GGA9 and the SCAN meta-GGA.10 We argue in  Appendix A that, within the exact Kohn–Sham theory, β1(N) = 0 regardless of degeneracies. Our proof does not hold for the perturbation series of the quantum chemical correlation energy for which it is well-known that the leading-order is Z0 for a non-degenerate ground state and Z1 for some degenerate ground states.28,29 The exact Hartree–Fock and exact Kohn–Sham exchange energies have different uniform scaling behaviors,7 and their wavefunctions differ substantially in the case of a degenerate ground-state. Thus we can expect different large-Z limits of their corresponding correlation energies.

The quantum chemical coefficients α1(N) and β2(N) for small N are reported in Tables I and II of Ref. 30. Table I shows that α1(N)/N2/3 and β2(N)/(N ln N) are nearly independent of N (not noticed in Ref. 30), suggesting that α1(N) ∼ N2/3 and β2(N) ∼ (N ln N). Here, we use ∼ to denote the leading-order behavior of a function. Thus, at least for NZ and for N such that β1(N) = 0,

Ex(N,Z)0.354N2/3Z,
(7)
Ec(N,Z)0.02NlnN.
(8)

The value −0.354 is the analytic large-Z limit of α1(N)/N2/3, as explained around Eq. (10). The value −0.02 is a roundoff of all the values of β2(N)/(N ln N) for N > 2 from Table I, accounting for the larger uncertainty in the numeric values of this coefficient. Now, setting N = Z in Eqs. (7) and (8) leads to a result that is qualitatively and semi-quantitatively like the leading terms of Eqs. (2) and (3). In particular, Ex(Z, Z) ∼ Z5/3 and Ec(Z, Z) ∼ Z ln Z. The primary difference is in the coefficient of Z5/3, in part, because the density functional for the exchange energy depends on the detailed shape of the electron density (e.g., hydrogenic vs self-consistent neutral Thomas–Fermi) but not so for correlation in the leading order.  Appendix B presents a simple derivation of the leading-order terms in the asymptotic series for the correlation energy and shows that they are identical for hydrogenic and self-consistent neutral Thomas–Fermi densities.

TABLE I.

The leading coefficients α1(N) and β2(N) (both in hartree) in Eqs. (5) and (6), from Ref. 30, divided by the displayed functions of electron number N.

Nα1(N)/N2/3β2(N)/(N ln N)
−0.3937 −0.0337 
−0.3471 −0.0163 
−0.3508 −0.0174 
−0.3569 −0.0184 
−0.3658 −0.0187 
10 −0.3766 −0.0186 
11 −0.3647 −0.0173 
Nα1(N)/N2/3β2(N)/(N ln N)
−0.3937 −0.0337 
−0.3471 −0.0163 
−0.3508 −0.0174 
−0.3569 −0.0184 
−0.3658 −0.0187 
10 −0.3766 −0.0186 
11 −0.3647 −0.0173 

The more important conclusion, which will be validated and applied in the rest of this article, is that the large-N asymptotics of the exchange and correlation energies for ZN are discernible even for small N. A corollary to this is that the large-N asymptotics can be roughly estimated from the small-N energetics.

Conlon31 showed that, in the limit Z = N, the Hartree–Fock exchange energy for any Coulomb system tends to the LDA exchange energy evaluated on the self-consistent Thomas–Fermi density. Thus, there is no inherent contradiction in the values of the exchange coefficient: hydrogenic densities are not the correct ZN limits. The Thomas–Fermi approximation to the hydrogenic density, nhTF, is known analytically [refer to Eq. (B3) in  Appendix B or Refs. 17, 32, and 33]. Evaluating the LDA on the Thomas–Fermi approximation for the hydrogenic density of a neutral atom yields, for N = Z,

ExLDA[nhTF]=231/34π2Z5/30.354Z5/3.
(9)

For heavy positive ions, where N < Z, the LDA exchange energy evaluated on the hydrogenic Thomas–Fermi density tends to, for NZ (see Refs.32 and 34),

ExLDA[nhTF]=0.354N2/3Z1+ONZ.
(10)

A numeric estimate of this coefficient from Ref. 35 on a neutral hydrogenic density of 2030 electrons agreed precisely with the analytic values of Eqs. (9) and (10).

The exact exchange energy for a given spin-unpolarized density n(r) is expected to be bounded by the conjectured tight lower bound,36,37

Ex[n]1.174ExLDA[n],
(11)

satisfied for all spin-unpolarized densities by LDA and SCAN,10 but not by PBE.9 This bound holds rigorously for a spin-unpolarized two-electron density.36 The LDA is expected to be relatively less accurate for such densities than for spin-unpolarized densities with N > 2, and no violation of the conjectured bound is known. The LDA typically becomes relatively exact as more electrons are packed into a given volume of space. We give an alternative derivation of the exactness of LDA in the large fixed N and Z limit in  Appendix B based on a scaling argument.

Equation (5) suggests that the leading correction BxZ of Eq. (2) arises, in part, from the two 1s electrons that are present in any atom of large Z and for which the density never becomes slowly varying. This is analogous to the Scott correction38 to Thomas–Fermi theory. Within Thomas–Fermi theory for atoms and ions, the majority of electrons are within a distance Z−1/3 from the nucleus, whereas the Scott correction applies39 to electrons within a distance Z−1. In the limit N = Z and at distances short compared to the Thomas–Fermi scale Z−1/3, the density is approximately independent of N, a Bohr-like atom.39 

In the limit 1 ≪ NZ, the Thomas–Fermi hydrogenic density becomes slowly varying on the scale of the exchange energy, but not on the scale of the correlation energy. This is demonstrated in  Appendix B using scaling arguments. The LDA is by definition exact for any uniform density; thus, the LDA exchange energy is exact as Z (where the density becomes high and locally uniform). This argument cannot determine the exactness of the LDA correlation energy but suggests that the convergence of the exact correlation energy to the LDA in the limit 1 ≪ NZ is slower than for the exchange energy. It is known19 from a direct semiclassical calculation that the LDA correlation energy is exact to the leading-order for heavy neutral atoms N = Z. At the next-to-leading order, the LDA coefficient is of the wrong sign,17 motivating the need for gradient corrections for real systems.40 

Figure 1 shows that the exact exchange energies23 of the neutral rare-gas atoms divided by Z5/3 vary from an almost-hydrogenic value at Z = 2 to a more Thomas–Fermi-like value at Z = 54. The more compact hydrogenic density has a more negative exchange energy for a given Z = N. The exact correlation energies17 divided by Z ln Z in Fig. 1 show an almost constant value of −0.014 from Z = 18 to 54. Interestingly, a much better coefficient of Z ln Z for neutral atoms of large Z can be found by approaching from the fixed N, large Z direction, as shown in Table I and  Appendix B. For the series of neutral atoms, one needs to extrapolate carefully17,18,20 to much larger Z to approach the asymptotic limit, but that limit is clearly pre-figured in the energies of real atoms. This suggests that widely predictive approximate functionals should be constrained to recover the correct large-Z asymptotics. In fact, LDA,3 PBE,9 SCAN,10 and acPBE18 sequentially improve10,18 the large-Z asymptotics. Functionals that satisfy sufficient exact constraints, including but not limited to Eqs. (2) and (3) with correct coefficients, can be predictive over the wide space of atoms, molecules, and solids, without being fitted to any bonded system. While that goal is not yet reached, the successes of SCAN41–46 suggest that it can be.

FIG. 1.

The scaled exact exchange23 (blue) and scaled exact correlation energies17 (orange) of the non-relativistic, neutral rare gas series. The smooth curves are a guide for the eye and not intended to be extrapolative.

FIG. 1.

The scaled exact exchange23 (blue) and scaled exact correlation energies17 (orange) of the non-relativistic, neutral rare gas series. The smooth curves are a guide for the eye and not intended to be extrapolative.

Close modal

Moreover, in the non-relativistic Z limit, the Periodic Table becomes perfectly periodic, with limiting first ionization energies for each column that increase across each row and for which the local density approximation for exchange and correlation may become relatively exact.47 

The explanations given here for the fundamental relevance of the LDA and its generalizations to real atoms and molecules were pioneered in Refs. 13, 14, 17, and (18), which focus on LDA’s correctness for neutral atoms, molecules, and solids of large atomic numbers. They constitute a third wave of such explanations. The first wave focused on LDA’s satisfaction of exact constraints on the exchange–correlation hole around an electron.48–50 The second wave, following Refs. 9 and 10, focused on the fact that the LDA inherits many (about nine) exact constraints on the density functional for the exchange–correlation energy from its appropriate norm, the uniform electron gas.

Our observation, that the large-N asymptotics of the exchange–correlation energy are “writ small” in the energies of atoms and atomic ions of very small N, is a contribution in support of this third wave of explanations. This prefiguration in small-Z neutral atoms is known,15 but not always emphasized.16,23,51 (The reasonable accuracy of LDA exchange for neutral atoms of small Z is also shown, for example, in Fig. 3 of Ref. 51.) The large-Z asymptotes by themselves are insufficient, since a pseudopotential method can remove them without significantly changing the energetics of the valence electrons, but Ref. 47 demonstrated its importance for energy differences.

Recent work has shown that the Perdew–Zunger self-interaction correction (PZ SIC), which makes any approximate density functional exact for non-overlapped one-electron densities, introduces errors of as much as 6% in the large-Z or slowly varying-density limit20 and that these errors can degrade the predicted equilibrium properties of molecules. Carefully scaling down the PZ SIC correction to zero in slowly varying regions yields much better results.21,22 The LYP correlation functional is highly accurate for first and second row molecules, but highly inaccurate for the uniform gas. To be appropriate for both solids and molecules, a density functional for exchange and correlation should recover the LDA in the uniform limit and correct the LDA for finite systems. Beyond this, Ref. 18 demonstrated that small modifications of the PBE functional (which, in its original form, recovers the LDA exactly in the appropriate limit and roughly predicts the next-order terms of the large-Z expansion), to recover exactly the next-order terms in Eqs. (2) and (3), moderately improve the PBE atomization energies of molecules. In future work, we hope to test the idea that the recovery of the next-order terms in Eqs. (2) and (3) can be used more generally to improve approximate density functionals.

The data that support the findings of this study are available within the article.

A.D.K. was supported by the Department of Energy (DOE), Office of Science (OS), Basic Energy Sciences (BES) through Grant No. DE-SC0012575 to the Energy Frontier Research Center: Center for Complex Materials from First Principles. B.S. was supported by DOE, OS, BES under Grant No. DE-SC0018331. P.B., K.W., and J.P.P. were supported by the U.S. National Science Foundation under Grant No. DMR-1939528. S.T.u.R.C. and H.T. were supported by DOE, OS, BES under Grant No. DE-SC0018194. K.B. was supported by DOE under Grant No. DE-FG02-08ER46496. K.B. thanks John Snyder, Jeremy Ovadia, Krishanu Ray, and Timothy Middlemas whose unpublished notes contributed to  Appendix B. We thank a referee for suggesting the possibility that β1(N) = 0 in Eq. (6) for all N within exact Kohn–Sham theory.

Consider the Hamiltonian Ĥ=Ĥ0+V^ee, where Ĥ0=T^+V^ext. Here, T^ is the kinetic energy operator, V^ext is the external potential or Coulomb attraction to the nuclei, and V^ee is the electron–electron Coulomb repulsion.

The exact Kohn–Sham correlation energy is defined as7 

Ec=Ψn|(T^+V^ee)|ΨnΦn|(T^+V^ee)|Φn.
(A1)

Here, Ψn, the ground-state wavefunction, minimizes the expectation value of Ĥ and defines the ground-state density n(r), while Φn, the exact Kohn–Sham wavefunction, is that ground eigenstate of the non-interacting or Kohn–Sham Hamiltonian ĤKS (a linear combination of at most a few Slater determinants) that yields the same ground-state density n(r). Consequently,

Ec=Ψn|Ĥ|ΨnΦn|Ĥ|Φn.
(A2)

Now, write Φn = Ψn + δΨ, and expand everything [including the ground state density n(r)] in powers of V^ee via degenerate perturbation theory. As V^ee0, Φn tends to the right linear combination of degenerate ground states of Ĥ0 and δΨ tends to zero like V^ee. Since Ψn minimizes the expectation value of Ĥ, the leading term of the correlation energy is of order (δΨ)2 or V^ee2. In the 1/Z expansion of Eq. (4), this is a contribution of order Z2(1/Z)2 = Z0.

Note that the exact Kohn–Sham exchange energy7Φn|V^ee|ΦnUH[n], where UH[n] is the Hartree electrostatic energy, can also differ substantially from its quantum chemical counterpart, the Hartree–Fock exchange energy.

Therefore, in the perturbation series for the total energy, the term of the correlation energy linear in Z (known precisely from quantum chemical calculations28) must appear in other components of the Kohn–Sham expansion (e.g., the exact Kohn–Sham exchange energy or exact Hartree potential).

In the standard degenerate or non-degenerate perturbation expansion of Eqs. (4)–(6), the density changes along with the coupling constant α that multiplies V^ee. A different perturbation expansion, in which the density is fixed at its physical or α = 1 value, was proposed by Görling and Levy (GL). In non-degenerate GL perturbation theory, the coefficient of the leading or α2 contribution to the correlation energy is given by Eq. (4) of Ref. 52.

In the heavy ion limit, with N fixed to a large value and Z, noninteracting Thomas–Fermi theory becomes relatively exact.25,33 One can show that the density of non-interacting Thomas–Fermi theory is given by

n(r)=13π22(μvext(r))3/2Θ(μvext(r)),
(B1)

where the chemical potential μ is a Lagrange multiplier determined by ∫n(r)d3r = N and Θ is a step function defining the turning surface. Let

ν=N/Z
(B2)

be the degree of ionization such that 0 < ν ≤ 1. We do not consider the case N > Z, as it has been found that NZ + 1 for real atoms.53 In the spherical case, one can rewrite the non-interacting Thomas–Fermi density in terms of the dimensionless variable x = r/rc, with rc being the turning surface radius. For hydrogenic densities, vext(r) = −Z/r, and thus, rc = −Z/μ [Note: μ < 0 in this case as vext(r) < 0 for all r]. Constraining the density to integrate to N determines rc=(18N2)1/3/Z, allowing the density to be rewritten as17,39

n(x)=Z2ν29π2(1/x1)3/2Θ(1x).
(B3)

In the limit ZN − 1, a hydrogenic density built up from Z-scaled hydrogen-atom orbitals obeying the hydrogen-atom aufbau principle becomes relatively exact. With the additional condition N ≫ 1, Eq. (B3) becomes relatively exact almost everywhere, so it usefully imitates a hydrogenic density with N ≫ 1. Here, “almost everywhere” means a region excluding electrons in the density tail or very near the nucleus, but including all of the electrons in the order of limits in which Z is followed by N. As we lack a simple, closed-form expression39 for a hydrogenic density constructed from hydrogen-atom orbitals in the large N limit, the non-interacting Thomas–Fermi density is needed to derive the asymptotic properties of hydrogenic densities.

From the scaling of n(x) with Z and ν, we see that the Wigner–Seitz radius rs(r) = [4πn(r)/3]−1/3 scales like

rs(x)=ν1/3Z2/33π1/32(1/x1)1/2Θ(1x).
(B4)

In the heavy ion limit, as N < Z with Z, rs tends to zero (the “high-density limit”), except near x = 1. In the heavy neutral atom limit, N = Z, rs → 0 as well. Note also that rc → 0 as Z, implying that the density localizes near the origin in this limit.

The LDA is exact for any uniform electron density and is relatively exact for any slowly varying electron density. We say that a density is slowly varying when its dimensionless gradients are less than order one. For exchange, the appropriate dimensionless density-gradient is s(r), defined on the scale of the Fermi wavevector kF(r)=[3π2n(r)]1/3,

s(r)=|n(r)|2(3π2)1/3n(r)4/3,
(B5)

and for correlation, the appropriate density-gradient is t(r), defined on the scale of the Thomas–Fermi screening wavevector ks=4kF/π,

t(r)=3π2161/3s(r)ϕ(ζ)rs(r)1/2.
(B6)

Here, ϕ(ζ) is a function of the spin-polarization ζ defined in Ref. 9. We know that the density is effectively spin-unpolarized in the limit of large N for which ϕ(ζ = 0) equals one and can be ignored throughout. The LDA exchange (correlation) energy will be relatively exact in the heavy atom limit if 0 < s(r) ≪ 1 [0 < t(r) ≪ 1] and will be exact if s(r) → 0 [t(r) → 0] as Z. These are presumptively sufficient conditions to determine the exactness of the LDA; the LDA may still be relatively exact even if they are violated.

From Eq. (B3),

|n(x)|=Z7/3ν5/34184/3π232x2(1/x1)1/2Θ(1x)+(1/x1)3/2δ(1x);
(B7)

the delta function is irrelevant, as x approaches 1 from below; thus, the dimensionless gradients scale as

s(x)=Z1/3ν1/33162/3x2(1/x1)3/2Θ(1x),
(B8)
t(x)=ν1/2183π21/2x2(1/x1)5/4Θ(1x).
(B9)

Both s and t are divergent as x → 0 and x → 1. However, for 0 < x < 1, s(x) → 0 as Z in both the heavy ion and neutral atom limits; thus, the LDA exchange energy becomes relatively exact as Z. This constitutes a simple argument for the exactness of the LDA exchange energy as found in Refs. 31 and 34. However, in the neutral atom limit, t(x)O(1), not characteristic of a slowly varying density. As the condition 0 ≤ t ≪ 1 is only sufficient, our scaling argument cannot determine if the LDA correlation energy is exact in the heavy neutral atom and heavy ion limits.

This scaling analysis, applied to the self-consistent Thomas–Fermi density of a neutral atom, shows that sZ−1/3, but tO(1) as Z. As noted in Ref. 13, the LDA still seems to get the correct leading term in Eq. (3), although PBE (which employs t) preserves the correct leading term and improves the next term in Eq. (3). Reference 13 also shows that s is already rather small in the closed-shell atoms Kr and Rn.

The coefficient of the leading-order terms in Eq. (3) can be determined by the scaling properties of a given density. Moreover, the scaling behavior of the density is sufficient to determine if the LDA exchange and correlation energies are separately exact in the high-density limit. In the high-density spin-unpolarized limit, the LDA correlation energy per particle tends to εcLDA=c0lnrsc1, where c0 = (1 − ln 2)/π2 ≈ 0.031 090 7 and c1 ≈ 0.046 644 are known from many-body perturbation theory.54 Consider a family of densities n(r/rc) that can be written as a function of the dimensionless position variable r/rc, with rc being the turning-surface radius defined by μ = vext(rc). The density may be expressed in powers of the nuclear charge Z, the ionization degree 0 < ν = N/Z ≤ 1, and a position-dependent function f(r/rc),

n(r/rc)=Zανβf(r/rc),
(B10)

with α ≥ 0 and β being real numbers. The high-density expansion of the LDA correlation energy evaluated on this family of densities yields an asymptotic series of the form

EcLDA=c03NαlnZ+βlnνKLDAN.
(B11)

KLDA is a constant dependent upon f(r/rc). For a self-consistent Thomas–Fermi density (α = 2 and β = 0) in the neutral limit (ν = 1), this series yields the coefficient Ac = −2c0/3 = −0.020 727 in Eq. (3). It was shown in Eq. (B3) that the non-interacting Thomas–Fermi hydrogenic density satisfies this scaling behavior with α = 2 and β = −1; thus, its leading order asymptotic series of Eq. (3) is identical to that of the self-consistent Thomas–Fermi density in the neutral limit.

A concern might come to mind: Eq. (B11) has explicit Z dependence, but the perturbation series of Eq. (6) suggests that the exact correlation energy should not depend upon Z, when there are no emerging degeneracies. Recasting Eq. (B11) as

EcLDA=c0α3NlnNc03(βα)lnν+KLDAN
(B12)

allows for a more direct comparison with Eq. (6). The dominant or N ln N term in Eq. (B12) is consistent with the numeric results of Table I. The next or N term is known to be incorrect even for the neutral case. The LDA correction applies at best only to the limit Z with N/Z fixed and not to the limit Z with N fixed. There exists an exact K that is not recovered accurately by the LDA.

Note that the scaling behavior shared by hydrogenic and self-consistent Thomas–Fermi densities, both derived from potentials that depend linearly upon Z, is not universal. For an isotropic harmonic oscillator potential vext(r) = kr2/2, where k > 0, the non-interacting Thomas–Fermi density is55 

n(r/rc)=k3/4N1/2241/23π2[1(r/rc)2]3/2Θ(1r/rc),rc=24Nk3/21/6,
(B13)

which yields a coefficient −5c0/12 of the N ln N term of Eq. (B12), when k = Z in the “neutral” limit (ν = 1). As a harmonic oscillator potential with linear dependence on Z is a mathematical artifice, no physical interpretation can be ascribed to the neutral limit. Setting k = Z4 would recover a 1/Z expansion for the total energy of the same form as Eq. (4).

1.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry
(
Macmillan
,
1982
).
2.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
3.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
4.
J. K.
Percus
,
Int. J. Quantum Chem.
13
,
89
(
1978
).
5.
M.
Levy
,
Proc. Natl. Acad. Sci. U. S. A.
76
,
6062
(
1979
).
6.
E. H.
Lieb
,
Int. J. Quantum Chem.
24
,
243
(
1983
).
7.
M.
Levy
and
J. P.
Perdew
,
Phys. Rev. A
32
,
2010
(
1985
).
8.
E. K. U.
Gross
,
M.
Petersilka
, and
T.
Grabo
,
Chemical Applications of Density-Functional Theory
, ACS Symposium Series Vol. 629 (
American Chemical Society
,
1996
), pp.
3
42
.
9.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
10.
J.
Sun
,
A.
Ruzsinszky
, and
J. P.
Perdew
,
Phys. Rev. Lett.
115
,
036402
(
2015
).
11.
P. A. M.
Dirac
,
Math. Proc. Cambridge Philos. Soc.
26
,
376
(
1930
).
12.
J.
Schwinger
,
Phys. Rev. A
24
,
2353
(
1981
).
13.
J. P.
Perdew
,
L. A.
Constantin
,
E.
Sagvolden
, and
K.
Burke
,
Phys. Rev. Lett.
97
,
223002
(
2006
).
14.
P.
Elliott
,
D.
Lee
,
A.
Cangi
, and
K.
Burke
,
Phys. Rev. Lett.
100
,
256406
(
2008
).
15.
A.
Becke
,
J. Chem. Phys.
84
,
4524
(
1986
).
16.
P.
Elliott
and
K.
Burke
,
Can. J. Chem.
87
,
1485
(
2009
).
17.
K.
Burke
,
A.
Cancio
,
T.
Gould
, and
S.
Pittalis
,
J. Chem. Phys.
145
,
054112
(
2016
).
18.
A.
Cancio
,
G. P.
Chen
,
B. T.
Krull
, and
K.
Burke
,
J. Chem. Phys.
149
,
084116
(
2018
).
19.
H.
Kunz
and
R.
Rueedi
,
Phys. Rev. A
81
,
032122
(
2010
).
20.
B.
Santra
and
J. P.
Perdew
,
J. Chem. Phys.
150
,
174106
(
2019
).
21.
R. R.
Zope
,
Y.
Yamamoto
,
C. M.
Diaz
,
T.
Baruah
,
J. E.
Peralta
,
K. A.
Jackson
,
B.
Santra
, and
J. P.
Perdew
,
J. Chem. Phys.
151
,
214108
(
2019
).
22.
P.
Bhattarai
,
K.
Wagle
,
C.
Shahi
,
Y.
Yamamoto
,
S.
Romero
,
B.
Santra
,
R. R.
Zope
,
K. A.
Jackson
,
J. E.
Peralta
, and
J. P.
Perdew
,
J. Chem. Phys.
152
,
214109
(
2020
).
23.
A. D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
).
24.
N. H.
March
and
R. J.
White
,
J. Phys. B
5
,
466
(
1972
).
25.
Y.
Tal
and
M.
Levy
,
J. Chem. Phys.
72
,
4009
(
1980
).
26.
A.
Savin
,
F.
Colonna
, and
R.
Pollet
,
Int. J. Quantum Chem.
93
,
166
(
2003
).
27.
M.
Levy
,
Phys. Rev. A
43
,
4637
(
1991
).
28.
J. P.
Perdew
,
M.
Ernzerhof
,
K.
Burke
, and
A.
Savin
,
Int. J. Quantum Chem.
61
,
197
(
1997
).
29.
J. W.
Hollett
and
P. M. W.
Gill
,
J. Chem. Phys.
134
,
114111
(
2011
).
30.
V. N.
Staroverov
,
G. E.
Scuseria
,
J. P.
Perdew
,
J.
Tao
, and
E. R.
Davidson
,
Phys. Rev. A
70
,
012502
(
2004
).
31.
J. G.
Conlon
,
Commun. Math. Phys.
88
,
133
(
1983
).
32.
N. H.
March
,
J. Phys. B
9
,
L73
(
1976
).
33.
B.-G.
Englert
,
Semiclassical Theory of Atoms
, Lecture Notes in Physics (
Springer-Verlag
,
Berlin; New York
,
1988
), pp.
vii; 401
.
34.
I. K.
Dmitrieva
and
G. I.
Plindov
,
Phys. Lett. A
55
,
3
(
1975
).
35.
J. C.
Snyder
,
J.
Ovadia
,
D.
Lee
,
K.
Ray
, and
K.
Burke
, in
Meeting Abstract 71-COMP
(
American Chemical Society
,
2011
).
36.
J. P.
Perdew
,
A.
Ruzsinszky
,
J.
Sun
, and
K.
Burke
,
J. Chem. Phys.
140
,
18A533
(
2014
).
37.
J.
Sun
,
J. P.
Perdew
, and
A.
Ruzsinszky
,
Proc. Natl. Acad. Sci. U. S. A.
112
,
685
(
2015
).
38.
J. M. C.
Scott
,
Philos. Mag.
43
,
859
(
1952
).
39.
O. J.
Heilmann
and
E. H.
Lieb
,
Phys. Rev. A
52
,
3628
(
1995
).
40.
S.-K.
Ma
and
K.
Brueckner
,
Phys. Rev.
165
,
18
(
1968
).
41.
J.
Sun
,
R.
Remsing
,
Y.
Zhang
,
Z.
Sun
,
A.
Ruzsinszky
,
H.
Peng
,
Z.
Yang
,
A.
Paul
,
U.
Waghmare
,
X.
Wu
,
M. L.
Klein
, and
J. P.
Perdew
,
Nat. Chem.
8
,
831
(
2016
).
42.
M.
Chen
,
H.-Y.
Ko
,
R. C.
Remsing
,
M. F. C.
Andrade
,
B.
Santra
,
Z.
Sun
,
A.
Selloni
,
R.
Car
,
M. L.
Klein
,
J. P.
Perdew
, and
X.
Wu
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
10846
(
2017
).
43.
L.
Goerigk
,
A.
Hansen
,
C.
Bauer
,
S.
Ehrlich
,
A.
Najibi
, and
S.
Grimme
,
Phys. Chem. Chem. Phys.
19
,
32184
(
2017
).
44.
J. W.
Furness
,
Y.
Zhang
,
C.
Lane
,
I. G.
Buda
,
B.
Barbiellini
,
R. S.
Markiewicz
,
A.
Bansil
, and
J.
Sun
,
Commun. Phys.
1
,
11
(
2018
).
45.
C.
Shahi
,
J.
Sun
, and
J. P.
Perdew
,
Phys. Rev. B
97
,
094111
(
2018
).
46.
Y.
Zhang
,
D. A.
Kitchaev
,
J.
Yang
,
T.
Chen
,
S. T.
Dacek
,
R. A.
Sarmiento-Pérez
,
M. A. L.
Marques
,
H.
Peng
,
G.
Ceder
,
J. P.
Perdew
, and
J.
Sun
,
npj Comput. Mater.
4
,
9
(
2018
).
47.
L. A.
Constantin
,
J.
Snyder
,
J. P.
Perdew
, and
K.
Burke
,
J. Chem. Phys.
133
,
241103
(
2010
).
48.
D. C.
Langreth
and
J. P.
Perdew
,
Solid State Commun.
17
,
1425
(
1975
).
49.
O.
Gunnarsson
and
B. I.
Lundqvist
,
Phys. Rev. B
13
,
4274
(
1976
).
50.
K.
Burke
,
J. P.
Perdew
, and
M.
Ernzerhof
,
J. Chem. Phys.
109
,
3760
(
1998
).
51.
A.
Pribram-Jones
,
D. A.
Gross
, and
K.
Burke
,
Annu. Rev. Phys. Chem.
66
,
283
(
2015
).
52.
A.
Görling
and
M.
Levy
,
Phys. Rev. B
47
,
13105
(
1993
).
53.
R.
Benguria
and
E. H.
Lieb
,
J. Phys. B: At. Mol. Phys.
18
,
1045
(
1985
).
54.
M.
Gell-Mann
and
K. A.
Brueckner
,
Phys. Rev.
106
,
364
(
1957
).
55.
M.
Brack
and
B. P.
van Zyl
,
Phys. Rev. Lett.
86
,
1574
(
2001
).