On the v-Representabilty Problem in Density Functional Theory: Application to Non-Interacting Systems

Based on a computational procedure for determining the functional derivative with respect to the density of any antisymmetric N-particle wave function for a non-interacting system that leads to the density, we devise a test as to whether or not a wave function known to lead to a given density corresponds to a solution of a Schrodinger equation for some potential. We examine explicitly the case of non-interacting systems described by Slater determinants. Numerical examples for the cases of a one-dimensional square-well potential with infinite walls and the harmonic oscillator potential illustrate the formalism.


Introduction
The Hohenberg-Kohn theorems [1] form the foundation of density functional theory (DFT), on which most ab initio electronic structure methods currently are based. The first theorem proves that the ground state of a many-electron system, described by an external potential, is uniquely determined by the ground state density (within an arbitrary constant). Therefore, the ground state density is the independent variable of DFT. In practice, it is still challenging to find the many-body wave function for the ground state of a potential. That task has been simplified by Kohn and Sham [2], who introduced a non-interacting system (Kohn-Sham system) with the same density as the interacting system that minimized the kinetic energy only. The existence of such a system is only assumed, but essential to the Kohn-Sham formulation of DFT.
In most of today's applications of DFT, only one direction of the Hohenberg-Kohn theorem is used, to find the ground state density for a given system described by an external potential. For the other direction, the inverse problem, it is a priori not clear if even a potential exists that leads to a given density through the solution of the Schrödinger equation, whether for the interacting or non-interacting system, the problem of v-representability. We refer to the literature [3][4][5][6][7][8][9][10][11] for discussions on the issue of v-representability. Specific examples of v-representable densities as well as non-v-representable densities have been identified.
An alternative formulation of the problem can be developed based on the N-representability property [3][4][5][6] of any density normalized to an integral number of particles, namely that a density can always be obtained from a wave function for a state of N-particles that, for fermions, is antisymmetric with respect to interchange of individual particle coordinates (and spins for electrons).
The v-representability problem can be stated as follows: given a density corresponding to a minimizing wave function, as obtained through the constrained search [3], determine whether or not the wave function is derived as the solution of the Schrödinger equation for a system of N electrons, either interacting, or non-interacting. A more general question is concerned with whether or not a wave function that leads to a given density is either interacting or non-interacting v-representable. A formal procedure that addresses this problem for general interacting systems has been published elsewhere [12]. In the following, a density defined with the following properties: it is everywhere non-negative, is normalized to the number of particles, and leads to a finite value of the kinetic energy.
In this work, we focus on the v-representability for non-interacting systems, even though all procedures presented here can be applied for the interacting system in the same way [12]. In previous work [13][14][15], we have derived closed expressions for the functional derivative of terms that only depend indirectly on the density, with respect to it. Using that formalism, we now show analytically that differentiating the non-interacting kinetic energy for a v-representable wave function leads to the negative potential (Kohn-Sham potential), up to a constant. This does not hold for arbitrary wave functions leading to the density. We obtain different wave functions leading to the same density by using a procedure by Cioslowski [16,17], who also showed that each of these wave functions is differentiable with respect to the density. This provides a test establishing whether a given wave function (non-interacting), leading to a density, is non-interacting v-representable. We illustrate these procedures on the one-dimensional particle in a box problem as well as the harmonic oscillator.

