Leading correction to the local density approximation of the kinetic energy in one dimension

In the tens of thousands of density functional calculations published annually,¹ most employ the gradient of the density to estimate the exchange-correlation energy of the Kohn–Sham equations.² Such approximations begin from the gradient expansion of a slowly varying electron gas,³ which is then “generalized” to an integral over an energy density with some function of the density gradient.⁴ The first such attempt came already in 1968 when Ma and Bruckner showed that severe problems applying this gradient expansion approximation (GEA) for the correlation energy to atoms could be overcome by this procedure.⁵ Since then, a variety of procedures and philosophies have been used to construct such generalized gradient approximations (GGAs)⁶ and other more sophisticated constructions.⁷ Some GGAs are more accurate and popular in chemistry,^8,9 while others work better for (weakly correlated) materials.¹⁰ This diversity reflects the ambiguity in their derivation. The older, simpler local density approximation (LDA)^2,11 is uniquely determined by the energy of the uniform electron gas.^12,13

Long ago, Lieb and Simon proved that, for any electronic system, the relative error in Thomas–Fermi (TF) theory vanishes in a well-defined semiclassical limit in which the particle number tends to infinity.^14–16 Much work since then studies corrections to this limit order-by-order, including extensions of Thomas–Fermi theory.¹⁷ Such work is sometimes limited to atoms where spherical symmetry simplifies the situation. Englert beautifully summarized the work with Schwinger on this subject.^17,18 However, this problem is complicated by the interaction between electrons, the Coulomb attraction to nuclei, and the complexities of semiclassics in three dimensions.

This work studies the origin of the errors in applying the gradient expansion in the simplest relevant case, namely, the kinetic energy of non-interacting electrons in one dimension. This is not of quantitative relevance to realistic electronic structure calculations. The primary purpose is the construction of a mathematical framework in which this question can be directly addressed, and the errors of the gradient expansion explicitly identified and calculated in a systematic expansion in powers of ℏ. We show that, for simple model cases, the formalism is exact, and we also calculate the order-by-order expansion, finding great quantitative improvements in energies when the corrections are accounted for. We discuss the nature of these corrections and how they might be incorporated in density functional approximations.

Consider a symmetric potential v(x), with zero chosen so that v(0) = 0 and which could tend to D, the well-depth, at large x. Let ε_j be the eigenvalues of the Schrödinger equation using (Hartree) atomic units (setting m = ℏ = 1), and let M be the highest bound state if there is one. The number staircase is

where Θ(x) is the Heaviside step function, i.e., this is the number of states with ε_j < ε. Next, consider a smooth monotonic function I(ε) such that

where $j¯=j−1/2$⁠. The 1/2 comes from the Maslov index for two turning points.¹⁹ As ℏ → 0, a possible I(ε) is the classical action across the well, divided by π. We also define ε(y) as the inverse of I so that $εj=ε(j¯)$⁠. Then,

where ⌊x⌋ is the highest positive integer less than x, so N is the nearest integer to I. Next, define the periodic function

so that

To invert N_e(ε), we turn on a temperature that is much smaller than any energy or difference,

where f⁻¹(x) = 1 + exp(−x) and β is inversely proportional to temperature. We then define $μβ(N)$ as the inverse of N_β, $N∈R$⁠. For any finite temperature, $μβ(N)$ exists and is well-defined. We take β → ∞ at the end of the derivations and stop mentioning the temperature explicitly.

Figure 1 illustrates these functions for a Poschl–Teller well with D = 5/2 (see below). The smooth I(ε) generates the staircase N_e(ε), whose steps are rounded by the temperature, making it invertible. The difference has a sawtooth shape, crossing zero at the eigenvalues, so that N_e = I when both are (half)-integers.

We wish to develop an expression for the sum of the eigenvalues, which would be the total energy of N same-spin fermions in the well. Define the energy staircase,

i.e., the sum of eigenvalues with energy below ε. It will be especially useful to consider

and a well-known semiclassical result is²⁰ 

i.e., N_e(ε) determines the energy staircase. However, we really want S as a function of continuous number, $N$⁠, which is

This expression is well-defined for continuous values of non-negative $N$ (for non-zero temperature). As the temperature goes to zero, it becomes piece-wise linear, with changes of slope at integer values of $N$⁠, so that knowledge at integer values is sufficient to determine the whole function. Both $μ(N)$ and $G(μ(N))$ have step-like features, yielding

Note that $μ(N)=ε(N)$ for integers, where ε = I⁻¹, i.e., the discontinuous contributions in $μ(N)$ vanish identically at integers, so they are not needed to find S_N. Moreover, we change variables in the integration in G. If y = I(ε), then

where I′ = dI/dε and only the last term requires an integral over oscillations. Equation (12) is a central result, providing the machinery to construct the sum of the eigenvalues directly from I(ε), in continuous and discontinuous contributions. We define the first continuous term in G as

