key: cord-0806155-fm6g9mwc
authors: Cernatic, Filip; Senjean, Bruno; Robert, Vincent; Fromager, Emmanuel
title: Ensemble Density Functional Theory of Neutral and Charged Excitations
date: 2021-09-10
journal: Topics in current chemistry
DOI: 10.1007/s41061-021-00359-1
sha: 4ec5e5975d5a120416e014094232f5f23ae6a8e0
doc_id: 806155
cord_uid: fm6g9mwc

Recent progress in the field of (time-independent) ensemble density-functional theory (DFT) for excited states are reviewed. Both Gross-Oliveira-Kohn (GOK) and $N$-centered ensemble formalisms, which are mathematically very similar and allow for an in-principle-exact description of neutral and charged electronic excitations, respectively, are discussed. Key exact results like, for example, the equivalence between the infamous derivative discontinuity problem and the description of weight dependencies in the ensemble exchange-correlation density functional, are highlighted. The variational evaluation of orbital-dependent ensemble Hartree-exchange (Hx) energies is discussed in detail. We show in passing that state-averaging individual exact Hx energies can lead to severe (solvable though) $v$-representability issues. Finally, we explore the possibility to use the concept of density-driven correlation, which has been recently introduced and does not exist in regular ground-state DFT, for improving state-of-the-art correlation density-functional approximations for ensembles. The present review reflects the efforts of a growing community to turn ensemble DFT into a rigorous and reliable low-cost computational method for excited states. We hope that, in the near future, this contribution will stimulate new formal and practical developments in the field.

Kohn-Sham density-functional theory (KS-DFT) [1] has become over the last two decades the method of choice for modeling the electronic structure of molecules and materials. This success originates from the relatively low computational cost of the method and its relatively good accuracy in the description of ground-state properties such as equilibrium structures and activation barriers. KS-DFT is an in-principle-exact ground-state theory. As such, it cannot be used straightforwardly for calculating excited-state properties. The formal beauty of KS-DFT lies in its universal description of electronic ground states. Indeed, in KS-DFT, all the quantum many-electron effects are encoded into a system-independent exchange-correlation (xc) density functional E xc [n] which, if it were known, would allow us to compute the exact ground-state energy of any electronic system, simply by solving self-consistent one-electron equations, instead of the many-electron Schrödinger equation. This universality is somehow lost, at least partially, when turning to the excited states. An obvious reason is that there are different types of electron excitations. First we should distinguish charged excitations, where the number of electrons in the system is modified, from neutral excitations, which occur at a fixed electron number [2] . In the context of DFT, these two processes are usually approached very differently. In the case of charged excitations, one traditionally refers to the extension of DFT to fractional electron numbers [3, 4, 5, 6, 7, 8, 9] . Its implementation at the simplest (semi-) local xc density functional level of approximation usually yields too small fundamental gaps for solids [10] . This can be related to the discontinuities that the density functional derivative of the xc energy should in principle exhibit, when crossing an integer electron number, but that are completely absent from standard approximations. This is the reason why hybrid functionals (where a fraction of orbital-dependent exchange energy is combined with density functionals) [11, 12, 8, 13] or even more involved frequency-dependent post KS-DFT approaches like GW [14, 15, 16, 17, 18, 19, 20, 21] are usually employed for improving the description of fundamental gaps, which implies a substantial increase in computational cost. If we now turn to neutral excitation processes, the most popular approach is (linear response) time-dependent-DFT (TD-DFT) [22, 23] . In TD-DFT, the excitation energies are determined by searching for the poles of the KS linear response function. The success of TD-DFT lies in the fact that in many (but not all) cases the (rather crude) adiabatic approximation performs relatively well with a moderate computational cost. In the latter approximation, the time-dependent density-functional xc energy (it is in fact an action [24, 25] , to be more precise) is evaluated, in the time range t 0 ≤ t ≤ t 1 , from the regular (time-independent) ground-state xc functional and the density n(t) ≡ n(r, t) at time t as follows, A xc [n] ≈ t1 t0 dt E xc [n(t)] [24, 25, 26, 23] . Nevertheless, the description of charge-transfer excitations becomes problematic when semi-local functionals are employed [23, 27, 28, 29] . Moreover, multiple excitations are completely absent from standard TD-DFT spectra, precisely because of the adiabatic approximation [30, 23, 29] . Let us finally mention that TD-DFT can be seen as a single-reference post-DFT method, because it relies on the single-configuration KS ground-state wave function. This can be problematic, for example, when the system under study has near-degenerate low-lying states, like in the vicinity of a conical intersection, for example [23] .

The various limitations (in terms of computational cost, accuracy, or physics) of the above-mentioned frequency-dependent post-DFT approaches explain why, in recent years, time-independent formulations of DFT for excited states have attracted an increasing attention. For charged excitations, the extended Koopman's theorem turns out to be an appealing alternative [31, 32] . DFT for fractional electron numbers has also been further developed in the last decade (see, for example, Refs. [33, 34, 35, 36, 37, 38, 39, 40] ). It has also been argued that restricting to integer electron numbers, in the calculation of charged excitations, is a completely valid alternative [41, 42] . For neutral excitations, various state-specific approaches have been explored at the formal [43, 44] but also practical levels. In the latter case, orbital optimizations must be performed under proper constraints in order to avoid variational collapses. This leads to various variational computational schemes such as the ∆ self-consistent field (∆SCF) approach [45, 46, 47] , the maximum overlap method (MOM) [48, 49, 50, 51, 52] , orthogonally-constrained DFT [53, 54, 55, 56] , and constricted DFT [57, 58, 59, 60, 61] . Interestingly, when a system has a Coulombic one-electron external local potential, which is the case for any real molecule, an excited state can be identified directly from its density [62, 63, 64, 65] . This fundamental property can be used for constructing an in-principle-exact DFT for individual excited states. The practical implementation of such a theory is not straightforward though, in particular because density functionals must be defined also for non-Coulombic densities, so that functional derivatives can be evaluated. Another strategy, which is the main topic of this review, is ensemble DFT (eDFT). The ensemble formalism is often referred to in DFT for mathematical purposes like, for example, extending the domain of definition of density functionals [66] or describing strict degeneracies [67, 68] . It has probably be underestimated as a potential alternative to standard time-dependent methods for the practical calculation of (charged or neutral) excitations [69, 70, 71, 72, 73, 74, 75] . A clear advantage of eDFT over time-dependent approaches is that its computational cost is essentially that of a standard KS-DFT calculation. The only difference is that, in an ensemble, orbitals can be fractionally occupied. Moreover, like in TD-DFT, regular ground-state xc functionals can be recycled in eDFT [71] . Note that, unlike in thermal DFT [76, 77, 78] , the fractional orbital occupation numbers are actually known before the eDFT calculation starts. They are determined by the ensemble weights that the user (arbitrarily) assigns to the M (M = 1, 2, . . .) lowest excited states she/he wants to study. In (say canonical) thermal DFT [78] , the ensemble weights are determined not only from the temperature (that is arbitrarily fixed by the user) but also from the (to-be-calculated) KS orbital energies. Another important feature of eDFT is that it can in principle describe any kind of excitation, including the double excitations [79, 80] that standard approximate TD-DFT misses. The eDFT formalism for neutral excitations is often referred to as Gross-Oliveira-Kohn DFT (GOK-DFT) because it relies on the GOK variational principle [81, 82] . A similar formalism, referred to as N -centered eDFT, has been derived recently by Senjean and Fromager [83] for the description of charged excitations. The present review aims at highlighting recent progress in eDFT, with a particular focus on the exact theory and the development of approximations from first principles.

The chapter is organized as follows. An introduction to GOK-DFT and N -centered eDFT is given in Sec. 2. Even though the two theories describe completely different physical processes, their mathematical formulations are very similar, as highlighted in the section. We also explain how individual energy levels (which give access to excited-state properties) can be extracted, in principle exactly, from these theories. Then, in Sec. 3, we discuss the equivalence between the xc derivative discontinuity, which is a fundamental concept in DFT, and the ensemble weight derivative of the xc density functional, which is central in eDFT. Strategies for developing weight-dependent xc densityfunctional approximations (DFAs) for ensembles, which is the most challenging task in eDFT, are then reviewed. Key concepts will be illustrated with the prototypical asymmetric Hubbard dimer model [84, 85] . The rigorous construction of hybrid functionals for ensembles is discussed in Sec. 4 . We reveal that using state-averaged exact exchange energies, which is common in computational eDFT studies [75] , can lead to severe (solvable though) v-representability issues. Finally, we discuss in Sec. 5 the concept of ensemble density-driven correlation, which was recently introduced by Gould and Pittalis [86] , and how it could be used in the design of correlation DFAs for ensembles. Conclusions and perspectives are given in Sec. 6 . Detailed derivations of some key equations are provided in the appendices. 

whereĤ =T +Ŵ ee +V ext (2) is the electronic Hamiltonian within the Born-Oppenheimer approximation andT

are the M -electron kinetic energy, Coulomb repulsion, and local (i.e., multiplicative) external potential operators, respectively. Both neutral and charged excitations can be described within a unified eDFT formalism. The calculations will simply differ in the type of excited states (charged or neutral) that is included into the ensemble. On the one hand, the GOK-DFT formalism [82] will be employed for neutral excitations while, for charged excitations, we will use the more recent N -centered eDFT formalism [83] . In this section, we derive key equations for each theory, and we show how individual excited-state properties (energy levels and densities) can be extracted, in principle exactly, from the KS ensemble. Real algebra will be used throughout this work. For the sake of clarity, derivations will be detailed only for ensembles consisting of non-degenerate states. The theory obviously applies to more general cases [82, 87] .

GOK-DFT has been formulated in the end of the 1980's by Gross, Oliveira, and Kohn [81, 82, 88] and is a generalization of the equiensemble DFT of Theophilou [89] . In contrast to standard DFT, which is a ground-state theory, GOK-DFT can describe both ground and (neutral) excited states. In this context, the ensemble density is used as a basic variable (in place of the ground-state density).

Before deriving the main equations of GOK-DFT, let us introduce the exact ensemble theory. We start with the ensemble GOK energy expression [81] 

which is simply a state-averaged energy where w = (w 1 , w 2 , . . .) denotes the collection of ensemble weights that are assigned to the excited states, and E I ≡ E N I are the energies of the N -electron ground (I = 0) and excited (I > 0) states Ψ N I . We assumed in Eq. (4) that the full set of weights (which includes the weight w 0 assigned to the ground state) is normalized, i.e., w 0 = 1 − I>0 w I , so that

For ordered weights w I ≥ w I+1 ≥ 0, with I ≥ 0, the following (so-called GOK) variational principle holds [81] ,

where Ψ I is a trial set of orthonormal N -electron wave functions. Note that the lower bound E w , which is the exact ensemble energy, is not an observable. It is just an (artificial) auxiliary quantity from which properties of interest, such as the excitation energies, can be extracted. Since it varies linearly with the ensemble weights, the extraction of individual energy levels is actually trivial. Indeed, combining the following two relations [see Eq. (5)],

and

leads to

Despite its simplicity the above expression has not been used until very recently for extracting excited-state energies from a GOK-DFT calculation [90, 91] . Further details will be given in the next section.

In GOK-DFT, the ensemble energy is obtained variationally as follows [82] ,

where the minimization is restricted to N -electron densities, i.e., dr n(r) = N , and the universal GOK density functional

which is evaluated from the density-functional eigenfunctions {Ψ w I [n]} that fulfill the density constraint I w I n Ψ w I [n] (r) = n(r), is the analog for GOK ensembles of the universal Hohenberg-Kohn functional. Its construction relies on a potential-ensemble-density map that is established for a given and fixed set w of ensemble weight values. Therefore, the universality of the functional implies, like in ground-state DFT, that it does not depend on the local external potential. However, it does not mean that it is ensemble-independent and therefore applicable to any excited state. As discussed in further detail in Secs. 4 and 5, encoding ensemble dependencies into density functionals is probably the most challenging task in eDFT.

In the standard KS formulation of GOK-DFT [82] , the GOK functional is split into non-interacting kinetic and Hartree-xc (Hxc) ensemble energy contributions, by analogy with regular KS-DFT:

