Mathematical Approach to N-Centered DFT and Associated Methods for the Study of Open Systems

A fundamental aspect of the study of N − electronic systems (systems containing N electrons) is to obtain information on the states in which these systems have minimal energy. In practice a numerical search of such states is impossible to carry out, so that alternative approaches have been developped, the one around which this work revolves being to consider electronic systems through their electronic density rather than their state. This approach, known today as Density Functional Theory (DFT), was formalised in Kohn and Sham’s seminal article [1] and its mathematical aspects were studied a few years later by Lieb [2]. Since then, the ideas leading to the construction of DFT have been adapted to the context of electronic systems with a fractionnal number of electrons (open systems), ﬁrst through PPLB DFT[3] and more recently through the deﬁnition of N − centered DFT[4, 5]. In both cases it is unclear wherether the mathematical properties established for classical DFT can be expected to hold true. This question is the main problematic of our work, in which we shall study the analogy between N − centered and classical DFT, from their construction to the methods that are derived from them. This will lead us to construct a Kohn-Sham scheme for N − centered DFT, investigate the links between this theory and optimal transport and present the Hubbard Dimer in this particular situation.


Introduction
The energy of an open system with fractionnal electron number N + ξ is expressed as a linear combination of energies associated to systems with integer electron number : It can be readily seen that variations in ξ of this quantity give access to the electron affinity (EA) E N +1 − E N , which could then be expressed after setting up a Kohn-Sham method in this context through approximations of the variation in ξ of an open-system exchange correlation energy E ξ xc . The modelisation of this dependance is problematic since E ξ xc is defined only for densities integrating to N + ξ, so that the study of its weight-dependance for a fixed density is impossible. As an alternative to this formulation we can choose to work with different weights. This leads to the definition of N -centered DFT, for which the domain of the exchange energy no longer depends on ξ. In both cases, the electronic problem is formulated as a minimization over weighted ensembles w N Γ N ⊕ w N +1 Γ N +1 of their ensemble energy. Although a constrained search can be set up in order to define a universal functional for such systems, it is unclear whether a strong analogy can be expected to hold with classical DFT theory since having modified the set of states defining the electronic problem could very well have modified the nature of the functional we obtain.
In Section 1 we justify that the work done by Lieb in order to propose a precise definition of DFT [2] can still be followed to obtain a similar construction of N −centered DFT. As in classical DFT we shall see that the main question is to prove lower semicontinuity of the N −centered functional. From this construction and the formalism we will have defined we derive the Kohn-Sham formulation of N −centered DFT in Section 2. There we shall see that the regularization method proposed in [6] in order to answer some of the problems inherent to this method can still be set up. This will be seen, in Section 3, not to be the case for the study of the semi-classical limit, which gives rise to a problem that doesn't take the excpected form of an optimal transport problem with Coulomb cost.
In the last section, we propose a comprehensive introduction to the Hubbard dimer model, an extremely simplified model of an electronic system whose major purpose is to be used as a toy model to "test" ideas proposed in DFT. In this context we shall study the associated N −centered functionals and show that this model may not correctly reproduce some fundamental aspects of N −centered DFT.

Mathematical formulation
We shall start by setting the notions and vocabulary which will enable us to define an N −centered universal density functional, the fundamental object of any DFT method, and show the expected duality between this functional and the N −centered energy functional holds.