Background
The Hamiltonian describing an interacting system of N electrons in an external potential takes the usual form, with the operatorsV,T N andÛ N corresponding, respectively, to the external field, the kinetic energy and the inter-particle interaction (Coulomb repulsion for electrons). The ground-state energy of the system is given by the expectation value: where Ψ N g denotes the many-particle ground state ofĤ N . We use the notation, Ψ N → n(r), and say Ψ N leads to the density n(r), to denote the property: where n(r) denotes the single-particle density function normalized to the total number of particles, N. We now write Equation (2) in the form, in terms of the constrained search functional [3,18], Given a density, n(r), the constrained search examines all antisymmetric N-particle wave functions that lead to the density and delivers the state (in the absence of degeneracy) that produces the minimum value of Ψ N T N +Û N Ψ N .
As in the initial formulation of DFT by Kohn and Sham [2], we postulate the existence of a fictitious non-interacting N-particle system described by the Hamiltonian, under the action of an external potential,V s , whose ground-state density is identical to the density of the interacting system described byĤ N . In analogy with Equation (5), we define the constrained search functional [3], implying a search over all Slater determinants leading to a given density and returning that determinant that minimizes the expectation value of the kinetic energy for a system of N particles. Now, assume that the density is non-interacting v-representable for some potential, v s (r), so that from the second theorem of Hohenberg and Kohn [1], it follows that: defined within an arbitrary constant. The potential, v s (r), is such that the orbitals entering the construction of Φ N (r (N) ) arise through the solutions of a single-particle Schrödinger equation with eigenenergies j . Clearly, the orbitals minimize the expectation value of the kinetic energy, and, as such, form the Slater determinant identified by the constrained search in establishing T s [n]. For the case of Slater determinants, Equation (7) reduces to the direct sum:

Functional Differentiation
The process of performing the functional differentiation of expectation values with respect to the density, however, is not immediately evident: The dependence of a given wave function on the density is implicit, rather than explicitly displayed in terms of the density (or analytic functions of it). This difficulty is resolved through the use of the equidensity orbitals [19][20][21] that are written as explicit functionals of the density and, as shown in [21], form an orthonormal and complete basis in three-dimensional space. This property has been exploited in previous work [13][14][15] to allow the functional differentiation of the Coulomb energy associated with a Slater determinant with respect to the density in order to determine the Coulomb potential including the contribution of the exchange term. A complete discussion is given in [14]. Here, we apply that formalism to derivatives of the non-interacting kinetic energy with respect to the density.
We can obtain the functional derivative of T s [n] by differentiating the wave functions under the integral signs in Equation (7), and determine the potential-like quantity, where the double bars prohibit the operators from acting on the quantity beyond them. The reason for this is that functional derivatives of wave functions with respect to the density are generally not elements of the Hilbert space. Because of this feature, the operators in the last two equations act on the wave functions, either to their left or their right. It is seen that the combination of the constrained search and functional differentiation allows in principle the determination of all Slater determinants associated with given density that may correspond to a potential. For the case of a Slater determinant, the last expression takes the form, where the del operator acts on arguments to its right, and the last line gives the functional derivative as a function of the coordinates that may or may not correspond to a potential. In this expression, the quantity, T s [n], denotes the expectation value of the non-interacting kinetic energy operator,T N , with respect to any wave function that leads to a density.
In previous work [13][14][15], we derived a closed-form expression of the functional derivative of an orbital f (r) with respect to the density, δ f (r) δn(r ) , in terms of ordinary (spatial) derivatives. In the one-dimensional case, excluding terms that would lead to a constant shift in the potential (constants and expressions solely depending on x but not x ), the derivative reads Primes on functions denote spatial derivatives with respect to x, and Θ is the Heaviside step function. Expressions for the three-dimensional case can be found in references [14,15] for both the cases of cartesian and spherical coordinates, respectively.
In Appendix A, we derive an expression for the differentiation of the expression for the kinetic energy (one-dimensional for simplicity) using Equations (12) and (13) (see Equation (A4)). If all orbitals contributing to the kinetic energy originate from the same potential through the Schrödinger or Kohn-Sham equations, we show analytically that the differentiation leads to the negative potential (within an arbitrary constant). The three-dimensional case follows analogously.

