On derivatives of the energy with respect to total electron number and orbital occupation numbers. A critique of Janak's theorem

Abstract

The relation between the derivative of the energy with respect to occupation number and the orbital energy, $\mathrm{\partial }E/\mathrm{\partial }{n}_i={ϵ}_i$, was first introduced by Slater for approximate total energy expressions such as Hartree–Fock and exchange-only LDA, and his derivation holds also for hybrid functionals. We argue that Janak's extension of this relation to (exact) Kohn–Sham density functional theory is not valid. The reason is the nonexistence of systems with noninteger electron number, and therefore of the derivative of the total energy with respect to electron number, $\mathrm{\partial }E/\mathrm{\partial }N$. How to handle the lack of a defined derivative $\mathrm{\partial }E/\mathrm{\partial }N$ at the integer point, is demonstrated using the Lagrange multiplier technique to enforce constraints. The well-known straight-line behaviour of the energy as derived from statistical physical considerations [J.P. Perdew, R.G. Parr, M. Levy and J.L. Balduz, Phys. Rev. Lett. 49, 1691 (1982)] for the average energy of a molecule in a macroscopic sample (‘dilute gas’) as a function of average electron number is not a property of a single molecule at T=0. One may choose to represent the energy of a molecule in the nonphysical domain of noninteger densities by a straight-line functional, but the arbitrariness of this choice precludes the drawing of physical conclusions from it.

GRAPHICAL ABSTRACT

Keywords:

• DFT
• Janak's theorem
• occupation numbers
• fractional electrons

1. Introduction

Fractional occupation numbers have remained popular ever since Slater introduced his expressions of the Hartree–Fock (HF) energy [1] and of the Xα energy (exchange-only LDA, XLDA) using explicit occupation numbers. His relation [2] (1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) holds for the mentioned model total energies. Slater applied his derivation explicitly to Hartree–Fock and to exchange-only LDA ($\mathrm{X}\alpha$), but his method only relies on the application of the chain rule to density and density matrix-dependent functionals, so is valid for all cited cases. Slater used relation (1(1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) ) to solve a problem that plagued his Xα scheme and all (semi-)local density functional approximations (LDFAs) since then: the orbital energy of the highest occupied orbital is not close to minus the ionisation energy $-I$. In the case of HF, it can be shown that the HOMO orbital energy ${ϵ}_H^{HF}$ is close to $-I$, according to Koopmans' frozen orbital approximation for the ion, ${ϵ}_H^{HF}=-I^{HF}(froz.orb.)$. The frozen orbital approximation leads to ${ϵ}_H^{HF}$ being a bit lower than the exact $-I$ (often about 1 eV). But for the local and semi-local density functional approximations (LDFAs) like LDA and GGAs, it is always some 4–6 eV too high. Slater introduced his so-called transition state (TS) method, which uses (1(1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) ), to obtain a good approximation to the ionisation energy as the HOMO orbital energy of a half-occupied orbital. See Ref. [3] for a recent review of orbital energies and (fractional) occupation numbers.

The interest in fractional occupation numbers has stimulated their introduction in density functional theory (DFT). Janak's theorem [4] (2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) expresses relation (1(1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) ) for the exact energy and Kohn–Sham orbital energies. The position of noninteger electron systems (implicit in (2(2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) )) in straightforward DFT for single molecules at T=0 is actually remarkable if one recognises that fractional electrons do not exist. One expects that everything that is ‘explained’ or ‘proven’ with fractional electrons, can also be treated with integer electron numbers, cf. Ref. [5]. In the case of molecular ionisation this is clear: On the one hand, Slater's relation (1(1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) ) only has meaning if occupation numbers can become fractional, and it can be used to prove Koopmans' theorem in HF using fractional electron numbers (see Section 5). It can also be used to calculate with the TS method the ionisation energies including orbital relaxation from orbital energies at fractional occupation number. On the other hand, Koopmans' theorem has first been formulated with just integer electron systems, and Slater's TS method just approximates total energy differences of integer electron systems. Slater's TS method is a mathematical device, it does not lay claim to any physical meaning of the fractional electron numbers. Janak's theorem in exact DFT, however, requires that physical meaning be given to (the energy of) noninteger electron systems (otherwise the derivative of Equation (2(2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) ) would not exist). The introduction of an ensemble description of noninteger electron systems by Perdew, Parr, Levy and Balduz (PPLB) [6] is very well known, but it should be kept in mind that this is based on statistical mechanics (using the grand canonical ensemble) [7]. It should be realised that the results hold for the average energy as a function of average electron number in, for instance, a macroscopic gas of molecules in the $T\to 0$ limit, not for a single atom or molecule.

We argue in this paper that Janak's theorem is in fact not a valid theorem in exact DFT (meaning DFT for single finite electron systems at T=0). Orbital occupation numbers do not properly belong to the KS theory. The KS system of noninteracting electrons will have the KS determinant as wavefunction. Orbitals either occur in that determinant (are occupied) or not. In order for the derivative in (2(2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) ) to be defined, the exact energy corresponding to a small increase and decrease of some $n_p$ (in a neighborhood of the integer occupation numbers) should be defined. It is not clear what that means in terms of the noninteracting electron system of the KS model (a KS orbital cannot ‘appear more’ or ‘appear less’ in the determinant). Increase of an orbital occupation above 1, would surely entail increase of the integrated density beyond the integer N -value. However, for such densities E is not defined in exact DFT, and therefore the derivative $\mathrm{\partial }E/\mathrm{\partial }N$ is not defined, and neither would be the derivative with respect to orbital occupation $\mathrm{\partial }E/\mathrm{\partial }{n}_i$.

We discuss in Section 2 the optimisation of the density according to the Hohenberg–Kohn theorem, where the constraint of fixed electron number N can, in principle, be treated with the Lagrange multiplier $\mu =\mathrm{\partial }E/\mathrm{\partial }N$. The derivative $\mathrm{\partial }E/\mathrm{\partial }N$, however, is not defined in DFT due to the domain on which the Hohenberg–Kohn functional (and Levy–Lieb and Lieb functionals) are defined being limited to N-conserving densities. The standard choice in functional analysis of $E\left[\rho \right]=+\mathrm{\infty }$ [8–11] for $\int \rho d\mathit{r}\ne N$ is problematic for definition of $\mathrm{\partial }E/\mathrm{\partial }N$. However, the essential arbitrariness of $\mathrm{\partial }E/\mathrm{\partial }N$ does not preclude the optimisation of the energy by density variation (at the integer point). It can be done either by constraining the density variation to normconserving ones, or by the Lagrange multiplier technique. The theory of constrained derivatives is discussed, in particular how the total functional derivative of the energy $\delta E/\delta \rho \left(\mathit{r}\right)$ can be split into the derivative under the constraint of constant N, written $\delta E/{\delta }_N\rho \left(\mathit{r}\right)$, and the derivative with respect to particle number N, $\mathrm{\partial }E/\mathrm{\partial }N$. Only the former is in principle needed, but is difficult to obtain. Therefore, a constraint like $\int \rho \mathrm{d}\mathit{r}=N$ is traditionally handled with the Lagrange multiplier technique. The latter, however, requires that $\mathrm{\partial }E/\mathrm{\partial }N$ (which enters the formalism as the force of constraint) is known. This raises the question how this technique should be applied if $\mathrm{\partial }E/\mathrm{\partial }N$ is not defined. Optimisation is possible, but requires that a suitable choice of the essentially arbitrary $\mathrm{\partial }E/\mathrm{\partial }N$ should be made (suitable means finite and continuous in the neighbourhood of the integer N point).

In Section 3, the consequences of $\mathrm{\partial }E/\mathrm{\partial }{n}_i$ not being defined are discussed. On the basis of the theory of constrained functional derivatives discussed in Section 2, a critique of Janak's theorem is given. The conclusion is that Janak's theorem has no validity in (exact) Kohn–Sham DFT. It is also investigated whether a KS-like formalism can be set up with possibly fractional occupation numbers (summing up to N) by dropping the KS physical model of N noninteracting electrons in N one-electron states in a local potential. This is shown to revert to the KS model.

Even if systems with a noninteger electron number are not physical, one may still define them as auxiliary systems to obtain useful results for integer electron systems. We have already mentioned Slater's TS method which uses ionisation by half an electron in order to derive an ionisation energy (physical quantity) from an orbital energy of an unphysical half-occupied orbital. This method defines, by a specific introduction of (fractional) occupation numbers, the energy of the auxiliary system as a (nearly) quadratic function of the occupation number. A different introduction of (fractional) occupation numbers, is the one based on statistical mechanical considerations of macroscopic systems at a finite temperature or judiciously chosen $T\to 0$ limit, as introduced by Mermin [12] and Perdew et al. [6, 7]. These occupation numbers describe the distribution of the electronic systems in the macroscopic sample over states with different occupations of the one-electron levels, see e.g. the derivation of the Fermi–Dirac distribution in [13]. So they are physical and describe averages. The search over ensembles of integer electron density matrices of PPLB [6] leads, in the $T\to 0$ limit, to linear interpolation (straight-line behaviour) of the average energy and the average electron densities for fractional electron number.

If an analogous procedure is followed for a single finite quantum system (atom or molecule), using an ensemble of density matrices with different electron numbers, density optimisation will again result in a specific (‘straight’) path through density space and a straight-line energy as a function of fractional electron number, see Section 4.1. However, a single system with noninteger electron number is nonphysical. This represents just one out of many possibilities to choose the arbitrary (because nonphysical) $\mathrm{\partial }E/\mathrm{\partial }N$ in this case. This choice should not be considered the only viable or even ‘exact’ definition [14–16] of $E\left[\rho \right]$ of a single quantum system (atom or molecule) for $\int \rho \mathrm{d}\mathit{r}\ne N$. The nonexistence of a noninteger electron system makes it impossible to ever obtain a benchmark value for the energy (wavefunction) of such a system.

It has been argued that the straight-line behaviour of the energy can be derived for an atom or molecule from size-consistency arguments, see Section 4.2. We caution that the underlying assumption is that the functional be local, which is not the case for the exact functional. A straight-line energy functional yields dissociation of any diatomic molecule into neutral atoms [6] (Section 4.3), which may be considered a strong point. However, we draw attention to the fact that this again relies on the unwarranted assumption of locality of the functional. Indeed, Section 4.4 discusses that the jump of the KS potential when an integer is passed, which is inherent in the straight-line functional, is a local property of a fragment in the dissociation and is not compatible with a correct description of dissociation of a heteronuclear diatomic.

Finally in Section 5, the possibility of a meaningful use of orbital occupation numbers in approximate energy expressions is highlighted. The freedom to introduce such additional variables in the total energy expressions is stressed, which may be done in Slater's way or in a different way. The Lagrange multiplier technique can be used to constrain the occupation numbers to physically meaningful values (typically 1 and 0). From this then emerges a physical meaning of the Lagrange multipliers $\left\{{\lambda }_i\right\}$ for this constraint, namely they become equal to the Lagrange multipliers $\left\{{ϵ}_i\right\}$ for the normalisation constraint, ${\lambda }_i={ϵ}_i$. It is shown that relation (1(1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) ) only holds if a particular choice is made for the dependence of the total energy on occupation numbers. Suppose one uses a model energy that allows to identify a derivative with respect to an occupation number with the related orbital energy according to the Slater relation (1(1) $\begin{array}{rl}& \text{Slater:}\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_i={ϵ}_i^{model}\\ & model=\mathrm{H}\mathrm{F},\mathrm{X}\alpha ,\mathrm{L}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{G}\mathrm{A},\mathrm{h}\mathrm{y}\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{d},\end{array}$(1) ), e.g. $\mathrm{\partial }{E}^{model}/\mathrm{\partial }{n}_{LUMO}={ϵ}_{LUMO}$. One may equate the LUMO orbital energy, in this example, to $(\mathrm{\partial }E/\mathrm{\partial }N{)}_{+}$, but the essential arbitrariness of the latter precludes any conclusion about the value of this orbital energy. Section 6 summarises.

We should stress that we are dealing here with systems in their ground state at T=0. The treatment of finite temperature effects and the statistical physics of large numbers of systems that can exchange energy (canonical ensembles) or also particles (grand canonical ensembles) with a bath are outside the scope of this treatment. The properties of individual quantum systems at T=0 are required in the first place, they determine the behaviour of ensembles. Individual systems can be in different states (excited, with different numbers of electrons) but the probabilities of occurrence of such states in statistical ensembles at finite temperatures is not the subject of Slater's or Janak's relations.

2. Energy derivatives and indeterminacy of the Lagrange multiplier for constant electron number N

2.1. The Lagrange multiplier technique and the force of constraint

