Systematic corrections to the Thomas-Fermi approximation without a gradient expansion

We improve on the Thomas-Fermi approximation for the single-particle density of fermions by introducing inhomogeneity corrections. Rather than invoking a gradient expansion, we relate the density to the unitary evolution operator for the given effective potential energy and approximate this operator by a Suzuki-Trotter factorization. This yields a hierarchy of approximations, one for each approximate factorization. For the purpose of a first benchmarking, we examine the approximate densities for a few cases with known exact densities and observe a very satisfactory, and encouraging, performance. As a bonus, we also obtain a simple fourth-order leapfrog algorithm for the symplectic integration of classical equations of motion.


Introduction
All practical applications of density-functional theory (DFT) to systems of interacting particles require trustworthy approximations to the functionals for the kinetic energy and the interaction energy or-more relevant for the set of equations that one needs to solve self-consistently-to their functional derivatives. While the Kohn-Sham (KS) scheme [1] avoids approximations for the kinetic-energy functional (KEF), this comes at the high price of a CPU-costly solution of the eigenvalue-eigenstate problem for the effective single-particle Hamiltonian. The popular alternative to KS-DFT proceeds from the KEF in Thomas-Fermi (TF) approximation [2,3] and improves on that by the inclusion of inhomogeneity corrections in the form of gradient terms, with the von Weizsäcker term [4] as the leading correction.
This gradient expansion is notorious for its lack of convergence (see, e.g., [5]), the wrong sign of the von Weizsäcker term for one-dimensional systems [6], and the vanishing of all corrections for two-dimensional systems [6][7][8]-or so it seems [9]. While cures have been suggested, such as the use of a Padé approximant rather than the power series (see, e.g., [10]), or the partial re-summation of the series with the aid of Airy-averaging techniques [11][12][13][14][15], the situation is hardly satisfactory.
Recently, however, Ribeiro et al. [16] demonstrated that one can improve very substantially on the TF approximation without any gradient expansion at all. True, the method employed in [16] and related papers [17,18] is designed for one-dimensional problems and has so far resisted all attempts of extending it to two and threedimensional situations. But this work strongly encourages the search for other approximation schemes that do not rely on gradient expansions and are not subject to the said limitations. We are here reporting one scheme of this kind.
Our approach is based on the reformulation of DFT in which the effective singleparticle potential energy V (r) is a variational variable on equal footing with the singleparticle density n(r) [14,19]. All functionals of V (r) can be stated in terms of the effective single-particle Hamiltonian, and we relate them to the corresponding unitary evolution operator. That, in turn, is then systematically approximated by products of simpler unitary operators of the Suzuki-Trotter (ST) kind-known as split-operator approximations; see, e.g., [20]. There is, in particular, a relatively simple five-factor approximation that is correct to fourth order. It promises a vast improvement over the TF approximation without the high costs of the KS method and without the dimensional limitation of the Ribeiro et al. method.
The unitary ST approximations directly provide symplectic approximants for classical time evolution. The here developed higher-order ST approximations therefore have the potential to significantly improve upon standard second-order leapfrog algorithms [21,22] or even fourth-order Runge-Kutta methods.