The non-interacting ensemble kinetic energy functional can be expressed more explicitly as follows within the constrained-search formalism [92] ,

where Tr denotes the trace,γ w = I w I |Φ I Φ I | is a trial ensemble density matrix operator that fulfills the density constraint nγw (r) ≡ Tr [γ wn (r)] = I w I n Φ I (r) = n(r), andn(r) ≡ N i=1 δ(r−r i ) is the electron density operator at position r. Combining Eqs. (10), (12) and (13) leads to the final GOK-DFT variational energy expression

where the minimization can be restricted to single-configuration wave functions (determinants or configuration state functions), hence the minimization over orbitals {ϕ p } on the first line of Eq. (15) . The minimizing KS orbitals {ϕ w p }, from which the KS wave functions {Φ w I [n w ] ≡ Φ w I } in the minimizing density matrix operatorγ w KS are constructed, fulfill the following self-consistent GOK-DFT equations,

where v w

is the ensemble Hxc density-functional potential. In the exact theory, the ensemble KS orbitals reproduce the exact (interacting) ensemble density, i.e.,

where the individual KS densities read as

and n I p is the (weight-independent) occupation number of the orbital ϕ w p in the single-configuration wave function Φ w I .

Let us now focus on the ensemble Hxc density functional. By analogy with regular KS-DFT, it can decomposed into Hx and correlation energy contribu-

In the original formulation of GOK-DFT [82] , the Hx functional is further decomposed as follows,

where

is the standard weight-independent Hartree functional, and

is the exact (complementary and weight-dependent) ensemble exchange functional. Note that Φ w I [n], which describes one of the configurations included into the ensemble, may not be a pure Slater determinant [93] . Other (weightdependent) definitions for the ensemble Hartree energy, where the explicit dependence on the ensemble density n is lost, have been explored [94] . In the most intuitive one, the ensemble Hartree energy is evaluated as the weighted sum of the individual KS Hartree energies:

For the sake of generality, we will keep in the following both Hartree and exchange energies into a single functional E w Hx [n] which is defined as

The remaining weight-dependent correlation energy can then be expressed as follows, according to Eqs. (11) and (14),

where the non-interacting KS {Φ w I [n]} and interacting {Ψ w I [n]} wave functions, which both reproduce the (weight-independent here) trial ensemble density n, whatever the ensemble weight values w, are in principle weight-dependent [95, 96] . Interestingly, the interacting density-functional wave functions lose their weight dependence when the trial density n matches the exact physical ensemble density n w , i.e., Ψ w I [n w ] = Ψ I ≡ Ψ N I . However, as shown in Sec. 5.2, the KS wave functions remain weight-dependent, even in this special case.

As readily seen from Eqs. (12) , (15) and (16), the only (but crucial) difference between regular ground-state KS-DFT and GOK-DFT is the weight dependence in the ensemble density-functional Hxc energy and potential. The computational cost should essentially be the same in both approaches. The challenge lies in the proper description of the weight-dependent ensemble Hxc density functional. Different approximations have been considered, such as the use of (weight-independent) regular ground-state functionals [78, 97] , or the use of an ensemble exact-exchange energy [93, 98, 99] with or without approximate weight-dependent correlation functionals [99, 91, 79] . Note that the expected linearity-in-weight of the ensemble energy is not always reproduced in (approximate) practical GOK-DFT calculations [97] . As a result, different weights can give different excitation energies, which is a serious issue. This lead to different computation strategies, such as trying to find an optimal value for the weights [100] , using Boltzmann weights instead [78] , restricting to equiensembles [91, 79] , or considering the ground-state w = 0 limit of the theory, like in the direct ensemble correction (DEC) scheme [87, 101] . A linear interpolation method has also been proposed [97, 102] .

Designing weight-dependent ensemble DFAs that systematically reduce the curvature in weight of the ensemble energy, while providing at the same time accurate excitation energies, is an important and challenging task. Recent progress in this matter will be extensively discussed in Secs. 4 and 5.

In Sec. 2.1.2, we have shown that both exact ensemble energy and density can be calculated, in principle exactly, within GOK-DFT. At this point we should stress that the KS and true physical densities are not expected to match individually, even though they both reproduce the same ensemble density [see Eq. (18) ]. This subtle point will be discussed in more detail in Sec. 5.2. Nevertheless, in complete analogy with Eq. (9), the exact individual densities can be extracted from the ensemble density as follows [103] ,

which, by inserting the expression in Eq. (18) , leads to the key result [103] 

where, according to Eq. (19) , the weight derivative of the individual KS densities

can be evaluated from the (static) linear response of the KS orbitals. This can be done, in practice, by solving an ensemble coupled-perturbed equation [73, 103] , for example.

Turning to the excitation energies, we obtain from the variational GOK-DFT ensemble energy expression and the Hellmann-Feynman theorem the following expression, where the derivatives of the minimizing (and therefore stationary) KS wave functions do not contribute,

with ∆γ w KS,I = |Φ w

This expression can be further simplified as follows [90] :

where E w I denotes the Ith (weight-dependent) KS energy which is obtained by summing up the energies {ε w p } of the KS orbitals that are occupied in Φ w I . Hence, the excitation energies can all be determined, in principle exactly, from a single GOK-DFT calculation.

As shown by Deur and Fromager [90] , individual energy levels can also be extracted (from the KS ensemble) and written in a compact form. For that purpose, we will use the exact expression of Eq. (9) where we see, in the light of Eq. (30) , that it is convenient to express the total ensemble energy [first term on the right-hand side of Eq. (9)] in terms of total KS energies. Levy and Zahariev (LZ) made such a suggestion in the context of regular ground-state DFT [104] . For that purpose, they introduced a shift in the Hxc potential that can be trivially generalized to GOK ensembles as follows [90] ,

Note that, if the exact LZ-shifted Hxc potential were known, we would be able to evaluate exact ensemble density-functional Hxc energies as follows,

Once the LZ shift has been applied to the ensemble Hxc potential, the (total N -electron) KS energies will be modified as follows,

and the true ensemble energy will simply read as a weighted sum of (LZshifted) KS energies:

Note that the KS excitation energies are not affected by the shift:

Thus, by combining Eqs. (9) , (30) , (34) and (35), we recover the exact expression of Ref. [90] for ground-and excited-state energy levels:

As readily seen from Eq. (36), applying the LZ shift is not sufficient for reaching an exact energy level. The ensemble weight derivatives of the Hxc density functional are also needed for that purpose.

Finally, it is instructive to consider the general expression of Eq. (36) in the ground-state w = 0 limit of the theory, which gives

where n Ψ0 is the exact ground-state density. As readily seen from Eq. (37), as we start from a pure I = 0 ground-state theory (we recover the energy expression of Levy and Zahariev in this case [104] ), the inclusion of a given I > 0 excited state into the ensemble induces an additional shift in the Hxc potential, which corresponds to the weight derivative ∂E w Hxc [n Ψ0 ]/∂w I | w=0 and can be interpreted as a derivative discontinuity, as shown in Ref. [105] and extensively discussed in Sec. 3, in the context of charged excitations.

A recent adaptation of GOK-DFT to charged excitations, which is referred to as N -centered eDFT [83] , is introduced in the present section.

The N -centered ensemble [83] can be seen as the "grand canonical" groundstate version of GOK ensembles. It is constructed from the M -electron ground states where the three possible values of the integer M ∈ {N − 1, N, N + 1} are, like the corresponding ensemble density (see below), centered in N , hence the name "N -centered". The exact N -centered ensemble energy is defined as follows [83] ,

where the two N -centered ensemble weights ξ − and ξ + , which describe the removal/addition of an electron from/to the N -electron system, respectively, are collected in

Similarly, the N -centered ensemble density reads as

Designed by analogy with GOK-DFT (which describes neutral excitations), the N -centered ensemble density integrates to the central integer number N of electrons:

In other words, even though we describe charged excitation processes, the number of electrons remains fixed and equal to the integer N whatever the value of the ensemble weights ξ. This major difference with the conventional DFT for fractional electron numbers [3] has important implications that will be discussed extensively in Sec. 3.

In this context, the ensemble energy can be determined variationally, as a direct consequence of the conventional Rayleigh-Ritz variational principle for a fixed number of electrons, i.e.,

where {Ψ M } are trial M -electron normalized wave functions, provided that the (so-called convexity) conditions ξ − ≥ 0, ξ + ≥ 0, and ξ − (N − 1)+ξ + (N + 1) ≤ N are fulfilled. Like in GOK-DFT, the ensemble energy E ξ 0 varies linearly with the ensemble weights. As a result, charged excitation energies can be extracted through differentiation with respect to the N -centered ensemble weights. For example, since

the exact fundamental gap can be determined as follows,

We can also extract the individual cationic, anionic, and neutral energies, respectively, as follows,

and

Eqs. (45)-(47) will be used in the following for deriving exact ionization potential and electron affinity theorems.

In complete analogy with GOK-DFT, the N -centered ensemble energy can be determined variationally as follows,

where, in the KS formulation of the theory [83] , the universal N -centered ensemble functional reads as

The non-interacting kinetic energy functional

is now determined through a minimization over N -centered density matrix operatorsγ

that fulfill the density constraint nγξ(r) = Tr γ ξn (r) = n(r). Combining Eqs. (48) , (49) and (50) leads to the final ensemble energy expression,

which is mathematically identical to its analog in GOK-DFT [see Eq. (15)], even though the physics it describes is completely different. The orbitals {ϕ ξ p }, from which the minimizing single-configuration KS wave functions Φ M,ξ 0 in γ ξ KS are constructed, fulfill self-consistent KS equations that are similar to those of regular (N -electron ground-state) KS-DFT:

The only difference is that the N -centered ensemble Hxc potential v ξ Hxc [n](r) = δE ξ Hxc [n]/δn(r) is now employed. In the exact theory, the ensemble KS orbitals are expected to reproduce the interacting N -centered ensemble density, i.e.,

or, equivalently [83] ,

Turning to the N -centered ensemble Hxc density functional, it can be decomposed as E ξ

, where, by analogy with GOK-DFT, the exact Hx energy is expressed in terms of the N -centered ensemble densityfunctional KS wave functions as follows,

and the complementary correlation functional reads as

+ξ + T +Ŵ ee

where {Ψ M,ξ 0

[n]} denotes the interacting density-functional N -centered ensemble.

When comparison is made with Sec. 2.1.2, it becomes clear that N -centered and GOK eDFTs are essentially the same theory (they only differ in the definition of the ensemble). From that point of view, we now have a unified eDFT for charged and neutral electronic excitations. As a result, N -centered eDFT can benefit from progress made in GOK-DFT, and vice versa.

We have shown in Sec. 2.2.1 that neutral, anionic, and cationic ground-state energies can be extracted exactly from the N -centered ensemble energy [see Eqs. (45) , (46) , and (47) ]. We can now use the variational density-functional expression of Eq. (52) to obtain expressions for the fundamental gap, the ionization potential (IP), and the electron affinity (EA). Note that these quantities are traditionally derived in the context of DFT for fractional electron numbers [3] (see Sec. 3 for a detailed comparison). According to the Hellmann-Feynman theorem, we can express the weight derivatives of the ensemble energy as follows,

where ∆ ±γ

are constructed from orbitals that fulfill the KS Eq. (53), the above energy derivative can be rewritten in terms of the KS orbital energies as [83] 

By plugging Eq. (60) into Eq. (44), we immediately obtain the following exact expression for the fundamental gap:

