Asymptotic behaviour of the electron density and the Kohn–Sham potential in case of a Kohn–Sham HOMO nodal plane

ABSTRACT

It is known that the asymptotic decay of the electron density $n\left(\mathbf{r}\right)$ outside a molecule is informative about its first ionisation potential I₀, $n\left(\right|\mathbf{r}|\to \infty )\sim \text{exp}(-2\sqrt{2I_0}r)$. This dictates the orbital energy of the highest occupied Kohn–Sham (KS) molecular orbital (HOMO) to be ε_H = −I₀, if the KS potential goes to zero at infinity. However, when the KS HOMO has a nodal plane, the KS density in that plane will decay as $exp(-2\sqrt{-2{ϵ}_{H-1}}r)$. Conflicting proposals exist for the KS potential: from exact exchange calculations it has been found that the KS potential approaches a positive constant in the plane, but from the assumption of isotropic decay of the exact (interacting) density, it has been concluded this constant needs to be negative. Here, we show that either (1) the exact density decays differently (according to the second ionisation potential I₁) in the HOMO nodal plane than elsewhere, and the KS potential has a regular asymptotic behaviour (going to zero everywhere) provided that ε_{H − 1} = −I₁; or (2) the density does decay like $\text{exp}(-2\sqrt{2I_0}r)$ everywhere but the KS potential exhibits strongly irregular if not divergent behaviour around (at) the nodal plane.

GRAPHICAL ABSTRACT

Keywords:

• DFT
• Kohn–Sham potential
• HOMO nodal plane
• asymptotic density

1. Introduction

In Kohn–Sham (KS) density functional theory (DFT), one of the most widely used techniques in electronic structure theory, ground-state properties are calculated via the KS system, consisting of non-interacting electrons moving in the local KS potential $v_s\left(\mathbf{r}\right)$ [1]. In principle, the KS potential $v_s\left(\mathbf{r}\right)$ ensures that the electron density $n\left(\mathbf{r}\right)$ of the non-interacting KS system is the same as that of the physical, interacting system. Exact properties of $n\left(\mathbf{r}\right)$ and $v_s\left(\mathbf{r}\right)$ have played – and continue to play – a crucial role in constructing and improving approximations.

Both the square root of the density and the KS orbitals ${\psi }_k\left(\mathbf{r}\right)$ obey Schrödinger-type equations, (1) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{eff}}\left(\mathbf{r}\right)\right)\sqrt{n\left(\mathbf{r}\right)}=-I_0\sqrt{n\left(\mathbf{r}\right)}\end{array}$$(1) (2) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{Hxc}}\left(\mathbf{r}\right)\right){\psi }_k\left(\mathbf{r}\right)={ϵ}_k{\psi }_k\left(\mathbf{r}\right),\end{array}$$(2)

where the sum of the external and the Hartree-exchange-correlation potentials constitutes the KS potential, v_s = v_ext + v_Hxc. The eigenvalue in Equation (1(1) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{eff}}\left(\mathbf{r}\right)\right)\sqrt{n\left(\mathbf{r}\right)}=-I_0\sqrt{n\left(\mathbf{r}\right)}\end{array}$$(1) ) [2,3] is the first ionisation potential, I₀ = E^{N − 1}₀ − E^N₀, and the occupied KS orbitals reproduce the density, ${\sum }_k^N{\left|{\psi }_k\left(\mathbf{r}\right)\right|}^2=n\left(\mathbf{r}\right)$.

Here, we are mainly concerned with molecules, where the external potential $v_{\mathrm{ext}}\left(\mathbf{r}\right)$ goes to zero at large distance like −Z/r, with Z representing the total charge of all nuclei and r the distance from the centre of nuclear charge. In this case, according to Equations (1(1) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{eff}}\left(\mathbf{r}\right)\right)\sqrt{n\left(\mathbf{r}\right)}=-I_0\sqrt{n\left(\mathbf{r}\right)}\end{array}$$(1) )–(2(2) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{Hxc}}\left(\mathbf{r}\right)\right){\psi }_k\left(\mathbf{r}\right)={ϵ}_k{\psi }_k\left(\mathbf{r}\right),\end{array}$$(2) ), the asymptotic (r → ∞) decay of $\sqrt{n\left(\mathbf{r}\right)}$ and ${\psi }_k\left(\mathbf{r}\right)$ is (3) $$\begin{array}{c}\hfill \sqrt{n\left(\mathbf{r}\right)}\sim e^{-\sqrt{2(I_0+v_{\mathrm{eff}}\left(\infty \right))}r}\end{array}$$(3) (4) $$\begin{array}{c}\hfill {\psi }_k\left(\mathbf{r}\right)\sim e^{-\sqrt{2(-{ϵ}_k+v_{\mathrm{Hxc}}\left(\infty \right))}r}.\end{array}$$(4) Both the effective potential $v_{\mathrm{eff}}\left(\mathbf{r}\right)$ for $\sqrt{n\left(\mathbf{r}\right)}$ and the Hartree-exchange-correlation potential $v_{\mathrm{Hxc}}\left(\mathbf{r}\right)$ had been thought, until recently, to go to zero asymptotically everywhere in space (in other words, it seemed always possible to fix the arbitrary constant in the functional derivatives for finite systems such that the effective potentials go to zero in all possible directions). In the exact KS model (and in exact generalised KS models as well [4]), the model density, which decays as the square of the highest occupied molecular orbital, should decay like the exact density, leading to the identification I₀ = −ε_H [3,5].

In case the KS highest occupied molecular orbital (HOMO) has a nodal plane (HNP) extending to infinity, a very straightforward argument was given by Wu et al. [6] that the KS potential should go to a negative constant for asymptotic points r_p → ∞ in that plane. The KS density in that plane is governed by the HOMO−1, assuming HOMO−1 does not have the same nodal plane. Wu et al. [6] made the common assumption that the exact interacting density has the same asymptotic behaviour everywhere, and they observed that then the decay of the HOMO−1 in the HNP must be equal to the decay of the total density, implying $-\sqrt{2(-{ϵ}_{H-1}+v_{\mathrm{Hxc}}(r_p\to \infty ))}{r}_p=-\sqrt{2I_0}{r}_p$, so that v_Hxc(r_p → ∞) should tend to the negative constant I₀ + ε_{H − 1} = −(ε_H − ε_{H − 1}). On the other hand, it had earlier been argued, and numerical evidence had been provided, that the optimised effective potential method for the exact exchange model (xOEP) of KS theory leads to an asymptotic constant in the HNP, but positive [7–10]. Since it was pointed out that this behaviour of the KS potential has a significant effect on the orbital energies of particularly the higher lying unoccupied orbitals, with large consequences for the excitation energies calculated with time-dependent DFT [7,8], the matter is also relevant for practical calculations.

In this work, we analyse this issue, starting from a very simple question: is it true that an asymptotic constant on the HNP in the KS potential changes the exponential decay of the HOMO−1 (and of all the other orbitals) on that plane? We will see in Section 2 that the argument of Wu et al. [6] did not consider the role of the angular part of the laplacian, which, instead, cannot be neglected.

We then turn to the key question: how does the exact density behave? We use (see Section 3) the expansion of the density in the squares of the Dyson orbitals. We analyze the energy-independent equations obeyed by the Dyson orbitals, in order to derive the behaviour of the two leading Dyson orbitals, the first with eigenvalue −I₀ and the second with eigenvalue −I₁, and of the density in and close to the plane. In Section 4, we discuss the case (Case 1) in which the KS potential can have a simple and regular behaviour (the expected isotropic −1/r asymptotics everywhere). In this case, it is necessary that the asymptotic decay of the exact interacting density be different in different directions: in the HOMO nodal plane, the density decay $n\left(\right|{\mathbf{r}}_p|\to \infty )$ should be faster (according to the second ionisation potential, $\text{exp}[-2\sqrt{2I_1}{r}_p]$) than outside the plane, and the KS orbital energy of the HOMO−1 should be equal to the second ionisation potential, −ε_{H − 1} = I₁. It depends on the properties of the system (e.g. spatial and spin symmetries of the ion states) if such a situation will occur, with for instance the prototype molecules ethylene and benzene [8] as good candidates. However, such behaviour of the density will not always occur. In Section 5, we discuss the case (Case 2) that the density has the same exponential decay, governed by I₀ = −ε_H, on the nodal plane as everywhere else (although it must then decay polynomially slower). In that case, the KS potential cannot decay uniformly like −1/r but has to exhibit very irregular behaviour near the plane (possibly even divergent behaviour in the plane) in order to impart on the KS orbitals the shapes that will make them reproduce the true density.

2. Asymptotic decay of the orbitals on a plane and their corresponding potential

The usual argument to derive the asymptotic decay of the KS orbitals for a finite system (or, in general, of a single-particle Schrödinger equation with a multiplicative potential) is based on the fact that if the potential $v_s\left(\mathbf{r}\right)$ goes to 0 when $\left|\mathbf{r}\right|\to \infty$, then asymptotically the single particle equation reads $-\frac{1}{2}{\nabla }^2\psi =ϵ\psi$ implying $\psi \left(\mathbf{r}\right)\sim e^{-\sqrt{-2ϵ}r}$. If, instead, the potential goes asymptotically to a constant v_s(∞), then we trivially obtain from the same equation $\psi \left(\mathbf{r}\right)\sim e^{-\sqrt{2(-ϵ+v_s\left(\infty \right))}r}$. This leads to the idea that if the exponential decay of $\psi \left(\mathbf{r}\right)$ is the same everywhere except on a plane, then in that plane the potential has to go asymptotically to a constant. And viceversa, one expects that an asymptotic constant in $v_s\left(\mathbf{r}\right)$ on a plane implies that the asymptotic decay of $\psi \left(\mathbf{r}\right)$ on the plane has a different exponent.