N −centered states and densities
An N −electronic state is a self-adjoint trace class linear operator on L 2 a (R 3N ) (the Hilbert space of square integrable antisymmetric functions of N spatial variables, endowed with the usual L 2 scalar product) such that Tr(Γ) = 1 and Tr((−∆) 1/2 Γ(−∆) 1/2 ) < ∞.
The set of N −electronic states is written S(N ).
To any two given real numbers 0 ≤ ξ − , ξ + such that N 0 ≤ w N (ξ − , ξ + ) : we associate the set S One of the notable properties of trace-class operators acting on a space of square integrable functions is that they can be described via kernels. More precisely, following [7], Proposition 7.13, we obtain the following result regarding the kernels of N −electronic states. where L 2 a×a (R 3N × R 3N ) is the subspace of functions γ ∈ L 2 (R 3N × R 3N ) such that γ(x, ·) and γ(·, y) lie in L 2 a (R 3N ), for all x, y ∈ R 3N . Furthermore, Γ N 2 = γ N L 2 a×a , where Γ 2 = (Tr L 2 a (R 3N ) (Γ * Γ)) 1/2 .
Let us emphasize that this property only results from the Hilbertian nature of electronic states, without having to take into consideration the (stronger) trace condition Tr(Γ) = 1 requested in their definition. The natural question is then to understand how this condition transfers to operator kernels. Again, basing ourselves on [7], we might expect to obtain the equality: but clearly this cannot make sense based only upon the definition of γ N : by changing its diagonal values it is possible to make the integral appearing in (2) take any real value. On the other hand, it is possible to construct for any Γ N ∈ S(N ) a kernel satisfying (2), for example by writing Γ N as the product of two operators with kernels r, s and considering the kernel γ N of Γ N defined by: In the rest of this work, in view of what has just been said, we shall say that a map is a kernel of Γ N if it is a kernel in the usual sense and satisfies (2). For notational simplicity we choose not to distinguish a state from its kernel : since a kernel and the operator it is associated with are very different objects, the contexts in which one or the other is used will indeed make it clear to which of the two we are referring to.
The density D N (Γ N ) of a state Γ N is defined by which in turn provides a definition of the density D The goal of DFT being to consider electronic systems through their densities rather than their states, a better understanding of what the densities we are working with are is a fundamental question which we shall answer to for the N −centered case after recalling the characterization obtained by Lieb for N −electronic densities [2].
Conversely, consider ρ ∈ I N . The result is yielded by the following alternate expression of ρ : , the ground-state energy of an N −electronic system subject to v is defined by where The N −centered ensemble ground state energy of a system subject to v is then given by which can be rewritten under a variational form : where Due to Corollary (3) and Expression (6), we may carry out the usual "constrained search" method to yield yet another expression of the N −centered energy: where for all density The major advantage coming from (7) is the simplification of the domain over which the search of the infimum is performed. Most notably, the size of the system (the number of electrons N ) has no influence on the complexity of I N (only the integral constraint on admissible densities is changed when considering bigger systems) whereas this is far from being true concerning S {ξ−,ξ+} N . Unsurprisingly, this simplification comes at a cost, "hidden" in the N −centered functional F {ξ−,ξ+} N (ρ) : as in the N −electronic case we do not have any straightforward expression of this functional at our disposal with which to perform the minimization (7) numerically. Interestingly though the core properties of the N −electronic functional transfer to the N −centered one, which is what we shall be discussing in the rest of this section.
For simplicity, we shall be focusing from now on on the so-called "right N -centered" functionals obtained via the particular choice of weights ξ − = 0. In this context, we shall use the notational simplifications Proposition 4 The functional F ξ N is proper, convex and its domain is equal to I N .
Proof Corollary (3) and Definition (8) directly imply that Dom(F ξ N ) ⊂ I N . The equality is obtained by making use of the N − and (N + 1)− electronic universal functionals F N and F N +1 to express F ξ N , namely under the form of an infimal convolution: and by recalling that Dom(F N ) = I N . Since convexity is not guaranteed for the infimal convolution of two convex functions when working in a non-reflexive Banach space, Formula (9) does not yield the convexity result immediately. Fortunately, it can easily be shown by a direct argument: consider t ∈ [0, 1], ρ 1 , ρ 2 ∈ I N and ρ = tρ 1 + (1 − t)ρ 2 . The set inclusion , from which convexity follows.
Proof First note that for all Γ ∈ S ξ N , Thus, it suffices to prove that there exists Γ ∈ S + (N, α) with Γ → ρ such that which prompts the introduction of several notations: Note that the operators h −1 k being bounded, the same is true for h −1 . Without loss we can assume that for all integer n, gn is finite and that g = limn gn exists. By definition of gn, for all integer n there exists a right N −centered state Γn ∈ S ξ N such that: Tr(hΓnh) = Tr(h 2 Γn) ≤ gn + 1 n D ξ N (Γn) = ρn. Setting yn = hΓnh, boundedness of this sequence in the trace norm and the Banach-Alaoglu theorem imply the existence of a subsequence, still written yn, and of some trace-class operator y such that for all compact operator A: lim n Tr(ynA) = Tr(yA).
Positivity of yn implies that of y (consider a Hilbert basis (e i ) of L 2 (R 3N ) ⊕ L 2 (R 3(N +1) ) and set A = ·, e i e i ). Therefore, weak-⋆ lower-semicontinuity of the trace norm implies: lim inf n Tr(yn) ≥ Tr(y).
Setting Γ = h −1 yh −1 , if we prove that Γ is in S ξ N and that Γ → ρ then equation (12) will yield the desired result (10). Since Γn = h −1 ynh −1 is in S ξ N there exists two states Γ N n , Γ N +1 n in S(N ), S(N + 1) such that: The operator sequences y N n := h N Γ N n h N and y N +1 n h N are bounded in the trace norm, so that we may -as before-successively extract subsequences (still written y N n and y N +1 n ) converging to positive trace-class operators y N , y N +1 in the weak-⋆ topology. This implies the convergence of (yn)n = (ω ξ N y N n ⊕ ξy N +1 n )n to (ω ξ N y N ⊕ ξy N +1 ) in the weak-⋆ topology and, by unicity of the limit: . The exact same reasoning as the one proposed in te original proof shows that ρ k n converges to ρ k weakly in L 1 . In particular, the wheighted sum ω ξ N ρ N n + ξρ N +1 n converges to ω ξ N ρ N + ξρ N +1 weakly in L 1 . All that is left to note is that ρn = D ξ N (Γn) = ω ξ N ρ N n + ξρ N +1 n D ξ N (Γ) = ω ξ N ρ N + ξρ N +1 and to recall that ρn converges to ρ weakly in L 1 by hypothesis, which enables us to conclude that D ξ N (Γ) = ρ.
This proposition immediately yields the two fundamental results that follow, the first which will later be seen to justify the fact that minimizing densities in (7) are precisely densities associated to minimizing ensembles in (5), and the second justifying that F ξ N satisfies all the properties required to study it through the lens of convex analysis, as we shall see in the next section.