Single-particle density and evolution operator
We consider a system of N unpolarized spin-1 2 fermions of mass m, subject to external forces that derive from the potential energy V ext (r). The energy functional which has n(r) and V (r) as well as the chemical potential µ as variables, is stationary at the actual values. Therefore, we have the three variational equations δµ : which we need to solve jointly. Here, E int [n] is the standard interaction-energy functional (IEF) of the Hohenberg-Kohn (HK) theorem [23], and the single-particle is related to the HK-KEF by a Legendre transformation, Equation (2a) yields n(r) for given V (r) and µ, whereas we obtain V (r) for given n(r) from equation (2b); and the correct normalization of n(r) to the particle number N is ensured by equation (2c) because it implies when combined with equation (2a).
In the KS formalism, we have with the eigenbras r| and eigenkets |r of the position operator R and the singleparticle Hamiltonian Here and throughout the paper, the upper-case letters P , R denote the quantum mechanical momentum and position operators, whereas the lower-case letters p, r stand for their classical counterparts. We can either accept equation (5) as determining the right-hand side in equation (2a), which is exact if we include the difference between the KEFs for interacting and non-interacting particles in the IEF, or we regard equation (5) as stating an approximation for the right-hand side in equation (2a)-an approximation that has a very good track record. Whichever point of view we adopt, we shall have to deal with equation (5). The task, then, is to evaluate, in a good approximation, the diagonal matrix element of the step function of µ − H without computing the eigenvalues of H and the corresponding single-particle wave functions ("orbitals"). In the tradition of the TF approximation and its refinements by gradient terms, and also in the spirit of the Ribeiro et al. work, we target an orbital-free (OF) formalism. In a first step toward this goal, we relate n(r) to the unitary evolution operator exp − it H(P , R) , where the integration path from t = −∞ to t = ∞ crosses the imaginary t axis in the lower half-plane [24][25][26]. Therefore, an approximation for exp − it H(P , R) provides a corresponding approximation for n(r).