This argument, however, does not take into account the fact that the laplacian also contains angular derivatives. This angular part has a 1/r² prefactor and one may think that this makes it negligible, when r → ∞, with respect to the constant terms (the eigenvalue and, if present, the constant v_s(∞)). However, if we have a different exponential decay on the plane than elsewhere, the relative difference between the orbital on the plane and very close to it increases exponentially with r, so that the angular derivative of the orbital very close to the plane increases much faster than 1/r². This diverging behaviour needs to be compensated by $v_s\left(\mathbf{r}\right)$. We illustrate this with two very simple examples.

First of all, let us consider an example that shows that an asymptotic constant on a plane in the potential does not necessarily imply a change in the exponential decay of the orbital on the plane. Suppose we have an orbital with the following asymptotic behaviour (5) $$\psi \left(\mathbf{r}\right)\sim e^{-r}\left(cos{\left(\theta \right)}^2+\frac{1}{r^2}\right),$$(5) where θ = π/2 defines the xy plane in spherical coordinates. This orbital has the same asymptotic exponential decay everywhere, but on the xy plane it decays 1/r² faster than elsewhere. We can compute the corresponding potential by inversion, $v\left(\mathbf{r}\right)={\nabla }^2\psi \left(\mathbf{r}\right)/2\psi \left(\mathbf{r}\right)$, and we find that v(r → ∞, θ) = 1/2, but v(r → ∞, π/2) = 3/2. The potential is shown in Figure 1 as a function of r and θ (we have subtracted 1/2 so that $v\left(\mathbf{r}\right)$ goes asymptotically to zero outside the plane). We clearly see that this potential has a ‘ridge’ on the xy plane, where it goes asymptotically to a constant. The ‘ridge’ shrinks as r gets larger and larger. This kind of behaviour in the potential has been usually associated in the literature to a change in the exponential decay of the orbital on the plane, but we see here a clear counterexample.

Figure 1. The large-r behaviour of the potential that generates the orbital with the asymptotic behaviour of Equation (5(5) $$\psi \left(\mathbf{r}\right)\sim e^{-r}\left(cos{\left(\theta \right)}^2+\frac{1}{r^2}\right),$$(5) ) as a function of r and the azimuthal angle θ = arccos(z/r). We clearly see that this potential has a ‘ridge’ on the xy plane (corresponding to θ = π/2), where it goes asymptotically to a constant. The ‘ridge’ shrinks as r gets larger and larger. Despite this constant, the orbital has the same exponential decay everywhere, it only decays 1/r² faster on the plane.

As a second example, we consider an orbital with an exponentially faster decay on the plane: (6) $$\psi \left(\mathbf{r}\right)\sim e^{-r}cos{\left(\theta \right)}^2+e^{-2r}.$$(6) Again, we compute the corresponding potential by inversion, and we find that this potential has also a ‘ridge’ on the plane, where it diverges exponentially, as shown in Figure 2. Again, the ‘ridge’ shrinks as r → ∞. This example shows that a different exponential decay of the orbital on the plane does not necessarily imply that the corresponding potential goes to a constant on the plane: as we see here, the potential might diverge, in order to compensate the derivative perpendicular to the plane, which increases exponentially as r increases.

Figure 2. The large-r behaviour of the potential that generates the orbital with the asymptotic behaviour of Equation (6(6) $$\psi \left(\mathbf{r}\right)\sim e^{-r}cos{\left(\theta \right)}^2+e^{-2r}.$$(6) ) as a function of r and the azimuthal angle θ = arccos(z/r). We clearly see that this potential has a ‘ridge’ on the xy plane (corresponding to θ = π/2), where it diverges exponentially. The ‘ridge’ shrinks as r gets larger and larger.

In Section 5, we will discuss the case that there is slower decay of an orbital on the plane than elsewhere, where it will be demonstrated that dependent on how this slower decay on the plane is approached (θ-dependence), the potential may exhibit asymptotically either a negative constant, or may diverge to +∞ or to −∞. Even in the case of a negative constant on the plane, we will see that the potential close to the plane displays very irregular behaviour.

3. Asymptotic behaviour of the exact density from the Dyson orbital expansion

Although it is commonly assumed in the literature on the asymptotic decay of the density that the exact density decays everywhere in the same way, it should be recognised that this is not always true. First of all, it is evident for a non-interacting electron system that if there is a nodal plane in the HOMO, the density decays differently in that plane, because the HOMO density does not contribute in the plane. However, the interacting case, where a configuration mixing will involve many configurations, is more subtle. While in a non-interacting system described by a single determinant, a HOMO nodal plane determines a nodal plane for the whole wavefunction (the wavefunction when all the electrons are on the nodal plane is zero), such a nodal plane in general does not survive in the many-configuration interacting wavefunction. It is also known that in the case of nodes in the wavefunction due to the fermionic character of the electrons (antisymmetry of the wavefunction under permutation), electron interaction substantially modifies the nodes [11,12].

To study the density decay in the general interacting case, we express the exact N-electron wavefunction Ψ^N₀ and the exact density in terms of the Dyson orbitals $d_i\left(\mathbf{x}\right)$, (7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) where the Ψ^{N − 1}_i are the exact (N − 1)-electron states and $\mathbf{x}=\mathbf{r},s$. The sum over i goes over both the spin-↑ and spin-↓ Dyson orbitals. If e.g. $\mathbf{x}=(\mathbf{r},s=\uparrow )$ then only the spin-↑ Dyson orbitals are non-zero at $\mathbf{x}$ and contribute to $n\left(\mathbf{x}\right)=n(\mathbf{r},\uparrow )$ ($=\frac{1}{2}n\left(\mathbf{r}\right)$ in closed shell systems). Each state of the ion is associated with a one-particle wavefunction, its Dyson orbital. These orbitals constitute a non-orthogonal non-normal, in general, linearly dependent set. We define the conditional amplitude $\Phi (2\cdots N;\mathbf{x})$ [13] and associated quantities, (8) $$\begin{array}{ccc}\hfill \Phi (2\cdots N;\mathbf{x})& =& \frac{{\Psi }_0^N(\mathbf{x},2\cdots N)}{\sqrt{n\left(\mathbf{x}\right)/N}},\hfill \\ \hfill n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})& =& {(N-1)\int \left|\Phi (2\cdots N|\mathbf{x})\right|}^2\mathrm{d}3\cdots \mathrm{d}N,\hfill \\ \hfill v^{\text{cond}}\left(\mathbf{x}\right)& =& \int \frac{n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})}{|\mathbf{r}-{\mathbf{r}}_2|}\mathrm{d}{\mathbf{x}}_2.\hfill \end{array}$$(8) $\Phi (2\cdots N;\mathbf{x})$ is a normalised (N − 1)-electron wavefunction depending parametrically on the position $\mathbf{x}$. Its square describes the probability distribution of electrons at positions 2⋅⋅⋅N when one electron is known to be at $\mathbf{x}$. Its associated one-electron density $n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})$ is the density of the other electrons at position ${\mathbf{x}}_2$ when one electron is at $\mathbf{x}$, which is the normal one-electron density $n\left({\mathbf{x}}_2\right)$ plus the full exchange-correlation hole surrounding position $\mathbf{x}$, $n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})=n\left({\mathbf{x}}_2\right)+n_{xc}^{\text{hole}}\left({\mathbf{x}}_2\right|\mathbf{x})$, $v^{\text{cond}}\left(\mathbf{x}\right)=\int d{\mathbf{x}}_2{n}^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})/|{\mathbf{x}}_2-\mathbf{x}|=v_{H\text{artree}}\left(\mathbf{x}\right)+v_{xc}^{\text{hole}}\left(\mathbf{x}\right)$. Projecting the Schrödinger equation ${\widehat{H}}^N{\Psi }_0^N=E_0^N{\Psi }_0^N$ against Ψ^{N − 1}_i(2⋅⋅⋅N) and using the expansion of Equation (7(7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) ), one obtains the energy-independent equations for the Dyson orbitals, (9) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)\right)d_i\left(\mathbf{x}\right)+\sum _{k=0}^{\infty }{X}_{ik}\left(\mathbf{x}\right)d_k\left(\mathbf{x}\right)=-I_i{d}_i\left(\mathbf{x}\right).\end{array}$$(9) Katriel and Davidson (KD) [14] pointed out that, due to the coupling integrals (10) $$X_{ik}\left(\mathbf{x}\right)\equiv {⟨{\Psi }_i^{N-1}\left|\sum _{j>1}^N\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}\right|{\Psi }_k^{N-1}⟩}_{2..N},$$(10) the exponential decay of the coupled Dyson orbitals will be the same. This was demonstrated by Handy et al. [15] for the analogous case of coupling of the Hartree–Fock orbitals by the exchange term. The first Dyson orbital d₀ will have exponential decay $\sim e^{-\sqrt{2I_0}r}$ multiplied by a factor r^β with $\beta =(Z-N+1)/\sqrt{2I_0}-1$, due to the −Z/r decay of v_ext and the (N − 1)/r decay of the coupling term. KD find that higher Dyson orbitals which have non-zero X_i0 with the first Dyson orbital will have decay $r^{\beta -L^{*}}{e}^{-\sqrt{2I_0}r}$ with L* ≥ 2. Dyson orbitals that are not connected to d₀ will have different exponential decay, governed by the eigenvalue of the first orbital in such a connected set (which is disjunct from other sets). Considering the expansion of the density in Dyson orbitals in Equation (7(7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) ), KD have concluded that, if the density decays for $\left|\mathbf{r}\right|\to \infty$ as the most slowly decaying term $|d_0{\left(\mathbf{x}\right)|}^2$, its exponential decay would be $\sim e^{-2\sqrt{2I_0}r}$. Levy, Perdew and Sahni (LPS) [3] proved this exponential decay in a different way, thereby proving that the leading term is not overruled by the infinite sum of the faster decaying terms in Equation (7(7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) ). The result (11) $$n\left(\right|\mathbf{r}|\to \infty )\sim |d_0{\left(\mathbf{r}\right)|}^2\sim e^{-2\sqrt{2I_0}r}$$(11) is considered well established.

