N-centered ensemble density-functional theory for open systems

Two (so-called left and right) variants of N-centered ensemble density-functional theory (DFT) [Senjean and Fromager, Phys. Rev. A 98, 022513 (2018)] are presented. Unlike the original formulation of the theory, these variants allow for the description of systems with a fractional electron number. While conventional DFT for open systems uses only the true electron density as basic variable, left/right N-centered ensemble DFT relies instead on (i) a fictitious ensemble density that integrates to a central (integral) number N of electrons, and (ii) a grand canonical ensemble weight $\alpha$ which is equal to the deviation of the true electron number from N. Within such a formalism, the infamous derivative discontinuity that appears when crossing an integral number of electrons is described exactly through the dependence in $\alpha$ of the left and right N-centered ensemble Hartree-exchange-correlation density functionals. Incorporating N-centered ensembles into existing density-functional embedding theories is expected to pave the way towards the in-principle-exact description of an open fragment by means of a pure-state N-electron many-body wavefunction. Work is currently in progress in this direction


I. INTRODUCTION
Density-functional theory (DFT) has become over the last two decades the method of choice for performing routine large-scale electronic structure calculations. This success is essentially due to the relatively low computational cost of the method and its (often but not always) good accuracy. In most applications, DFT is applied to closed electronic systems, i.e. systems with an integral number N of electrons. However, at the formal level, there is no such a restriction. In other words, DFT is in principle able to describe also systems with a fractional electron number, as shown in the pioneering work of Perdew, Parr, Levy and Balduz (PPLB) [1]. The extension of DFT to open systems plays a crucial role in the description of charged electronic excitations [2]. More recently, it became a key ingredient in the derivation of DFT-based embedding approaches such as partition DFT [3][4][5][6][7][8][9][10], potential-functional embedding theory [11], and frozen density embedding theory for non-integer subsystems' particle numbers [12].
The density n(r) of an open system integrates to a fractional number N of electrons. As shown in Refs. [1,2], the ground-state energy of such a system varies linearly with the deviation α = N − ⌊N ⌋ of the true electron number from its floor integral value. In the language of ensemble DFT [6,13], the density is nothing but a weighted sum of ground-state densities n N −1 (r) and n N (r) integrating to N − 1 = ⌊N ⌋ and N = ⌈N ⌉ electrons, respectively: n(r) = (1 − α)n N −1 (r) + α n N (r), (1) * Corresponding author; bsenjean@gmail.com where, as readily seen, the ensemble weights are α and (1 − α), respectively. As shown in Ref. [13], the explicit expression in Eq. (1) is convenient for constructing density functional approximations to the open-system exchange-correlation (xc) energy E xc [n] as functions of α.
Quite recently, the authors have proposed an inprinciple-exact reformulation of the fundamental gap problem in DFT [14]. For that purpose, they introduced the concept of N -centered ensemble where, unlike in the true physical ensemble described in Eq. (1), the ensemble density integrates to an integral number of electrons. The practical advantage of such a reformulation is that the infamous derivative discontinuity contribution to the true gap [2] can be expressed as an ensemble Hartree (H) xc energy derivative with respect to the ensemble weight [14], exactly like in DFT for canonical ensembles [15]. In the original formulation of the theory [14], the ensemble weights have no physical meaning. They are just auxiliary variables which, through their variations, enable the extraction of quantities of interest like the ionization potential (IP), the electron affinity (EA) or, when a single weight is used, the fundamental gap. In fact, N -centered ensemble DFT allows for a direct extraction of individual energies for the neutral, anionic and cationic systems [14], which are all closed systems, in a single calculation. In this work, we show how the theory can be adapted to open systems. By considering ionization and affinity processes separately, we obtain two variants (referred to as left and right) of N -centered ensemble DFT where the ensemble weight is now connected to α [see Eq. (1)].
The paper is organized as follows. After a brief review on the original formulation of N -centered ensemble arXiv:1912.07125v3 [cond-mat.str-el] 1 Feb 2020 DFT in Sec. II A, the concept of left and right N -centered ensembles is introduced (Sec. II B). The left/right variants of N -centered ensemble DFT are then derived in Sec. II C. In this context we obtain an in-principle-exact reformulation of the IP/EA theorem, as discussed in Sec. II D. An explicit connection (through scaling relations) between left and right ensemble density functionals is then made in Sec. II E. Finally, a brief discussion on density-driven correlation effects, whose importance in canonical ensembles was recently revealed [16][17][18], is proposed in Sec. II F. As a conclusion to the "Theory" Sec. II, we compare left/right N -centered ensemble DFT with conventional PPLB-DFT for open systems (Sec. II G). A proof of concept study of the Hubbard dimer model is presented in Sec. III, with a particular emphasis on the weight-dependence of the Hxc left/right ensemble functionals. Conclusions and perspectives are finally given in Sec. IV.