Theorem 6
The following properties hold true : The sequence (ρ n ) n defined by ρ n := ρ converges to ρ so that, according to Proposition (5), there exists an N −centered state Γ ∈ S ξ N such that D ξ N (Γ) = ρ and:

Duality relations
The usefulness of the properties of the N −centered functional may seem underwhelming at first glance, since they do not directly contribute towards solving the minimization problem of interest (7). As in the classical context, their utility resides elsewhere : we shall see in this Section that they enable us to justify (convex) duality relations linking F ξ N and E ξ N , which will then make it possible to set N −centered DFT in the context of convex analysis.
These duality relations are obtained by making use of tools of convex analysis, regarding the properties of the Legendre transform (for an in-depth introduction one may refer to [8]). A slight formal modification has to be made in order to take into account the fact that the energy and universal functionals are not both convex: one of them is whereas the other is concave.
Definition 7 Consider a Banach space X, a convex functional f on X and a concave functional g on X * .
The Legendre transform f * , respectively double Legendre transform f * * , of f is the functional on X * , respectively X, defined by The concave conjugate f ∧ of f is the functional on X * defined by : The convex conjugate g ∨ of g is the functional on X defined by : Lemma 8 Legendre transforms and concave/convex conjugates are linked by the following relations : The following result clarifies the usefulness of these definitions.