This picture changes if the KS HOMO has a nodal plane. The common thinking is that always the interacting density decays everywhere in the same way due to correlation effects [6,16]. Instead, by analysing the Dyson orbitals, we can see that a KS HNP may also imply a special behaviour of the interacting density. A nodal plane extending to infinity is typically related to a symmetry plane (consider, e.g. ethylene or benzene [7,9]), but a HOMO nodal surface can occur also in more general situations. In the case of a symmetry plane, the exact interacting states of the molecule are either symmetric or antisymmetric with respect to the plane. For example, the ground-state wavefunction corresponding to a closed-shell configuration is totally symmetric with respect to that plane, while the first ion state Ψ^{N − 1}₀ will be antisymmetric (the KS first ion state surely will be so, and we will consider the usual case that the same holds for the exact ion state). For points ${\mathbf{r}}_p$ in the HNP, the conditional amplitude $\Phi (2\cdots N|{\mathbf{x}}_p)$ will be symmetric with respect to the plane. Therefore, the matrix element ${⟨{\Psi }_0^{N-1}|\Phi (2\cdots N|{\mathbf{x}}_p)⟩}_{2..N}$ will vanish, so that the first Dyson orbital is zero in the plane: (12) $$\begin{array}{cc}& d_0\left({\mathbf{x}}_p\right)=\sqrt{n\left({\mathbf{x}}_p\right)}{⟨{\Psi }_0^{N-1}|\Phi (2\cdots N|{\mathbf{x}}_p)⟩}_{2..N}=0.\hfill \end{array}$$(12) In fact, d₀ is antisymmetric with respect to the plane. To obtain the asymptotic behaviour of higher Dyson orbitals in the plane, we have to solve Equation (9(9) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)\right)d_i\left(\mathbf{x}\right)+\sum _{k=0}^{\infty }{X}_{ik}\left(\mathbf{x}\right)d_k\left(\mathbf{x}\right)=-I_i{d}_i\left(\mathbf{x}\right).\end{array}$$(9) ) for points ${\mathbf{x}}_p$ in the plane. Since $d_0\left({\mathbf{x}}_p\right)=0$, it looks as if the coupling to d₀ will be absent for any higher Dyson orbital d_{i > 0}. The decay in the HNP of the second Dyson orbital (and, thus, of the density) is then not governed by d₀, but will be according to the second ionisation potential, I₁ (if d₁ is not also zero on the plane). We will in Section 4 discuss this situation (Case 1), which would allow for a smooth KS potential that decays uniformly as −1/r. However, a closer scrutiny of Equation (9(9) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)\right)d_i\left(\mathbf{x}\right)+\sum _{k=0}^{\infty }{X}_{ik}\left(\mathbf{x}\right)d_k\left(\mathbf{x}\right)=-I_i{d}_i\left(\mathbf{x}\right).\end{array}$$(9) ) reveals that in general d₁ and the density will inherit in the HNP the slower decay according to I₀ from d₀. In the KS case, the HOMO−1 does not couple to another orbital, and the asymptotic behaviour of the density in the HNP is determined by the HOMO−1. Given its eigenvalue close to −I₁, and the corresponding ‘fast’ decay everywhere else, the slow decay according to I₀ on only the HNP requires very irregular features in the KS potential, as will be discussed in Section 5.

4. Case 1: The density decay on the HNP is governed by the second ionisation potential

4.1. Asymptotic behaviour of the density

Since d₀ = 0 in the HNP, we have to turn to $d_1\left({\mathbf{x}}_p\right)$ to determine the asymptotic behaviour of the density on the plane, $n\left(\right|{\mathbf{r}}_p|\to \infty )$. If the corresponding excited ion state Ψ^{N − 1}₁ has the same symmetry with respect to the HNP as Ψ^N₀, this Dyson orbital will not be zero in the plane (see Equation (7(7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) )). In this section, we explore the case that there would be no coupling between d₁ and d₀, which would be the case if X₁₀ ≡ 0. For simplicity, we only consider the coupling of d₁ to d₀ and consider the asymptotic terms in the equation for d₁ (13) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{ext}\right)d_1+X_{11}{d}_1+X_{10}{d}_0=-I_1{d}_1\end{array}$$(13) The term with X₁₁ is not problematic since X₁₁ ∼ (N − 1)/r, so it can be combined with v_ext with the same asymptotic behaviour to give a Q/r term with Q = −Z + N − 1. On the other hand, X₁₀ can play a large role (14) $$\begin{array}{cc}\hfill X_{10}& \left(\mathbf{x}\right)=\int {\Psi }_1^{N-1}{({\mathbf{x}}_2\cdots {\mathbf{x}}_N)}^{*}\hfill \\ & \sum _{j>1}\frac{1}{|r_j-r|}{\Psi }_0^{N-1}({\mathbf{x}}_2\cdots {\mathbf{x}}_N)d{\mathbf{x}}_2\cdots d{\mathbf{x}}_n\hfill \end{array}$$(14) We are dealing with the situation that the ion ground state Ψ^{N − 1}₀ is antisymmetric under reflection of all electronic coordinates with respect to the HNP, while Ψ^{N − 1}₁ is symmetric. The operator in Equation (14(14) $$\begin{array}{cc}\hfill X_{10}& \left(\mathbf{x}\right)=\int {\Psi }_1^{N-1}{({\mathbf{x}}_2\cdots {\mathbf{x}}_N)}^{*}\hfill \\ & \sum _{j>1}\frac{1}{|r_j-r|}{\Psi }_0^{N-1}({\mathbf{x}}_2\cdots {\mathbf{x}}_N)d{\mathbf{x}}_2\cdots d{\mathbf{x}}_n\hfill \end{array}$$(14) ) can be split in a symmetric and an antisymmetric part with respect to reflection of the coordinates ${\mathbf{r}}_j$ in the (xy) plane, (15) $$\begin{array}{cc}\hfill \sum _{j>1}\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}& =\frac{1}{2}\sum _{j>1}\left(\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}-\frac{1}{|{\mathbf{r}}_j-\widehat{{\sigma }_h}\mathbf{r}|}\right)\hfill \\ & +\frac{1}{2}\sum _{j>1}\left(\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}+\frac{1}{|{\mathbf{r}}_j-\widehat{{\sigma }_h}\mathbf{r}|}\right)\hfill \end{array}$$(15) where $\widehat{{\sigma }_h}$ is the operator for reflection in the horizontal (xy) plane, $\widehat{{\sigma }_h}(x\widehat{i}+y\widehat{j}+z\widehat{k})=x\widehat{i}+y\widehat{j}-z\widehat{k}$, where $\widehat{i},\widehat{j},\widehat{k}$ are unit vectors along the x, y and z axes. Only the antisymmetric part will survive, and one may consider the asymptotic behaviour for points $\mathbf{r}$ far beyond the positions ${\mathbf{r}}_j$ occurring in the integral, which range over the limited domain where the ion wavefunctions Ψ^{N − 1}₀ and Ψ^{N − 1}₁ have appreciable values. One may deduce that the leading term in the asymptotic behaviour of X₁₀ will be (16) $$\begin{array}{cc}& X_{10}(\mathbf{r}\to \infty )=\frac{cos\theta }{r^2}\int {\Psi }_1^{N-1}{({\mathbf{x}}_2\cdots {\mathbf{x}}_N)}^{*}\hfill \\ & \left(\sum _{j>1}{z}_j\right){\Psi }_0^{N-1}({\mathbf{x}}_2\cdots {\mathbf{x}}_N)d{\mathbf{x}}_2\cdots d{\mathbf{x}}_n=\frac{kcos\theta }{r^2}\hfill \end{array}$$(16) This leads to a term (kcos θ/r²)d₀ in Equation (13(13) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{ext}\right)d_1+X_{11}{d}_1+X_{10}{d}_0=-I_1{d}_1\end{array}$$(13) ), where k is the integral in (16(16) $$\begin{array}{cc}& X_{10}(\mathbf{r}\to \infty )=\frac{cos\theta }{r^2}\int {\Psi }_1^{N-1}{({\mathbf{x}}_2\cdots {\mathbf{x}}_N)}^{*}\hfill \\ & \left(\sum _{j>1}{z}_j\right){\Psi }_0^{N-1}({\mathbf{x}}_2\cdots {\mathbf{x}}_N)d{\mathbf{x}}_2\cdots d{\mathbf{x}}_n=\frac{kcos\theta }{r^2}\hfill \end{array}$$(16) ). The symmetries with respect to the nodal plane (both Ψ^{N − 1}₀ and the operator in the matrix element k are antisymmetric with respect to the HNP, while Ψ^{N − 1}₁ is symmetric) do not force the k integral to be zero. However, in special cases k will still be zero for symmetry reasons, if for instance Ψ^{N − 1}₁, Ψ^{N − 1}₀ and the operator ∑_{j > 1}z_j belong to such irreducible representations of the molecular point group that the integral is zero. (The same argument then holds for further terms in the expansion of (15(15) $$\begin{array}{cc}\hfill \sum _{j>1}\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}& =\frac{1}{2}\sum _{j>1}\left(\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}-\frac{1}{|{\mathbf{r}}_j-\widehat{{\sigma }_h}\mathbf{r}|}\right)\hfill \\ & +\frac{1}{2}\sum _{j>1}\left(\frac{1}{|{\mathbf{r}}_j-\mathbf{r}|}+\frac{1}{|{\mathbf{r}}_j-\widehat{{\sigma }_h}\mathbf{r}|}\right)\hfill \end{array}$$(15) ) with odd powers of z_j.) This is not an esoteric possibility. It holds, for instance, for the prototype molecules ethylene and benzene. Suppose that d₁ also does not couple to d₀ indirectly (through coupling to higher Dyson orbitals which themselves might couple to d₀), then d₁ will decay in the plane as $e^{-\sqrt{2I_1}{r}_p}$ like everywhere else, since in the eigenvalue Equation (9(9) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)\right)d_i\left(\mathbf{x}\right)+\sum _{k=0}^{\infty }{X}_{ik}\left(\mathbf{x}\right)d_k\left(\mathbf{x}\right)=-I_i{d}_i\left(\mathbf{x}\right).\end{array}$$(9) ) all the terms except those from ∇² can be neglected in the asymptotic region compared to I₁. According to Equation (7(7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) ), the density will, for points in the HNP, decay as $|d_1\left({\mathbf{x}}_p\right){|}^2\sim e^{-2\sqrt{2I_1}{r}_p}$, i.e. different in the plane than outside the plane. We have to consider the possibility that this special case occurs. As will be seen, it is the only case where the KS potential can have the simple, generally assumed uniform asymptotic −1/r behaviour. This will be discussed as Case 1 in this section, while in Section 5, we will discuss Case 2 where the density has the slow decay according to I₀ also on the plane, either because the k integral is not zero and d₁ couples to d₀, or because higher Dyson orbitals couple to d₀. We emphasise that it does not appear to be likely that no higher d_i(i > 1) would couple to d₀, so Case 1 with its regular −1/r asymptotic KS potential must be exceptional, if it exists at all.