Hierarchy of ST approximations-TF and beyond
In sections 3.1-3.5, we derive explicit approximations of the time-evolution operator and the according densities based on equation (7). These approximations successively incorporate exact constraints, eventually leading to a five-term factorization U 7 →0 that captures the exact evolution operator up to the fourth order in time t, is exact for the constant-force potential, and is consistent with the leading gradient correction. Table 1 summarizes our results. Table 1. Summary of the hierarchy of systematic ST approximations of the timeevolution operator discussed in sections 3.1-3.5. We require these approximations to be exact for the constant-force potential, and be at least of third order in time t for consistency with the leading gradient correction. The fourth-order approximation U 7 | →0 derived in section 3.4 meets all these requirements.
Exact for constant-force potential no no yes yes yes yes Correct up to order 2nd 2nd 2nd 2nd 2nd 4th Consistent with leading gradient correction no no no no no yes
As one verifies easily, this is recovered by the three-factor ST approximation (ST3) for which We note that the symmetrized approximation of equation (9) is not needed here, the two-factor approximations exp − it V (R) exp − it 1 2m P 2 or exp − it 1 2m P 2 exp − it V (R) would work just as well, and so would other asymmetric ways of sandwiching the kinetic-energy factor by potential-energy factors. Yet, the symmetric version does have an advantage: An approximation U (t), with the property U (t)U (−t) = 1, as is the case for U 3 , cannot have a leading error proportional to an even power of t. Since ST3 is obviously correct to first order in t, the terms ∝ t 2 must also be correct; indeed, the leading error is of order t 3 . In this sense, ST3 is a second-order approximation. If we interchange the roles of the kinetic and the potential energy in equation (9), we obtain another three-factor approximation of second order (ST3'), This gives n 3 (r) = 2 (dp 1 )(dp 2 )(dr 1 ) with for the corresponding approximation for n(r).

Beyond the TF approximation-ST5
Since U 3 and U 3 of equations (9) and (12) have errors for a constant-force situation, ∇V (r) = constant, they do not incorporate Langer's correction [27] which improves the description enormously at the border between the classically allowed and forbidden regions. Langer's correction is a key ingredient in the Ribeiro et al. method [16]. With this in mind we now consider a five-factor approximation (ST5) of the form and require that it is exact for vanishing dyadic ∇∇V (r). This identifies a set of admissible coefficients, parameterized by z. While U 3 and U 3 are particular cases of U 5 , they are not in the z-parameterized set. The error in U 5 is of third order; the leading error term is the double commutator Since ∆ 5 = 0 for all real z values (z = i/ √ 3 is not an option; see below), we need weighted sums such as for a third-order approximation, but then the unitary evolution operator is approximated by the non-unitary sum of two unitary operators, and this should better be avoided [20]. The resulting approximation for n(r) is where is the sum of a kinetic-energy term and a potential-energy term. For this kinetic energy to be positive, we need y 0 , y 1 , y 2 ≥ 0, hence − 1 3 ≤ z ≤ 1 3 . Although the z values are outside of this range in the z → ∞ limit in equation (17), this does provide an acceptable approximation, namely which is a three-factor approximation for the evolution operator. The corresponding expression for the density is with the Airy function Ai( ) and We emphasize that there is no gradient expansion in equation (21); the appearance of ∇V is a consequence of

Beyond the TF approximation-ST7
Since U 5 has a nonzero third-order error coefficient, ∆ 5 = 0, it does not account in full for the leading gradient correction, unless we resort to the non-unitary approximations of equation (18). ST5 also fails to reproduce the leading correction to the Wigner function of the evolution operator (see, e.g., [15]). Let us, therefore, consider a sevenfactor approximation (ST7) of the form for which with Here, we need x 1 + x 2 + x 3 + x 4 = 1 and and the cyclic analogs for y 1 and y 2 with to ensure that Langer's correction is incorporated. The third-order error coefficient is now We return to ST5 for When insisting on y 0 , y 1 , y 2 ≥ 0 to ensure positive kinetic-energy terms in H 7 , the minimal value of ∆ 7 is positive, ∆ 7 ≥ 1 12 − 1 18 √ 2 = 0.0048. While this lower bound is below that of ∆ 5 , ∆ 5 ≥ 1 72 = 0.0139, the improvement is quantitative, not qualitative. However, in the following section we demonstrate that ∆ 7 = 0 in the limit y 0 → 0 − , of the kind in equation (21).

A fourth-order approximation
The symmetric version of this limit uses for which and . For → 0, then, the error coefficient vanishes, ∆ 7 → 0, and we obtain as well as where The right-hand side in equation (33) has the symmetry property discussed at equation (11) and, as a consequence, is a fourth-order approximation to the evolution operator. It also reproduces the leading correction to the Wigner function of U (t); see appendix A. The simplicity of the expression is striking if one compares it with the fourth-order approximations in equations (44) and (144) in [20]; both have eleven factors.

Approximate densities
In this section, we present numerically tractable expressions for the approximate densities in equations (13), (19), (22), (26), and (34). Upon performing the p integrations in equation (13), we have where J D ( ) is the Bessel function of order D and is the effective Fermi wave number, with [x] Similarly, the single-particle density obtained from U 5 is where and The right-hand side of equation (38) applies also for n 7 (r) with k 5 replaced by The density for the fourth-order approximation is  (36) and (38), an equivalent expression for more efficient numerical implementation is derived in appendices B and C. When we replace k (x) 7 on the right-hand side of equation (42) by we obtain the z → ∞ limit of n 5 (r). We recover the TF density in equation (8)

Examples
In the following we assess the quality of the approximate densities presented in section 3.5. Figure 1 illustrates the systematic improvement of the ST approximation in equation (34) over the TF approximation in equation (8). The plot also shows the densities associated with the three-factor and five-factor ST approximations n 3 and n 5 ; see equation (13) and equation (19) with z = 0. The isotropic single-particle densities n(r) of spin-1 2 fermions with harmonic potential energy in D = 3 dimensions are obtained via equations (2a) and (2b) for a fixed chemical potential µ = 5 ω. In contrast to the quasi-classical n TF , n 3 decays smoothly across the boundary between classically allowed and forbidden regions, and the higher-order ST approximations approach the exact density n ex with successively higher accuracy. The results for n 7 are in line with the Airy-averaged density n Ai that also captures the leading quantum corrections, albeit via a systematic gradient expansion (the densities n Ai in figures 1 and 2 are calculated for a tempurature of one pico-Kelvin, i.e., very close to the ground state); see [15]. r [osc]ñ TF n Aĩ n 5 n 7 n ex n TF n Ai n 5 n 7 n ex Figure 2. Overview of particle densities with the same parameters as in figure 1, together with the scaled densitiesñ(r) = 1 2 r 2 n(r) that emphasize spatial regions of importance for integrated observables.
trap center (see figure 2), while their deviations from the exact density at small radii are less essential for global quantities like energy or particle number; cf. the weighted densitiesñ(r) = 1 2 r 2 n(r) in figure 2. A fair comparison of the performance of the various ST approximations in the context of equations (2a)-(2c) requires fixed µ and V , although this yields (slightly) different particle numbers for n 3 , n 5 , and n 7 ; from the physical point of view, µ is an auxilliary quantity whose sole purpose is to fix the particle number N for an actual physical system. In figure 3, we find excellent agreement between n 5 and the exact density for the Hooke atom of two Coulomb-interacting electrons in external harmonic confinement [28,29], although this semiclassical approximation should not be expected to excel for systems with such small particle numbers. We shall also consider n 7 for the Hooke atom, among other systems, as soon as the numerics is optimized. When approximating the interaction contribution to V (r) by the Hartree term which is void of exchange and correlation contributions, we obtain the densities in the inset of figure 3 via the self-consistent solution of equations (2a)-(2c). The disparity in results between the inset and the main figure serves as a reminder that a truly systematic description requires expressions for particle density and interaction functional at the same level of approximation.  [n] whose approximate nature is responsible for the differences to the main plot. We set = m = 2ω = 1.

Interaction energy functional
The various approximations for n(r) in terms of V (r) enter on the right-hand side in equation (2a). The partner equation (2b) requires a corresponding approximation for the IEF, where V pair (a) is the potential energy of an interacting pair of particles separated by distance a. For simplicity we are here content with considering spin-independent forces between the particles; spin-dependent forces, such as those between dipolar atoms, can be dealt with as well [30,31]. We can use any available approximate IEF for this purpose or, better, search for an IEF in an approximation that is consistent with that for the KEF (or, rather, its functional derivative). We are thus assigned the task of expressing the two-particle density n (2) (r, r ; r, r ) in terms of the singleparticle density n(r) = n (1) (r; r). Although this matter is beyond the scope of this article, let us briefly mention a possible strategy. Following Dirac [32], we accept the approximation n (2) (r, r ; r, r ) n(r)n(r ) − 1 2 n (1) (r; r )n (1) (r ; r) , which corresponds to splitting the IEF into a direct energy and an exchange energy contribution. For the single-particle density matrix n (1) (r; r ) we employ what equation (5) suggests, that is n (1) (r; r ) 2 r|η µ − H(P , R) |r , then relate the step function of µ − H to the evolution operator, to which we apply a suitable ST factorization. This yields n (1) (r; r ) in terms of V (r), and to complete the job we have to invert the mapping V (r) → n(r) in equation (2a) in a consistent approximation. For the TF approximation, for example, this results in δE int [n] δn(r) = (dr ) V pair (r )n(r + r ) D , for use in equation (2b); the first integral is the Hartree term of equation (45), there stated for the Coulomb interaction of the Hooke atom. We shall return to these matters in due course.

Digression: Beyond the leapfrog approximation
Besides supplying approximate solutions to the time-dependent Schrödinger equation in quantum mechanics, the split-operator approximations for the unitary evolution operator such as U 3 and U 7 →0 have counter parts in the context of symplectic approximations in solving Hamilton's equations of motion in classical mechanics. One obtains an approximate solution for d dt by the so-called leapfrog algorithm [21] (see also, e.g., [22]): (iii) repeat step (i); which takes one from time t to time t + τ (see the green dashed lines in figure 4). Clearly, this second-order leapfrog is the analog of the ST3 factorization in equation (9). There is also a leapfrog analog of ST3'. The fourth-order approximation in equation (33) yields a corresponding fourthorder leapfrog (see the blue solid lines in figure 4): where in step (iii). This fourth-order leapfrog promises substantial improvements over the second-order leapfrog. We are currently exploring this territory with very encouraging initial results [33].

Summary and outlook
We presented a method for incorporating inhomogeneity corrections into the approximate single-particle density without resorting to a gradient expansion. By insisting on no error for a constant force, we ensure that Langer's correction is properly accounted for. This improves the density very much at the transition from the classically allowed to the classically forbidden region, and beyond, as is confirmed in the examples of our benchmarking exercise. In order to advance the agenda, we shall next develop corresponding approximations for the interaction-energy functional, thereby setting the stage for the computation of self-consistent approximate solutions for a variety of fermion systems, among them the electron gas in atoms, molecules, and solids, as well as ultracold neutral atoms in optical traps. Generally, we expect the computational effort for our orbital-free approach to be linear in system size -in contrast to the cubic scaling of Kohn-Sham calculations in the self-consistent scheme. Some of the applications are likely to require the use of pseudopotentials [34,35]. Further, it will be interesting to investigate the analogous approximations for momentum-space densities [36,19,37,38], which are directly accessible in experiments with degenerate atom gases.
The fourth-order approximation in equation (33) gives rise to an approximate Wigner function of the time-evolution operator One reverts to the leading-order gradient correction of the Wigner function by expanding V (r ± 1 2 s) and V t (r + 1 2 a) around r up to second order in ∇ and then evaluating the resulting Gaussian integrals. Keeping only terms through second order in ∇, we so arrive at This agrees with the approximate Wigner function of U (t) in appendix B of [15], there obtained from a gradient expansion. We now regard this as an approximation to equation (A.1), which improves on the gradient expansion.

Appendix B. Derivation of the approximate densities
One way of deriving n 7 proceeds from the combination of equations (25) and (7) × (dp 1 )(dp 0 )(dp 2 ) exp − it 2 m y 1 p 2 1 + y 0 p 2 0 + y 2 p 2 2 × exp i p 1 · (r − r 1 ) + p 0 · (r 1 − r 2 ) + p 2 · (r 2 − r) , (B.1) where the evaluation of the Gaussian p integrals results in with a 2 as in equation (39). We use the identity to make the exponent linear in t, after which the resulting t integral yields a step function and equation (B.2) becomes The elementary p integration now establishes equation (38) with k 5 replaced by k 7 as in equation (41). Quite similarly, we get the analog of equation (B.2) for the density corresponding to the fourth-order approximation U 7 →0 in equation (33)    we can close the contour of the s integration by a large-radius semicircle in the lower half-plane when A − xB < 0, and by a semicircle in the upper half-plane when A − xB > 0. Since the only singularity of the integrand is at s = 0, we get a null result for K D (A, B) if A − xB < 0 and pick up the residue if A − xB > 0, A generating function of the Bessel functions, as the k = 0 term. It follows that as in equation (43). We obtain equation (42) The two positive solutions of ϕ(s) = φ > ϕ 0 are denoted s 1 (φ) and s 2 (φ), respectively, with 0 < s 1 (φ) < s 0 < s 2 (φ) and s 1 (φ) < 0 < s 2 (φ). For symmetry reasons, we have if the integration path crosses the imaginary axis at s = −i with > 0. We note that iϕ(−i ) is real and also that s n e iϕ(s) s=−i → 0 as 0 < → 0 for all powers n, negative or positive. It is expedient to integrate from s = −i to s = s 0 and then along the real axis from s 0 to ∞. This decomposes K D (A, B) into two pieces, . (C.10) Together with equation (C.7), this yields It follows that η(φ − ϕ 0 ) sin φ , (C.14) and (C.15) as well as Aφ η(φ) (C. 16) Equations (C.14)-(C.16) are well suited for a numerical evaluation.