Asymptotic nodal planes in the electron density and the potential in the effective equation for the square root of the density

It is known that the asymptotic decay (|r|→∞) of the electron density n(r) outside a molecule is informative about its first ionization potential I0. It has recently become clear that the special circumstance that the Kohn–Sham (KS) highest-occupied molecular orbital (HOMO) has a nodal plane that extends to infinity may give rise to different cases for the asymptotic behavior of the exact density and of the exact KS potential [P. Gori-Giorgi et al., Mol. Phys. 114, 1086 (2016)]. Here we investigate the consequences of such a HOMO nodal plane for the effective potential in the Schrödinger-like equation for the square root of the density, showing that for atoms and molecules it will usually diverge asymptotically on the plane, either exponentially or polynomially, depending on the coupling between Dyson orbitals. We also analyze the issue in the external harmonic potential, reporting an example of an exact analytic density for a fully interacting system that exhibits a different asymptotic behavior on the nodal plane.


Introduction
Both the (square root of the) density and the KS orbitals ψ k (r) obey Schrödinger-type equations, − 1 2 ∇ 2 + v ext (r) + v eff (r) n(r) = −I 0 n(r), (1) where the sum of the external and the Hartree-exchangecorrelation potentials constitutes the KS potential, v s = v ext + v Hxc . The eigenvalue in equation (1) [1,2] is the first ionization potential, I 0 = E N −1 0 − E N 0 , and the occupied KS orbitals reproduce the density, N k |ψ k (r)| 2 = n(r). For the derivation of equation (1) it is essential to assume that the ground-state, interacting, N -electron wavefunction is real [2]. When this is not the case, an additional vector potential appears in the left-hand-side of equation (1), as in the exact factorization approach put forward by Gross and coworkers [3,4]. In what follows we only focus on the case in which the ground-state wavefunction is real, leaving the interesting investigation of the complex case to future work. In molecules, the external potential v ext (r) goes to zero at large distance like −Z/r, with Z representing the total charge of all nuclei and r the distance from the barycenter of nuclear charge. In this case, according to equations (1) and (2), the asymptotic (|r| → ∞) decay of n(r) and ψ k (r) is (with r = |r|) Both the effective potential v eff (r) for n(r) and the Hartree-exchange-correlation potential v Hxc (r) had been thought, until recently, to go to zero asymptotically. However, if a nodal plane extending to infinity is present in the Kohn-Sham highest-occupied molecular orbital (HOMO), a special behavior on this plane may result. Proposals for the asymptotic behavior of the exact v Hxc (r) have ranged from tending to a positive constant in the HOMO nodal plane (HNP) [5][6][7][8] to a negative constant [9]. The main question is whether the exact, interacting, density has a different decay on the HNP or whether this different decay is only a feature of a single-particle description, as this is what ultimately determines the exact v Hxc (r). A comprehensive investigation of this question has recently been performed for the exact density and the corresponding Hartree-exchange-correlation potential [10], as well as for approximations like GGA's and metaGGA's [11]. For the exact density and KS potential, the analysis of the coupled equations for the Dyson orbitals d i (r) did not allow for a final unique answer, but showed that two different scenarios are compatible with the structure of these equations. The simplest case might be a KS potential uniformly decaying like −1/r, even when the electron density has a different asymptotic decay in the HNP (namely as ∼ exp[−2 √ 2I 1 r p ], with I 1 the second vertical ionization potential and r p a direction belonging to the HNP) than elsewhere (where it is known to have asymptotics ∼ exp[−2 √ 2I 0 r]). However, there are also cases where the density has exponential decay in the HNP according to I 0 , like everywhere else. This will arise if the second Dyson orbital d 1 inherits this asymptotic behavior from the first Dyson orbital d 0 through angular coupling, and then necessarily the KS HOMO−1 will have to provide this same asymptotics in the HNP, since the KS HOMO does not contribute there. In that case the KS potential will exhibit rather special behavior [10] in order to induce the asymptotic decay in the HOMO−1 orbital different from the one according to its eigenvalue (which is typically close to the second ionization potential I 1 ).
In the present paper we further investigate this issue focussing on the effective potential v eff (r) appearing in equation (1), which is related to the functional derivative of the von Weizsäcker kinetic energy functional [12], via the relation The functional T W [n] is also often used in the construction of orbital-free kinetic energy functionals (see, e.g., Refs. [13][14][15][16]). The corresponding energy density, τ W (r) = |∇n(r)| 2 8n(r) , plays a crucial role for metaGGA functionals, where it is used to detect one-electron and iso-orbital regions (see, e.g., Refs. [17][18][19][20][21][22]). Interestingly, even in regions that one traditionally would consider to be oneelectron regions, the ratio τ W /τ may differ from 1 when the system has an orbital nodal plane, see discussion in Section 4 of reference [23].
As mentioned, the main open issue that ultimately determines the behavior of the exact KS and effective potential v eff (r) is whether a density coming from the wavefunction of a fully interacting system has or does not have a different asymptotic decay on the HNP. To further shed light on this open question, here we also report an interacting case that can be solved analytically (two spinpolarized electrons in the harmonic external potential), showing that its density displays different asymptotic decay on the HNP, and thus a different behavior of v eff (r), and discussing the implications for systems bound by the Coulomb potential.
The paper is organized as follows. In Section 2 we will review some quantities needed in the discussion of the asymptotic behavior of the exact density. In particular, the definition of Dyson orbitals and the expansion of the exact density in a sum over the squares of the Dyson orbitals are relevant. Importantly, the asymptotic decay of the exact interacting density may be different in different directions: if there is a KS HOMO nodal plane, it may be inferred that also the first Dyson orbital (at eigenvalue −I 0 ) will have asymptotically the same nodal plane. In that case the decay n(|r p | → ∞) in that plane can be different (faster, according to the eigenvalue I 1 of the second Dyson orbital) than the decay outside the plane. But it might also happen that the second Dyson orbital inherits the slower decay on the plane from the first one, through the coupled equations (9) for Dyson orbitals. In Section 3 we give theoretical expressions for v eff (r), basically relating v eff (r) to wavefunction related quantities (such as the Dyson orbitals). We will recall that the KS potential can be expressed with the help of similar ingredients, with input from the KS independent particle wavefunction. The behavior of v eff (r) is then highlighted in Sections 4 and 5 using examples from both a Coulombic external potential −Z/r and a harmonic external potential 1 2 ω 2 r 2 . The latter affords exact solutions, including electron correlation, for specific values of ω [24], which will be used in Section 5.2. It is established that the potential v eff (r) of equations (1) and (6), while normally going to zero asymptotically, exhibits a different asymptotic behavior in directions where the density decays differently: in the KS HNP it will usually diverge either exponentially or polynomially. Conclusions are given in Section 6.