Supposing then that neither d₁ nor the higher Dyson orbitals d_i(i > 1) have the slow $e^{-\sqrt{2I_0}{r}_p}$ decay, we proceed to show that a consistent picture of the density decay in the HNP can be given, with the regular −1/r asymptotic behaviour of the KS potential. First, we consider the question if the decay $e^{-2\sqrt{2I_1}{r}_p}$ of $|d_1\left({\mathbf{r}}_p\right){|}^2$ will actually be the decay of the total density, i.e. whether the infinite summation over the other Dyson orbitals squared does not overrule the decay of the first term, cf. Equation (7(7) $$\begin{array}{ccc}\hfill {\Psi }_0^N& =& N^{-1/2}\sum _{i=0}^{\infty }{d}_i\left(\mathbf{x}\right){\Psi }_i^{N-1}(2\cdots N),\hfill \\ \hfill d_i\left(\mathbf{x}\right)& =& \sqrt{N}\int {\Psi }_i^{N-1}{(2\cdots N)}^{*}{\Psi }_0^N(\mathbf{x},2\cdots N)\mathrm{d}2\cdots \mathrm{d}N,\hfill \\ \hfill n\left(\mathbf{x}\right)& =& \sum _{i=0}^{\infty }{\left|d_i\left(\mathbf{x}\right)\right|}^2,\hfill \end{array}$$(7) ). The relation between d₀ and the density as in Equation (12(12) $$\begin{array}{cc}& d_0\left({\mathbf{x}}_p\right)=\sqrt{n\left({\mathbf{x}}_p\right)}{⟨{\Psi }_0^{N-1}|\Phi (2\cdots N|{\mathbf{x}}_p)⟩}_{2..N}=0.\hfill \end{array}$$(12) ) (but in general directions) has been used by LPS [3] to show that indeed the decay of $|d_0\left({\mathbf{r}}_p\right){|}^2$ and the density are the same. In the HNP, we now have to use the analogous equation for d₁, (17) $$\begin{array}{cc}& d_1\left({\mathbf{x}}_p\right)=\sqrt{n\left({\mathbf{x}}_p\right)}⟨{\Psi }_1^{N-1}|\Phi (2\cdots N|{\mathbf{x}}_p){⟩}_{2..N}.\hfill \end{array}$$(17) We show that the special properties of Φ in this case afford the required relation between $|d_1\left({\mathbf{x}}_p\right){|}^2$ and $n\left({\mathbf{x}}_p\right)$ for $|{\mathbf{r}}_p|\to \infty$. The conditional amplitude Φ is a normalised (N − 1)-electron wavefunction that describes, when $|{\mathbf{r}}_p|\to \infty$, the probability distribution of the electrons that remain behind when one electron is infinitely far away (in this case, in the plane). The integral in Equation (17(17) $$\begin{array}{cc}& d_1\left({\mathbf{x}}_p\right)=\sqrt{n\left({\mathbf{x}}_p\right)}⟨{\Psi }_1^{N-1}|\Phi (2\cdots N|{\mathbf{x}}_p){⟩}_{2..N}.\hfill \end{array}$$(17) ) must be ≤1, since Φ and Ψ^{N − 1}₁ are both normalised. We can also show that it cannot decay to zero, but will go to a finite constant. Consider the expansion of the ground-state wavefunction in the leading KS-independent particle determinantal wavefunction Ψ^N_{s, 0}, plus all its excitations (which are orthogonal to the leading term). It is elementary to show from properties of the determinant that for ${\mathbf{r}}_p$ in the HNP, the conditional amplitude ${\Phi }_s(2\cdots N|{\mathbf{x}}_p)$ of the KS determinantal wavefunction reduces to the second ion state Ψ^{N − 1}_{s, 1} of the non-interacting ion (the determinant with a hole in HOMO−1) for r_p → ∞ (see Appendix). Therefore, the full Φ will consist in large part of this KS ion state (the contribution of the HF or KS determinant in the wavefunction is typically substantial, 80%–90% is not uncommon). The overlap of the exact second ion state with the second KS ion state, ⟨Ψ^{N − 1}₁|Ψ_{s, 1}^{N − 1}⟩, will be a finite constant. The integral in Equation (17(17) $$\begin{array}{cc}& d_1\left({\mathbf{x}}_p\right)=\sqrt{n\left({\mathbf{x}}_p\right)}⟨{\Psi }_1^{N-1}|\Phi (2\cdots N|{\mathbf{x}}_p){⟩}_{2..N}.\hfill \end{array}$$(17) ), therefore, remains finite. We will show below that the integral is actually 1, since Φ ‘collapses’ to Ψ^{N − 1}₁ for r_p → ∞, but at this point the fact that the integral goes to some finite constant is sufficient to see that the exponential decay of $\sqrt{n\left({\mathbf{x}}_p\right)}$ must be the same as that of $d_1\left({\mathbf{x}}_p\right)$. We can thus conclude that in the present Case 1, the exact density will have a different (faster) decay in the HNP ($e^{-2\sqrt{2I_1}{r}_p}$) than in general directions ($e^{-2\sqrt{2I_0}{r}_p}$).

4.2. Asymptotic behaviour of the Kohn–Sham potential

What happens in this case with the KS potential for $|{\mathbf{r}}_p|\to \infty$?

The argument of Ref. 6 for a negative constant does not apply here since it was based on uniform decay of the density in all directions like $\text{exp}[-2\sqrt{2I_0}r]$ while in Case 1 the exact density has a different decay in the HNP than elsewhere. If the exact density decays like $e^{-2\sqrt{2I_1}{r}_p}$ on the HNP, the decay $\text{exp}[-2\sqrt{2(-{ϵ}_{H-1}+v_s\left(\infty \right))}{r}_p]$ of the HOMO−1 density can represent the density decay with a simple KS potential that goes uniformly like −1/r (v_s(∞) = 0) if ε_{H − 1} = −I₁. Since the decay of both HOMO−1 and HOMO are uniform (the same outside and in the plane), the problems related to the angular derivatives (Section 2) do not appear and the uniform −1/r asymptotic behaviour of the KS potential is consistent with the asymptotic behaviour of these solutions.

It is interesting to observe that the few accurate (but not exact) calculations that are available for the KS orbital energies [17,18] show that ε_{H − 1} ≈ −I₁. The equality has not been established to better than ca. 0.05 eV, since the KS orbital energies are only obtained to ca. 0.05 eV accuracy (the calculations of the ‘exact’ orbital energies use KS potentials that are generated with the criterium that they reproduce accurate CI densities).

A positive asymptotic constant for the KS potential in the HNP has been found in Refs. 7–10 for the exact-exchange model. This does not imply that a positive constant on the plane is also present in the asymptotic behaviour of the full KS potential [16,19]. Moreover, it is not clear whether this positive constant would modify the (exponential) decay of HOMO−1 in the nodal plane, see Section 2, but if it did and the exponential decay of the HOMO−1 density would be $\text{exp}[-2\sqrt{2(-{ϵ}_{H-1}+C)}{r}_p]\approx \text{exp}[-2\sqrt{2(I_1+C)}{r}_p]$, different from the exact density, this would not have any consequence. If the local potential of the non-interacting electron system is not explicitly required to reproduce the exact density, but is determined by some other criterium such as minimum EXX energy, we may encounter essential differences between the non-interacting density and the exact one.