We recall a few salient features of the Lagrange multiplier method [17, 18] for the imposition of constraints when finding an extremum of a function. Figure 1 shows an example in the case of an objective function f (the function for which an extremum has to be found) in two dimensions, under the constraint $h(x_1,x_2)=c$. At an arbitrary point, depicted to the left in the domain of feasible points (those obeying the constraint, i.e. on the curve $h(x_1,x_2)=c$), there is a component $\mathrm{\nabla }{f}_{//}$ of $\mathrm{\nabla }f$ along the curve of feasible points, which implies that the optimal point has not been reached. The perpendicular component is compensated by the ‘force of constraint’ exerted by the constraint, $-\lambda \mathrm{\nabla }h=-\mathrm{\nabla }{f}_{\mathrm{\perp }}$. [It is common practice to use the language of mechanical problems and denote f as a potential in which the minimum is searched and where the probe particle moves in the feasible domain until no forces are any more exerted.] The optimal point $(x_1^{\ast },x_2^{\ast })$ is characterised by the absence of a parallel force, $\mathrm{\nabla }{f}_{//}(x_1^{\ast },x_2^{\ast })=0$, and the Lagrange multiplier ${\lambda }^{\ast }$ at that point is a measure for the strength of the necessary force of constraint.

Figure 1. A two-dimensional example of application of the Lagrange multiplier: finding the minimum value of the function $f(x_1,x_2)$ under the constraint $h(x_1,x_2)=c$. The crossing point of h=c and $f=c_5$ represents an arbitrary point of the feasible set of points (those obeying the constraint) at which the component $\mathrm{\nabla }{f}_{//}$ in the domain of feasible points is still nonzero. $(x_1^{\ast },x_2^{\ast })$ is the optimal point, where $\mathrm{\nabla }{f}_{//}(x_1^{\ast },x_2^{\ast })=0$ and $\mathrm{\nabla }{f}_{\mathrm{\perp }}(x_1^{\ast },x_2^{\ast })$ is compensated by the force of constraint $-{\lambda }^{\ast }\mathrm{\nabla }h(x_1^{\ast },x_2^{\ast })$.

In textbook examples, the objective function is typically defined also outside the domain of feasible points, so that $\mathrm{\nabla }{f}_{\mathrm{\perp }}$ is defined. This is also necessary for straightforward application of the Lagrange multiplier technique, which requires all derivatives of the Lagrangean, without constraints, to be zero: (3) $\begin{array}{rl}L(x_1,x_2)& =f(x_1,x_2)-\lambda \left(h\right(x_1,x_2)-c)\\ \frac{\mathrm{\partial }L}{\mathrm{\partial }{x}_i}& =\frac{\mathrm{\partial }f}{\mathrm{\partial }{x}_i}-\lambda \frac{\mathrm{\partial }h}{\mathrm{\partial }{x}_i}=0\mathrm{\forall }i\end{array}$(3) It may happen of course that the objective function is not defined outside the feasible domain, for instance for physical reasons (as in our case of the Hohenberg–Kohn functional). That is not a problem if $\mathrm{\nabla }{f}_{//}(x_1,x_2)$ can be determined: the point $(x_1^{\ast },x_2^{\ast })$ where $\mathrm{\nabla }{f}_{//}(x_1^{\ast },x_2^{\ast })=0$ will still be the same. Application of the Lagrange multiplier technique, however, requires that $f(x_1,x_2)$ is defined outside the feasible domain, so that the derivative $\mathrm{\nabla }{f}_{\mathrm{\perp }}(x_1,x_2)$ can be determined. As a way out, one can simply extend in an arbitrary but suitable way the definition of f into the nonfeasible domain, so the derivatives of f exist and $\mathrm{\nabla }{f}_{\mathrm{\perp }}$ is defined. The extension is arbitrary in the sense that the actual value of $\mathrm{\nabla }{f}_{\mathrm{\perp }}$ is unimportant, it only affects the magnitude of the force of constraint which in this case is only an auxiliary quantity (chosen by the user), not a physical quantity. The extension must be suitable in the sense that the derivative must exist. A nonexistent derivative, for instance different left and right derivatives at the point $(x_1^{\ast },x_2^{\ast })$, precludes the application of the Lagrange multiplier method; the whole purpose of this method is to enable the calculation of the full derivatives, without any possibly very difficult restrictions on those derivatives due to constraints.

2.2. Differentiation of $E\left[\rho \right]$ and the arbitrariness of $\mathrm{\partial }E/\mathrm{\partial }N$ and of the Lagrange multiplier for the constraint of fixed electron number

In the Hohenberg–Kohn theory, we have the situation that the objective functional $E_v\left[\rho \right]$ is not defined outside the feasible domain of densities that integrate to N electrons, $\int \rho \mathrm{d}\mathit{r}=N$. Part of the derivative $\delta E/\delta \rho$ is therefore undefined (we omit the subscript v when there is no risk of confusion). Following Parr and Bartolotti [19] (see also Gál [20]), we make this explicit by writing the density as the product of N and a ‘shape function’ $\sigma \left(\mathit{r}\right)$ which is normalised to 1: (4) $\rho \left(\mathit{r}\right)=N\frac{g\left(\mathit{r}\right)}{\int g\left({\mathit{r}}^{\prime }\right)\mathrm{d}{\mathit{r}}^{\prime }}\equiv N\sigma \left(\mathit{r}\right),\int \sigma \left(\mathit{r}\right)\mathrm{d}\mathit{r}=1$(4) where $g\left(\mathit{r}\right)$ can be any decent nonnegative function. The density variation $\delta \rho \left(\mathit{r}\right)$ can be split in N-conserving and shape-conserving components, (5) $\delta \rho \left(\mathit{r}\right)=N\delta \sigma \left(\mathit{r}\right)+\sigma \left(\mathit{r}\right)\delta N\equiv {\delta }_N\rho \left(\mathit{r}\right)+{\delta }_{\sigma }\rho \left(\mathit{r}\right),$(5) where the subscripts N and σ denote variations $\delta \rho \left(\mathit{r}\right)$ where N or σ are kept constant, respectively. $\delta \sigma$ should obviously integrate to zero (and should be orthogonal to σ, $\int \sigma \left(\mathit{r}\right)\delta \sigma \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$, see Appendix 1). The energy variation becomes (6) $\begin{array}{rl}\delta E\left[\rho \right]& =\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \left[\frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\delta E\left[\rho \right]}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}\right]\left({\delta }_N\rho \left(\mathit{r}\right)+{\delta }_{\sigma }\rho \left(\mathit{r}\right)\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+\int \frac{\delta E\left[\rho \right]}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\end{array}$(6) where we have introduced the N-conserving and shape-conserving constrained derivatives, indicated with subscript N and σ in the functional derivative, following Gál [20]. The constrained derivatives for functionals are like partial derivatives for functions, with of course special adaptations. In going from the second to the third line in (6(6) $\begin{array}{rl}\delta E\left[\rho \right]& =\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \left[\frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\delta E\left[\rho \right]}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}\right]\left({\delta }_N\rho \left(\mathit{r}\right)+{\delta }_{\sigma }\rho \left(\mathit{r}\right)\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+\int \frac{\delta E\left[\rho \right]}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\end{array}$(6) ), we have used the fact that the cross terms (N-conserving derivative working on shape-conserving density change, and vice versa) are zero. See Appendix 1 for further elaboration and definitions. Note the analogy with the partial derivatives of functions of more variables, $\delta f=\left(\mathrm{\partial }f/\mathrm{\partial }{x}_1\right)\delta x_1+\left(\mathrm{\partial }f/\mathrm{\partial }{x}_2\right)\delta x_2$, where there are also no cross terms. Now (7) $\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\sigma \left(\mathit{r}\right)\delta N\mathrm{d}\mathit{r}=\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r}\delta N=\frac{\mathrm{\partial }E\left[\rho \right]}{\mathrm{\partial }N}\delta N$(7) (cf. Equation (A8(A8) $\frac{\mathrm{\partial }E}{\mathrm{\partial }N}=\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\frac{\mathrm{\partial }\rho \left(\mathit{r}\right)}{\mathrm{\partial }N}\mathrm{d}\mathit{r}=\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r}.$(A8) )) which defines the partial derivative of $E\left[\rho \right]$ with respect to N, $\mathrm{\partial }E/\mathrm{\partial }N$. This is $\mathit{r}$-independent, a constant. However, it is undefined in HK theory, since $E\left[\rho \right]$ is not defined for noninteger N's in the neighborhood of an integer N.

The problem of the value of $E_v\left[\rho \right]$ outside the domain of normalised densities has been addressed in various ways, usually in discussions of the Hohenberg–Kohn functional F, defined according to (8) $\begin{array}{rl}{E}_v\left[\rho \right]=F\left[\rho \right]& +\int \rho v\mathrm{d}\mathit{r};\\ & F\left[\rho \right]=⟨\mathrm{\Psi }\left[\rho \right]|\hat{T}+\hat{W}|\mathrm{\Psi }\left[\rho \right]⟩\end{array}$(8) The expectation values of the kinetic energy $\left(\hat{T}\right)$ and the electron–electron Coulomb interaction energy $\left(\hat{W}\right)$ have been defined for fixed electron number (N-electron) systems according to the Hohenberg–Kohn, Levy–Lieb or Lieb prescriptions [9, 21]. The Lieb functional $F_L\left[\rho \right]$, which leads to a convex functional, is always used in mathematical treatments (see earlier work by Levy [22] and Valone [23] on searches over density matrices of N-electron states) (9) $\begin{array}{rl}\text{Lieb:}{F}_L\left[\rho \right]& =\underset{\hat{D}\to \rho }{inf}Tr\hat{D}\left(\hat{T}+\hat{W}\right)\end{array}$(9) (10) $\begin{array}{rl}\hat{D}=\sum _{i=1}^n{\lambda }_i|{\mathrm{\Psi }}_i^N⟩⟨{\mathrm{\Psi }}_i^N|;0& \le {\lambda }_i\le 1,\sum _{i=1}^n{\lambda }_i=1\end{array}$(10) Lieb [8] in his study of functional analysis aspects of DFT proposed to define $F_L\left[\rho \right]$ as $F_L\left[\rho \right]=+\mathrm{\infty }$ for densities that fall outside the domain of positive functions that integrate to N. Infinite values are well defined in the theory of convex functionals, and they are usually introduced to deal in a simple way with domain issues [9]. This choice has usually been followed in analyses of the differentiability of $F_L$ [9–11, 24, 25]. It is however very different from the choice of Ref. [6] to define $F_L\left[\rho \right]$ for such densities by an extension of the density matrices over which the search is performed to ensembles of density matrices with different integer N values, see Section 4.1. That search then leads to linear interpolation of the energy between the N- and $(N+1)$-electron ground state energies at the electron-rich side (with derivative $-A$) and between the N- and $(N-1)$-electron ground state energies at the electron deficient side (with derivative $-I$). But an infinite derivative or the lack of differentiability at the integer point (different left and right derivatives) invalidates the Euler–Lagrange equation (14(14) $\int \frac{\delta }{\delta \rho \left(\mathit{r}\right)}\left[E\left[\rho \right]-\mu \left(\int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}-N\right)\right]\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$(14) ) (the force of constraint cannot be determined). Lieb has also indicated that more regular choices can be made, such as proportionality of $F_L$ to $\int \rho \mathrm{d}\mathit{r}=N$ (N noninteger), e.g. $F_L\left[\rho \right]=(\int \rho \mathrm{d}\mathit{r})F_L\left[\sigma \right(\mathit{r}\left)\right]$ with straightforward derivative (see remarks at Theorem 3.3 [8]). We stress here the arbitrariness of the definition of $F_L$ in the noninteger N domain, hence the freedom to make a suitable choice.