Proposition 9
Energy and universal functionals are set in duality via concave and convex conjugates: Proof The first relation (13) is simply the one obtained via the constrained-search method (7). According to Lemma (8), it follows that: (E ξ N ) ∨ = (F ξ N ) * * . Since F ξ N is convex, proper and lower semicontinuous, it is equal to its double Legendre transform and the second relation (14) is obtained.
With this result at hand, a strong analogy between N −centered DFT and N −electronic DFT can be expected, since most of the results established in the latter case are derived from the duality relations between the universal functional and the energy.
The following proposition establishes the physical meaning of ground state densities, i.e. densities ρ such that Let us write GSD(v) the set of such densities, and GS(v) the set of (N −centered) ground states for v.
Proposition 10 Consider some potential v. There exists a ground state density ρ for this potential if and only if there exists an N −centered ground state Γ for v. Furthermore, the following set equality holds true : Proof Assume that ρ ∈ GSD(v). There exists a state Γ → ρ such that F ξ N (ρ) = Tr(H 0 Γ), which yields : Conversely, consider Γ 0 ∈ GS(v) and let Γ 0 → ρ. Then for all Γ → ρ, . Furthermore, for all ρ ∈ L 1 ∩ L 3 , the following equivalence holds true :

Investigating the weight-dependance of the universal functionals
We have established that for a fixed parameter ξ, N −centered density functional theory can be constructed and treated in complete analogy with N −electronic DFT. In this last paragraph dedicated to these functionals we investigate their weight dependance for a given fixed density ρ. To do so let us define for all ξ ∈ [0, N/(N + 1)]: and set W ρ N (ξ) = +∞ elsewhere. The importance of this parameter-dependent functional is expanded upon in [9], where it is used to express some important features of the electronic system under consideration, namely its electron affinity and ionization potential. More precisely, these expressions are obtained by making use of the derivative of W ρ N . Here we shall be presenting the conclusions we have come to regarding the existence of this derivative.
First and foremost, let us point out that this derivative does indeed exist at all but possibly countably many points of the interval ]0, N/(N + 1)[. This comes from the fact that W ρ N is convex, since by definition it is the point-wise supremum of the family of convex functions Another advantage of this convexity property is that we may study the differentiability of W ρ N through its subdifferential. This leads us to the following proposition.
Further assume that the Hohenberg-Kohn property is satisfied : each v ξ is unique up to an additive constant. Then for all ξ ∈]0, N/(N + 1)[, the two scalars are well defined and Proof First off, let us verify that m ξ N and M ξ N are well defined, i.e that it is independent of the choice of the family v ξ . The proof is straightforward : consider another family −w ξ ∈ ∂F ξ N (ρ). By the Hohenberg-Kohn assumption, for all ξ there exists a constant C ξ such that w ξ = v ξ + C ξ , which in turn yields : Then for all β ∈ [0, N/(N + 1)] and v ∈ L 3/2 + L ∞ : , this can be rewritten under the form : The mindset that lead us to the previous result was to obtain a caracterisation of ∂W ρ N . Since this has not been achieved (we only obtained an inclusion), we might as well omit the Hohenberg-Kohn hypothesis in order to write the following more general result : . For all β > ξ and µ < ξ such that −v β ∈ ∂F β N (ρ) and −vµ ∈ ∂F µ N (ρ) then : Proof The same proof as that of the previous result yields the inclusion ∂W ρ We are then left to note that for all γ :

N −centered generalized Kohn-Sham equations
In this section we shall be focusing on a particular aspect of DFT theory : the Kohn-Sham method. More precisely, our goal is to establish the heuristics leading to N −centered Kohn-Sham equations presented in [9] within the mathematical formalism presented previously, thus identifying that this approach leads to similar mathematical problems as the classical Kohn-Sham method for N −electronic DFT. We will then be taking a quick look at a method proposed in [6] to solve some of these problems, which will further establish the analogy between N −centered and N −electronic DFT.

General idea
The ultimate goal of the Kohn-Sham method is to find an element of GSD(v) for a given v ∈ L ∞ +L 3/2 . As seen previously, this amounts to finding a density ρ such that Before trying to do so, it is necessary to construct a noninteracting (N −centered) DFT. This is done by adapting the previous construction to the situation where electron-electron interactions are not taken into account, leading to the definition of the following noninteracting energy and universal functionals : All results established in the previous section still hold true for these noninteracting functionals, with identical justifications.
The Kohn-Sham approach then amounts to finding a way to transform equation (16) into a noninteracting ground-state density search problem, ideally by finding a functional v s : L 1 ∩ L 3 → L ∞ + L 3/2 such that : Most of the literature on this subject (in the N −electronic case) assumes differentiability of the exchange energy : in order to yield the following Kohn-Sham equation : Unfortunately, the set I N on which the exchange energy is defined has an empty interior in L 1 ∩ L 3 , so that the meaning of the differential notation used in (18) is rather opaque : the exchange energy is obviously not differentiable in the usual sense.

Generalized Kohn-Sham equations
In order to overcome the problem underlined in the previous subsection, a regularization method has been proposed in [6], which leads to a less restrictive Kohn-Sham approach than the one embodied by the equivalence (17) : rather than trying to find a Kohn-Sham potential for a fixed density ρ, the density itself is allowed to be modified. Formally, the goal is now to find two functionals ρ s : L 1 ∩ L 3 → L ∞ + L 3/2 and v s : The reader familiar with the publication in which this method is presented will certainly have no doubt that it can be applied in the N −centered context without adding any theoretical difficulty. Nevertheless, for the sake of completeness and to make things clear for readers discovering this subject, we choose to expose the details of this approach for N −centered DFT. The major prerequisite necessary to understand it is a convex analysis tool : the Moreau-Yosida regularization. We recall its definition along with the fundamental properties to which we shall be appealing and refer to [10] for details. Definition 14 Let (H, ·, · ) be a Hilbert space, f a functional on H and ε > 0. The ε−Moreau-Yosida regularization M Yε(f ) of f is the functional on H defined by : Proposition 15 Let f ∈ Γ 0 (H), ε > 0. Then : 1. M Y ε (f ) is convex, real valued and exact. (18) is uniquely attained at some point Prox εf (ρ).

The infimum in
One might wonder how this tool can be of any use for DFT since it features functionals on a Hilbert space when the universal functionals are defined on the non-reflexive Banach space L 1 ∩ L 3 . This indeed needs to be discussed : even though the Moreau-Yosida transform can be defined on reflexive Banach space without sacrificing much of its properties, it doesn't fare as well in the non-reflexive context [11]. A first step towards working in a reflexive context is made by recalling the inclusion L 1 ∩ L 3 ⊂ L 2 , which combined with (7) yields : Unfortunately, this inclusion does not imply that Γ 0 (L 1 ∩ L 3 ) ⊂ Γ 0 (L 2 ) (quite the contrary), so that the duality linking the N −centered energy and functional cannot be guaranteed to hold in L 2 : we cannot directly assert that so that it becomes necessary to define a new "universal functional" in duality with the energy on L 2 (R 3 ). This is made possible by the continuity of E ξ N on L 2 (R 3 ), leading to : . Furthermore, for all ρ ∈ I N the following equivalence holds : Proof The first statement is a direct consequence of the inclusion . We may then write : Conversely, assume that −v ∈ ∂F ξ N (ρ) and F ξ N (ρ) = F ξ N (ρ). We immediately see that v ∈ L 2 (R 3 ). It then follows that for all n ∈ L 1 ∩ L 3 : One of the motivations behind N −centered density functional theory as opposed to ensemble DFT is the fact that the domains of the universal functionals F ξ N do not depend on the choice of the parameter ξ. The following result shows that this still holds true for the functionals F Proof From (22) and lemma 16 we get : Thus, letting G ξ N be the functional equal to F ξ N on I N and +∞ elsewhere, the functionals G ξ N and F ξ N are both in Γ 0 (L 2 (R 3 )) and have same concave conjugate. This implies that these functionals are equal and yields the result.
Of course it is again possible to adapt the construction of F ξ N to the noninteracting case, yielding a noninteracting functional F ξ N,∅ in duality with the noninteracting N −centered energy in L 2 . Let us now show that the method proposed in [6] in order to link the subdifferentials of F ξ N and F ξ N,∅ can be carried out in the N −centered context with no significant modifications.

Proposition 18 (Generalized Kohn-Sham equations for N −centered DFT)
The following equivalence holds for all ε > 0 : The key idea is to make use of the properties of the Moreau-Yosida transform in order to construct the path linking the noninteracting and interacting systems. This being said, the result can be obtained by a simple series of equivalent statements :

Remark 19
The resolution of this generalized Kohn-Sham equation can, in theory, be handled in the same way as the usual one by setting up fixed point iterations : If the convergence of this scheme for the usual Kohn-Sham method was already a problem to be considered, having set up this generalised scheme does not simplify matters. Indeed, it is now necessary to evaluate a regularized exchange potential, with no clear idea of the physical meaning of this object, which itself is defined through a generalised universal functional itself carrying a less clear physical sense than the true functional F ξ N .

Semi-classical limits for N −centered DFT
This far, all properties of conventional DFT we have investigated have been shown to hold in the N −centered DFT setting and the mathematical problems arising to be the same (existence of a Kohn-Sham configuration, applicability of the Hellmann-Feynmann theorem). In order to expand on this we shall study in this Section the behavior of another result, established on firm mathematical grounds for N −electronic DFT, which goes by the name of "semiclassical limit of the universal functional", when written in the N −centered DFT context.

The
Let us begin by recalling some results established in [12] regarding this semiclassical limit.
Definition 20 For all N and ρ ∈ I N we can define the Coulomb multi-marginal problem : Theorem 21 For all ρ ∈ I N : Thereafter, we write : Theorem 22 Consider χ ∈ C ∞ c (R 3 , R) a radial function with support in the unit ball such that χ 2 = 1 and P a symmetric probability measure on R 3N with support in D c N,α := (x 1 , . . . , x N )|x i − x j | ≥ α, ∀i = j , for some α > 0 and such that For all α/4 > ε > 0 there exists Γε ∈ S(N ) such that : and for all symmetric Φ ∈ C 2 (D c α−4ε ) : Let us now go back to the N −centered case. For simplicity we limit ourselves to the left N −centered case. Let us define ). The proof of the following result is based on the following statement : for all Proposition 23 For all ρ ∈ I N , Proof For all η write : It follows that ≥ E N,ξ OT (ρ). Conversely, consider ρ N ∈ I N and ρ N +1 ∈ I N +1 such that ρ = ω ξ N ρ N + ξρ N +1 . Let P N (respectively P N +1 ) be a minimizer for E N OT (ρ N ) (respectively for E N +1 OT (ρ N +1 )) and Γ N ε (respectively Γ N +1 ε ) the state yielded by theorem 22 for P N (respectively P N +1 ). There exists β, δ such that P N has its support in D c N,β , P N +1 in D c N +1,δ . Defining Γε = (1 − α)Γ N ε + α + Γ N +1 ε we may then write : The properties of Γ N ε and Γ N +1 ε lead to : Thus, taking η small enough and ε = η 1/4 , the previous inequalities yield : , and letting η → 0 we obtain : This being true for all ρ N , ρ N +1 such that ρ = ω ξ N ρ N + ξρ N +1 , the result follows.