In the remainder of this section, we investigate if this simple Case 1 situation (different decay of the exact density on HNP than elsewhere and a regular −1/r asymptotics of the KS potential in all directions) is consistent with what is known analytically about v_s. We consider the KS potential in the convenient form in which it can be written [20]: (18) $$\begin{array}{c}\hfill v_s=v_{\mathrm{ext}}+v^{\text{cond}}+\left(v^{\text{kin}}-v_s^{\text{kin}}\right)+(v^{N-1}-v_s^{N-1}).\end{array}$$(18) v^cond has been defined in and below Equation (8(8) $$\begin{array}{ccc}\hfill \Phi (2\cdots N;\mathbf{x})& =& \frac{{\Psi }_0^N(\mathbf{x},2\cdots N)}{\sqrt{n\left(\mathbf{x}\right)/N}},\hfill \\ \hfill n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})& =& {(N-1)\int \left|\Phi (2\cdots N|\mathbf{x})\right|}^2\mathrm{d}3\cdots \mathrm{d}N,\hfill \\ \hfill v^{\text{cond}}\left(\mathbf{x}\right)& =& \int \frac{n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})}{|\mathbf{r}-{\mathbf{r}}_2|}\mathrm{d}{\mathbf{x}}_2.\hfill \end{array}$$(8) ) and the definitions of v_kin and v^{N − 1} appear in the second and third lines, respectively, of the expression for the effective potential $v_{\mathrm{eff}}\left(\mathbf{r}\right)$ for $\sqrt{n\left(\mathbf{r}\right)}$, Equation (1(1) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{eff}}\left(\mathbf{r}\right)\right)\sqrt{n\left(\mathbf{r}\right)}=-I_0\sqrt{n\left(\mathbf{r}\right)}\end{array}$$(1) ), see [3,20] (19) $$\begin{array}{ccc}\hfill v_{\mathrm{eff}}\left(\mathbf{x}\right)& =& v^{\text{cond}}\left(\mathbf{x}\right)+\frac{1}{2}⟨{\nabla }_x\Phi (2\cdots N|\mathbf{x})|{\nabla }_x\Phi (2\cdots N|\mathbf{x})⟩\hfill \\ & & +⟨\Phi (2\cdots N|\mathbf{x})\left|{\widehat{H}}^{N-1}-E_0^{N-1}\right|\Phi (2\cdots N|\mathbf{x})⟩\hfill \\ & \equiv & v^{\text{cond}}\left(\mathbf{x}\right)+v^{\text{kin}}\left(\mathbf{x}\right)+v^{N-1}\left(\mathbf{x}\right).\hfill \end{array}$$(19) The potentials $v_s^{\text{kin}}$ and v^{N − 1}_s in v_s are defined by replacing the exact conditional amplitude Φ in the definitions of v^kin and v^{N − 1} (see (19(19) $$\begin{array}{ccc}\hfill v_{\mathrm{eff}}\left(\mathbf{x}\right)& =& v^{\text{cond}}\left(\mathbf{x}\right)+\frac{1}{2}⟨{\nabla }_x\Phi (2\cdots N|\mathbf{x})|{\nabla }_x\Phi (2\cdots N|\mathbf{x})⟩\hfill \\ & & +⟨\Phi (2\cdots N|\mathbf{x})\left|{\widehat{H}}^{N-1}-E_0^{N-1}\right|\Phi (2\cdots N|\mathbf{x})⟩\hfill \\ & \equiv & v^{\text{cond}}\left(\mathbf{x}\right)+v^{\text{kin}}\left(\mathbf{x}\right)+v^{N-1}\left(\mathbf{x}\right).\hfill \end{array}$$(19) )) with the conditional amplitude for the KS determinantal wavefunction, Φ_s.

LPS [3] noted that each term in Equation (19(19) $$\begin{array}{ccc}\hfill v_{\mathrm{eff}}\left(\mathbf{x}\right)& =& v^{\text{cond}}\left(\mathbf{x}\right)+\frac{1}{2}⟨{\nabla }_x\Phi (2\cdots N|\mathbf{x})|{\nabla }_x\Phi (2\cdots N|\mathbf{x})⟩\hfill \\ & & +⟨\Phi (2\cdots N|\mathbf{x})\left|{\widehat{H}}^{N-1}-E_0^{N-1}\right|\Phi (2\cdots N|\mathbf{x})⟩\hfill \\ & \equiv & v^{\text{cond}}\left(\mathbf{x}\right)+v^{\text{kin}}\left(\mathbf{x}\right)+v^{N-1}\left(\mathbf{x}\right).\hfill \end{array}$$(19) ) is everywhere non-negative and should tend to zero asymptotically. In fact, $v^{\text{cond}}\left(\mathbf{x}\right)$ (see Equation (8(8) $$\begin{array}{ccc}\hfill \Phi (2\cdots N;\mathbf{x})& =& \frac{{\Psi }_0^N(\mathbf{x},2\cdots N)}{\sqrt{n\left(\mathbf{x}\right)/N}},\hfill \\ \hfill n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})& =& {(N-1)\int \left|\Phi (2\cdots N|\mathbf{x})\right|}^2\mathrm{d}3\cdots \mathrm{d}N,\hfill \\ \hfill v^{\text{cond}}\left(\mathbf{x}\right)& =& \int \frac{n^{\text{cond}}\left({\mathbf{x}}_2\right|\mathbf{x})}{|\mathbf{r}-{\mathbf{r}}_2|}\mathrm{d}{\mathbf{x}}_2.\hfill \end{array}$$(8) )), being the repulsive Coulomb potential of a localised charge distribution of (N − 1) electrons, decays like (N − 1)/r. With v_ext = −Z/r and Z = N for neutral systems, the expected −1/r behaviour emerges in v_s and v_eff if the remaining terms are asymptotically zero. However, in the presence of a HNP, the asymptotic behaviour of the other terms is more complicated. Clearly, we will find −1/r behaviour for v_s if both $v^{\text{kin}}-v_s^{\text{kin}}$, and v^{N − 1} − v^{N − 1}_s are asymptotically zero.

Considering first v^{N − 1} − v^{N − 1}_s, also called the response potential v^resp [17,20,21], we note that v^{N − 1} is positive since in general Φ will not be the ground-state wavefunction of the ion, so its expectation value will be larger than E^{N − 1}₀. When $\left|\mathbf{r}\right|\to \infty ,$ it has been inferred that, in general, the conditional amplitude collapses to the ion ground state Ψ^{N − 1}₀ [14] (when s = ↑, then Φ will collapse to the M_S = −1/2 state of the doublet ion), so that $v^{N-1}\left(\right|\mathbf{r}|\to \infty )\to 0$. But in the HNP, this changes. By expanding the conditional amplitude $\Phi (2\cdots N|\mathbf{x})$ in terms of the exact (N − 1)-electron states, (20) $$\Phi (2\cdots N|\mathbf{x})=\sum _{i=0}^{\infty }\frac{d_i\left(\mathbf{x}\right)}{\sqrt{n\left(\mathbf{x}\right)}}{\Psi }_i^{N-1}(2\cdots N),$$(20) we see that, since on the HNP d₀ = 0 and $|d_1\left({\mathbf{x}}_p\right)|(r_p\to \infty )\sim \sqrt{n\left({\mathbf{x}}_p\right)}$, while all higher d_i decay a factor $r^{-L^{*}}$ faster, with L* ≥ 2, the conditional amplitude tends asymptotically for ${\mathbf{x}}_p\to \infty$ on the plane to the first-excited ion state, Φ → Ψ^{N − 1}₁ (note that Φ is normalised for any position $\mathbf{x}$). This implies that (21) $$v^{N-1}\left(\right|{\mathbf{r}}_p|\to \infty )=E_1^{N-1}-E_0^{N-1}=I_1-I_0.$$(21) This is a positive constant. It would appear in the asymptotics of v_eff only on the HNP. It can be shown that v^{N − 1}_s also goes to the constant I₁ − I₀ for asymptotic points ${\mathbf{x}}_p$ in the nodal plane, and therefore cancels v^{N − 1}, see Equation (21(21) $$v^{N-1}\left(\right|{\mathbf{r}}_p|\to \infty )=E_1^{N-1}-E_0^{N-1}=I_1-I_0.$$(21) ). We have already noticed that the conditional amplitude of the non-interacting KS system with determinantal ground state collapses to the second ion state of the non-interacting system for points in the HOMO nodal plane, so (22) $$\begin{array}{cc}& v_s^{N-1}\left({\mathbf{x}}_p\right)\hfill \\ & =⟨{\Phi }_s(2\cdots N|{\mathbf{x}}_p)\left|H_s^{N-1}\right|{\Phi }_s(2\cdots N|{\mathbf{x}}_p)⟩-E_{s,0}^{N-1}\hfill \\ & =⟨{\Psi }_{s,1}^{N-1}(2\cdots N)\left|H_s^{N-1}\right|{\Psi }_{s,1}^{N-1}(2\cdots N)⟩-E_{s,0}^{N-1}\hfill \\ & =E_{s,1}^{N-1}-E_{s,0}^{N-1}={ϵ}_H-{ϵ}_{H-1}=I_1-I_0\hfill \end{array}$$(22) The response potential v^{N − 1} − v^{N − 1}_s, therefore, goes to zero.

Turning next to $v_c^{\text{kin}}=v^{\text{kin}}-v_s^{\text{kin}}$, we observe that v^kin can also be non-zero at infinity: when crossing the HNP, the asymptotic conditional amplitude changes from Ψ^{N − 1}₀ to Ψ^{N − 1}₁, so that the $\mathbf{r}$-derivative of Φ perpendicular to the plane can be non-zero on the HNP also when $\left|\mathbf{r}\right|\to \infty$. The behaviour of v^kin then depends on how $d_0(\mathbf{r}\to {\mathbf{r}}_p)$ goes to zero when approaching the nodal plane. We note that for the determinantal wavefunction of a non-interacting system, the Dyson orbitals are precisely the occupied independent particle orbitals. In the interacting system, the first Dyson orbitals for primary ion states (those corresponding to a simple orbital ionisation) still are very similar to the KS orbitals: overlaps are typically >0.999 [22]. This agrees with our finding in this paper that when the KS HOMO is antisymmetric with respect to a plane, the corresponding Dyson orbital also is antisymmetric with respect to that plane. Let us then take as example that asymptotically, in spherical coordinates, $d_0\sim f(cos\theta )R\left(r\right)e^{-\sqrt{2I_0}r}$, with f(0) = 0, and f′(0) ≠ 0, as would be the case for a π orbital, which has fR = rcos θ = z. By writing v^kin in the form [17,20,21] (23) $$v^{\text{kin}}\left(\mathbf{r}\right)=\sum _{i=1}^{\infty }\frac{|\nabla d_i{\left(\mathbf{x}\right)|}^2}{n\left(\mathbf{r}\right)}-\frac{{\left|\nabla n\left(\mathbf{r}\right)\right|}^2}{8n{\left(\mathbf{r}\right)}^2},$$(23) and using $d_1\sim fR{e}^{-\sqrt{2I_1}r}$, it is easy to see that (24) $$v^{\text{kin}}(r_p\to \infty )\to \frac{1}{2}{f}^{\text{'}}{\left(0\right)}^2\frac{R^2}{r^2}{e}^{2(\sqrt{2I_1}-\sqrt{2I_0})r},$$(24) showing that v^kin can go asymptotically to infinity on the HNP. (This is not detrimental for the solution of $\sqrt{n}$ with Equation (1(1) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)+v_{\mathrm{eff}}\left(\mathbf{r}\right)\right)\sqrt{n\left(\mathbf{r}\right)}=-I_0\sqrt{n\left(\mathbf{r}\right)}\end{array}$$(1) ) since it can be shown that terms coming from ∇² cancel this divergence of v^kin.) The complete kinetic term $v_c^{\text{kin}}=v^{\text{kin}}-v_s^{\text{kin}}$ can be evaluated using the expression for $v_s^{\text{kin}}$ analogous to Equation (23(23) $$v^{\text{kin}}\left(\mathbf{r}\right)=\sum _{i=1}^{\infty }\frac{|\nabla d_i{\left(\mathbf{x}\right)|}^2}{n\left(\mathbf{r}\right)}-\frac{{\left|\nabla n\left(\mathbf{r}\right)\right|}^2}{8n{\left(\mathbf{r}\right)}^2},$$(23) ), but now written for the KS wavefunction (note that the H KS orbitals ψ_i, i = 1..H, are the exact Dyson orbitals of the non-interacting KS system, which are a finite number in this case), (25) $$v_s^{\text{kin}}\left(\mathbf{r}\right)=\sum _{i=1}^H\frac{|\nabla {\psi }_i{\left(\mathbf{x}\right)|}^2}{n\left(\mathbf{r}\right)}-\frac{{\left|\nabla n\left(\mathbf{r}\right)\right|}^2}{8n{\left(\mathbf{r}\right)}^2}.$$(25) On the plane, the KS HOMO ψ_H has a node. It has the same behaviour in the neighbourhood of the nodal plane as the first Dyson orbital d₀, ${\psi }_H\sim f_H(cos\theta )R_H\left(r\right)e^{-\sqrt{2I_0}r}$, with f_H(0) = 0, and f′_H(0) ≠ 0. Note that asymptotically the density very close to the HNP is determined by its slowest decaying part $|d_0\left(\right|\mathbf{r}|\to \infty ,z={\delta \left)\right|}^2$. But in the KS representation, it is determined by the HOMO, $|{\psi }_H\left(\right|\mathbf{r}|\to \infty ,z={\delta \left)\right|}^2$. We, therefore, expect the first Dyson orbital and the KS HOMO to have identical behaviour at the nodal plane, i.e. f_H(0) = f(0) = 0, f′_H(0) = f′(0) and R_H(r) = R(r) for r → ∞. For the HOMO−1, we have the asymptotic behaviour ${\psi }_{H-1}\sim e^{-\sqrt{2I_1}{r}_p}$. We then obtain for the asymptotic behaviour of $v_s^{\text{kin}}$ an expression analogous to Equation (24(24) $$v^{\text{kin}}(r_p\to \infty )\to \frac{1}{2}{f}^{\text{'}}{\left(0\right)}^2\frac{R^2}{r^2}{e}^{2(\sqrt{2I_1}-\sqrt{2I_0})r},$$(24) ), and (26) $$\begin{array}{c}\hfill v_c^{\text{kin}}(r_p\to \infty )=v^{\text{kin}}(r_p\to \infty )-v_s^{\text{kin}}(r_p\to \infty )\\ \hfill \to (f^{\text{'}}{\left(0\right)}^2-f_H^{\text{'}}{\left(0\right)}^2)\frac{R^2}{2r^2}{e}^{2(\sqrt{2I_1}-\sqrt{2I_0})r_p},\end{array}$$(26) Divergence of $v_c^{\text{kin}}$ does not occur if, as anticipated, f′(0)² − f′_H(0)² = 0, which requires perfect similarity between the first Dyson orbital and the KS HOMO at the nodal plane.

We conclude that in Case 1, under rather mild conditions on similar behaviour of the KS HOMO and the first Dyson orbital d₀ at the HNP, the KS potential indeed has the simple, uniform −1/r asymptotic behaviour that is generally assumed. We have not rigorously proven that a small positive or negative asymptotic constant cannot exist in the potential. At this point, however, we feel that postulating such a constant in the present Case 1 is not plausible and would require convincing proof.

5. Case 2: The density decay on HNP is exponentially the same as everywhere, although polynomially slower

5.1. Asymptotic behaviour of the density

We now investigate the possibility that the Dyson orbital d₁ inherits, through Equation (9(9) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)\right)d_i\left(\mathbf{x}\right)+\sum _{k=0}^{\infty }{X}_{ik}\left(\mathbf{x}\right)d_k\left(\mathbf{x}\right)=-I_i{d}_i\left(\mathbf{x}\right).\end{array}$$(9) ), the slow decay $e^{-\sqrt{2I_0}{r}_p}$ from the first Dyson orbital d₀. In that case, the exact density would not have slower exponential decay on the HNP than elsewhere. We have observed that the term coupling d₁ to d₀ in the eigenvalue Equation (13(13) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{ext}\right)d_1+X_{11}{d}_1+X_{10}{d}_0=-I_1{d}_1\end{array}$$(13) ) for d₁ can be written as kcos ²θ/r². It may happen that k ≠ 0, which is Case 2 discussed in this section. The fact that X₁₀d₀ is zero in the HNP (because d₀ is zero there and the prefactor as well) does not preclude coupling of d₁ to d₀. With $d_0\sim rcos\theta e^{-\sqrt{2I_0}r}$, we have an inhomogeneous term $(k{cos}^2\theta /r)e^{-\sqrt{2I_0}r}$ in (13(13) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{ext}\right)d_1+X_{11}{d}_1+X_{10}{d}_0=-I_1{d}_1\end{array}$$(13) ). Clearly, (13(13) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{ext}\right)d_1+X_{11}{d}_1+X_{10}{d}_0=-I_1{d}_1\end{array}$$(13) ) can only be obeyed if this term is cancelled by an equal term coming from −(1/2)∇²d₁. This can be provided by a term in d₁ proportional to $(k{cos}^2\theta /r)e^{-\sqrt{2I_0}r}\equiv f\left(\theta \right)e^{-\sqrt{2I_0}r}/r$. Equation (13(13) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{ext}\right)d_1+X_{11}{d}_1+X_{10}{d}_0=-I_1{d}_1\end{array}$$(13) ) becomes an identity in terms $e^{-\sqrt{2I_0}r}/r$ if (27) $$\begin{array}{c}\hfill \frac{e^{-\sqrt{2I_0}r}}{r}\left(-I_{0}f\left(\theta \right)+k{cos}^2\theta \right)=-I_{1}f\left(\theta \right)\frac{e^{-\sqrt{2I_0}r}}{r}\end{array}$$(27) yielding (28) $$f\left(\theta \right)=\frac{k}{I_0-I_1}{cos}^2\theta$$(28) The presence of this $e^{-\sqrt{2I_0}r}/r$ term does not yet change the behaviour in the HNP because f(π/2) = 0. However, this term in turn necessitates a $e^{-\sqrt{2I_0}r}/r^2$ term, which again has a zero prefactor in the HNP. But the $e^{-\sqrt{2I_0}r}/r$ and $e^{-\sqrt{2I_0}r}/r^2$ terms necessitate next a $e^{-\sqrt{2I_0}r}/r^3$ term, with a prefactor that does not become zero in the HNP and changes the asymptotic behaviour of d₁ and of the density in that plane, (29) $$\begin{array}{ccc}\hfill d_1& =& e^{-\sqrt{2I_1}r}-\frac{k}{I_1-I_0}{cos}^2\theta \frac{e^{-\sqrt{2I_0}r}}{r}\hfill \\ & & +\frac{Qk}{{(I_1-I_0)}^2}{cos}^2\theta \frac{e^{-\sqrt{2I_0}r}}{r^2}+C\left(\theta \right)\frac{e^{-\sqrt{2I_0}r}}{r^3}+\cdots \hfill \\ \hfill C\left(\theta \right)& =& \frac{Qk(\sqrt{2I_0}+Q)}{{(I_1-I_0)}^3}{cos}^2\theta \hfill \\ & & +\frac{k}{{(I_1-I_0)}^2}(2{cos}^2\theta -{sin}^2\theta )\hfill \end{array}$$(29) We note that the $e^{-\sqrt{2I_0}r}/r^3$ term in (29(29) $$\begin{array}{ccc}\hfill d_1& =& e^{-\sqrt{2I_1}r}-\frac{k}{I_1-I_0}{cos}^2\theta \frac{e^{-\sqrt{2I_0}r}}{r}\hfill \\ & & +\frac{Qk}{{(I_1-I_0)}^2}{cos}^2\theta \frac{e^{-\sqrt{2I_0}r}}{r^2}+C\left(\theta \right)\frac{e^{-\sqrt{2I_0}r}}{r^3}+\cdots \hfill \\ \hfill C\left(\theta \right)& =& \frac{Qk(\sqrt{2I_0}+Q)}{{(I_1-I_0)}^3}{cos}^2\theta \hfill \\ & & +\frac{k}{{(I_1-I_0)}^2}(2{cos}^2\theta -{sin}^2\theta )\hfill \end{array}$$(29) ) has a non-zero part in the plane (the sin ²θ term). Therefore, when k ≠ 0, the density in the plane has the same exponential decay as elsewhere. The special circumstance of a HNP shows up in an asymptotic decay of the density by a factor 1/r⁸ faster than the decay of the leading contribution |d₀(r → ∞)|² in other directions.

5.2. Asymptotic behaviour of the Kohn–Sham potential

If the density has the slow decay according to $e^{-2\sqrt{2I_0}r}$ on the HNP, this has significant consequences for the KS potential. Since the HOMO is zero in the HNP, the slow decay must come from HOMO−1. The HOMO−1 KS orbital has eigenvalue ε_{H − 1} which is rather different from −I₀ (actually ≈−I₁). In every other direction than HNP, the KS potential is assumed to go asymptotically to zero like −1/r, so that the HOMO (and the total density dominated by |ψ_{s, H}|²) will decay correctly according to its eigenvalue ε_H = −I₀. The HOMO−1 will then have a decay $e^{-\sqrt{-2{ϵ}_{H-1}}r}\approx e^{-\sqrt{2I_1}r}$, different from $e^{-\sqrt{2I_0}r}$, in every other direction, even arbitrarily close to the HNP. In the one-electron Schrödinger equation for the KS orbitals, there is only a local potential, there is no coupling to other orbitals that could modify the asymptotic behaviour. The asymptotic $e^{-2\sqrt{2I_0}r}$ density decay in HNP must come from HOMO−1, which then has to switch its ‘fast’ decay $e^{-\sqrt{-2{ϵ}_{H-1}}r}$ outside the HNP to the ‘slow’ decay $e^{-\sqrt{2I_0}r}$ on the HNP.

This raises the following question: which properties must the KS potential have in order to make this special behaviour of HOMO−1 possible? It is commonly assumed that the asymptotic decay is just governed by a constant in the potential. Thus, Wu et al. [6] proposed that the decay $e^{-\sqrt{2I_0}r}$ of HOMO−1 in the HNP will be achieved by an appropriate constant asymptotic value of the KS potential in the HNP such that $-\sqrt{2(-{ϵ}_{H-1}+v_{\mathrm{Hxc}}(r_p\to \infty ))}{r}_p=-\sqrt{2I_0}{r}_p$, i.e. v_Hxc(r_p → ∞) should tend to the negative constant I₀ + ε_{H − 1} = −(ε_H − ε_{H − 1}). However, the switching of asymptotic behaviour generates important derivatives perpendicular to the plane, which cannot be neglected, see Section 2. Calculations on real molecular systems that are numerically exact or sufficiently accurate in the asymptotic region to demonstrate the asymptotic behaviour of the KS potential are not possible. Some insight in the peculiar features that this requirement might introduce in the KS potential can be gleaned from simple analytical models. A simple analytical function which has the required decay (fast outside the plane, slow on the plane) may be written like (30) $$\begin{array}{cc}& {\psi }_{s,H-1}(r\to \infty ,\theta )\to \frac{e^{-\sqrt{2I_1}r}}{f(cos\theta )+C{e}^{-(\sqrt{2I_1}-\sqrt{2I_0})r}}\hfill \\ & =\frac{1}{e^{\sqrt{2I_1}r}f(cos\theta )+C{e}^{-\sqrt{2I_0}r}}\hfill \end{array}$$(30) With f(0) = 0, we see that as long as θ ≠ π/2 (31) $${\psi }_{s,H-1}\left(r\to \infty ,\theta \ne \frac{\pi }{2}\right)\to \frac{e^{-\sqrt{2I_1}r}}{f(cos\theta )}$$(31) but on the plane (32) $${\psi }_{s,H-1}\left(r\to \infty ,\theta =\frac{\pi }{2}\right)\to \frac{1}{C}{e}^{-\sqrt{2I_0}r}$$(32) The asymptotics of the KS potential can be calculated directly from the eigenvalue equation obeyed by ψ_{s, H − 1} (33) $$v_s\left(\mathbf{r}\right)={ϵ}_{H-1}+\frac{1}{2}\frac{{\nabla }^2{\psi }_{s,H-1}}{{\psi }_{s,H-1}}$$(33) Of course, the results for the asymptotic behaviour of v_s depend on the chosen function f(cos θ) in our example. Different choices for the function f that determines the θ-dependence of the switching from the decay outside the plane to the different decay on the HNP lead to totally different behaviour for the asymptotics of the potential in HNP, all with strong irregularities at or close to the nodal plane:

• for f = cos θ: $v({\mathbf{r}}_p\to \infty )\to +\infty$ (but not applicable since ψ_{s, H − 1} should be symmetric with respect to the plane);

• for f = cos ²θ: $v({\mathbf{r}}_p\to \infty )\to -\infty$ (and has also large positive ridges around the plane, at θ = π/2 ± Δ, with Δ decreasing as r → ∞);

• for f = cos ⁴θ: $v({\mathbf{r}}_p\to \infty )\to -(I_1-I_0)$ (but has large positive ridges around the plane, at θ = π/2 ± Δ, with Δ decreasing as r → ∞);

• for $f=exp\left(\frac{-1}{{cos}^2\theta }\right)$ : $v({\mathbf{r}}_p\to \infty )\to -(I_1-I_0)$ (but has positive ridges around the plane, at θ = π/2 ± Δ, with Δ decreasing as r → ∞).

It can be seen analytically that the negative constant in the plane can arise if the first and second derivatives of f (with respect to θ) are zero, which is the case with the last two choices for f. The behaviour of v_s for f = cos ⁴θ is illustrated in Figure 3 where v_s is plotted for the choice $\sqrt{2I_1}=3.0$ and $\sqrt{2I_0}=2.0$, (34) $${\psi }_{s,H-1}(r,\theta )=\frac{e^{-3r}}{{cos}^4\theta +e^{-2r}}$$(34) We see in Figure 3 that although on the plane the potential goes to a negative constant, its most prominent features are actually very high positive peaks close to the plane. The region delimited by the peaks becomes very narrow (the peaks become true spikes) when r → ∞, while their height grows exponentially with r.

Figure 3. The potential v(r, θ) that generates the orbital with the asymptotic behaviour of Equation (34(34) $${\psi }_{s,H-1}(r,\theta )=\frac{e^{-3r}}{{cos}^4\theta +e^{-2r}}$$(34) ), shown as a function of the azimuthal angle $\theta =arccos(z/r)$, for different values of r. We see that although on the plane the potential goes to a negative constant, its most prominent features are actually very high positive peaks close to the plane. The region delimited by the peaks becomes very narrow (the peaks become true spikes) when r → ∞, while their height grows exponentially with r.

With the choice f(cos θ) = cos ²θ, the slowly decaying behaviour of ψ_{s, H − 1} in the nodal pane is approached somewhat more steeply than with f(cos θ) = cos ⁴θ. The KS potential now goes to −∞ on the plane for r → ∞, as clearly shown in Figure 4. Notice that this diverging behaviour of the potential at infinity in the HNP is compatible with a fairly regular analytic form of the eigenfunction, as in Equation (30(30) $$\begin{array}{cc}& {\psi }_{s,H-1}(r\to \infty ,\theta )\to \frac{e^{-\sqrt{2I_1}r}}{f(cos\theta )+C{e}^{-(\sqrt{2I_1}-\sqrt{2I_0})r}}\hfill \\ & =\frac{1}{e^{\sqrt{2I_1}r}f(cos\theta )+C{e}^{-\sqrt{2I_0}r}}\hfill \end{array}$$(30) ). In fact, the negative exponential divergence of the potential is cancelled in the KS one-electron equation by a positive divergence from the angular part of −(1/2)∇²ψ_{s, H − 1}, as also discussed in Section 2.