In applications of the Lagrange multiplier method, it is required that $\mathrm{\partial }E\left[\rho \right]/\mathrm{\partial }N$ is a defined quantity. Then we should make sure our results are independent of the choice we have made. In principle (when trying to determine the N-electron ground state density ${\rho }_0^N$), we do not need $\mathrm{\partial }E\left[\rho \right]/\mathrm{\partial }N$ because we could only search in the domain of N-electron densities, so at ${\rho }_0^N$ we must have, using the density shape factor g of Equation (4(4) $\rho \left(\mathit{r}\right)=N\frac{g\left(\mathit{r}\right)}{\int g\left({\mathit{r}}^{\prime }\right)\mathrm{d}{\mathit{r}}^{\prime }}\equiv N\sigma \left(\mathit{r}\right),\int \sigma \left(\mathit{r}\right)\mathrm{d}\mathit{r}=1$(4) ) (cf. [19, 26]) (11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) where we have used the chain rule. Since [26] (12) $\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}=\frac{N}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\left[\delta (\mathit{r}-{\mathit{r}}^{\prime })-\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\right]$(12) we have from (11(11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) ) and (12(12) $\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}=\frac{N}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\left[\delta (\mathit{r}-{\mathit{r}}^{\prime })-\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\right]$(12) ) at ${\rho }_0^N$ (13) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(13) This shows that at ${\rho }_0^N$ the functional derivative $\delta E/\delta \rho \left(\mathit{r}\right)$ is not a function of $\mathit{r}$ but is a constant. This fits in with the fact that at ${\rho }_0^N$ the functional derivative for N-conserving variations of ρ, $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$, is zero, see next paragraph (and in agreement with Equation (11(11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) )). So the full derivative $\delta E\left[{\rho }_0^N\right]/\delta \rho \left(\mathit{r}\right)$ reduces to the constant (independent of $\mathit{r}$) $\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N$, as is evident by comparing (13(13) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(13) ) with (7(7) $\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\sigma \left(\mathit{r}\right)\delta N\mathrm{d}\mathit{r}=\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r}\delta N=\frac{\mathrm{\partial }E\left[\rho \right]}{\mathrm{\partial }N}\delta N$(7) ). We cannot derive this constant from Equation (13(13) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(13) ), since any constant value is compatible with that equation. That is in order, since we have required the functional derivative to be zero only with respect to shape variation, Equation (11(11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) ). We can expect to obtain a result that is compatible with any value of $\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N$. The fact that $\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N$ is not defined in DFT, since HK theory does not provide values of the functional $E\left[\rho \right]$ for N-nonconserving ρ, and therefore not of the derivative, is not a problem. It should suffice to only require the derivative to be zero for N-conserving variations of ρ, $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)=0$, or $\delta E\left[{\rho }_0^N\right]/\delta g\left(\mathit{r}\right)=0$. There is a risk that it is not recognised that the undefined constant is part of the total derivative, but not of the constrained derivative, see comments on Equation (17(17) $\begin{array}{rl}\delta E& =\int \mathrm{d}\mathit{r}\frac{\delta E}{\delta \rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\\ \mathrm{i}\mathrm{s}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{o}\delta E& =\int \mathrm{d}\mathit{r}\left(\frac{\delta E}{\delta \rho \left(\mathit{r}\right)}+C\right){\delta }_N\rho \left(\mathit{r}\right)\end{array}$(17) ) and see Appendix 2.

In practice, constrained derivatives are not used since they are hard or impossible to obtain, and the standard Lagrange multiplier method is used because it has the great advantage that it allows to work with unconstrained derivatives, which are usually straightforwardly obtained. We try to find the minimising density from (14) $\int \frac{\delta }{\delta \rho \left(\mathit{r}\right)}\left[E\left[\rho \right]-\mu \left(\int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}-N\right)\right]\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$(14) (N is integer). This should hold for arbitrary $\delta \rho$. We may first choose for $\delta \rho \left(\mathit{r}\right)$ just ${\delta }_{\sigma }\rho \left(\mathit{r}\right)$ and must have at the optimum density ${\rho }_0^N$ (15) $\int \left[\frac{\delta E\left[{\rho }_0^N\right]}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\mathrm{\partial }E\left[{\rho }_0^N\right]}{\mathrm{\partial }N}-\mu \right]{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0.$(15) Since $\int \left[\delta E\right[\rho \left]/{\delta }_N\rho \right(\mathit{r}\left)\right]{\delta }_{\sigma }\rho \left(\mathit{r}\right)=0$ (see Appendix 1) this yields ${\mu }^{\ast }=\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N$: the force of constraint is determined by the choice of $\mathrm{\partial }E/\mathrm{\partial }N$. Using next the remaining $\delta \rho \left(\mathit{r}\right)$ space of norm-conserving $\left\{{\delta }_N\rho \right(\mathit{r}\left)\right\}$, and substituting ${\mu }^{\ast }=\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N$, it is clear that at ${\rho }_0^N$ we must have (16) $\begin{array}{rl}& \int \frac{\delta E\left[{\rho }_0^N\right]}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0\\ & \to \int \frac{\delta E\left[{\rho }_0^N\right]}{{\delta }_N\rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0\to \frac{\delta E\left[{\rho }_0^N\right]}{{\delta }_N\rho \left(\mathit{r}\right)}=0\end{array}$(16) where the second line follows since we already had $\int \delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right){\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$, as used above. Here it does not suffice to use only the subspace $\left\{{\delta }_N\rho \right(\mathit{r}\left)\right\}$ and require that the N-conserving derivative is only zero for application on that subspace. That requirement alone does not completely fix $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$.

There is a pitfall lurking here. It has often been noted in the literature [5, 9, 19, 26–28] that, since only norm conserving density variations are allowed in HK theory, derivatives like $\delta E/\delta \rho \left(\mathit{r}\right)$ (or $\delta F\left[\rho \right]/\delta \rho \left(\mathit{r}\right)$) are only defined up to a constant: if $\int {\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$, then (17) $\begin{array}{rl}\delta E& =\int \mathrm{d}\mathit{r}\frac{\delta E}{\delta \rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\\ \mathrm{i}\mathrm{s}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{o}\delta E& =\int \mathrm{d}\mathit{r}\left(\frac{\delta E}{\delta \rho \left(\mathit{r}\right)}+C\right){\delta }_N\rho \left(\mathit{r}\right)\end{array}$(17) It then appears that the same argument can be applied to the constrained derivative $\delta E/{\delta }_N\rho \left(\mathit{r}\right)$, which is then also believed to be only defined up to a constant. However, full definition of the functional derivatives requires their ‘operation’ to be defined on arbitrary $\delta \rho$, not only normconserving ${\delta }_N\rho$. As we have seen, this prohibits a free constant in $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$. In the Hohenberg–Kohn theory, there is such an undefined constant in the full derivative $\delta E\left[{\rho }_0^N\right]/\delta \rho \left(\mathit{r}\right)$, which stems from the HK restriction to integer electron systems, causing $\mathrm{\partial }E/\mathrm{\partial }N$ to be an undefined quantity in the theory. A constant is not allowed in the constrained derivative $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$ since it would give $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$ a component in the subspace $\left\{{\delta }_{\sigma }\rho \right(\mathit{r}\left)\right\}$, to which it should be ‘orthogonal’ (note the analogy of $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$ with $\mathrm{\nabla }{f}_{//}$ of Figure 1). See also discussion in Appendix 2.

We have seen that it is advisable to define $E\left[N\sigma \right]$ as a continuous function of N in a neighbourhood of the integer value N, so that the derivative $\mathrm{\partial }E\left[N\sigma \right]/\mathrm{\partial }N$ exists. We mentioned such a choice, derived from Lieb, earlier. Parr and Yang [26, p. 84] suggest to take a parabolic fit over the three points $E(N-1),E\left(N\right)$ and $E(N+1)$. This yields a continuous $E\left(N\right)$ in the neighbourhood of N, with $\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N=-(I+A)/2$, i.e. Mulliken's electronegativity. Görling [5] has analysed the arbitrary constant in the total derivative of the energy (or equivalently the Hohenberg–Kohn functional $F\left[\rho \right]$) along similar lines, and fixes it by adhering to the energy behaviour at noninteger N according to the PPLB straight lines picture of Ref. [6]. This has the difficulty of different left and right derivatives at the integer point. It would also be possible to choose $E\left[\rho \right]$ independent of N in the neighbourhood of N, which would yield $\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N=0$. All of this is compatible with the fact that $\mathrm{\partial }E/\mathrm{\partial }N$ is arbitrary.

3. Janak's theorem: not a DFT theorem

Relation (2(2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) ) is widely quoted as Janak's theorem and is often considered a fundamental relation in Kohn–Sham DFT. However, there are difficulties with this relation. One can try to formulate relation (2(2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) ) for a mathematically defined energy where occupation numbers (that supposedly can vary continuously) are introduced. Janak [4] introduced occupation numbers in a total energy expression (called $\tilde{E}$), which we call $E^{KS1}$ because it has occupation numbers to the first power (see Section 5 for an alternative) (18) $\begin{array}{rl}{E}^{KS1}& =\sum _{\sigma =\alpha ,\beta }\sum _p{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}p,q\\ \sigma ,\tau =\alpha ,\beta \end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\\ & +E_{xc}[{\rho }^{\alpha },{\rho }^{\beta }],\\ {\rho }^{\sigma }\left(\mathit{r}\right)& =\sum _p{n}_p^{\sigma }|{\phi }_p^{\sigma }\left(\mathit{r}\right){|}^2\end{array}$(18) The linear dependency on occupation numbers in the expression for the density and for the one-electron terms, which had earlier been introduced in the Xα and Hartree–Fock total energy expressions, is arbitrary (cf. Section 5) but plausible. For integer values (1 and 0) of the occupation numbers, there is no effect. In order for $E^{KS1}$ to have a defined derivative with respect to an occupation number, it is necessary that $E_{xc}[{\rho }^{\alpha },{\rho }^{\beta }]$ is defined for densities with fractional electron numbers (a neighbourhood of an integer $n_p$). But there is a problem with the value of $E_{xc}\left[\rho \right]$ in Equation (18(18) $\begin{array}{rl}{E}^{KS1}& =\sum _{\sigma =\alpha ,\beta }\sum _p{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}p,q\\ \sigma ,\tau =\alpha ,\beta \end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\\ & +E_{xc}[{\rho }^{\alpha },{\rho }^{\beta }],\\ {\rho }^{\sigma }\left(\mathit{r}\right)& =\sum _p{n}_p^{\sigma }|{\phi }_p^{\sigma }\left(\mathit{r}\right){|}^2\end{array}$(18) ) at noninteger electron numbers. Equation (18(18) $\begin{array}{rl}{E}^{KS1}& =\sum _{\sigma =\alpha ,\beta }\sum _p{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}p,q\\ \sigma ,\tau =\alpha ,\beta \end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\\ & +E_{xc}[{\rho }^{\alpha },{\rho }^{\beta }],\\ {\rho }^{\sigma }\left(\mathit{r}\right)& =\sum _p{n}_p^{\sigma }|{\phi }_p^{\sigma }\left(\mathit{r}\right){|}^2\end{array}$(18) ) is in fact the defining equation for $E_{xc}\left[\rho \right]$, its value is determined by all other terms (in particular $E^{KS1}\left[\rho \right]$) which need to be defined. But at noninteger electron numbers, the energy $E^{KS1}\left[\rho \right]$ (which should be $E_v\left[\rho \right]$) is not defined, and therefore $E_{xc}\left[\rho \right]$ is not defined. There are simply no physical systems with a noninteger number of electrons, for which an energy $E_v\left[\rho \right]$ would be obtainable. Also the Kohn–Sham construction of a noninteracting system of electrons is not possible for fractional electron numbers, so $E_{xc}\left[\rho \right]$ at noninteger electron densities is not a physical quantity. It is of course possible to choose an energy for noninteger electron densities in some way or another, for instance using ensembles of integer electron systems (see Section 4.1), but then it has to be kept in mind that no physical information can follow from such a man-made choice. Janak [4] did not specify $E^{KS1}\left[\rho \right]$ for noninteger densities.

We may wonder if Janak's theorem (2(2) $\text{Janak:}\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_p^{KS}}={ϵ}_p^{KS}$(2) ) may at least have meaning at the integer point, with occupation numbers introduced according to $E^{KS1}$ (18(18) $\begin{array}{rl}{E}^{KS1}& =\sum _{\sigma =\alpha ,\beta }\sum _p{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}p,q\\ \sigma ,\tau =\alpha ,\beta \end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\\ & +E_{xc}[{\rho }^{\alpha },{\rho }^{\beta }],\\ {\rho }^{\sigma }\left(\mathit{r}\right)& =\sum _p{n}_p^{\sigma }|{\phi }_p^{\sigma }\left(\mathit{r}\right){|}^2\end{array}$(18) ). Let us consider explicitly the derivative (19) $\begin{array}{rl}{\left.\frac{\mathrm{\partial }{E}^{KS1}}{\mathrm{\partial }{n}_p^{\sigma }}\right|}_{{\rho }_0^N}& =⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}+v_{Coul}|{\phi }_p^{\sigma }⟩\\ & +\int {\left.\frac{\delta E_{xc}}{\delta \rho \left(\mathit{r}\right)}\right|}_{{\rho }_0^N}{\left.\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{n}_p^{\sigma }}\right|}_{{\rho }_0^N}\mathrm{d}\mathit{r}\\ & =⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}+v_{Coul}\\ & +(v_{xc}+C)|{\phi }_p^{\sigma }⟩={ϵ}_p^{\sigma }+C\end{array}$(19) (we consider for simplicity the spin-compensated case where $v_{xc}$ is the same for spin-up and spin-down orbitals). We have indicated the occurrence of an unknown constant $C=\mathrm{\partial }{E}_{xc}/\mathrm{\partial }N$. It is unknown because the derivative for a change $\delta N$ in the number of electrons is unknown, see Section 2. In the total derivative of $E_{xc}$, (20) $\frac{\delta E_{xc}\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}=\frac{\delta E_{xc}\left[{\rho }_0^N\right]}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\mathrm{\partial }{E}_{xc}\left[{\rho }_0^N\right]}{\mathrm{\partial }N},$(20) the N-conserving derivative $\delta E_{xc}\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$ taken at the ${\rho }_0^N$ point may be identified with the exchange–correlation potential $v_{xc}$ (see below), but still there is the undefined constant $\mathrm{\partial }{E}_{xc}/\mathrm{\partial }N$. So the derivative $\mathrm{\partial }E/\mathrm{\partial }{n}_p$ is not a defined quantity in Kohn–Sham DFT, not even at the ground-state density. It certainly is not defined at fractional electron densities, where not even the $\tilde{N}$-conserving derivative $\delta E\left[\rho \right]/{\delta }_{\tilde{N}}\rho \left(\mathit{r}\right),\int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=\tilde{N}=N+\omega$ is defined. At the N-electron ground-state density ${\rho }_0^N$, we are free to choose the constant $\mathrm{\partial }{E}_{xc}\left[{\rho }_0^N\right]/\mathrm{\partial }N$, but then we cannot draw conclusions about the physics of a system from that choice.

The undetermined constant in the functional derivative, i.e. in the xc potential, also arises in the derivation of the Kohn–Sham equations, but it is harmless there. In the traditional derivation of the KS equations [29], one minimises the KS expression for the energy with respect to variations in the density via the orbitals, from which the one-electron Kohn–Sham equations follow. The application of the chain rule in the variation of the $E_{xc}$ term in the energy leads to the exchange–correlation potential $v_{xc}\left(\mathit{r}\right)=\delta E_{xc}/\delta \rho \left(\mathit{r}\right)$, without any precaution being taken that the total derivative is not defined (i.e. not for N-nonconserving density variation). However, the chain rule for the $E_{xc}$ term gives (21) $\begin{array}{rl}\delta E_{xc}& =\sum _p^N\int \frac{\delta E_{xc}}{\delta {\phi }_p\left(\mathit{r}\right)}\delta {\phi }_p\left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\sum _p^N\int \int \frac{\delta E_{xc}}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta {\phi }_p\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }\delta {\phi }_p\left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\sum _p^N\int \left(\frac{\delta E_{xc}}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\mathrm{\partial }{E}_{xc}}{\mathrm{\partial }N}\right){\phi }_p^{\ast }\left(\mathit{r}\right)\delta {\phi }_p\left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\sum _p^N\int \left(v_{xc}\right(\mathit{r})+C){\phi }_p^{\ast }\left(\mathit{r}\right)\delta {\phi }_p\left(\mathit{r}\right)\mathrm{d}\mathit{r}\end{array}$(21) where $v_{xc}\left(\mathit{r}\right)$ is the defined derivative $\delta E_{xc}/{\delta }_N\rho \left(\mathit{r}\right)$ (because $E_{xc}$ is defined for N-conserving densities) and one is left with the undefined constant $C=\mathrm{\partial }{E}_{xc}/\mathrm{\partial }N$. But in this case, this arbitrary constant in the KS potential does not have any physical consequences. A constant in the potential (which extends over all space, including asymptotic regions) just shifts the whole eigenvalue spectrum up or down. This is the well-known gauge freedom of a local potential, which can be eliminated always by choosing the potential to go to zero at infinity, so that all calculations work with the same gauge and the orbital energies become comparable. This is commonly done. Although the commonly used definition $v_{xc}\left(\mathit{r}\right)=\delta E_{xc}/\delta \rho \left(\mathit{r}\right)$ is strictly speaking not correct, or not complete, since it ignores the fact that the HK theorem only allows N-conserving densities, we now see that this omission does not have any consequences. There is not a similar saving grace in the derivation of Janak's theorem, since there the undefined constant exactly expresses that it is undefined what the theorem claims to tell, namely how the energy changes under a change $\delta N=\delta n_p$.

One may wonder if fractional occupations could be meaningfully introduced by abandoning the Kohn–Sham model system of N noninteracting electrons in a local potential. Parr and Yang [26, § 7.6] have raised the question if then the Janak theorem could be put on a secure footing. They consider a generalisation of the noninteracting kinetic energy $T_s\left[\rho \right]$ to a Janak kinetic energy (22) $T_J\left[\rho \right]=\underset{n_i,{\psi }_i\to \rho }{min}\sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩$(22) where the search is over all possible $n_i$ $(0\le n_i\le 1)$ and orthonormal orbitals ${\psi }_i$ yielding the given density (constraining it to N electrons) (23) $\rho \left(\mathit{r}\right)=\sum _i^{\mathrm{\infty }}{n}_i|{\psi }_i\left(\mathit{r}\right){|}^2\sum _i{n}_i=N$(23) The exact total energy functional can now be written (24) $E\left[\rho \right]=T_J\left[\rho \right]+\int v\left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+W_{Hartree}+E_{xc}^J\left[\rho \right]$(24) where $W_{Hartree}=\left(1/2\right)\int \rho \left({\mathit{r}}_1\right)\rho \left({\mathit{r}}_2/r_{12}\right)\mathrm{d}{\mathit{r}}_1\mathrm{d}{\mathit{r}}_2$ and $E_{xc}^J\left[\rho \right]$ is defined by this equation and is different from the standard $E_{xc}\left[\rho \right]$ if $T_J\left[\rho \right]$ would be different from $T_s\left[\rho \right]$. The question is if this set-up might lead to noninteger $n_i$. Then occupation numbers would have been introduced, possibly fractional, without the physical model of Kohn and Sham of noninteracting electrons in a local potential, which is defined with integer (1 and 0) occupation numbers. With integer $n_i$, the model of Equations (22(22) $T_J\left[\rho \right]=\underset{n_i,{\psi }_i\to \rho }{min}\sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩$(22) ) and (24(24) $E\left[\rho \right]=T_J\left[\rho \right]+\int v\left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+W_{Hartree}+E_{xc}^J\left[\rho \right]$(24) ) would reduce to the standard noninteracting kinetic energy of the KS determinantal wavefunction with the given ρ and lowest kinetic energy.

For $E\left[\rho \right]$ of (24(24) $E\left[\rho \right]=T_J\left[\rho \right]+\int v\left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+W_{Hartree}+E_{xc}^J\left[\rho \right]$(24) ) to be defined, ρ must be a density belonging to a ground state wavefunction (HK) or at least an N-electron wavefunction or density matrix (Levy–Lieb). Then $E_{xc}^J\left[\rho \right]$ is only defined for such integer-N densities. Parr and Yang [26] derive, at a given set of $\left\{n_i\right\}$, the KS-like equations for the ${\psi }_i$ which have $n_i\ne 0$ (25) $\left[-\frac{1}{2}{\mathrm{\nabla }}^2+v_{eff}^J\left(\mathit{r}\right)\right]{\psi }_i={ϵ}_i{\psi }_i$(25) where $v_{eff}^J=v\left(\mathit{r}\right)+v_{Hartree}\left(\mathit{r}\right)+v_{xc}^J\left(\mathit{r}\right)$ with $v_{xc}^J\left(\mathit{r}\right)=\delta E_{xc}^J/\delta \rho \left(\mathit{r}\right)$. Differentiating the total energy with respect to an occupation number $n_i$ is again problematic since the energy $E[\rho +\delta \rho ]$ is not defined for a small N-nonconserving density change $\delta \rho =\delta n_i|{\psi }_i{|}^2$, and neither is $E_{xc}^J[\rho +\delta \rho ]$. Ignoring this problem and again applying the chain rule for the derivative of $E_{xc}^J$ yields (26) $\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_i}={ϵ}_i$(26) The behaviour of the energy (24(24) $E\left[\rho \right]=T_J\left[\rho \right]+\int v\left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+W_{Hartree}+E_{xc}^J\left[\rho \right]$(24) ) under occupation number changes according to Equation (26(26) $\frac{\mathrm{\partial }E}{\mathrm{\partial }{n}_i}={ϵ}_i$(26) ) is well known to lead to Aufbau, i.e. $n_i=1$ for the lowest orbitals emerging from the optimisation, and $n_a=0$ for the remaining ‘virtual’ orbitals. This is easily seen by considering infinitesimal occupation number changes. The Janak–Kohn–Sham system with $T_J\left[\rho \right]$ instead of $T_s\left[\rho \right]$ then simply reverts to the Kohn–Sham system.

This raises the question if the Janak kinetic energy is actually different from $T_s\left[\rho \right]$. We can prove that this is not the case, i.e. the orbitals resulting from the minimisation of Equation (22(22) $T_J\left[\rho \right]=\underset{n_i,{\psi }_i\to \rho }{min}\sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩$(22) ) are the KS orbitals belonging to ρ and the optimal occupation numbers obey Aufbau, so (27) $T_J\left[\rho \right]=T_s\left[\rho \right].$(27) Let us minimise $T_J\left[\rho \right]$ by varying $\left\{n_i\right\}$ and $\left\{{\psi }_i\right\}$ under the constraint that $\sum _i^{\mathrm{\infty }}{n}_i|{\psi }_i\left(\mathit{r}\right){|}^2=\rho \left(\mathit{r}\right)$ at each $\mathit{r}$, for which purpose we introduce the $\mathit{r}$-dependent Lagrange multiplier $\mu \left(\mathit{r}\right)$. Normalisation of the ${\psi }_i$ is maintained with the usual Lagrange multipliers ${ϵ}_i^{\prime }$. The Lagrangean $L_J$ becomes (28) $\begin{array}{rl}{L}_J=& \sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩\\ & +\int \mu \left(\mathit{r}\right)\left(\sum _i^{\mathrm{\infty }}{n}_i|{\psi }_i\left(\mathit{r}\right){|}^2-\rho \left(\mathit{r}\right)\right)\mathrm{d}\mathit{r}\\ & -\sum _i^{\mathrm{\infty }}{ϵ}_i^{\prime }\left(\int |{\psi }_i\left(\mathit{r}\right){|}^2\mathrm{d}\mathit{r}-1\right)\end{array}$(28) For a fixed set of $\left\{n_i\right\}$, optimisation of the orbitals [26] yields the equations (29) $\left(-\frac{1}{2}{\mathrm{\nabla }}^2+\mu \left(\mathit{r}\right)\right){\psi }_i={ϵ}_i{\psi }_i$(29) for the occupied orbitals, where ${ϵ}_i={ϵ}_i^{\prime }/n_i$. The Lagrange multiplier $\mu \left(\mathit{r}\right)$ acts as local potential in this one-electron equation. Orthogonality of the orbitals then follows from the hermiticity of the operator and need not be enforced separately [26]. Before determining $\mu \left(\mathit{r}\right)$ in the usual way from the constraints, we note that with a fixed potential $\mu \left(\mathit{r}\right)$ the energy of the noninteracting electron system with Janak kinetic energy would be (30) $E_s^J=\sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩+\int \mu \left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}$(30) Minimisation of this energy by variation of the $\left\{n_i\right\}$ and $\left\{{\psi }_i\right\}$ under the constraint $\sum _i^{\mathrm{\infty }}{n}_i|{\psi }_i\left(\mathit{r}\right){|}^2=\rho \left(\mathit{r}\right)$ leaves the $\int \mu \left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}$ term invariant and is therefore equivalent to (should lead to the same $\left\{n_i\right\}$ and $\left\{{\psi }_i\right\}$ as) the minimisation in the definition of Janak kinetic energy according to (22(22) $T_J\left[\rho \right]=\underset{n_i,{\psi }_i\to \rho }{min}\sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩$(22) ). Considering then variation of the $\left\{n_i\right\}$, we note that now there is no problem with taking the derivative with respect to $n_i$ in (30(30) $E_s^J=\sum _i^{\mathrm{\infty }}{n}_i⟨{\psi }_i|-\frac{1}{2}{\mathrm{\nabla }}^2|{\psi }_i⟩+\int \mu \left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}$(30) ) since an exact energy for a noninteger electron system is not required for its definition. The derivative will just be (31) $\frac{\mathrm{\partial }{E}_s^J}{\mathrm{\partial }{n}_i}=\int {\psi }_i^{\ast }\left(\mathit{r}\right)\left(-\frac{1}{2}{\mathrm{\nabla }}^2+\mu \left(\mathit{r}\right)\right){\psi }_i\left(\mathit{r}\right)\mathrm{d}\mathit{r}={ϵ}_i$(31) So the minimum will be obtained for Aufbau. If the fixed set of $\left\{n_i\right\}$ with which we started was inadvertently not the Aufbau choice of $\left\{n_i\right\}$, which we may call $\left\{n_i^A\right\}$, the condition $\sum _i{n}_i^A|{\psi }_i\left(\mathit{r}\right){|}^2=\rho \left(\mathit{r}\right)$ will not be obeyed with the present set of $\left\{{\psi }_i\right\}$ and $\left\{n_i\right\}$. So then we have to repeat the process with occupation numbers according to Aufbau. The potential ${\mu }^A\left(\mathit{r}\right)$ will become exactly the KS potential $v_s\left[\rho \right]\left(\mathit{r}\right)$ belonging to density $\rho \left(\mathit{r}\right)$, since that KS potential is unique. So $T_J\left[\rho \right]$ is not a new ‘Janak’ kinetic energy but it is just $T_s\left[\rho \right]$. There is no Janak kinetic energy and no Janak–Kohn–Sham model.