Remark 24
The result proven there-before is weaker than the one shown by Lewin, which stated the existence of a constant C(ρ) such that If such a result holds in the N −centered context a further analysis of C(ρ N , ρ N +1 ) shall be necessary to yield it, in particular it would suffice to findC such that Regarding the analogy between N −centered and N −electronic DFT, the result obtained in the previous proposition is not entirely satisfying since the N −centered semiclassical limit does not explicitly appear to be a standard optimal transport problem.

An upper bound of E N +1
OT which is an N −electronic optimal transport problem In this subsection we present the construction of an upper bound for E N,ξ OT , which we derived when trying to express this functional as an optimal transport problem. The idea we were investigating was that of expressing the (N + 1)electronic semiclassical limit as an N −dimensional optimal transport problem, hoping to make the link between E N OT and E N +1 OT more clear and uncover an alternative expression of E N,ξ OT . The initial thought was to construct a bijection between the probability sets N (ρ/N ) and N +1 (ρ/N ) for a given density ρ ∈ I N . This was quickly set aside since it can already be seen to be false when N equals one (in which case the first set is a singleton). On the other hand, it should be possible to define injections/surjections between these sets, which could then possibly yield inequalities linking E N OT and E N +1 OT . Let us begin by investigating the construction of an injective mapping, while trying to guarantee "as much surjectivity as possible". Defining the probability measure µ ρ (A) = A ρ(x)dx for borelian sets A ∈ B(R 3 ), an obvious way to map a probability P in N (ρ/N ) to a probability measure on R 3(N +1) whose marginals are all equal to ρ/N is to define Of course this is far from satisfying since λ(P) does not belong to N +1 (ρ/N ) in general : the only case where λ(P) would be symmetric is if P were given by P = N i=1 µρ N . This can be remedied by defining the map Λ obtained by mixing probability measures similar to λ(P) : Proposition 25 The map Λ : Π N (ρ/N ) → Π N +1 (ρ/N ) is injective.
Proof Let us consider P 1 , P 2 ∈ Π N (ρ/N ) such that Λ(P 1 ) = Λ(P 2 ). For integers 1 ≤ k ≤ N we propose a recursive proof of the statement : Since P 1 and P 2 are both in Π N (ρ/N ), (H 1 ) is true (they have same marginals). Let us assume that (H k ) holds true for some k < N and write : Doing the same for P 2 it follows that : which yields (H k+1 ).