If we now apply the LZ shift-in-potential procedure [104] , by analogy with GOK-DFT (see Sec. 

we can express both the ensemble energy and its derivatives in terms of the LZ-shifted KS orbital energies ε ξ p , thus leading to the following compact expressions for the ensemble and individual energies [83] , respectively:

and

By subtraction, we immediately obtain in-principle-exact IP and EA theorems:

and

Interestingly, in the regular ground-state N -electron limit (i.e., when ξ = 0), the expression of Levy and Zahariev [104] is recovered for the IP,

where the asymptotic value of the LZ-shifted Hxc potential away from the system [see Ref. [104] and Eq. (132) ] can now be expressed explicitly, within the N -centered ensemble formalism, as ∂E ξ

. Similarly, we obtain the following expression for the EA:

As readily seen from the above expressions, neutral and charged systems cannot be described with the same (LZ-shifted) Hxc potential. As shown in Sec. 3, the additional ensemble weight derivative correction [second term on the righthand side of Eqs. (69) and (70)] is actually connected to the concept of derivative discontinuity which manifests in conventional DFT for fractional electron numbers, when crossing an integer [106] .

The concept of derivative discontinuity originally appeared in the context of DFT for fractional electron numbers [3] , which is the conventional theoretical framework for the description of charged excitations. The (xc functional) derivative discontinuities play a crucial role in the evaluation of fundamental gaps [106] . More specifically, they correct the bare KS gap which is only an approximation to the true interacting gap. It is well known that standard (semi)-local DFAs do not contain such discontinuities, which explains why post-DFT methods based on Green functions, for example, are preferred for the computation of accurate gaps [14, 15, 16, 17, 18, 19, 20, 21] . Their substantially higher computational cost is a motivation for exploring simpler (frequencyindependent) strategies. The recently proposed N -centered ensemble formalism [83, 107] , which has been introduced in Sec. 2.2, is (among others [34, 36, 42, 108, 109, 110, 111, 112] ) promising in this respect.

From a more fundamental point of view, it is important to clarify the similarities and differences between N -centered eDFT and the standard formulation of DFT for charged excitations, which is often referred to as Perdew-Parr-Levy-Balduz (PPLB) DFT [3] . More specifically, we should explain what the derivative discontinuity, which is central in PPLB, becomes when switching to the N -centered formalism. This is the purpose of this section.

After a brief review in Sec. 3.1 of the PPLB formalism and its implications, we will show (in Sec. 3.2), on the basis of Ref. [113] , that derivative discontinuities exist also in N -centered eDFT and that they are directly connected to the ensemble weight derivatives of the xc functional, like in GOK-DFT [105] . Finally, we will explain in Sec. 3.3 why these discontinuities can essentially be removed from the theory, unlike in PPLB, and discuss the practical implications.

3.1 Review of the regular PPLB approach to charged excitations

The key idea in PPLB is to describe electron ionization or affinity processes through a continuous variation of the electron number, hence the need for an extension of DFT to fractional electron numbers. For that purpose, the energy of an artificial (zero-temperature) grand-canonical-type ensemble, which should not be confused with physical finite-temperature grand-canonical ensembles of statistical physics [114] , is constructed as follows,

where we minimize over integer numbers M of electrons and E M 0 denotes the exact M -electron ground-state energy of the system. In this formalism, the number of electrons in the system can be arbitrarily fixed by tuning the chemical potential µ. For example, if the following inequalities are fulfilled,

or, equivalently,

then the system contains an integer number N of electrons (it is assumed that the N -electron fundamental gap E N g = I N 0 − A N 0 is positive, so that Eq. (73) can be fulfilled). In the special case where one of the inequality becomes a strict equality, say

which means that the chemical potential is exactly equal to minus the Nelectron ionization potential,

the N -and (N − 1)-electron solutions are degenerate (grand-canonical energy wise). Therefore, they can be mixed as follows,

where 0 ≤ α ≤ 1, thus allowing for a continuous variation of the electron number N (which now becomes fractional) from N − 1 to N :

This is the central idea in PPLB for describing the ionization of an N -electron system. Ionizing the (N + 1)-electron system gives access to the N -electron affinity. Interestingly, we recover from Eqs. (77) , (78) , and (79) the well-known piecewise linearity of the energy with respect to the electron number [3] :

In order to establish a clearer connection between the PPLB and N -centered formalisms, we follow the approach of Kraisler and Kronik [34] where the ensemble weight α is used as a variable, in place of the electron number N . Therefore, in PPLB, the ensemble energy reads as

where

is the exact ground-state ensemble density. Note that, if we introduce the ensemble density matrix operator

the ensemble energy and density can be expressed in a compact way as follows,

and

respectively.

On the basis of the "grand-canonical" ensemble formalism introduced in the previous section, we can extend the domain of definition of the universal Hohenberg-Kohn functional F [n] to densities n that integrate to fractional electron numbers, i.e.,

as follows,

where the ground-state density-functional wave functions fulfill the density constraint

From now on we will take α in the range

so that the integer electron number case systematically corresponds to α = 1. Therefore, in the present density-functional PPLB ensemble, the N -electron state will always contribute (even infinitesimally), and

At this point it is essential to realize that, unlike in N -centered eDFT, the ensemble weight α is not an independent variable. Indeed, according to Eqs. (87) and (91), it is an explicit functional of the density:

Therefore, in PPLB, the ensemble is fully determined from the density. The latter remains, like in regular DFT for integer electron numbers, the sole basic variable in the theory. Following Levy and Lieb [92, 115, 66] , the extended universal functional of Eq. (88) can be expressed in a compact way as follows,

where we minimize over grand-canonical ensemble density matrix operatorŝ

that fulfill the following density contraint:

The commonly used KS formulation of PPLB is recovered when introducing the non-interacting kinetic energy functional

Tr γ αT (96) and the in-principle-exact decomposition

where the Hxc functional now applies to fractional electron numbers. Let us stress that, unlike in N -centered eDFT, the Hxc functional has no ensemble weight dependence because the weight is determined from the density n. Any dependence in α is incorporated into the functional through the density. This is a major difference with N -centered eDFT where the ensemble weight and the density are independent variables, like in GOK-DFT. This subtle point will be central later on when comparing the two theories.

According to the variational principle, the exact ensemble energy can be determined, for a given and fixed value of α, as follows,

where we minimize over densities that integrate to the desired number N −1+α of electrons. According to Eqs. (96) and (97), the ensemble energy can be rewritten as

where the minimizing KS density matrix operator

reproduces the exact ensemble density of Eq. (83):

The orbitals from which Φ N −1,α

where, as readily seen from the following ensemble density expression,

the highest occupied molecular orbital (HOMO) [i.e., ϕ α N ] is fractionally occupied. This is the main difference with conventional DFT calculations for integer electron numbers.

Once the ensemble energy E α has been determined (variationally), we can evaluate the IP, which is the quantity we are interested in, by differentiation with respect to the ensemble weight α [see Eq. (81)], i.e.,

which, according to the Hellmann-Feynman theorem and Eqs. (99) , (100) , and (102), can be written more explicitly as follows,

thus leading to the famous Janak's theorem [116] :

As readily seen from Eq. (107), the energy ε α N of the KS HOMO does not vary with the fraction α > 0 of electron that is introduced into the (N − 1)electron system. Therefore, it matches the N -electron KS HOMO energy that we simply denote ε N N :

At this point it is important to mention that, unlike in N -centered eDFT [see Eq. (69)], there is no ensemble weight derivative of the Hxc functional involved in Janak's theorem. Such a quantity does not exist in PPLB, simply because the ensemble weight α and the density n cannot vary independently. However, while the number of electrons is artificially held constant in the N -centered formalism, it is not the case in PPLB. Indeed, variations in α induce a change in density [see the third contribution on the right-hand side of Eq. (105)] that does not integrate to zero:

Therefore, it is crucial, when evaluating the functional derivative of the Hxc energy δE Hxc [nγα KS ]/δn(r) (i.e., the Hxc potential), to consider variations of the density δn(r) that do not integrate to zero. This is unnecessary in Ncentered eDFT. In PPLB, however, the proper modeling of the xc potential is essential for describing charged excitations. This is clearly illustrated by the fact that the exact xc potential exhibits derivative discontinuities when crossing an integer electron number, as discussed further in Sec. 3.2.

Let us finally discuss the unicity of the xc potential. We recall that, in the present review, the external potential is simply the (Coulomb) nuclear potential of the molecule under study. It is fixed and it vanishes away from the system:

which we simply denote v ext (∞) = 0 in the following. As readily seen from Eq. (107), when describing a continuous variation of the electron number N in the range N − 1 < N < N , the KS potential becomes truly unique, not anymore up to a constant. This can be related to the unicity of the chemical potential which allows for fractional electron numbers, as discussed previously in the interacting case [see Eq. (75)]. As a result, the xc potential is truly unique. More precisely, as illustrated in Appendix A for a one-dimensional (1D) system, Janak's theorem implies that [117] 

According to Janak's theorem, the fundamental gap can be evaluated in PPLB, in principle exactly, from the HOMO energies as follows,

What is truly challenging in practice, in particular in solids [10] , is the extraction of this gap from a single N -electron calculation. Indeed, the HOMO energy ε N +1 N +1 of the (N + 1)-electron system has no reason to match the lowest unoccupied molecular orbital (LUMO) energy ε N N +1 of the N -electron system, simply because the infinitesimal addition of an electron to the latter system will affect the density [see Eq. (A.10)] and, consequently, the xc potential. The impact of an electron addition on the xc potential will be scrutinized in Sec. 3.2, in the context of N -centered eDFT. If we denote

the deviation in energy between the above-mentioned HOMO and LUMO, we recover the usual expression [106] 

where ∆ N xc can now be interpreted as the difference in gap between the physical and KS systems. As readily seen from the key Eq. (61) of N -centered eDFT, that we take in the regular N -electron ground-state DFT limit (i.e., ξ + = ξ − = 0), ∆ N xc is indeed a nonzero correction to the KS gap that can be expressed more explicitly as follows,

Note that we used in Eq. (115) the in-principle-exact decomposition

where the regular (weight-independent) Hartree functional is employed. The practical disadvantage of such a decomposition will be extensively discussed in Sec. 4. We focus here on the exact theory.

In the language of N -centered eDFT, ∆ N xc describes the variation in ensemble weights (while holding the ensemble density fixed and equal to the N -electron ground-state density n Ψ N 0 ) of the N -centered ensemble xc energy due to the infinitesimal removal/addition of an electron from/to the N -electron system. Evidently, standard (local or semi-local) DFAs do not incorporate such a weight dependence because they were not designed for N -centered eDFT calculations (we recall that the concept of N -centered ensemble has been proposed quite recently [83, 107, 113] ). Therefore, when such DFAs are used, the physical gap is systematically approximated by the (also approximate) KS one. Note that the resulting underestimation of the fundamental gap is highly problematic, for example, when computing transport properties [118, 119] . The interpretation that is given in PPLB for ∆ N xc is completely different. The latter actually originates from the discontinuity that the xc potential (which is the functional derivative of the xc energy) exhibits when crossing an integer electron number, hence the name derivative discontinuity. In the language of PPLB, ∆ N xc is not described at all when (semi-) local xc functionals are employed, simply because the latter do not incorporate functional derivative discontinuities. The connection between these two very different interpretations will be made in Sec. 3.2.

Let us finish the previous discussion with a detailed comment on the use of (orbital-dependent) exact exchange energies, which is often recommended for improving the description of fundamental gaps [10] . From the perspective of N -centered eDFT, using an exact exchange energy (or a fraction of it) is a way to incorporate weight dependencies into the ensemble exchange density functional. Indeed, according to Eq. (116), the exact exchange-only derivative discontinuity can be rewritten as

or, equivalently,

where we have introduced the double-weight N -centered ensemble KS density (119) which reduces to n Ψ N 0 when α = ξ = 0. Thus, we can remove all the contributions involving the derivative of the ensemble density [see the second line of Eq. (118) ] that are erroneously introduced by the first term on the right-hand side of Eq. (118)]. We recall that, in the evaluation of ∆ N x , we must differentiate with respect to the weight for a fixed density, as readily seen from Eq. (117) . At first sight, Eq. (118) is uselessly complicated, when compared with Eq. (117), but it will actually enable us to obtain simpler expressions. This trick has been introduced in the context of GOK-DFT [91, 103] . First we need to realize that, according to the exact expression of the N -centered ensemble Hx functional in Eq. (57),

Moreover, we have

where, in the α = ξ = 0 limit, the derivative of the density

can be evaluated from the regular N -electron KS frontier orbitals. By combining Eqs. (120)-(122) we finally obtain a simple (orbital-dependent) expression for the exchange-only derivative discontinuity:

Adding this correction to the exact KS gap leads to the following approximate fundamental gap expression,