Recently Li et al. [16] studied the Janak construction with a different purpose, namely to see if in case of a noninteger ρ this would not lead to occupation of higher virtual KS orbitals but just to the fractional occupation of frontier orbitals (Aufbau). They concluded to this Aufbau behaviour of the occupation numbers. However, we have indicated that $E\left[\rho \right]$ is problematic in case ρ is a noninteger density. Even if a $T_J\left[\rho \right]$ exists for noninteger densities, it cannot be used in Equation (24(24) $E\left[\rho \right]=T_J\left[\rho \right]+\int v\left(\mathit{r}\right)\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+W_{Hartree}+E_{xc}^J\left[\rho \right]$(24) ), since with a noninteger density, $E\left[\rho \right]$ in that equation is not defined.

Valiev and Fernando [30] have objected against the Janak theorem on somewhat different grounds. They point out that the results of the work by Englisch and Englisch [31] imply that the energy of a noninteracting system at fractional values of the occupation numbers (except for the degenerate levels at the Fermi energy) would not be differentiable with respect to the density.

We are not considering here fractional occupations of degenerate levels at the Fermi energy, which is theoretically [8, 32] and practically [33, 34] a well understood case, and neither do we consider density matrix functional theory (DMFT). Using the spectral resolution of the 1-matrix $\gamma ({\mathit{x}}_1,{\mathit{x}}_1^{\prime })=\sum _p{n}_p{\psi }_p\left({\mathit{x}}_1\right){\psi }_p^{\ast }\left({\mathit{x}}_1^{\prime }\right),$ the optimisation in DMFT is with respect to both the natural orbitals and the occupation numbers. In DMFT, the latter are an intrinsic part of the theory, they typically become all fractional. The inequality constraints $0\le n_p\le 1$ have to be applied with the Karush–Kuhn–Tucker method, see Giesbertz and Baerends [35]. Orbital energies start to play a very different role, cf. Gilbert's [36] famous finding of degeneracy for all fractionally occupied natural orbitals, which usually means all natural orbitals.

4. Energy behaviour at noninteger electron numbers

Ever since the seminal PPLB paper [6], the straight-line behaviour of the energy over the $(N,N+1)$ and $(N,N-1)$ intervals has received much attention in the literature. It should be emphasised that PPLB addresses the issue of fluctuating particle number as can occur in macroscopic samples at finite temperature, and use the grand canonical ensemble of statistical mechanics to treat this variable particle number. They derive straight-line behaviour in the $T\to 0$ limit (see the full details treated in [7]). On the basis of this work, it is sometimes assumed that even for a single quantum system (atom or molecule) $\mathrm{\partial }E/\mathrm{\partial }N$ is NOT arbitrary, but has to be $-I$ on the $(N,N-1)$ interval and $-A$ on the $(N,N+1)$ interval, with discontinuous derivative at N. We do not feel that it is a correct interpretation of the results of Ref. [6] to consider this behaviour as mandatory or ‘exact’ for a single molecule. For the common case of DFT calculation on a single molecule at T=0, the interpolation of the density between the integer N and $N\pm 1$ ground-state densities and the corresponding straight-line behaviour of the energy constitute just one choice for the essentially arbitrary continuation (see Section 2) to nonphysical fractional electron number. We will discuss this in Section 4.1.

One can also argue in favour of linear energy behaviour on the basis of size-consistency requirements for dissociation of molecules, as has been done by Yang, Zhang and Ayers [37, 38]. This is the subject of Section 4.2, where we argue that it is not size-consistency but the (unwarranted) requirement of locality of the functional that leads to the linear behaviour. In Section 4.3, we recall that application of the straight-line energy (assuming it holds for single quantum systems) has the well-known success of describing dissociation of a molecule into integer electron fragments, but we caution that again the unwarranted local approximation is invoked. In Section 4.4, it is shown that the local approximation, which is inherent in the application of the straight-line energy to dissociation, also causes it to fail: the derivative discontinuity jump in the KS potential is not quantitatively correct, in that it does not lead to proper dissociation of a heterogeneous electron pair bond.

4.1. Straight-line energy behaviour from grand canonical ensemble considerations and for single molecules

The statistical mechanical approach of [6, 7] implies that the straight-line behaviour applies to the average energy of the molecules in a dilute macroscopic gas of molecules, that can exchange electrons with a reservoir, at the properly taken $T\to 0$ limit. Also a very simple consideration at just T=0 makes that clear. We can then take μ as a tunable parameter of the reservoir, regulating the energy involved in electron transfer to and from the reservoir. A physical realisation is discussed by Perdew in [7] taking the reservoir to be a metal with workfunction Φ having negligible coupling integrals with a far away molecule. Then at T=0 when $\mu =-\mathrm{\Phi }$ drops below minus the molecular ionisation energy, evidently the ground-state energy corresponds to all molecules giving up an electron to the reservoir, i.e. become ionised. Similarly, as soon as μ would rise above minus the electron affinity, all molecules would turn into negative ions. At the points $\mu =-I$ and $\mu =-A$, the number of electrons on the molecules is undetermined. This is also what Perdew observes (see Equation (27) of [7]) for a molecule and a metal with work function Φ as reservoir. At those specific values of μ, many situations are possible, each characterised by specific probabilities to find N and N−1 electrons on a molecule (for $\mu =-I$) or N and N+1 electrons (for $\mu =-A$). Perdew [7] discusses the wavefunctions one may construct for the reservoir-molecule system that yield specific probabilities for the possible integer number of electrons on a molecule. When there is negligible coupling between molecule and reservoir, the wavefunction construction can yield the same probabilities as a mixture of density matrices [6, 7]. This yields a clear interpretation of the noninteger electron number and the corresponding energy as the averages of these quantities over all molecules of the dilute gas. For instance, if at $\mu =-A$ m electrons (say) are taken up from the reservoir by the gas with a very large number M of molecules, they have to go (at zero temperature) to m molecules, the change of the average number of electrons per molecule is $\mathrm{\Delta }\overline{N}=m/M$ and the average energy change is $\mathrm{\Delta }\overline{E}=-mA/M$. One obtains the derivative of the average energy $\mathrm{\partial }\overline{E}/\mathrm{\partial }\overline{N}=\mathrm{\Delta }\overline{E}/\mathrm{\Delta }\overline{N}=(-mA/M)/\left(m/M\right)=-A$ at $\overline{N}=N+\omega$ electrons. Going to $\mu =-I$ and $\overline{N}=N-\omega$ electrons, the derivative makes a quantum jump at integer N to $-I$ because of the quantum nature of the molecules. At integer N, there is not a defined derivative $\mathrm{\partial }\overline{E}/\mathrm{\partial }\overline{N}$ (because these considerations are at T=0 and do not take into account the proper $T\to 0$ limit). It is to be emphasised that this behaviour of the average energy and its derivative with respect to the average number of electrons $\overline{N}$ do not imply the same properties for any individual atom or molecule. No actual system can have a fractional number of electrons. Such systems are fictitious and we cannot calculate a wavefunction for them or get any information such as the energy for them.

The arbitrariness of the energy of a single atom or molecule at noninteger density implies that one may construct a functional which has some convenient behaviour. It is important then, of course, to refrain from conclusions about the physical (integer electron) system from the arbitrary choice. As example, let us construct a straight-line behaviour of the energy for a single molecule at noninteger electron number along similar lines as in Ref. [6], but now we relinquish any statistical mechanical underpinning and are dealing frankly with the nonphysical system of a single molecule with fractional electron number. We first recall that following work by Levy [22] and Valone [39], Lieb [8] formulated the Lieb functional in terms of a constrained search over an ensemble of N-electron pure state density matrices, Equation (9(9) $\begin{array}{rl}\text{Lieb:}{F}_L\left[\rho \right]& =\underset{\hat{D}\to \rho }{inf}Tr\hat{D}\left(\hat{T}+\hat{W}\right)\end{array}$(9) ). The Lieb functional has the important property of being convex. The further step can be taken of considering a density that integrates to a noninteger number of electrons, where the choice is made that the density be produced by taking an ensemble of density matrices of states ${\mathrm{\Psi }}_i$ of different electron numbers $N_i$, (32) $\begin{array}{rl}{E}_v\left[\rho \right]& =\int \rho \left(\mathit{r}\right)v\left(\mathit{r}\right)\mathrm{d}\mathit{r}+F\left[\rho \right]\\ F\left[\rho \right]& \equiv \underset{\hat{\mathcal{D}}\to \rho }{min}Tr\left[\hat{\mathcal{D}}\right(\hat{T}+\hat{W}\left)\right]\\ \hat{\mathcal{D}}& =\sum _i{\lambda }_i|{\mathrm{\Psi }}_i⟩⟨{\mathrm{\Psi }}_i|,\sum _i{\lambda }_i=1,\\ Tr\left(\hat{\mathcal{D}}\hat{\rho }\right(\mathit{r}\left)\right)& =\sum _i{\lambda }_i{\rho }_i\left(\mathit{r}\right)=\rho \left(\mathit{r}\right),\\ \int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}& =\tilde{N}=\sum _i{\lambda }_i{N}_i,\end{array}$(32) where we indicate general possibly noninteger N with $\tilde{N}$ (see Eschrig [24] for a comprehensive treatment). With the generally assumed convexity of the energy as a function of number of electrons $\left(E\right(N)\le (E(N+1)+E(N-1)\left)/2\right)$, the minimum energy at $\tilde{N}$ in between, say, the integers N and N+1, will be straight line interpolation between $E_0^N$ and $E_0^{N+1}$. So at $\tilde{N}=N+\omega$ (33) $\begin{array}{rl}{\hat{\mathcal{D}}}^{min}& =(1-\omega )|{\mathrm{\Psi }}_0^N⟩⟨{\mathrm{\Psi }}_0^N|+\omega |{\mathrm{\Psi }}_0^{N+1}⟩⟨{\mathrm{\Psi }}_0^{N+1}|\\ \tilde{E}\left(\tilde{N}\right)& =(1-\omega )E_0^N+\omega E_0^{N+1}\\ \rho \left(\tilde{N}\right)& =(1-\omega ){\rho }_0^N+\omega {\rho }_0^{N+1}\end{array}$(33) However, $\tilde{E}\left(\tilde{N}\right)$ cannot be taken to be the exact energy of a noninteger electron system at T=0. Such a system does not exist, so it is not possible to determine its exact energy. The density $\rho \left(\tilde{N}\right)$ and energy $\tilde{E}\left(\tilde{N}\right)$ represent just one out of the many possible ways to continue the density and energy of an N-electron system into the nonphysical fractional electron domain. This choice leads to a $\tilde{N}$-derivative with the constant value $-A$ on the $(N,N+1)$ interval, and $-I$ on the $(N-1,N)$ interval. Given the arbitrariness of $\mathrm{\partial }E/\mathrm{\partial }N$ according to Section 2, this choice is possible, but it would not work in a density optimisation with a Lagrange multiplier for the $\int \rho \mathrm{d}\mathit{r}=N$ constraint, since the derivative is not continuous at N. Equations (32(32) $\begin{array}{rl}{E}_v\left[\rho \right]& =\int \rho \left(\mathit{r}\right)v\left(\mathit{r}\right)\mathrm{d}\mathit{r}+F\left[\rho \right]\\ F\left[\rho \right]& \equiv \underset{\hat{\mathcal{D}}\to \rho }{min}Tr\left[\hat{\mathcal{D}}\right(\hat{T}+\hat{W}\left)\right]\\ \hat{\mathcal{D}}& =\sum _i{\lambda }_i|{\mathrm{\Psi }}_i⟩⟨{\mathrm{\Psi }}_i|,\sum _i{\lambda }_i=1,\\ Tr\left(\hat{\mathcal{D}}\hat{\rho }\right(\mathit{r}\left)\right)& =\sum _i{\lambda }_i{\rho }_i\left(\mathit{r}\right)=\rho \left(\mathit{r}\right),\\ \int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}& =\tilde{N}=\sum _i{\lambda }_i{N}_i,\end{array}$(32) ) and (33(33) $\begin{array}{rl}{\hat{\mathcal{D}}}^{min}& =(1-\omega )|{\mathrm{\Psi }}_0^N⟩⟨{\mathrm{\Psi }}_0^N|+\omega |{\mathrm{\Psi }}_0^{N+1}⟩⟨{\mathrm{\Psi }}_0^{N+1}|\\ \tilde{E}\left(\tilde{N}\right)& =(1-\omega )E_0^N+\omega E_0^{N+1}\\ \rho \left(\tilde{N}\right)& =(1-\omega ){\rho }_0^N+\omega {\rho }_0^{N+1}\end{array}$(33) ) for a single $\tilde{N}$-electron system lack the (statistical) physical interpretation of the straight-line behaviour of Refs. [6, 7].

