{'fraction_non_alphanumeric': 0.03935118255941711, 'fraction_numerical': 0.009165558318394726, 'mean_word_length': 4.476959619952494, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Excitations and benchmark ensemble density functional theory for two electrons

A new method for extracting ensemble Kohn-Sham potentials from accurate excited state densities is applied to a variety of two electron systems, exploring the behavior of exact ensemble density functional theory. The issue of separating the Hartree energy and the choice of degenerate eigenstates is explored. A new approximation, spin eigenstate Hartree-exchange (SEHX), is derived. Exact conditions that are proven include the signs of the correlation energy components, the virial theorem for both exchange and correlation, and the asymptotic behavior of the potential for small weights of the excited states. Many energy components are given as a function of the weights for two electrons in a one-dimensional flat box, in a box with a large barrier to create charge transfer excitations, in a three-dimensional harmonic well (Hooke's atom), and for the He atom singlet-triplet ensemble, singlet-triplet-singlet ensemble, and triplet bi-ensemble.

Ensemble DFT (EDFT) [18][19][20][21] is one such alternative approach. Unlike TDDFT, it is based on an energy variational principle. An ensemble of monotonically decreasing weights is constructed from the M + 1 lowest levels of the system, and the expectation value of the Hamiltonian over orthogonal trial wavefunctions is minimized by the M + 1 exact lowest eigenfunctions. 19 A one-toone correspondence can be established between ensemble densities and potentials for a given set of weights, providing a Hohenberg-Kohn theorem, and application to noninteracting electrons of the same ensemble density yields a Kohn-Sham scheme with corresponding equations. 20 In principle, this yields the exact ensemble energy, from which individual excitations may be extracted.
But to make a practical scheme, approximations must be used. [22][23][24][25][26] These have been less successful for EDFT than those of ground-state [27][28][29][30][31] and TDDFT, 6,32 and their accuracy is not yet competitive with TDDFT transition frequencies from standard approximations. Some progress has been made in identifying some major sources of error. [33][34][35] To help speed up that progress, we have developed a numerical algorithm to calculate ensemble Kohn-Sham (KS) quantities (orbital energies, energy components, potentials, etc.) essentially exactly, 36 from highly accurate excited-state densities. In the present paper, we provide reference KS calculations and results for two-electron systems under a variety of conditions. The potentials we find differ in significant ways from the approximations suggested so far, hopefully leading to new and better approximations.
To illustrate the essential idea, we perform calculations on simple model systems. For example, Sec. VI A presents two 'electrons' in a one-dimensional box, repelling one another via a (slightly softened) Coulomb repulsion. In Fig. 1, we show their ground-and excitedstate densities, with I indicating the specific ground or excited state. We also plot the ensemble exchangecorrelation potentials for equally weighted mixtures of the ground and excited states, which result from our inversion scheme. In this lower plot, I = 1 denotes the ground-state exchange-correlation potential, and I > 1 indicates the potential corresponding to an equal mixture of the ground state and all multiplets up to and including the I-th state. Excitation energies for all these states are extracted using the EDFT methods described below.
The paper is laid out as follows. In the next section, we briefly review the state-of-the-art for EDFT, introducing our notation. Then we give some formal considerations about how to define the Hartree energy. The naive definition, taken directly from ground-state DFT, introduces spurious unphysical contributions (which then must be corrected-for) called 'ghost' corrections. 33 We also consider how to make choices among KS eigenstates when they are degenerate, and show that such choices mat- ter to the accuracy of the approximations. We close that section by showing how to construct symmetry-projected ensembles.
In the following section, we prove a variety of exact conditions within EDFT. Such conditions have been vital in constructing useful approximations in groundstate DFT. 30,37 Following that, we describe our numerical methods in some detail.
The results section consists of calculations for quite distinct systems, but all with just two electrons. The onedimensional flat box was used for the illustration here, which also gives rise to double excitations. A box with a high, asymmetric barrier produces charge-transfer excitations. Hooke's atom is a three-dimensional system, containing two Coulomb-repelling electrons in a harmonic oscillator external potential. 38 It has proven useful in the past to test ideas and approximations in both groundstate and TDDFT calculations. 39 We close the section reporting several new results for the He atom, using ensembles that include low-lying triplet states. Atomic units [e = = m e = 1/(4πǫ 0 ) = 1] are used throughout unless otherwise specified.

A. Basic theory
The ensemble variational principle 19 states that, for an ensemble of the lowest M + 1 eigenstates Ψ 0 , . . . , Ψ M of the HamiltonianĤ and a set of orthonormal trial func-tionsΨ 0 , . . . ,Ψ M , when the set of weights w m satisfies and E m is the eigenvalue of the mth eigenstate ofĤ. Equality holds only forΨ m = Ψ m . The density matrix of such an ensemble is defined bŷ where W denotes the entire set of weight parameters. Properties of the ensemble are then defined as traces of the corresponding operators with the density matrix. The ensemble density n W (r) is and the ensemble energy E W is n W (r) is normalized to the number of electrons, implying M m=0 w m = 1. A Hohenberg-Kohn (HK) 1 type theorem for the oneto-one correspondence between n W (r) and the potential inĤ has been proven, 18,20 so all ensemble properties are functionals of n W (r), includingD W . The ensemble HK theorem allows the definition of a non-interacting KS system, which reproduces the exact n W (r). The existence of an ensemble KS system assumes ensemble v-representability. EDFT itself, however, only requires ensemble non-interacting N -representability, since a constrained-search formalism is available. 20 Ensemble N -and v-representability are not yet proven, only assumed.
As in the ground-state case, only the ensemble energy functional is formally known, which is where v(r) is the external potential. The ensemble universal functional F W is defined as whereT andV ee are the kinetic and electron-electron interaction potential operators, respectively. The ensemble variational principle ensures that the ensemble energy functional evaluated at the exact ensemble density associated with v(r) is the minimum of this functional, Eq. (5). The ensemble KS system is defined as the noninteracting system that reproduces n W (r) and satisfies the following non-interacting Schrödinger equation: The ensemble KS system has the same set of w m as the interacting system. There is no formal proof for this consistency, and it has non-trivial implications even for simple systems. This will be explored more in Sec. II B. The KS density matrixD s,W iŝ where Φ m are non-interacting N -particle wavefunctions, usually assumed to be single Slater determinants formed by KS orbitals φ j,W . We find that this choice can be problematic, and it will be discussed in Sec. III A. The ensemble density n w (r) is reproduced by the KS system, meaning where n m (r) = Ψ m |n(r)| Ψ m , and n s,m (r) = Φ m |n(r)| Φ m . The KS densities of the individual states are generally not related to those of the interacting system; only their weighted sums are equal, as in Eq. (10). = tr{D S,WT } + d 3 r n(r)v(r) where only the ensemble exchange-correlation (XC) energy E XC,W is unknown. The form of v S,W (r) is then determined according to the variational principle by re- where where w ∈ [0, 1/M I ]. In this ensemble, here called GOK, I denotes the set of degenerate states (or 'multiplet') with the highest energy in the ensemble, g I is the multiplicity of the I-th multiplet, and M I is the total number of states up to the I-th multiplet. GOK ensembles must contain full sets of degenerate states to be welldefined. The weight parameter w interpolates between two ensembles: the equiensemble up to the I-th multiplet (w = 1/M I ) and the equiensemble up to the (I − 1)-th multiplet (w = 0). All previous studies of EDFT have been based on this type of ensemble. The purpose of EDFT is to calculate excited-state properties, not ensemble properties. With the GOK ensemble, the excitation energy of multiplet I from the ground state, ω I , is obtained using ensembles up to the I-th multiplet as which simplifies to for finite r, where ∆v XC is the change in the KS highestoccupied-molecular-orbital (HOMO) energy between w = 0 (ground state) and w → 0 + . 41 ∆v XC is a property of electron-number-neutral excitations, and should not be confused with the ground-state derivative discontinuity ∆ XC , which is related to ionization energies and electron affinities. 42

B. Degeneracies in the Kohn-Sham system
Taking the He atom as our example, the interacting system has a non-degenerate ground state, triply degenerate first excited state, and a non-degenerate second excited state. However, the KS system has a four-fold degenerate first excited state (corresponding to four Slater determinants), due to the KS singlet and triplet being degenerate (Fig. 2). Consider an ensemble of these states with arbitrary weights. Represent the ensemble energy functional Eq. (5) as the KS ensemble energy, E S,W , plus a correction, G W . The correction then must encode the switch from depending only on the sum of the weights of the excited states as a whole in the KS case to depending on the sum of triplet weights and the singlet weight separately.˜ For the interacting system, the ensemble energy and density take the forms where ω i = E i −E 0 , and so on, w T is the sum of the triplet weights, and w S is the singlet weight. On the other hand, for the KS system we have The functional showing that the exact ensemble energy functional (which can also be decomposed as in Eq. (11)) has to encode the change in the multiplet structure between noninteracting and interacting systems, even for a simple system like the He atom. Such information is unknown a priori for general systems, and can be very difficult to incorporate into approximations. This problem can be alleviated if the degeneracies are the result of symmetry. This will be discussed in Sec. III C.

C. Approximations
Available approximations to the ensemble E XC include the quasi-local-density approximation (qLDA) functional 21,43 and the 'ghost'-corrected exact exchange (EXX) functional. 23, 33 The qLDA functional is based on the equiensemble qLDA, 43 and it interpolates between two consecutive equiensembles: 21 (20) where v eqLDA XC (r) is the equiensemble qLDA functional defined in terms of finite-temperature LDA in Ref. 43.
The ensemble Hartree energy is defined analogously to the ground-state Hartree energy as shown in Eq. (11).
Similarly, Nagy provides a definition of the exchange energy for bi-ensembles: 23 where n σ (r, r ′ ) is the reduced density matrix defined analogously to its ground-state counterpart, assuming a spin-up electron is excited in the first excited state: , L↑= N ↑ + 1 and H↑= N ↑ , the spin-up lowest-unoccupied-molecularorbital (LUMO) and HOMO, respectively. Both E H in Eq. (11) and (21) contain 'ghost' terms, 33 which are cross-terms between different states in the ensemble due to the summation form of n w (r) in Eq. (4) and n w (r, r ′ ) in Eq. (22). An EXX functional is obtained after such spurious terms are corrected. For two-state ensembles, the GPG XC energy functional is then where w = 1−w. These 'ghost' corrections are small compared to the Hartree and exchange energies. However, they are large corrections to the excitation energies, as Eq. (14) contains energy derivatives instead of energies. Table I shows a few examples.
With the help of the exact ensemble KS systems to be presented in this paper, we notice inconsistencies with the GPG functional. These problems and our proposed solutions will be explained in Secs. III A and III B.

III. THEORETICAL CONSIDERATIONS
In this section, we extend EDFT to improve the consistency and generality of the theory.

A. Choice of Hartree energy
The energy decomposition in Eq. (11) is analogous to its ground-state counterpart. However, unlike T S and V , the choices for E H and E X and E C are ambiguous; only their sum is uniquely determined. As shown in Eq. (11) and (21), definitions for E H and E X can introduce 'ghost' terms. Corrections can be considered either a part of E H and E X or a part of E C . Such correction terms also take a complicated form when generalized to multi-state ensembles.
A more natural way of defining E H and E XC for ensembles can be achieved by considering the purpose of this otherwise arbitrary energy decomposition. In the ground-state case, the electron-electron repulsion reduces 44 to the Hartree energy for large Z, which is a simple functional of the density. The remaining unknown, E XC (and its components E X and E C ), is a small portion of the total energy, so errors introduced by approximations to it are small.
For ensembles, we propose a slightly different energy decomposition. Instead of defining E H and E X in analogy to their ground-state counterparts, we first define the combined Hartree-exchange energy E HX , which is the more fundamental object in EDFT. E HX can be explicitly represented as the trace of the KS density matrix: For the ground state, both Hartree and exchange contributions are first-order in the adiabatic coupling constant, while correlation consists of all higher-order terms. According to the definition above, we retain this property in the ensemble. Eq. (24) contains no 'ghost' terms by definition, eliminating the need to correct them. As a consequence, the correlation energy, E C , is defined and decomposed as where This form of E HX reveals a deeper problem in EDFT. As demonstrated in Sec. II B, the multiplet structure of real and KS He atoms is different. Real He has a triplet state and a singlet state as the first and second excited states, but KS He has four degenerate single Slater determinants as the first excited states. Worse, the KS single Slater determinants are not eigenstates of the total spin operatorŜ 2 , so their ordering is completely arbitrary. The KS system is constructed to yield only the real spin densities, not other quantities. KS wavefunctions that are not eigenstates ofŜ 2 do not generally affect commonly calculated ground-state DFT properties, 45 but things are clearly different in EDFT. Consider the biensemble of the ground state and the triplet excited state of He. Then E HX,w [n w ] depends on which three of the four KS excited-state Slater determinants are chosen, though it must be uniquely defined. Therefore, we choose the KS wavefunctions in EDFT to be linear combinations of the degenerate KS Slater determinants, preserving spatial and spin symmetries and eliminating ambiguity in E HX . Such multi-determinant, spin eigenstates are also required for construction of symmetry-projected ensembles, as described in Sec. III C.
The multi-determinant KS eigenstates and ensemble E HX proposed here avoid the errors in 'ghost'-corrected EXX, 33 which introduces spurious spin-polarization in closed-shell systems and inherent ambiguity in the treatment of triplet states. We observe considerable improvement in the first excitation energies of some atoms, as reported in the fourth line of Table I.
With E HX fixed, the definitions of E H and E X depend on one another, but E C does not. Defining a Hartree functional in the same form as the ground-state we can examine different definitions for the GOK ensemble. We can define a 'ghost'-free ensemble Hartree, E ens H , as i.e., the ensemble sum of the Hartree energies of the interacting densities, or the slightly different i.e., the ensemble sum of the Hartree energies of the KS densities. The traditional Hartree definition, introduces 'ghost' terms through the fictitious interaction of ground-and excited-state densities. Traditional and ensemble definitions differ in their production of 'ghosts,' as well as in their w-dependence. The 'ghost'-corrected As mentioned previously, the 'ghost'-corrected EXX of Ref. 33 introduces spurious spin-polarization even for closed-shell systems. The root of this problem is the use of single-Slater-determinant wavefunctions, which are not symmetry eigenstates. We have now identified E HX as being more consistent with the EDFT formalism than E H and E X . Having also justified multideterminant ensemble KS wavefunctions, we now derive a spin-consistent EXX potential, the symmetry-eigenstate Hartree-exchange (SEHX). Define the two-electron repulsion integral and L µνκλ = (µν | κλ) δ σµ,σκ δ σν ,σ λ .
φ µ (r) denotes the µ-th KS orbital and σ µ its spin state. If the occupation of the p-th Slater determinant of the µ-th KS orbital of theĩ-th multiplet of the exact system is f whereg(ĩ) is the KS multiplicity of the i-th multiplet, and C's are the coefficients of the multi-determinant wavefunctions defined by in order to write Then, if the Hartree-exchange energy for up to the I-th multiplet is where g(i) is the exact multiplicity of the i-th multiplet.
The v HX,W potential is then v SEHX HX,W ,σ (r) = δE HX,W δn W,σ (r) which yields an OEP-type equation for v HX,W (r).
The v HX,W (r) of Eq. (39) produces neither 'ghost' terms nor spurious spin-polarizations. For closed-shell systems, Eq. (39) yields v HX,W,↑ (r) = v HX,W ,↓ (r), unlike Ref. 33. An explicit v HX,W (r) can be obtained by applying the usual Krieger-Li-Iafrate(KLI) 48 approximation. Here we provide the example of the singlet bi-ensemble studied in our previous paper. 36 E HX for a closed-shell, singlet ensemble is where n orb j (r) = |φ j (r)| 2 is the KS orbital density. Spin is not explicitly written out because the system is closedshell. After applying the KLI approximation, we obtain v HX,w (r) = 1 2n w (r) with Eq. (41) is an integral equation for v HX (r) that can be easily solved.
To fully understand the performance of v HX (r), selfconsistent EDFT calculations would be needed at different values of w, which is beyond the scope of this paper. We demonstrate the performance at w = 0 later in Sec. IV D.

C. Symmetry-projected Hamiltonian
The ensemble variational principle holds for any Hamiltonian. If the HamiltonianĤ commutes with another operatorÔ, one can apply toĤ a projection operator formed by the eigenvectors ofÔ. One obtains a new Hamiltonian, and the ensemble variational principle holds for this subspace ofĤ, allowing an EDFT to be formulated.
An example would be the total spin operatorŜ 2 , where and |S are its eigenvectors. Define a new Hamiltonian H 1 asĤ H 1 has the same set of eigenvectors asĤ, but the eigenvalues are 0 for the eigenvectors not having spin S. Since one can change the additive constant inĤ arbitrarily, it is always possible to make the eigenvalues of any set of spin-S eigenvectors negative and thus ensure that they are the lowest energy states ofĤ 1 . The ensemble variational principle holds for ensembles of spin-S states. We have employed this symmetry argument in our previous paper 36 for a purely singlet two-state ensemble of the He atom.
A similar statement is available in ground-state DFT, allowing direct calculation of the lowest state of a certain symmetry. 49,50 The differences between the subspace and full treatments are encoded in the differences in their corresponding E XC . Thus the lowest two states within each spatial and spin symmetry category can be treated in EDFT in a two-state-ensemble fashion, which is vastly simpler than the multi-state formalism.
Since the multiplet structures of the interacting system and the KS system must be compatible, a symmetryprojected ensemble also requires a symmetry-projected KS system, which is impossible if KS wavefunctions are single Slater determinants, as discussed in Sec. III A.

IV. EXACT CONDITIONS
Here we prove some basic relations for the signs of various components of the KS scheme and a virial for the potentials. We describe a feature of the ensemble derivative discontinuity and extraction of excited properties from the ground state.

A. Inequalities
Simple exact inequalities of the energy components (such as E C < 0) have been proven in ground-state DFT. 42 If these are true in EDFT, experiences designing approximated E XC in ground-state DFT may be transferrable to EDFT. Here we show that inequalities related to the correlation energy are still valid in EDFT.
Due to the variational principle, 19 the wavefunctions that minimize the ensemble energy Eq. (5) are the interacting wavefunctions Ψ m . Thus The existence of a non-interacting KS system 20 means T S,W is the smallest possible kinetic energy for a given density n W (r), resulting in From Eq. (47) and (48) we immediately obtain and These inequalities are later verified with exact ensemble KS calculations.
According to the variational principle, the exact interacting wavefunctions Ψ m minimize the ensemble energy E W for a given set W . First-order variations of the ensemble energy therefore vanish. Thus yielding We can write the energy of an ensemble for a given set of weights W and a given coupling constant λ as where v λ W (r) is the W -dependent, one-body potential maintaining a constant density for all degrees of interaction. 37 The functional derivative of Eq. (55) with respect to n(r) is Thus − d 3 r n(r)r·∇ δF λ W δn(r) = d 3 r n(r)r·∇v λ W (r). (57) Applying the ensemble virial theorem of Eq. (54) to Eq.
Considering the KS quantities, we can insert the energy in terms of Hartree-exchange-correlation and the λdependent one-body potential into the general virial: which can easily be converted to an "unambiguous" energy density. 52

C. Asymptotic behavior
Ref. 40 derived the ensemble derivative discontinuity of Eq. (16) for bi-ensembles, in the limit of w → 0. For finite w of an atomic system, as shown in our previous paper, 36 ∆v XC is close to a finite constant for small r, and jumps to 0 at some position denoted by r C . We provide the derivation of the location of r C as a function of w here.
For atoms, the HOMO wavefunction and LUMO wavefunctions have the following behavior: with α ≥ α ′ . For the bi-ensemble of the ground state and the first excited state, the ensemble density is assuming that the HOMO is doubly-occupied. The behavior of the density at large r is dominated by the density of the doubly-occupied HOMO and the second term. In order to see where the density decay switches from that of the HOMO to the LUMO, we find the r-value at which the two differently decaying contributions are equal: As w → 0, r C is then with ∆α = α − α ′ . The ionization energies are available for the He ground state and singlet excited state. Since we obtain for the He singlet bi-ensemble with w close to 0.

D. Connection to ground-state DFT
With weights as in Eq. (13), calculation of the excitation energies is done recursively: for the M th excited state, one needs to perform an EDFT calculation with the M th state highest in the ensemble, and another EDFT calculation with the (M − 1)th as the highest state, and so on. Thus for the M th state, one needs to perform M separate EDFT calculations for its excitation energy.
For bi-ensembles, however, the calculation of the excitation energy can be greatly simplified. Eq. (14) holds for w = 0, so one can work with ground-state data only and obtain the first-excited state energy, without the need for an explicit EDFT calculation of the two-state ensemble.
We calculate the first excitation energies of various atoms and ions with Eq. (14) at w = 0 with both qLDA 21,43 (based on LDA ground states), EXX, 23 GPG, 33 and SEHX, with the last three based on OEP-EXX (KLI) ground states. 48 In order to ensure the correct symmetry in the end result, SEHX must be performed on spin-restricted ground states. However, for closed-shell systems, these results coincide with those of spin-unrestricted calculations. We use these readily available results when possible in this paper. The wderivatives of the E XC 's for qLDA and EXX required in Eq. (14) are (considering Eq. (76)) −n H (r)n L (r ′ ) + n H (r, r ′ )n L (r, r ′ ) where j sums over the spin-up densities. Only ground state properties are needed to evaluate Eq. (73). The results are listed in Table I. SEHX improves calculated excitation energies for systems where GPG has large errors, such as Be and Mg atoms.

V. NUMERICAL PROCEDURE
We invert the ensemble KS equation with exact densities to obtain the exact KS potential. We describe the numerical inversion procedure in Ref. 36 I. First non-triplet excitation energies (in eV) of various atoms and ions calculated with qLDA, EXX, GPG, and SEHX functionals. qLDA calculations were performed upon LDA (PW92) 53 ground states; EXX 23 ground states were used for the rest. Asterisks indicate use of spin-restricted ground states. qLDA relies on ground-state LDA orbital energy differences; it cannot be used with the single bound orbital of LDA He.
We implemented the numerical procedure on a realspace grid. The ensemble KS equation (8) is solved by direct diagonalization of the discrete Hamiltonian. The grid is in general nonuniform, which complicates the discretization of the KS kinetic energy operator. We tested two discretization schemes, details of which are available in the Supplemental Material. Based on these tests, all results presented in this paper have been obtained using the finite-difference representation . (74)

A. Derivative Corrections
Exactness of the inversion process can be verified by calculating the excitation energies with Eq. (14) at different w values. Eq. (14) requires calculating E XC,w of the exact ensemble KS system, Since we do not have a closed-form expression for the exact E XC , its derivative can only be calculated numerically. However, the numerical derivative of E XC , ∂E XC,w [n w ]/∂w, is not the quantity required in Eq. (14). It is related to the true derivative through The correction to the numerical derivative of E XC,w adjusts for the w-dependence of the ensemble density, which is not inherent to E XC,w . All our calculations show that the two terms on the right hand side of Eq. (76) are of the same order of magnitude. This shows that the exact E XC,w [n] changes more slowly than n w (r) as w changes. Though the calculations of E XC,w and ∂E XC,w [n]/∂w| n=nw both involve integrations containing v XC,w (r), they are independent of the additive constant.

VI. RESULTS
We apply the numerical procedure described in Sect. V to both 1D and 3D model systems in order to further demonstrate our method for inverting ensemble densities.

A. 1D flat box
The external potential of the 1D flat box is The exact wavefunctions can be solved numerically for two electrons with the following soft-Coulomb interaction: where we choose a = 0.1.  Table II shows the total and kinetic energies of the exact ground state and first four excited states for L = 1 a.u., calculated on a 2D uniform grid with 1000 points for each position variable. The third excited state is a doubly-excited state corresponding to both electrons occupying the second orbital of the box. Fig. 1 shows the exact densities of the ground state and first four excited states, together with the XC potential of equiensembles containing 1 to 5 multiplets. Table III lists calculated excitation energies, showing that the excitation energy is independent of w, no matter how many states are included in the ensemble. This is a non-trivial exact condition for the ensemble E XC .
Double excitations are generally difficult to calculate. It has been shown that adiabatic TDDFT cannot treat double or multiple excitations. 12 Table III shows that there is no fundamental difficulty in treating double excitations with EDFT. Fig. 1 shows that v XC,w (r) for the 4multiplet equiensemble resembles the potentials of other ensembles. The exact two-multiplet ensemble XC potentials at different w are plotted in Fig. 4. The bump up near the center of the box in these potentials ensures that the ensemble KS density matches that of the real ensemble density. Increasing the proportion of the excited state density (see Fig. 1) included in the ensemble density requires a corresponding increase in the height of this bump (see Supplemental Material). With no asymptotic region, there is no derivative discontinuity for the box, and v XC,w→0 (r) is equal to the ground-state v XC (r). Energy components for the bi-ensemble of the 1D box satisfy the inequalities shown in Sec. IV A and are reported in the Supplemental Material. there is no fundamental difficulty in treating CT excitations with EDFT. Since EDFT calculations do not involve transition densities, they do not suffer from the lack-of-orbital-overlap problem in TDDFT.
The external potential for the CT box is with the barrier dimensions chosen for numerical stability of the inversion process. The lowest two eigenstate densities are given in the top of Fig. 5. The groundstate and first-excited-state total and kinetic energies of the CT system described are This significant increase in kinetic energy together with a small total energy change designate the CT character of the first excited state. The electrons become distributed between the two wells of the potential, instead of being confined in one well.
The ground-and first-excited-state densities and ensemble XC potentials are plotted in Fig. 5. The potentials show the characteristic step-like structures of charge-transfer excitations, which align the chemical potentials of the two wells. 55,56

C. Hooke's atom
Hooke's atom is a popular model system 57,58 with the following external potential: For our calculation, k = 1/4. Though the first excited state has cylindrical symmetry, we use a spherical grid, as it has been shown that the error due to spherical averaging is small. 59 As a closed-shell system, the spatial parts, and therefore the densities, of the spin-up and spin-down ensemble KS orbitals have to be the same, so we treat this system as a bi-ensemble. The magnitude of the external potential of the Hooke's atom is smallest at r = 0, and becomes larger as r increases. This is completely different from the Coulomb potential of real atoms. Since the electron-electron interaction is still coulombic, v XC (r) can be expected to have a −1/r behavior as r → ∞, which is negligibly small compared to v(r). Combined with a density that decays faster than real atomic densities, n(r) ∼ exp(−ar 2 ) versus n(r) ∼ exp(−br), convergence of the Hooke's atom v XC (r) is difficult in the asymptotic region. Additionally, v XC (r) ≫ v(r) for small r, so larger discretization errors in this region also contribute to poorer inversion performance. Despite these challenges, we still obtain highly accurate excitation energies.
A logarithmic grid with 550 points ranging from r = 10 −5 a.u. to 10 a.u. is used for all the Hooke's atom calculations. On this grid, the exact ground-and first excitedstate energies are Calculated ω 2 was 9.786 eV for all values of w tested (see Supplemental Material). Unlike the He atom and the 1D flat box, the n w (r) and v XC,w (r) show little variation with w (see Supplemental Material). The second KS orbital of the Hooke's atom is a p-type orbital, which has no radial node and a radial shape similar to that of the first KS orbital. Consequently, the changes in the KS and xc potentials are also smaller.

D. He
Using the methods in Ref. 36, we employ a Hylleraas expansion of the many-body wavefunction 60 to calculate highly accurate densities of the first few states of the He atom. We report the exact ensemble XC potentials for He singlet ensemble in that paper. Table V shows accurate excitation energies calculated from mixed symmetry, three-multiplet, and strictly triplet ensembles, demonstrating the versatility of EDFT. Fig. 6 compares v XC,w (r) for four types of He equiensembles, highlighting their different features. The characteristic bump up in these potentials is shifted left in the 2-multiplet case, relative to the others shown. This shift has little impact on the first "shell" of the ensemble density's shell-like structure, but the second is shifted left and has sharper decay, noticeably different from that of the singlet ensemble. The inequalities shown in Sec. IV A and the virial theorem Eq. (64) are verified by the exact results. Behaviors of the energy components for the singlet ensemble versus w are plotted in Fig. 7. Correlation energies show strong non-linear behavior in w. According to Eq. (14), the excitation energies are related to the derivative of E XC versus w. Therefore, E C is crucial for accurate excitation energies, even though its absolute magnitude is small.

VII. CONCLUSION
This paper is an in-depth exploration of ensemble DFT, an alternative to TDDFT for extracting excitations from DFT methodology. Unlike TDDFT, EDFT is based on a variational principle, and so one can expect that the failures and successes of approximate functionals should occur in different systems than those of TDDFT.
Apart from exploring the formalism and showing several new results, the main result of this work is to apply a new algorithm to highly-accurate densities of eigenstates to explore the exact EDFT XC potential. We find intriguing characteristic features of the exact potentials that can be compared against the performance of old and new approximations. We also extract the weightdependence of the KS eigenvalues, which are needed to extract accurate transition frequencies, and find that a large cancellation of weight-dependence occurs in the exact ensemble. Many details of these calculations are reported in the supplemental information.
From the original works of Gross, Oliviera, and Kohn, ensemble DFT has been slowly developed over three decades by a few brave pioneering groups, most prominently that of Nagy. We hope that the insight these exact results bring will lead to a plethora of new ensemble approximations and calculations and, just possibly, a competitive method to treating excitations within DFT.