A. N -centered ensembles
In a previous work [14], we introduced the concept of N -centered ensemble for the purpose of calculating fundamental gaps, in principle exactly, within DFT. In standard approaches, the number of electrons is considered to vary continuously between two integers, thus allowing for the description of ionization or affinity processes [13,. The situation is completely different in an N -centered ensemble where, by construction, the number of electrons (which is obtained by integration of the N -centered ensemble density) is not affected by the charged excitation. It remains equal to the so-called central number N of electrons which is an integer. Formally, an N -centered ensemble is described by the following density matrix operator [14]: where ξ − and ξ + are ensemble weights assigned to the cationic (N − 1)-electron and anionic (N + 1)-electron systems, respectively. Each individual density matrix operatorΓ Ne , where N e = N, N ± 1, is used to describe an N e -electron ground state. If the latter can be described with a pure-state wavefunction Ψ Ne then Γ Ne ≡ Ψ Ne ⟩ ⟨Ψ Ne . In case of degeneracy, it may be written as a convex combination of pure N e -electron statesΓ Ne ≡ ∑ I λ I Ψ Ne I ⟩ ⟨Ψ Ne I , where ∑ I λ I = 1. As readily seen from Eq. (2), the number of electrons within the N -centered ensemble is Tr Γ {N,ξ−,ξ+}N = N where Tr [. . .] denotes the trace,N = ∫ drn(r) is the electron counting operator, andn(r) is the density operator at position r.
As shown in Ref. [14], using two independent weights ξ − and ξ + allows for a separate extraction of the IP and EA from the N -centered ensemble energy Tr Γ {N,ξ−,ξ+}Ĥ , whereĤ is the (second-quantized) Hamiltonian. In the particular case where ξ − = ξ + = ξ, the N -centered ensemble density matrix operator simplifies as follows:Γ so that the fundamental gap can be extracted directly by differentiating the corresponding N -centered ensemble energy Tr Γ {N,ξ}Ĥ with respect to ξ [14]. The potential advantage of such a formalism over conventional approaches is that the infamous derivative discontinuity contribution to the gap can be written as the derivative of the N -centered ensemble Hxc density functional with respect to the ensemble weight ξ, exactly like in DFT for canonical ensembles [15,41].
Let us stress that, in the current formulation of Ncentered ensemble DFT, the ensemble weights have no physical meaning. They are just convenient auxiliary variables. Actually, the properties of interest [namely the ground-state energies of the neutral (N -electron), cationic and anionic systems] should not depend on the value of the ensemble weights. In the rest of this work, we explore two variants of the theory which are directly applicable to open systems. For such systems, a convenient and physical choice for the ensemble weight value is the deviation of the (fractional) electron number N from either its integral floor value ⌊N ⌋ or the ceiling one ⌈N ⌉.

B. Left and right N -centered ensembles
Starting from Eq. (2), we explore in this section different choices of N -centered ensemble weights ξ − and ξ + for the purpose of describing an open N -electron system. In order to treat both ⌈N ⌉ = N and ⌊N ⌋ = N scenarios, thus allowing for the description of fluctuations around the central (integral) number N of electrons, we introduce what we will refer to as left (subscript '−') and right (subscript '+') N -centered ensemble density matrix operators, respectively: where 0 ≤ α ≤ 1. Note that, by construction, the number of electrons within these two ensembles is the central number N : Moreover, left and right ensembles can be connected to each other as follows: These ensembles are just special cases of the original Ncentered one introduced in Ref. [14]. Indeed, with the notations of Eq. (2), we havê The original (two-weight) N -centered ensemble can actually be recovered from the left and right ones as follows: Following the standard description of open systems in DFT [1], we propose to use for the value of α in Eqs. (4) and (5) the deviation of the true (fractional) electron number N , that we assume to vary in the range N − 1 < N < N + 1, from the central number N , i.e. α = N − N .
If ⌈N ⌉ = N , then α = N − N and the open system will be described by the left N -centered ensemble of Eq. (4). If ⌊N ⌋ = N , then α = N − N and the right N -centered ensemble of Eq. (5) will be used instead. This is illustrated graphically in Fig. 1.
Let us stress that these ensembles are not the physical ones, which are described by the following density matrix operators [13]: thus leading to the final expressionŝ andΓ (N ) A direct consequence of Eqs. (15) and (16) is that the physical energy E (N ) = Tr Γ (N )Ĥ of an open system can be extracted from a left or right N -centered ensemble as follows: and where E and Let us finally point out that, according to Eqs. (13) and (14), Consequently, the IP and EA can be extracted from the left and right N -centered ensemble energies as follows: and As shown in the following, by deriving a DFT for left and right N -centered ensembles, we obtain from Eqs. (17), (18), (19), and (20) a novel and in-principleexact density-functional description of open systems. In this section we derive a variational Kohn-Sham (KS) DFT expression for the left/right N -centered ensemble energies. As a result, the open system problem will be mapped [via Eqs. (17) and (18)] onto an auxiliary N -electron non-interacting ensemble, where the (central) number N is an integer. The derivation follows closely the one presented in Ref. [14] for two-weight N -centered ensembles [see Eq. (2)].
Let us consider the usual ab-initio quantum chemical HamiltonianĤ =T +Ŵ ee +V ext whereT andŴ ee are the kinetic energy and two-electron repulsion operators, respectively, andV ext = ∫ dr v ext (r)n(r) is the external potential operator. For an isolated molecule, the external local potential v ext (r) is simply the nuclear potential. If we introduce the analog for left/right N -centered ensembles of Levy's constrained-search universal functional, where the minimization is performed over left/right Ncentered ensemble density matrix operatorsγ then it becomes possible to express the exact left/right N -centered ensemble energies variationally as follows: The minimum is reached when the density equals the exact left/right N -centered ensemble one nΓ{N,α} ± . Note that the true (i.e. interacting) density matrix operatorŝ Γ {N,α} ± in Eqs. (4) and (5) would be obtained by solving the ground-state Schrödinger equationĤΨ Ne = E Ne Ψ Ne for N e = N and N e = N ± 1 electrons. Like in regular DFT, we bypass this complicated task by introducing the following KS decomposition, where the non-interacting density-functional kinetic energy contribution reads and E

{N,α}
Hxc± [n] is the left/right N -centered ensemble Hxc functional. Note that this functional is, for a fixed density n, α-dependent. It differs from the conventional (Nelectron) ground-state Hxc functional E Hxc [n], which is recovered when α = 0: As shown in the following, modeling the dependence in α of the Hxc functional in this context is analogous to modeling the derivative discontinuity in conventional DFT for open systems. In the light of Ref. [42], we define the exact Hx contribution as follows: whereγ

{N,α} s±
[n] is the minimizing left/right N -centered ensemble density matrix operator in Eq. (29). According to Eqs. (25), (28), and (29), the correlation contribution can then be expressed as Let us now return to the variational ensemble energy expression of Eq. (27). By inserting the exact decomposition of Eq. (28) we obtain, according to Eq. (29), the final expressions whereĥ =T +V ext . Like in regular ensemble DFTs, the minimizing (non-interacting) density matrix operatorsγ and nγ{N,α} The KS orbitals are obtained by solving the following self-consistent equations, . (40) Note that like in the original version of N -centered ensemble DFT [14], since and Moreover, by applying the Hellmann-Feynman theorem to the variational energy expressions in Eq. (33), it comes and, similarly, By combining Eqs. (17), (18), and (42)-(45) we can finally express the exact open-system energy in terms of the LZ-shifted KS orbital energies as follows: and Eqs. (46) and (47) are the central result of this work.
As readily seen from these equations, the exact physical energies are recovered by summing up the LZ-shifted occupied KS orbital energies not only in the N -electron system (i.e. when α = 0), as expected from Ref. [43], but also in the (N ± 1)-electron ones (i.e. when α = 1). This is a non-trivial result since the left/right N -centered ensemble densities do not reduce exactly to the physical ones when α = 1 [see the prefactors 1±(α N ) in Eqs. (19) and (20)]. In addition, keeping in mind that the physical energies E(N ± α) vary linearly with α, we expect the total LZ-shifted KS energies for N − α and N + α electrons, respectively, to vary (at least) quadratically with α. Linearity should then be recovered when adding the ensemble Hxc first-orderderivative corrections [third terms on the right-hand side of Eqs. (46) and (47)]. This is illustrated in Fig. 2 with the Hubbard dimer model [see Sec. III for further details].
D. Reformulation of the IP/EA theorem As shown in Sec. II B, the left/right N -centered ensemble energies (and their first-order derivatives in α) give directly access not only to the physical energy for any fractional electron number N (i.e. any value of α), but also to the IP and the EA [see Eqs. (23) and (24)]. From the KS-DFT expressions in Eqs. (42)-(45) we obtain the following reformulation of Janak's theorem [45], which holds for any value of N or, equivalently, for any value of α in the range 0 ≤ α ≤ 1: and This is the second key result of the paper.
In the limits N → N ± η where η → 0 + , which consists in taking the α = 0 limit in both left and right N -centered ensembles, the density becomes simply the conventional N -electron ground-state density: Consequently, we have and, in particular, where ε N and ε N +1 are, respectively, the energies of the HOMO and LUMO obtained from a regular N -electron KS-DFT calculation. Note that the LZ shift does not affect the HOMO-LUMO gap. Therefore, according to the IP/EA expressions in Eqs. (50) and (51), the exact fundamental gap can be written as follows: where the last two terms on the right-hand side play all together the role of the derivative discontinuity correction to the gap in conventional KS-DFT. Note that, according to Eqs. (8) and (9) and E

{N,α}
Thus we recover the compact expression [14]: Interestingly, designing density-functional approximations for the original (two-weight) N -centered ensembles benefits automatically to the left/right variants of the theory through Eqs. (57) and (58).

E. Left or right ?
For clarity and convenience, we have derived explicitly both left and right variants of N -centered ensemble DFT. This is in principle not necessary as we may opt for one of the ensemble (left or right) and simply increment or decrement the central number of electrons. As illustrated in Fig. 1, an open N -electron system with ⌈N ⌉ = N can be described either by a left ensemble centered on N or a right ensemble centered on (N − 1). In other words, one may express the number of electrons as or, equivalently, where α ′ = 1 − α. The two descriptions should of course lead to the same physical density and energy. It is actually quite simple to obtain the left ensemble Hxc density functional from the right one, and vice versa. Indeed, according to Eqs. (7) and (25), where n is an N -electron density. Similarly, we can show that, if n is an (N − 1)-electron density, then .
As these scaling relations apply also to the noninteracting kinetic energy, we conclude that .
In summary, a single (left or right) N -centered Hxc functional is needed for describing the variation of the electron number between two integers. One functional is obtained from the other via the scaling relations of Eqs. (64) and (65).

F. Density-driven correlations in N -centered ensembles
Very recently, Gould and Pittalis [16,17] revealed that direct approximations to Gross-Oliveira-Kohn (GOK) ensemble DFT, where the ensemble consists of N -electron ground and excited states [15], miss an important correlation effect that they refer to as density-driven correlation. The latter originates from the deviation in density of the individual KS states within the ensemble from the true interacting ones, as shown explicitly by one of the author [18]. One may naturally wonder if this kind of correlation exists also in N -centered ensemble DFT. In the following, we will briefly explain how density-driven correlations can be defined in this context. We leave their detailed study for future work.
Let us follow the derivation of Ref. [18] and adapt it to (say left) N -centered ensembles. We start by extracting the individual M -electron (M = N or N −1) total energies E(M ) = Tr Γ MĤ . For that purpose we use Eqs. (13) and (14) respectively. One should then realize that the density constraint in Eq. (29) does not necessarily imply that the individual KS densities match the true interacting ones: This can be illustrated easily with the asymmetric Hubbard dimer model, which will be studied in detail in Sec. III. In this model, the one-electron ground-state density (which is an atomic site occupation) reads where t is the strength of the kinetic energy and ∆v is the analog of the local potential (here it is just a number that controls the asymmetry in the dimer). When the dimer is asymmetric (i.e. ∆v ≠ 0), the Hxc potential [which is the difference in potential between the KS noninteracting N -centered ensemble and the interacting one] is nonzero. In the weakly correlated regime, this can be readily seen from Eq. (79). Therefore, according to Eq. (71), the true and KS one-electron states [each of them being a component of a (left) N -centered ensemble with N = 2] do not have the same density. If we now return to the general case, it is then relevant to decompose each individual correlation energy into state-driven (SD) and density-driven (DD) contributions, by analogy with GOK ensembles [18]: While we do not expect the conventional ground-state functional E Hxc [n] to be a good approximation to (neutrally) excited-state Hxc energies, one may wonder how good the approximation E

{N,α}
Hxc−,M n, n M ≈ E Hxc n M is since an N -centered ensemble consists of ground states only. Answering this question is expected to pave the way toward the rationalization and development of densityfunctional approximations that incorporate derivative discontinuity corrections. Work is currently in progress in this direction.

G. Comparison with conventional DFT for open systems
In standard PPLB-DFT for open systems [1], the density n of the true physical open system is used as basic variable. It gives, by integration, the total (fractional) number N of electrons: The Hxc energy is then obtained from the universal ground-state functional E Hxc [n] which is defined over the domain of densities that can integrate to any integral or fractional number of electrons. Despite its formal beauty, this grand canonical formulation of DFT is not necessarily the most appealing one when it comes to perform practical calculations. Indeed, modeling E Hxc [n] accurately for any number of electrons (including fractional ones) is not straightforward. Most importantly, the model should be able to reproduce the derivative discontinuity that the functional is expected to exhibit when crossing an integer [2]. Note that, for open systems, the KS potential is truly unique (not up to a constant anymore). Indeed, the chemical potential has to be adjusted such that the grand canonical ground-state energies of the (N − 1)-and N -electron systems [we assume that N − 1 < N < N ] are equal.
The situation is completely different in (say) left N -centered ensemble DFT, where we use two variables instead. The first one is the left N -centered ensemble density that integrates to the so-called central number N = ⌈N ⌉ of electrons. The second variable α = N − N is the deviation of the central number from the true electron number N . Even though alternative ensemble DFT approaches [13,36] do not rely on a centered density, they use information about the (N − 1) and N -electron systems (i.e. systems with an integral number of electrons) as N -centered ensemble DFT does. They are similar from this point of view.
Since a density n that integrates to N can be reproduced by either a pure N -electron ground state or a left N -centered ensemble, it is essential to construct a functional E

{N,α}
Hxc− [n] that is both density-and αdependent. At first sight, this version of DFT for open systems looks much more complicated than the PPLB one. What makes the N -centered formulation appealing is that the infamous derivative discontinuity is now obtained through the derivative in α of the N -centered ensemble density-functional Hxc energy. As a result, if one can model the dependence in α of the latter functional, one can in principle solve the derivative discontinuity problem.

A. Hamiltonian and density functionals
A proof of concept study of the (not necessarily symmetric) Hubbard dimer is presented in this section. This is one of the simplest solvable models exhibiting a nontrivial interplay between electron-electron interaction and inhomogeneity. It is often used as a lab for testing new ideas in DFT, gaining more insight into those and their approximate formulations [14,44,[46][47][48][49][50][51][52][53]. In this model, the ab initio HamiltonianĤ is simplified as follows:T wheren i = ∑ σ=↑↓niσ is the density operator on site i (i = 0, 1). Note that the external potential reduces to a single number ∆v ext 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 = N − n for an N -centered ensemble. In the following, the central number of electrons will be fixed to N = 2 and the hopping parameter will be set to t = 1. All density functionals can be derived analytically except the correlation ones that can be computed to arbitrary accuracy via Lieb maximizations [14,44,48].
As shown in Appendix A, within left/right N -centered ensemble DFT, the exact ground-state energy of an open system with N ± = 2 ± α electrons reads . (76) According to Eqs. (8), (9), (29), and (31), the expressions for the left/right N -centered ensemble non-interacting kinetic and Hx energy functionals can be obtained from the two-weight N -centered ensemble functionals of our previous work [14]. Indeed, the left functionals correspond to ξ − = N α N −1 and ξ + = 0, thus leading to where ∆v ∂n is the KS potential for the left ensemble, and Similarly for the right functionals, we have ξ − = 0 and ξ + = N α N +1 , thus leading to and As readily seen from Eq. (81), the non-interacting vrepresentability condition for the right N -centered ensemble is α-dependent: Various density-functional approximations are tested in the following. For simplicity, we will restrict our study to functional-driven errors [54] which means that all density-functional energies will be computed with the exact N -centered ensemble densities. The investigation of density-driven errors [54] is left for future work. Let us start with two approximations that would be applicable also to ab initio problems. The simplest one consists (n) + Ec(n) Eqs. (79), (82), Ref. [47] in neglecting correlation effects and using a regular (αindependent) Ground-State exchange functional. In this model, we simply have to consider the α = 0 limit of Eqs. (79) and (82). The approximation will be referred to as GS-Hx. A conventional (α-independent) Ground-State correlation functional might then be added, thus leading to a second approximation referred to as GS-Hxc. In this context, we use the correlation functional of Carrascal et al. [47] which has been parameterized on the two-electron Hubbard dimer. Finally, in order to investigate the importance of the α-dependence in both exchange and correlation ensemble energies, we consider two additional approximations. In the first one, which is referred to as ensemble exact exchange (EEXX), Eqs. (79) and (82)  As shown in Eqs. (46) and (47), the true physical energy of an open system cannot be expressed solely as a function of the shifted KS orbital energies. An additional ensemble Hxc first-order-derivative correction has to be accounted for as soon as the true number of electrons deviates from an integer. In Fig. 2, we plot the exact expressions for the true physical energy using the left and right (N =2)-centered ensembles, with and without the ensemble derivative-in-α Hxc corrections. As expected, the exact energy is a piecewise linear function of the electron number. Removing the derivative-in-α correction induces curvature. Hence, in the language of left/right N -centered ensemble DFT, describing the piecewise linearity of the energy consists in modeling the dependence in α of the ensemble Hxc density functional. As clearly shown in Fig. 2, no derivative-in-α correction is needed when the number of electrons is an integer. This is a direct consequence of the scaling relations in Eqs. (64) and (65), and the fact that, by construction, summing up the LZ-shifted occupied KS orbital energies gives the exact  48) and (49), and the text that follows for further details.
energy for an integral number of electrons [43].

C. Symmetric case
In this section, we investigate the performance of the various approximations (see Tab. I) in the symmetric dimer (∆v ext t = 0) for different values of U t, as shown in Fig. 3. As expected, the GS-Hxc energy is always below the GS-Hx one, because the correlation energy functional is always negative (and there is no additional derivative-in-α correction). This statement also holds when comparing GS-c with EEXX. Since we use the highly accurate parametrization of Carrascal et al. [47] for the two-electron GS correlation energy functional, the GS-Hxc and GS-c approximations are on top of the exact curve for N = 2. This is definitely not the case for the GS-Hx approximation which describes neither the correlation effects nor the α-dependence. This is the reason why GS-Hx performs so poorly for any N value, especially when U t is large. Note that it also gives a wrong energy for N = 1. Indeed, even if there is only one electron (and therefore no correlation), the Hx part of the functional is still α-independent within GS-Hx, and therefore only meaningful when N = 2 (i.e. α = 0).
Returning to the GS-Hxc and GS-c approximations, they are both exact for N = 2. When deviating from this central electron number, the α-independent correlation functional of Carrascal et al. becomes an approximation, as one cannot expect the correlation energy to be the same for systems with different numbers of electrons. In contrast, EEXX is exact for N = 1, by construction. It is also exact for N = 3, which is due to the particle-hole symmetry of the model. In the range 1 < N < 3, EEXX is not exact anymore due to the lack of the α-dependent correlation functional.
The effect of the latter on the true physical energy can be directly evaluated by Interestingly, all the approximations give physical energies that vary linearly with N (or α). As shown in Appendix B, this is due to the fact that, in the symmetric dimer, the left and right (N = 2)-centered ensemble densities are equal to n = 1.

D. Asymmetric case
Let us now investigate the asymmetric case (we choose ∆v ext t = 5) for different values of U t. As shown in Fig. 4, unlike in the symmetric case, approximate energies exhibit curvature in N . In the particular case of GS-Hx and GS-Hxc (which are α-independent functionals), the curvature is due to the α-dependence in both the non-interacting kinetic energy density functional and the interacting N -centered ensemble density. In the weaklycorrelated regime (U t = 1 and U ∆v ext = 1 5) one can- not distinguish the exact curve from the EEXX one, showing that the ensemble (α-dependent) correlation density functional is equal to 0 in this regime, for any α. This is due to the fact that the ensemble densities reach the border of their non-interacting v-representability domain when ∆v ext ≫ U and ∆v ext ≫ t. For this particular value of the ensemble density, the ensemble correlation density functional vanishes. Note that the same behaviour occurred in the original N -centered ensemble DFT, for which the explanation can be found in Appendix C.3 of our previous work [14]. Because the vrepresentability domain of the left N -centered ensemble density is α-independent, the two-electron (α = 0) correlation density functional is also very close to 0 in between 1 ≤ N ≤ 2, such that GS-c is now on top of the exact curve just like EEXX. It also explains why GS-Hxc and GS-Hx give the same true physical energy here but, in contrast to EEXX and GS-c, they are not exact thus highlighting the importance of the (α-dependent) ensemble Hx functional in this regime. For 2 ≤ N ≤ 3, the two-electron correlation energy is not equal to zero such that EEXX and GS-Hx differ slightly from GS-c and GS-Hxc, respectively. Hence, GS-c is not exact anymore, in contrast to EEXX. As shown in the lower panels of Fig. 4, the curvature of the approximate energies becomes more pronounced as U t increases. As expected, EEXX is not accurate anymore, especially around N = 2, where the correlation energy becomes important. Unlike in the stronglycorrelated symmetric case, GS-c reproduces essentially the exact energy when N = 1 (i.e. α = 1). In this case, the left N -centered ensemble density is actually equal to n Thus we conclude that, for large ∆v ext t values, the two ensemble densities reach the border of their noninteracting v-representability domains [see Eqs. (78) and (81)], such that the exact (α-dependent) ensemble density-functional correlation energy becomes zero. For the left ensemble, it is clear that the α-independent correlation energy used in the GS-c approximation is also equal to zero, as the border of the v-representability domain does not depend on α. Therefore, both GS-c and EEXX become exact when N = 1. However, when N = 3, GS-c introduces a spurious correlation energy contribution obtained by inserting the right ensemble density value n = 2 − 2 3 into the two-electron (α = 0) groundstate correlation functional.

IV. CONCLUSIONS AND PERSPECTIVES
So-called left and right variants of N -centered ensemble DFT have been explored. Unlike in the original formulation of the theory [14], both open and closed electronic systems can be described.
In left/right N -centered ensemble DFT, the key variables are the left/right ensemble density, which integrates to the central integral number N of electrons, and the (absolute) deviation α of the true electron number from N . Within such a formalism, the infamous derivative discontinuity is taken care of by the α-dependent N -centered ensemble Hxc functional. Its α-dependence plays also a key role in reproducing the correct piecewise linearity of the energy, as illustrated in this work with the Hubbard dimer model. What we learn from this model is that conventional (α-independent) xc functionals are not sufficient for describing open systems in the context of N -centered ensemble DFT. Developing ab initio α-dependent density-functional approximations is a challenging but necessary task. As a starting point, a (semi-)local approximation could be obtained, for example, by applying the theory to finite uniform electron gases [55]. Work is currently in progress in this direction.
Various applications of left/right N -centered ensemble DFT can already be foreseen. Regarding densityfunctional embedding techniques, it could be used in place of conventional DFT for grand canonical ensembles when describing open fragments. The practical advantage would come from the explicit description of derivative discontinuity contributions to the energy through the dependence in α of the Hxc functional. Another even more appealing feature of the N -centered formalism is the possibility to use a (pure-state) N -electron many-body wavefunction for describing an open fragment. Even though such an idea seems counterintuitive and difficult to implement, it can be made formally exact if appropriate density-functional corrections are introduced. The basic idea is to consider the following decomposition: where nΓ{N,α} ± is the left/right N -centered ensemble density [from which the density of the open system can be extracted, according to Eqs. (19) and (20) where Ψ

{N,α} ±
is an auxiliary many-body wavefunction with density nΓ{N,α} The expressions for the ground-state energy of the twoelectron Hubbard dimer and its derivatives are given in the Appendix of Ref. [48].
for GS-Hxc, and finally for GS-c, Note that in the case N = 2 − α, the slope with respect to N in Figs. 3 and 4 are the same as the above but with opposite sign. As readily seen from Eqs. (B14) and (B17) in the atomic limit t → 0, the energy of the open system within GS-Hxc has the same slope (equal to zero) as the exact one for 1 ≤ N ≤ 2. As GS-Hxc is exact at N = 2, GS-Hxc is also exact in the range 1 ≤ N ≤ 2 in the atomic limit. Another interesting behaviour, in the symmetric case and in the atomic limit, is that none of the approximations considered in this article can describe the derivative discontinuity of the grand canonical energy when crossing N = N = 2. Indeed, all the approximations feature a constant slope in between 1 ≤ N ≤ 3, making no difference between the ionization potential and the electronic affinity, and thus no fundamental gap. This is obviously not correct, as the exact fundamental gap is equal to U in this limit. Therefore, in order to reproduce the correct behaviour, weight-dependent correlation functionals are required. Such developments are left for future work.