where the physical interacting wave functions have been replaced by the KS ones, and correlation is introduced only through the correlation potential.

Note that a consistent implementation of Eq. (124) on the basis of Eq. (117) would in principle require using optimized effective potentials (OEPs) [120, 121] . Indeed, in the present formulation of N -centered eDFT, the KS orbitals are expected to be generated from a local (i.e., multiplicative) xc potential [see Eq. (53)], unlike in Hartree-Fock (HF)-based methods where the exchange potential is nonlocal. In practice, we would have to consider a trial local potential v and determine the corresponding KS orbitals,

thus ensuring that the KS wave functions in Eq. (120) are eigenfunctions of a non-interacting Hamiltonian, like in the exact theory. The ensemble energy would then be minimized with respect to v rather than the orbitals (hence the name OEP given to the method). Obviously, such a procedure induces a substantial increase in computational complexity, even though simplifications can be made in the optimization process [122] . Alternative variational evaluations of orbital-dependent ensemble exchange energies exist [120] . Their practical advantages and drawbacks will be discussed in detail in Sec. 4.

Crossing an integer electron number, which is a key concept in PPLB, can be described in the context of N -centered eDFT by considering so-called left and right N -centered ensembles [107] . These ensembles are recovered when ξ + = 0 (electron removal only) and ξ − = 0 (electron addition only), respectively. In the following, we will use the following shorthand notations,

right N -centered ensemble: (0, ξ + )

for convenience. For example, the right N -centered ensemble Hxc functional and KS orbital energies will be denoted as E 

and

respectively. Note that, with these notations, we have the following equivalence relation,

as readily seen from Eqs. (128) and (129). When Eq. (130) is fulfilled, the system is in its pure N -electron ground state which means, in the language of PPLB, that it contains exactly the integer number N of electrons.

A clearer connection between the two theories can be established by comparing the two limits ξ + → 0 + (which describes the infinitesimal addition of an electron to the N -electron system) and ξ + = 0 (or, equivalently, ξ − = 0). For that purpose, we first need to realize that, by analogy with PPLB (see Appendix A for the proof in the simpler 1D case), the exact IP/EA theorems of N -centered eDFT in Eqs. (67) and (68) can be alternatively written as follows [113] ,

and

where we recall that v 

and

We recall that the decomposition of Eq. (116) is employed for a direct comparison with PPLB. Let us now consider the ξ + → 0 + and ξ − = 0 limits in Eqs. (133) and (134), respectively. Since

it comes, by subtraction,

or, equivalently,

where we used Eq. (115) and the relation v

(r), according to Eq. (130) . Note that, as readily seen from Eq. (139) , ∆ N xc is insensitive to constant shifts in the xc potential, as expected from Eq. (114) .

The connection that is made explicit in Eq. (139) between the N -centered ensemble weight derivative ∆ N xc of the xc density-functional energy [see Eq. (115)] and the xc potential is an important result that was highlighted very recently in Ref. [113] . It proves that weight derivatives and derivative discontinuities are equivalent, thus extending to charged excitations what was already known for neutral excitations [105] . Indeed, if we systematically choose (but we do not have to in the N -centered formalism, unlike in PPLB) the xc potential that asymptotically goes to zero, i.e.,

then we recover what looks like a Janak's theorem [see Eqs. (131) and (132)] and, according to Eq. (138),