The value of ε(0), negative in Fig. 1, is irrelevant to G(N), as the step function vanishes for arguments less than 1/2, but not to J(N) or G_disc(N). Because J′ = I, S_N is fully determined by J(ε). Changing variables to ε(y) and integrating by parts yield the more succinct²¹ 

where II denotes integer-interpolation, i.e., a quantity that matches S at integers, but may differ in between. Figure 2 plots quantities vs $N$ for a PT well, showing that Eq. (14) agrees with $S(N)$ only at (half)-integers.

A harmonic oscillator is instructive. Here, I(ε) = ε/ω, so ε(y) = yω, I′ = 1/ω, J = ωN²/2, and G_disc vanishes because the average of ⟨y⟩ over one period vanishes if I′ is constant. For a particle in an infinite well,

Then, μ(N) = π²(N + 1/2)²/(2L²) and 1/I′ = π²(y + 1/2)/L². Then, $J(N)$ is trivial to integrate, but, because I′ varies, G_disc does not vanish,

The constant term gives no contribution, while the integral over y⟨y⟩ is −N/24, producing the exact answer S_N = π²N(N² + 3N/2 + 1/2)/(6L²). A less trivial example is the Poschl–Teller well of depth D,

Writing $αe=2D+1/4$⁠, then

yielding the eigenvalues

The simple result I′ = 1/(α_e − I) makes the calculation easy, using the same integral over ⟨y⟩ as before, giving

So far, this result might be considered a simple tautology. Its real use comes when a semiclassical expansion is performed. We multiply ℏ by a dimensionless number η and consider the limit as η → 0. Elementary analysis shows $εj(η)[v]=η2εjv/η2$ and, as η → 0,

where the expansion is known from WKB theory.²² Here,

is the classical action divided by π, and p(x) is the real part of the local classical momentum, $2(ε−v(x))$⁠, which yields the (zero-order) WKB eigenvalues. As ε⁽⁰⁾(0) = 0,

The WKB expansion is an expansion of individual eigenvalues in powers of η, keeping $ηj¯$ fixed, but the expansion of $S(η)(N)$ keeps $ηN$ fixed. One can both sum the WKB values to compare with S_N and compare S_N − S_N−1 with the Nth WKB eigenvalue. In general, these differ order-by-order (but infinite sums are identical). For example, summing WKB eigenvalues for the PT well produces an additional N/24 relative to S⁽⁰⁾(N). For N > 1, we expect the expansion of S(N) to outperform the sum of WKB eigenvalues to the same order, as the semiclassical approximation is used only at N and not at each individual j up to N, where it should be less accurate.

The leading correction is trickier to evaluate, due to a singularity as the turning points are approached. Define

where b(ε) is the turning point at energy ε, 0 < a < b(ε). As a → 0, a singularity develops, which must be canceled,

and ΔI⁽²⁾(ε) is found by taking a → 0. This cumbersome procedure can be elegantly avoided by an integral over a contour surrounding the turning points.²² Higher-order terms involve even stronger singularities. For a harmonic oscillator, Dunham²³ showed that all higher-order terms are identically zero so that WKB yields the exact answers. Likewise for a particle in a box, as all derivatives of v vanish, but $ΔI(2)=1/(82D)$ for PT.

However, we can instead evaluate the expansion for J(ε) and perform the energy integration before the spatial integral.²⁴ Consider

where b(ε) is the turning point at energy ε. For positive δ, this has no singularities and the order of integration can be reversed. As δ → 0, a singular term appears (of order $1/δ$⁠), which cancels that of Eq. (25). Thus,

Because of the oscillation, G_disc is already of higher-order than the continuous terms. Thus,

Because of the periodicity of ⟨y⟩, only the endpoints contribute to the integral as η → 0, yielding

Expanding ε(N) to second order in J⁽⁰⁾ yields

which is $N2/(16(2D))−N/12$ for the PT well.

Finally, we are ready to connect with density functional theory (DFT). For 1D same-spin non-interacting fermions in a slowly varying potential in an extended system, there is a well-known expansion of both the density n(x) and kinetic energy T in gradients of the potential,²⁵ 

and

where p_F(x) = p(ε_F, x) and ∫dxn(x) = N determines ε_F. Zero-order is Thomas–Fermi (TF) theory, and 2nd order is the gradient expansion. The combination $NεF(N)−S(N)$ yields T⁽⁰⁾ − 2ΔT⁽²⁾/3, agreeing with the semiclassical expansion of J, showing

where GEA denotes the (2nd-order) gradient expansion approximation. The semiclassical expansion for S_N reduces to the gradient expansion for extended systems, where the sawtooth contribution vanishes. However, $Gdisc(2)$ produces the exact leading-order correction to the local approximation for finite systems.