The representation of a fractional system by an ensemble of density matrices corresponding to different electron numbers leads to a jump in the KS potential when passing the integer N, as first noted in [6, 7]. It is known that when we represent the ground-state density of a system, ${\rho }_0^N$ say, with a KS determinant, the exact asymptotic decay of the density according to $e^{-2\sqrt{2I}r}$ should be matched by the decay of the KS density, which is that of the slowest decaying orbital density, the HOMO density, $e^{-2\sqrt{-2{ϵ}_H}r}$. From this follows ${ϵ}_H=-I$ [40–42]. If we consider the density $\rho \left(\tilde{N}\right)$, it is clear that as soon as $\omega >0$ there is a contribution of the ground state $(N+1)$-electron density to $\rho \left(\tilde{N}\right)$. One can represent the density $\rho \left(\tilde{N}\right)$ with an ensemble of two KS densities, with the same KS potential and the same orbitals, but with different occupations of the orbitals. One KS density will have weight $(1-\omega )$ and N occupied orbitals (ignoring spin), and the other one will have weight ω and $(N+1)$ occupied orbitals. So the $(N+1)$-th orbital (the former LUMO) becomes occupied with ω electrons, and the slowest decay of the ensemble density will be according to $e^{-2\sqrt{-2{\widetilde{ϵ}}_L}r}$. Here the orbital energy of the LUMO must be ${\widetilde{ϵ}}_L=-A$, at any $\omega >0$, since the ionisation energy of the negative ion, which dictates the decay of ${\rho }_0^{N+1}$, is A. Now in general the LUMO level of the N-electron system, ${ϵ}_L$, will be lower than $-A$, therefore the KS potential must shift up by a constant Δ over the atom or molecule (but not asymptotically) so that the LUMO level is raised to $-A:\mathrm{\Delta }=-A-{ϵ}_L$. If the N-electron system was open shell, then HOMO and LUMO would be at the same energy $-I$ (ignoring spin polarisation effects), and the upshift Δ of the potential in this case is I−A, equal to the postulated discontinuity I−A in the derivative of the energy at $\tilde{N}=N$, (34) ${\left.\frac{\mathrm{\partial }\tilde{E}}{\mathrm{\partial }\tilde{N}}\right|}_{\tilde{N}\downarrow N}-{\left.\frac{\mathrm{\partial }\tilde{E}}{\mathrm{\partial }\tilde{N}}\right|}_{\tilde{N}\uparrow N}=I-A=\mathrm{\Delta }$(34) The upshift is therefore also often denoted as ‘the derivative discontinuity of the potential’, even in the case it is $\mathrm{\Delta }=-A-{ϵ}_L$. If we start with a closed shell system, the LUMO level will be above ${ϵ}_H=-I$ and the upshift $\mathrm{\Delta }={\widetilde{ϵ}}_L-{ϵ}_L<I-A$. It is still often denoted as derivative discontinuity. We return to the jump behaviour in Section 4.4

4.2. Locality and size consistency

We have encountered the straight lines (SL) energy as either the average energy according to statistical mechanics of the systems in a macroscopic sample at a properly taken $T\to 0$ limit according to PPLB, or as an arbitrary behaviour for a single noninteger electron system at T=0 following from a particular construction for the energy in this nonphysical domain, see Section 4.1. We call the energy of (33(33) $\begin{array}{rl}{\hat{\mathcal{D}}}^{min}& =(1-\omega )|{\mathrm{\Psi }}_0^N⟩⟨{\mathrm{\Psi }}_0^N|+\omega |{\mathrm{\Psi }}_0^{N+1}⟩⟨{\mathrm{\Psi }}_0^{N+1}|\\ \tilde{E}\left(\tilde{N}\right)& =(1-\omega )E_0^N+\omega E_0^{N+1}\\ \rho \left(\tilde{N}\right)& =(1-\omega ){\rho }_0^N+\omega {\rho }_0^{N+1}\end{array}$(33) ) the straight-line energy $E^{SL}$. An interesting attempt to prove the linearity of such an energy (without recourse to statistical physics and a grand canonical ensemble) has been made by Yang, Zhang and Ayers (YZA) [37, 38]. These authors considered the dissociation behaviour of molecules, extending similar arguments by PPLB [6] (see next subsection), and find that linearity follows when the functional is required to be local. Since the property we call ‘locality’ is called ‘size-consistency’ by YZA, it is important to define terms here.

Size-consistency is a requirement on any correct theoretical treatment of a set of noninteracting (sub)systems (35) $\text{Size consistency:}E(A\cdots B)=E_A+E_B$(35) (the dots indicate long enough distance to make the interaction negligible). Size-consistency is a physical requirement formulated for a separation of a system into physical subsystems A and B with defined energies $E_A$ and $E_B$. It does not apply for unphysical fragments lacking a defined energy, such as fragments with a noninteger number of electrons (which are necessarily entangled with other fragments). Size consistency should be obeyed by any proper functional (any proper quantum chemical method) since the energy is an extensive property.

The property of locality of a functional is something very different. It is not about (physical) systems but about functionals. A DFA is (semi-)local if it computes the energy with an energy density $ϵ\left({\mathit{r}}_1\right)$ that uses at ${\mathit{r}}_1$ only $\rho \left({\mathit{r}}_1\right)$ or (in the semilocal case) only nearby densities, for instance when derivatives $\mathrm{\nabla }\rho \left({\mathit{r}}_1\right)$, ${\mathrm{\nabla }}^2\rho \left({\mathit{r}}_1\right),\dots$ are used. LDA, GGA and meta-GGA are examples. We define also a property we call domain-locality. Domain-locality arises when a functional does not determine the energy density at point $\mathit{r}$ from the electron density at (and in a neighbourhood of) $\mathit{r}$, but determines the energy contributions from separate (nonoverlapping) densities, in cases where the total density is built up from such disjunct pieces. A domain-local functional is for instance the $E^{SL}$ functional for a fractional electron density: with a density that is separated from the rest of the system, as for instance the atomic densities in dissociated H${}_2^{+}$, it needs integration over that local domain to find the fractional number of electrons $\tilde{N}$, which is an ingredient in the $E^{SL}$ energy determination. For a local functional (either strictly local, semi-local or domain-local), the following relation holds (the dots indicate nonoverlapping densities) (36) $\begin{array}{rl}& \text{(strict, semi--, domain--)locality:}\\ & E^{DFA}[{\rho }_A\cdots {\rho }_B]=E^{DFA}\left[{\rho }_A\right]+E^{DFA}\left[{\rho }_B\right]\end{array}$(36) So locality implies that energy values are assigned even to possibly unphysical (noninteger) densities ${\rho }_A$ and ${\rho }_B$. Nothing in the proof of the HK theorem makes us expect the HK functional to be local. The prime example of nonlocality of the exact functional occurs for a system of two open shell atoms X and Y at large separation in its singlet ground state. The $E_{xc}$ functional cannot be local since its derivative, the xc potential, must be upshifted by a constant over the atom with highest ionisation energy, the constant being determined by the ionisation energy of the other atom, however remote [7, 14, 43–51]. Another example is found when a separated fragment density ${\rho }_A$ does not correspond to a physical system, for instance when it does not integrate to an integer number of electrons. In that case the exact functional cannot be local. The simplest and very well-known case of unphysical fragments is H${}_2^{+}$ at long distance, with electron charge densities of $\left(1/2\right)e$ at each site, or an array of P protons with very large distances and 1 electron, with charge densities of $\left(1/P\right)e$ at each site. For such a one-electron system, the $E_{xc}$ energy has to cancel the Hartree term, i.e. $\int \rho \left(\mathit{r}\right){ϵ}_{xc}\left(\mathit{r}\right)\mathrm{d}\mathit{r}=-\left(1/2\right)\int \rho \left({\mathit{r}}_1\right)\rho \left({\mathit{r}}_2\right)/r_{12}\mathrm{d}{\mathit{r}}_1\mathrm{d}{\mathit{r}}_2$, so ${ϵ}_{xc}\left({\mathit{r}}_1\right)=-\left(1/2\right)\int \rho \left({\mathit{r}}_2\right)/r_{12}\mathrm{d}{\mathit{r}}_2$, manifestly nonlocal.

Now one could try to approximate the correct result, for instance for a two-site system, with a (domain-)local functional, (37) $F\left[\rho \right]\stackrel{?}{=}{F}^{local}\left[\rho \right]=F^{local}\left[\rho \right(A\left)\right]+F^{local}\left[\rho \right(B\left)\right]$(37) This local functional delivers the correct number (but only at infinite distance, outside the range of the Coulomb potential) if it is defined, for fractional electron number, to be a linear interpolation between the defined functionals for integer electron systems. This is how the $E^{SL}$ functional could be applied. First it is to be considered as a domain-local functional, which is to be applied for each of the two nonoverlapping densities ${\rho }_A$ and ${\rho }_B$ separately. These are, in the H${}_2^{+}$ example, both 1/2-electron densities. For example, ${\rho }_A$ is the average of a one-electron density ${\rho }_H\left(A\right)=|1s_H(\mathit{r}-{\mathit{R}}_A){|}^2$ and a zero electron density ${\rho }_0\left(A\right)=0$. For such a local noninteger electron density, $E^{SL}$ defines the energy as the average of the one-electron and zero electron energies, (38) $\begin{array}{rl}{F}^{SL}[{\rho }_A\cdots {\rho }_B]& =F^{SL}\left[\rho \right(A\left)\right]+F^{SL}\left[\rho \right(B\left)\right]\left(\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\right)\\ F^{SL}\left[\rho \right(A\left)\right]& =F^{SL}\left[\frac{1}{2}{\rho }_H\left(A\right)+\frac{1}{2}{\rho }_0\left(A\right)\right]\\ & \equiv \frac{1}{2}F\left[{\rho }_H\right(A\left)\right]+\frac{1}{2}F\left[{\rho }_0\right(A\left)\right]\end{array}$(38) where the first line expresses the requirement that the SL functional is to be applied as a domain-local functional. In the case of P proton sites with one delocalised electron (39) $\begin{array}{rl}{F}^{SL}\left[\rho \right(p\left)\right]& =F^{SL}\left[\frac{1}{P}{\rho }_H\left(p\right)+\frac{P-1}{P}{\rho }_0\left(p\right)\right]\\ & \equiv \frac{1}{P}F\left[{\rho }_H\right(p\left)\right]+\frac{P-1}{P}F\left[{\rho }_0\right(p\left)\right]\end{array}$(39) (with ${\rho }_H\left(p\right)=|1s_H(\mathit{r}-{\mathit{R}}_p){|}^2$ and ${\rho }_0\left(p\right)=0$). Equation(39(39) $\begin{array}{rl}{F}^{SL}\left[\rho \right(p\left)\right]& =F^{SL}\left[\frac{1}{P}{\rho }_H\left(p\right)+\frac{P-1}{P}{\rho }_0\left(p\right)\right]\\ & \equiv \frac{1}{P}F\left[{\rho }_H\right(p\left)\right]+\frac{P-1}{P}F\left[{\rho }_0\right(p\left)\right]\end{array}$(39) ), where P can be arbitrarily large, shows that the interpolation has to be linear. This does not imply that a fractional electron system is physical. It only shows that we can, in the limit of noninteracting subsystems, try to ignore the nonlocal nature of the exact functional and introduce another functional, not exact, which is taken to be local (as defined in Equation (36(36) $\begin{array}{rl}& \text{(strict, semi--, domain--)locality:}\\ & E^{DFA}[{\rho }_A\cdots {\rho }_B]=E^{DFA}\left[{\rho }_A\right]+E^{DFA}\left[{\rho }_B\right]\end{array}$(36) )) and which delivers the correct energy for the total (fully dissociated) system. We have to endow it with the following properties: (a) it accepts unphysical fractional electron densities (since such fragment densities may occur); (b) it interpolates linearly between the functionals of the nearby integer densities (Equation (39(39) $\begin{array}{rl}{F}^{SL}\left[\rho \right(p\left)\right]& =F^{SL}\left[\frac{1}{P}{\rho }_H\left(p\right)+\frac{P-1}{P}{\rho }_0\left(p\right)\right]\\ & \equiv \frac{1}{P}F\left[{\rho }_H\right(p\left)\right]+\frac{P-1}{P}F\left[{\rho }_0\right(p\left)\right]\end{array}$(39) ) with variable P shows that linearity is required [37]).

YZA [37] have demonstrated more generally (for the case of degenerate subsystems) that the property of locality of a functional (in the case of disjoint subsystems) leads to the requirement that the functional has to be defined for noninteger densities (which the noninteracting subsystems may possess). It is then deduced, as above, that the behaviour of the energy for fractional electron densities should be linear. However, we wish to caution that locality is not a property of the exact functional. The requirement of locality (36(36) $\begin{array}{rl}& \text{(strict, semi--, domain--)locality:}\\ & E^{DFA}[{\rho }_A\cdots {\rho }_B]=E^{DFA}\left[{\rho }_A\right]+E^{DFA}\left[{\rho }_B\right]\end{array}$(36) ) is called size consistency in Ref. [37], or functional size consistency. But the two should not be confused: size consistency is the property (35(35) $\text{Size consistency:}E(A\cdots B)=E_A+E_B$(35) ) (this is called energy size consistency in Ref. [37]). Size consistency is a property the exact functional will possess (as any bona fide theory must), and it should be required of approximate functionals (recall the requirement that A and B of Equation (35(35) $\text{Size consistency:}E(A\cdots B)=E_A+E_B$(35) ) be physical systems). Locality (36(36) $\begin{array}{rl}& \text{(strict, semi--, domain--)locality:}\\ & E^{DFA}[{\rho }_A\cdots {\rho }_B]=E^{DFA}\left[{\rho }_A\right]+E^{DFA}\left[{\rho }_B\right]\end{array}$(36) ) is not a property of the exact functional, and can only be a property of approximate functionals. The requirement of size-consistency, in particular in the case of degenerate subsystems, and the difficulties involved in obeying this requirement for approximate (local) functionals, have been analysed by Savin [52, 53] and Gori-Giorgi and Savin [54]. Locality is not a requirement for size consistency. The exact functional will be nonlocal and must be size consistent.