It then becomes clear that, in the region of the system under study (i.e., where the density n Ψ N 0 (r) is nonzero), the xc potentials obtained in the ξ + → 0 + and ξ + = 0 limits, respectively, cannot match. In order to fulfill the arbitrary constraint of Eq. (140), while still reproducing for ξ + > 0 the correct density in all regions of space [which includes a proper description of the density's asymptotic behavior (see Appendix A)], the xc potential must be shifted in the region the system, thus ensuring that the ground-state density n Ψ N 0 (r) is also correctly reproduced in that region. This has been nicely illustrated in Ref. [113] for an atom in 1D. Thus, we recover a well-known result of PPLB: When an electron is infinitesimally added (i.e., ξ + → 0 + ) to a system with an integer number of electrons (ξ + = 0 case), the xc potential exhibits a jump (in the region of the system) which, according to Eqs. (114) and (141), corresponds exactly to the deviation in fundamental gap of the true system from the KS one.

The fundamental gap expression of Eq. (61), which has been derived within the N -centered eDFT formalism, may intrigue PPLB practitioners. Indeed, it makes it possible to describe charged excitations, in principle exactly, without invoking explicitly the concept of derivative discontinuity. Instead, all our attention should be focused on the weight dependence of the N -centered ensemble xc density functional. Note that, despite this major difference between N -centered eDFT and PPLB, the xc potential exhibits derivative discontinuities in both theories, as we have seen in Sec. 3.2. One may argue that modeling weight dependencies in ensemble xc density functionals is actually easier than modeling functional derivative discontinuities. Nevertheless, as discussed in further detail in Secs. 4 and 5, designing weight-dependent exchange and correlation DFAs from first principles raises several fundamental questions to which, up to now, no definitive answers have been given.

From a conceptual point of view, the fact that we do not need anymore to put efforts into the explicit description of derivative discontinuities, once we have moved from the standard PPLB picture to the N -centered one, can be interpreted as follows. Unlike in PPLB, the constraint in Eq. (140) is arbitrary because the KS potential remains unique up to a constant when charged excitations occur in N -centered eDFT, by construction. If, for simplicity, we keep this constraint for ξ + = 0, i.e., we setṽ ξ+=0 xc

(∞) = 0, which is likely to be fulfilled in a practical N -electron DFT calculation, it can be relaxed as follows, when ξ

thus leading toṽ 

In other words, via the shifting procedure of Eq. (142), we can simply move the derivative discontinuity away from the system, i.e., in regions where the density is essentially equal to zero. Consequently, with this change of paradigm, the absence of derivative discontinuities in standard semi-local DFAs should not be considered as an issue anymore. The ability of the local density approximation (LDA) to reproduce relatively accurate LZ-shifted KS orbital energies, as shown in a 1D atomic model [113] , is actually encouraging since the latter are central in the evaluation of both the IP and the EA [see Eqs. (69) and (70)]. On the other hand, the resulting charged excitation energies are rather poor because weight dependencies are completely absent from standard LDA [113] . We hope that, in the near future, (much) more efforts will be put into the design of weight-dependent DFAs. Recent developments based on uniform electron gas models [91, 79] are a first and important step in this direction.

We have shown in Sec. 3 that the infamous derivative discontinuity problem, which appears in DFT when describing electronic excitations, can be bypassed, in principle exactly, via a proper modeling of the ensemble weight dependence in the xc density functional. We focus in this section on the design of weight-dependent exchange DFAs. Despite several (scarce though) attempts [100, 123, 124, 79] , it is still unclear how weight dependencies can be introduced into standard (semi-) local exchange functionals in a general and systematically improvable way. The use of orbital-dependent exchange functionals seems much more promising in this respect [71, 91, 120] .

Combining (a fraction of) orbital-dependent HF-like exchange energies with semi-local DFAs has been a key ingredient in the success of regular ground-state DFT in chemistry. This procedure finds its rigorous foundation in the generalized KS theory of Seidl et al. [125] . As we will see in the following, its extension to ensembles is nontrivial because different formulations that have pros and cons are possible. The resulting dilemma is nicely summarized by the title "Ensemble generalized Kohn-Sham theory: The good, the bad, and the ugly" of a recent paper by Gould and Kronik [120] . Their discussion of the current situation will serve as a guideline for this section. Following Lieb [66] , we will show how an in-principle-exact (OEP-free) hybrid eDFT approach can be derived simply by exploiting the concavity (in potential) of the state-averaged HF energy. For simplicity, we will focus on GOK ensembles but the discussion applies to other eDFTs like, for example, PPLB or N -centered eDFT (see Sec. 3).

The reason why extending generalized KS-DFT [125] to ensembles leads to a dilemma has actually nothing to do with DFT. It is more a wave function theory problem that arises at the HF level of approximation. Therefore, for clarity, we will first discuss the extension of HF theory to ensembles.

We start with a brief review of the procedure that is usually followed by DFT practitioners for extending HF to (GOK in the present case) ensembles. In the regular scheme, the ensemble HF energy is evaluated variationally by inserting the ensemble (spin-summed one-electron reduced) density matrix (eDM) into the ground-state DM-functional HF energy. For that reason, we refer to the approach as eDMHF. The corresponding potential-functional ensemble energy can be expressed as follows,

whereV = dr v(r)n(r) is a local potential operator (in practice it would correspond to the nuclear potential). The eDM is evaluated from the trial orthonormal set {Φ I } of single-configuration wave functions (i.e., Slater determinants or configuration state functions). In an arbitrary orthonormal orbital basis {ϕ p }, the eDM reads in second quantization as

In this context, the ensemble Hx energy is evaluated as follows,

where the ground-state HF interaction functional reads as

and

As shown in Appendix B, the orbitals ϕ w p , from which the minimizing wave functions Φ w I in Eq. (144) are constructed, fulfill the following stationarity condition:

The (possibly fractional) occupation number θ w p of the orbital ϕ w p within the ensemble is determined from the ensemble weights and the (fixed) integer occupation numbers n I p of ϕ w p in each Φ w I as follows,

The ensemble Fock matrix elements f w rs ≡ f rs (D w ) in Eq. (149) are functionals of the eDM D w ≡ {D w nl = δ nl θ w l }:

where h rs ≡ ϕ w r |ĥ|ϕ w s [withĥ ≡ − 1 2 ∇ 2 r + v(r)] and rn|sl ≡ dr dr ϕ w r (r)ϕ w n (r )ϕ w s (r)ϕ w l (r )/|r − r |

are regular one-and two-electron integrals, respectively.

By analogy with the complete active space SCF (CASSCF) method [126] , we can distinguish the doubly occupied (so-called inactive) ϕ w i , ϕ w j orbitals from the partially occupied (so-called active) ϕ w u , ϕ w v and unoccupied (virtual) ϕ w a , ϕ w b ones. Thus, the stationarity condition of Eq. (149) can be detailed as follows:

and

Since θ w i = θ w j = 2 and θ w a = θ w b = 0, there is no specific condition for the inactive-inactive and virtual-virtual blocks of the Fock matrix, like in a regular ground-state HF calculation. Therefore, we can freely rotate the orbitals within the inactive and virtual orbital subspaces. As readily seen from Eq. (154), this statement holds also for active orbital subspaces in which the orbitals have the same fractional occupation (θ w u = θ w v ). However, if θ w u = θ w v , then f w uv = 0. In conclusion, an optimal set of orbitals can be determined by diagonalizing the ensemble Fock matrix, i.e., by solving the (self-consistent) eigenvalue equation

The fact that standard SCF routines can be trivially recycled in this context is the main reason why ensemble HF and, more generally, hybrid eDFT calculations are performed this way. However, as discussed further in the following, the eDMHF energy is unphysical in many ways. For example, by construction, it varies quadratically with the ensemble weights [see Eqs. (146) and (147)] while the true physical ensemble energy is expected to vary linearly.

The most severe issue with the eDMHF energy expression of Eq. (144) is that it incorporates unphysical interactions between the states of the ensemble. These are known as ghost interactions (GIs) [127] . The error, which is inherent to the eDMHF method, simply originates from the fact that, at the HF level of approximation, the ground-state interaction energy is a quadratic functional of the density matrix. More explicitly, we have

where, as readily seen, GI terms arise from all "I = J" pairs. Even though GI corrections can be applied on top of the converged eDMHF energies [91, 128] , the procedure is not variational, thus making the evaluation of energy derivatives (and therefore, according to Eq. (9), of excited-state properties) less straightforward. Let us stress that, in the original formulation of GOK-DFT [82] , the ensemble Hartree energy, which is evaluated from the standard (ground-state) Hartree functional [see Eq. (20) ], includes GI errors [see the first term on the right-hand side of Eq. (156) ]. In the exact theory, the latter are supposed to be removed by the weight-dependent ensemble exchange functional. It is not necessarily the case in practice when, for example, standard (weight-independent) local or semi-local DFAs are employed [128, 97] .

For wave function theory practitioners, using eDMHF with ad hoc GI corrections would probably seem uselessly complicated. Indeed, substituting the (GI-free) weighted sum of individual (single-configuration) interaction energies, which are evaluated from the individual density matrices, for the HF interaction eDM functional of Eq. (144) looks, at least at first sight, like a simple and straightforward solution to the problem:

Note that, in Eq. (157), we assumed that individual interaction energies can be written as functionals of the individual (one-electron reduced) density matrices, for simplicity. The variational evaluation of state-averaged interaction energies is discussed in detail in Sec. 4.1.3 on that basis. Such a simplification is always valid for single Slater determinants. For more general multideterminant (single configuration though) wave functions, the simplification in Eq. (157) might be used as an (additional) approximation. Alternatively, one may evaluate exactly multideterminant interaction energies from the individual two-electron reduced density matrices, by analogy with the state-averaged CASSCF (SA-CASSCF) method [126] . The latter approach is not described further in the present review. Both (approximate anyway) strategies lead ultimately to an in-principle-exact eDFT, once a proper complementary correlation ensemble density functional has been introduced (see Sec. 4.4).

As discussed in the next section, the reason why the state-averaging of interaction energies is not as popular as one would expect is that, as we switch from eDMHF to the state-averaged energy paradigm, the orbital optimization cannot be performed anymore with standard SCF routines. An additional implementation work is needed in this case [71] . We stress that this statement holds even when individual one-electron reduced density matrix-functional interaction energies are employed [see Eq. (157)], as assumed in the rest of the review.

The paradigm on the right-hand side of Eq. (157) can be seen as an adaptation of the SA-CASSCF method [126] to single-configuration wave functions. While, in SA-CASSCF, (correlated) multiconfigurational wave functions are employed, we simply restrict, in the present case, the energy minimization to sets of (uncorrelated) single-configuration wave functions. The resulting (GIfree) total ensemble energy, which is obtained from the eDMHF energy expression and the substitution in Eq. (157), will be referred to as state-averaged HF (SAHF) energy in the following. It reads as follows,

where, as already mentioned after Eq. (157), we assume for simplicity that interaction energies can be evaluated from the individual (one-electron reduced) density matrices. Let us denote Φ w I (with a tilde symbol) the minimizing single-configuration wave functions so that they can be clearly distinguished from the eDMHF ones. As further discussed in Appendix C, these minimizing SAHF wave functions are all constructed from the same set of orthonormal molecular orbitals which can be optimized variationally through orbital rotations. The minimizing orbitals φ w p fulfill the following stationarity conditions [see Appendix C],

where the individual (non-local) density-matrix-functional Hx operators read asv 

so that the stationarity condition of Eq. (159) can be rewritten as follows,

Let us stress that, unlike in GOK-DFT (see Sec. 2.1.2) or eDMHF, the minimizing orbitals are a priori not eigenfunctions of their associated Fock operator. Indeed, if we consider two fractionally occupied orbitalsφ w u andφ w v , and assume that φ w u |F w u |φ w v = 0, Eq. (164) would immediately imply that φ w u |F w v |φ w v = 0, which is unlikely due to the orbital dependence of the Fock operator [see Eq. (163)]. As a result, the self-consistent SAHF equations will have the following general structure, where off-diagonal one-electron energy couplings ε w qp p =q cannot be removed. In practice, Eq. (165) can be solved with the coupling operator technique [129, 75] which consists in diagonalizing repeatedly (θ w pF w p −θ w qF w q )/(θ w p −θ w q ) until convergence is reached. In the latter case, the matrixε w ≡ ε w qp becomes hermitian, as a consequence of Eqs. (164) and (165), and the hermiticity of the Fock operators:

Let us summarize what we have learned from the previous subsections. While the orbital optimization in eDMHF is relatively straightforward, because standard SCF routines can be recycled in this context, it is more involved in SAHF because no to-be-diagonalized Fock operator emerges from the stationarity condition. On the other hand, the commonly used eDMHF energy expression suffers from severe GI errors, while SAHF is completely GI-free. The latter point is probably the strongest argument for promoting SAHF over eDMHF. Note that both schemes would be good starting points for turning the recently formulated ensemble reduced density matrix functional theory (w-RDMFT) [130] into a practical method for the computation of low-lying excited states. In the following, we will show how eDMHF and SAHF can be merged rigorously with eDFT.

In order to derive in-principle-exact hybrid eDFT schemes, where (a fraction of) orbital-dependent exchange energies are combined with ensemble density functionals, we need to "exactify" the eDMHF and SAHF approximations reviewed previously. In a DFT perspective, an approximation becomes exact when it reproduces the exact density of the system under study. This is how KS-DFT transforms an approximate non-interacting problem into an exact one. In the present case, we want to extend the Hohenberg-Kohn theorem to the more advanced eDMHF and SAHF approximations. For that purpose, convex analysis [66, 131, 132] turns out to be a powerful mathematical tool because, as we will see, it allows for the derivation of several exact eDFTs within the same (unified) formalism.

For the sake of generality, we will express the various approximate ensemble energies discussed previously as follows,

where κ denotes the collection of variational parameters (in the present case, the latter will be orbital rotation parameters). Each approximation (noninteracting (KS), eDMHF, or SAHF) is characterized by a specific potentialindependent function of κ:

where single-configuration ground-and excited-state wave functions {Φ I (κ)} are employed. The potential-dependent contribution to the energy expression of Eq. (168) is determined from the ensemble density n w (κ, r) = I w I n Φ I (κ) (r).

From a mathematical point of view, the fact that the approximate ensemble energies are evaluated variationally has important implications. Even though they are approximate, these energies still share a fundamental property with the exact ensemble energy, namely the concavity with respect to the local potential v. Indeed, for two potentials v a and v b , α in the range 0 ≤ α ≤ 1, and any set κ of variational parameters, we have

thus leading to

Following Lieb [66] , we can now construct (thanks to this concavity property) an approximation to the universal GOK density functional as follows,

where we assume, for simplicity, that a maximum is reached. An even more rigorous definition (from a mathematical point of view [66] ) would actually be obtained by using a "sup" instead of a "max" [and a "inf" instead of a "min", in Eq. (168)]. The maximizing potential v w approx.

[n] in Eq. (174) fulfills the following stationarity condition:

= n(r).

We conclude from the Hellmann-Feynman theorem that, when the approximate ensemble energy of Eq. (168) is calculated with v = v w approx.

[n], the minimizing single-configuration wave functions (which are determined from the minimizing κ) reproduce the desired ensemble density n. Thus, we automatically extend the Hohenberg-Kohn theorem to eDMHF and SAHF ensembles. If we choose for n the true physical ensemble density of a given system, both approximations become exact density wise, because they reproduce the correct density. Exact ensemble energies can then be recovered from the approximate ensembles once a complementary ensemble (x)c density functional has been introduced. This final step will be discussed in Sec. 4.4 Let us finally focus on the OEP-and GI-free SAHF approximation. The corresponding universal density functional reads more explicitly as [see Eqs. (168), (171), and (174)]

where v w SAHF [n] denotes the stationary (maximizing) density-functional potential of Eq. (175) in the particular case of the SAHF approximation. By analogy with the constrained-search formalism of Levy [92] , the SAHF functional can be rewritten as follows,

where the density constraint {Φ I } w → n imposed on the single-configuration wave functions {Φ I } reads as I w I n Φ I (r) = n(r). Note that, as illustrated in Sec. 4.3, some densities may not be v-representable by a single SAHF ensemble. In this case, a more general Levy-Lieb-like [66] functional, where an ensemble of ensembles is considered, should be employed:

where the density constraint reads as 

We stress that, in the above more general definition, the additional ensemble weights α (i) are not given, unlike the GOK ensemble weights w. They are determined from the density constraint of Eq. (179). An interesting feature of the constrained-search formalism is that it applies to densities that might be ensemble N -representable but not SAHF v-representable, i.e., densities that cannot be generated from a given potential v according to Eq. (168). For clarity, we will use in Sec. 4.4 the simpler density functional expression of Eq. (177) [rather than the one in Eq. (178)].

In order to compare SAHF with eDMHF, both methods are applied in this section to the (two-electron) Hubbard dimer model [84, 85, 96] . Despite its simplicity, it is nontrivial and has become in recent years the model of choice for analyzing and understanding failures of DFT or TD-DFT, but also for exploring new concepts [133, 134, 135, 99, 101, 136] . In this model, the ab initio Hamiltonian is simplified as follows,

where operators are written in second quantization andn i = τ =↑↓n iτ is the density operator on site i (i = 0, 1). Note that the local potential reduces to a single number ∆v which controls the asymmetry of the dimer. The density also reduces to a single number n = n 0 , which is the occupation of site 0, given that n 1 = 2 − n.

We consider in the following a two-electron singlet biensemble consisting of the ground state and the singly-excited state.

In (restricted) SAHF theory, both ground-and excited-state wave functions are approximated by configuration state functions. Following this view, trial singlet ground-state Φ 0 and first excited-state Φ 1 wave functions of the twoelectron Hubbard dimer read as follows, in second quantization,

and

respectively, where σ 0 and σ 1 stand for the molecular orbitals (MOs) written in the basis of the orthonormal local atomic orbitals (AOs) a and b. In this context, the latter simply correspond to site 0 (ĉ † aτ ≡ĉ † 0τ ) and 1 (ĉ † bτ ≡ĉ † 1τ ), respectively. Bonding σ 0 and anti-bonding σ 1 MOs can be determined through orbital rotation as follows,

where the angle α is the sole variational parameter in the model. In SAHF, the energy that must be minimized, for a fixed ensemble weight w in the range 0 ≤ w ≤ 1/2, is constructed as follows,

where

and

We use standard notations for Coulomb and exchange two-electron integrals,

both expressed in the MO basis {σ 0 , σ 1 } of Eq. (184). At this point, let us stress that the SAHF energy substantially differs from that of a (truncated) configuration interaction calculation. Indeed, in the present case, the weight w is fixed. Moreover, the configurations Φ 0 and Φ 1 are never coupled explicitly, whether the dimer is symmetric or not. They are just mixed through the ensemble formalism. Note that, in the symmetric ∆v = 0 case, symmetry can in principle be artificially broken, like in (spin) unrestricted calculations. This feature will be discussed in further detail in the next section. For convenience, we denote

so that the symmetric solution corresponds to θ = 0. Consequently, the SAHF energy expression of Eq. (185) reduces to

As readily seen from the above expression, the SAHF energy is π-periodic, which means that it is sufficient to vary θ in the range − π 2 ≤ θ ≤ π 2 .

Similarly, from the general expression in Eq. (144), we obtain the following analytical expression for the eDMHF energy:

Note that, according to Eq. (184), the ensemble density (on site 0) varies with the trial angle θ as follows,

(193)

Let us concentrate on the symmetric dimer, for which ∆v = 0. In this case, the dimer is a prototype for the H 2 molecule. If we denote

then the SAHF energy reads as

where

By taking its first derivative with respect to θ, we obtain the following stationarity condition,

Therefore, θ = 0 is systematically an extremum where the traditional in-phase and out-of-phase linear combinations for σ 0 and σ 1 are recovered. The nature (maximum or minimum) of this stationary point, which is discussed in the following, emerges from a straightforward evaluation of the energy curvature:

Remembering that |ρ(θ)| ≤ 1, the stationarity condition of Eq. (197) is also fulfilled for two additional (opposite) θ values given by

as long as

Note that, with this notation, the successive energy derivatives can be expressed as follows,

and

The symmetric solution θ = 0 will not be the absolute minimum anymore when the above curvature becomes strictly negative, which implies 1/ρ 0 > 1, as readily seen from Eq. (202). Obviously, this constraint can only be fulfilled if w > 1/3, since 1/ρ 0 must be strictly positive. This is a necessary but not sufficient condition. More precisely, for weights in the range 1/3 < w ≤ 1/2, electron correlation should be strong enough such that ρ 0 < 1 or, equivalently,

Interestingly, if we introduce effective weight-dependent hoppingt = t(1 − w) and on-site interactionŨ = U (3w − 1) parameters, the condition in Eq. (203) can be rewritten as

which resembles the usual definition of moderate (up to strong) electron correlation in lattices. Actually, in the commonly used equiensemble case (w = 1/2), the effective ratio matches the physical one 2t/U .

Finally, when w ≤ 1/3, the SAHF energy becomes convex at θ = 0 and, since it can have two additional (say θ + > 0 and θ − = −θ + ) stationary points at most in the range −π/2 ≤ θ ≤ π/2 [see Eqs. (194) and (201)], the symmetric solution has to be the absolute minimum. The different possible scenarios are summarized in Table 1 . Note that a connection can be made with the (spin) 

unrestricted energy of a symmetric dimer (e.g., the H 2 molecule in the minimal basis) [137] . For sufficiently large bond distances, the restricted solution becomes a saddle point and two unrestricted lower-in-energy solutions emerge. Accordingly, the constraint in Eq. (203) is compatible with a reduction of the t value featuring an increasing bond length. Still, even in the strictly correlated t → 0 limit, θ = 0 remains the global minimum for weights in the range 0 ≤ w ≤ 1/3.

In order to illustate the above discussion, trial SAHF energies are plotted as functions of the rotation angle θ in the top panel of Fig. 1 for the strongly correlated U/t = 3.5 dimer and various biensemble weight values. As expected, the symmetric θ = 0 solution gives systematically the lowest ensemble energy as long as w ≤ 1/3. When w > 1/3, two scenarios can be observed. For example, when w = 0.35, which gives 2(1 − w)/(3w − 1) = 26 U/t, the dimer is not strongly correlated enough to break the symmetry and θ = 0 is still the global minimum. However, for the larger w = 0.475 weight value [2(1 − w)/(3w − 1) = 2.47 < U/t in this case] or in the commonly used equiensemble case [w = 0.5 and 2(1 − w)/(3w − 1) = 2 < U/t], the energy has the expected double-well shape, thus leading to two degenerate (non-zero) minima, both corresponding to asymmetric solutions. Interestingly, in such situations, the popular eDMHF approach always favors the symmetric solution, as shown in the bottom panel of Fig. 1. 

We show in Fig. 2 the potential-ensemble-density maps ∆v → n w (θ min (∆v)) (205) generated from Eq. (193) and

at both eDMHF and SAHF levels of approximation for the fixed U/t = 3.5 interaction strength value. Various scenarios are illustrated, in particular those where broken symmetry SAHF solutions are obtained when the dimer is symmetric (see Sec. 4.3.2). As we will see, what happens in the symmetric case can play a crucial role in the density-functional description of the asymmetric dimer. In cases where w ≤ 1/3, or w > 1/3 and 2(1 − w)/(3w − 1) > U/t, both approximations give smooth density profiles. We note in passing the one-to-one correspondence between potentials and ensemble densities, as expected from the concavity of the eDMHF and SAHF energies (see Sec. 4.2). However, when w > 1/3 and 2(1 − w)/(3w − 1) < U/t, the SAHF density profile exhibits a discontinuity at ∆v = 0, unlike the eDMHF one. This step in density can be interpreted as follows. If w > 1/3 and the constraint of Eq. (203) is fulfilled, as ∆v → 0 ± , we will recover the SAHF biensemble so-

and θ ± are the minimizing angles associated to the broken-symmetry orbitals. Any slight deviation from ∆v = 0 will favor one of these solutions, depending on its sign, as illustrated in the top panel of Fig. 1 (see the "∆v = +0.15" curve which exhibits, in the equiensemble case, a single absolute minimum in the vicinity of θ + ). Note that, in eDMHF, the minimizing angle simply passes through θ = 0 when the potential changes from ∆v = 0 − to ∆v = 0 + [see the "∆v = +0.15" curve in the bottom panel of Fig. 1 ], and no discontinuity is observed in the density profile. The step in density observed in SAHF covers the density range n − ≤ n ≤ n + , where

In the equiensemble case (w = 0.5), we have n± = 1.0 ± 0.410326, as readily seen from the top panel of Fig. 2 . We keep many digits for analysis purposes (see the bottom panel of Fig. 3 ). It is important to stress that none of the densities in the range n − < n < n + , which includes the a priori simple symmetric n = 1 case, can be represented by a single SAHF ensemble. This severe v-representability issue, which would deserve further investigation at the ab initio level, for example, in the stretched H 2 molecule, might be used as an argument for promoting the eDMHF approach over the SAHF one in practical calculations. A counter-argument is of course the presence of GI errors in eDMHF. At the formal level, the v-representability issue can be solved easily as follows. As illustrated in Fig. 3 , Lieb's maximization of Eq. (174) systematically returns ∆v = 0 for input densities in the range n − ≤ n ≤ n + . As soon as the input density leaves this interval, a non-zero maximizing potential is obtained [see the bottom panel of Fig. 3 ]. If we exploit the strict degeneracy of the broken-symmetry solutions, we can write, for ∆v = 0 [see Eq. (181)],

where 0 ≤ α ≤ 1 andγ

is the convex combination of the two degenerate SAHF ensemble density matrix operators. The density (on site 0) of the resulting "ensemble of ensembles" reads as

and, as readily seen, it can vary continuously from n − to n + . The generalization of this approach to the ab initio theory is provided in Eqs. (178) former misses all correlation effects. These effects can actually be introduced into the theory as a density-functional complement:

As a result, according to Eqs. (10) and (177), the exact ensemble energy can be calculated variationally as follows,

or, equivalently,

thus leading to the final expression:

By differentiating the density-functional correlation energy with respect to any (orbital rotation) variational parameter κ pq [see Appendix C],

we realize that the minimizing orbitals in Eq. (214), from which the singleconfiguration wave functions that reproduce the exact ensemble density n w are constructed (we denote themΦ w I for convenience), fulfill SAHF-like selfconsistent equations [see Eq. (166)] where the density-functional correlation potential δẼ w c [n]/δn(r)

is simply added to the local external one.

In practice, the quantities of interest are usually the excitation energies and, more generally, the ground-and excited-state energy levels (which are needed, for example, for geometry optimizations). According to Eqs. (9) and (214), and the Hellmann-Feynman theorem, the latter can be evaluated, in principle exactly, as follows,

where the following relation has been used,

andṼ w c ≡Ṽ w c [n w ] is the LZ-shifted [104, 91] ensemble correlation-only densityfunctional potential:

δn(r) n(r) dr n(r)

.

Interestingly, an expression similar to that of Eq. (216) has been derived in the context of GOK-DFT [91, 103] (see also Sec. 5) where, unlike in the present case, a local ensemble exchange potential was used.

By analogy with conventional (ground-state) hybrid functionals, we can modify as follows the SAHF-based functional of Eq. (176), in order to combine rigorously a fraction λ of SAHF exchange energy with a weighted sum of (approximate) individual density-functional exchange energies, thus leading to another (Hx-only) approximation to the GOK functional:

where

and E Hx [n] is the regular ground-state Hx density functional of KS-DFT. The above approximate ensemble energy can be computed from SAHF routines simply by scaling the individual non-local exchange potentials and then adding to each of them the complementary fraction (1 − λ) of individual local densityfunctional exchange potential. As readily seen from Eq. (220), in the present scheme, the total ensemble Hx energy remains GI-free, even in the limiting λ = 0 case. On that basis, an exact hybrid eDFT can be derived along the lines of Sec. 4.4 by considering the following in-principle-exact decomposition of the universal GOK functional (where correlation is described with an ensemble density functional):

A density-functional correction ∆Ẽ w,λ x [n] to the ensemble exchange energy must in principle be introduced since each individual exchange energy E x [n Φ I ] is evaluated, for both ground-and excited-state densities, as a ground-state exchange energy. This correction is usually neglected in practical calculations [75] . This observation would actually hold also for approximate ensemble correlation energies that are constructed from the regular ground-state correlation functional E c [n] of KS-DFT (see below and Sec. 5.1). Note finally the λ-dependence of both ∆Ẽ w,λ x [n] andẼ w,λ c [n] functionals in Eq. (221). It originates from the fact that the ensemble xc energy is now evaluated from the (SAHF-like) λ-dependent single-configuration wave functions that reproduce the desired density n.

Finally, as discussed in further detail in Sec. 5, it is quite common to construct weight-dependent ensemble density-functional correlation energies by recycling the ground-state correlation functional as follows [71] ,

This standard approximation, which is referred to as ground-state individual correlations (GS-ic) scheme in the following, can also be made formally exact within the present formalism. Indeed, once we have introduced the following (approximate) potential-functional energy

and the subsequent density functional

we only need to consider the alternative (but still exact) partitioning of the GOK functional

where complementary density-functional corrections to both exchange and (GS-ic) correlation energies have been introduced. Note that, in practice, these corrections are simply neglected [75] . While SAHF is expected to provide, through its orbital dependence, a proper description of the ensemble exchange energy, modeling the correlation density-functional correction ∆Ě w,λ c [n] remains a necessary and challenging task that has attracted too little attention until now. For that purpose, it is essential to have a deeper understanding of how individual correlation energies are connected to the ensemble one. This is the main focus of the next section.

While the previous section was dedicated to the description of orbital-and weight-dependent ensemble Hx energies, this last section deals with correlation effects in many-body ensembles. For convenience, we continue focusing on GOK ensembles but the discussion applies to other types of ensembles like, for example, N -centered [83] or thermal ones [78, 76, 77, 136] . We will work within the original GOK-DFT formalism [82] , where a local multiplicative ensembledensity-functional Hxc potential is employed, but the discussion holds also when orbital-dependent exchange energies are employed (see Sec. 4).

To the best of our knowledge, very few works have addressed the construction of weight-dependent ensemble correlation DFAs from first principles. We can essentially distinguish two different general strategies. In the first and most straightforward one, which was introduced in Eq. (222) and that we referred to as GS-ic, the (weight-independent) ground-state correlation functional is recycled as follows,

where, in the exact theory, the KS wave functions {Φ w I } are expected to reproduce the true ensemble density n w .

More recently [91, 79] , Loos and coworkers explored another path. They designed a first generation of weight-dependent ensemble LDA (eLDA) correlation functionals where the regular ground-state LDA functional E LDA c [n] = dr n(r) c (n(r)), which is based on the infinite uniform electron gas (UEG) model, is combined with the density-functional correlation excitation energies of a finite UEG (hence the acronym f LDA used below) as follows,

The individual correlation functional E f LDA c,I

[n] = dr n(r) f c,I (n(r)) is constructed from the Ith state of the finite UEG:

where N f is the number of electrons in the finite gas. If N f is fixed, a parameterization of the correlation energy per particle f c,I (n) as a function of the uniform density n is obtained by varying the volume of the gas [91] . Note that, in a uniform system, there is no need to introduce weight dependencies into the interacting and non-interacting (ground-or excited-state) densityfunctional wave functions, unlike in the general definition of Eq. (25) . Indeed, all the eigenstates of the (interacting or non-interacting) uniform gas have the same (uniform) density n, which then becomes the ensemble density of the gas, whatever the value of the ensemble weights:

Note also that, while the finite UEG allows for the incorporation of weight dependencies into the correlation functional, the use of a regular LDA correlation functional reduces finite-size errors. Refinements are possible, for example, by including a dependence in the Fermi hole curvature [138] .

The strategies depicted in Eqs. (226) and (227) miss various correlation effects that we briefly review below. More insight will be given in the next subsections. Let us start with the GS-ic approximation. From the exact expression,

where Ψ 0 [n] and Φ 0 [n] are the interacting and KS non-interacting ground-state density-functional wave functions of regular KS-DFT, respectively, we immediately identify two sources of errors. The first one is related to the fact that, as already mentioned in Sec. 2.1.3, the individual KS density n Φ w I does not necessarily match the interacting individual one n Ψ I . This subtle point was recently highlighted by Gould and Pittalis [86, 139] . It induces what the authors referred to as density-driven (DD) correlation effects. Even if the true individual densities {n Ψ I } (which can be extracted in principle exactly from the KS ensemble, as shown in Eq. (27) and Ref. [103] ) were inserted into the expression of Eq. (230), we would still not recover the correct individual excited-state correlation energies simply because Ψ 0 [n Ψ I ] will always be a ground-state wave function, even when n Ψ I is an excited-state density. The missing energy contribution is connected to the concept of state-driven (SD) correlation [86] . Interestingly, eLDA describes (approximately) SD correlations, as readily seen from Eq. (228). However, it completely misses DD ones, simply because KS and interacting (ground-or excited-state) wave functions have the same density in a uniform system.

Even though the physical meaning of DD and SD correlations is rather clear, it is less obvious how their contributions to the total exact ensemble correlation energy should be defined mathematically [86, 103, 139, 140] . Addressing this fundamental question is of primary importance for the design of more accurate and systematically improvable ensemble correlation DFAs, which is probably the most challenging task in GOK-DFT. Up to now, we have discussed the concept of DD and SD correlations in the light of the GSic approximation [see Eqs. (226) and (230)]. We may actually wonder if a proper definition can be (or should be) given without referring explicitly to the GS correlation functional of KS-DFT. Indeed, the latter appears naturally in GOK-DFT only in the limiting w = 0 case. Gould and Pittalis [86] , and then Fromager [103] , recently addressed this SD/DD ensemble correlation energy decomposition issue from that perspective. A detailed and complemented review of the two approaches is presented in the following.

Before proceeding with the extraction of individual correlation energies from the GOK-DFT ensemble energy, which is convenient for deriving in-principleexact SD/DD decompositions [103] , we would like to highlight the importance of weight dependencies in the KS wave functions. It might be surprizing at first sight because the true ground and excited states of the system under study are of course weight-independent. We explain below, with a simple argument, why it cannot be the case in the KS ensemble.

Since the KS and true ensemble densities match for any set of weights w, their derivatives with respect to the weights also match. Therefore, if we consider the ground-state w = 0 limit of GOK-DFT, it comes

or, equivalently,

where {Φ I } denote here the ground-and excited-state KS wave functions generated from a regular ground-state KS-DFT calculation. In KS-DFT, the density constraint applies to the ground state only, i.e., n Φ0 (r) = n Ψ0 (r), not to the excited states. Thus, we obtain the exact individual excited-state density expression, which can be recovered from Eq. (27) when w = 0,

where we readily see that, in GOK-DFT, the KS wave functions (the groundstate one in the present case) are necessarily weight-dependent. This feature is central in the design of DD correlation energies [103] , as discussed further in the following. We refer to Eq. (275) for an illustrative example (based on the prototypical Hubbard dimer) of weight-dependent KS ground-state density.

In this section we revisit the derivation of the individual energy levels in Eq. (36) in order to construct individual correlation energies within the ensemble under study. For that purpose, we start from the exact relation between individual and ensemble energies in Eq. (9) , and the variational GOK-DFT ensemble energy expression of Eq. (15), thus leading to, according to the Hellmann-Feynman theorem,

or, equivalently,

where, in analogy with Eq. (119), the following double-weight ensemble KS density has been introduced:

The last contribution (that is subtracted) on the right-hand side of Eq. (235) originates from the Hellmann-Feynman theorem. In other words, derivatives of the KS wave functions (and, therefore, of their densities) do not contribute to the derivatives of the total ensemble energy, because the latter is variational.

As shown in Refs. [91, 103] , the Hx contribution to the individual Jth energy level reduces to the expectation value of the two-electron repulsion operator evaluated for the Jth KS state, as one would guess. Indeed, once we have realized that, for given weight values ξ, the ensemble KS potential that reproduces n ξ,w is simply the one that reproduces the true ensemble density n w , we deduce from Eq. (24) that

and, consequently,

simply because the KS density n Φ w J does not match, in general, the true physical density n Ψ J . The concept of DD correlation, which was introduced recently by Gould and Pittalis [86] , originates from this observation. The important property that the true individual correlation energies share with the individual correlation components is that both of them can be used to construct the total ensemble correlation energy, i.e.,

The above expression can be deduced from Eq. (242) and the fact that, for any

From now on we will substitute the decomposition of Eq. (246) for the more conventional one of Eq. (243). As shown in the following, with this change of paradigm, DD-type correlation energy contributions will naturally emerge from the derivation of a more explicit expression. Unlike in Ref. [86] , the approach of Ref. [103] , which is reviewed in the next section, does not require additional (state-specific) KS systems, thus avoiding formal issues such as the non-uniqueness of KS potentials for individual excited states or vrepresentability issues [103] .

Let us now derive a more explicit expression for the ensemble correlation energy, on the basis of Eq. (246). We start with a simplification of the true individual correlation energy expression of Eq. (242), where the standard decomposition into components [see Eq. (243)] of the ensemble correlation energy will be employed. On the one hand, we will have

where, according to Eq. (244), 

On the other hand, according to Eq. (236),

thus leading to [see Eq. (27)]

Finally, by combining Eqs. (242), (250), and (252), we recover the expression of Ref. [103] for the deviation of the true Jth individual correlation energy

thus leading [see Eq. (244)] to the following exact expression for individual correlation energies:

The above expression is a key result of Ref. [103] which, as we will see, allows us to explore in-principle-exact SD/DD correlation energy decompositions.

Let us now analyze the different contributions on the right-hand side of Eq. (254). While, on the first line, the bare Jth correlation energy component is recovered, the additional terms on the second and third lines ensure that the external potential energy is evaluated with the correct true density (see Eqs. (241) and (245), and the supplementary material of Ref. [103] ). Interestingly, in the summation (in K) over all the states that belong to the ensemble [see the second line of Eq. (254)], one may separate the contribution of the state under consideration (i.e., the Jth state) from the others, thus defining an individual SD correlation energy:

The above definition, which was denoted SD in Ref. [103] (the "overline" notation is dropped in the present work, for simplicity), differs substantially from the definition of Gould and Pittalis [86] . In the latter, an additional statespecific KS wave function, which is expected to reproduce the true individual density of the state under consideration, is introduced. In this case, the name "state-driven" means that the correlation energy is evaluated from interacting and non-interacting wave functions which share the same density. Here, no additional KS wave function is introduced, which is obviously appealing from a computational point of view. One possible criticism about the definition in Eq. (255) is its arbitrariness. Indeed, we may opt for a more density-based definition, in the spirit of what Gould and Pittalis proposed, by introducing, for example, the following auxiliary wave functions:

Note that, in the ground-state w = 0 limit, Φ w 0 reduces to the conventional KS wave function Φ w=0 0 of KS-DFT. What might be interesting in the (somehow artificial) construction of the above individual auxiliary KS states is the possibility it gives to recover, like in the Gould-Pittalis approach [86] , all the KS contributions (to the individual correlation energy) that appear on the first two lines of Eq. (254) from a single expectation value, thus generating, on the other hand, (several) additional terms that should ultimately be removed:

Moreover, according to Eq. (27) , we recover (among other terms) the correct physical density:

On that basis, we could argue that the first two lines on the right-hand side of Eq. (254) should be interpreted as a SD correlation energy, while the third line would correspond to the missing DD correlation energy. The issue with such a decomposition is that the individual DD correlation energies would then cancel out in the weighted sum:

which means that the ensemble DD correlation energy would be zero. As a result, with such an interpretation, the concept of DD correlation would not be of any help in the development of correlation DFAs for ensembles. This is of course not what we want [86] . In this respect, the definition in Eq. (255) is much more appealing. We will stick to this definition from now on. Consequently, the complementary ensemble DD correlation energy will read as [see Eqs. (243), (244), and (255)]

Thus, we recover another key result of Ref. [103] . As readily seen from Eq. (261), the exact evaluation of the DD correlation energy only requires computing the static linear response of the KS wave functions that belong to the ensemble, which is computationally affordable.

We can even establish a direct connection with the total ensemble functional f w [n], by analogy with Eq. (237) . Indeed, since

it comes

The importance of DD correlations, which was revealed in Ref. [86] , has been confirmed in Ref. [103] , in the weakly asymmetric and stronly correlated regime of the two-electron Hubbard dimer. We propose in the following to complete the study of Ref. [103] by exploring all asymmetry and correlation regimes, and comparing exact results with that of standard approximations.

The Hubbard dimer has been introduced in Sec. 4.3. In this simple model system, exact (two-electron and singlet) biensemble density-functional correlation energies E w c (n) can be evaluated through Lieb maximizations [96, 99] from the following analytical expressions for the exact potential-functional interacting energies [136, 84, 85] :

where r = 3(4t 2 + ∆v 2 ) + U 2 (267) and θ = 1

Exact ground-and excited-state densities are then obtained from the Hellmann-Feynman theorem [see Eq. (181)],

and the cubic polynomial equation that the energies fulfill (see the Appendix of Ref. [96] ). The resulting biensemble density reads as n w = (1 − w)n Ψ0 + w n Ψ1 . The Hx-only GOK functional introduced in Eq. (262) can be expressed analytically as follows [96] ,

so that, as shown in Appendix D, the exact DD ensemble correlation energy reads explicitly as

Since the KS excited-state density is always equal to 1 in this model [96] , the prefactor (n Ψ1 − 1) matches the deviation in density of the true physical excited state from the KS one:

Then it becomes clear that E w,DD c (n w ) is a DD correlation energy. We also see from the expression of Eq. (271) that, in the regular ground-state DFT limit (w = 0), this type of correlation disappears.

In the following we test two common DFAs: A (weight-independent) groundstate density-functional description of the ensemble correlation energy (GSec) [97, 96] ,

where E c (n) = E w=0 c (n), and the GS-ic approximation introduced in Sec. 5.1 which, in the present case, gives

Note that the KS ground-state density n Φ w 0 fulfills the constraint (1 − w)n Φ w 0 + wn Φ w 1 = n w , thus leading to

where we readily see that, in general, n Φ w 0 = n Ψ0 . In the following, the local potential will be fixed. It is then analogous to the external potential of ab initio calculations, hence the notation ∆v = ∆v ext .

Let us first discuss the strictly symmetric (∆v ext = 0) dimer in which simple analytical expressions can be derived for both exact and approximate ensemble correlation energies. In this special case, ground-and excited-state densities are equal to 1, in both KS and physical interacting systems. Consequently, the ensemble DD correlation energy vanishes [see Eq. (271)]. Total and SD ensemble correlation energies are equal and vary linearly with the ensemble weight [99] as

with the positive slope −E c (n = 1). As readily seen, the excited state exhibits no correlation effects in this density regime. Turning to the approximations, GS-ic erroneously assigns a (ground-state) correlation energy to the excited state [see Eq. (274)], thus leading to a total ensemble correlation energy that is wrong and equal to E c (n = 1), like in GS-ec [see Eq. (273)]. In conclusion, in the symmetric case, both GS-ic and GS-ec approximations completely miss the weight dependence of the ensemble correlation energy.

We now discuss the performance of GS-ec and GS-ic in the asymmetric dimer. Results are shown in Fig. 4 . The features described in the symmetric case are preserved in the weakly asymmetric regime (see the top left panel of Fig. 4) . When the asymmetry is more pronounced, both exact and approximate ensemble correlation energies exhibit curvature. By construction, these energies all reduce to the same (ground-state) correlation energy when w = 0. They differ substantially by their slope in the ground-state limit (w = 0). Further insight into GS-ic, for example, is obtained from the following analytical expression,

where E c (n = 1)−E c (n = n Ψ0 ) ≤ 0, as readily seen from Fig. 4 of Ref. [99] . Interestingly, in the strongly asymmetric ∆v ext /U 1 regime, the true ground-(n Ψ0 ) and excited-state (n Ψ1 ) densities tend to 2 and 1, respectively (see Fig.  1 of Ref. [96] ). This is the situation where the slope in weight expressed in Eq. (277) reaches its maximum (in absolute value), thus inducing a large deviation from the exact slope, as shown in the bottom left panel of Fig. 4 . At the GS-ec level of approximation, the situation is less critical, at least for small weight values. As readily seen from the following expression [see Eq. (273)],

when ∆v ext ∼ U , the slope (at w = 0) is relatively small since n Ψ1 ∼ n Ψ0 (see Fig. 1 of Ref. [96] ). This is in agreement with the right panels of Fig. 4 . Note that, in this regime, GS-ic can exhibit positive slopes (see the top right panel of Fig. 4) . In this case, the density derivative contribution [second line of Eq. (277)], which is positive [96, 99] , is not negligible anymore and it (more than) compensates the negative correlation energy difference [first line of Eq. (277)]. When the asymmetry of the dimer is more pronounced (i.e., ∆v ext U ), the GS-ec slope (in weight) remains negligible, as shown in the bottom left panel of Fig. 4 . Indeed, in this case, n Ψ0 tends to 2. Moreover, since the ground-state correlation functional expands as follows in the weakly and strongly correlated regimes [99] ,

and E c (n)

respectively, we immediately see that ∂E c (n)/∂n ≈ 0 when n approaches 2, whether U/t is large or small. Note finally that, as already mentioned, in regimes where the asymmetry is weaker than the correlation, i.e., ∆v ext /t < t/U < 1 (see the top left panel of Fig. 4) , the slopes obtained at w = 0 with GS-ec and GS-ic are identical and relatively weak. This can now be understood from the expressions in Eqs. (277) and (278), and the fact that ∂E c (n)/∂n| n=1 = 0 [99] , knowing that n Ψ0 ≈ 1 [96] in this case.

We see in Fig. 4 that the overall weight dependence of the exact ensemble correlation energy differs substantially from that of the GS-ec and GS-ic approximations. It is again instructive to look at the slope at w = 0. It can be expressed exactly as follows,

where we readily see from Eq. (278) that GS-ec neglects the derivative in weight of the ensemble correlation density functional. As highlighted in Eq. (37) [see also Sec. 3.2 for a more detailed discussion in the context of charged excitations], the latter contribution is connected to the derivative discontinuity that the xc potential exhibits when an excited state is incorporated into the ensemble. Since, in the weakly correlated regime [99] , 

which means that, when approaching the equiensemble w = 1/2 case, it systematically decreases with the ensemble weight (because of the term (1−w) 2 in the denominator), unlike the total ensemble correlation energy [see Eq. (284)].

As long as the dimer remains close to symmetric, which requires that ∆v ext reduces as U increases (see Fig. 1 of Ref. [96] ), the numerator in Eq. (289) will be small enough such that DD correlations are at most equal to the total correlation energy. This feature is actually observed in the moderately correlated U/t = 1 regime (see the top left panel of Fig. 5 ). However, in asymmetric and strongly correlated regimes where 0 < ∆v ext U (n Ψ0 ≈ 1 and n Ψ1 > 1 in this case) or ∆v ext ≈ U (i.e., n Ψ0 ≈ n Ψ1 ≈ 1.5) [96] , the numerator is not negligible anymore and, consequently, the DD ensemble correlation energy is significantly lower than the total one (see the bottom left panel of Fig. 5 ; see also Ref. [103] ). In such cases, the complementary SD ensemble correlation energy can be large and positive. This may look unphysical at first sight but, if we return to the definition of Eq. (255), we see that the individual SD correlation energies are not guaranteed to be negative. The reason is that, unlike the total ensemble correlation energy, they are not evaluated variationally. Note finally that, when ∆v ext > U t (n Ψ0 ≈ 2 and n Ψ1 ≈ 1 in this case), the numerator in Eq. (289) will be relatively small, because of the (n Ψ1 − 1) prefactor, thus reducing the energy difference between total and DD correlations (see the bottom right panel of Fig. 5 ).

In summary, with the present SD/DD decomposition [see Eqs. (260) and (261)], both SD and DD correlation energies become relatively large (when compared to the total ensemble correlation energy), especially in the commonly used equiensemble case, and they mostly compensate when the Hubbard dimer has a pronounced asymmetry. This is clearly not a favorable situation for the development of DFAs, which was the initial motivation for introducing the SD/DD decomposition [86, 103] . The latter should definitely be implemented for atoms and diatomics, for example, in order to get further insight. In the case of stretched diatomics, the present study of the Hubbard dimer might be enlightening [141] . We should also stress that, in the asymmetric ∆v ext = U case, standard GS-ic and GS-ec approximations give ensemble correlation energies that are of the same order of magnitude as the exact one, unlike the SD and DD correlation energies. As briefly mentioned in Sec. 5.1, exploring alternative SD/DD decompositions that rely explicitly on GS-ic, which is maybe a better starting point, would be relevant in this respect. Work is currently in progress in this direction. 

Despite the success of time-dependent post-DFT approaches for the description of charged and neutral electronic excitations, various limitations (in terms of accuracy, computational cost, or physics) have motivated in recent years the development of time-independent alternatives. In the present review, we focused on GOK-DFT and N -centered eDFT, which are two flavors of eDFT for neutral and charged excitations, respectively. Their computational cost is essentially that of a standard KS-DFT calculation, because they both rely on self-consistent one-electron KS equations [see Eqs. (16) and (53)]. A major difference though is that, in eDFT, the Hxc density functional is ensemble weight-dependent. This weight dependence is central in eDFT. It allows for the in-principle-exact extraction, from the KS ensemble, of individual (groundand excited-state) energy levels [see Eqs. (36) , (64) and (65)] and densities [Eq. (27) ]. We have also shown that the infamous derivative discontinuity problem that must be addressed when computing fundamental (or optical) gaps [see Eqs. (30) , (61) and (141) functional. Recent progress in the design of weight-dependent xc DFAs have been reviewed. The pros and cons of using an (orbital-dependent) ensemble density matrix functional exchange energy or state-averaging individual exact exchange energies have been discussed in detail. We reveal in passing that, in the latter case, severe (solvable though) v-representability issues can occur when electron correlation becomes strong. Turning to the design of DFAs for ensemble correlation energies, state-of-the-art strategies have been discussed, in particular the combination of finite and infinite uniform electron gas models as well as the recycling of standard (ground-state) correlation DFAs through state-averaging. In the latter case, further improvements may emerge from the concept of density-driven correlation, which does not exist in ground-state KS-DFT. How to define mathematically the corresponding correlation energy is an open question to which we provided a tentative answer [see Eqs. (261) and (263)]. Test calculations on the Hubbard dimer reveal how difficult it is to have a definition that is both rigorous and useful for the development of approximations. Work is currently in progress in other (briefly discussed) directions. Even though it was not mentioned explicitly in the review, we would like to stress that current formulations of eDFT do not give a direct access to exact response properties such as oscillator strengths, Dyson orbitals, or nonadiabatic couplings. Extending Görling-Levy perturbation theory [142, 143, 144] to ensembles might be enlightening in this respect. We recently became aware of such an extension [145] for the computation of excitation energies within the DEC scheme [87, 101] , which is an important first step. Nevertheless, a general quasi-degenerate density-functional perturbation theory based on ensembles, where individual energy levels and properties can be evaluated, is still highly desirable. Work is currently in progress in this direction.

In conclusion, we have summarized in the present review recent efforts of a growing community to put eDFT to the front of the scene. We highlighted several formal and practical aspects of the theory that should be investigated further in the near future in order to turn eDFT into a reliable and low-cost computational method for excited states.

Let us consider the simpler one-dimensional (1D) case in which the KS-PPLB equations read as

thus leading to

where we used the limits vext(∞) = v α H (∞) = 0. Note that ϕ α i (x) is expected to decay as |x| → +∞, which implies −2 ε α i − v α xc (∞) > 0. Therefore,

and

In the true interacting system, the N -electron ground-state wave function Ψ N 0 fulfills

where wee(|x i − x j |) is a well-behaved two-electron repulsion energy in 1D. Let us consider the situation where |x 1 | → +∞ while x 2 , . . . , x N remain in the region of the system, which corresponds to an ionization process in the ground state. Since wee(|x 1 − x j |) → 0, the (to-be-antisymmetrized) wave function and its density can be rewritten as We now turn to the left and right formulations of N -centered eDFT. We recall the shorthand notations (ξ − , 0) notation ≡ ξ − and (0, ξ + ) notation ≡ ξ + . When ξ + > 0, the right N -centered ensemble density, which is mapped onto a non-interacting KS ensemble, has the following asymptotic behavior [we just need to substitute N + 1 for Similarly, for ξ − ≥ 0, we have

Thus, we conclude that 

For convenience, we use the following exponential parameterization of the single-configuration wave functions [126] ,

where κ ≡ {κpq} p<q are the variational orbital rotation parameters andκ is the corresponding real singlet rotation quantum operator. The latter reads as follows in second quantization,

where the index p refers to the orbital ϕ w p andÊpq = τ =↑,↓â † pτâqτ . Therefore, the eDMHF energy becomes a function of κ, where D w (κ) = I w I D Φ I (κ) is a trial ensemble density matrix, and E HF (D) is the conventional ground-state HF density matrix functional energy:

The individual Hx-only GOK energies are extracted from the ensemble one,

as follows, where n ξ,w = (1 − ξ)n Φ w 0 + ξn Φ w 1 .

(D.5)

Since n Φ w 1 = 1 and (1 − w)n Φ w 0 + wn Φ w 1 = n w , (D. 6) or, equivalently,

it comes

and

From the weight derivative expression 

(D. 13) or, equivalently,

it comes

As a result, × 2t

(D. 19) 

Proc. Natl. Acad. Sci

Molecular electronic-structure theory

We conclude from Eq. (235) that the energy levels can be evaluated exactly within GOK-DFT as follows,where the individual correlation energy of the Jth state is determined from the ensemble correlation density functional as follows,In the following section, we will see how the concept of DD correlation emerges from Eq. (242), once it has been rewritten more explicitly in terms of individual densities.

According to the definition of the ensemble correlation functional in GOK-DFT [see Eq. (25) ], the exact ensemble correlation energy can be decomposed as follows,where the individual components read asLet us stress that these components do not match the individual correlation energies of Eq. (242). Indeed, unlike the latter [see Eq. (241)], they do not give access to the exact individual energy levels,Therefore, as long as the ground state does not deviate too much from the symmetric n Ψ0 = 1 density profile, which is the case when ∆v ext U , the exact slope is not negligible, and it is positive. This is in agreement with the top left panel of Fig. 4 . Interestingly, in this density regime, this feature is preserved when the strength of electron correlation increases (not shown). Indeed, in this case, the ensemble correlation functional reads as [99] E w c (n)thus leading to ∂E w c (n Ψ0 )/∂w ≈ U/2. When the dimer is strongly asymmetric, the ground-state density approaches 2 and, in this case [99] ,so thatthus leading to ∂E w c (n Ψ0 )/∂w| w=0 ≈ 0. As readily seen from Eq. (283), the same result is obtained in the weakly correlated regime. This is in complete agreement with the bottom left panel of Fig. 4 . It also explains why the deviation of GS-ec from the exact result drastically reduces when ∆v ext increases for a fixed interaction strength U and relatively small weight values. Finally, in the particular case where ∆v ext = U , the computed ground-state densities equal n Ψ0 ≈ 1.30 and n Ψ0 ≈ 1.46 in the moderately U/t = 1 and strongly U/t = 5 correlated regimes, respectively. As expected from Eqs. (283) and (287), the exact slope will be substantial and negative, which agrees with the right panels of Fig. 4 .We now focus on the exact SD/DD decomposition of the ensemble correlation energy. Results are shown in Fig. 5 for various correlation and asymmetry regimes. In the ground-state limit, the slope of the DD ensemble correlation energy reads as [see Eq. (271)]

The authors declare that they have no conflict of interest.By construction, the minimum is reached when κ = 0, and we denote D w = D w (κ = 0). Note thatwhere the occupation number n I p is an integer, andwhere θ w p can be fractional. The stationarity condition that is fulfilled by the minimizing eDMHF orbitals can now be written explicitly as follows,whereis the conventional density matrix functional Fock operator matrix element, and 

We use the same parameterization as in Appendix B, i.e.,where the indices {p} in creation/annihilation operators (as well as in one-and two-electron integrals) now refer to the minimizing SAHF orbitals φ w p . The to-be-minimized SAHF energy can be expressed as follows,so that the stationarity condition reads aswhere v w Therefore, if we use the notation φ w r |ĥ +v w Hx,I |φ wEq. (C.3) can be rewritten in a compact form as follows, 

For convenience, we will use the following exact expression for the ensemble DD correlation energy: