Algebraic Density Functionals

A systematic strategy for the calculation of density functionals (DFs) consists in coding informations about the density and the energy into polynomials of the degrees of freedom of wave functions. DFs and Kohn-Sham potentials (KSPs) are then obtained by standard elimination procedures of such degrees of freedom between the polynomials. Numerical examples illustrate the formalism.

Existence theorems [1] for DFs do not provide directly constructive algorithms. Fortunately, the Kohn-Sham (KS) method [2] spares the construction of a "kinetic functional" and reduces energy and density calculations to the tuning of a local potential, v KS (r). Hence, a considerable amount of work has been dedicated to detailed estimates of electronic correlation energies and the corresponding KSPs, see for instance [3][4][5]. Many authors were also concerned with representability and stability questions, see for instance [6] and, for calculations in subspaces, see [7] and [8]. For cases where the mapping between potential and density shows singularities, see [9]. For reviews of the rich multiplicity of derivations of DFs and KS solutions and their properties, we refer to [10] and [11], and, for nuclear physics, to [12].
Local or quasi local approximations use the continuous infinity of values ρ(r), ∀ r, as the parameters of the problem. However, whether for atoms, molecules or nuclei, a finite number of parameters is enough to describe physical situations. For instance, Woods-Saxon nuclear profiles notoriously make good approximations, depending only on a handful of parameters, and it is easy to add a few parameters describing, for example, long tails and/or moderate oscillations of the density. High frequency oscillations are unlikely, for they might cost large excitation energies. The number of parameters to describe a nuclear density, therefore, can be restricted to maybe ∼ 10 at most; situations with ∼ 20 parameters are a luxury. For molecules, shapes are much more numerous, but a finite, while large number of parameters, truncating a list of multipoles for instance, still makes a reasonable frame. Practical DFs, therefore, can boil down to functions of a finite number of parameters. Integro-differential variations can then be replaced by simple derivatives.
This Letter shows how informations about both the density and the energy can be recast into polynomials. This allows eliminations of part of the parameters. Further polynomial manipulations locate energy extrema. Only density parameters are left. The same method gives KSPs. Finally we offer a discussion and conclusion.
A toy model of two particles illustrates the strategy. Consider the first 4 wave functions of the one-dimensional harmonic oscillator, ϕ 0 (r),. . . ,ϕ 3 (r), and a Slater determinant Φ made of positive and negative parity orbitals, where all 6 parameters, u, . . . , z are real numbers. Parity ensures orthogonality. Normalizations read, The density is a Gaussian, multiplied by a polynomial, with 3 independent coefficients, since the integral, equates to the particle number. Given a Hamiltonian H, the DF is defined by the constrained minimization [13,14], F [ρ] = Min Φ⇒ρ Φ |H| Φ , where the constraint, Φ ⇒ ρ, means that, given a 6 , a 4 , a 2 , (and hence a 0 ,) the parameters, u, v, . . . , z, must ensure that, Insert Eqs. (1) into Eq. (5) and take advantage of the harmonic oscillator basis states. The density constraint, Φ ⇒ ρ, then reduces into the 3 quadratic constraints, Recall there are 6 parameters, namely u, v, . . . , z to handle such 3 constraints, Eqs. (6), and also the 2 normalization constraints, Eqs. (2), and the energy. The latter, η = Φ |H| Φ , is also a polynomial. Indeed, if H is a one-body operator, with given matrix elements, ϕ i |H| ϕ j , the function η is quadratic with respect to the initial 6 parameters. If H contains two-body operators, then η has degree 4 with respect to the same 6 parameters, because of matrix elements, For a model of maximum simplicity, consider the sum of 2 single-particle harmonic oscillators, without their zero-point energy, A trivial elimination of 5 among the parameters, by means of the 5 constraints, gives, There remains to minimize η with respect to that last degree of freedom, v, left after the elimination. Let a "precursor" be a polynomial resulting from the definition of the energy, Eq. (7), after all constraints have eliminated a maximum number of degrees of freedom. In this toy model, a precursor is the left-hand side, P (η, v), of Eq. (8). Extrema of η correspond to the cancellation of, That solution which directly relates η to a 6 only is clearly spurious. The correct solution is, therefore, v = 0. Accordingly, the rigorous link between η and the parameters of the density is given by the polynomial equation, This is not an open formula of the form, η = F [a 2 , a 4 , a 6 ], but it provides roots for η at any realistic degree of numerical accuracy. Incidentally, it gives excited energies as well as minima, a well known property [15] of DFs. For numerical tests we show in Fig. 1 a plot of a density, ρ 1 (r) = (0.37861 + 3.64114r 2 − 3.30773r 4 + 1.21686r 6 )e −r 2 / √ π. Fig. 2 shows the η roots of Eq. (8) when v takes on all values that induce non-complex roots. We took great care to verify that the same values of η are obtained if v is run as a parameter of Eqs. (6) and Eqs. (2) to obtain u, w, x, y, and z by a brute force numerical algorithm without the elimination process. Figs. 3 and 4 show another density, ρ 2 (r) = (1.11314 − 0.0853536r 2 + 0.10563r 4 + 0.453499r 6 )e −r 2 / √ π, and the η roots obtained therefrom. In both cases of the density,  and in every other case that we tested at random, the minimum of η for v = 0 is confirmed. It may be noticed that, because of a better centering and milder oscillations for the density seen in Fig. 3, its pattern of energies, see Fig. 4, lies lower than that seen in Fig. 2.
A side issue arises here, and will also arise in all future models using this polynomial method, because a final minimization of η must be performed within a convex domain of densities: what conditions the coefficients a i (and their resulting a 0 ) must satisfy to maintain ρ(r) positive? This question was recently [16] solved by means of the Sturm criterion, for positive functions having positive Fourier transforms. The criterion gives the number of real roots of a polynomial, and can be used to ensure that a polynomial has no real roots. In the present model, the minimum energy among all densities corresponds to that determinant made of ϕ 0 and ϕ 1 , with density, ρ(r) = (1 + 2r 2 )e −r 2 / √ π. If we set, therefore, a 2 = 2 and a 4 = a 6 = 0 in Eq. (10), it is easy to verify that the root, η = 1, is found, as expected from Eq. (7) with u = x = 1 and v = w = y = z = 0.
The relation, τ = (η + 1) − ω, that connects the kinetic energy to the potential one, and to the full energy, transforms Eq. (13) into Eq. (10). Among generalizations of some interest, consider the radial density functional (RDF) theory [17], valid for nuclei and/or atoms, isolated, described by rotation invariant Hamiltonians. The constrained density minimization of energy [13,14] returns isotropic densities, with radial profiles, ρ(r), 0 r < ∞. Choose a positive measure, µ(r), having some physical relevance, and find a basis {p ℓn } of orthonornal polynomials, to define a set of "raw orbitals", ϕ nℓm (r) = Y ℓm (Ω r ) µ(r) p ℓn (r). Physical orbitals, φ α = nℓm c α nℓm ϕ nℓm , can be expanded on this basis, truncated at some maximum polynomial order. Furthermore, whatever the proton, neutron or electron numbers, states Φ q made of such orbitals can be mixed to generate isotropic densities ρ(r). Any such ρ is a polynomial of r multiplying [µ(r)] 2 . Simultaneously, the density, and the energy, will be polynomials with respect to the coefficients c α nℓm and/or the coefficients of configuration mixing, Ψ = q C q Φ q . There will therefore remain free parameters once a maximum number of parameters has been eliminated between the energy polynomial and the polynomials describing the density fit constraints. These free parameters stay in what we called a precursor polynomial. Energy minimization by manipulations of such a precursor requires again standard algebraic operations only. In a forthcoming paper, we shall describe a model more realistic than that used for this Letter.
This strategy gives directly the full DF, hence making the KS detour useless in principle. If spurious solutions [7,8] pop up, an analysis for their detection is easy, because the polynomial equation for the energy can only create a finite number of candidate solution branches. The KS approach may remain useful to verify n-and v-representability. A constructive derivation of KSPs is available, indeed. For instance, truncate again the harmonic oscillator basis at a maximum polynomial order n for single particle wave functions and let P be the projector upon the resulting, finite dimensional subspace for a system of N fermions, with their Hamiltonian H, or rather now, PHP. Any density in the subspace shows, besides its exponential factor, e −νr 2 , a polynomial of maximum order 2n. Given the kinetic energy operator T, choose a local potential v 0 (r), hence a one-body operator V 0 = N i=1 v 0 (r i ), hence again a one-body Hamiltonian H 0 = T + V 0 , so that the ground state of PH 0 P, a Slater determinant Φ 0 , be non degenerate and providing an approximate density ρ 0 for the system. For any density ρ in the subspace, the difference, ∆ρ = ρ − ρ 0 , with ∆ρ = 0, remains the product of e −νr 2 by a polynomial of maximum order 2n. (Here and in the following, the integral sign, , means r d−2 dr depending on the d-dimensional problem under consideration.) Expand ∆ρ in a basis of orthonormal polynomials S β (r), "constrained by vanishing averages" [18][19][20], ∆ρ(r) = 2n β=1 b β ξ β (r), where ξ β (r) = e −νr 2 S β (r). Again, given a determinant Φ with the parameters c α nℓm of its orbitals, or given a correlated state, Ψ = q C q Φ q , the constraints, Φ ⇒ b β or Ψ ⇒ b β , are polynomials of the parameters. Given H 0 , the polynomial method provides a polynomial K(κ; b 1 , . . . , b 2n ) for a reference functional, such that the lowest root of the equation, K = 0, represents the constrained minimum, κ ′ = Min Φ⇒b 1 ,...,b 2n Φ|H 0 |Φ , for the determinants in the subspace. In the same way, given the full H, the method provides a polynomial E(η; b 1 , . . . , b 2n ), the lowest η root of which is the constrained minimum, η ′ = Min Ψ⇒b 1 ,...,b 2n Ψ|H|Ψ , for correlated states in the subspace. Then it is trivial to derive from K and E a polynomial, Ω(ω; b 1 , . . . , b 2n ), for the difference, ω = η − κ. The diagonalization of PHP then reads, With the ratio, v β = −(∂Ω/∂b β )/(∂Ω/∂ω), representing ∂ω ′ /∂b β , define the one-body, local potential, v ∆ (r) = 2n β=1 v β ξ β (r). Let Φ be the ground state of P H 0 + N i=1 v ∆ (r i ) P. Notice that Φ|Pξ β P|Φ = Φ|ξ β |Φ . Then the energy E of Φ has derivatives, because of the orthornomality of the S β 's under the weight, e −2νr 2 . The numbers, b β0 = ρ 0 e −νr 2 S β , are easily pretabulated. The quantities, v β and (b β + b β0 ), are Legendre conjugates, and, moreover, ∂/∂(b β + b β0 ) = ∂/∂b β . The conditions, Eqs. (15), read as the diagonalization for a determinant Φ with the same density ρ as that of the eigenstate Ψ of PHP. The potential, P(v 0 + v ∆ )P, is a KSP valid for the subspace.
Rather than coefficients, we also fitted local values, ρ(r i ), i = 1, . . . , n, for even profiles. Such fits are also polynomial constraints with respect to the parameters. An optimization of the choice of such points r i makes a side issue, of some interest. But, at present, we feel that the simplest algebra results from fitting density polynomial coefficients. This is a highly nonlocal parametrization of ρ that deviates from the quasi-local tradition of the field. Other parametrizations, such as by moments of ρ from a Fourier transform, deserve consideration. In every case, our unconventional parametrization of ρ, in the frame of a polynomial algebra, creates a new zoology of DFs. Nothing of this zoology is known to us, but its interest is obvious, since manipulations of polynomials and properties of their roots, including bounds, are basic subjects. The number of available, exactly solvable models is huge, both at zero and finite temperature. It is limited only by computational power. For nuclei or atoms, the models will be "radial" [17], hence somewhat simple, although density tails adapted to Coulomb potentials, even screened ones, might demand many components. For electrons in molecules or extended systems (metals, thin layers, etc.), however, a necessary algebra of polynomials of 2 or 3 variables will burden the models. Anyhow, one can always test whether functionals from "smaller" models may remain valid for "larger" ones. Asymptotic properties might guide towards more traditional derivations of DFs. The models allow comparisons between the KS and the true kinetic energies of correlated systems. They also provide DF terms for those correlation energies due to interactions. SK acknowledges support from the National Research Foundation of South Africa and thanks CEA/Saclay for hospitality during part of this work. BG thanks Rhodes University for its hospitality during part of this work.