4.3. Dissociation into integer electron fragments and (non)locality

A pleasing property of the SL energy $E^{SL}$ (i.e. $\tilde{E}\left(\tilde{N}\right)$ of (33(33) $\begin{array}{rl}{\hat{\mathcal{D}}}^{min}& =(1-\omega )|{\mathrm{\Psi }}_0^N⟩⟨{\mathrm{\Psi }}_0^N|+\omega |{\mathrm{\Psi }}_0^{N+1}⟩⟨{\mathrm{\Psi }}_0^{N+1}|\\ \tilde{E}\left(\tilde{N}\right)& =(1-\omega )E_0^N+\omega E_0^{N+1}\\ \rho \left(\tilde{N}\right)& =(1-\omega ){\rho }_0^N+\omega {\rho }_0^{N+1}\end{array}$(33) )) is that it affords correct dissociation of molecules into integer electron fragments (atoms or molecular fragments), while many LDFA's yield fractionally charged fragments. This problem was originally raised by Slater with the dissociation of NaCl as example, cf. Ref. [2, Ch. 4], and was shown to be solved by a SL energy in Ref. [6]. However, it should be recognised that this correct dissociation behaviour only follows if the assumption is made that the functional is applicable in local fragments separately (domain-locality). We will review in this section how an $E^{SL}$ leads to dissociation of a diatomic into integer electron atoms. In the next subsection, we will then demonstrate that the straight line behaviour of the energy still cannot be ‘the exact functional’. The derivative discontinuity jump of the KS potential discussed at the end of Section 4.1 for the case of a single molecule, does not have the right magnitude in the case of dissociation of a heteronuclear diatomic. The right jump of the KS potential over an atom in a dissociated system has a nonlocal origin, which is why a locality Ansatz has to fail.

Let us denote the energy for fractional electron number $\tilde{E}\left(\tilde{N}\right)$, where $\tilde{N}=N+\omega ,0<\omega <1$ or $-1<\omega <0$ is the fractional electron number and $\tilde{E}$ is linear on the two intervals. When a system X−Y separates into two parts X and Y with non-overlapping densities $\rho \left(X\right)$ and $\rho \left(Y\right)$, not necessarily integer, it is assumed the energy can be obtained as the sum of subsystem energies determined by the subsystem densities: $E\left[\rho \right(X\cdots Y\left)\right]=\tilde{E}\left[\rho \right(X\left)\right]+\tilde{E}\left[\rho \right(Y\left)\right]$ (assumption of locality of the functional). If the linear $\tilde{E}$ is applied to each subsystem X and Y separately, dissociation into integer electron systems has the lowest energy. The argument runs as follows. Suppose the starting electron distribution is ‘wrong’ in the sense that $A\left(Y\right)>I\left(X\right)$ (e.g. Y =Na${}^{+}$ and X=Cl⁻). With the linearity for the energy at fractional electron number $\tilde{N}$ and constant $\mathrm{\partial }\tilde{E}/\mathrm{\partial }\tilde{N}=-I$ or $-A$ for $\tilde{N}<N$ and $\tilde{N}>N$, respectively, we will have that a fractional charge transfer $\delta \tilde{N}$ will occur from X to Y with energy change (40) $\begin{array}{rl}\mathrm{\Delta }E& =\frac{\mathrm{\partial }\tilde{E}\left[X\right]}{\mathrm{\partial }\tilde{N}}(-\delta \tilde{N})+\frac{\mathrm{\partial }\tilde{E}\left[Y\right]}{\mathrm{\partial }\tilde{N}}\left(\delta \tilde{N}\right)\\ & =I\left(X\right)\delta \tilde{N}-A\left(Y\right)\delta \tilde{N}<0\end{array}$(40) If e.g. Y =Na${}^{+}$, the $3s$ level of Na${}^{+}$ will, due to the derivative discontinuity jump of the KS potential, jump up to the $3s$ level of neutral Na, at $-I$(Na), for an arbitrary amount ω $(0<\omega <1)$ of electron transfer, while the $3p_{\sigma }$ of Cl${}^{-1+\omega }$ remains at $-I$(Cl⁻), see Figure 2(a). Since the energy derivatives remain constant at $-A\left(Y\right)$ and $-I\left(X\right)$ all along the fractional charge on X on the $(N_X,N_X-1)$ interval and on Y on the $(N_Y,N_Y+1)$ interval, transfer will continue until a complete electron has transferred from X to Y. So separation of X and Y will occur into two integer electron systems. However, this is not a stable end result for the SL energy model, see next section.

Figure 2. (a) Na${}^{+}$Cl⁻: The constant frontier orbital energies of Na${}^{+}$ and Cl⁻ for a straight-line energy during the charge equalisation process from charged fragments to neutral atoms. After the discontinuity jump of the Na${}^{+}$, 3s level has occurred at transfer of any amount ω $(0<\omega \le 1)$ electrons to the 3s level of neutral Na at $-I$(Na), it is still below the $3p_{\sigma }$ levels of Cl⁻ (α and β level indicated with up and down arrows). (b) NaCl: The dissociated system of two neutral atoms Na and Cl and the derivative discontinuity jump upon a small electron transfer δ from the Na-3$s\alpha$ orbital to the lower lying Cl-3$p_{\sigma }\alpha$ orbital. The discontinuity jump puts Cl-3$p_{\sigma }\alpha$ at $-A$(Cl)=$-I$(Cl⁻), above the Na-$3s$ level. Degeneracy of Na-3$s\alpha$ and Cl-3$p_{\sigma }\beta$ does not result.

4.4. Nonlocality and dissociation of an electron pair bond

In spite of the correct dissociation into integer electron fragments, the linear energy behaviour and the inherent derivative discontinuity of the KS potential do not lead to completely correct dissociation. The problem is that the linear energy $\tilde{E}$ has to be applied in combination with a (domain-)locality Ansatz. The systems X and Y are treated as independent subsystems to which locally the linear energy behaviour can be applied if they have noninteger electron number. The typical result of the dissociation discussed so far will be two neutral atoms ($Y=\mathrm{N}\mathrm{a}\left[\right(3s\alpha {)}^1]$, $X=\mathrm{C}\mathrm{l}\left[\right(3p_{\pi }{)}^4\left(3p_{\sigma }\beta {)}^1\right]$ in the given example). But this is not a stable situation. For the neutral atoms $Y=\mathrm{N}\mathrm{a}$ and $X=\mathrm{C}\mathrm{l}$, we now have that $I\left(X\right)>I\left(Y\right)$, see Figure 2(b). For the atom X, the presence of another open shell atom Y far away in the universe, with a smaller ionisation potential $I\left(Y\right)$ than $I\left(X\right)$, makes it necessary that the KS potential exhibits an upshift all over the region of X (but not asymptotically) by the constant $C=I\left(X\right)-I\left(Y\right)$. This upshift has been identified a long time ago [7, 43] and has received considerable interest [14, 44–51]. Its origin is in the response part of the KS potential [45]. That is understandable: the conditional amplitude that is underlying the response potential describes the strong correlation in this case: when one electron is on Y the other electron stays away on X and does not (not even partially) delocalise towards Y. The potential step following from the conditional amplitude arises because of the long-range correlation, a manifestly nonlocal effect. It is obvious that this step of the KS potential over atom X requires a nonlocal functional – the presence of another atom $\left(Y\right)$ very far away and the magnitude of its ionisation energy $I\left(Y\right)$ cannot be described with only the local density $\rho \left(X\right)$ available.

The potential jump $I\left(X\right)-I\left(Y\right)$ over atom X makes the unpaired spin levels of both atoms degenerate. A closed-shell KS molecular orbital solution can result, with doubly occupied HOMO having 50–50 mixture of X and Y character. This KS solution will have equal probability of spin α or β on Y and X, in agreement with the exact wavefunction.

This correct solution does not result with the standard derivative discontinuity jump of the KS potential. Suppose that the dissociation to neutral atoms described in Section 4.3 has been completed, see Figure 2(b). As soon as a full electron has been transferred from ${\mathrm{C}\mathrm{l}}^{-}$, the orbital energies of the generated neutral atoms revert to their ‘normal’ values ($-I$(Cl) for the $3p_{\sigma }\beta$ on Cl, $-I$(Na) for $3s\alpha$ on Na). So the empty Cl-$3p_{\sigma }\alpha$ will be below the occupied Na-$3s\alpha$. Suppose that on a further iteration in an SCF calculation some number δ of α electrons is transferred from the higher lying Na-$3s\alpha$ to the lower lying unoccupied Cl-$3p_{\sigma }\alpha$ (it could be 1 electron if the transfer is not damped). The ‘derivative discontinuity’ jump of the KS potential will put Cl-$3p_{\sigma }\alpha$ at $-A$(Cl) for δ electrons. But this will typically be a larger upshift than the required I(Cl)-I(Na), so this is not a stable situation: the electrons will on the following iteration be sent back to the now lower lying Na-$3s\alpha$ level. If actually carried out, this calculation will lead to infinite oscillation rather than a self-consistent stable solution because the desired situation with 50–50 mixing of X(=Cl) and Y (=Na) orbitals has to come from perfect degeneracy of X and Y orbitals, which is not achieved if the KS potential plateau of $C=I\left(X\right)-I\left(Y\right)$ over atom X is not generated. And it is not generated because the $\tilde{E}\left(\tilde{N}\right)$ functional lacks the non-locality (i.e. knowledge of Y ) that is required to determine the right upshift of the KS potential over the X atom. The ‘derivative discontinuity’ jump of the KS potential over the X atom of magnitude $I\left(X\right)-A\left(X\right)$ is too large and prevents a stable proper dissociation situation, see Figure 2(b). The jump should not be $I\left(X\right)-A\left(X\right)$ (the locally determined derivative discontinuity jump) but $I\left(X\right)-I\left(Y\right)$ (dictated nonlocally) in order to obtain proper delocalisation of the α and β electrons over X and Y.

The problems discussed in this section and Section 4.2 arise from attempts to treat dissociated systems with neglect of the fact they are still entangled. Entanglement does not show up (or very weakly) in the local densities, but it does in the two-particle density matrix (pair density, conditional amplitude and conditional probabilities, for instance), and manifests itself in the nonlocality of the functional. Nonlocal functionals leading to proper dissociation can be formulated, but they typically require orbital dependency. Refs. [55–58] give examples in a KS context, and in a DMFT context they are commonplace [59–64]. A genuine density functional treatment of strongly correlated electrons, which is exact in the limiting case of dissociated H${}_2$, is provided by the strictly correlated electrons (SCE) functional [65–68].

5. Meaning and applications of the relation $\mathrm{\partial }E/\mathrm{\partial }{n}_i={ϵ}_i$ for Hartree–Fock, Xα and approximate density functional energies

While orbital occupation numbers do not have a place in exact Kohn–Sham DFT, in particular not when noninteger values are considered, this is different in many approximate total energy expressions, corresponding to approximate electronic structure models such as Hartree–Fock, or density functional approximations such as Xα, LDA, GGAs (LDFAs in general) and hybrid functionals. Even if fractional occupation numbers do not have a physical meaning, they still yield mathematically defined energy values, which allow mathematical manipulations to extract physically meaningful quantities. An example is provided by the use of occupation numbers in the Hartree–Fock energy, as originally done by Slater [1, 2, 69],

