From Kohn-Sham to many-electron energies via step structures in the exchange-correlation potential

Accurately describing excited states within Kohn-Sham (KS) density functional theory (DFT), particularly those which induce ionization and charge transfer, remains a great challenge. Common exchange-correlation (xc) approximations are unreliable for excited states owing, in part, to the absence of a derivative discontinuity in the xc energy ($\Delta$), which relates a many-electron energy difference to the corresponding KS energy difference. We demonstrate, analytically and numerically, how the relationship between KS and many-electron energies leads to the step structures observed in the exact xc potential, in four scenarios: electron addition, molecular dissociation, excitation of a finite system, and charge transfer. We further show that steps in the potential can be obtained also with common xc approximations, as simple as the LDA, when addressed from the ensemble perspective. The article therefore highlights how capturing the relationship between KS and many-electron energies with advanced xc approximations is crucial for accurately calculating excitations, as well as the ground-state density and energy of systems which consist of distinct subsystems. (* - equal contribution)


I. INTRODUCTION
Describing many-electron excited states at an affordable computational cost remains an important goal within solid state physics, quantum chemistry and materials science [1]. In principle, this is possible within density functional theory (DFT) [2][3][4][5][6][7][8] as the groundstate density, n(r), contains all the information about the many-electron system's ground and excited states according to the first Hohenberg-Kohn (HK) theorem [9]. However, in practice such a description is extremely challenging. The excitation spectrum, the fundamental gap (the difference between the ionization potential (IP), I, and the electron affinity (EA), A) and charge-transfer energies (the difference between the IP of the donor, I d , and the EA of the acceptor, A a ) are of particular importance [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. The unreliable performance of standard exchange-correlation (xc) approximations for these quantities is in contrast to the remarkable success of Kohn-Sham (KS) DFT for various applications to ground state properties of materials [26,[29][30][31][32][33][34][35][36][37]. In this article we explore the exact relationship between KS excitation energies and the corresponding many-electron quantities with standard and ensemble DFT. We study the consequences of this relationship on the exact KS potential and its importance for the advancement of approximate xc density functionals.
Unlike other commonly used methods for electronic structure calculations, e.g., many-body perturbation theory [38][39][40], within KS-DFT the relationship between the KS energy levels, {ε i }, and the many-electron energies, {E i }, is not generally straightforward. For example, while for the exact KS potential the highest occupied (ho) * These authors contributed equally KS energy level, ε ho , equals minus the IP [41][42][43][44][45][46][47], −I, the fundamental gap, E g = I − A, does not simply equal the KS gap, E KS g = ε lu − ε ho (i.e., the difference between the lowest unoccupied (lu) and the ho KS energies), even for the exact KS potential. Instead, the KS gap differs from the fundamental gap by ∆, known as the derivative discontinuity [41,42,[48][49][50][51][52][53][54][55][56][57][58][59][60][61]: E g = I − A = ε lu − ε ho + ∆. ∆ manifests in the exact xc potential as a uniform shift when the number of electrons within the system infinitesimally surpasses an integer. It occurs because the xc energy of the system has discontinuities in its derivative as a function of electron number, N , at integer values of N .
Similarly, it has been shown recently [62] that the charge-transfer energy in stretched systems differs from the corresponding KS energy difference by the chargetransfer derivative discontinuity (CTDD), ∆ CT , which occurs when a fraction of charge is transferred from one subsystem to another within the whole system. The CTDD proved to be an important concept for accurately modeling charge transfer within KS theory in practice [63].
In 1995 Levy proposed that the optical (uncharged) gap, i.e., the energy to excite an electron from the ground to its first excited state ( ω og ), is related to the corresponding KS gap (ε lu − ε ho = ω KS og ) via a derivative discontinuity [64], as such ω og = ω KS og + ∆ og . All the discontinuities mentioned above -∆, ∆ CT and ∆ og -are important and rather delicate properties of the exact xc functional. Their existence gives rise to step structures in the exact xc potential -sudden changes in the magnitude of the potential over a short region of space. These steps have a strong nonlocal dependence on the electron density, which partly explains why they are not captured by most existing approximations.
In Ref. 62 the relationship between the derivative discontinuity ∆ in the xc energy and the spatial step S that appears in the exact xc potential of stretched diatomics was established. In this article we further study the step structure of the exact xc potential and relate it to the excitation energies of the interacting many-electron system. Particularly, we show how the steps are crucial in the prediction of the fundamental gap, excitation energies, such as charge transfer, and the correct distribution of charge in stretched systems.
This article is organized as follows. Section II gives a detailed introduction to the interatomic step S within a stretched diatomic molecule in its ground state. Then the derivative discontinuity, ∆, which occurs for groundstate systems with a fractional electron number, is discussed. Finally the CTDD, ∆ CT , is analytically studied for both a stretched diatomic molecule with a fractional N and for a stretched diatomic molecule that experiences charge transfer upon excitation. Section III provides the numerical details of the calculations performed in this work. Section IV discusses the relationship between ∆ and S, numerically addressing finite and stretched systems. Section V presents the exact KS potential obtained from an excited-state calculation of a one-dimensional (1D) stretched diatomic molecule, which undergoes charge transfer. Then, in Sec. VI an excited atom is analyzed to show that steps and plateaus in the KS potential appear not only for a stretched, but also for a finite system, upon excitation within ensemble DFT. In Sec. VII we show that steps can be found not only in the usually unreachable exact KS potential, but also in approximate potentials, as simple as the one that stems from the local density approximation (LDA), by means of numerical inversion of the LDA ensemble density. Finally, in Sec. VIII we summarize our work.

II. PROPERTIES OF THE EXACT EXCHANGE-CORRELATION POTENTIAL
A. The spatial step S In general, the sharp spatial steps may occur in the exact xc potential [51,65,66] at any point where the electron density decays at a rate which abruptly changes. One scenario is an atom with spatially distinct electron shells (see, e.g., Refs. 67 and 68). In this case, approaching the atom inwards from infinity, the decay of the outermost shell is substituted by the decay of the next, inner shell. The potential then experiences a step, which can be revealed [69][70][71], particularly when using orbital-dependent, exact-exchange-based approximations, within the optimized effective potential (OEP) method [72][73][74][75][76][77]; however, this approach has known numerical difficulties which arise from the use of a finite basis set [74,75,78,79].
Another, very important scenario is a complex system, which consists of several spatially distinct subsystems, e.g., atoms within a molecule. For such systems one can introduce the local effective ionization potential (LEIP) [68], which stems from the decay rate of a given subsystem. Moving from one subsystem to another leads to a change in the LEIP, which causes a sharp spatial step in the xc potential. The height of the step is analytically derived below from the density decay rate, following Refs. 62 and 68.
A simple and instructive example of a system with a step in the xc potential is a stretched diatomic molecule L · · · R sketched in Fig. 1. In this case each atom within the system can be considered a subsystem. Additional, more complicate examples include donor-acceptor pairs, which are important in photovoltaics [80][81][82][83][84][85] and a molecule between two metallic contacts in a transport experiment [21,28,[86][87][88][89][90]. Therefore, understanding the steps in the exact KS potential is crucial, as it allows one to accurately describe various scenarios in real materials of high practical importance with KS DFT.
In the diatomic molecule L · · · R with interatomic distance d = |d|, Atom L is located at x = − 1 2 d and Atom R at x = 1 2 d with x being the interatomic axis (see Fig. 1). In the limit d → ∞, the energy of the molecule equals the sum of the energies of the constituent atoms (the subsystems), as such and the density is the sum of the (shifted) atomic densities: see Fig. 2 (top). The equilibrium number of electrons in the molecule is thus N 0 L···R = N 0 L + N 0 R . Now we ask what form the exact KS potential of the whole molecule, v KS L···R (r), takes for large d [91], and how it relates to the atomic potentials, v KS L (r) and v KS R (r). Is it that, similarly to the molecular density, There is reason to think that the limit above holds, at least in the vicinity of each atom, because near, say, Atom L, the molecular potential v KS L···R (r) has to reproduce the atomic density n L (r + 1 2 d). From the HK theorem [9] we know that this potential is unique, up to a constant, and equals v KS L (r + 1 2 d) [92]. The same is, of course, also true for Atom R. However, the simple superposition of the atomic KS potentials can create the following problem (see Fig. 2 (middle)): the lu KS energy level of one of the atoms (say, Atom R), ε lu R , can lie below the ho Top: A sketch of the density nL···R(r) in a stretched diatomic molecule, L · · · R. Middle: The atomic potentials, v KS L (r + 1 2 d) and v KS R (r − 1 2 d), and their ho and lu KS energy levels. The problem of ε lu R lying below ε ho L is illustrated. Bottom: The molecular KS potential, v KS L···R (r) (blue), compared to the atomic potentials (gray). v KS L···R (r) possesses a step S in between the atoms (as well as a complementary step (down), −S to the right of Atom R, not shown). The molecular ho and lu KS levels are marked.
level of the other atom (Atom L), ε ho L . Then, from the perspective of the KS system, the electron which should localize on L will spuriously do so on R, resulting in the wrong number of electrons on each atom [93]. In the case in which the atoms within the molecule are bonded, the molecular ho levels of Atoms L and R ought to be aligned; this does not always happen if the two atomic potentials are simply superimposed.
What must the exact KS potential do to maintain the correct atomic densities in the vicinity of each atom whilst yielding the correct distribution of charge within the molecule? The answer is to raise the level of the potential around one of the atoms, in our case Atom R, forming a plateau, which results in a spatially abrupt step in the KS potential between the atoms (and a comple-mentary step far to the right of Atom R) [44,65]. In the vicinity of Atom R the molecular potential equals v KS R (r − 1 2 d), up to a constant, hence no violation of the HK theorem occurs. The density in this vicinity equals n R (r − 1 2 d), as required. Following Ref. 68, we now show how the height of the step in the KS potential of a stretched diatomic molecule is related, in the general case, to the IPs of the constituent atoms, I L and I R , and the molecular orbital energies of the system as a whole (see also Ref. 62 and references therein). We consider, therefore, a diatomic molecule L · · · R with a large, but finite separation d and assume that it has been solved within KS DFT, and the molecular KS potential, v KS L···R (r), as well as the molecular energy levels are known; see Fig. 2 (bottom). We denote here the molecular KS energy levels by {η i } to clearly distinguish them from the atomic KS energy levels, {ε i }. We also explicitly indicate whether the molecular orbitals localize on one of the atoms by the subscripts L and R. Generally, in the vicinity of Atom L the molecular KS potential v KS L···R (r) is identical to the atomic potential, v KS L (r + 1 2 d), up to a constant, v , and in the vicinity of R [94]. Furthermore, in the vicinity of Atom L the molecular density n L···R (r), which equals the (shifted) atomic density, n L (r + 1 2 d) (see Eq. (2)), and decays as ∼ exp(−2 √ 2I L |r + 1 2 d|) [95] [43,45,[96][97][98][99]. From the KS perspective, the decay of the atomic density is governed by the square of the ho KS orbital, which is localized on L, |ϕ ho L (r)| 2 . This orbital decays as [100] |ϕ ho As the exact KS density equals the many-electron density, the two decay rates are equal and hence v = η ho L +I L . Similar analysis for Atom R yields v = η ho R + I R . Combining these two results, and recalling that S = v − v , we arrive at an expression for the interatomic step [68]: Importantly, the constraint that the multiplicative KS potential must yield a single-particle density which exactly equals the many-electron density leads to the step S in the potential [101]. The step is generally nonzero, because the KS energy differences do not equal the manyelectron energy differences, as mentioned in the Introduction. In the particular case here, I R − I L = η ho L − η ho R . The step forms at the point in the electron density where the decay from the left meets the decay from the right, and the LEIP abruptly changes.
We wish to emphasize that the right-hand side of Eq. (4) includes the molecular energy levels, {η i }, and not the atomic levels, {ε i }. Therefore, in general, Eq. (4) does not allow one to directly obtain the step height in the molecular potential, S, relying only on atomic calculations. This equation rather shows the relationship between S, the molecular KS energies and the manyelectron energies, I L and I R , associated with each atom.
Equation (4) refers to the general case, where L and R can be any atoms, and therefore the energies η ho L and η ho R need not be assumed equal. The latter is true when L and/or R are closed-shell atoms. In the particular case that L and R are bonded, the ho KS orbital stretches over both atoms and therefore it follows that, in the notation adopted here, η ho R = η ho L . As a result, Eq. (4) reduces to the famous result S = I R − I L by Almbladh and von Barth [65] , [102].
Depending on the atoms L and R, either I L or I R is the overall IP of the molecule; in the case depicted in Fig. 2 it is I L . Thus, the overall highest occupied molecular orbital (HOMO) energy is η ho L and is equal to the atomic orbital ε ho L when v = 0. Furthermore, due to the IP theorem in DFT [41,43,46,47,103,104], which we discuss in detail below, ε ho L = −I L . It then follows that Eq. (4) reduces to S = I R + η ho R . It does not necessarily follow, however, that S vanishes. A generally nonzero S stems from the inclusion of the molecular energy, η ho R , opposed to the atomic energy, ε ho R , in Eq. (4). The atomic energy ε ho R equals −I R , whereas the molecular energy η ho R does not, as it is elevated relative to the atomic energy by the step height S: η ho R = ε ho R + S Our decomposition of this molecule into fragments is reminiscent of Partition DFT (PDFT) [105] in which the exact KS potential is separated into the KS potential for each individual subsystem plus the 'partition potential'. In the limit that the subsystems are completely separated -in our case the two atoms -the partition potential consists of the interatomic step described above [106]. The partition potential is a functional of the density of each fragment of the system [107] and hence is nonlocal in character [108]. In addition, the exact partition potential is known to contain derivative discontinuities [109]. The perspective allowed by PDFT offers an approach to developing approximations which capture these discontinuous features, yield accurate binding energies of disassociated diatomics [109][110][111] or a reliable description of charge transfer [112,113]. The partition potential has also been shown to be a chemically significant reactivity potential [114,115].

B. The uniform jump ∆
The uniform jump ∆ occurs in the KS potential when the number of electrons, N , varies continuously, and infinitesimally surpasses an integer value. A fractional number of electrons in our systems of interest may be considered as a time average of the number of electrons in an open system, namely in a system which is free to exchange electrons with its surroundings (see, e.g., Ref. 116, §14). The ground state of such a system can no longer be described by a pure quantum-mechanical state. Instead, it is a statistical mixture, or ensemble, of pure (integerelectron) states [41].
In the following we consider three types of manyelectron systems. First, in this section, we describe in detail a finite system that is connected to an electron reservoir, which allows N to change continuously. Second, in Sec. II C we consider a stretched diatomic molecule L · · · R, whose total number of electrons can vary continuously, and for which any additional charge localizes on Atom R, whereas any charge deficiency results in decrease of charge around Atom L. Third, in Sec. II C we consider a stretched diatomic molecule L · · · R, whose total number of electrons is fixed at a given integer value, but the number of electrons on each atom can become fractional by transferring charge between the atoms.
We start with a finite system, like an atom or a molecule, with N = N 0 + α electrons, where N 0 is an integer number, and 0 α 1. As mentioned above, the ground state of such a system is an ensemble, which combines states each with a different integer number of electrons. For systems with Coulomb interaction at zero temperature, this ensemble consists only of states for N 0 and N 0 + 1 electrons, |Ψ N0 and |Ψ N0+1 : with the statistical weights of (1 − α) and α, respectively [2,41,[117][118][119]. As a direct consequence of Eq. (5), the expectation value of any operatorÔ in the ensemble state is [41]. In particular, the average density of a system with N electrons is n(r; N ) = (1 − α)n(r; N 0 ) + αn(r; N 0 + 1), where n(r; N 0 ) is the ground-state density for the N 0electron system and n(r; N 0 + 1) is the ground-state density for the (N 0 + 1)-electron system. Furthermore, the total energy as a function of N equals As can be seen in Fig. 3(top), E(N ) is piecewise-linear in N : for any fractional N , the energy is linear, but it can change its slope when N passes an integer. Consequently, the chemical potential, µ = ∂E/∂N , is a stair-step function of N . For example, in the ground state: where is the EA of the system. Clearly, the chemical potential is generally discontinuous at integer N ; the height of this discontinuity equals the fundamental gap of the system, E g = I − A. Furthermore, from combination of piecewise-linearity of the energy and Janak's theorem [120], which states that the i th KS eigenenergy, ε i = ∂E/∂f i -the derivative of the total energy with respect to the occupation of the i th level, f i -we find that the ho KS energy level, ε ho (N ), equals the chemical potential, µ(N ), and is also discontinuous at integer N (see Fig. 3(middle)). This is the content of the IP theorem in DFT [41,43,46,47,103,104]: for the exact xc potential, infinitesimally below an integer, ε ho (N − 0 ) = −I and infinitesimally above ε ho (N + 0 ) = −A. The IP theorem in KS DFT is an exact result, for the exact xc potential.
Satisfying the aforementioned IP theorem creates a challenge for the exact xc potential, v xc (r). From the perspective of the KS system, increasing N above an integer means occupying the next KS level, ε lu (N − 0 ). As ε lu (N − 0 ) does not necessarily equal −A, even for the exact KS potential (see Fig. 3(middle)), the only thing the exact potential can do in order to satisfy the IP theorem is to discontinuously change as N infinitesimally surpasses an integer. However, due to the continuity of the density with N (see Eq. (6)) and the HK theorem, the discontinuity of the KS potential can change only by a spatially uniform constant (see Fig. 3(bottom)), which is usually denoted ∆. This discontinuity in the KS potential, v KS (r), can only come from v xc (r), because the Hartree potential is continuous and the external potential is N -independent. Therefore, The value of ∆ is easy to deduce from the arguments above: it is the difference between the value that ε ho (N + 0 ) ought to have, namely −A, and the value it has in absence of discontinuity, . Together with ε ho (N − 0 ) + I = 0, and dropping here the argument N − 0 for brevity, we arrive at the following familiar form for ∆: where ∆ is expressed as the difference between the fundamental gap of the system, E g = I −A, and the KS gap, The derivative discontinuity is a topic of great importance and has received much attention over the years [41, 42, 48-54, 56-60, 104, 121, 122]. Yet, many common approximate xc functionals lack this important feature; advanced approximations are being developed to reconstruct it (see, e.g., [11, 12, 24, 27, 28, 36, 55, 57-60, 104, 123-146].

C. Charge-transfer derivative discontinuity
Let us now consider a stretched diatomic molecule L · · · R, where the separation between the atoms is large enough for the energy and density of the molecule to satisfy Eqs. (1) and (2). At first, the molecule possesses N 0 L electrons on Atom L and N 0 R electrons on Atom R, so the total number of electrons equals N 0 L···R = N 0 L +N 0 R . Next, we allow the total number of electrons to vary continuously: We consider the specific case for which any additional charge localizes on Atom R, whereas any charge deficiency results in decrease of charge around Atom L. This case is indeed specific but not esoteric -it is the prototype case for a donor-acceptor pair.
Combining Eqs. (1) and (7) we can conclude that the total energy of the molecule is piecewise-linear with the number of electrons (see Fig. 4): Dependence of the total energy of a stretched diatomic molecule L · · · R on α -the deviation of the total number of electrons from its integer value, N 0 L···R . The slopes of the graph are associated with the IP of the left atom and the EA of the right atom.
The chemical potential of the molecule as a whole, being the derivative of its energy with respect to N L···R , or equivalently to α, is a stair-step function discontinuous at integers, qualitatively similar to the chemical potential depicted in Fig. 3(middle): Notably, here the height of the discontinuity in µ L···R is the left-to-right charge-transfer energy, E CT L→R = I L −A R , namely the energy required to remove one electron from Atom L minus the energy gained by adding an electron to an infinitely distant Atom R. As for the finite system discussed above, the stretched molecule L · · · R also obeys the IP theorem. Namely, the overall HOMO energy, η ho (N L···R ), has to equal µ L···R (N L···R ). For N L···R slightly below N 0 L···R the overall ho energy equals η ho L (N 0− L···R ), which in our case, as explained in Sec. II A, equals −I L . As the overall number of electrons increases above N 0 L···R , the overall ho level is localized around Atom R and has to equal −A R . As a result, the molecular potential v KS L···R (r) jumps by the constant (cf. Eq. (10)). This quantity was first introduced in Ref. 62, where it has been termed charge-transfer derivative discontinuity. ∆ CT L→R is the difference between the charge-transfer energy, E CT L→R = I L − A R and the corresponding quantity in the KS system (η lu R − η ho L ) (cf. Eq. (10)).
Finally, we consider a stretched but finite diatomic molecule in which the atomic separation is large enough to define individual atoms within the molecule but in which the electrons localized on the left atom experience the Coulomb repulsion of the electron localized on the right atom and vice versa. The total number of electrons within the molecule, N 0 L···R , is constant and inte- Dependence of the total energy of a stretched diatomic molecule L · · · R on q -a fraction of an electron transferred from Atom L to Atom R. The slopes of the graph are associated with the IPs and the EAs of the constituent atoms.
ger. When the molecule is excited a fraction of q electrons is transferred from Atom L to Atom R. We define an ensemble consisting of the ground state, |Ψ 0 , of the molecule and the first excited state |Ψ 1 where the latter has charge-transfer character, i.e. the nature of |Ψ 1 is such that compared to the ground state, one electron is transferred from Atom L to Atom R. The statistical operator describing this ensemble is give bŷ Both states, |Ψ 0 and |Ψ 1 have fixed (integer) particle number N 0 . The ensemble expectation value of any operatorÔ, by virtue of Eq. (14), In particular, the ensemble density is given by where n 0 (r) and n 1 (r) are the densities of the ground state and the first excited state, respectively. Likewise the total ensemble energy as a function of q equals where the subscript 0 corresponds to the ground state, whereas the subscript 1 corresponds to the first excited state. Therefore, E 1 = E 0 + E CT ; for this system with a large but finite atomic separation, E CT =Ĩ L −Ã R for q > 0 and E CT =Ĩ R −Ã L for q < 0, whereĨ L is the ionization energy of the whole molecule which corresponds to an electron localized to the left atom whileÃ R is the molecule's affinity and corresponds to the addition of an electron to the right atom once the electron on the left atom has been ionized -this is the nature of a chargetransfer excitation. Consequently, bothĨ L andÃ R are influenced by the Coulomb interaction between the left and right atoms; this effect has previously been emitted because the atoms were assumed to be infinitely separated. lim d→∞ĨL −Ã R = I L − A R (as defined above). By modeling the system with a finite separation we more closely model a real donor-acceptor pair for short-to mediumrange charge transfer. The difference betweenĨ L −Ã R and I L − A R is the electron-hole electrostatic interaction. For large separation between the donor and acceptor, it is usually approximated as −1/d [41,147,148]. Plugging this definition for E CT in this system into Eq. (16), we obtain Analogously, for a charge transfer from R to L, we obtain Hence the total energy is piecewise-linear with respect to q (see Fig. 5). Therefore, its derivative, m(q) = ∂E L···R /∂q, which is the change in energy as a result of transfer of charge, is a stair-step function: From the Gross-Oliveira-Kohn (GOK) theorem [149][150][151], we can express the charge-transfer energy as such where η q i is the i th KS energy of the ensemble system.
Therefore, recalling that in the limit of infinite atomic separation Eq. (20) is equivalent to Eq. (13), we arrive at an expression for the CTDD for the ensemble system, defined in terms of the derivative of the ensemble xc energy: This expression allows one to calculate the CTDD from any explicit q-dependent xc functional [152][153][154]. In Ref. 155 ∆ CT L→R -as it is defined by Eq. (21) -was evaluated experimentally for donor-acceptor pairs.
Note that in the limit that Atom L and Atom R become infinitely separated, m(q) equals the difference between the chemical potentials of the constituent atoms with the atomic chemical potentials given by Eq. (8).
The discontinuity in m(q) around 0, denoted here D = lim q→0 + m(q) − m(−q), equals being the sum of the left-to-right and the right-to-left charge-transfer energies. It can also be expressed as the sum of the atomic fundamental gaps: D = E g,L + E g,R .
Using Eq. (13), D can be also expressed in terms of the KS quantities: (24) in direct analogy with results presented above. D may also be expressed solely in terms of the KS gaps and ∆'s of the constituent atoms using Eq. (10): Hence, for this stretched system the derivative discontinuity, D, can equally be expressed in the KS system in terms of the derivative discontinuities of the individual atoms and also in terms of the charge-transfer derivative discontinuities of the system as a whole. We shall see below in Sec. V that the interatomic step, S, derived in Sec. II A is related to both the derivative discontinuity of the individual atoms and to the CTDDs. Finally, we emphasize two additional results. From Eqs. (24) and (25) we arrive at the following relation for the CTDDs, which shows the close relationship between them to the atomic ∆'s. Furthermore, we wish to draw attention to the following relation, which emerges from Eq. (23): meaning that the sum of the left-to-right and the rightto-left charge-transfer energies, between any two distant subsystems, equals the sum of the fundamental gaps of these subsystems.

III. NUMERICAL DETAILS
We use a 1D model to investigate the structure of the exact KS potential. Our 1D models -in Secs. IV A, V and VI -employ the iDEA code [156] in which the exact, fully-correlated many-electron wavefunction may be calculated for an arbitrary external potential. In addition to the ground state, the many-electron excited states are calculated by solving the many-electron Schrödinger equation [157]. As a result we have access to the exact many-electron ground-state and excited-state electron densities, from which the exact corresponding KS xc potentials can be calculated by a numerical inversion of the KS equations. Our inversion algorithm employs that of Ref. 156 to calculate the KS potential, with parameters p = 0.05 and µ = 0.1. For the 1D systems the KS potential is considered converged when the mean absolute error between the many-electron and KS densities is < 10 −9 Bohr −3 .
Results for Sec. VII were obtained using the ORCHID program [158], version 3.1, on a natural logarithmic radial grid, r ∈ [e c /Z, L], with c = −13 and L = 35 Bohr.
The total energy and the eigenvalues are converged below 10 −6 Hartree. The inversion procedure [156] used the parameters p = 0.1 and µ = 0.72. The convergence criterion for the inversion procedure is ln (n(r)/n target (r)) < 10 −4 , enforced for r ∈ [e c /Z, L ], with L = 30 Bohr. Finally, the parameters a and b required for the alignment of the KS potentials, which show the asymptotic behavior of ∼ a/r + b (see details in Sec. VII), have been obtained by a linear fit of the potential vs. 1/r at 20 and 30 Bohr.

IV. THE RELATIONSHIP BETWEEN S AND ∆
The properties S and ∆ of the exact xc potential discussed in Secs. II A and II B, respectively, have been known for a long time [41,44,49,65,67,[159][160][161], but whether these two are completely independent or related properties, remained elusive until recently [62]. Indeed, S and ∆ are not one and the same: first, they can be derived from two different perspectives, as performed in Sec. II. Second, the EA and the lu energy, which contribute to ∆ (Eq. (10)), are absent from the expression for S (Eq. (4)). Finally, the shift ∆ occurs when varying the charge of the system, whereas S occurs at a fixed, integer number of electrons. However, it was realized early on that both S and ∆ occur for a finite system when the decay rate of the electron density abruptly changes [44,65]. This suggests a close relationship between the two properties. In the following we characterize this relationship in detail, by formulating and subsequently resolving two paradoxes that arise from the combination of the concepts presented in Secs. II A and II B.

A. Uniform jump paradox
Paradox 1 -The spatial uniformity of the jump in the KS potential implies ∆ = 0.
In Sec. II B we described a finite system with a varying number of electrons N and concluded that as N passes an integer the KS potential jumps by a spatially uniform constant ∆. Here we address a finite system again, like in Sec. II B, but now we are applying the approach from Sec. II A. In other words, we find ∆ by examining the exponential decay of the density.
If the number of electrons in the system equals an integer N 0 or a little bit less, the density decay is determined by the IP of the system, i.e., n(r; N 0 ) ∝ exp −2 √ 2I|r| (denoted I-decay). From the KS perspective, the density decay is governed by the ho orbital squared, |ϕ ho (r)| 2 ∝ exp −2 −2ε ho (N − 0 )|r| . As the exact KS density equals the many-electron density, ε ho (N − 0 ) = −I. If the number of electrons is now slightly increased above N 0 by a small fraction of an electron, α, the density becomes a linear combination of n(r; N 0 ) and n(r; N 0 + 1), as in Eq. (6). The term n(r; N 0 +1) decays ∝ exp −2 √ 2A|r| (A-decay) which is slower than the decay of n(r; N 0 ) because I > A for all known systems (known as the convexity conjecture [2,41,53,117]). Therefore the A-decay asymptotically dominates the density decay. From the KS perspective, the decay of the density is dominated by the now highest, partially occupied orbital (the former lu orbital). The problem arises when taking Fig. 3(bottom) at face value, namely assuming that the KS potential indeed jumps by a completely uniform constant ∆. Then, one may think that the decay of the highest, partially occupied orbital is ∝ exp −2 −2(ε ho (N + 0 ) − ∆)|r| , i.e., the decay rate is governed by the ho energy, ε ho (N + 0 ), relative to the overall potential shift, ∆ (cf. Eq. (3)). Recalling that ε ho (N + 0 ) = ε lu (N − 0 ) + ∆, one may further infer that the density decays ∝ exp −2 −2ε lu (N − 0 )|r| . This leads to the paradoxical conclusion that ε lu (N − 0 ) = −A and hence ∆ = 0. In other words, if the jump ∆ is uniform, its height is zero.
To resolve this paradox we look more closely at Eq. (6), keeping in mind that in our case α → 0 + . Although n(r; N 0 + 1) decays slower and is thus the asymptotically dominant term, it is multiplied by the small coefficient, α. As a result, we have a competition between the two decay rates: when we reduce α to 0 while looking at a fixed and large r, the region in which the A-decay is dominant moves away from the nucleus as the term αn(r; N 0 + 1) vanishes and the term (1 − α)n(r; N 0 ) prevails. The process is illustrated in Fig. 6 for an exactly solved 1D model of an atom with v ext (x) = −2.0/(0.4 · |x| + 1), with 1 + α same-spin electrons. It is useful to look at the natural logarithm of the density to clearly see the decay rates, as such a region of an exponential decay appears as a linear line of negative slope. Indeed, in Fig. 6(a) we clearly observe the I-and A-regions of exponential decay. As α decreases, the A-decay region appears further away from the nucleus. Next, recalling our conclusion from Sec. II A that a change in the decay rate of the density (no matter what the reason) leads to a step in the KS potential, we indeed find in Fig. 6(b) that for all positive α the KS potential is elevated near the origin, comparing to the (α = 0)-case, and presents steps far from the origin, at the point where the decay rate changes and hence where the LEIP changes. In Fig. 6(c), subtracting the (α = 0)-potential from all the potentials of Fig. 6(b), we clearly see a plateau around the origin, in agreement with previous studies (see, e.g., Refs. [44,51,62,65,162]). As α vanishes, the width of the plateau increases, approaching infinity. However, at any finite α the plateau width is finite and asymptotically the KS potential approaches the value of 0 (and not ∆), i.e., the shift for finite α is not uniform. This resolves our paradox: the correct decay rate of the density in the region of A-decay is ∝ exp −2 −2ε ho (N + clusion that ε ho (N + 0 ) = ε lu (N − 0 ) + ∆ = −A, as required; whereas in the region of I-decay the potential is elevated by ∆. As a result, steps form in the potential as shown in Fig. 6. Thus in this case, for a finite system with varying N , the quantities ∆ and S have the following relationship: lim α→0 + S = ∆. For the system presented in Fig. 6 this has been numerically verified as ∆ was obtained also from total-energy differences.
Finally, we wish to add several comments on plateaus in finite systems. First, the shape of the steps observed includes characteristic dips clearly seen in Fig. 6(c) (cf. Refs. 62,[162][163][164][165]. These features are numerically robust and are required to yield the exact KS density, although their significance in our context is not as high as that of the step. Second, the value of the KS potential of a finite system far from its center is an example for an order-of-limits problem, namely lim |r|→∞ lim α→0 + v KS (r, N 0 + α) = ∆, whereas lim α→0 + lim |r|→∞ v KS (r, N 0 +α) = 0. In words, if we examine the value of the KS potential at some finite point |r| while continuously decreasing α to zero, for a certain α the plateau will be wide enough to reach |r| and elevate the potential there. Taking then |r| to infinity will result with the height ∆ for the KS potential. Conversely, taking |r| to infinity first while keeping α finite, ensures that for any finite α, no matter how small, we will reach the edge of the plateau and the potential value will drop to 0.

B. Charge transfer paradox
Paradox 2 -The transfer of charge in a diatomic molecule results in a plateau, ∆, around the acceptor atom. Yet, the overall interatomic step height must remain S.
To further explore the relationship between ∆ and S we study the stretched diatomic molecule presented in Sec. II A, but now taking into account also the results of Sec. II B. We consider two scenarios that model charge transfer (cf. Sec. II C): (i) The overall number of electrons in the stretched molecule is increased; the additional charge localizes on one of the atoms, say, Atom R. (ii) When we increase the number of electrons on Atom R, we decrease the number of electrons on Atom L by means of charge-transfer excitation of the molecule so that the overall number of electrons is constant. From the results shown in Fig. 6(c), we would expect a plateau of height ∆ R to emerge around the acceptor atom, in our case Atom R (with no significant change around L). But this is contrary to the results of Sec. II A: there exists a plateau of height S around Atom R, irrespective of any infinitesimal transfer of charge, to ensure the correct distribution of charge in the ground-state KS system. As S = ∆ R , and (thinking of the complimentary scenario of right-to-left charge transfer) S = ∆ L either, there appears to be a contradiction.
To resolve this paradox, we refer again to the density of the system. For both Cases (i) and (ii) the natural logarithm of the density in between the two atoms is sketched in Fig. 7(a). We expect three regions of exponential decay between the atoms: going from right to left, the density decay is first governed by I R and then by A R (changing at point (2); cf. Fig. 6(a)), due to the extra charge on Atom R. Then, the A R -decay meets the I L -decay at point (1), simply due to the fact that the two atoms form one molecule. As a result, we expect not one, but two steps in the KS potential between the atoms in this diatomic molecule (Fig. 7(b)). The height of the steps can be deduced analytically [62], similarly to the derivation of Eq. (4): the step S (2) , which depends solely on quantities related to Atom R, equals ∆ R , whereas the step S (1) equals −∆ CT L→R . Importantly, the steps S (1) and S (2) combine to yield the overall step S of Eq. (4). This resolves the paradox raised above: indeed, a plateau of height ∆ R is expected to form on the receiving Atom R FIG. 7. (a) A diagram of ln (n) far from, and in between, the atoms of a molecule L · · · R. Three regions of density decay are present: IR-, AR-and IL-regions. Transition from the IRto the AR-region occurs at point (2) and from the AR-to the IL-region at point (1). The changes in the density give rise to two steps in the KS potential (b). upon charge transfer or addition. But in conjunction, in the region of Atom L, the KS potential shifts when the 'local electron number' decreases below an integer. The combination of these two plateaus yields an overall interatomic step of height S.
The internal structure of the step S in Case (i) has been illustrated and extensively discussed in Ref. 62. The two steps, S (1) and S (2) , have been identified both in a 1D model of a stretched diatomic molecule and in a 3D (Li · · · Be) 3+ ion. Case (ii) is numerically illustrated in Sec. V below for a charge transfer in a stretched 1D diatomic molecule induced by exciting the system.

V. CHARGE TRANSFER IN A DIATOMIC MOLECULE
Simulation of a charge transfer process, and particularly obtaining the exact KS potential that describes the process is by no means a trivial task [166]. To this end it is necessary to exactly obtain not only the ground state of the system, but also its first excited state that corresponds to a charge transfer.
In this section we present a prototypical 1D stretched diatomic molecule L · · · R, which we excite in order to transfer charge from Atom L to Atom R. Our system consists of an integer number of same-spin electrons, in this case N 0 L···R = 2. Figure 8 illustrates the charge-transfer process: the external potential, v ext (x) = −4/(0.6 · |x − 7.5| + 1) − 2/(0.4 · |x + 7.5| + 1) is asymmetric, chosen such that the ground-state electron density corresponds to a system with one electron localized on Atom L and one electron on Atom R, whereas in the first excited state both electrons are localized on Atom R. Hence, by exciting this system we can initiate a transfer of charge from L to R. We first find the exact many-electron groundstate density n 0 (x) and the first excited-state density n 1 (x). Then we construct an ensemble electron density, which corresponds to a transfer of a fraction of q electrons   15) where 0 q 1 [166]. We emphasize that all the densities present in Eq. (15) integrate to an integer number of electrons.
The GOK theorem ensures a one-to-one mapping between the density and the local potential for this excited system, provided that 0 q 0.5. Hence there exists a KS system, which exactly reproduces the electron density of Eq. (15), and thus we can obtain this KS potential from the density n(x; q) by numerical inversion (Sec. III). In our case, where N 0 L···R = 2, the density is given in terms of the KS orbitals by [167] n(x; q) = |φ 0 (x)| 2 + (1 − q) · |φ 1 (x)| 2 + q · |φ 2 (x)| 2 . When the system is excited, a fraction of the electron (q) initially occupying the first excited KS orbital is transferred into the second excited KS orbital localized in our case on Atom R, while the overall number of electrons stays constant and integer; in this sense this type of excitation is uncharged (the number of electrons within the overall system is unchanged) but in the vicinity of each atom, this excitation corresponds to a charged excitation (the number of electrons changes locally). This observation may explain why approximate KS theories, such as linear response time-dependent DFT (TDDFT), struggle to accurately describe charge transfer [21,168,169]. Figure 9(a) shows the natural logarithm of the exact ground-state electron density, ln (n 0 (x)), for our diatomic molecule: each electron occupies its own potential well, and far from the well the density decays exponentially. There are two regions of decay between the atoms -the I L -and the I R -decay -and hence one step at the point where the decay of the density changes yielding a change in the LEIP; see Sec. II A. The height of this step is given by Eq. (4). Figure 9(b) shows the KS potential corresponding to this ground-state density. The potential has an interatomic step which acts to localize one electron on each atom in the KS system, as required. Another step of height −S is expected far to the right of Atom R, when the I L -decay will prevail over the I R -decay (not shown on the figure). Both steps together form a plateau of height S around Atom R. Figure 9(c) shows ln (n(x; q)), the natural logarithm of the exact excited many-electron density, given by Eq. (15), with q = 5 × 10 −4 . For reference, ln (n 0 (x)) is also shown. There are now three regions of exponential decay in the density n(x; q): the I L -, A R -and the I R -decay, as we expected (cf. Sec. IV B). These three regions of decay give rise to two steps in the corresponding exact KS potential, at the points in the density where the decay rate changes. Figure 9(d) shows the corresponding exact KS potential of our excited system with the two steps apparent, S (1) and S (2) (arrows). The right (acceptor) atom experiences the jump in the KS potential characteristic of the derivative discontinuity, i.e., S (2) = ∆ R , owing to the local number of electrons of Atom R surpassing an integer by a small amount (q). Simultaneously a plateau forms in the vicinity of the left donor atom. The height of the plateau is ∆ CT L→R , i.e., the CTDD associated with transferring an electron from left to right atom (Eq. (13)). S (1) is therefore equal to −∆ CT L→R (the minus sign describes the fact that S (1) is a step down between the atoms; whereas ∆ R is a step up). The sum of the two steps equals the overall step of Eq. (4). Figure 9 is notably similar to Fig. 2 in Ref. 62, where the same 1D diatomic molecule is modeled but for a system with a fractional number of electrons N L···R = 2.0005 in the ground state. This means that the approach chosen in Ref. 62 to reveal the internal structure of the interatomic step S and find the CTDD, employing calculation which are much cheaper numerically, is appropriate. Therefore, there is reason to assume that modeling of full charge transfer for 3D systems, as the one analyzed in Ref. 62 and others mentioned in Sec. IV B, will also yield extremely similar results to those already obtained by varying the total number of electrons.
To summarize, simulation of a charge transfer by means of excitation of a 1D diatomic molecule showed that the interatomic step hence it has an internal structure, as expected: it consists of the ∆ of the acceptor atom, in our case Atom R, and the (negative of the) relevant CTDD, ∆ CT L→R . If charge is transferred from right to left a similar picture is expected: the overall step will split as S = −∆ L +∆ CT R→L (cf. Eq. (26)). Therefore, also in the case of a stretched diatomic molecule the relationship between the interatomic step S and the ∆'s of the constituent atoms is established (Eq. (28)) via the CTDD (Eq. (13)).

VI. DISCONTINUITIES IN EXCITED FINITE SYSTEMS WITH INTEGER ELECTRON NUMBER
In Sec. V we demonstrated that derivative discontinuities arise upon excitation of a stretched system, which in-  (2)). (d) The exact KS potential corresponding to the excited density (solid blue). Two plateaus are present: one corresponding to the derivative discontinuity of the right atom, ∆R, the other corresponds to the CTDD, ∆ CT L→R . These steps combine to give an overall step whose height is given by Eq. (4). duces charge transfer. But what happens to a finite (and not stretched) system, upon excitation from its ground to first excited state, not necessarily related to a transfer of charge? To explore this question we model a single atom with an integer N in its ground and excited states, to find whether its KS potential forms any plateaus upon excitation. This concept was first proposed by Levy [64]. Below we analyze Levy's concept numerically and study how the change in the exact KS potential of the excited ensemble state varies with the ensemble weight, β.
We model a single atom in 1D with the external potential v ext (x) = −2.0/(0.4 |x| + 1) with N 0 = 2 (again, same-spin electrons). We calculate the exact groundstate and the first excited-state density. We then find the ensemble electron density employing the 1D version of Eq. (15), where q = β in this case, for β = 10 −4 , 10 −3 and 10 −2 , and invert the KS equations to find the corresponding exact KS potential associated with each den-sity.  The excited density has two regions of decay in each case: closer to the origin the I-decay region is present, but then the decay rate changes and the density decays slower. The rate of decay of this excited density is determined by I − ω 01 , where ω 01 is the energy required to excite the many-electron system from the ground to the first excited state. Due to this change in the density decay rate we expect steps in the potential of the corresponding KS system.
The steps are clearly seen in Fig. 10(b), which shows ∆v xc -the difference between the xc potential of the excited system and the ground-state xc potential. In the central region of the system the excited KS potential is elevated by a plateau of height ∆ 01 . The height of the plateau can be analytically deduced, as before: where ε β i is the i th KS energy of the ensemble system. As β → 0 + , ε β N0+1 − ε β N0 = ε lu − ε ho = ω KS og , is the energy required to excite a KS electron from the ho to the lu KS orbital. ω 01 = ω og , the many-electron optical gap. Thus, ∆ 01 = ∆ og and This equation is the same result found by Levy [64]. It is a time-independent way of calculating exact excitation energies [152], similar to the calculation of the fundamental gap (discussed above) [170].
We find that ∆ og is always relatively small, below 0.03 Hartree (< 1 eV), for different 1D atoms with a slightly less or more confining external potential, e.g., v ext (x) = −8.0/(|x| + 1), with N = 2 -this implies that exciting one KS electron for this system is indeed a good model for the many-electron excitation of the twoelectron system. Hence, for this system, ω KS og ≈ ω og which implies that as long as ∆ og is small, the groundstate KS energy levels are reasonably good approximations to the many-electron excitation energies in their own right, i.e., neglecting the contribution of the Hartreexc (Hxc) kernel within TDDFT, which has been observed by others [54,157,171,172]. For more strongly correlated systems, or indeed the charge-transfer system above, this is not the case, and the role of the Hxc kernel or the corresponding ∆ becomes crucial [173][174][175].
From the analysis in the sections above, we conclude that any electron donor experiences a discontinuous shift in its xc potential despite the local number of electrons decreasing below an integer. This discontinuity emerges because a truly isolated system with a fractional number of electrons cannot exist in reality; there must be a source of electrons, e.g., an electron reservoir (the donor), with which a finite system, like an atom or molecule, can exchange electrons (the acceptor). Imagine that the chemical potential of the reservoir is adjusted such that an infinitesimal amount of charge is transferred to the finite system. The xc potential of the system as a whole (reservoir plus the finite system) experiences a uniform shift of height ∆ CT , which is the CTDD associated with transferring an electron from the reservoir to the finite system; see Sec. II C. This shift in the potential is truly uniform as it manifests as a result of an excitation experienced by the whole system, like the atom in this section As the amount of charge transferred from the reservoir is steadily increased, a plateau localizes in the vicinity of the acceptor which is associated with the derivative discontinuity of that finite system, ∆. In conjunction, the shift in the xc potential associated with the CTDD localizes to the donor. This occurs for the diatomic molecule of Fig. 9; in this case the donor atom acts as the electron reservoir. The charge-transfer derivative discontinuity, ∆ CT L→R , manifests as a uniform shift in the xc potential of the donor-acceptor when the transferred (excited) charge is infinitesimal. As the amount of charge is increased a plateau of height ∆ R localizes to the acceptor atom which in the vicinity of just the acceptor looks to be uniform -S (2) = ∆ R in Fig. 9. In conjunction, a complementary plateau forms around the donor atom of height ∆ CT L→R because the donor and acceptor form one system -S (1) = ∆ CT L→R . Consequently, the shift to the xc potential associated with the derivative discontinuity of the finite system when the local number of electron increases above an integer, ∆, can never be truly uniform.

VII. PLATEAUS IN APPROXIMATE XC POTENTIALS
So far we have addressed exact many-electron densities and the corresponding exact KS potentials obtained from the densities by means of numerical inversion. But what happens when working within one of the common approximations to the xc functional, like the local density approximation (LDA) or a generalized gradient approximation (GGA)? Do the resultant KS potentials possess any steps or form any plateaus in the various scenarios discussed above?
The immediate answer to this question is negative. It is well-known that if one addresses a finite system with a varying number of electrons, N = N 0 + α, with, e.g., the LDA in its standard implementation (i.e. constructing the density for fractional N by occupying the last KS level with α electrons), one obtains a gradually changing xc potential, without any plateau of the sort presented in Fig. 6(c).
However, in the spirit of the present work, it is possible to obtain the KS potential for fractional N , relying on LDA densities, also in a different way: First, one solves the system self-consistently for N 0 and separately for N 0 + 1 electrons, within a given xc approximation. Second, one creates the ensemble density, n(r; N ), using Eq. (6), thus assuring piecewise-linearity of the density. Third, one obtains the KS potential, up to a constant, via numerical inversion of the ensemble density.
We obtained this 'inverted LDA' (invLDA) potential for the Li ion (N 0 = 2) for varying α. Remarkably, the potentials show a clear asymptotic behavior of ∼ a/r + b far from the nucleus (with a, b being α-dependent parameters), rather than the exponential decay of the standard LDA. This allows us to align each potential such that it decays to 0 (and not to some finite constant, b) and subsequently subtract from it the KS potential for N = 2. The resultant differences are shown in Fig. 11. We can clearly see that the invLDA KS potential does form a plateau of height S LDA = 0.134 Hartree in the vicinity of the nucleus. As α → 0 + , the height of the plateau converges and its width logarithmically approaches infinity (cf. [62,176]) A qualitative understanding of the emergence of plateaus in the invLDA can be gained by looking at the density decay rates, presented in Fig. 12. Surely, the decay rate of the ensemble densities obtained via Eq. (6) is slower than the decay rate of the density obtained from a standard LDA calculation with fractional occupations. Then, clearly, whereas the change in the decay rates of the latter yields a plateau of height zero, a density with a slower decay will yield a non-zero plateau.
Next, we establish the quantitative relationship between S LDA found with invLDA and ∆ LDA for Li + ob-  Hence, according to Eq. (10), ∆ LDA = 0.7288 Hartree. Alternatively, ∆ LDA = −A LDA − ε lu LDA = 0.0475 Hartree. For the exact xc functional, ∆ = ∆, but for an approximate one, like the LDA, the above equality is not necessarily true, because the IP theorem is not obeyed. In any case, neither ∆ LDA nor ∆ LDA seem to equal S LDA .
We resolve the above conundrum by realignment of the KS potentials to satisfy the IP theorem [177]. This means that for each α the KS potential is shifted by the amount required for the ho level to equal the IP. For 0868 Hartree is required. We denote ∆v = v 1 − v 0 = −0.5945 Hartree and recall that ε ho (N + 0 ) = ε lu + lim α→0 + S, to find that For the exact potential ∆v = 0, and we return to the basic relationship between ∆ and S derived in Sec. IV B. To summarize, within approximate KS DFT calculations for finite systems with a fractional N there are two, equally legitimate approaches to obtain the KS potential. They lead to two qualitatively different results: the standard approach yields a smoothly varying potential, without steps, which exponentially decays at infinity. The invLDA approach yields steps in the KS potential, and the asymptotic decay is ∼ a/r. We relate these improvements to the piecewise-linearity in the density, which is enforced in the invLDA approach. This internal inconsistency within semi-local xc approximations closely relates to another inconsistency: the IP of finite systems, like atoms and small molecules can be obtained with common xc approximations from total-energy differences with high accuracy of a few percent, whereas obtaining the same quantity directly from the ho energy level results in discrepancies of ∼ 50% (see, e.g., Refs. 135 and 158 and references therein); whereas when the associated ∆ is added to the KS energy difference the exact manyelectron energy difference is obtained for the exact xc potential (as shown above).

VIII. CONCLUSIONS
In this article we studied the relationship between the Kohn-Sham energies and the many-electron energies of various systems, such as atoms and diatomic molecules, and related them to the step structures that appear in the exchange-correlation (xc) potential.
Steps can occur in the exact potential in different scenarios. In this article we address four: (i) a finite system (an atom) in the ground state with a varying number of electrons (Sec. II B and IV A); (ii) a finite, excited system with a constant number of electrons (Sec. VI); (iii) a system comprised of subsystems (stretched diatomic molecule) in the ground state with a varying overall number of electrons (Sec. II A and V); (iv) a system comprised of subsystems that experiences a charge transfer upon excitation (Sec. II C and V). By these examples we address the processes of ionisation, excitation, dissociation and charge transfer.
As a general rule, steps in the potential occur at points where the exponential decay rate of the density changes, and hence changes the 'local effective ionization potential' (LEIP) [68]. This rule is true irrespectively of the specific physical or chemical process the system undergoes, be it adding a small fraction of an electron to the system, exciting the system, inducing transfer of charge or even bringing two subsystems together. In a sense, the complex step structure of the potential is the price one pays for the decision to describe an interacting manyelectron system via a non-interacting system with a multiplicative potential [101]. An expression for the height of the step in the exact KS potential can be derived from the changes in the LEIP.
By analyzing the exact KS potential, we show the general relationship between the step structures in the potential and derivative discontinuities in the xc energy: in the cases discussed here, the many-electron energy difference equals the corresponding KS energy difference plus the associated derivative discontinuity.
The well-known derivative discontinuity of the xc energy (∆) of a system with a varying number of electrons relates the fundamental gap and the Kohn-Sham (KS) gap: E g = E KS g + ∆. This relationship manifests in the potential as a uniform shift as the system's electron number infinitesimally surpasses an integer value. For a small finite fraction of an additional electron, spatial step structures form in the exact xc potential on the periphery of the system in order to elevate the level of the potential in the center by ∆; as this additional amount of electron tends to zero the plateau created by the steps becomes the uniform shift.
The relationship between a particular step structure in the xc potential and derivative discontinuities is not always straightforward. The infamous interatomic step, S, which forms in a stretched diatomic molecule in order to correctly distribute the electron density throughout the system has usually been regarded as unrelated to the derivative discontinuity because the system typically consists of a fixed number of electrons and the height of the step is seemingly unrelated to the ∆'s of any of the constituent atoms. We demonstrate that upon the transfer of charge from one atom to another within the diatomic molecule, the acceptor atom experiences a shift which corresponds to ∆ of that atom owing to the 'local number of electrons' on that atom surpasses an integer, ∆ a . Simultaneously the donor atom experiences a shift which corresponds to the charge-transfer derivative discontinuity (CTDD) [62], ∆ CT d→a .
We demonstrate that this discontinuity occurs within the exact KS potential within ensemble DFT of a system which undergoes charge transfer when excited. Analysis of this potential, which is studied here for the first time, can offer valuable insight for the development of advanced approximations to the xc energy within ensemble DFT. In this case, we show that S = ∆ a − ∆ CT d→a and hence the interatomic step is comprised of two derivative discontinuities, which are revealed when charge transfer occurs. In addition, this derivative discontinuity occurs when a fraction of an electron is added to the overall system while the additional charge localizes on one of the atoms. In both cases ∆ CT is related to the discontinuity of the derivative of the xc energy of the stretched molecule.
We also show that the many-electron excitation energy from the ground to the first excited state is related to the KS energy difference plus the associated derivative discontinuity [64]. We demonstrate this numerically for a single atom and show that this excitation is well approximated by the ground-state KS energy differences for this system alone, i.e., in this case the ∆ is small. This implies that the Hartree-xc kernel plays a small role in yielding accurate spectra for our single atom. This is not the case for the charge-transfer system, however, as we typically find the CTDD to be large. Hence in this case the Hxc kernel must have important features which, at least in part, correspond to the CTDD in the potential. Capturing these features in approximations to the ground-state and excited xc potential of DFT and ensemble DFT respectively, as well as the xc kernel of time-dependent DFT, is crucial for accurately obtaining many-electron excitation energies from KS theory.
Finally, we demonstrate that step structures are obtainable also from approximate xc functionals, as simple as the LDA. With the 'inverted LDA' (invLDA) approach introduced here, we construct an ensemble of LDA densities with integer number of electrons for each. Upon 'reverse-engineering' these densities we find that the corresponding potential possesses step structures, which resemble those present in the exact potential. Ensuring that our invLDA potentials obey the IP theorem, we establish the relationship between the step height and the derivative discontinuity in approximate xc functionals.
Supplemental Material for the article: "From Kohn-Sham to many-electron energies via step structures in the exchange-correlation potential" In this document we provide more details as to the inversion procedure of densities obtained for atoms and ions within common exchange-correlation (xc) approximations. The key results are presented in the main text, Section VII. The data provided here complements the main results, giving more technical details and peripheral information.
For the Li ion with N = N 0 + α electrons, where N 0 = 2 and α ∈ [0, 1], one obtains the Kohn-Sham (KS) potentials depicted in Fig. 1, by means of numerical inversion. The figure shows raw data, prior to any alignment procedure, as detailed below. The inversion is performed on the ensemble density n(r; N ) = (1−α)n(r; N 0 )+αn(r; N 0 +1), where the integer-number densities n(r; N 0 ) and n(r; N 0 + 1) were obtained within the local density approximation (LDA). As mentioned in the main text, conversion of the inversion procedure was required for r < 30 Bohr, hence the numerical artefact that is observed for higher values of r.
As it always happens in cases of numerical inversion, the resultant potentials are obtained up to a constant. The first choice of such a constant (for each value of α) presented in the main text is such that the KS potential approaches 0 at infinity. The way this require- * These authors contributed equally ment is enforced here is not by increasing L (the high limit of the variable r) to very high values, but rather by analysing the asymptotic behaviour of the resultant KS potentials. Fortunately, in our case the asymptotic behaviour is rather clear. Figure 2 shows the Hartree-exchange-correlation (Hxc) potentials (obtained by subtracting the external potential, −Z/r, from all the potentials of Fig. 1), at far distances, as a function of 1/r. It is easy to see that in the converged region, r < 30 Bohr, i.e., 1/r > 0.033 Bohr −1 , the potential v Hxc (r) ≈ a · 1 r + b, with a and b being different for each α. The values of a and b were found by linear fitting of the potentials to a straight line (in terms of 1/r), at two points: r = 20 and 30 Bohrs (equivaent of 1/r = 0.033 and 0.05 Bohr −1 ), denoted in Fig. 2 by two vertical lines. These values of a and b are given in Table I.
The fact that the Hxc potentials obtained by inversion (invLDA) decay with a power law, and not expoenetially, as is usually expected in (semi-)local functionals, such as the LDA, is a significant improvement, and is attributed by us to the fact that the piecewise-linearity criterion for the density is enforced by the above inversion procedure.
Alignment of the KS potentials of Fig. 1 by subtracting from each of the potentials the corresponding constant b     given in Table I and further subtracting the KS potential for N = 2, provides Fig. 11 of the main text. As α → 0 + , the invLDA KS potential forms a plateau around the origin, whose height approaches a constant value and whose width broadens, for the plateau to eventually fill all space and appear as a uniform shift. The dependence of the plateau height S and width R 0 on α are given in Figs. 3 and 4. The step width is determined manually from the graphs for the KS potentials, observing the position of the dip in the potential, with the accuracy of ∼ 0.03 Bohr. From Fig. 3 it is clear that the width of the plateau grows logarithmically, as the number of electrons approaches an integer from above. The convergence of the step height S seems to be proportionate to 1/ ln(α).
Back to the alignment of the KS potentials, the second choice for alignment presented in the main text is such that for each value of α the system satisfies the ionization potential (IP) theorem, namely that the highest occupied (ho) KS energy level, ε ho (α), equals the negative of the LDA IP calculated from total energy differences, for all α. In the case of Li, I LDA = E(Li + ) − E(Li) = −7.142178 − (−7.334610) = 0.192432 Hartree. Figure 5 shows with green x's the ho energy levels for the invLDA potentials given in Fig. 11 of the main text (and described here as the first alignment choice) along with the ho energy levels obtained from the standard LDA runs for systems with fractional N (dark red circles) and the reference level of (the negative of) the LDA IP (solid blue line). Interestingly, unlike the standard LDA results, which change significantly with α, the invLDA results show a strictly flat behaviour, as required, although of the wrong height. Therefore, to satisfy the IP theorem, all KS potentials for N ∈ (2, 3) are realigned (after the first alignment procedure detailed above) by just subtracting 0.086832 Hartree. Similarly, for N ∈ (1, 2) the constant 0.681309 Hartree has to be subtracted (see Fig. 6). The resultant graph is given in the main text, Fig. 13. Performing calculations with the local spin-density approximation (LSDA), namely, associating 1 + α electrons with the ↑-channel and 1 electron with the ↓-channel, we find that v ↑ KS (r) experiences a plateau, whereas no sharp features are found in v ↓ KS (r). Aligning the ↑-potentials such that they all approach zero at infinity, in the same manner detailed above for the LDA, one obtains Fig. 7. The optimal inversion parameters for the LSDA runs were found to be: p ↑ = 0.5, µ ↑ = 0.05, p ↓ = 9.0, µ ↓ = 0.03.
Similar is the situation for the Perdew-Burke-Ernzerhof (PBE) Generalized Gradient Approximation (GGA): for the ↑-channel a plateau is found and shown in Fig. 8; no plateaus are found for the ↓-channel. The optimal inversion parameters are the same as in LSDA runs.