This result leads to an upper bound for E N +1
OT taking the form of an N −dimentional optimal transport problem, namely : which the following technical lemma will help make more explicit.
Lemma 26 For all symmetric measurable function f , Proof The definition of Λ(P) and the symmetry of f and P yield : Corollary 27 For all ρ ∈ I N , We may further precise the right hand side of (24) by noting that it is equal to the following expression : We hope to sharpen this estimation in future work in order to get closer to a more standard optimal transport formulation of the semiclassical limit in the N -centered setting.

Application to the Hubbard Dimer
In this Section our goal is to provide a self-contained presentation of the Hubbard dimer, a toy model used by chemists in order to test and provide proof-of-concept for DFT-based calculations with at its core an extremely simplified Hamiltonian, and to apply this model in the N -centered setting. In order to fully understand the meaning of this operator, which we write down explicitly in formula (26) there-below, a prerequisite concept we need to introduce is that of the second quantization. For further detail on the subject one may refer to [13].

CAR algebra and second quantization
In this section, (H, ·, · L 2 (R 3 ) ) ⊂ L 2 (R 3 ) is a Hilbert space and {f i } a basis of H. Let us begin by presenting a brief overview of the CAR relations and second quantization, two tools necessary to understand the definition of the Hubbard Dimer which are discussed for example in [12].
Definition 30 The N-particle Fock space associated to H is the Hilbert space : The Fock space associated to H is The vacuum state is defined by Ω = (1, 0, . . . ) ∈ Fa.
Definition 31 For all f ∈ H, the creation operator The annihilation operator a(f ) is the adjoint of a † (f ) and can be explicitly defined by : Annihilation and creation operators satisfy the canonical anticommutation relations (CAR) : Let A be a self-adjoint, bounded from below operator on H of domain can be re-written in terms of creation and annihilation operators : Similarly, given an operator W acting on H 2 a , the operator can be rewritten : In particular, the hamiltonian H N (v) can be written :

The Hubbard dimer and its functionals
Based on expression 25, a toy Hamiltonian known as the Hubbard Dimer has been proposed in order to, amongst other things, apply DFT methods in a simplified context. The main assumptions leading to this model are the following : rather than considering an orthonormal family for the whole space L 2 (R 3N ) we limit ourselves to four sites {f 1↑ , f 1↓ , f 2↑ , f 2↓ } (numbers correspond to "spatial" position and arrows to spin direction) and assume that inter-electronic interaction between two spatial sites are negligible. This leads to the study of the operator : Proposition 32 The matricial representations of the operators HD k (v 1 , v 2 ) in the basis span(f i1 ∧ · · · ∧ f i k ) i1<···<i k are given by : In the context of the Hubbard dimer, the universal functionals are defined for n = (n 1 , n 2 ) ∈ R 2 by : . Explicit expressions of these eigenvalues are known and are recalled in the appendix. From these analytical expressions and an asymptotic study we derive several results relative to the functionals, gathered in the following result.

Proof See appendix.
Corollary 34 The domain of the right 2−centered electronic functional The functional n → f ξ 2 (n, 2 − n) for various choices of ξ.
Proof Simply note that f ξ 2 can be written as an infimal convolution of dilatations of f 2 and f 3 : This result shows that, in opposition to physical N −centered DFT, N −centered DFT applied to the Hubbard dimer does not yield a functional who's domain is independent of α. A graphical representation of this assertion can be found in figure 1, where we plotted the functional n → f ξ 2 (n, 2 − n) for various weights. As such, this model does not seem adequate in order to illustrate N −centered DFT methods and some alternatives should be considered.
Before completely dismissing the Hubbard model, let us note the following : even though the domain of the 2−centered functional does indeed depend on the weight, this dependance does not lie in the "integral condition" n 1 + n 2 = 2 (analog to ρ = 2) but rather seems to come from the fact that the 3−electronic functional itself diverges from the physical model : the condition n ≥ 1 on the occupation numbers for this functional's domain has no analog for physical 3-electronic densities I 3 . As such, the problem we are faced with seems to come not from the fact that we are trying to set up N −centered DFT, but rather that the domains of N −electronic density functionals for the Hubbard model (dimer or more generally considering K sites) themselves carry nonphysical constraints when K < N, due to the fact that in this situation the discrete nature of this model all sites have to be occupied by at least one electron.
As such, the Hubbard model will correctly represent the domains of N −centered universal functionals in cases where all sites are free to be occupied by 0, 1 or 2 electrons, in other words for Hubbard models with a number of sites K ≥ N. This is illustrated for the 1-centered functional associated to the 2-site model in figure 2. Figure 2 The functional n → f ξ 1 (n, 1 − n) for various choices of ξ.