(41) $\begin{array}{rl}{E}^{HF1}& =\sum _{\begin{array}{c}\sigma =\alpha ,\beta \\ p\end{array}}{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ p,q\end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(41)

We use the more general case of unrestricted Hartree–Fock. The summations over σ and τ run over the spin functions $(\alpha ,\beta )$. The Hartree–Fock model is defined as the lowest energy determinantal wavefunction. Orbitals occur in the determinant or not. Fractional occupation numbers have no meaning in the theory, an orbital cannot be partly in the determinant. The introduction of occupation numbers can be seen as a convenient notational device, the occurrence of a spinorbital being determined by the occupation number being 0 or 1. The summations are over all orbitals (we use the conventional notation of indexing occupied orbitals with $i,j,k,l,\dots$, unoccupied orbitals with $a,b,\dots$ and general orbitals with $p,q,r,\dots$). The presence or the absence of an orbital is governed by $n_p=1$ or 0. In order to indicate that we extended the Hartree–Fock energy expression with occupation numbers as additional variables, for which we choose linear dependency (power 1), we denote this energy expression as ‘HF1’. Alternative expressions with occupation numbers are possible, see Equation (46(46) $\begin{array}{rl}{E}^{HF2}& =\sum _{i,\sigma }(n_i^{\sigma }{)}^2⟨{\phi }_i^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_i^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ i,j\end{array}}(n_i^{\sigma }{n}_j^{\tau }{)}^2\\ & \times \int \frac{{\phi }_i^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_j^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_i^{\sigma }\left({\mathit{x}}_1\right){\phi }_j^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(46) ). The energy is usually optimised under orbital variation in $E^{HF1}$ with fixed occupation numbers, where the orthonormality of the orbitals is treated with Lagrange multipliers. Treating the occupation numbers as variables, the optimisation has to be carried out under the constraints $\{n_i^{\alpha }=1,i\le H^{\alpha };n_a^{\alpha }=0,a>H^{\alpha }\},\{n_j^{\beta }=1,j\le H^{\beta };n_b^{\beta }=0,b>H^{\beta }\}$, where $H^{\alpha }$ and $H^{\beta }$ are the highest occupied α and β spin orbitals, respectively. It is elementary to derive from Equation (41(41) $\begin{array}{rl}{E}^{HF1}& =\sum _{\begin{array}{c}\sigma =\alpha ,\beta \\ p\end{array}}{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ p,q\end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(41) ) the Slater relation (42) $\frac{\mathrm{\partial }{E}^{HF1}}{\mathrm{\partial }{n}_p^{\sigma }}={ϵ}_p^{\sigma }$(42) The constraints on the occupation numbers can be treated with Lagrange multipliers because Equation (41(41) $\begin{array}{rl}{E}^{HF1}& =\sum _{\begin{array}{c}\sigma =\alpha ,\beta \\ p\end{array}}{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ p,q\end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(41) ) defines the objective function $E^{HF1}$ also outside the feasible values of 0 and 1 for the occupation numbers (in contrast to the energy according to HK functional of the density). The Hartree–Fock model itself and the Hartree–Fock energy have no physical meaning outside the feasible values $n_i^{\sigma }=1,n_a^{\sigma }=0$. The choice for the value of the objective function $E^{HF}$ outside the feasible domain could be made differently than in $E^{HF1}$, see below. Optimisation with Lagrange multipliers for the constraints (including the well-known constraints of orbital orthonormalisation) requires stationarity of the Lagrangean (43) $\begin{array}{rl}{L}^{HF1}& =E^{HF1}-\sum _{\sigma }\left(\sum _{p,q}{ϵ}_{pq}^{\sigma }\left(⟨{\phi }_p^{\sigma }|{\phi }_q^{\sigma }⟩-{\delta }_{pq}\right)\right.\\ & \left.-\sum _{i\le H^{\sigma }}{\lambda }_i^{\sigma }(n_i^{\sigma }-1)-\sum _{a>H^{\sigma }}{\lambda }_a^{\sigma }(n_a^{\sigma }-0)\right)\end{array}$(43) The values of the Lagrange multipliers for the occupation number constraints follow from the conditions $\mathrm{\partial }{L}^{HF1}/\mathrm{\partial }{n}_p^{\sigma }=0$, (44) $\mathrm{\partial }{L}^{HF1}/\mathrm{\partial }{n}_p^{\sigma }=0\Rightarrow {\lambda }_p^{\sigma }=\frac{\mathrm{\partial }{E}^{HF1}}{\mathrm{\partial }{n}_p^{\sigma }}$(44) The force of constraint, needed to keep the occupation number at its prescribed value, is proportional to the corresponding Lagrange multiplier ${\lambda }_p^{\sigma }$ (in this case equal to it, since $\mathrm{\partial }(n_p^{\sigma }-1)/\mathrm{\partial }{n}_p^{\sigma }=\mathrm{\partial }(n_p^{\sigma }-0)/\mathrm{\partial }{n}_p^{\sigma }=1$). From Equations (44(44) $\mathrm{\partial }{L}^{HF1}/\mathrm{\partial }{n}_p^{\sigma }=0\Rightarrow {\lambda }_p^{\sigma }=\frac{\mathrm{\partial }{E}^{HF1}}{\mathrm{\partial }{n}_p^{\sigma }}$(44) ) and (42(42) $\frac{\mathrm{\partial }{E}^{HF1}}{\mathrm{\partial }{n}_p^{\sigma }}={ϵ}_p^{\sigma }$(42) ) follows that (45) ${\lambda }_p^{\sigma }={ϵ}_p^{\sigma }$(45) The fact that the Lagrange multiplier for the constraint of orbital normalisation, ${ϵ}_p^{\sigma }$, is equal in magnitude to the force of constraint for constant occupation number, makes physically sense: regardless of the fact that increase of $n_p^{\sigma }$ to a value >1 does not have physical meaning in the Hartree–Fock model, such increase in (41(41) $\begin{array}{rl}{E}^{HF1}& =\sum _{\begin{array}{c}\sigma =\alpha ,\beta \\ p\end{array}}{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ p,q\end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(41) ) would lower $E^{HF1}$, in particular for large negative ${ϵ}_p^{\sigma }$ (core orbitals).

We may recall [3] that a different choice could be made for the dependence of the energy on occupation numbers. Since the value of the objective function, $E^{HF}$ outside the domain to which the occupation numbers are constrained is arbitrary, one can choose the function behaviour there at will. Of course the force of constraint at the optimum point will depend on the chosen behaviour outside the domain, see Section 2. We can illustrate this by using for Hartree–Fock the energy (46) $\begin{array}{rl}{E}^{HF2}& =\sum _{i,\sigma }(n_i^{\sigma }{)}^2⟨{\phi }_i^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_i^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ i,j\end{array}}(n_i^{\sigma }{n}_j^{\tau }{)}^2\\ & \times \int \frac{{\phi }_i^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_j^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_i^{\sigma }\left({\mathit{x}}_1\right){\phi }_j^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(46) The total energies $E^{HF1}$ and $E^{HF2}$ are exactly the same at the values $n_i^{\sigma }=1$ and $n_a^{\sigma }=0$ (i.e. on the feasible domain), only the behaviour outside that domain differs. Optimisation with the same constraint on the occupation numbers being 0 and 1 leads to exactly the same orbitals and orbital energies and the same total energy. Of course, the derivatives with respect to $n_p^{\sigma }$ are different (47) $\mathrm{\partial }{E}^{HF2}/\mathrm{\partial }{n}_p^{\sigma }=2n_p^{\sigma }{ϵ}_p^{\sigma }$(47) Obviously, the force of constraint ${\lambda }_i^{\sigma }=2n_i^{\sigma }{ϵ}_i^{\sigma }$ for occupied orbitals (to keep $n_i^{\sigma }=1$) is larger, since the energy now changes with the square of the occupation number (in the unphysical domain outside $n_i^{\sigma }=1$). The derivative is now larger and then a larger force of constraint is needed to cancel it. The unoccupied orbitals on the other hand now have zero force of constraint since the required zero occupation number is at the minimum of the quadratic $(n_a^{\sigma }{)}^2$ behaviour for $a\in unocc$.

It has become popular to move to hybrid functionals with a certain percentage of exact (Hartree–Fock) exchange, in addition to a pure density functional for correlation and the rest of the exchange. When such an energy expression is optimised by orbital variation, as in the traditional Hartree–Fock approach, one will obtain a partly nonlocal potential (a fraction of the Hartree–Fock exchange operator). This procedure is called Hartree–Fock–Kohn–Sham [26] or generalised Kohn–Sham (gKS) [70, 71]. Suppose now that we enter occupation numbers in the exact-exchange part of the energy in the same way as in Hartree–Fock, and also write the electron density and the kinetic energy and electron-nuclear energies as linear expressions in the occupation numbers, all as in the energy $E^{HF1}$. It is elementary to show that in that case the Slater relation $\mathrm{\partial }{E}^{hybrid1}/\mathrm{\partial }{n}_i={ϵ}_i$ also holds (we use again $hybrid1$ to signal that linear dependence on occupation numbers has been introduced). We recall that the derivation of the Slater relation only uses the chain rule for taking the derivative of the density-dependent part $E_{xc}^{DFA}\left[\rho \right]$ and, if present, of the one-particle reduced density matrix (1RDM) part $E_{xc}^{DFA}\left[\gamma \right]$ of the total energy, cf. [3]. This relation has also been called the generalised Janak theorem [71, 72]. Assume further that a ground-state calculation is performed, constraining the occupation numbers to 1 and 0. [If we would relieve that constraint and allow fractional occupation numbers, the functional would no longer be distinct from a functional of the one-electron reduced density matrix (1RDM). The exact 1RDM functional has been shown by Gilbert [36] to have complete degeneracy of the orbital energies for all partially occupied orbitals.] If we now consider the derivative with respect to the occupation number of the LUMO, we are dealing with an infinitesimal increase of the total number of electrons. This derivative with respect to occupation number then appears to be the same as the derivative with respect to total electron number, $\mathrm{\partial }{E}^{hybrid1}/\mathrm{\partial }{n}_L=(\mathrm{\partial }E/\mathrm{\partial }N{)}_{+}$ [71, 73]. With the straight-line energy, the latter is $-A$. This has led to the conclusion that the LUMO orbital energy for an accurate hybrid functional must be close to $-A$. However, this is questionable. In the first place, the derivative $\mathrm{\partial }E/\mathrm{\partial }N$ is not defined in exact DFT, as discussed in Section 2. It is only $-A$ with a particular choice of the energy behaviour for noninteger N, namely the straight-line energy $\tilde{E}\left(\tilde{N}\right)$, Equation (33(33) $\begin{array}{rl}{\hat{\mathcal{D}}}^{min}& =(1-\omega )|{\mathrm{\Psi }}_0^N⟩⟨{\mathrm{\Psi }}_0^N|+\omega |{\mathrm{\Psi }}_0^{N+1}⟩⟨{\mathrm{\Psi }}_0^{N+1}|\\ \tilde{E}\left(\tilde{N}\right)& =(1-\omega )E_0^N+\omega E_0^{N+1}\\ \rho \left(\tilde{N}\right)& =(1-\omega ){\rho }_0^N+\omega {\rho }_0^{N+1}\end{array}$(33) ). In the second place, the derivative with respect to occupation number depends on the way occupation numbers are introduced in a DFA. The energy at fractional occupation numbers is a mathematical device, it does not correspond to a true physical system and therefore using it to infer properties of the actual physical system (at integer electron number) is questionable. A LUMO close to $-A$ is typical for the Hartree–Fock model and can for a hybrid be obtained with a very high percentage of exact exchange [15]. For common hybrids, with exact exchange percentages in the order of 20–30%, the LUMO is much below $-A$ [3]. It is remarkable that in solids (but not in molecules) the ${ϵ}_L^{DFA}\left(solid\right)$ is equal to the DFA calculated (total energy difference) $-A^{DFA}\left(solid\right)$ for DFAs that obey the Slater relation [72]. This difference between solids and molecules is not self-evident. It requires for its proof that the ‘solid state limit’ to arbitrarily large size of the system can be taken. This is discussed at some length in Ref. [3], see also Ref. [72].

6. Summary

We have applied the theory of functional derivatives with constraints and the Lagrange multiplier technique to the problem of optimisation of the density in the Hohenberg–Kohn functional with the $\int \rho \mathrm{d}\mathit{r}=N$ constraint. We have highlighted the arbitrariness of a part of the total functional derivative, namely the derivative of the total energy with respect to electron number, $\mathrm{\partial }E/\mathrm{\partial }N$. It stems from the restricted domain on which $F\left[\rho \right]$ is defined. The Hohenberg–Kohn functional $F\left[\rho \right]$ (and the Levy–Lieb and Lieb functionals) are only defined on the domain of N-electron densities. As a corollary, Janak's theorem does not exist in (exact) Kohn–Sham DFT. The essential difficulty is with the nonexistence of fractional electron systems. Any change of an orbital occupation number from integer value (1 or 0) entails a change of electron number from the integer value. There are no wavefunctions for such systems, not for interacting electrons and not for the noninteracting electrons of the Kohn–Sham system. Therefore, there are no exact energies and no exchange–correlation energies $E_{xc}$ for noninteger N (at least not for single, finite electron systems (atoms and molecules) at T=0).

For application of the Euler–Lagrange equation to find the optimum N-electron density ${\rho }_0^N$ for which the Hohenberg–Kohn energy functional $E_v\left[\rho \right]$ attains its minimum value, the arbitrariness of $\mathrm{\partial }E/\mathrm{\partial }N$ can be solved in two ways. The functional derivative $\delta E_v\left[\rho \right]/\delta \rho \left(\mathit{r}\right)$ should be taken with the constraint $\int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=N$ (N integer). The constraint derivative $\delta E_v\left[\rho \right]/{\delta }_N\rho \left(\mathit{r}\right)$ should be zero at the ground-state density ${\rho }_0^N$. Alternative one can make a suitable choice of the $\mathrm{\partial }E/\mathrm{\partial }N$ derivative and apply the Lagrange multiplier technique. The choice should be for a continuous and finite derivative, otherwise it would impede the application of the Lagrange multiplier technique for density optimisation. Choices in the literature, such as $F_L\left[\rho \right]=+\mathrm{\infty }$ in the functional analysis literature of DFT [8–11, 24] or a discontinuous derivative (different derivatives to the electron-rich and electron-deficient sides of the integer), are not suitable in this sense.

It is important that we are dealing here with single quantum systems at T=0 (for which the majority of DFT calculations in chemistry are being done). Perdew et al. [6] have studied the statistical mechanics of systems with fluctuating electron number with a constrained search over ensembles of integer electron density matrices, leading to linear energy behaviour in between integers for the average energy at average electron number. It is an altogether different matter to define the energy of a single molecule with a noninteger electron number at T=0 by a mixture of density matrices of ground states with different electron numbers, see Equations (32(32) $\begin{array}{rl}{E}_v\left[\rho \right]& =\int \rho \left(\mathit{r}\right)v\left(\mathit{r}\right)\mathrm{d}\mathit{r}+F\left[\rho \right]\\ F\left[\rho \right]& \equiv \underset{\hat{\mathcal{D}}\to \rho }{min}Tr\left[\hat{\mathcal{D}}\right(\hat{T}+\hat{W}\left)\right]\\ \hat{\mathcal{D}}& =\sum _i{\lambda }_i|{\mathrm{\Psi }}_i⟩⟨{\mathrm{\Psi }}_i|,\sum _i{\lambda }_i=1,\\ Tr\left(\hat{\mathcal{D}}\hat{\rho }\right(\mathit{r}\left)\right)& =\sum _i{\lambda }_i{\rho }_i\left(\mathit{r}\right)=\rho \left(\mathit{r}\right),\\ \int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}& =\tilde{N}=\sum _i{\lambda }_i{N}_i,\end{array}$(32) ) and (33(33) $\begin{array}{rl}{\hat{\mathcal{D}}}^{min}& =(1-\omega )|{\mathrm{\Psi }}_0^N⟩⟨{\mathrm{\Psi }}_0^N|+\omega |{\mathrm{\Psi }}_0^{N+1}⟩⟨{\mathrm{\Psi }}_0^{N+1}|\\ \tilde{E}\left(\tilde{N}\right)& =(1-\omega )E_0^N+\omega E_0^{N+1}\\ \rho \left(\tilde{N}\right)& =(1-\omega ){\rho }_0^N+\omega {\rho }_0^{N+1}\end{array}$(33) ). That procedure represents just one possible extension of the energy functional into the nonphysical fractional electron density domain. Such an extension is essentially arbitrary, cf. Section 2. The linear choice, which we have denoted $\tilde{E}\left(\tilde{N}\right)$, has the interesting property that it leads to the proper dissociation of molecules into neutral atoms, see Section 4.3. We have also emphasised that such success is achieved by treating $\tilde{E}\left(\tilde{N}\right)$ as a local functional, applicable in the local domains of the densities of the fragments separately. This at the same time prevents such a functional from being ‘exact’ because fully correct dissociation requires a nonlocal functional, see Section 4.4. In order to obtain the proper KS wavefunction for a heteronuclear dissociated diatomic molecule $X\cdots Y$ with $I\left(X\right)>I\left(Y\right)$, the derivative discontinuity of the potential, which is an inherent property of $\tilde{E}\left(\tilde{N}\right)$, and which would be $I\left(X\right)-A\left(X\right)$ on atom X (i.e. only dependent on X), is wrong; because of nonlocality effects it should be $I\left(X\right)-I\left(Y\right)$, i.e. be dependent on both X and Y.

Derivatives of the energy with respect to occupation numbers do not exist in exact KS DFT. However, if occupation numbers are introduced as additional variables in the total energy such derivatives can be defined for approximate energy expressions, such as $E^{HF1}$ (41(41) $\begin{array}{rl}{E}^{HF1}& =\sum _{\begin{array}{c}\sigma =\alpha ,\beta \\ p\end{array}}{n}_p^{\sigma }⟨{\phi }_p^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_p^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ p,q\end{array}}{n}_p^{\sigma }{n}_q^{\tau }\\ & \times \int \frac{{\phi }_p^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_q^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_p^{\sigma }\left({\mathit{x}}_1\right){\phi }_q^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(41) ) and $E^{HF2}$ (46(46) $\begin{array}{rl}{E}^{HF2}& =\sum _{i,\sigma }(n_i^{\sigma }{)}^2⟨{\phi }_i^{\sigma }|-\frac{1}{2}{\mathrm{\nabla }}^2+v_{ext}|{\phi }_i^{\sigma }⟩\\ & +\frac{1}{2}\sum _{\begin{array}{c}\sigma ,\tau =\alpha ,\beta \\ i,j\end{array}}(n_i^{\sigma }{n}_j^{\tau }{)}^2\\ & \times \int \frac{{\phi }_i^{\sigma \star }\left({\mathit{x}}_1\right){\phi }_j^{\tau \star }\left({\mathit{x}}_2\right)(1-P_{12}){\phi }_i^{\sigma }\left({\mathit{x}}_1\right){\phi }_j^{\tau }\left({\mathit{x}}_2\right)}{|{\mathit{r}}_1-{\mathit{r}}_2|}\\ & \times \mathrm{d}{\mathit{x}}_1\mathrm{d}{\mathit{x}}_2\end{array}$(46) ), and similarly for Xα, LDA, GGA and hybrid functionals. At fractional values of these occupation numbers, one obtains mathematically defined energies, although without physical meaning. For such energies, Slater had already obtained the relation between derivative with respect to orbital occupation number and orbital energy which bears his name and which Janak has tried to extend to (exact) DFT. The Slater relation is particularly useful when total energies are hard to calculate but orbital energies are readily available [3]. These results do not imply that an extension to exact Kohn–Sham DFT is possible. In order to avoid confusion, application of the Slater relation with the approximate DFT energy expressions mentioned above should be denoted as such, and not be called application of Janak's theorem. Janak's theorem is an attempt to extend Slater's relation to exact Kohn–Sham DFT, which has been questioned in Section 3.

Acknowledgements

I am grateful to Paola Gori-Giorgi and Andreas Savin for helpful discussions and keen interest.

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix 1.

N-conserving functional derivatives

A few aspects of constraints on functional derivatives are reviewed. We refer to Gál [74, 75] for detailed discussion. In the present case, the constraint of conservation of particle number, $\int \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=N$ is specifically addressed. We follow Parr and Bartolotti [19] in separating the N-dependency and shape dependency of the density by writing it as N times a shape function $\sigma \left(\mathit{r}\right)$ which integrates to 1, (48) $\begin{array}{rl}\rho \left(\mathit{r}\right)& =N\frac{g\left(\mathit{r}\right)}{\int g\left({\mathit{r}}^{\prime }\right)\mathrm{d}{\mathit{r}}^{\prime }}=N\sigma \left(\mathit{r}\right);\int \sigma \left(\mathit{r}\right)\mathrm{d}\mathit{r}=1\\ \delta \rho \left(\mathit{r}\right)& =N\delta \sigma \left(\mathit{r}\right)+\sigma \left(\mathit{r}\right)\delta N={\delta }_N\rho \left(\mathit{r}\right)+{\delta }_{\sigma }\rho \left(\mathit{r}\right)\\ \int {\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}& =0,\int {\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=\int \sigma \left(\mathit{r}\right)\mathrm{d}\mathit{r}\delta N=\delta N\end{array}$(48) where $g\left(\mathit{r}\right)$ is defined up to a scaling constant; it can be $g\left(\mathit{r}\right)=\rho \left(\mathit{r}\right)$. Subscripts N and σ denote constant electron number N and constant shape function $\sigma \left(\mathit{r}\right)$, respectively. The constrained derivatives for functionals are analogous to partial derivatives for functions. $\delta \sigma$ should be orthogonal to σ otherwise the $N\delta \sigma$ term contains part of the electron number variation, which should follow from $\sigma \delta N$ (if $\delta \sigma$ has a component $a\sigma$ along σ, $\int N\delta \sigma \mathrm{d}\mathit{r}=\int Na\sigma \mathrm{d}\mathit{r}=Na$).

Define the changes ${\delta }_{N}E\left[\rho \right]$ of the functional $E\left[\rho \right]$ as those variations of the energy that are induced by N-conserving density changes ${\delta }_N\rho \left(\mathit{r}\right)$, $\int \left(\rho \right(\mathit{r})+{\delta }_N\rho (\mathit{r}\left)\right)\mathrm{d}\mathit{r}=N$, (A2) ${\delta }_{N}E\left[\rho \right]\equiv \int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}.$(A2) Application of the Lagrange multiplier technique requires variation of E to be defined for arbitrary variations $\delta \rho$, so if we want to break down the total derivative $\delta E/\delta \rho$ into norm-conserving and shape-conserving parts, these have to be defined as operators acting on arbitrary $\delta \rho$, (A3) $\int \frac{\delta E}{\delta \rho }\delta \rho \mathrm{d}\mathit{r}=\int \left(\frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}\right)\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}$(A3) The constraint derivative $\delta E\left[\rho \right]/{\delta }_N\rho \left(\mathit{r}\right)$ is defined as that operator (kernel of an integral operator) that delivers ${\delta }_{N}E$ from a general $\delta \rho$, by picking out the ${\delta }_N\rho \left(\mathit{r}\right)$ part of $\delta \rho \left(\mathit{r}\right)$ (51) $\begin{array}{rl}{\delta }_{N}E& =\int \frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=\frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}\left({\delta }_N\rho \right(\mathit{r})+{\delta }_{\sigma }\rho (\mathit{r}\left)\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+\int \frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+0.\end{array}$(51) The fact that the last term on the second line is zero is an essential element in the definition of $\delta E/{\delta }_N\rho$ – the total ${\delta }_{N}E$ stems from and is determined only by the N-conserving part of $\delta \rho$. The shape-conserving part, which only allows a change of N, cannot contribute to ${\delta }_{N}E$. This is an important point. In the literature [19, 26–28], it has been an issue that $\delta E\left[\rho \right]/{\delta }_N\rho \left(\mathit{r}\right)$ would be ambiguous, containing an undefined constant. One must then be careful to note that this stems from restricting the domain on which the derivative is defined to the set of N-conserving density changes ${\delta }_N\rho \left(\mathit{r}\right)$ and the observation that $\int C{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$ [9, 20]. However, in the application of the Lagrange multiplier method (see text), we need derivatives (operators) that are defined for action on an arbitrary $\delta \rho$, not just on the domain of N-conserving ${\delta }_N\rho$. So ‘orthogonality’ of $\delta E/{\delta }_N\rho \left(\mathit{r}\right)$ on ${\delta }_{\sigma }\rho \left(\mathit{r}\right)$, $\int \left(\delta E/{\delta }_N\rho \right(\mathit{r}\left)\right){\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$, is a necessary part of the definition of $\delta E/{\delta }_N\rho \left(\mathit{r}\right)$. This property would be destroyed by an additional constant $C:\int C{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=C\delta N$. A constant would give $\delta E/{\delta }_N\rho \left(\mathit{r}\right)$ an unwarranted component in the space ‘orthogonal’ to the set $\left\{{\delta }_N\rho \right(\mathit{r}\left)\right\}$. This is analogous to adding to ${\mathbf{\nabla }}_{//}f$ in Figure 1 a component $C{\mathit{n}}_{\mathrm{\perp }}$ in the direction ${\mathit{n}}_{\mathrm{\perp }}$ perpendicular to ${\mathbf{\nabla }}_{//}f$, so it would clutter $df=({\mathbf{\nabla }}_{//}f+C{\mathit{n}}_{\mathrm{\perp }})¢erdot\delta \mathit{l}$ for a general displacement $\delta \mathit{l}$ with contributions $C{\mathit{n}}_{\mathrm{\perp }}¢erdot\delta {\mathit{l}}_{\mathrm{\perp }}$ from the perpendicular component $\delta {\mathit{l}}_{\mathrm{\perp }}$.

In the same way, for a complete definition of $\delta E/{\delta }_{\sigma }\rho$ we must have $\int \left(\delta E/{\delta }_{\sigma }\rho \right(\mathit{r}\left)\right){\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$. Now for arbitrary $\delta \rho$ (A5) $\begin{array}{rl}& \int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+\int \frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}{\delta }_{\sigma }\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}+\int \frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}\\ & =\int \left(\frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}\right)\delta \rho \left(\mathit{r}\right)\mathrm{d}\mathit{r},\end{array}$(A5) consistent with the identity (A6) $\frac{\delta E}{\delta \rho \left(\mathit{r}\right)}=\frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}.$(A6) Now a derivative with constant shape function σ can only involve a change of N, so we expect (A7) $\frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}=\frac{\mathrm{\partial }E}{\mathrm{\partial }N}$(A7) and with the chain rule (A8) $\frac{\mathrm{\partial }E}{\mathrm{\partial }N}=\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\frac{\mathrm{\partial }\rho \left(\mathit{r}\right)}{\mathrm{\partial }N}\mathrm{d}\mathit{r}=\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r}.$(A8) A derivation can be given for $\delta E\left[\rho \right]/{\delta }_N\rho \left(\mathit{r}\right)$ employing the factorisation of $\rho \left(\mathit{r}\right)$ as $Ng\left(\mathit{r}\right)/\int g\left({\mathit{r}}^{\prime }\right)\mathrm{d}{\mathit{r}}^{\prime }$ (see [20]) which leads to (A9) $\frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}=\frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}-\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r},$(A9) in perfect agreement with Equation (A8(A8) $\frac{\mathrm{\partial }E}{\mathrm{\partial }N}=\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\frac{\mathrm{\partial }\rho \left(\mathit{r}\right)}{\mathrm{\partial }N}\mathrm{d}\mathit{r}=\int \frac{\delta E}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r}.$(A8) ), so that with Equations (A6(A6) $\frac{\delta E}{\delta \rho \left(\mathit{r}\right)}=\frac{\delta E}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\delta E}{{\delta }_{\sigma }\rho \left(\mathit{r}\right)}.$(A6) )–(A9(A9) $\frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}=\frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}-\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\frac{\rho \left(\mathit{r}\right)}{N}\mathrm{d}\mathit{r},$(A9) ) (A10) $\frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}=\frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\mathrm{\partial }E}{\mathrm{\partial }N}.$(A10) Note that indeed $\int \left(\delta E/{\delta }_{\sigma }\rho \right(\mathit{r}\left)\right){\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=\int \left(\mathrm{\partial }E/\mathrm{\partial }N\right){\delta }_N\rho \left(\mathit{r}\right)\mathrm{d}\mathit{r}=\left(\mathrm{\partial }E/\mathrm{\partial }N\right)N\int \delta \sigma \left(\mathit{r}\right)\mathrm{d}\mathit{r}=0$. The fact that $\mathrm{\partial }E/\mathrm{\partial }N$ is not defined by the Hohenberg–Kohn theorem, is not a problem in DFT, since only the energy changes ${\delta }_{N}E\left[\rho \right]$ are required. The fact that any constant C has no effect when acting on a feasible ${\delta }_N\rho \left(\mathit{r}\right)$ makes sure that the essential change ${\delta }_{N}E\left[\rho \right]$ is defined whatever the (chosen) value of $\mathrm{\partial }E/\mathrm{\partial }N$: (A11) $\begin{array}{rl}\int \frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)& =\int \left(\frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}+\frac{\mathrm{\partial }E}{\mathrm{\partial }N}\right){\delta }_N\rho \left(\mathit{r}\right)\\ & =\int \frac{\delta E\left[\rho \right]}{{\delta }_N\rho \left(\mathit{r}\right)}{\delta }_N\rho \left(\mathit{r}\right)={\delta }_{N}E\left[\rho \right].\end{array}$(A11)

