One-electron self-interaction and the asymptotics of the Kohn-Sham potential: an impaired relation

One-electron self-interaction and an incorrect asymptotic behavior of the Kohn-Sham exchange-correlation potential are among the most prominent limitations of many present-day density functionals. However, a one-electron self-interaction-free energy does not necessarily lead to the correct long-range potential. This is here shown explicitly for local hybrid functionals. Furthermore, carefully studying the ratio of the von Weizs\"acker kinetic energy density to the (positive) Kohn-Sham kinetic energy density, $\tau_\mathrm{W}/\tau$, reveals that this ratio, which frequently serves as an iso-orbital indicator and is used to eliminate one-electron self-interaction effects in meta-generalized-gradient approximations and local hybrid functionals, can fail to approach its expected value in the vicinity of orbital nodal planes. This perspective article suggests that the nature and consequences of one-electron self-interaction and some of the strategies for its correction need to be reconsidered.


Density functional approximations and their Kohn-Sham potentials
During the past few decades, Kohn-Sham density-functional theory (DFT) 1,2 has evolved into a standard tool for electronic structure calculations of atoms, molecules and solids. The decisive quantity of DFT is the exchange-correlation (xc) energy functional, E xc , which contains all electronic interactions beyond the classical electrostatic Hartree contribution, E H . E xc in practice has to be approximated, and the approximation used governs the accuracy of a DFT calculation. 3,4 It is one of the puzzles of DFT that explicit density functionals such as the generalized gradient approximations (GGAs) can predict binding energies and bond lengths of complex many-electron systems reliably, but make substantial errors in describing simple one-electron systems. The underlying problem is wellknown as the one-electron ''self-interaction problem'': 5 for the exact functional, E xc + E H will vanish for any one-electron ground-state density because one electron does not interact with itself -but most approximate functionals yield a spurious finite value for this case. Following ref. 5 a functional is considered to be one-electron self-interaction free if it fulfills the condition where n is = |j is (r)| 2 designates a single spin-orbital density.
Self-interaction plays a decisive (although not the only) role in the (un)reliability of density functional theory calculations, and its consequences are particularly pronounced, e.g., in questions of orbital localization, 5-8 ionization processes, [9][10][11][12] charge transfer, [13][14][15] and the interpretability of eigenvalues and orbitals, e.g., as photoemission observables. [16][17][18][19][20][21][22][23] Many of these observables can also be directly related to properties of the Kohn-Sham exchange-correlation potential, which is defined as the functional derivative of the xc energy with respect to the ground-state density n(r), i.e., v xc ðrÞ ¼ dE xc ½n dnðrÞ . It is generally expected that there is a close relation between freedom from self-interaction and xc potential features. The field-counteracting term that is important for obtaining correct response properties is one such feature. 24,25 Another example, and probably the most prominent one, is the long-range asymptotic behavior of the xc potential, 26,27 v xc ðrÞ À ! jrj!1 À 1 jrj : (Hartree units are used here and throughout.) In this perspective article we focus exclusively on Kohn-Sham theory, i.e., on a local multiplicative xc potential, as opposed to orbitalspecific (non-multiplicative) potentials that arise in generalized Kohn-Sham theory 28 and are used in the standard application of hybrid functionals. In the Kohn-Sham approach, the local xc potential models the interaction of one particle with all others and it therefore appears intuitively plausible that a functional that is not self-interaction-free cannot show the correct À1/r long-range asymptotic behavior: as one particle of a finite, overall electrically neutral system ventures out to infinity, it will ''feel'' the hole of charge 1 that it left behind in the total charge. This gives rise to the À1/r potential asymptotics (see, e.g., ref. 4, p. 242 for a more detailed argument along these lines). However, a particle that spuriously self-interacts will ''feel'' itself, and thus not the proper hole. Consequently, the potential will not have the proper long-range decay.
The correct asymptotics of the xc potential has proven to be important for a variety of physical quantities. It plays a prominent role for obtaining stable anions in DFT, it leads to a Rydberg series in the Kohn-Sham eigenvalues and generally to unoccupied eigenvalues of improved interpretability, and as a consequence allows for improved accuracy in the prediction of various response properties. [29][30][31][32][33] The correct asymptotic behavior is also important for the ionization potential (IP) theorem, 26,27,34,35 which states that the negative of the highest occupied Kohn-Sham eigenvalue Àe ho should correspond to the vertical IP, and for developing functionals that allow for approximately predicting IPs from ground-state eigenvalues. 36,37 There have been fruitful attempts to incorporate the correct behavior in the limit |r| -N directly into the xc potential, [38][39][40][41] leading to improvements in the description of some of the aforementioned properties. However, since directly designed potential expressions are typically not functional derivatives of any energy functional, the use of such ''potential only'' approximations is necessarily limited, as discussed, e.g., in detail in ref. 42-44. A functional that combines freedom from self-interaction and the correct asymptotics of the potential is exact exchange (EXX), being defined as the Fock integral evaluated using Kohn-Sham orbitals j is (r), where i labels orbitals in spin channel s: Here, N s is the number of electrons with spin s. Treating exchange exactly with a local Kohn-Sham potential leads to a significant improvement in the quality of Kohn-Sham eigenvalues when compared to (semi-)local functionals. 20,45,46 EXX also tends to increase Kohn-Sham gaps, 8,47-50 leads to a desired particle number discontinuity in static 51 and time-dependent 52 situations, and improves the description of charge transfer, 24,25 dissociation 53 and ionization processes. 52 However, using bare EXX is known for its rather poor description of binding energies and structural properties (see, e.g., ref. 54 and 55, and ref. 4, chapter 2). Adding a (semi-)local correlation term to EXX hardly improves the situation and typically leads to results that are inferior to the ones from (semi-)local functionals. The reason for this failure is the wellknown incompatibility of the fully non-local Fock exchange with a purely (semi-)local correlation term. 56 A class of approximations which has been designed to remedy this incompatibility is that of local hybrid functionals, 57,58 sometimes also called hyper-GGAs. 56 Whereas global hybrid functionals 59-64 mix a constant, fixed fraction of Fock exchange with (semi-)local exchange and correlation, local hybrids replace the fixed fraction by a density dependent local mixing function (LMF). Both types of hybrids originate from the concept of the coupling-constant integration, i.e., the adiabatic connection scheme. 59,65 Global hybrids are successful in modeling the coupling-constant averaged, integrated energy. Local hybrids can go one step further and aim to model the coupling-constant curve itself 66 instead of just the integral. Thus, in contrast to the global hybrid functionals which are used in practical applications of DFT and combine GGA components with a fixed fraction of exact exchange, local hybrids can incorporate full exact exchange and can be fully one-electron self-interaction-free.
An early local hybrid with reduced one-electron selfinteraction error showed promising results for dissociation curves and reaction barriers, but its accuracy for binding energies was limited. 58 A self-consistent implementation of a local hybrid functional was given in ref. 67, and over the years several local hybrids were constructed, using different LMFs and (semi-)local exchange and correlation functionals, [68][69][70][71][72][73][74][75][76] striving to reach greater accuracy by refining the positiondependent mixing of nonlocal and local components. Many of these functionals rely on the concept of an iso-orbital indicator, i.e., a functional that allows one to distinguish regions of space in which the density is dominated by one orbital shape from regions of space where several orbitals of different shapes contribute to the density. The most prominent iso-orbital indicator, which goes back to a long tradition of using kinetic energy densities in density functional construction, [77][78][79] is the ratio of the von Weizsäcker kinetic energy density t W to the positive (as opposed to other possible definitions, see, e.g., ref. 80) Kohn-Sham kinetic energy density t, discussed in detail below.
By using full EXX and an iso-orbital indicator, local hybrids aim at being one electron self-interaction-free and producing a Kohn-Sham potential with the proper long-range asymptotic decay. They are a paradigm class of functionals designed for simultaneously curing both of these two prominent problems of (semi-)local density functionals. In the following, we therefore use the example of a local hybrid functional to shed light on the relation between a functional's self-interaction and its potential asymptotics, as well as the properties of the t W /t indicator. We argue that quite generally a one-electron selfinteraction-free energy does not guarantee the correct longrange potential, and that t W /t loses its indicator ability in the vicinity of nodal planes of the highest-occupied molecular orbital (HOMO).
with e xc ([n];r) denoting the xc energy density per particle. The definition of e xc (r) is not unique and subject to a gauge-dependence. 72 Yet, for local hybrid functionals it has become common to define this energy in the form e lh xc (r) = e ex x (r) + f (r)(e sl x (r) À e ex x (r)) + e sl c (r).
Here, e ex x (r) marks the exchange energy density per particle corresponding to the EXX energy of eqn (3). This non-local term is mixed with (semi-)local exchange and correlation energy densities e sl x (r) and e sl c (r), respectively. The position dependent mixing ratio f (r), which is itself a density functional, marks the LMF.
Often, the LMF is designed in a way that aims at eliminating the one-electron self-interaction error of eqn (1) that is inherent in most density functionals. An established method for reducing self-interaction effects is to detect regions of space where a single Kohn-Sham orbital shape dominates the density (''iso-orbital regions''), and then enforce eqn (1) in these regions. One of the most popular 58,67,70,71,74,76,[80][81][82][83] indicator functions for detecting iso-orbital regions is where t W (r) = |rn(r)| 2 /(8n(r)) denotes the von Weizsäcker kinetic energy density and tðrÞ ¼ 1 2 rj is ðrÞ j j 2 is the positive Kohn-Sham kinetic energy density. In iso-orbital regions, t(r)t W (r) and therefore g(r) -1. In the case of a slowly varying density, t W (r) -0 and, since t(r) remains finite, g(r) -0. This indicator function is typically a decisive ingredient in the LMF, f (r), of local hybrids. With its help one can construct f (r) such that eqn (5) reduces to correct limiting cases, e.g., e sl x (r) + e sl c (r) for slowly varying densities, and e ex x (r) for single orbital regions. The latter case additionally requires that e sl c (r) vanishes in single-orbital regions, a condition that we discuss below.
In the asymptotic limit, |r| -N, the xc energy density for a finite system should be dominated by e ex x (r). When e sl c (r) vanishes sufficiently fast in the asymptotic region (a condition that is usually fulfilled), then is the requirement that one aims at, because it leads to the correct asymptotic limit of the xc energy density per particle e lh xc ðrÞ $ e ex x ðrÞ À ! jrj!1 À 1 2r : (Note the difference to the asymptotic limit of the xc potential, see ref. 38).
Since for a finite system each Kohn-Sham orbital decays exponentially with an exponent set by its eigenvalue, 84 the density is asymptotically dominated by the HOMO density, i.e., it becomes of iso-orbital character. Therefore, g(r) can be used in the construction of the LMF to realize eqn (8).
Considerations of the type discussed above are inherent to many density functional constructions. As a particular example for a local hybrid functional we here use a recently proposed, physically motivated LMF, 83 which reads The function g(r) in the numerator is multiplied by the squared spin polarization z(r) = (n m (r) À n k (r))/(n m (r) + n k (r)), which lets the LMF not only identify iso-orbital regions, but also correctly distinguish between true one-orbital regions, and regions with two identical spin-orbitals. The function g(r) is used in such a way that f t (r) vanishes for one-orbital regions, as required. The use of the reduced density gradient where a 0 is the Bohr radius and F(z(r)) = 1 2 ((1 + z) 2/3 + (1 À z) 2/3 ), in the denominator of f t (r), ensures the correct behavior of E xc under uniform coordinate scaling rgr. 85,86 The density transforms as n g (r) = g 3 n(gr) and as a consequence eqn (9) uses full exact exchange in the sense of ref. 72 The function t 2 (r) is multiplied by a parameter c that we cannot determine, at least presently, from fundamental constraints. It allows for adjustments in the functional ansatz. In the case of slowly varying densities, f t (r) -1 and eqn (5) reduces to its purely (semi-)local components. As an aside we note that this LMF comprises the one of ref. 58 as the special case when c = 0 and z(r) = 1 8 r. We denote this case by f 0 (r), i.e., f 0 ðrÞ ¼ 1 À t W ðrÞ tðrÞ .
For the semi-local exchange we use the LSDA, 87 i.e., e sl x (r) = e LSDA x (r), whereas e sl c ðrÞ ¼ 1 À t W ðrÞ tðrÞ z 2 ðrÞ e LSDA c ðrÞ. The additional multiplication with the numerator of eqn (9) consistently reduces eqn (5) to pure EXX in the one-spin-orbital case, where e LSDA c (r) alone does not vanish. The general questions that we discuss in this perspective article, i.e., whether there is a relation between self-interaction and the xc potential asymptotics and how far the iso-orbital indicator t W /t can be used to enforce freedom from selfinteraction, can be scrutinzed with the local hybrid of eqn (9) as an instructive example.