Conclusion
The construction and study of the N -centered universal functional presented in this work provides a clear mathematical formalism for the study of this theory, which can easily be adapted to other ensemble states with different choices of weights, such as the states used in PPLB DFT for example. Thanks to this formalism we were able to show that the Kohn-Sham approach can be obtained in an identical fashion as in the N -electronic context, whereas other branches of the classical theory cannot immediately be adapted to the N -centered context, as is the case for the semiclassical limit or some specific Hubbard dimer models.
Amongst the theoretical problems related to ensemble DFT, other than the ones that already appear for single state DFT, the understanding of the weightdependance of the universal functional is a major issue of which we proposed an interpretation through the lens of convex analysis, which led us to precise the link between weight-dependance and electronic affinity and that we hope to investigate further in future work.
Let us first note that the columns of HD 2 (v 1 , v 2 ) containing only v 1 + v 2 do not impact on ǫ 2 .
Lemma 35 For all v 1 , v 2 ∈ R, ǫ 2 (v 1 , v 2 ) is the lowest eigenvalue of Proof It suffices to prove that the lowest eigenvalue of HD r 2 is smaller than v 1 + v 2 (which is the only eigenvalue of HD 2 that may not be in the spectrum of HD 2 ), i.e. that : This is easily seen by considering x = (0, 1, 1, 0) for which we obtain : Lemma 36 For all v 1 , v 2 ∈ R, With the previous lemma now in hand, we recall the expressions of E k (δv) = ǫ k (δv/2, −δv/2) (see for example [15] to understand where these come from). We limit ourselves in what follows to the particular case U = U 1 = U 2 .
With the analytic expressions of ǫ k now at our disposal, we can prove proposition (33).