Illustration of the Formalism
The calculations quoted here illustrate the formalism developed above for cases where the potential is given and the wave functions are known analytically. The goal is to reproduce the potential from the functional derivative of the kinetic energy with respect to the density.
For simplicity, we discuss one-dimensional systems in the non-interacting case where the wave function can be represented as a Slater determinant. Analogous procedures will apply in higher dimensions and the interacting framework.
For our first example, we consider N non-interacting particles in a box of length L bounded by potential walls of infinite height. The wave functions of the particles confined in the box are determined through the single-particle Schrödinger equation, (14) and vanish at the edges of the box. The n denote the eigenvalues. When v(x) = 0 (or a constant), the normalized single-particle wave functions are given by the expressions, with quantum numbers n = 1, 2, 3, . . . . The energies take the form, with the ground state density for a system of N Fermions given by n(x) = ∑ N n=1 |φ n (x)| 2 . Figure 1 illustrates the six lowest in energy eigenstates (L = 1) (left panel) and their moduli squared (corresponding to single-particle densities) on the right.  Figure 1. The six lowest in energy wave functions (left) and the corresponding single particle densities (right) for the particle-in-a-box model (L = 1). For illustrative purposes, all curves are shifted by the corresponding energies (see Equation (16)) in atomic unitsh = m = 1.
The second example discussed here is the one-dimensional harmonic oscillator. The potential is given by We make the usual transformation y = √ αx and α = mω h . The normalized solutions of the Schrödinger equation with the potential (17) read as where the H n (y) are Hermite polynomials, which are solutions of the differential equations H n (y) − 2yH n (y) + 2ny = 0, and n are integer quantum numbers ≥ 0. The corresponding energies are given by In the following discussion, we set m = ω =h = 1, which corresponds to α = 1 and V(x) = 1 2 x 2 . The six lowest in energy single particle wave functions and their corresponding moduli squared are shown in Figure 2.
For both examples, we use Equation (A4) to compute the potential. Since the wave functions are given analytically, the potentials can be calculated analytically. We use Mathematica R [22] to accomplish that task. We tested cases up to N = 6, and obtained every time the corresponding potential up to a constant. That is not surprising because we have shown in Appendix A that the formalism to calculate those functional derivatives leads to the potential. Equation (A9) also obtains a constant that does not have a physical meaning. Applying l'Hopital's rule twice at the right boundary for the particles in a box example we obtain the constant: x Figure 2. The six lowest in energy wave functions (left) and the corresponding single particle densities (right) for the Harmonic oscillator. All curves are shifted by the corresponding energies (see Equation (19)). The quadratic potential is shown in black.
Due to the exponential decay of the wave function and the density in the harmonic oscillator example, only the highest in energy of the occupied orbitals contributes to the constant. We find the constant to be equal to the eigenvalue of the highest occupied orbital.
So far, those tests have been performed occupying the lowest N orbitals. We find that any other combination, mimicking a non-ground state, e.g., by occupying orbitals 1, 2, 3, 4, 5 and 7, also leads to the same potential. The constants for the particles in the box example contains then the sum over occupied orbitals instead a sum from 1 to N, for the harmonic oscillator the constant is still determined by the highest occupied orbital.
The finding that the formalism yields the potential from the kinetic energy is not obvious because the condition δT s [n] δn(r) = −v s (r) implies the condition of the ground state. On the other hand, those states correspond to a local extremum so that the expression still holds.
Obtaining the potential by a functional differentiation for these examples may sound trivial, but it demonstrates the power of the procedure and the ability to calculate the functional derivative analytically. On the other hand, the procedure can be applied to cases where the the single particle orbitals are given numerically.
So far, we have shown the formalism determines the potential exact (up to a constant) from wave functions by functional differentiation of the kinetic energy, where it is known that all orbitals originate from the same potential, the condition for v-representability. The next step is to determine if a given set of orbitals, forming the density, is v-representable. That can be accomplished by calculating the potential v s (r) as shown above, and then checking if each of the orbitals fulfills the Kohn-Sham equation by comparing the right and left-hand sides of the expression: The determinant is v-representable when the equality holds, (for some α j ), for all orbitals forming the determinant. It is convinient for the comparison to normalize the left and right-hand side separately in order to be independent of the α j . Now, we apply this procedure to the case when the density is given. Generally, the task would be to find the wave function that minimizes the kinetic energy through the procedure of the constrained search [3]. Cioslowski [16,17] has established the formal procedure for generating the complete set of antisymmetric wave functions that lead to a given density. For a Slater determinant, the procedure is outlined in Appendix B. Based on his procedure, for each density, we construct a number of different Slater determinants that leads to it. For each Slater determinant, we differentiate the expectation value of kinetic energy with respect to the the density, obtaining v s (r). Finally, we apply Equation (21) for every orbital of the determinant. The results of the test for the case of the particle in the box and harmonic oscillator are shown in Figures 3 and 4, respectively, for a three particle system.
In the figures, the top rows exhibit the density, the second shows the set of three mutually orthonormal orbitals, leading to the density. The third row contains v s , obtained using Equation (A4), for each set of orbitals and the remaining rows compare each of the orbitals with the normalized output corresponding to Equation (21), where, in some, the output is multiplied by −1.
In both examples, the target density (top row) is obtained using the three lowest in energy orbitals (see Equations (15) and (18); L = 1) for panels A,B and C, the density in panel D is constructed by using orbitals 2, 3 and 4. We applied Cioslowski's procedure [16] for a single Slater determinant using auxiliary functions φ given in Tables 1 and 2 to obtain new Slater determinants leading to the given density. Those new orbitals are shown in the second row. If the exact orbitals are used, Cioslowski's procedure returns the same orbitals (column A); otherwise, different orbitals are obtained. We note that the auxiliary functions φ do not need to be orthogonal or normalized, as long as they are linearly independent. The new orbitals, although mutually orthonormal, can have quite strange forms, as seen in the second row. They do not even have to obey the symmetry of the problem, nor have a specific number of nodes. Table 1. Initial test functions used in Cioslowski's procedure [16] for the square well example corresponding to Figure 3. The functions in the first row are the exact orbitals (see Equation (15)). Table 2. Initial test functions used in Cioslowski's procedure [16] for the harmonic oscillator example corresponding to Figure 4. The functions in the first row are the exact orbitals (see Equation (18)). vs(x) Figure 3. Results for the square well example, as described in the text. The given density is shown in the top row, where the new orbitals obtained though Cioslowski's procedure [16] using initial functions from Table 1. The potential, obtained using functional differentiation of the expression of the kinetic energy, is plotted in the third row. The three rows at the bottom show the results of the test (see Equation (21) vs(x) Figure 4. Results analogous to those of Figure 3, but for the harmonic oscillator example. New orbitals were obtained through Cioslowski's procedure [16] using initial functions from Examining these cases, where different sets of orbitals form the same density (columns B and C), we see that the left-and right-hand side of Equation (21) are not the same (see bottom three rows), even though we know the density to be v-representable. This suggests that most of the sets generated by Cioslowski's procedure [16] may not correspond to a potential, yet all lead to the same density by construction. This behavior is expected from the first theorem of Hohenberg and Kohn: if v-representable, the ground state density determines uniquely the potential and hence the corresponding ground state wave function.