The Kohn-Sham exchangecorrelation potential of local hybrid functionals
In order to implement local hybrids self-consistently within the Kohn-Sham scheme, one has to find the local multiplicative xc potential corresponding to the energy of eqn (4) and (5). The fact that local hybrids use EXX and typically also t(r) makes them explicitly orbital-dependent. Therefore, the local xc potential must be obtained from the optimized effective potential (OEP) equation (see, e.g., ref. 45, 55 and 88-91). The computational effort can be reduced significantly by employing the approximation of Krieger, Li and Iafrate (KLI). 92 For the local hybrid of eqn (9) it has been shown that the total energy E tot and the highest occupied Kohn-Sham eigenvalue e ho obtained using the KLI approximation agree quite well with the ones from the full OEP. 83 Furthermore, it is a general finding 84 that the KLI approximation does not affect the potential asymptotics to leading order. In the actual calculations presented in the following we therefore always use the KLI approximation.
In the OEP (and KLI) scheme the chain rule for functional derivatives 45 relates the derivative with respect to the density to the derivatives with respect to the orbitals, From the structure of the OEP equation it further follows that to first order i.e., the functional derivative with respect to the HOMO in general determines the potential asymptotics. 84 Therefore, investigating the HOMO functional derivative is the key to determining the asymptotic behavior of an orbital dependent functional's xc potential. When one takes the functional derivative (with respect to the orbital) of a local hybrid one obtains three terms, corresponding to the three addends in eqn (5): Evaluating the asymptotical behavior of each of these three terms for the highest occupied orbital allows one to predict the potential asymptotics.
The first term can be derived directly from eqn (3) and reads This term evaluated for the HOMO indeed provides the correct asymptotic behavior 45 The third term u c-sl is (r), on the other hand, does not contribute to the asymptotics of eqn (16) as it decays exponentially due to its purely (semi-)local nature.
Evaluating the second term on the right-hand side of eqn (14) requires careful consideration. Intuitively, one might expect that an asymptotically vanishing LMF will surpress any asymptotic contribution of this term to the potential. In the following we check this expectation. Details of the underlying calculation for both LMFs used in this work, i.e., f t (r) and f 0 (r), can be found in ref. 83 and in Appendix B, eqn (27), respectively.
By defining P(r) = n(r)e sl x (r) and Q(r) = n(r)e ex x (r) one can write The first two terms consist of (semi-)local components and thus vanish exponentially. Evaluating the third term on the other hand is not as trivial as it contains the non-local quantity Q(r) as well the functional derivative of the LMF with respect to the corresponding Kohn-Sham orbital. For the LMFs addressed in this perspective we find that this term does not contribute to the asymptotics either (see ref. 83 and Appendix B for details).
Thus, only the fourth term in eqn (17) is relevant in the asymptotic limit and therefore The first term in this equation equals Àu exx is (r) of eqn (15), locally multiplied by f (r)/2. Due to eqn (7) it vanishes faster than the leading term of u lh is , which is given in eqn (16). The second term, however, is of a different structure, as it evaluates the LMF under the integral. By considering the HOMO level, its asymptotic limit is This corresponds to a Hartree-like potential caused by the spin-orbital density of the HOMO averaged over all space, with the LMF as a weighting function. Thus, this term gives a finite contribution in the asymptotic limit despite eqn (7). Now, when adding the asymptotically significant components, eqn (16) and (19), for the evaluation of eqn (13), we arrive at v xcs ðrÞ À ! jrj!1 À g s jrj : Here, the parameter g s denotes the reduced slope of the potential asymptotics, which can numerically be extracted from a self-consistent Kohn-Sham calculation via Eqn (21) is a central result of this work, as it demonstrates that a local hybrid of the form of eqn (5) does not lead to the exact asymptotic behavior of the xc potential. Eqn (20) holds for all f (r) that vanish in the asymptotic limit and for which the third term of eqn (17) does not contribute to the asymptotics of the functional derivative u c-nl is (r), i.e., under very general conditions.
Further details of the calculation, specifically regarding the question of the xc potential asymptotics in different spin-channels, are given in Appendix B. The LMF is limited between 0 r f (r) r 1 and therefore the asymptote is bound between 1 2 o g s r 1. Consequently, the exact value g s = 1 can only be reached by setting f (r) = 0 8 r, which corresponds to the trivial case of using EXX, ''as is'' or combined with a purely (semi-)local correlation functional.
A different extreme case, f (r) = 1 8 r, does not, as one could naïvely believe due to eqn (21), lead to g s = 1 2 . Here, we have to take the neglected first term of eqn (18) into account again, and from this we see that g s actually vanishes. This is to be expected, since in this case the local hybrid reduces to a purely (semi-)local functional. Fig. 1 shows a numerical verification of the above analytical considerations (see Appendix A for numerical details). It depicts the xc (KLI) potential corresponding to the local hybrid of eqn (9) in comparison with the asymptotic decay according to eqn (20) and (21) for the carbon atom. An additonal curve indicates the exact À1/r decay, which is clearly not reached. The xc potential, instead of decaying with g m = 1 as one would intuitively expect, 69 approaches the predicition of eqn (21) (g m = 0.716) quite rapidly. Fig. 2 shows the xc energy density e xc (r) for the same system in comparison to its correct asymptotic of À1/(2r). Clearly the xc energy density shows the correct asymptotic, cf. eqn (8).
We thus see that while the behavior of e xc can directly be controlled via the LMF in eqn (5), the process of finding the local xc potential via functional differentiation leads to non-local evaluations of the LMF that decisively impact the potential's asymptotics.
A physically meaningful quantity closely related to the asymptotics of the xc potential is the highest occupied eigenvalue e ho . Table 1 shows Àe ho compared to the experimental IP for the carbon and fluorine atoms for different functionals, together with the corresponding value of g s of the xc potential from eqn (21).
The LSDA, as generally known, significantly underestimates the IP due to the wrong potential asymptotics and the inherent self-interaction error. Using pure EXX with the correct asymptotic decay and no self-interaction error leads to a much better prediction of the IP. When employing a local hybrid with the LMF f t (r), the explicit dependence on the parameter c becomes evident: with growing c, the asymptotic value g s grows and the description of the IP improves. Fig. 3 sheds further light on the situation. It shows potentials of local hybrids which are all based on eqn (9) but use different values of c. Growing values of c increase the amount of EXX and lead to an overall deeper potential. This explains that the eigenvalues become more negative.
We thus see that while all of the local hybrids used here can (so far, see caveat in the next section) be thought of as being one-electron self-interaction free, they show different potential asymptotics and their highest occupied eigenvalues predict the IP with significantly different reliability. The relation between freedom from self-interaction, potential asymptotics, and physical Fig. 1 The xc potential v xcm (r) for the C atom along the x-axis (see Appendix A for definition), computed using f t (r) with c = 0.5. Also displayed is the asymptotic curve according to eqn (20) with g m (c = 0.5) = 0.716 and the correct asymptotic À1/r. The inset shows the potential plotted along the complete axis. , computed using f t (r) with c = 0.5. Also displayed is the corresponding asymptotic slope of À1/(2r). The inset shows the energy density plotted along the complete axis.  interpretability of the highest occupied eigenvalue as the negative IP is therefore much less clear than intuitively believed. This observation also calls for taking a closer look at the iso-orbital indicator g(r) that is used in enforcing freedom from selfinteraction. This is the topic of the next section.