Some results for a generic Poschl–Teller well with α_e = 4.4 are given in Table I. The left side gives errors for the sum of eigenvalues, and the right side gives errors for the individual levels. The second column of errors on the left shows the result of the 2nd-order GEA, which sometimes worsens results relative to TF theory. However, inclusion of the correction reduces those errors by two orders of magnitude. Addition of the next order reduces the errors to the microHartree range. Deeper wells are even more favorable. Switching to the right side, for this well, even the eigenvalues are better approximated by differences in the sums within TF theory, but comparison of their second-order contributions shows that, in the asymptotic limit, WKB will have smaller errors than ΔTF for the top 1/6th of the levels. Table II repeats this calculation for D = 1, which binds one particle with energy −1/2 relative to the outside and has a second level right at threshold, far from the semiclassical limit. The trends are the same, but errors are greater, and WKB does better for the eigenvalue at the top of the well.

Finally, consider density functionals. Simply invert Eq. (31) and insert the result into Eq. (32) to find²⁵ 

It is straightforward to convert Eq. (29) into a functional of the TF density for the present circumstances. For potentials with a parabolic minimum at the origin,

where

Thus, $Gdisc(2)$ contains both highly local and non-local contributions (integrals over local functionals). Inserting the TF density for the PT well correctly yields N/24. Equation (35) is a correction unlike any in the existing literature; it has been derived, not devised.

The local density approximation applies to almost all situations. The potential functional correction to GEA of Eq. (29) applies to many circumstances, such as semi-infinite systems with surfaces, where the Maslov index differs, but must be generalized for multiple wells. On the other hand, when converted to a density functional [Eq. (35)], the form of the functional depends even further on the general class of problem. For example, the form differs from Eq. (35) for v(x) = |x|.

No general prescription was given for finding I(ε). For the DFT results, one needs only its well-defined asymptotic expansion. For simple model systems, the formulas used here suffice. However, adding any other function that vanishes at the eigenvalues generates equally viable candidates. Different I(ε) yield different continuous and discontinuous contributions but still yield the exact sums.

Many phenomena in DFT have a simple analog within this 1D world, as shown by two examples. The first is the well-known inaccuracy of functional derivatives of reasonably accurate semilocal approximations for the energy.²⁶ This infamous misbehavior of the LDA XC potential leads to highly inaccurate KS orbital eigenvalues. The analog here is N-particle density,

The archetype in 1D is the harmonic oscillator in TF theory, which yields the exact eigenenergies (and their sums), but whose density is highly inaccurate. The local approximation is exact for the harmonic potential, but not when small point-wise changes are made, as in Eq. (37). Only smooth changes in the potential should be expected to be correct in a local theory [the first four moments (0–3) of the TF density of the oscillator are exact]. Including the second-order correction yields densities that are singular at the turning points. This simply reflects the incompatibility of the order of limits by expanding in ℏ before differentiating.

The second is the well-known difficulty of semilocal functionals when bonds are stretched, a specific type of strong correlation.²⁷ Their failure has been traced to a delocalization error and related to curvatures of E vs $N$⁠. The same error shows up more strongly for the 1D kinetic energy. For one particle in two well-separated identical potentials, half the density ends up in each, leading to a factor of 4 reduction in the kinetic energy relative to the one-well result. However, a model for the double well is

where β is now a fixed large number, chosen to mimic the energy splitting between even and odd levels. The local approximation is much smoother

and produces a huge overestimate for j = 1. In fact, the equivalence of TF and WKB approximations breaks down, as there are now four turning points, leading to ambiguities analogous to the symmetry dilemma²⁸ for stretched H₂.

The expansion considered here is the analog of the standard Lieb–Simon (or Thomas–Fermi) expansion as Z → ∞ of neutral atoms.¹⁷ (Interacting) TF theory yields the dominant term, just as here, but the next two terms (Scott correction and Z^5/3 contribution) differ, due to the Coulomb potential, Coulomb interaction, and three dimensions. However, like those terms, one contribution to the leading correction is given by the gradient expansion and the other involves a difference between a term evaluated at the highest occupied level and at the bottom of the well.

This work represents a culmination of a series of earlier works,^29–31 which focused on finding the density as a functional of the potential but failed to yield systematically improved kinetic energies.³² The earlier results should prove useful when understood in the present context. While model results are not directly relevant to realistic calculations, the understanding achieved from previous studies has already had significant practical impact: the derivation of the parameter in the B88 functional,³³ an exact condition in PBEsol³⁴ three exact conditions in the SCAN meta-GGA,⁷ and the recently improved GGA correlation energy.³⁵ 

This paper aimed at the implications of this framework for DFT. Related work focuses on higher-order asymptotics of these sum formulas.²¹ The power of modern asymptotic analysis has recently been demonstrated³⁶ on v = |x|, finding contributions that are beyond all finite powers in the semiclassical expansion, with errors of 0.1 mH for just the second level. It is of tremendous interest to apply this machinery in three dimensions and to interacting systems.

I acknowledge funding from NSF (Grant No. CHE 1856165) and the University of Bristol for a Benjamin Meaker Professorship. I thank Attila Cangi and Raphael Ribeiro, whose earlier work inspired this advance, Michael Berry for instruction in semiclassics, and Chris Hughes for useful discussions.