Asymptotic behavior of the exact density
We briefly review a few aspects of the asymptotic behavior of the exact density based on the analysis of reference [10]. It is known that the asymptotic behavior of the exact density of a molecule is related to its ionization energy [25]. The relation of the density to ion states can be made explicit with the so-called Dyson expansion of the wavefunction where the Ψ N −1 i are the exact (N − 1)-electron states and x = r, s, with s =↑ or ↓. Each state of the ion is associated with a one-particle wavefunction, its Dyson orbital. The sum over i goes over both the spin-↑ and spin-↓ Dyson orbitals. If e.g. x = (r, s =↑) then only the spin-↑ Dyson orbitals are nonzero at x and contribute to n(x) = n(r, ↑) (= 1 2 n(r) in closed shell systems). The Dyson orbitals constitute a nonorthogonal nonnormal, in general linearly dependent set. The Dyson orbitals are, however, not completely esoteric objects. In an independent particle system with a determinantal ground state wavefunction, such as the KS electrons, it follows from (7) that the Dyson orbitals are just the occupied orbitals (in this case there is only a finite number of nonzero Dyson orbitals). The expression of the density in terms of squares of Dyson orbitals is then equivalent to the KS expansion of the density in squares of KS orbitals.
We define the conditional amplitude Φ(2 · · · N ; x) [26] and associated quantities, the conditional density n cond (x 2 |x) and conditional potential v cond (x), Φ(2 · · · N ; x) is a normalized (N − 1)-electron wavefunction depending parametrically on the position x. Its square describes the probability distribution of electrons at positions 2 · · · N when one electron is known to be at x. Its associated one-electron density n cond (x 2 |x) is the density of the other electrons at position x 2 when one electron is at x, which is the normal one-electron density n(x 2 ) plus the full exchange-correlation hole surrounding position x. Projecting the Schrödinger equation (2 · · · N ) and using the expansion of equation (7) one obtains the usual equations for the Dyson orbitals, Katriel and Davidson (KD) [25] pointed out that, due to the coupling integrals the exponential decay of the coupled Dyson orbitals will be the same, cf. reference [27] for the case of Hartree-Fock orbitals. The first Dyson orbital d 0 will have exponential decay ∼e − √ 2I0 r multiplied by a factor r β with β = (Z − N + 1)/ √ 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 nonzero X i0 with the first Dyson orbital will have decay r β−L * e − √ 2I0 r with L * ≥ 2. Dyson orbitals that are not connected to d 0 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), KD have concluded that, if the density decays for |r| → ∞ as the most slowly decaying term |d 0 (x)| 2 , its exponential decay (we do not write here the polynomial prefactor) would be ∼e −2 √ 2I0 r . Levy et al. (LPS) [2] proved this exponential decay in a different way, thereby showing that the leading term is not overruled by the infinite sum of the faster decaying terms in equation (7). The result is generally accepted. It has been realized [5,7,9,28] that there may be special cases where the asymptotic behavior of the density is different in some directions than the one of equation (11). This may happen, for instance, when there is a symmetry plane in the system (as in many π systems, like ethylene and benzene), but 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, if the KS HOMO has a nodal plane, the ground state KS wavefunction (and very likely also the exact ground state wavefunction) corresponding to a closed shell configuration is totally symmetric with respect to that plane, while the first ion state 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 r p in the HNP the conditional amplitude Φ(2 · · · N |x p ) will be symmetric with respect to the plane. Therefore, the matrix element Ψ N −1 0 |Φ(2 · · · N |x p ) 2...N will vanish, so that the first Dyson orbital is zero in the plane: In fact, d 0 is antisymmetric with respect to the plane. When we consider the asymptotic behavior of higher Dyson orbitals, it is clear that with d 0 (x p ) = 0, the coupling to d 0 in equation (9) for points in the HNP at first sight seems to be zero for any higher Dyson orbital d i>0 (but see below). The decay in the HNP of the second Dyson orbital (and thus of the density) is then not governed by d 0 but by d 1 with asymptotic behavior . This is what has been called Case 1 in reference [10]. It is exemplified by the minimal model for a density employed by Aschebrock et al. [11], with a p z type orbital with (outside the HNP z = 0) slow decay ∼ exp[−α p r] and a lower lying s-type orbital with faster decay ∼ exp[−α s r], α s > α p (cf. our − √ 2I 0 r and − √ 2I 1 r for the exponents of HOMO and HOMO−1 respectively). Reference [11] gives a comprehensive discussion of the shape of exchange potentials obtained as functional derivatives of GGA exchange energy approximations (Armiento and Kümmel [29] and Becke [30]), as well as potential functionals like Becke and Johnson [31] and van Leeuwen and Baerends [32]. In that investigation the minimal model density is fed into the density functionals for the various potentials. Often an exponentially diverging behavior is obtained of the form exp[k(α s − α p )r] (k = 1 or 1/2). Remarkably, the same exponential divergence has been observed [10] for the effective potential for the square root of the density for an exact density like the minimal model. Such a density with different decay in a particular plane than elsewhere has been called Case 1 [10] (fast decay according to I 1 in the plane, slower decay according to I 0 everywhere else). However, it is an important issue whether a true density of Coulombically interacting electrons can have such different exponential decay in different directions. The present authors have argued that an exact density will typically not exhibit this different exponential decay in different directions (although the polynomial prefactor may differ). This has been called Case 2 in [10]. Decay of the HOMO−1 according to I 0 in the HNP leads to rather intricate consequences for the KS potential, which requires very special features to generate a decay of HOMO−1 in HNP according to I 0 and not according to its eigenvalue (which is equal to (or close to) I 1 ). The situation for v eff (r) is, however, simpler than for the KS potential, since it can be related directly to wavefunction quantities, as discussed in Section 3.
3 Asymptotic behavior of the effective potential for √ n The potential v eff (r) of equation (1) can be written in the form [2,33] v For future reference we note that for the exact KS potential an analogous expression holds, The potentials v kin s and v N −1 s depend on the KS independent particle wavefunction Ψ N s and notably its associated conditional amplitude Φ s in exactly the same way as v kin and v N −1 depend on the exact wavefunction and conditional amplitude. In equation (13) v eff (r) is expressed in terms of only wavefunction quantities. LPS [2] stressed that each term in equation (13) is everywhere nonnegative and should tend to zero asymptotically. In fact, v cond (x) [see equation (8)], being the repulsive Coulomb potential of a localized charge distribution of (N − 1) electrons, decays like (N − 1)/r. The third term of v eff , 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 0 . When |r| → ∞ it has been inferred [25] that the conditional amplitude collapses to the ion ground state Ψ N −1 0 (when s =↑ then Φ will collapse to the M S = −1/2 state of the doublet ion), so that v N −1 (|r| → ∞) → 0. The second term, v kin , is manifestly positive and is expected to go to zero asymptotically since the derivative of Φ with respectto r when the reference electron is very far becomes zero (Φ remains constantthe ion ground state -under small change of r at ∞). These expectations are not borne out if there is a HNP, see below.
3.1 Case 1: The density decay on the HNP is governed by the HOMO−1 For points r p in the HNP d 0 (x p ) = 0 because of spatial symmetry. By expanding the conditional amplitude Φ(2 · · · N |x) in terms of the exact N − 1 states, we see that, since on the HNP d 0 = 0 and |d 1 (x p )|(r p → ∞) ∼ n(x p ), while all higher d i decay a factor r −L * faster [25], with L * ≥ 2, the conditional amplitude tends asymptotically on the plane to the first-excited ion state, This is a positive constant appearing in the asymptotics of v eff only on the HNP. The exponential decay of √ n is governed according to equation (1) by exp(− 2(I 0 + v eff (∞))r). The positive asymptotic constant in v eff (r → ∞) looks perfectly in order: this value for v eff (∞) gives precisely the asympotic decay exp(− √ 2I 1 r p ) we have assumed for √ n on the HNP in Case 1. Also the second term in equation (13), v kin , can be nonzero at infinity: when crossing the HNP, the asymptotic conditional amplitude changes from Ψ N −1 0 , to which it collapses for asymptotic points in general directions, to Ψ N −1 1 , to which it collapses for asymptotic points in HNP, see above. So the r-derivative of Φ perpendicular to the plane can be nonzero on the HNP also when |r| → ∞. Its actual value depends on how d 0 (r → r p ) goes to zero when approaching the nodal plane. We have noted that for the determinantal wavefunction of a noninteracting 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 ionization) still are very similar to the Kohn-Sham orbitals: overlaps are typically >0.999 [34]. This agrees with our finding [10] 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 ∼ f (cos θ)R(r)e − √ 2I0 r , with f (0) = 0, and f (0) = 0, as would be the case for a π orbital, which has f R = r cos θ = z. By writing v kin in the form [33,35,36] v kin (r) = 1 2 and using d 1 ∼ e − √ 2I1 r , it is found after some manipulation that showing that v kin can go asymptotically to infinity on the HNP. A simple illustration of this fact is given in the next Section 4 for non-interacting electrons (a Case 1 density). The asymptotically diverging behavior of equation (18) is perfectly compatible with an analytical, well-behaved density. It induces in the density the special behavior in the HNP of Case 1 which is certainly realizable by noninteracting electrons in one-electron states (orbitals): fast decay according to I 1 in HNP coming from HOMO−1, slow decay everywhere else according to I 0 from HOMO. The key point is that when we project equation (1) on the plane, we have to take into account also ∇ 2 √ n in the direction perpendicular to the plane. Usually, the θ and φ derivatives in ∇ 2 are zero when r → ∞, but when there is a HNP this is not the case. Indeed in our example, using spherical coordinates, the − 1 r 2 sin θ ∂ ∂θ sin θ ∂ ∂θ operation on √ n exactly cancels the diverging behavior coming from v kin . This leaves for r p → ∞ just the radial part of the one-electron Schrödinger equation. With the remaining potential −1/r from v ext + v cond , and the I 1 − I 0 constant of v N −1 combined with the eigenvalue −I 0 , √ n acquires the asymptotic decay in the HNP according to I 1 .
For an interacting electron system the density is described by the leading terms in the Dyson expansion, n(x) = |d 0 (x)| 2 + |d 1 (x)| 2 + · · · . Only if there is no coupling of d i>0 to d 0 in equation (9) will d 1 (and the orbitals in the same set) have asymptotics according to I 1 and will this picture for noninteracting electrons also prevail for the interacting electron system. We discuss in the next subsection the Case 2 where such coupling does occur.
Considering the kinetic correlation potential v kin c = v kin − v kin s in the KS potential, we have argued in reference [10] that if the exact density is like the noninteracting (KS) density with a HOMO nodal plane, the behavior of the Dyson orbitals d 0 and d 1 close to HNP should be identical to that of HOMO and HOMO−1, and in v kin c the divergence of v kin is canceled by an equal divergence of v kin s . Then in Case 1 the KS potential will have asymptotically the simple uniform −1/r behavior, compatible with solutions of the KS equations with a HOMO with a nodal plane and a HOMO−1 with uniformly faster decay (the density of the minimal model of Aschebrock et al. [11] is compatible with such a regular KS potential). This is then a consistent picture. However, we have also indicated that the situation where coupling of d 1 to d 0 in equation (9) generates slow decay in d 1 will be prevalent in interacting electron systems, see discussion of Case 2 in next section.
As recalled in equation (6), the potential v eff (r) essentially gives the functional derivative of the von Weizsäcker kinetic energy functional, which, thus, also exhibits the same diverging behavior on the nodal plane in Case 1. Notice that all the spherical harmonics Y m (θ, φ) approach their nodal planes linearly with cos θ, so that the divergence predicted by equation (18) is expected to occur in the general case. This might have consequences for calculations using orbital-free kinetic energy functionals and metaGGA functionals, probably in a way qualitatively similar to the one reported by Aschebrock et al. [11].
3.2 Case 2: The density decay on HNP is exponentially the same as everywhere, although polynomially faster In reference [10] it has been shown that, in Coulombically interacting systems, coupling of some of the d i>0 to d 0 in equation (9) will usually occur (Case 2), and will lead to d 1 (and other orbitals) acquiring on the nodal plane the same slow exponential asymptotic decay (dictated by I 0 ) as d 0 in general directions. The decay of d 1 will then still be polynomially faster on the plane (by 1/r 4 ) (and the correlated density a factor 1/r 8 faster). We have investigated what the asymptotic behavior of v kin will be for such a Case 2 density. Let us consider the essential terms in d 1 responsable for the slow e − √ 2I0 r behavior [10], notably also the term C 1 e − √ 2I0 r /r 3 yielding this slow decay of d 1 on the plane where the constant C 1 is non-zero if f (0) = 0, which usually will occur, since, as discussed in reference [10], in the vast majority of cases we will have f (x) = x 2 . After some manipulation one obtains for the asymptotics of v kin v kin (r → ∞, θ = π 2 ) = q 2 1 r 6 2( So v kin will still be diverging on the nodal plane of the first Dyson orbital d 0 , but the divergence is no longer exponential, but becomes polynomial (like r 6 ). Note that if C 1 → 0, i.e. when d 1 does not have the slow decay on the plane, then we are back in Case 1 above, where v kin diverges more rapidly on the plane, in fact exponentially as in equation (18). One may again verify that the r 6 divergent behavior does not pose any problem in equation (1) for √ n since it is canceled by an opposite divergent term coming from the Laplacian of √ n. A more detailed analysis including all the Dyson orbitals d i>0 that inherit the slow decay on the HNP from d 0 through the same kind of angular coupling does not change qualitatively the conclusion of equation (20), with the caveat that one should always be careful with asymptotic expansions expressed as infinite sums. Equation (20), where i ∈ G denotes the set of all the Dyson orbitals having the same coupling with d 0 as d 1 . This set includes all the Dyson orbitals for which the matrix element k i appearing in equation (16) of reference [10] is nonzero, the constants C i being determined by the matrix element k i divided by (I i − I 0 ) 2 . In highly symmetrical systems, like ethylene and benzene, many k i will be zero by symmetry [10]. Also, one should keep in mind that now the conditional probablity on the HNP does not collapse asymptotically anymore to the first excited state of the ion, but to a superposition of all the ion states with Dyson orbitals that are non-zero on the plane (belonging to the set G having As a consequence, equation (16) does not hold anymore and we have, instead, which reduces to equation (16) when only i = 1 is considered.

External harmonic potential
We also analyze how these conclusions may become different for a different external potential. We take the case of an harmonic external potential, v ext = 1 2 ω 2 r 2 , which has the interesting property that it affords analytic solutions for two Coulombically interacting electrons for specific values of ω [24]. Moreover, just as the Coulombic external potential, it has many applications in physics (quantum dots, cold atoms, plasmas, etc.). Equation (1) is generally valid for any binding external potential, including the harmonic confinement. From the decomposition of v eff (r) of equation (13) we clearly see that also in this case v eff (r) is expected to go asymptotically to zero. In other words, it is only the Coulombic nature of the electron-electron repulsion, determining the conditional amplitude and related quantities, that matters for the asymptotic behavior of v eff (r). It is again a special behavior of the density in certain directions, such as in a HNP, that may induce special behavior of v eff (r) in these directions.
Obviously, the way the information on I 0 is embodied into the asymptotics of the density is different with the harmonic external potential than with the Coulombic external potential of the molecular case. In d dimensions, with v eff (|r| → ∞) → 0, equation (1) for the square root of the density of N electrons confined in an harmonic trap reads, when r → ∞, As well known, the solution has the form e − ω 2 r 2 u(r), where for large r, u(r) ∼ r q , with q ∈ R + . While the gaussian decay only depends on ω, it is now the polynomial prefactor (q) that carries the information on with and hence in three dimensions (d = 3) Notice that, due to the unbounded nature of the potential, the ionization energy cannot be defined as removing the particle to infinity with zero kinetic energy. Removing a particle (ionization) is equivalent to putting the particle with zero kinetic energy at the zero of the harmonic potential well (without interaction with the other particles). The energy of a (N − 1)-particle state in the harmonic potential minus the energy of a N particle state is then negative with −I 0 ≥ 3 2 ω, so that q is positive. The negative ionization energy does not change the derivation of equation (1). The exact solutions that are possible in this case afford an analytical study of the asymptotic behavior of v eff (r) for both a noninteracting and an interacting correlated electron system in the presence of a KS HOMO nodal plane, as reported in Section 5.
We should note, however, that because the ionization information is now only in the polynomial prefactor, Cases 1 and 2 discussed for the external Coulomb potential in Sections 3.1 and 3.2, respectively, can get mixed in the harmonic external potential. The reason is that, as explained in Section 3.2, when the angular coupling between the Dyson orbitals makes d 1 inherit the slower behavior of d 0 on the plane, such behavior is damped by a factor 1/r 4 (which becomes 1/r 8 in the density). This polynomial damping does not prevent d 1 from getting a slower decay on the plane if the difference in the two asymptotic behaviors is exponential, as it happens for the Coulomb external potential case, but can have an important effect when the difference is only polynomial.

Non-interacting electrons in the Coulomb external potential
In order to illustrate the diverging behavior predicted by equation (18), we consider N = 3 non-interacting electrons in the Coulomb external potential v ext (r) = −Z/r in the configuration 1s 2 2p 1 z , which is one of the possible degenerate ground-states. The corresponding wavefunction is antisymmetric with respect to the nodal plane z = 0. We construct the corresponding density and we calculate v eff (r) by inversion, v eff (r) = ∇ 2 n(r) The result is shown in Figure 1 for Z = 3 as a function of r for different values of θ = arccos( z r ) (top panel) and as a function of θ for different values of r (bottom panel). We see that v eff (r) is a smooth function going to zero asymptotically everywhere except on the plane z = 0, where it diverges exponentially according to equation (18).
The diverging behavior we illustrate here for v eff (r) is precisely the same as the one obtained for various approximations (like GGAs) to the KS potential by Aschebrock et al. [11], if a density like the present one is inserted in the corresponding exchange potential expressions. However, such a density (called minimal model in [11]) is generated here by a purely Coulombic potential in the Schrödinger equation (2). It is not clear if the similarities between the present exact v eff (r) potential for this type of √ n (with faster decay on the HNP) and the approximate (GGA) KS potentials obtained with this density are more than accidental.

Electrons in a harmonic external potential
We illustrate here the behavior of √ n and v eff (r) in the harmonic external potential in the two cases of non-interacting and interacting electrons, where we use one of the analytic solutions of Taut [24] for the spin-polarized (triplet) case.

Non-interacting electrons
We first consider again N = 3 non-interacting electrons and we put them in the harmonic potential v ext (r) = 1 2 ω 2 r 2 . With the lowest totally symmetric orbital (s type) doubly occupied, and one electron available for the degenerate p orbitals (configuration s 2 p 1 ), we select among the three degenerate ground states the one with m = 0 so that the HOMO is again a p z orbital. The corresponding v eff (r) is calculated as in the previous section, see (28), and is reported in Figure 2. We see that, on the HNP, v eff (r p → ∞) again does not go to zero, but this time it tends to a constant. It is easy to verify analytically that if, asymptotically, ψ H ∼ r q0 f (cos θ)e − ω 2 r 2 with f (0) = 0, and ψ H−1 ∼ r q1 e − ω 2 r 2 , we have v kin (r p → ∞) ∼ f (0) 2 r 2(q0−q1−1) , For non-interacting electrons we have always E N −1 = ω so that the potential goes to a constant in the plane, For N ≥ 3 interacting electrons, depending on how correlated is the system, we could have E N −1 > ω, and thus a polynomially diverging behavior of v eff on the HNP. Comparison of equations (29) and (18) shows that in the presence of a HOMO nodal plane that extends to infinity the asymptotic behavior of v eff (r) on the plane can depend dramatically on the kind of binding external potential.

N = 2 spin-polarized interacting electrons
We consider now N = 2 spin-polarized interacting electrons (with standard Coulomb 1/r 12 interaction). As well known, the corresponding hamiltonian is separable into center-of-mass R = 1 2 (r 1 + r 2 ) and relative r 12 = r 2 − r 1 coordinates, so that its exact wavefunction reads Ψ N 0 (r 1 , r 2 ) = ξ(R)φ(r 12 ). With spin-polarized electrons, the spatial wavefunction must satisfy Ψ N 0 (r 1 , r 2 ) = −Ψ N 0 (r 2 , r 1 ), which implies that the ground state corresponds to the 12 = 1 spherical harmonic for the relative vector r 12 . We have then 3 degenerate ground-state wavefunctions, and we choose one of them by fixing m 12 = 0: this way, we obtain an interacting density with a symmetry plane like the one encountered in molecules.
For N = 2 there is an infinite set of special values of ω for which q in equation (27) is integer: they correspond to analytical solutions of the interacting hamiltonian [24]. For 12 = 1, ω = 1 4 is one of those. The interacting wave function for the case m 12 = 0 then is equal to with C a normalization constant. The associated density is given by n(r) = C n e − r 2 4 π 3/2 2 2 26 + r 2 − z 2 + z 2 32 + r 2 + 4π r 5 e − r 2 4 −24 r z 2 +8r 3 1+z 2 +2r 5 2+z 2 + √ π erf r 2 24z 2 + r 6 (2 + z 2 ) − 4r 2 2 + 3z 2 with the normalization constant Inserting equation (32) into equation (28) we find that the corresponding v eff (r) has the same asymptotic behavior as observed for non-interacting electrons in the harmonic potential in Figure 2. This is shown in Figure 3, where we compare our v eff (r) for interacting electrons with the one for two non-interacting spin-polarized electrons in the same external potential. We clearly see that, in this case, the asymptotic behavior close to the HNP is exactly the same. This is in agreement with our findings of Sections 3 and 5.1: the behavior close to the HNP is entirely determined by the differences between the ground and the first excited states of the N − 1 state. For N = 2, the N − 1 states are the same for both interacting and non-interacting electrons.
In this case we can also compute analytically the first two Dyson orbitals that can be obtained from