The implications of orbital nodal planes
As explained in the preceding sections, many local hybrids and other functionals such as meta-GGAs rely on the function g(r) tending to 1 to detect regions of space in which a single orbital shape dominates the density, and then, e.g., correct for selfinteraction in such regions. However, a first caveat that one has to take note of is that g(r) -1 holds for one-particle densities of ground-state character. This is a possibly far reaching restriction for the use of g(r) because electron orbital densities typically have nodes, i.e., are not of ground-state character. As a specific, illustrative example, consider an atom where the HOMO has an azimuthal quantum number m i.e., is expressed in spherical coordinates as j ho (r) = R(r,y)e imf . In the region where the density is HOMO-dominated, t W (r) is therefore given approximately by 1 2 |rR(r,y)| 2 . However, in the same region t(r) is given approximately by 1 2 |rj ho (r)| 2 = 1 2 (|rR(r,y)| 2 + m 2 R 2 (r,y)/(r sin y) 2 ). Thus, if m = 0, then t(r) = t W (r), but for m a 0 this is no longer the case. Instead, t(r) and t W only approach each other asymptotically, as the m 2 -dependent-term of t(r) decays to zero with large r.
One may counter-argue that this restriction is not so severe because in density functional construction, the condition t(r)t W (r) is mostly used to detect those regions of space in a finite system which are far from all nuclei, and where the density decays nodelessly. However, we here show that even in such regions the condition t(r)t W (r) can be violated. This leads to a second caveat about the reliability of the g(r) indicator. It is rooted in the existence of orbital densities that have nodal planes or nodal axes. Fig. 4 illustrates this case. It shows g(r) evaluated in the (xz)-plane for the carbon atom (see Appendix A for numerical details, including grid setup). The density here was obtained using f t (r)(c = 0.5), but the density features relevant here are not sensitive to functional details. The important observation is that g(r) approaches 1 in the asymptotic limit in every direction -except for in the vicinity of x = 0.
The first step towards an understanding of this finding is to note that the z-axis is a nodal axis, being the intersection of the nodal planes of the two HOMOs of carbon, which are degenerate and of p-orbital character.
The consequences of the existence of nodal planes can be studied analytically. To this end we look at a schematic density that is dominated by the HOMO j ho (r), but also take the next lower lying orbital j hoÀ1 (r) into account, i.e. n(r) B |j ho (r)| 2 + |j hoÀ1 (r)| 2 . With this ansatz one finds t W $ rjj ho j 2 þ rjj hoÀ1 j 2 2 8 jj ho j 2 þ jj hoÀ1 j 2 (22) and t B 1 2 (|rj ho | 2 + |rj hoÀ1 | 2 ).
These two terms combined and evaluated on or close to a nodal plane (denoted by À! n:p: Even though j ho vanishes on the nodal plane, its gradient still yields a finite value and keeps the function g(r) from approaching 1. Fig. 4 shows that the deviation from 1 has a noticeable spatial extension of a few a.u. This raises the question of how well the use of the iso-orbital indicator g(r) leads to freedom from self-interaction, as in some regions that so far have been considered as iso-orbital ones, e.g., all space far from the system's center, self-interaction effects may not be eliminated fully when the indicator aberrates due to the presence of a nodal plane or axis. A different interpretation of Fig. 4 would be Fig. 4 The function gðrÞ ¼ t W ðrÞ tðrÞ on the numerical grid for the C atom. Fig. 3 Asymptotics of the xc potential v xcm (r) for the C atom along the x-axis, computed with pure EXX and local hybrids using f t (r) from eqn (9) with parameters c = 0.5 and c = 5.0. Also displayed are the complete potentials in the inset.
to reconsider one's expectation of where iso-orbital regions are, or what they are. The traditional point of view has been that all space far from a finite system's center is of iso-orbital nature. Fig. 4 and eqn (24) may be interpreted to show that this is not the case when the HOMO has a nodal plane/axis extending to infinity. From this perspective one might say that g(r) does exactly what it is supposed to be doing, i.e., it indicates that the nodal plane region is not of iso-orbital character. Yet, also from this perspective Fig. 4 reveals a surprising finding, namely that even infinitely far from a finite system's center, the density may not be of iso-orbital character.
The nodal plane observation also forces us to take a yet closer look at the central topic of this perspective, the potential asymptotics. Nodal planes can influence the asymptotics of a local hybrid's xc potential in two ways. First, it has been argued that all orbital-dependent functionals show non-vanishing asymptotic constants in their xc potential along nodal planes of the highest occupied Kohn-Sham orbital that extend to infinity. This was first discussed in ref. 94 and 95 for the case of pure EXX, and the occurring shift was determined to be with % v xcis = Ð j is *(r)v xcs (r)j is (r)d 3 r and % u xcis = Ð j is *(r)u xcis (r)j is (r)d 3 r. The index M s denotes the highest lying Kohn-Sham orbital that does not show a vanishing spin-orbital density along the nodal plane of the HOMO. Since eqn (25) follows from the KLI (OEP) equation without referring to a specific functional, non-vanishing asymptotic constants on nodal planes of the HOMO are expected on rather general grounds.
Second, the fact that g(r) -1 is not guaranteed on a nodal plane can also affect the potential. For the sake of clarity, we again discuss this effect for the specific example of local hybrids. When the LMF tends to zero on the nodal plane, i.e., f ðrÞ À! n:p: 0 and eqn (7) is obeyed, then the non-vanishing constant of eqn (25) is the only effect. An example for this case is the LMF f t (r) with a finite value of the parameter c. It is depicted in Fig. 5 for the C atom density in the (xz)-plane, and one sees that there are no asymptotic features. This is because the reduced density gradient in the denominator causes f t (r) to vanish in the asymptotic limit, regardless of the occurrence of a nodal plane. The potential decays like Àg s /r in all directions, but along the z-axis a non-vanishing constant v xcs ðrÞ À! n:p: appears. This is shown in Fig. 6 for f t (r)(c = 0.5) and the F atom.
One can clearly see how v xcm (r) decays with g m = 0.6650, but, instead of zero, approaches a constant of C m = 0.0244, in agreement with eqn (25). A different situation occurs when f ðrÞ ! = n:p: 0, i.e., the behavior of the indicator function along a nodal plane/axis of the HOMO prevents the LMF from reaching its intended limit. This happens, e.g., for f 0 (r) or f t (r)(c = 0) and is depicted in Fig. 7, again for the C atom. The occurrence of a nodal axis here very clearly affects the LMF. Since in this case eqn (7) is violated in the direction of the z-axis, the previous derivations cannot be used to predict the potential's asymptotic behavior. However, we have numerically checked the xc potential's behavior. On the   nodal axis it neither tends to À1/r, nor to C s À g s jrj , but rather tends to some other value. Thus, the nodal axis in this case has a very noticeable influence on the potential asymptotics, which is hard to predict a priori.

Conclusions
With local hybrid functionals serving as an explicit example we have argued that freedom from self-interaction in the sense of eqn (1) does not necessarily lead to the expected À1/r decay of the local Kohn-Sham xc potential. We have further argued that the ratio of the von Weizsäcker kinetic energy density to the positive Kohn-Sham kinetic energy density, which is frequently used in functional construction for indicating iso-orbital regions and eliminating selfinteraction effects in these, may not serve its intended purpose because it is very sensitive to excited state features such as orbital nodal planes that are present in Kohn-Sham orbitals that construct ground-state densities of many-electron systems. These findings have immediate and somewhat discomforting consequences for the local hybrid approach. For a large class of functionals one has to accept that the correct long-range xc potential simply cannot be obtained. This observation plays a role in explaining why it is very hard to construct a local hybrid that yields good binding energetics and physically meaningful eigenvalues with the same functional form and set of parameters. 83 However, the impaired relation between self-interaction and the xc potential's asymptotics, and also the impact of nodal planes, stand in a context that is much larger than the local hybrid one. The iso-orbital indicator g(r) has been used in many functionals, not only local hybrids. Nodal planes are known to impact the exact exchange potential in surprising ways. 94,95 They have appeared here as a prominent feature in kinetic energy ratios, and we expect 96 that they play a much larger role in the exchange potential than has been realized so far. The observation that a one-electron self-interactionfree energy can go together with a potential that does not fall of like À1/r is not only a feature of local hybrids, but has also been reported for a ''scaled down'' version of the Perdew-Zunger self-interaction correction. 97 One may therefore wonder whether semi-local indicator functionals are in some sense incompatible with the fully nonlocal self-interaction correction that is achieved by EXX or full Perdew-Zunger-type correction approaches. It has also been pointed out recently 98 that eqn (1) itself, which is the basis of the present definition of one-electron self-interaction, leads to questions when evaluated for orbital densities, because E xc [n] is intended to be used with ground state densities, whereas orbital densities are excited state densities. Further conceptual questions about eqn (1) relate to its inherent identification of orbitals with electrons and its unitary variance. 5,99,100 The success of self-interaction corrections schemes that rely on eqn (1) tells us that the equation is meaningful. However, the sum of the insights into its limitations that emerged over the years suggests that there is more to the question of selfinteraction in density functional theory.
While the above considerations point out areas that require further thought and work, one should also note that there have been developments in DFT that shine a bright light into the future.
The concept of many-electron self-interaction 101,102 is not as straightforward to use as eqn (1), but it avoids the conceptual questions that are associated with this equation. Rangeseparated hybrids yield the correct asymptotic potential and have proven to be a very successful concept, without being selfinteraction-free. 23,36,[103][104][105][106][107][108][109][110][111][112][113][114][115][116][117] There have been successful functional constructions that can be seen as combinations of the local hybrid and the range-separation idea. 118,119 Ensemble corrections 120 allow to extract information from functionals in an unexpected way, and can, e.g., further improve IP prediction. Finally, it has recently been shown 121 that a new type of a generalized gradient approximation can show features that were so far thought of as being associated only with exact exchange, such as step structures and surprising nodal plane features, 96 and understanding potentials in terms of xc charges has provided new insights. [122][123][124] Therefore, the battle against DFT's old foe, the selfinteraction error, and its surprisingly independent side-kick, the wrong potential fall-off, is far from being lost.

Appendix A: numerical details
We used the all-electron code DARSEC 125 for all calculations presented in this perspective. This code exploits the rotational symmetry of diatomic molecules along the interatomic axis z, treating the azimuthal angle f analytically and thus effectively reducing the problem of solving the Kohn-Sham equations in two dimensions. The equations are represented on a real-space grid of prolate-spheroidal coordinates. In such a coordinate system, the nuclear position(s) coincide with the focal point(s) of the grid located at z = AER/2, with R being the bond length of the diatomic molecule. This is the case also for calculations of single atoms: the position of the nucleus is not equivalent to the origin of the coordinate system, but it is located at z = ÀR/2, where R was set to 0.5 a.u. e.g., the C atom in our plots is centered at r C = (x,z) = (0, À0.25). The x-axis is defined as perpendicular to the z-axis, crossing the latter at z = 0, i.e., at a point being equidistant from the focal points of the grid (see ref. 125 for details).
In order to avoid numerical instabilities due to singularities in the Laplacian, the grid was chosen such that it does not include the actual z-axis, i.e. the interatomic axis. As a consequence, in this direction all quantities can only be plotted along a projected z*-axis, which takes into account all grid points that are closest to the actual z-axis. Since the discrepancy between the projected and the real z-axis decreases with increasing number of grid points, we made sure that the difference between z and z* is small by choosing sufficiently dense and large grids. sufficient rate in the derivation of eqn (20). In the present work, we investigated two possibilities for the decay of the LMF.
First, f t (r) for a finite value of the parameter c vanishes exponentially because f t ðrÞ $ t À2 ðrÞ $ e À 2 3 ffiffiffiffiffiffiffiffiffi À2e ho p r . In this case, all individual terms in each functional derivative, u c-nl is (see eqn (17)), vanish exponentially in the asymptotic limit as well, except for the second term in eqn (18). Eventually, this remaining term is responsible for the reduced asymptotic decay of eqn (20) due to the non-local evaluation of f t (r). Consequently, the xc potential in both spin channels decays with Àg s /r. However, a different picture emerges when evaluating f 0 ðrÞ ¼ 1 À t W ðrÞ tðrÞ . This function decays much more slowly than f t (r) with finite c, as Fig. 8 shows for the carbon atom. Consequently, not all terms in the functional derivative originating from f 0 (r) vanish individually and more detailed investigations are necessary. Defining K(r) = n(r)(e sl x (r) À e ex x (r)), the functional derivative in this case reads u c-nl is ðrÞ ¼ À Therefore, both df 0 ðrÞ dtðrÞ and df 0 ðrÞ dt W ðrÞ reach the same absolute value in the asymptotic limit, but show opposite signs. Now, we have to distinguish between the spin channels: If one looks at u c-nl Ns ho s ho ðrÞ in the spin channel that has the global HOMO, i.e. n(r) B |j N sho s ho (r)| 2 , then one can see from eqn (27) that the fourth and fifth terms are equivalent in the asymptotic limit except for the sign. Therefore, they cancel each other and, since the first and second terms decay fast enough, only the third term remains, leading to the limit of Àg s ho /r. In the other spin channel however, the fourth and fifth terms do not cancel anymore, since the density is still dominated by j N sho s ho (r), whereas the fourth term features j N sho s ho (r). Therefore, in the other spin channel yet another asymptotic limit is obtained, again strictly following from the evaluation of the functional derivative.
This feature can be corrected by using a spin-polarized ansatz with an indicator function that is a spin-polarized LMF of the form g s ðrÞ ¼ t Ws ðrÞ t s ðrÞ , with t Ws (r) and t s (r) being the kinetic energy spin densities. In this case, a functional derivative that does not feature the total density n(r) follows and therefore the aforementioned effect does not occur. However, since for the spin channel s ho all derivations made are valid independently of the form of the LMF and since this spin channel features the physical meanigful quantity Àe ho , it suffices for this work to consider the more simple LMFs instead of their spin-polarized counterparts.