Column Test Orbital 1 Test Orbital 2 Test Orbital 3
We also applied the procedure to a different density to mimic an exited state (panel D). As seen in panels D4, D5 and D6, the set of obtained wave functions fail to satisfy Equation (21) and hence are not v-representable. The choice of the exact orbitals passes the test of v-representability, as already pointed out above.
In order to quantify the deviation of the left-and right-hand side of Equation (21), the L 2 norm of the difference is shown in Tables 3 and 4. L 2 norms smaller than the machine epsilon are marked as 0. Table 3. L 2 norm of the difference of left and right hand side of Equation (21) for the square well example, corresponding to rows 4-6 in Figure 1.

Conclusions
For a non-interacting system, we have shown how v-representability of a given wave function can be determined. This is accomplished by functional differentiation of the kinetic energy with respect to the density, and checking whether that derivative is a potential leading to the wave function. We have shown analytically (in one dimension), that, in the case of v-representability, the previously derived expressions for the functional derivative of orbitals lead to the correct result, the negative of the potential, when applied to terms depending on the diagonal part of the single-particle density matrix like the expression for the kinetic energy.
then the expression simplifies to: δT s δn(x ) = − V(x ) + c (A9) For a finite system, the x = ∞ in Equation (A8) signifies the right boundary. To derive Equation (A8), integration by parts is used:

Appendix B. Cioslowski's Procedure
For a single Slater determinant, we briefly describe how orthogonal orbitals, ψ(x), can be obtained, leading to a given density n(x) = ∑ N i=1 |ψ i | 2 , using Cioslowski's precedure [16]. From a set of linearly independent auxiliary functions φ(x), the ψ i (x) are given by the expression: The matrix s is given by The matrix S exists and has to be obtained self-consistently. That can be accomplished by iterating the equations above, e.g., using an identity matrix as a starting point.