Figure 4. The potential v(r, θ) that generates the orbital with the asymptotic behaviour of Equation (30(30) $$\begin{array}{cc}& {\psi }_{s,H-1}(r\to \infty ,\theta )\to \frac{e^{-\sqrt{2I_1}r}}{f(cos\theta )+C{e}^{-(\sqrt{2I_1}-\sqrt{2I_0})r}}\hfill \\ & =\frac{1}{e^{\sqrt{2I_1}r}f(cos\theta )+C{e}^{-\sqrt{2I_0}r}}\hfill \end{array}$$(30) ) with f(cos θ) = cos²θ, shown as a function of r and of the azimuthal angle $\theta =arccos(z/r)$. We see that as r increases, the potential tends exponentially fast to −∞ on the nodal plane, corresponding to θ = π/2.

We next choose a function that is flatter around θ = π/2 than cos ⁴θ. Consider the function f(cos θ) = exp(− 1/cos ²θ) with all derivatives tending to zero when the plane is approached, so that the ridge in ψ_{s, H − 1} at any finite r is rather broad, narrowing only very slowly with increasing r. Now, at finite values of r, a more extended region is obtained where the potential goes to minus the constant (I₁ − I₀), but the pattern of diverging positive and negative peaks is again observed at both sides of the plane, see Figure 5.

Figure 5. The potential v(r, θ) that generates the orbital with the asymptotic behaviour of Equation (30(30) $$\begin{array}{cc}& {\psi }_{s,H-1}(r\to \infty ,\theta )\to \frac{e^{-\sqrt{2I_1}r}}{f(cos\theta )+C{e}^{-(\sqrt{2I_1}-\sqrt{2I_0})r}}\hfill \\ & =\frac{1}{e^{\sqrt{2I_1}r}f(cos\theta )+C{e}^{-\sqrt{2I_0}r}}\hfill \end{array}$$(30) ) with f(cos θ) = exp(− 1/cos ²θ), shown as a function of the azimuthal angle $\theta =arccos(z/r)$, for different values of r. We see that now the region close to the plane where the potential goes to a negative constant is larger than the one of Figure 3. Yet, this region eventually shrinks as r → ∞, similarly to the case of Figure 3.

We conclude that the situation that the HOMO−1 KS orbital has fast decay according to $e^{-\sqrt{2I_1}r}$ in every direction, however, close to the HNP, but a slow decay $e^{-\sqrt{2I_0}r}$ on the plane, can be handled by special features of the KS potential. These features are, however, very irregular and will be very hard to represent properly in numerical approaches. Of course, these features will also affect the other KS orbitals, so the requirement that the exact total density be reproduced as the sum of squares of all occupied KS orbitals may induce further changes in the behaviour of the potential. Note that we have assumed here that the exponential decay of the HOMO is not changed by the precise behaviour of the KS potential in HNP because this is a nodal plane for the HOMO.

6. Summary and conclusions

It is known that there is an intimate relation between the asymptotic decay of the electron density and the first ionisation potential of an atom or molecule. We have investigated what this relation could be when there is a nodal plane in the KS HOMO. Since the behaviour of the exact density is crucial, we have invoked the expansion of the exact density in terms of squares of Dyson orbitals. The Dyson orbitals of primary ion states (those that can be associated with a simple orbital ionisation) are close to KS orbitals [22] (the orbitals of a non-interacting electron system, like the KS system, are the Dyson orbitals of that system). The analysis shows that one can distinguish two cases. When there would be, in the eigenvalue Equation (9(9) $$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}\left(\mathbf{r}\right)\right)d_i\left(\mathbf{x}\right)+\sum _{k=0}^{\infty }{X}_{ik}\left(\mathbf{x}\right)d_k\left(\mathbf{x}\right)=-I_i{d}_i\left(\mathbf{x}\right).\end{array}$$(9) ) obeyed by the second Dyson orbital, no coupling with the first Dyson orbital (with the same nodal plane as the KS HOMO), the second Dyson orbital would have asymptotic decay according to its eigenvalue, the second ionisation potential I₁. The exact density would decay in the nodal plane like the square of the second Dyson orbital, $e^{-2\sqrt{2I_1}r}$, although everywhere else according to the first ionisation energy, $e^{-2\sqrt{2I_0}r}$. This situation would allow for a perfectly regular KS potential which would decay uniformly (in all directions) like −1/r. However, lack of coupling of the second Dyson orbital with the first requires the integral k of Equation (16(16) $$\begin{array}{cc}& X_{10}(\mathbf{r}\to \infty )=\frac{cos\theta }{r^2}\int {\Psi }_1^{N-1}{({\mathbf{x}}_2\cdots {\mathbf{x}}_N)}^{*}\hfill \\ & \left(\sum _{j>1}{z}_j\right){\Psi }_0^{N-1}({\mathbf{x}}_2\cdots {\mathbf{x}}_N)d{\mathbf{x}}_2\cdots d{\mathbf{x}}_n=\frac{kcos\theta }{r^2}\hfill \end{array}$$(16) ) to be zero. This may happen in special cases, for instance, for symmetry reasons, but will not be true in general. Then, (Case 2) the density decay will be exponentially the same everywhere, although in the nodal plane it would be polynomially slower by a factor r⁻⁸. The KS HOMO−1 has to take care of this slow exponential decay in the HNP (the HOMO being zero there), while the HOMO−1 has decay according to its eigenvalue ε_{H − 1} ≈ −I₁ in all other directions. Such behaviour of the HOMO−1 puts very special demands on the KS potential. Apart from the negative constant (or perhaps −∞) to which it should tend asymptotically in the nodal plane, it will have to have strong, asymptotically increasing, oscillations around that plane.

We should stress that many of our arguments are based on plausibility and examples, rather than on rigorous mathematical proofs. Thus, we hope that the comprehensive study presented here will trigger interest in developing a true rigorous basis for our findings.

Acknowledgments

It is our pleasure to dedicate this paper to Andreas Savin, who always enjoyed to discover strange features in the Kohn–Sham potential and has been a pioneer in asking fundamental questions in exact DFT.

We thank the Netherlands Science Foundation NWO for a visitors grant for Tamás Gál and a Vidi grant for Paola Gori-Giorgi, and the WCU (World Class University) program of the Korea Science and Engineering Foundation (Project No. R32-2008-000-10180-0) for support. Paola Gori-Giorgi acknowledges useful discussions with A. Görling and S. Kümmel.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The Netherlands Science Foundation NWO [grant number 700.59.428]; The WCU (World Class University) program: Korea Science and Engineering Foundation [Project Number R32-2008-000-10180-0].

Appendix. Collapse of the KS conditional amplitude ${\Phi }_s(2\cdots N|{\mathbf{x}}_p)$ to the first excited KS ion state for r_p → ∞

For a non-interacting particle system, the conditional amplitude

(A1) $$\begin{array}{c}\hfill {\Phi }_s({\mathbf{x}}_2\cdots {\mathbf{x}}_N|{\mathbf{x}}_1)=\frac{{\Psi }_{s,0}({\mathbf{x}}_1\cdots {\mathbf{x}}_N)}{\sqrt{n\left({\mathbf{x}}_1\right)}}\end{array}$$(A1) ‘collapses to the second ion state for the reference point ${\mathbf{x}}_1$ going to infinity in the HNP. To see this, we consider the determinantal ground state Ψ_{s, 0} and first let the points ${\mathbf{x}}_1$ go to infinity in a direction where the HOMO is non-zero. The only important contribution to $\sqrt{n\left({\mathbf{x}}_1\right)}$ for ${\mathbf{x}}_1\to \infty$ will be ${\phi }_H\left({\mathbf{x}}_1\right)$. In the numerator of the conditional amplitude, Equation (A1(A1) $$\begin{array}{c}\hfill {\Phi }_s({\mathbf{x}}_2\cdots {\mathbf{x}}_N|{\mathbf{x}}_1)=\frac{{\Psi }_{s,0}({\mathbf{x}}_1\cdots {\mathbf{x}}_N)}{\sqrt{n\left({\mathbf{x}}_1\right)}}\end{array}$$(A1) ), we can expand the determinant. Every term where ${\mathbf{x}}_1$ is in another orbital than φ_H will be negligible. In the remaining terms, the factor ${\phi }_H\left({\mathbf{x}}_1\right)$ in the numerator cancels against the same factor in the denominator. We are left with a determinantal wavefunction with only the other orbitals, i.e. with the first ion state where the HOMO has been removed.

Now, suppose the points ${\mathbf{x}}_1$ are in the HOMO nodal plane. Expanding the determinant, all the terms with ${\mathbf{x}}_1$ in the HOMO are zero. The terms with ${\mathbf{x}}_1$ in another orbital than the HOMO−1 will be negligible. So, we retain the ones with ${\mathbf{x}}_1$ in HOMO−1. The decay of the density in HNP is governed by HOMO−1, i.e. $\sqrt{n\left({\mathbf{x}}_1\right)}$ for ${\mathbf{x}}_1\to \infty$ will be ${\phi }_{H-1}\left({\mathbf{x}}_1\right)$. So, we will have cancellation of φ_{H − 1} on the numerator and denominator. We are left with a determinant in which φ_{H − 1} has been crossed out, i.e. the conditional amplitude has collapsed to the second ion state.