Appendix 2.

An arbitrary constant in the constrained derivative?

Starting with the contribution by Parr and Bartolotti [19] remarks have been made in the literature that the constrained derivative (A12) ${\left[\frac{\delta E\left[\rho \right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N$(A12) which we have written for the sake of concise notation as $\delta E\left[\rho \right]/{\delta }_N\rho \left(\mathit{r}\right)$, following Gál [20], would only be defined up to an arbitrary constant. The origin is the same pitfall we mentioned in Section 2. Parr and Yang [26, p. 83] write Equation (11(11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) ) with an N-conserving derivative (A13) $\frac{\delta E\left[\rho \right]}{\delta g\left(\mathit{r}\right)}=\int {\left[\frac{\delta E\left[\rho \right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(A13) where $\left[\delta E\right[\rho \left]/\delta \rho \right(\mathit{r}){]}_N$ is ‘a functional differentiation with N fixed’ (what we denote $\delta E/{\delta }_N\rho \left(\mathit{r}\right)$). Application of the chain rule means that Equation (A13(A13) $\frac{\delta E\left[\rho \right]}{\delta g\left(\mathit{r}\right)}=\int {\left[\frac{\delta E\left[\rho \right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(A13) ) should in principle be written with the full derivative $\delta E\left[\rho \right]/\delta \rho \left({\mathit{r}}^{\prime }\right)$, as has been done in Equation (11(11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) ), but since $\delta \rho \left({\mathit{r}}^{\prime }\right)/\delta g\left(\mathit{r}\right)$ will deliver N-conserving variations ${\delta }_N\rho \left({\mathit{r}}^{\prime }\right)$, the use of the constraint derivative in this case is also correct. Using Equations (12(12) $\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}=\frac{N}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\left[\delta (\mathit{r}-{\mathit{r}}^{\prime })-\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\right]$(12) ) and (A13(A13) $\frac{\delta E\left[\rho \right]}{\delta g\left(\mathit{r}\right)}=\int {\left[\frac{\delta E\left[\rho \right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(A13) ) at the solution point $\left({\rho }_0^N\right)$, one then obtains Equation (4.4.7) of [26, p. 83], (A14) ${\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N=\int {\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(A14) from which it was concluded that $\left[\delta E\right[{\rho }_0^N\left]/\delta \rho \right(\mathit{r}){]}_N$ must be an undefined constant (the constant cannot be fixed from (A14(A14) ${\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N=\int {\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(A14) )). However, we have seen in Section 2 that at the solution point this N-conserving derivative is not an arbitrary constant (the pitfall mentioned there) but is zero. There is not a problem with Equation (A14(A14) ${\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N=\int {\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(A14) ), it is correct, but one should be careful with conclusions drawn from it. It is possible to conclude that the N-conserving derivative at the solution point is a constant, but it may be that on other grounds it is known that this constant is not arbitrary but has a specific value. The constrained derivative has been found to be zero at the solution point ${\rho }_0^N$, see Section 2. That is of course perfectly consistent with (A14(A14) ${\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N=\int {\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(A14) ). If on the other hand the full derivative would have been used in Equation (A13(A13) $\frac{\delta E\left[\rho \right]}{\delta g\left(\mathit{r}\right)}=\int {\left[\frac{\delta E\left[\rho \right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(A13) ), as we did in Equation (11(11) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta g\left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{\delta \rho \left({\mathit{r}}^{\prime }\right)}{\delta g\left(\mathit{r}\right)}\mathrm{d}{\mathit{r}}^{\prime }=0$(11) ), this would also have emerged in Equation (A14(A14) ${\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N=\int {\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(A14) ). We have noted in Section 2 that it is this full derivative that is equal to a constant at ${\rho }_0^N$, which is again compatible with equation (A14(A14) ${\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}\right]}_N=\int {\left[\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\right]}_N\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(A14) ) (for the full derivative $\delta E\left[{\rho }_0^N\right]/\delta \rho \left(\mathit{r}\right)$ instead of the constrained $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$). As we have seen, the constant in that case is indeed arbitrary, being $\mathrm{\partial }E\left[{\rho }_0^N\right]/\mathrm{\partial }N$. This difference – the constrained derivative $\delta E\left[{\rho }_0^N\right]/{\delta }_N\rho \left(\mathit{r}\right)$ being zero according to our present analysis, see discussion below Equation (13(13) $\frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left(\mathit{r}\right)}=\int \frac{\delta E\left[{\rho }_0^N\right]}{\delta \rho \left({\mathit{r}}^{\prime }\right)}\frac{g\left({\mathit{r}}^{\prime }\right)}{\int g\left({\mathit{r}}^{″}\right)\mathrm{d}{\mathit{r}}^{″}}\mathrm{d}{\mathit{r}}^{\prime }$(13) ), but equal to an undefined constant according to Ref. [26] – goes back to Parr and Bartolotti [19], who first proposed there to be an arbitrary constant in $\left[\delta E\right[{\rho }_0^N\left]/\delta \rho \right(\mathit{r}){]}_N$, which has since often surfaced in the DFT literature.