# Mathematical approach to N -centered ensemble DFT and associated methods for the study of open systems

Guillaume Grente, Zakaria Belhachmi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Guillaume Grente, Zakaria Belhachmi. Mathematical approach to N -centered ensemble DFT and associated methods for the study of open systems. 2023. hal-04512728

HAL Id: hal-04512728

https://hal.science/hal-04512728v1

Preprint submitted on 20 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mathematical approach to N-centered ensemble
DFT and associated methods for the study of
open systems

Guillaume Grente* Zakaria Belhachmi†

February 22, 2024

Abstract

A fundamental aspect of the study of N-electronic systems is to ob-
tain 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 developed, the one around which
this work revolves being to consider electronic systems through their elec-
tronic density rather than their state. A first implementation of this idea
was proposed by Thomas [21] and Fermi [1] and was later revisited by
Hohenberg and Kohn [12] who set the basis for the method known today
as Density Functional Theory (DFT) for which an in-depth mathematical
study has been established by Lieb [16]. Since then, the ideas leading to
the construction of DFT have been adapted to the context of electronic
systems with a fractional number of electrons (open systems), first through
PPLB DFT [17] and more recently through the definition of N-centered
ensemble DFT [19, 20]. In both cases it is unclear whether the mathe-
matical properties established for regular 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 ensemble and regular 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 ensemble
DFT, investigate the links between this theory and optimal transport and
present the Hubbard Dimer in this particular situation.

**Keywords.** Density Functional Theory, Convex Analysis, Open Systems,

Semiclassical limit, optimal transport

<footnote>* ggrente@unistra.fr, IRMA, Univ. Strasbourg and IRIMAS, Univ. Haute Alsace</footnote>

<footnote>†zakaria.belhachmi@uha.fr, IRIMAS, Univ. Haute Alsace</footnote>
<table><thead><tr><th>Notation</th><th>Description</th><th>Page List</th></tr></thead><tbody><tr><td>S<sub>N</sub></td><td>set of N-electronic states.</td><td>4</td></tr><tr><td>D<sub>N</sub>(Γ)</td><td>Density of an N-electronic state Γ.</td><td>5</td></tr><tr><td>D<sup>Ξ</sup>(Γ)</td><td>Density of an N-centered ensemble Γ.</td><td>5</td></tr><tr><td>I<sub>N</sub></td><td>{ρ ∈ L<sup>1</sup>(R<sup>3</sup>) ∩ L<sup>3</sup>(R<sup>3</sup>) | ∇ρ<sup>1/2</sup> ∈ L<sup>2</sup>, ρ ≥ 0, ||ρ||<sub>L<sup>1</sup></sub> = N}.</td><td>5</td></tr><tr><td>E<sub>N</sub>(v)</td><td>N-electronic ground-state energy for a potential v.</td><td>6</td></tr><tr><td>H<sup>N</sup><sub>0</sub></td><td>-∑<sub>i=1</sub><sup>N</sup> Δ<sub>i</sub> + ∑<sub>1≤i<j≤N</sub> <sup>1</sup>&frasl;<sub>|x<sub>i</sub>-x<sub>j</sub>|</sub>.</td><td>6</td></tr><tr><td>E<sup>Ξ</sup><sub>N</sub>(v)</td><td>N-centered ensemble ground-state energy.</td><td>6</td></tr><tr><td>H<sub>0</sub></td><td>⊕<sub>k≥1</sub> H<sup>k</sup><sub>0</sub>.</td><td>6</td></tr><tr><td>F<sup>Ξ</sup><sub>N</sub></td><td>N-centered ensemble universal functional.</td><td>7</td></tr><tr><td>F<sup>ξ</sup><sub>N</sub></td><td>F<sub>N</sub><sup>{0,ξ}</sup>.</td><td>7</td></tr><tr><td>w<sup>ξ</sup><sub>N</sub></td><td>w<sub>N</sub><sup>{0,ξ}</sup>.</td><td>7</td></tr><tr><td>lim inf<sub>n→∞</sub> f(x<sub>n</sub>)</td><td>lim<sub>n→∞</sub> (inf<sub>m≥n</sub> f(x<sub>m</sub>)).</td><td>9</td></tr><tr><td>f<sup>*</sup></td><td>Legendre transform of f.</td><td>11</td></tr><tr><td>f<sup>^</sup></td><td>Concave conjugate of f.</td><td>11</td></tr><tr><td>g<sup>√</sup></td><td>Convex conjugate of g.</td><td>11</td></tr><tr><td>E<sup>ξ</sup><sub>N,0</sub></td><td>Noninteracting right N-centered ground-state energy.</td><td>16</td></tr><tr><td>F<sup>ξ</sup><sub>N,0</sub></td><td>Noninteracting right N-centered universal functional.</td><td>16</td></tr><tr><td>E<sup>ξ</sup><sub>xc,N</sub></td><td>F<sup>ξ</sup><sub>N</sub> - F<sup>ξ</sup><sub>N,0</sub>.</td><td>16</td></tr><tr><td>MY<sub>ε</sub>(f)</td><td>ε-Moreau-Yosida regularization of f.</td><td>17</td></tr></tbody></table>
# Introduction

By extending the Hohenberg-Kohn approach to ensemble states with fractional electron number Perdew, Parr, Levy and Balduz [17] defined the standard DFT framework for the description of electron ionization and affinity processes, referred to as PPLB DFT. In this context charged excitations are accessed through a continuous variation of the electron number $N$ which, by construction of this theory, is entirely defined by the ensemble density which integrates to $N$. This implies that in an effort to construct model density functionals compatible with this approach one has to simultaneously take into account the influence of the variation of the density and of the electron number, which can lead to major difficulties. An alternative framework for the use of DFT methods for the description of charged excitations, $N$-centered ensemble DFT, was recently proposed by Senjean and Fromager [19, 20] and will be referred to as $N$-centered eDFT. This method contrasts with the PPLB approach in the sense that it is formulated in such a way that the fractional part $\xi$ of the electron number $N = N+\xi$ and the ensemble density can be considered as independent variables, making it possible to vary them separately.

In both frameworks, systems with fractional electron number $N = N + \xi$ are represented by weighted ensembles of the form $w_N^\xi \Gamma^N \oplus w_{N+1}^\xi \Gamma^{N+1}$ where $\Gamma^k$ are $k$-electronic states and:

* $(w_N^\xi, w_{N+1}^\xi) = (\xi, 1-\xi)$ for PPLB DFT.

* $N w_N^\xi + (N+1) w_{N+1}^\xi = N$ for $N$-centered DFT (so that the ensemble density integrates to $N$).

For a fixed choice of weights, the ensemble energy of an ensemble $w_N^\xi \Gamma^N \oplus w_{N+1}^\xi \Gamma^{N+1}$ is expressed as a linear combination of the energies of the integer-electronic states $\Gamma^N$ and $\Gamma^{N+1}$:

$$w_N^\xi E^N(\Gamma^N) + w_{N+1}^\xi E^{N+1}(\Gamma^{N+1}). \qquad (1)$$

The associated electronic problem (i.e. the search of the ensemble ground-state) in this context is the minimization of expression (1) over all ensembles $w_N \Gamma^N \oplus w_{N+1} \Gamma^{N+1}$. Although the usual constrained search method may be set up in order to reformulate this problem in a DFT setting by defining a universal functional of the density it is unclear whether a strong analogy can be expected to hold with regular DFT: having modified the set of states defining the electronic problem did we modify the nature of the universal functional that is obtained? Such questions are the main motivation of this work, in which we aim to clarify the mathematical formulation of $N$-centered DFT in order to precise its similarities and discrepancies with the regular theory.

Let us emphasize that, even though we chose to focus exclusively on $N$-centered DFT, ensemble DFT as a whole is an active research domain [9, 6, 8]. In this regard we wish to insist on the fact that all the results established in this work which do not involve the variation of the electron number do not rely on
any arguments that would not hold true for PPLB DFT or any other ensemble DFT method for charged excitations defined based on a choice of ensembles of the form $w_N^\xi \Gamma^N \oplus w_{N+1}^{\xi} \Gamma^{N+1}$. On the other hand the approach we propose cannot be adapted as such to ensemble DFT methods for neutral excitations, such as GOK DFT [10], which require additional orthogonality constraints on the states used to define the ensembles under consideration.

The paper is organized as follows. In Section 1 we justify that the work done by Lieb in order to propose a precise definition of DFT [16] can still be followed to obtain a similar construction of N-centered eDFT. 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 [13] 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 expected 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 eDFT.

# 1 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.

## 1.1 N-centered ensembles and densities

An N-electronic state is a self-adjoint trace class linear operator on $L_a^2(\mathbb{R}^{3N})$ (the Hilbert space of square integrable antisymmetric functions of N spatial variables, endowed with the usual $L^2$ scalar product) such that $\mathrm{Tr}(\Gamma) = 1$ and

$$ \mathrm{Tr}((-\Delta)^{1/2}\Gamma(-\Delta)^{1/2}) < \infty. $$

The set of N-electronic states is written $S_N$.

The main idea of N-centered DFT is the following: rather than working with the set of ensemble states of the form $(1-\xi)S_N \oplus \xi S_{N+1}$ that are used in PPLB DFT, we wish to consider non-physical ensembles (whose traces are not equal to one) depending on a certain number of parameters that can vary independently of the electron number which will remained fixed, equal to an integer N, for any parameter choice. We begin by providing a general formulation of this idea
and shall then for clarity focus on the study of a particular setting that was considered in the initial formulation of this theory [19, 20].

Begin by choosing the number $p$ of parameters and the finite subset $J \subset \mathbb{N}$ of integer-electronic components we wish to consider. These parameters will impact the $k$-electronic component of the ensemble considered in a way that is predetermined by the choice of real-valued positive maps $(w_k)_{k \in J}$, defined on a subset $\Omega \subset \mathbb{R}^p$, satisfying the condition $\sum_{k \in J} kw_k = N$. For a given parameter choice $\Xi \in \Omega$, the associated set of $N$-centered ensembles $S_N^\Xi$ is then defined by:

$$S_N^\Xi = \bigoplus_{k \in J} w_k(\Xi)S_k. \qquad (2)$$

For a given kernel $\gamma$ of an $N$-electronic state $\Gamma$ satisfying the integral condition $\int_{\mathbb{R}^{3N}} \gamma(x, x) = 1$ we define its density $D_N(\Gamma)$ by

$$D_N(\Gamma)(x) = N \int_{\mathbb{R}^{3(N-1)}} \gamma(x, x_2, \dots, x_N; x, x_2, \dots, x_N) dx_2 \dots dx_N,$$

which in turn provides a definition of the density $D^\Xi(\Gamma)$ of an $N$-centered ensemble $\Gamma = \bigoplus_{k \in J} w_k(\Xi)\Gamma_k$:

$$D^\Xi(\Gamma) = \sum_{k \in J} w_k(\Xi)D_k(\Gamma_k). \qquad (3)$$

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 [16].

**Proposition 1** For all integer $N$, the range of the density map $D_N$ is the set

$$\mathcal{I}_N = \{ \rho \in L^1(\mathbb{R}^3) \cap L^3(\mathbb{R}^3) | \nabla \rho^{1/2} \in L^2, \rho \ge 0, \|\rho\|_{L^1} = N \}.$$

**Corollary 2** For all integer $N$, any choice of weight maps and parameters $\Xi \in \Omega$, the range of $D^\Xi$ is the set of $N$-electronic densities $\mathcal{I}_N$.

**Proof.** Making use of (3) and proposition (1), this amounts to show that:

$$\mathcal{I}_N = \sum_{k \in J} w_k(\Xi)\mathcal{I}_k. \qquad (4)$$

Consider a density $\rho = \sum_k w_k(\Xi)\rho_k$ in the right-hand-side of (4). The density $\rho$ is in $L^1(\mathbb{R}^3) \cap L^3(\mathbb{R}^3)$ and since the weight maps $w_k$ are positive it is also positive. The convexity and homogeneity of $G(\rho) \mapsto \int (\nabla \rho^{1/2})^2$ yield:

$$G(\rho) \le \sum_{k \in J} w_k(\Xi)G(\rho_k) < +\infty,$$
and the last condition is derived by recalling that the weight maps were chosen so that $\sum_k kw_k(\Xi) = N$.

Conversely, consider $\rho \in \mathcal{I}_N$. The result is yielded by the following alternate expression of $\rho$:

$$ \rho = \sum_{k \in J} w_k(\Xi) \frac{k}{N} \rho, $$

after noting that the density $k\rho/N$ is in $\mathcal{I}_k$. $\square$

## 1.2 N-centered universal functionals

Given a potential $v \in L^{3/2}(\mathbb{R}^3) + L^\infty(\mathbb{R}^3)$, the ground-state energy $E_N(v)$ of an N-electronic system subject to v is defined by

$$ E_N(v) = \inf \left\{ \mathrm{Tr}(H_0^N \Gamma) + \int \mathcal{D}_N(\Gamma) v |\Gamma \in \mathcal{S}_N \right\}, \quad (5) $$

where $H_0^N$ is the operator

$$ H_0^N = - \sum_{i=1}^{N} \Delta_i + \sum_{1 \le i < j \le N} \frac{1}{|x_i - x_j|}. $$

The N-centered ensemble ground state energy $E_N^{\Xi}(v)$ of a system subject to v is then given by

$$ E_N^{\Xi}(v) = \sum_{k \in J} w_k(\Xi) E_k(v) $$

which can be rewritten under a variational form:

$$ E_N^{\Xi}(v) = \inf \left\{ \mathrm{Tr}(H_0 \Gamma) + \int \mathcal{D}^{\Xi}(\Gamma) v |\Gamma \in \mathcal{S}_N^{\Xi} \right\}, \quad (6) $$

where $H_0$ is the operator

$$ H_0 = \bigoplus_{k \ge 1} H_0^k. $$

Due to Corollary (2) and Expression (6), we may carry out the usual “constrained search” method to yield yet another expression of the N-centered energy:

$$ E_N^{\Xi}(v) = \inf \left\{ F_N^{\Xi}(\rho) + \int v \rho |\rho \in \mathcal{I}_N \right\}, \quad (7) $$

where for all density $\rho \in L^1(\mathbb{R}^3) \cap L^3(\mathbb{R}^3)$:

$$ F_N^{\Xi}(\rho) = \inf \{ \mathrm{Tr}(H_0 \Gamma) | \Gamma \in \mathcal{S}_N^{\Xi}, \mathcal{D}^{\Xi}(\Gamma) = \rho \}. \quad (8) $$

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 $\mathcal{I}_N$ (only the integral constraint on admissible densities is changed when considering bigger systems) whereas this is far from being true concerning $S_{\Xi_N}$. Unsurprisingly, this simplification comes at a cost, “hidden” in the N-centered functional $F_{\Xi_N}^\xi$: 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 for the particular choice of weight $\xi_-=0$. From now on the functional $F_N^{\{0,\xi\}}$ and the weight $w_N^{\{0,\xi\}}$ are written $F_N^\xi$ and $w_N^\xi$ and other notations relative to N-centered ensembles are adapted in a similar fashion.

In order to study the properties of the functional $F_N^\xi$ one may start by wondering if an explicit link can be found between this functional and the integer-electronic ones $F_N$ and $F_{N+1}$. The following result establishes that such a link does indeed exist: $F_N^\xi$ is shown to take the form of an infimal convolution between $F_N$ and $F_{N+1}$.

**Proposition 3** For all $\rho \in L^1(\mathbb{R}^3) \cap L^2(\mathbb{R}^3)$ the following equality holds true:

$$F_{\Xi_N}^{\Xi}(\rho) = \inf \left\{ \sum_{k \in J} w_k(\Xi) F_k \left( \frac{\rho_k}{w_k(\Xi)} \right) \mid \rho_k \in L^1(\mathbb{R}^3) \cap L^3(\mathbb{R}^3), \sum_{k \in J} \rho_k = \rho \right\}, \quad (9)$$

**Proof.** By separating the minimization defining it (8), the universal functional $F_{\Xi_N}^\xi$ reads:

$$F_{\Xi_N}^{\Xi}(\rho) = \inf \left\{ \sum_{k \in J} w_k(\Xi) F_k(\rho_k) \mid \rho_k \in \mathcal{I}_k, \sum_{k \in J} w_k(\Xi) \rho_k = \rho \right\}, \quad (10)$$

Noting that the constraint $\rho_k \in \mathcal{I}_k$ can be replaced by $\rho_k \in L^1 \cap L^3$ with no impact on this minimization and after replacing $\rho_k$ by $\rho_k/w_k(\Xi)$ equation (10) then yields the result.

□

**Remark 4** By strict analogy and setting $w_{N+1}^+(\xi) = w_{N+1}(\xi, 0) = 1 - \xi \frac{N}{N+1}$ the following holds for the left $(N+1)$-centered functional:

$$F_{N+1}^{\{\xi,0\}}(\rho) = \inf_{n \in L^1 \cap L^3} \left\{ \xi F_N \left( \frac{\rho - n}{\xi} \right) + w_{N+1}^+(\xi) F_{N+1} \left( \frac{n}{w_{N+1}^-(\xi)} \right) \right\}. \quad (11)$$

Since for all $\alpha \in [0, 1]$ the left and right centered weights are related by

$$\left( \frac{N+1}{N} \alpha, w_{N+1}^{+} \left( \frac{N+1}{N} \alpha \right) \right) = \frac{N+1}{N} \left( w_{N}^{\frac{N}{N+1}(1-\alpha)}, \frac{N}{N+1}(1-\alpha) \right),$$
we may note that equations (9) and (11) enable the recovery of the following link between left and right centered universal functionals, which was previously established in [20]:

$$F_{N+1}^{\{\frac{N+1}{N}\alpha,0\}}(\rho) = \frac{N+1}{N} F_N^{\{0, \frac{N}{N+1}(1-\alpha)\}}\left(\frac{N}{N+1}\rho\right).$$

The following proposition collects some properties of the integer-electronic functionals that may readily be seen to transfer to $F_N^\xi$ thanks to the expression (9) established in the previous lemma.

**Proposition 5** Independently of the choice of the weight $\xi$, the functional $F_N^\xi$ is proper, convex and its domain $\text{Dom}(F_N^\xi)$ is equal to $\mathcal{I}_N$.

**Proof.** Since $F_N^\xi$ is the infimal convolution (for more detail on this operation one may refer to Chapter 12 of [4]) of the two proper convex functionals

$$f_N^\xi : n \mapsto w_N^\xi F_N\left(\frac{n}{w_N^\xi}\right)$$

and $f_{N+1}^\xi : n \mapsto \xi F_{N+1}\left(\frac{n}{\xi}\right)$,

it is proper, convex and satisfies:

$$\text{Dom}(F_N^\xi) = \text{Dom}\left(f_N^\xi\right) + \text{Dom}\left(f_{N+1}^\xi\right) = w_N^\xi \mathcal{I}_N + \xi \mathcal{I}_{N+1} = \mathcal{I}_N,$$

where the last equality follows from 4. □

**Remark 6** Convexity of $F_N^\xi$ can also be established by a direct argument: consider $t \in [0, 1]$, $\rho_1, \rho_2 \in \mathcal{I}_N$ and $\rho = t\rho_1 + (1-t)\rho_2$. The set inclusion

$$\{\Gamma = t\Gamma_1 + (1-t)\Gamma_2 | \Gamma_k \in S_N^\xi, D_N^\xi(\Gamma_k) = \rho_k, k \in \{1, 2\}\} \subset \{\Gamma \in S_N^\xi | D_N^\xi(\Gamma) = \rho\}$$

yields the inequality:

$$F_N^\xi(\rho) \le \inf \left\{ t \operatorname{Tr}(H_0 \Gamma_1) + (1-t) \operatorname{Tr}(H_0 \Gamma_2) \middle| \Gamma_k \in S_N^\xi, D_N^\xi(\Gamma_k) = \rho_k, k \in \{1, 2\} \right\},$$

from which convexity follows.

In addition to convexity of the functional under study most results established in the framework of convex analysis require an additional property: lower semi-continuity, which we recall there-below.

**Definition 7** Given a Banach space $X$, a functional $f: X \to \mathbb{R} \cup \{+\infty\}$ is lower semi-continuous if for all $x \in X$ and any sequence $(x_n)$ converging to $x$ the following inequality holds:

$$f(x_0) \le \liminf_{n \to \infty} f(x_n), \qquad (12)$$
where $\liminf_{n \to \infty} f(x_n)$ is the standard notation for the limit inferior of the sequence $((f(x_n)))$:

$$ \liminf_{n \to \infty} f(x_n) = \lim_{n \to \infty} \left( \inf_{m \ge n} f(x_m) \right). \qquad (13) $$

Unfortunately even though $F_N$ and $F_{N+1}$ satisfy this property the expression (9) does not guarantee that it will transfer to $F_{N+1}^\xi$ since the Banach space $L^1(\mathbb{R}^3) \cap L^3(\mathbb{R}^3)$ in which this discussion takes place is non-reflexive. Following the procedure proposed in [16], in order to prove the lower semi-continuity of the functional $F_N^\xi$ we establish the following, slightly stronger, result. The proof can be found in Appendix A.

**Proposition 8** Consider a density $\rho \in \mathcal{I}_N$ and any sequence $(\rho_n)$ that converges to $\rho$ weakly in $L^1$. There exists an ensemble state $\Gamma \in S_N^\xi$ such that

$$ \begin{aligned} D_N^\xi(\Gamma) &= \rho \\ \mathrm{Tr}(H_0\Gamma) &\le \liminf_{n \to +\infty} F_N^\xi(\rho_n). \end{aligned} $$

This proposition 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^\xi$ satisfies all the properties required to study it through the lens of convex analysis, as we shall see in the next section.

**Theorem 9** The following properties hold true:

1. The infimum defining the functional $F_N^\xi$ (8) is a minimum.

2. The functional $F_N^\xi$ is lower semicontinuous.

**Proof.**

1. Let $\rho$ be a density in $\mathcal{I}_N$. The sequence $(\rho_n)_n$ defined by $\rho_n := \rho$ converges to $\rho$ so that, according to Proposition (8), there exists an ensemble $\Gamma \in S_N^\xi$ such that $D_N^\xi(\Gamma) = \rho$ and:

$$ \mathrm{Tr}(H_0\Gamma) \le F_N^\xi(\rho) = \inf \{ \mathrm{Tr}(H_0\Gamma) | \Gamma \in S_N^\xi, D_N^\xi(\Gamma) = \rho \} $$

which implies the result.

2. Consider a sequence $(\rho_n)$ that converges to $\rho$ in $L^1(\mathbb{R}^3) \cap L^3(\mathbb{R}^3)$. Then $(\rho_n)$ converges to $\rho$ weakly in $L^1(\mathbb{R}^3)$. It follows from proposition (8) that there exists $\Gamma \in S_N^\xi$ such that

$$ D_N^\xi(\Gamma) = \rho \qquad (14) $$

$$ \mathrm{Tr}(H_0\Gamma) \le \liminf_{n \to +\infty} F_N^\xi(\rho_n). \qquad (15) $$
The definition of $F_N^\xi$ and equation (14) lead to the inequality

$$F_N^\xi(\rho) \le \mathrm{Tr}(H_0\Gamma)$$

and equation (15) then yields the semicontinuity result.

□

## 1.3 N-centered DFT as a state space theory

The grand canonical setting generalizes the notion of states by defining them as being any nonnegative trace class self-adjoint operator $\Gamma$ on the Fock space

$$\mathcal{F}_a = \mathbb{C} \oplus \bigoplus_{n \ge 1} L_a^2(\mathbb{R}^{3N}) \quad (16)$$

such that $\mathrm{Tr}(\Gamma) = 1$. The set of grand canonical states is written $S_{\mathcal{F}_a}$. We may then easily identify the sets integer-electronic and PPLB states to subspaces of $S_{\mathcal{F}_a}$ and it might feel uneasing that the ensembles $S_N^\xi$ used to define N-centered DFT steer away from this unifying setting. In this section we propose an equivalent construction of N-centered DFT arising from a slightly different point of view in order to establish it as a physical state theory.

We shall keep the N-centered energy unchanged:

$$E_N^\xi(v) = \inf \left\{ (1 - \xi) \mathrm{Tr}(H_v^N \Gamma_N) + \frac{N\xi}{N+1} \mathrm{Tr}(H_v^{N+1} \Gamma_{N+1}) \mid \Gamma_k \in \mathcal{S}_k \right\}, \quad (17)$$

but now choose to write it as the ground-state energy of a modified hamiltonian over $(N + \xi)$-electronic states:

$$E_N^\xi(v) = \inf \{ \mathrm{Tr}(H_v^{Nc}\Gamma) \mid \Gamma \in \mathcal{S}_{N+\xi} \}, \quad (18)$$

where

$$H_v^{Nc} = H_v^N \oplus \frac{N}{N+1} H_v^{N+1}. \quad (19)$$

Of course, the constrained-search method then yields:

$$E_N^\xi(v) = \inf_{\rho \in \mathcal{I}_N} \left\{ F_N^\xi(\rho) + \int v \rho \right\}, \quad (20)$$

where $F_N^\xi$ is the same N-centered universal functional as was defined in the previous paragraph, but whose expression adapted to the setting considered here will be

$$F_N^\xi(\rho) = \inf \{ \mathrm{Tr}(H_0^{Nc}\Gamma) \mid \Gamma \in \mathcal{S}_{N+\xi}, \rho_\Gamma^N = \rho \}, \quad (21)$$
where we defined the N-centered density of a state $\Gamma = (1-w)\Gamma_N \oplus w\Gamma_{N+1}$ by:

$$ \rho_{\Gamma}^{Nc} = (1 - \xi)\rho_{\Gamma_N} + \frac{N\xi}{N+1}\rho_{\Gamma_{N+1}}. \qquad (22) $$

Before letting our focus go back to the original formulation of N-centered DFT in the rest of this paper, we may note here that by extending the previous hamiltonian to the whole Fock space:

$$ H_v^{Nc} = \bigoplus_{k \ge 1} \frac{N}{k} H_v^k, \qquad (23) $$

we could define a "grand canonical N-centered DFT", whose relation with right N-centered DFT would be the same as the one between grand-canonical DFT and PPLB DFT.

## 1.4 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 regular DFT context, their utility resides elsewhere: we shall see in this Section that they enable us to justify (convex) duality relations linking $F_N^\xi$ and $E_N^\xi$, which will then make it possible to set N-centered eDFT 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 [7]). 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 10** Consider a Banach space $X$, a convex functional $f$ on $X$ and a concave functional $g$ on $X^*$. Letting $\langle \cdot, \cdot \rangle$ be the duality pairing between $X$ and $X^*$ we define the following objects.

The Legendre transform $f^*$, respectively double Legendre transform $f^{**}$, of $f$ is the functional on $X^*$, respectively $X$, defined by

$$ f^*(v) = \sup \{ \langle v, \rho \rangle - f(\rho) | \rho \in X \} $$

$$ f^{**}(\rho) \sup \{ \langle v, \rho \rangle - f^{*}(v) | v \in X^{*} \}. $$

The concave conjugate $f^\wedge$ of $f$ is the functional on $X^*$ defined by:

$$ f^\wedge(v) = \inf \{ f(\rho) + \langle v, \rho \rangle | \rho \in X \}. $$

The convex conjugate $g^\vee$ of $g$ is the functional on $X$ defined by:

$$ g^\vee(\rho) = \sup \{ g(v) - \langle v, \rho \rangle | v \in X^* \}. $$
**Lemma 11** Legendre transforms and concave/convex conjugates are linked by the following relations:

$$f^{\wedge}(v) = -f^{*}(-v) \\ (f^{\wedge})^{\vee} = f^{**}$$

The following result clarifies the usefulness of these definitions.

**Proposition 12** Energy and universal functionals are set in duality via concave and convex conjugates:

$$E_N^{\xi} = (F_N^{\xi})^{\wedge} \qquad (24)$$

$$F_N^{\xi} = (E_N^{\xi})^{\vee}. \qquad (25)$$

**Proof.** The first relation (24) is simply the one obtained via the constrained-search method (7). According to Lemma (11), it follows that: $(E_N^{\xi})^{\vee} = (F_N^{\xi})^{**}$. Since $F_N^{\xi}$ is convex, proper and lower semicontinuous, it is equal to its double Legendre transform and the second relation (25) is obtained. □ With

this result at hand, a strong analogy between $N$-electronic and $N$-centered eDFT 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 $\rho$ such that

$$E_N^{\xi}(v) = F_N^{\xi}(\rho) + \int v\rho.$$

Let us write $GSD^{\xi}(v)$ the set of such densities, and $GS^{\xi}(v)$ the set of $(N$-centered) ground states for $v$, i.e. states $\Gamma \in S_N^{\xi}$ such that

$$E_N^{\xi}(v) = \mathrm{Tr}(H_0\Gamma) + \int \mathcal{D}_N^{\xi}(\Gamma)v.$$

**Proposition 13** Consider some potential $v$. There exists a ground state density $\rho$ for this potential if and only if there exists an $N$-centered ground state $\Gamma$ for $v$. Furthermore, the following set equality holds true:

$$GSD^{\xi}(v) = \{\rho \mid \exists \Gamma \in GS^{\xi}(v), \mathcal{D}_N^{\xi}(\Gamma) = \rho\}.$$

**Proof.** Assume that $\rho \in GSD^{\xi}(v)$. From the first item of Theorem (9) it follows that there exists a state $\Gamma \in S_N^{\xi}$ such that $\mathcal{D}_N^{\xi}(\Gamma) = \rho$ and $F_N^{\xi}(\rho) = \mathrm{Tr}(H_0\Gamma)$, which yields:

$$E_N^{\xi}(v) = F_N^{\xi}(\rho) + \int v\rho = \mathrm{Tr}(H_0\Gamma) + \int \mathcal{D}_N^{\xi}(\Gamma)v,$$
so that $\Gamma \in GS^\xi(v)$.

Conversely, consider $\rho = D_N^\xi(\Gamma_0)$ for some $\Gamma_0 \in GS^\xi(v)$. Then for all $\Gamma \in S_N^\xi$ such that $D_N^\xi(\Gamma) = \rho$:

$$ \mathrm{Tr}(H_0\Gamma_0) = E_N^\xi(v) - \int v\rho \le \mathrm{Tr}(H_0\Gamma) + \int D_N^\xi(\Gamma)v - \int v\rho = \mathrm{Tr}(H_0\Gamma), $$

so that $\mathrm{Tr}(H_0\Gamma_0) = F_N^\xi(\rho)$ and

$$ E_N^\xi(v) = F_N^\xi(\rho) + \int v\rho. $$

□

**Remark 14** Letting $\partial^\wedge E_N^\xi(v)$ denote the superdifferential of $E_N^\xi$ at $v$ we can see that

$$ GSD^\xi(v) = \partial^\wedge E_N^\xi(v) \cap (L^1 \cap L^3). $$

Furthermore, for all $\rho \in L^1 \cap L^3$, the following equivalence holds true:

$$ -v \in \partial^\vee F_N^\xi(\rho) \Leftrightarrow \rho \in \partial^\wedge E_N^\xi(v). \qquad (26) $$

## 1.5 Investigating the weight-dependence of the universal functionals

We have established that for a fixed parameter $\xi$, $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-dependence for a given fixed density $\rho$. To do so let us define for all $\xi \in [0, N/(N+1)]$:

$$ W_N^\rho(\xi) := F_N^\xi(\rho), $$

and set $W_N^\rho(\xi) = +\infty$ elsewhere.

The importance of this parameter-dependent functional is expanded upon in [20, 11], 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^\rho$. 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^\rho$ is convex, since by definition it is the point-wise supremum of the family of convex functions

$$ \xi \in [0, N/(N+1)] \mapsto E_N^\xi(v) - \int v\rho, \quad v \in L^{3/2} + L^\infty. $$

Another advantage of this convexity property is that we may study the differentiability of $W_N^\rho$ through its subdifferential. This leads us to the following proposition.
**Proposition 15** Consider $\rho \in L^1 \cap L^3$. Assume that $\rho \in \text{Dom}(\partial F_N^\xi)$ for all $\xi \in [0, N/(N+1)]$, choose $-v_\xi \in \partial F_N^\xi(\rho)$ and define

$$ \Omega_N(v) = E^{N+1}(v) - \frac{N+1}{N} E^N(v). $$

Further assume that the Hohenberg-Kohn property is satisfied: each $v_\xi$ is unique up to an additive constant. Then for all $\xi \in ]0, N/(N+1)]$, the two scalars

$$ m_N^\rho(\xi) = \sup_{\beta < \xi} \Omega_N(v_\beta) $$

$$ M_N^\rho(\xi) = \inf_{\beta > \xi} \Omega_N(v_\beta) $$

are well defined and

$$ \partial W_N^\rho(\xi) \subset [m_N^\rho(\xi), M_N^\rho(\xi)]. $$

**Proof.** First off, let us verify that the quantities $m_N^\rho(\xi)$ and $M_N^\rho(\xi)$ are well defined, i.e that they are independent of the choice of the family $v_\xi$. The proof is straightforward: we show that for any other family $-v'_\xi \in \partial F_N^\xi(\rho)$ we have $\Omega_N(v'_\xi) = \Omega_N(v_\xi)$. By the Hohenberg-Kohn assumption, for all $\xi$ there exists a constant $C_\xi$ such that $v'_\xi = v_\xi + C_\xi$, which in turn yields:

$$ \Omega_N(v'_\xi) = [E^{N+1}(v_\xi) + (N+1)C_\xi] - \frac{N+1}{N}[E^N(v_\xi) + NC_\xi] = \Omega_N(v_\xi). $$

Now consider some $\mu \in \partial W_N^\rho(\xi)$:

$$ \forall \beta \in [0, N/(N+1)], W_N^\rho(\beta) - W_N^\rho(\xi) \geq \mu(\beta - \xi). $$

Then for all $\beta \in [0, N/(N+1)]$ and $v \in L^{3/2} + L^\infty$:

$$ E_N^\xi(v) - \int v\rho \leq W_N^\rho(\beta) - \mu(\beta - \xi), $$

and since

$$ E_N^\xi(v) = E_N^\beta(v) + (\xi - \beta)\Omega_N(v), $$

this can be rewritten under the form:

$$ E_N^\beta(v) - \int v\rho \leq F_N^\beta(\rho) + (\beta - \xi)[\Omega_N(v) - \mu], \forall \beta \in [0, N/(N+1)], \forall v. $$

Since $\rho \in \text{Dom}(\partial F_N^\beta)$ for all $\beta \in [0, N/(N+1)]$, by setting $-v = -v_\beta \in \partial F_N^\beta(\rho)$ in the previous equality we get:

$$ 0 \leq (\beta - \xi)[\Omega_N(v_\beta) - \mu], \forall \beta \in [0, 1], $$

which yields:

$$ \sup_{\beta < \xi} \Omega_N(v_\beta) \leq \mu \leq \inf_{\beta > \xi} \Omega_N(v_\beta). $$
□

The mindset that led us to the previous result was to obtain a characterization of $\partial W_N^\rho$. 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:

**Proposition 16** Consider $-v_\xi \in \partial F_N^\xi(\rho)$. For all $\beta > \xi$ and $\mu < \xi$ such that $-v_\beta \in \partial F_N^\beta(\rho)$ and $-v_\mu \in \partial F_N^\mu(\rho)$ the following holds:

$$ \Omega_N(v_\xi) \in \partial W_N^\rho(\xi) \subset [\Omega_N(v_\mu), \Omega_N(v_\beta)]. $$

**Proof.** The same proof as that of the previous result yields that $\partial W_N^\rho(\xi)$ is a subset of $[\Omega_N(v_\mu), \Omega_N(v_\beta)]$. We are then left to note that for all $\beta \in \mathbb{R}$:

$$ \begin{aligned} W_N^\rho(\beta) - W_N^\rho(\xi) &= W_N^\rho(\beta) - E_N^\xi(v_\xi) + \int v_\xi \rho \\ &\ge E_N^\beta(v_\xi) - E_N^\xi(v_\xi) \\ &= (\beta - \xi)\Omega_N(v_\xi), \end{aligned} $$

which precisely means that $\Omega_N(v_\xi)$ is in $\partial W_N^\rho(\xi)$. □

Although not completely explicit, the bounds of the subdifferential of $W_N^\rho$ provided in the previous proposition do contribute to justify that weight variations of the N-centered functional are closely linked to energy gaps.

## 2 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 [20] within the mathematical formalism presented previously, thus identifying that this approach leads to similar mathematical problems as the Kohn-Sham method for N-electronic DFT. We will then be taking a quick look at a method proposed in [13] to solve some of these problems, which will further establish the analogy between N-centered and N-electronic DFT.

### 2.1 General idea

The ultimate goal of the Kohn-Sham method is to find an element of $GSD^\xi(v)$ for a given $v \in L^\infty + L^{3/2}$. As seen previously, this amounts to finding a density $\rho$ such that

$$ -v \in \partial^\vee F_N^\xi(\rho). \qquad (27) $$

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 $E_{N,\emptyset}^{\xi}$ and $F_{N,\emptyset}^{\xi}$:

$$E_{N,\emptyset}^{\xi}(v) = \inf \left\{ \mathrm{Tr}(K\Gamma) + \int \mathcal{D}_{N}^{\xi}(\Gamma)v|\Gamma \in \mathcal{S}_{N}^{\xi} \right\},$$

where

$$K = \bigoplus_{k=N-1}^{N+1} K^k,$$

$$K^k = -\sum_{i=1}^{k} \Delta_i,$$

and

$$F_{N,\emptyset}^{\xi}(\rho) = \inf \left\{ \mathrm{Tr}(K\Gamma)|\Gamma \in \mathcal{S}_{N}^{\xi}, \mathcal{D}_{N}^{\xi}(\Gamma) = \rho \right\}.$$

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 (27) into a noninteracting ground-state density search problem, ideally by finding a functional $v_s: L^1 \cap L^3 \to L^\infty + L^{3/2}$ such that:

$$-v \in \partial^{\vee} F_{N}^{\xi}(\rho) \Leftrightarrow -v_s(\rho) \in \partial F_{N,\emptyset}^{\xi}(\rho). \qquad (28)$$

Most of the literature on this subject (in the $N$-electronic case) assumes differentiability of the exchange-correlation energy $E_{xc,N}^{\xi}$:

$$E_{xc,N}^{\xi} = F_{N}^{\xi} - F_{N,\emptyset}^{\xi}$$

in order to yield the following Kohn-Sham equation:

$$-v \in \partial^{\vee} F_{N}^{\xi}(\rho) \Leftrightarrow -(v + v_{KS,N}^{\xi}(\rho)) \in \partial F_{N,\emptyset}^{\xi}(\rho),$$

where

$$v_{KS,N}^{\xi} = \frac{\delta E_{xc,N}^{\xi}}{\delta \rho}. \qquad (29)$$

Unfortunately, as was noted in [13], the set $\mathcal{I}_N$ on which the exchange-correlation energy is defined has an empty interior in $L^1 \cap L^3$, so that the meaning of the differential notation used in (29) is rather opaque: the exchange-correlation energy is obviously not differentiable in the usual sense.

## 2.2 Generalized Kohn-Sham equations

In order to overcome the problem underlined in the previous subsection, a regularization method has been proposed in [13], which leads to a less restrictive Kohn-Sham approach than the one embodied by the equivalence (28): rather than trying to find a Kohn-Sham potential for a fixed density $\rho$, the density
itself is allowed to be modified. Formally, the goal is now to find two functionals $\rho_s: L^1 \cap L^3 \to L^\infty + L^{3/2}$ and $v_s: L^1 \cap L^3 \to L^1 \cap L^3$ such that:

$$ -v \in \partial^\vee F_N^\xi(\rho) \Leftrightarrow -v_s(\rho) \in \partial F_{N, \emptyset}^\xi(\rho_s(\rho)). $$

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 eDFT. 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 [4] for details.

**Definition 17** Let $(\mathcal{H}, \langle \cdot, \cdot \rangle)$ be a Hilbert space, $f$ a functional on $\mathcal{H}$ and $\varepsilon > 0$. The $\varepsilon$-Moreau-Yosida regularization $MY_\varepsilon(f)$ of $f$ is the functional on $\mathcal{H}$ defined by:

$$ MY_\varepsilon(f)(\rho) = \left( f \Box \frac{1}{2\varepsilon} \| \cdot \|^2 \right) (\rho) = \inf_{n \in \mathcal{H}} \left\{ f(n) + \frac{1}{2\varepsilon} \|\rho - n\|^2 \right\}. \quad (30) $$

**Proposition 18** Let $f \in \Gamma_0(\mathcal{H})$, $\varepsilon > 0$. Then:

1. $MY_\varepsilon(f)$ is convex, real valued and exact.

2. The infimum in (29) is uniquely attained at some point $\text{Prox}_{\varepsilon f}(\rho)$.

3. $MY_\varepsilon(f)$ is Fréchet-differentiable on $\mathcal{H}$ and its gradient

$$ \nabla MY_\varepsilon(f) = \frac{Id - \text{Prox}_{\varepsilon f}}{\varepsilon} $$

is $\varepsilon^{-1}$-Lipschitz continuous.

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 \cap 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 [2].

A first step towards working in a reflexive context is made by recalling the inclusion $L^1 \cap L^3 \subset L^2$, which combined with (7) yields:

$$ E_N^\xi(v) = \inf \left\{ F_N^\xi(\rho) + \int v \rho |\rho| \in L^2(\mathbb{R}^3) \right\}. \quad (31) $$

Unfortunately, this inclusion does not imply that $\Gamma_0(L^1 \cap L^3) \subset \Gamma_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

$$ F_N^\xi(\rho) = \sup \left\{ E_N^\xi(v) - \int v \rho |v| \in L^2(\mathbb{R}^3) \right\}, \forall \rho \in L^2(\mathbb{R}^3), \quad (32) $$
so that it becomes necessary to define a new “universal functional” in duality with the energy on $L^2(\mathbb{R}^3)$. This is made possible by the continuity of $E_N^\xi$ on $L^2(\mathbb{R}^3)$, leading to:

$$E_N^\xi(v) = \inf \left\{ \overline{F}_N^\xi(\rho) + \int v\rho |\rho| \in L^2(\mathbb{R}) \right\}, \quad (33)$$

where

$$\overline{F}_N^\xi(\rho) = \sup \left\{ E_N^\xi(v) - \int v\rho |v| \in L^2(\mathbb{R}^3) \right\}.$$

**Lemma 19** For all $\rho \in L^2(\mathbb{R}^3)$, $\overline{F}_N^\xi(\rho) \le F_N^\xi(\rho)$. Furthermore, for all $\rho \in \mathcal{I}_N$ the following equivalence holds:

$$\overline{F}_N^\xi(\rho) = F_N^\xi(\rho) \text{ and } -v \in \partial \overline{F}_N^\xi(\rho) \Leftrightarrow -v \in \partial F_N^\xi(\rho) \cap L^2(\mathbb{R}^3)$$

**Proof.** The first statement is a direct consequence of the inclusion $L^2 \subset L^\infty + L^{3/2}$.

Next consider $-v \in \partial F_N^\xi(\rho) \cap L^2(\mathbb{R}^3)$. We may then write:

$$F_N^\xi(\rho) = E_N^\xi(v) - \int v\rho \le \overline{F}_N^\xi(\rho) \le F_N^\xi(\rho) \Rightarrow -v \in \partial \overline{F}_N^\xi(\rho).$$

Conversely, assume that $-v \in \partial \overline{F}_N^\xi(\rho)$ and $\overline{F}_N^\xi(\rho) = F_N^\xi(\rho)$. We immediately see that $v \in L^2(\mathbb{R}^3)$. It then follows that for all $n \in L^1 \cap L^3$:

$$F_N^\xi(n) \ge \overline{F}_N^\xi(n) \ge \overline{F}_N^\xi(\rho) + \int v(n-\rho) \\ = F_N^\xi(\rho) + \int v(n-\rho),$$

so that $-v \in \partial F_N^\xi(\rho)$. $\square$

One of the motivations behind $N$-centered density functional theory was to construct universal functionals with the same domain for any choice of weights, which was not the case for preexisting ensemble approaches such as PPLB DFT. The following result shows that this still holds true for the functionals $\overline{F}_N^\xi$.

**Proposition 20** For any admissible choice of parameter $\xi$, the domain of $\overline{F}_N^\xi$ is $\mathcal{I}_N$.

**Proof.** From (33) and lemma 19 we get:

$$E_N^\xi(v) \le \inf \left\{ \overline{F}_N^\xi(\rho) + \int v\rho |\rho| \in \mathcal{I}_N \right\} \le \inf \left\{ F_N^\xi(\rho) + \int v\rho |\rho| \in \mathcal{I}_N \right\} = E_N^\xi(v).$$
Thus, letting $\overline{G}_N^\xi$ be the functional equal to $\overline{F}_N^\xi$ on $\mathcal{I}_N$ and $+\infty$ elsewhere, the functionals $\overline{G}_N^\xi$ and $\overline{F}_N^\xi$ are both in $\Gamma_0(L^2(\mathbb{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 $\overline{F}_N^\xi$ to the noninteracting case, yielding a noninteracting functional $\overline{F}_{N,\emptyset}^\xi$ in duality with the noninteracting $N$-centered energy in $L^2$. Let us now show that the method proposed in [13] in order to link the subdifferentials of $\overline{F}_N^\xi$ and $\overline{F}_{N,\emptyset}^\xi$ can be carried out in the $N$-centered context with no significant modifications.

**Proposition 21** (Generalized Kohn-Sham equations for $N$-centered eDFT)
The following equivalence holds for all $\varepsilon > 0$:

$$ -v \in \partial \overline{F}_N^\xi(\rho) \Leftrightarrow -(v + v_{KS,N}^\xi(\rho - \varepsilon v)) \in \partial \overline{F}_{N,\emptyset}^\xi(\rho + \varepsilon v_{KS,N}^\xi(\rho - \varepsilon v)), $$

where

$$ v_{KS,N}^{\xi}(\rho) = \nabla E_{xc,N}^{\xi}(\rho) \\ E_{xc,N}^{\xi} = M Y_{\varepsilon}(\overline{F}_{N}^{\xi}) - M Y_{\varepsilon}(\overline{F}_{N,\emptyset}^{\xi}) $$

**Proof.** 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:

$$ \begin{align*} -v \in \partial \overline{F}_N^\xi(\rho) &\Leftrightarrow \rho \in \partial E_N^\xi(v) \\ &\Leftrightarrow \rho - \varepsilon v \in \partial (M Y_\varepsilon (\overline{F}_N^\xi))^\wedge (v) \\ &\Leftrightarrow -v = \nabla M Y_\varepsilon (\overline{F}_N^\xi) (\rho - \varepsilon v) \\ &\Leftrightarrow -(v + \nabla E_{xc,N}^\xi(\rho - \varepsilon v)) = \nabla M Y_\varepsilon (\overline{F}_{N,\emptyset}^\xi) (\rho - \varepsilon v) \\ &\Leftrightarrow \rho - \varepsilon v \in \partial (M Y_\varepsilon (\overline{F}_N^\xi))^\wedge (v + \nabla E_{xc,N}^\xi(\rho - \varepsilon v)) \\ &\Leftrightarrow \rho + \varepsilon \nabla E_{xc,N}^\xi(\rho - \varepsilon v) \in \partial E_N^\xi(v + \nabla E_{xc,N}^\xi(\rho - \varepsilon v)) \\ &\Leftrightarrow -(v + \nabla E_{xc,N}^\xi(\rho - \varepsilon v)) \in \partial \overline{F}_{N,\emptyset}^\xi(\rho + \varepsilon \nabla E_{xc,N}^\xi(\rho - \varepsilon v)). \end{align*} $$

□

**Remark 22** 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:
$$
\begin{align*}
\rho_0 &\in L^2 \\
v_i &= v_{KS,N}^{\xi}(\rho_i - \epsilon v) \\
\text{solve } &-(v + v_i) \in \partial \bar{F}_{N,\emptyset}^{\xi}(\rho) \\
\text{set } \rho_{i+1} &= \rho - \epsilon v_i.
\end{align*}
$$

If the convergence of this scheme for the usual Kohn-Sham method was already a problem to be considered, having set up this generalized scheme does not simplify matters. Indeed, it is now necessary to evaluate a regularized exchange-correlation potential, with no clear idea of the physical meaning of this object, which itself is defined through a generalized universal functional itself carrying a less clear physical sense than the true functional $F_N^\xi$. However carrying on the convergence analysis with this (relaxed) functional might be more appropriate at least as it would be set in a Hilbertian framework.

3 Semi-classical limits for N − centered eDFT

This far, all properties of conventional DFT we have investigated have been shown to hold in the $N$-centered eDFT 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 eDFT context.

3.1 The $N$-centered semiclassical limit $\mathcal{E}_{0T}^{N,\xi}$

Let us begin by recalling some results established in [15] regarding this semi-
classical limit.

**Definition 23** For all $N$ and $\rho \in \mathcal{I}_N$ we can define the Coulomb multi-marginal problem:

$$
\mathcal{E}_{OT}^{N}(\rho) = \min_{\mathbb{P} \in \Pi_N(\rho/N)} \int_{\mathbb{R}^{3N}} \sum_{1 \le i < j \le N} \frac{1}{|x_i - x_j|} d\mathbb{P},
$$

where $\prod_N(\rho / \int \rho)$ is the set of symmetric probabilities on $\mathbb{R}^{3N}$ with specified marginal:

$$
\frac{\rho(x)}{N} = \int_{\mathbb{R}^{3(N-1)}} d\mathbb{P}(x, x_2, \dots, x_N)
$$

**Theorem 24** For all $\rho \in \mathcal{I}_N$:

$$
\lim_{\eta \to 0} \min \left\{ \mathrm{Tr}(H_0^{N,\eta} \Gamma) \mid \Gamma \in \mathcal{S}_N, \mathcal{D}_N(\Gamma) = \rho \right\} = \mathcal{E}_{\mathrm{OT}}(\rho),
$$

where

$$
H_0^{N,\eta} = \eta K^N + V_{ee}^N.
$$
Thereafter, we write:

$$ \min \left\{ \mathrm{Tr}(H_0^{N,\eta}\Gamma) |\Gamma \in \mathcal{S}_N, \mathcal{D}_N(\Gamma) = \rho \right\} = F_{\eta}^N(\rho). $$

**Theorem 25** Consider $\chi \in C_c^\infty(\mathbb{R}^3, \mathbb{R})$ a radial function with support in the unit ball such that $\int \chi^2 = 1$ and $\mathbb{P}$ a symmetric probability measure on $\mathbb{R}^{3N}$ with support in

$$ D_{N,\alpha}^c := \{(x_1, \dots, x_N) | x_i - x_j| \ge \alpha, \forall i \ne j\}, $$

for some $\alpha > 0$ and such that

$$ \rho_{\mathbb{P}}(x) = N \int d\mathbb{P}(x, x_2, \dots, x_N) \in \mathcal{I}_N. $$

For all $\alpha/4 > \varepsilon > 0$ there exists $\Gamma_\varepsilon \in \mathcal{S}_N$ such that:

$$ \mathcal{D}_N(\Gamma) = \rho_{\mathbb{P}} $$

$$ \mathrm{Tr}(K^N \Gamma_\varepsilon) = \int |\nabla \sqrt{\rho_{\mathbb{P}}}|^2 + \frac{N}{\varepsilon^2} \int |\nabla \chi|^2 $$

and for all symmetric $\Phi \in C^2(D_{\alpha-4\varepsilon}^c)$:

$$ \begin{aligned} |\mathrm{Tr}(\Phi\Gamma_\varepsilon) - \int \Phi d\mathbb{P}| \le \varepsilon^2 \frac{1}{N} \sum_{j=1}^N \|\nabla_j \Phi\|_{L^\infty(D_{\alpha-4\varepsilon}^c)} \left( \int |\nabla \rho_{\mathbb{P}}| \right) \left( \int |u|^2 \chi(u)^2 du \right) \\ + 2 \sum_{j,k=1}^N \|\nabla_j \nabla_k \Phi\|_{L^\infty(D_{\alpha-4\varepsilon}^c)} \end{aligned} $$

Let us now go back to the $N$-centered case. For simplicity we limit ourselves to the left $N$-centered case. Let us define

$$ H_0^\eta = H_0^{\eta,N} \oplus H_0^{\eta,N+1}, $$

$$ F_{N,\eta}^{\xi}(\rho) = \min \left\{ \mathrm{Tr}(H_0^{\eta}\Gamma) |\Gamma \in \mathcal{S}_N^{\xi}, \mathcal{D}_N^{\xi}(\Gamma) = \rho \right\} = \mathrm{Tr}(H_0^{\eta}\Gamma_{N,\eta}^{\xi}(\rho)). $$

The proof of the following result is based on the following statement : for all $\Gamma = w_N^\xi \Gamma^N \oplus \xi \Gamma^{N+1} \in \mathcal{S}_N^\xi$ with $\mathcal{D}_N^\xi(\Gamma) = \rho$ we have:

$$ F_{N,\eta}^{\xi}(\rho) \le \eta \mathrm{Tr}(K\Gamma) + w_N^{\xi} \int V_{ee}^N(x) \Gamma^N(x, x) dx + \xi \int V_{ee}^{N+1}(x) \Gamma^{N+1}(x, x) dx. $$

**Proposition 26** For all $\rho \in \mathcal{I}_N$,

$$ \lim_{\eta \to 0} F_{N,\eta}^{\xi}(\rho) = \mathcal{E}_{OT}^{N,\xi}(\rho) $$

where

$$ \begin{aligned} \mathcal{E}_{OT}^{N,\xi}(\rho) = \inf \left\{ w_N^{\xi} \mathcal{E}_{OT}^N(\rho_N) + \xi \mathcal{E}_{OT}^{N+1}(\rho_{N+1}) | \rho_N \in \mathcal{I}_N, \rho_{N+1} \in \mathcal{I}_{N+1}, \right. \\ \left. \rho = w_N^{\xi} \rho_N + \xi \rho_{N+1} \right\} \end{aligned} $$
**Proof.** For all $\eta$ write:

$$ \Gamma_{N,\eta}^{\xi}(\rho) = w_{N}^{\xi} \Gamma_{N,\eta} \oplus \xi \Gamma_{N+1,\eta}. $$

Since $\mathcal{D}_N^\xi(\Gamma_{N,\eta}^\xi(\rho)) = w_N^\xi \rho_\eta^N + \xi \rho_\eta^{N+1} = \rho$ it follows that

$$ 
\begin{aligned}
F_{N,\eta}^{\xi}(\rho) &= w_{N}^{\xi} \mathrm{Tr}(H_{0}^{N,\eta} \Gamma_{N,\eta}) + \xi \mathrm{Tr}(H_{0}^{N+1,\eta} \Gamma_{N+1,\eta}) \\
&\geq w_{N}^{\xi} F_{\eta}^{N}(\rho_{\eta}^{N}) + \xi F_{\eta}^{N+1}(\rho_{\eta}^{N+1}) \\
&\geq w_{N}^{\xi} \mathcal{E}_{OT}(\rho_{\eta}^{N}) + \xi \mathcal{E}_{OT}(\rho_{\eta}^{N+1}) \\
&\geq \mathcal{E}_{OT}^{N,\xi}(\rho).
\end{aligned}
\quad (34) $$

Conversely, consider $\rho_N \in \mathcal{I}_N$ and $\rho_{N+1} \in \mathcal{I}_{N+1}$ such that $\rho = w_N^\xi \rho_N + \xi \rho_{N+1}$. Let $\mathbb{P}_N$ (respectively $\mathbb{P}_{N+1}$) be a minimizer for $\mathcal{E}_{OT}^N(\rho_N)$ (respectively for $\mathcal{E}_{OT}^{N+1}(\rho_{N+1})$) and $\Gamma_\epsilon^N$ (respectively $\Gamma_\epsilon^{N+1}$) the state yielded by theorem 25 for $\mathbb{P}_N$ (respectively $\mathbb{P}_{N+1}$). There exists $\beta, \delta$ such that $\mathbb{P}_N$ has its support in $D_{N,\beta}^c, \mathbb{P}_{N+1}$ in $D_{N+1,\delta}^c$. Defining $\Gamma_\epsilon = (1-\alpha)\Gamma_\epsilon^N + \alpha_+ \Gamma_\epsilon^{N+1}$ we may then write:

$$ 
\begin{aligned}
F_{N,\eta}^{\xi}(\rho) &\leq \eta \mathrm{Tr}(K\Gamma_{\epsilon}) + \mathrm{Tr}(V_{ee}\Gamma_{\epsilon}) \\
&= \eta \left[ w_{N}^{\xi} \mathrm{Tr}(K^{N}\Gamma_{\epsilon}^{N}) + \xi \mathrm{Tr}(K^{N+1}\Gamma_{\epsilon}^{N+1}) \right] + w_{N}^{\xi} \mathrm{Tr}(V_{ee}^{N}\Gamma_{\epsilon}^{N}) + \xi \mathrm{Tr}(V_{ee}^{N+1}\Gamma_{\epsilon}^{N+1}).
\end{aligned}
$$

The properties of $\Gamma_\epsilon^N$ and $\Gamma_\epsilon^{N+1}$ lead to:

$$ 
\begin{aligned}
\eta \mathrm{Tr}(K^N \Gamma_\epsilon) + \mathrm{Tr}(V_{ee}^N \Gamma^N) &\leq \eta \int |\nabla \sqrt{\rho_N}|^2 + \frac{N}{\epsilon^2} \int |\nabla \chi|^2 + \mathcal{E}_{OT}^N(\rho_N) \\
&\quad + C_N \epsilon^2 \left( \frac{1}{(\beta - 4\epsilon)^2} \int |\nabla \rho_N| \int |u|^2 \chi(u)^2 + \frac{1}{(\beta - 4\epsilon)^3} \right)
\end{aligned}
$$

and

$$ 
\begin{aligned}
\eta \mathrm{Tr}(K^{N+1} \Gamma_\epsilon) + \mathrm{Tr}(V_{ee}^N \Gamma^N) &\leq \eta \int |\nabla \sqrt{\rho_{N+1}}|^2 + \frac{N+1}{\epsilon^2} \int |\nabla \chi|^2 + \mathcal{E}_{OT}^{N+1}(\rho_{N+1}) \\
&\quad + C_{N+1} \epsilon^2 \left( \frac{1}{(\delta - 4\epsilon)^2} \int |\nabla \rho_{N+1}| \int |u|^2 \chi(u)^2 + \frac{1}{(\delta - 4\epsilon)^3} \right).
\end{aligned}
$$

Thus, taking $\eta$ small enough and $\epsilon = \eta^{1/4}$, the previous inequalities yield:

$$ F_{N,\eta}^{\xi}(\rho) \leq C(\rho_N, \rho_{N+1})(\eta + \sqrt{\eta}) + w_N^{\xi} \mathcal{E}_{OT}(\rho_N) + \xi \mathcal{E}_{OT}(\rho_{N+1}), \quad (35) $$

Consider any sequence $(\eta_n)$ that converges to 0. It follows from inequalities (34) and (35) that

$$ \mathcal{E}_{OT}^{N,\xi}(\rho) \leq \liminf_{n \to \infty} F_{N,\eta_n}^{\xi}(\rho) \leq \limsup_{n \to \infty} F_{N,\eta_n}^{\xi}(\rho) \leq w_N^{\xi} \mathcal{E}_{OT}(\rho_N) + \xi \mathcal{E}_{OT}(\rho_{N+1}). $$

This being true for all $\rho_N, \rho_{N+1}$ such that $\rho = w_N^\xi \rho_N + \xi \rho_{N+1}$, we infer that

$$ \mathcal{E}_{OT}^{N,\xi}(\rho) \leq \liminf_{n \to \infty} F_{N,\eta_n}^{\xi}(\rho) \leq \limsup_{n \to \infty} F_{N,\eta_n}^{\xi}(\rho) \leq \mathcal{E}_{OT}^{N,\xi}(\rho), $$
which then yields the result $\lim_{\eta \to 0} F_{N,\eta}^{\xi}(\rho) = \mathcal{E}_{OT}^{N,\xi}(\rho)$. $\square$

**Remark 27** The result proven there-before is weaker than the one shown by Lewin, which stated the existence of a constant $C(\rho)$ such that

$$F_{\eta}^{N}(\rho) \leq C(\rho)(\eta + \sqrt{\eta}) + \mathcal{E}_{OT}(\rho).$$

If such a result holds in the $N$-centered context a further analysis of $C(\rho_N, \rho_{N+1})$ shall be necessary to yield it, in particular it would suffice to find $\tilde{C}$ such that

$$C(\rho_N, \rho_{N+1}) \leq \tilde{C}(\rho), \forall \rho = (1 - \alpha)\rho_N + \alpha_+\rho_{N+1}.$$

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 classical optimal transport problem.

## 3.2 An upper bound of $\mathcal{E}_{OT}^{N+1}$ which is an $N$-electronic optimal transport problem

In this subsection we present the construction of an upper bound for $\mathcal{E}_{OT}^{N,\xi}$, 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 $\mathcal{E}_{OT}^N$ and $\mathcal{E}_{OT}^{N+1}$ more clear and uncover an alternative expression of $\mathcal{E}_{OT}^{N,\xi}$. The initial thought was to construct a bijection between the probability sets $\prod_N(\rho/N)$ and $\prod_{N+1}(\rho/N)$ for a given density $\rho \in \mathcal{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 $\mathcal{E}_{OT}^N$ and $\mathcal{E}_{OT}^{N+1}$.

Let us begin by investigating the construction of an injective mapping, while trying to guarantee “as much surjectivity as possible”. Defining the probability measure $\mu_{\rho}(A) = \int_{A} \rho(x)dx$ for borelian sets $A \in \mathcal{B}(\mathbb{R}^3)$, an obvious way to map a probability $P$ in $\prod_N(\rho/N)$ to a probability measure on $\mathbb{R}^{3(N+1)}$ whose marginals are all equal to $\rho/N$ is to define

$$\lambda(P) = P \otimes \frac{\mu_{\rho}}{N}. \qquad (36)$$

Of course this is far from satisfying since $\lambda(P)$ does not belong to $\prod_{N+1}(\rho/N)$ in general: the only case where $\lambda(P)$ would be symmetric is if $P$ were given by $P = \bigotimes_{i=1}^{N} \frac{\mu_{\rho}}{N}$. This can be remedied by defining the map $\Lambda$ obtained by mixing probability measures similar to $\lambda(P)$:

$$\Lambda(P) = \frac{1}{N+1} \sum_{i=1}^{N+1} \frac{\mu_{\rho}(x_i)}{N} \otimes P(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_{N+1}). \in \Pi_{N+1}(\rho/N).$$
**Proposition 28** The map $\Lambda : \Pi_N(\rho/N) \to \Pi_{N+1}(\rho/N)$ is injective.

**Proof.** Let us consider $\mathbb{P}_1, \mathbb{P}_2 \in \Pi_N(\rho/N)$ such that $\Lambda(\mathbb{P}_1) = \Lambda(\mathbb{P}_2)$. For integers $1 \le k \le N$ we propose a recursive proof of the statement:

$$ (H_k) : \text{"}\forall A_1, \dots, A_k \in \mathcal{B}(\mathbb{R}^{3N}), \mathbb{P}_1\left(\bigotimes_{i=1}^{k} A_i \bigotimes_{i=k+1}^{N} \mathbb{R}^3\right) = \mathbb{P}_2\left(\bigotimes_{i=1}^{k} A_i \bigotimes_{i=k+1}^{N} \mathbb{R}^3\right). \text{"} $$

Since $\mathbb{P}_1$ and $\mathbb{P}_2$ are both in $\Pi_N(\rho/N)$, $(H_1)$ is true (they have same marginals). Let us assume that $(H_k)$ holds true for some $k < N$ and write:

$$ \begin{aligned} \Lambda(\mathbb{P}_1) \left( \bigotimes_{i=1}^{k+1} A_i \bigotimes_{i=k+2}^{N+1} \mathbb{R}^3 \right) &= \frac{1}{N+1} \left( \sum_{i=1}^{k+1} \frac{\mu_\rho(A_i)}{N} \mathbb{P}_1 \left( \bigotimes_{j \neq i}^{k+1} A_j \bigotimes \mathbb{R}^3 \right) \right. \\ &\qquad \left. + \sum_{i=k+2}^{N+1} \frac{\mu_\rho(\mathbb{R})}{N} \mathbb{P}_1 \left( \bigotimes_{j=1}^{k+1} A_j \bigotimes \mathbb{R}^3 \right) \right) \\ &= \frac{1}{N+1} \left( \sum_{i=1}^{k+1} \frac{\mu_\rho(A_i)}{N} \mathbb{P}_2 \left( \bigotimes_{j \neq i}^{k+1} A_j \bigotimes \mathbb{R}^3 \right) \right. \\ &\qquad \left. + (N-k) \mathbb{P}_1 \left( \bigotimes_{j=1}^{k+1} A_j \bigotimes \mathbb{R}^3 \right) \right). \end{aligned} $$

Doing the same for $\mathbb{P}_2$ it follows that:

$$ \begin{aligned} 0 &= \Lambda(\mathbb{P}_1) \left( \bigotimes_{i=1}^{k+1} A_i \bigotimes_{i=k+2}^{N+1} \mathbb{R}^3 \right) - \Lambda(\mathbb{P}_2) \left( \bigotimes_{i=1}^{k+1} A_i \bigotimes_{i=k+2}^{N+1} \mathbb{R}^3 \right) \\ &= \frac{N-k}{N+1} \left( \mathbb{P}_1(\bigotimes_{j=1}^{k+1} A_j \bigotimes \mathbb{R}^3) - \mathbb{P}_2(\bigotimes_{j=1}^{k+1} A_j \bigotimes \mathbb{R}^3) \right), \end{aligned} $$

which yields $(H_{k+1})$. $\square$

<figure><img src="image_3.png" /></figure>

This result leads to an upper bound for $\mathcal{E}_{\text{OT}}^{N+1}$ taking the form of an $N$-dimensional optimal transport problem, namely:

$$ \mathcal{E}_{\text{OT}}^{N+1} \left( \frac{N+1}{N} \rho \right) \le \inf_{\mathbb{P} \in \Pi_N(\rho/N)} \int_{\mathbb{R}^{3N+1}} \sum_{1 \le i < j \le N+1} \frac{1}{|x_i - x_j|} d\Lambda(\mathbb{P}) $$

which the following technical lemma will help make more explicit.

**Lemma 29** For all symmetric measurable function $f$,

$$ \int_{\mathbb{R}^{3(N+1)}} f(x_1, \dots, x_{N+1}) d\Lambda(\mathbb{P})(x_1, \dots, x_{N+1}) \\ = \int_{\mathbb{R}^{3N}} \left[ \int_{\mathbb{R}^3} f(x_1, \dots, x_{N+1}) \frac{\rho(x_{N+1})}{N} dx_{N+1} \right] d\mathbb{P}(x_1, \dots, x_N) $$
**Proof.** The definition of $\Lambda(\mathbb{P})$ and the symmetry of $f$ and $\mathbb{P}$ yield:

$$ 
\begin{aligned}
& \int_{\mathbb{R}^{3(N+1)}} f(x_1, \dots, x_N) d\Lambda(\mathbb{P})(x_1, \dots, x_{N+1}) \\
&= \frac{1}{N+1} \sum_{i=1}^{N+1} \int_{\mathbb{R}^{3N}} \left[ \int_{\mathbb{R}^3} f(x_1, \dots, x_{N+1}) \frac{\rho(x_i)}{N} dx_i \right] d\mathbb{P}(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_{N+1}) \\
&= \frac{1}{N+1} \sum_{i=1}^{N+1} \int_{\mathbb{R}^{3N}} \left[ \int_{\mathbb{R}^3} f(x_1, \dots, x_{N+1}) \frac{\rho(x_{N+1})}{N} dx_{N+1} \right] d\mathbb{P}(x_1, \dots, x_N) \\
&= \int_{\mathbb{R}^{3N}} \left[ \int_{\mathbb{R}^3} f(x_1, \dots, x_{N+1}) \frac{\rho(x_{N+1})}{N} dx_{N+1} \right] d\mathbb{P}(x_1, \dots, x_N)
\end{aligned}
$$

□

**Corollary 30** For all $\rho \in \mathcal{I}_N$,

$$ \mathcal{E}_{OT}^{N+1}(\frac{N+1}{N}\rho) \leq \inf_{\mathbb{P} \in \Pi_N(\rho/N)} \int_{\mathbb{R}^{3N}} \left[ \int_{\mathbb{R}^3} \frac{\rho(x_{N+1})}{N} \sum_{1 \leq i < j \leq N+1} \frac{1}{|x_i - x_j|} dx_{N+1} \right] d\mathbb{P}. \quad (37) $$

We may further precise the right hand side of (37) by noting that it is equal to the following expression:

$$ \inf_{\mathbb{P} \in \Pi_N(\rho/N)} \left\{ \int_{\mathbb{R}^{3N}} \left[ \int_{\mathbb{R}^3} \frac{\rho(x_{N+1})}{N} \sum_{1 \leq i \leq N} \frac{1}{|x_i - x_{N+1}|} dx_{N+1} \right] d\mathbb{P}(x_1, \dots, x_N) + \int_{\mathbb{R}^{3N}} \sum_{1 \leq i < j \leq N} \frac{1}{|x_i - x_j|} d\mathbb{P}(x_1, \dots, x_N) \right\}. $$

We hope to sharpen this estimation in future work in order to get closer to a “more classical optimal transport formulation” of the semiclassical limit in the $N$-centered setting. In fact, with regard to the previous proposition, such optimal transport formulation could hold in the framework of weak (nonlinear) optimal transport problems.

# 4 Application to the Hubbard Dimer

In this Section we propose a presentation of the Hubbard dimer, a toy model used in order to test and provide proof-of-concept for DFT-based calculations with at its core an extremely simplified Hamiltonian (which still describes electron-electron interaction), and study this model in the $N$-centered setting. The prerequisite concepts of CAR algebra and second quantization used to define this model can be found in [14].
## 4.1 The Hubbard dimer and its functionals

Based on the second quantization of the hamiltonian obtained for a given orthonormal family, 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(\mathbb{R}^{3N})$ we limit ourselves to four sites $\{f_{1\uparrow}, f_{1\downarrow}, f_{2\uparrow}, f_{2\downarrow}\}$ (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:

$$HD(v_1, v_2) = -t \sum_{s=\uparrow,\downarrow} [a^\dagger(f_{1s})a(f_{2s}) + a^\dagger(f_{2s})a(f_{1s})] + \sum_{i=1}^{2} U_i n_{i\uparrow} n_{i\downarrow} + \sum_{i=1}^{2} v_i n_i, \quad (38)$$

where

$$n_{is} = a^{\dagger}(f_{is})a(f_{is}), \quad s \in \{\uparrow, \downarrow\}$$

$$n_i = n_{i\uparrow} + n_{i\downarrow}.$$

Letting $\mathcal{H} = \text{span}(f_{1\uparrow}, f_{1\downarrow}, f_{2\uparrow}, f_{2\downarrow})$ and $\mathcal{H}_a^N = \text{span}(f_{i_1} \wedge \cdots \wedge f_{i_N})_{i_1 < \cdots < i_N}$, where we have used the ordering convention: $1 \uparrow < 1 \downarrow < 2 \uparrow < 2 \downarrow$. It follows that what we are really interested in here is the action of the second quantized hamiltonian on the (16-dimensional) Fock space $\mathcal{F}_a^4$ (since $\mathcal{H}_a^N = \{0\}$ for all $N \ge 5$). Letting $HD_k$ denote the restriction of $HD$ to $\mathcal{H}_a^k$, the CAR relations yield their matricial representations.

**Proposition 31** *The matricial representations of the operators $HD_k(v_1, v_2)$ in the basis $\mathcal{H}_a^k$ are given by:*

$$HD_1(v_1, v_2) = \begin{pmatrix} v_1 & 0 & -t & 0 \\ 0 & v_1 & 0 & -t \\ -t & 0 & v_2 & 0 \\ 0 & -t & 0 & v_2 \end{pmatrix}$$

$$HD_2(v_1, v_2) = \begin{pmatrix} U_1 + 2v_1 & 0 & -t & t & 0 & 0 \\ 0 & v_1 + v_2 & 0 & 0 & 0 & 0 \\ -t & 0 & v_1 + v_2 & 0 & 0 & -t \\ t & 0 & 0 & v_1 + v_2 & 0 & t \\ 0 & 0 & 0 & 0 & v_1 + v_2 & 0 \\ 0 & 0 & -t & t & 0 & U_2 + 2v_2 \end{pmatrix}$$

$$HD_3(v_1, v_2) = \begin{pmatrix} U_1 + 2v_1 + v_2 & 0 & t & 0 \\ 0 & U_1 + 2v_1 + v_2 & 0 & t \\ t & 0 & U_2 + 2v_2 + v_1 & 0 \\ 0 & t & 0 & U_2 + 2v_2 + v_1 \end{pmatrix}$$

$$HD_4(v_1, v_2) = (2v_1 + 2v_2 + U_1 + U_2).$$
**Definition 32** For a given potential $(v_1, v_2) \in \mathbb{R}^2$ and an integer $k$ the lowest eigenvalue of $HD_k(v_1, v_2)$ (i.e. the k-electronic ground-state energy for the Hubbard dimer) is written $\epsilon_k(v_1, v_2)$.

Given $n = (n_1, n_2) \in \mathbb{R}^2$, the k-electronic universal functional $f_k$ and right k-centered universal functional $f_k^\xi$ for this model are defined by duality with the corresponding ground-state energies:

$$f_k(n) = \sup_{(v_1, v_2) \in \mathbb{R}^2} \{\epsilon_k(v_1, v_2) - n_1 v_1 - n_2 v_2\},$$

$$f_k^\xi(n) = \sup_{(v_1, v_2) \in \mathbb{R}^2} \{w_k^\xi \epsilon_k(v_1, v_2) + \xi \epsilon_{k+1}(v_1, v_2) - n_1 v_1 - n_2 v_2\}.$$

Explicit expressions of the eigenvalues $\epsilon_k(v_1, v_2)$ are known and are recalled in Appendix B. From these analytical expressions and an asymptotic study we derive several results relative to the functionals, gathered in the following result.

**Proposition 33** 1. $\text{Dom}(f_2) = \{(n_1, n_2) \in [0, 2]^2 | n_1 + n_2 = 2\}$.

2. $\text{Dom}(f_3) = \{(n_1, n_2) \in [1, 2]^2 | n_1 + n_2 = 3\}$. Furthermore, an explicit expression of $f_3$ can be found:

$$f_3(n_1, n_2) = -t\sqrt{1 - 4(n_1 - 3/2)^2} + U, \forall n_1 + n_2 = 3.$$

**Proof.** See appendix B. $\square$

**Corollary 34** The domain of the right 2-centered electronic functional $f_2^\xi$ is

$$\text{Dom}(f_2^\xi) = \{(n_1, n_2) \in [\xi, 2 - \xi] \mid n_1 + n_2 = 2\}$$

**Proof.** By definition, $f_2^\xi$ reads:

$$f_2^\xi(n_1, n_2) = \sup_{(v_1, v_2) \in \mathbb{R}^2} \{w_2^\xi \epsilon_2(v_1, v_2) + \xi \epsilon_3(v_1, v_2) - n_1 v_1 - n_2 v_2\}. \quad (39)$$

The functions $w_2^\xi \epsilon_2$ and $\xi \epsilon_3$ being proper convex and defined on all $\mathbb{R}^2$, the convex conjugate of their sum (which is precisely $f_2^\xi$ according to (39)) is equal to the infimal convolution of their convex conjugates (see Theorem 16.4 of [18]). Explicitly:

$$f_2^\xi(n_1, n_2) = \inf_{(\rho_1, \rho_2) \in \mathbb{R}^2} \{(w_2^\xi \epsilon_2)^\vee(n_1 - \rho_1, n_2 - \rho_2) + (\xi \epsilon_3)^\vee(\rho_1, \rho_2)\}. \quad (40)$$

Note that multiplying a function by a constant modifies its convex conjugate in the following fashion: $(af)^\vee(x) = af^\vee(x/a)$, so that the equality (40) yields:

$$f_2^\xi(n_1, n_2) = \inf_{(\rho_1, \rho_2)} \left\{ w_2^\xi f_2 \left( \frac{n_1 - \rho_1}{w_2^\xi}, \frac{n_2 - \rho_2}{w_2^\xi} \right) + \xi f_3 \left( \frac{\rho_1}{\xi}, \frac{\rho_2}{\xi} \right) \right\}.$$
<figure><img src="image_4.png" /><figcaption>Figure 1: The functional $n \mapsto f_2^\xi(n, 2-n)$ for various choices of $\xi$.</figcaption></figure>

Having expressed $f_2^\xi$ as an infimal convolution of two convex functionals this, along with Proposition (33), yields the result:

$$ \begin{aligned} \mathrm{Dom}(f_2^\xi) &= \mathrm{Dom}(w_2^\xi f_2(\cdot/w_2^\xi)) + \mathrm{Dom}(\xi f_3(\cdot/\xi)) \\ &= \left(1 - \frac{3}{2}\xi\right) \mathrm{Dom}(f_2) + \xi \mathrm{Dom}(f_3) \\ &= \{(n_1, n_2) \in [\xi, 2-\xi] \mid n_1 + n_2 = 2\} \end{aligned} $$

□

This result shows that, in opposition to physical $N$-centered eDFT, $N$-centered eDFT applied to the Hubbard dimer does not yield a functional who's domain is independent of $\xi$. A graphical representation of this assertion can be found in figure 1, where we plotted the functional $n \mapsto f_2^\xi(n, 2-n)$ for various weights. As such, this model does not seem adequate in order to illustrate $N$-centered eDFT 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 dependence does not lie in the “integral condition” $n_1 + n_2 = 2$ (analog to $\int \rho = 2$) but rather seems to come from the fact that the $3$-electronic functional itself diverges from the physical model: the condition $n \ge 1$ on the occupation numbers for this functional's domain has no analog for physical $3$-electronic densities $\mathcal{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 eDFT, 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 \ge N$.
This is illustrated for the 1-centered functional associated to the 2-site model
in figure 2.

<figure><img src="image_5.png" /><figcaption>Figure 2: The functional $n \mapsto f_1^\xi(n, 1-n)$ for various choices of $\xi$.</figcaption></figure>

# 5 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 regular 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 weight-dependence 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-dependence and electronic affinity and that we hope to investigate further in future work.

# Declarations

Data sharing not applicable to this article as no datasets were generated or
analysed during the current study.
A Proof of Proposition 8

This result can be established by following similar steps as the ones used in order to prove its N-electronic analog [16]. For completeness we recall these steps, while insisting on the technicalities arising from the use of states in $S_N^\xi$ rather than in $S_N$.

**Proof.** First note that for all $\Gamma \in S_N^\xi$,

$$ \mathrm{Tr}(H_0\Gamma) = \mathrm{Tr}((H_0+1)\Gamma) - 1 + \frac{\xi}{N}. $$

Thus, it suffices to prove that there exists $\Gamma \in S_N^\xi$ with $\mathcal{D}_N^\xi(\Gamma) = \rho$ such that

$$ \mathrm{Tr}((H_0 + 1)\Gamma) \le \liminf_{n \to \infty} \left\{ \inf \left\{ \mathrm{Tr}((H_0 + 1)\Gamma) \mid \mathcal{D}_N^\xi(\Gamma) = \rho_n, \Gamma \in S_N^\xi \right\} \right\}, \quad (41) $$

which prompts the introduction of several notations:

$$ \begin{aligned} h^2 &= H_0 + 1 \\ h_k^2 &= H_0^k + 1, \quad k \in \{N, N+1\} \\ g_n &= \inf \left\{ \mathrm{Tr}((H_0+1)\Gamma) \mid \Gamma \to \rho_n, \Gamma \in S_N^\xi \right\}, \quad n \in \mathbb{N}. \end{aligned} $$

Note that the operators $h_k^{-1}$ being bounded, the same is true for $h^{-1}$.

Without loss we can assume that for all integer $n$, $g_n$ is finite and that $g = \lim_n g_n$ exists. By definition of $g_n$, for all integer $n$ there exists a right $N$-centered ensemble $\Gamma_n \in S_N^\xi$ such that:

$$ \begin{aligned} \mathrm{Tr}(h\Gamma_n h) &= \mathrm{Tr}(h^2\Gamma_n) \le g_n + \frac{1}{n} \\ \mathcal{D}_N^\xi(\Gamma_n) &= \rho_n. \end{aligned} $$

Setting $y_n = h\Gamma_n h$, boundedness of this sequence in the trace norm and the Banach-Alaoglu theorem imply the existence of a subsequence, still written $y_n$, and of some trace-class operator $y$ such that for all compact operator $A$:

$$ \lim_n \mathrm{Tr}(y_n A) = \mathrm{Tr}(y A). \quad (42) $$

Positivity of $y_n$ implies that of $y$ (consider $A = \langle \cdot, e_i \rangle e_i$ for a given Hilbert basis ($e_i$) of $L^2(\mathbb{R}^{3N}) \oplus L^2(\mathbb{R}^{3(N+1)})$). Therefore, weak-$\star$ lower-semicontinuity of the trace norm implies:

$$ \liminf_n \mathrm{Tr}(y_n) \ge \mathrm{Tr}(y). \quad (43) $$

Setting $\Gamma = h^{-1}yh^{-1}$, if we prove that $\Gamma$ is in $S_N^\xi$ and that $\mathcal{D}_N^\xi(\Gamma) = \rho$ then equation (43) will yield the desired result (41). Since $\Gamma_n = h^{-1}y_n h^{-1}$ is in $S_N^\xi$ there exists two states $\Gamma_n^N, \Gamma_n^{N+1}$ in $S_N, S(N+1)$ such that:

$$ \Gamma_n = w_N^\xi \Gamma_n^N \oplus \xi \Gamma_n^{N+1}. $$
The operator sequences $y_n^N := h_N \Gamma_n^N h_N$ and $y_n^{N+1} := h_N \Gamma_n^{N+1} 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^{N+1}$) converging to positive trace-class operators $y^N, y^{N+1}$ in the weak-⋆ topology. This implies the convergence of $(y_n)_n = (w_N^\xi y_n^N \oplus \xi y_n^{N+1})_n$ to $(w_N^\xi y^N \oplus \xi y^{N+1})$ in the weak-⋆ topology and, by unicity of the limit:

$$y = w_N^\xi y^N \oplus \xi y^{N+1},$$ $$\Gamma = w_N^\xi h_N^{-1} y^N h_N^{-1} \oplus \xi h_{N+1}^{-1} y^{N+1} h_{N+1}^{-1}.$$

Let us write for $k \in \{N, N+1\}$:

$$\Gamma^k = h_k^{-1} y^k h_k^{-1}, \quad \rho^k = \mathcal{D}_k^\xi(\Gamma^k), \quad \rho_n^k = \mathcal{D}_k(\Gamma_n^k).$$

The exact same reasoning as the one proposed in the original proof shows that $\rho_n^k$ converges to $\rho^k$ weakly in $L^1$. In particular, the weighted sum $w_N^\xi \rho_n^N + \xi \rho_n^{N+1}$ converges to $w_N^\xi \rho^N + \xi \rho^{N+1}$ weakly in $L^1$. All that is left to note is that

$$\rho_n = \mathcal{D}_N^\xi(\Gamma_n) = w_N^\xi \rho_n^N + \xi \rho_n^{N+1}$$ $$\mathcal{D}_N^\xi(\Gamma) = w_N^\xi \rho^N + \xi \rho^{N+1}$$

and to recall that $\rho_n$ converges to $\rho$ weakly in $L^1$ by hypothesis, which enables us to conclude that $\mathcal{D}_N^\xi(\Gamma) = \rho$. $\square$

## B Hubbard dimer calculations

Explicit expressions of the lowest eigenvalues of $HD_k(v_1, v_2)$ can be found in multiple sources [5], so that we will not further discuss them. On the other hand, the expressions derived in these sources are given for potentials satisfying $v_1 = -v_2$, as functions of $\delta v = v_1 - v_2$. For the sake of completeness we show how this special case is in fact sufficient to yield the lowest eigenvalues in the general case.

Let us first note that the columns of $HD_2(v_1, v_2)$ containing only $v_1 + v_2$ do not impact on $\epsilon_2$.

**Lemma 35** For all $v_1, v_2 \in \mathbb{R}$, $\epsilon_2(v_1, v_2)$ is the lowest eigenvalue of

$$HD_2^r(v_1, v_2) = \begin{pmatrix} U_1 + 2v_1 & -t & t & 0 \\ -t & v_1 + v_2 & 0 & -t \\ t & 0 & v_1 + v_2 & t \\ 0 & -t & t & U_2 + 2v_2 \end{pmatrix}.$$

**Proof.** It suffices to prove that the lowest eigenvalue of $HD_2^r$ 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:

$$\inf_{x \in \mathbb{R}^4} \frac{x^T HD_2^r x}{x^T x} \le v_1 + v_2.$$
This is easily seen by considering $x = (0, 1, 1, 0)$ for which we obtain :

$$ \frac{x^T H D_2 x}{x^T x} = v_1 + v_2. $$

□

**Lemma 36** For all $v_1, v_2 \in \mathbb{R}$,

$$ \epsilon_2(v_1, v_2) = \epsilon_2 \left( \frac{\delta v}{2}, -\frac{\delta v}{2} \right) + v_1 + v_2 $$

$$ \epsilon_3(v_1, v_2) = \epsilon_3 \left( \frac{\delta v}{2}, -\frac{\delta v}{2} \right) + \frac{3}{2}(v_1 + v_2) \quad (44) $$

With the previous lemma now in hand, we recall the expressions of $E_k(\delta v) = \epsilon_k(\delta v/2, -\delta v/2)$ (see for example [3] to understand where these come from). We limit ourselves in what follows to the particular case $U = U_1 = U_2$.

**Proposition 37** The 2-electronic and 3-electronic energies of the Hubbard dimer are given by:

$$ E_3(\delta v) - \frac{1}{2}\sqrt{4t^2 + \delta v^2} + U $$

$$ E_2(\delta v) = \frac{4t}{3} \left( u - w(\delta v) \sin \left( \theta(\delta v) + \frac{\pi}{6} \right) \right), $$

where

$$ u = U/2t $$

$$ w(\delta v) = \sqrt{3(1 + \nu^2) + u^2} $$

$$ \nu = \delta v/2t $$

$$ \cos(3\theta(\delta v)) = \frac{(9(\nu^2 - 1/2) - u^2)u}{w(\delta v)^3}. $$

With the analytic expressions of $\epsilon_k$ now at our disposal, we can prove proposition (33).

**Proof of point 1.:** A simple limit study yields:

$$ \lim_{+\infty} E'_2(\delta v) = -1 $$

$$ \lim_{-\infty} E'_2(\delta v) = 1. $$

Since $E'_2$ is increasing (by concavity of $E_2$), this also yields:

$$ -1 \le E'_2 \le 1, $$
and by continuity $E'_2(\mathbb{R}) =]-1, 1[$. It follows that if $(n_1, n_2) \in \{(n_1, n_2) \in ]0, 2[|n_1 + n_2 = 2\}$, we can choose $v_1, v_2$ such that the Jacobian

$$ \begin{pmatrix} E'_2(v_1 - v_2) + 1 - n_1 \\ -E'_2(v_1 - v_2) - 1 + n_1 \end{pmatrix} $$

of the concave function we seek to maximize is zero, so that $f_2(n_1, n_2)$ is finite.

If $(n_1, n_2) = (0, 2)$,

$$ E_2(v_1 - v_2) + v_1(1 - n_1) + v_2(1 - n_2) = E_2(v_1 - v_2) + v_1 - v_2. $$

Since $E'_2(\delta v) \ge -1$, the function $\delta v \mapsto E_2(\delta v) + \delta v$ is increasing and

$$ \begin{aligned} \sup_{\delta v \in \mathbb{R}} E_2(\delta v) + \delta v &= \lim_{\delta v \to \infty} E_2(\delta v) + \delta v \\ &= \lim_{\delta v \to \infty} \frac{4t}{3} \left( u - \sqrt{3(1 + \delta v^2/4t^2) + u^2} \sin \left( \theta(\delta v) + \frac{\pi}{6} \right) \right) + \delta v \\ &= \lim_{\nu \to \infty} \frac{4t}{3} u + 2t\nu - \frac{4t}{3} \sqrt{3(1 + \nu^2) + u^2} \sin \left( \theta(2t\nu) + \frac{\pi}{6} \right). \end{aligned} $$

A small asymptotic study is needed. First let us note that

$$ \begin{aligned} u(9\nu^2 - 9/2 - u^2)(3\nu^2 + 3 + u^2)^{-3/2} &= u(9\nu^2 - 9/2 - u^2) \frac{\left(1 + \frac{3+u^2}{3\nu^2}\right)^{-3/2}}{3^{3/2}v^3} \\ &=_{+\infty} u \left( \frac{\sqrt{3}}{\nu} + o(1/\nu^2) \right) \left( 1 - \frac{3}{2} \frac{3+u^2}{3\nu^2} + o(1/\nu^2) \right) \\ &= \frac{\sqrt{3}u}{v} + o(1/\nu^2). \end{aligned} $$

Thus,

$$ \begin{aligned} \theta(2t\nu) &= \frac{1}{3}\operatorname{acos}\left[u(9\nu^2 - 9/2 - u^2)(3\nu^2 + 3 + u^2)^{-3/2}\right] \\ &=_{+\infty} \frac{1}{3}\left(\frac{\pi}{2} - \frac{\sqrt{3}u}{v} + o(1/\nu^2)\right). \end{aligned} $$

Since

$$ \begin{aligned} \sin x &=_{x \to \pi/3} \frac{\sqrt{3}}{2} + \frac{1}{2}(x - \pi/3) + o(x), \\ \sin \left(\theta(2t\nu) + \frac{\pi}{6}\right) &=_{+\infty} \frac{\sqrt{3}}{2} - \frac{u}{2\sqrt{3}\nu} + o(1/\nu). \end{aligned} $$

We then obtain:

$$ \begin{aligned} 2t\nu - \frac{4t}{3}\sqrt{3(1+\nu^2)+u^2}\sin\left(\theta(2t\nu)+\frac{\pi}{6}\right) &= 2t\nu - 2t\sqrt{1+\nu^2+\frac{u^2}{3}}\left(1-\frac{u}{3\nu}+o(1/\nu)\right) \\ &= 2t\nu\left[1-\sqrt{1+\frac{3+u^2}{3\nu^2}}\left(1-\frac{1}{3\nu}+o(1/\nu)\right)\right] \\ &= \frac{2t}{3}u+o(1). \end{aligned} $$
and may then conclude that

$$f(n_1, n_2) = \sup E_2(\delta v) = \frac{6t}{3}u = U < \infty.$$

For $(n_1, n_2) = (2, 0)$, note that

$$E_2(\delta v) = E_2(-\delta v),$$

so that

$$\begin{aligned} f(n_1, n_2) &= \sup_{\delta v} E_2(\delta v) - \delta v \\ &= \sup_{\delta v} E_2(\delta v) + \delta v \\ &= \frac{6t}{3}u. \end{aligned}$$

Conversely, considering $(n_1, n_2)$ such that $n_1 \notin [0, 2]$, we get:

$$\begin{aligned} f(n_1, n_2) &= \sup_{v_2} \left\{ \sup_{v_1} \left\{ E_2(v_1 - v_2) + (1 - n_1)v_1 \right\} + (1 - n_2)v_2 \right\} \\ &= \sup_{v_2} \left\{ \sup_{\delta v} \left\{ E_2(\delta v) + (1 - n_1)\delta v \right\} + (2 - n_1 - n_2)v_2 \right\}. \end{aligned}$$

If $n_1 < 0$, $\delta v \mapsto E_2(\delta v) + (1 - n_1)\delta v$ is increasing and

$$\lim_{\delta v \to \infty} E_2(\delta v) + (1 - n_1)\delta v = \lim_{\delta v \to \infty} E_2(\delta v) + \delta v - n_1\delta v = +\infty.$$

If $n_1 > 2$, $\delta v \mapsto E_2(\delta v) + (1 - n_1)\delta v$ is decreasing and

$$\lim_{\delta v \to -\infty} E_2(\delta v) + (1 - n_1)\delta v = \lim_{\delta v \to -\infty} E_2(\delta v) - \delta v + (2 - n_1)\delta v = +\infty.$$

The $n_2 \notin [0, 2]$ case is treated in the same way.

Finally, let us consider $(n_1, n_2)$ such that $n_1 + n_2 \neq 2$. In this case,

$$\begin{aligned} f(n_1, n_2) &\ge \sup_{v_1=v_2} \{E_2(v_1-v_2) + v_1(1-n_1) + v_2(1-n_2)\} \\ &= \sup_v \{E_2(0) + v(2-n_1-n_2)\} \\ &= +\infty. \end{aligned}$$

**Proof of point 2:** Since

$$\frac{dE^3}{d\delta v}(\delta v) = - \frac{\delta v}{2\sqrt{4t^2 + \delta v^2}},$$

it is clear that

$$\frac{dE^3}{d\delta v}(\mathbb{R}) = ]-1/2, 1/2[.$$
Thus, if $(n_1, n_2) \in \{(n_1, n_2) \in ]1, 2[ | n_1 + n_2 = 3\}$, we can choose $v_1, v_2$ such that

$$ (E^3)'(v_1 - v_2) = n_1 - 3/2 \in ]-1/2, 1/2[, \quad (45) $$

so that the Jacobian

$$ \begin{pmatrix} (E^3)'(v_1 - v_2) + 3/2 - n_1 \\ -(E^3)'(v_1 - v_2) - 3/2 + n_1 \end{pmatrix} $$

of the function to minimize equals 0. Equation 45 has an explicit solution:

$$ \text{if } n_1 > 3/2, \quad v_1 - v_2 = \frac{-4t(n_1 - 3/2)}{\sqrt{1 - 4(n_1 - 3/2)^2}} $$

$$ \text{if } n_1 = 3/2, \quad v_1 - v_2 = 0 $$

$$ \text{if } n_1 < 3/2, \quad v_1 - v_2 = \frac{-4t(n_1 - 3/2)}{\sqrt{1 - 4(n_1 - 3/2)^2}}, $$

so that

$$ f^3(n_1, n_2) = \begin{cases} -\frac{1}{2}\sqrt{4t^2 \left(1 + \frac{4(n_1 - 3/2)^2}{1 - 4(n_1 - 3/2)^2}\right)} + U + 4t \frac{(n_1 - 3/2)^2}{\sqrt{1 - 4(n_1 - 3/2)^2}}, & n_1 > 3/2 \\ -t + U, & n_1 = 3/2, \\ -\frac{1}{2}\sqrt{4t^2 \left(1 + \frac{4(n_1 - 3/2)^2}{1 - 4(n_1 - 3/2)^2}\right)} + U + 4t \frac{(n_1 - 3/2)^2}{\sqrt{1 - 4(n_1 - 3/2)^2}}, & n_1 < 3/2 \end{cases} $$

which can be further simplified:

$$ f^3(n_1, n_2) = -t\sqrt{1 - 4(n_1 - 3/2)^2} + U \quad (46) $$

If $(n_1, n_2) = (1, 2)$,

$$ \begin{aligned} f^3(n_1, n_2) &= \sup_{\delta v} U - \frac{1}{2}\sqrt{4t^2 + \delta v^2} - \frac{1}{2}\delta v \\ &= \lim_{\delta v \to -\infty} U - \frac{1}{2}\sqrt{4t^2 + \delta v^2} - \frac{1}{2}\delta v = U. \end{aligned} $$

Similarly,

$$ f^3(2, 1) = \lim_{\delta v \to +\infty} U - \frac{1}{2}\sqrt{4t^2 + \delta v^2} + \frac{1}{2}\delta v = U. $$

Note that the expression 46 holds true for these situations also.

## References

[1] Un metodo statistico per la determinazione di alcune proprieta dell'atomo, Rendiconti Academia Dei Lincei 6 (1927), 602-607.
[2] Maicon Marques Alves and Benar Fux Svaiter, *Moreau-Yosida Regularization of Maximal Monotone Operators of Type (D)*, Set-Valued and Variational Analysis **19** (2011), no. 1, 97–106.

[3] Heinz H. Bauschke and Patrick L. Combettes, *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*, vol. 77, CMS Books in Mathematics, no. 23, Springer International Publishing, Cham, 2017.

[4] ————, *Correction to: Convex Analysis and Monotone Operator Theory in Hilbert Spaces*, 2019, pp. C1–C4.

[5] D. J. Carrascal, J. Ferrer, J. C. Smith, and K. Burke, *The Hubbard dimer: A density functional case study of a many-body problem*, Journal of Physics Condensed Matter **27** (2015), no. 39.

[6] Filip Cernatic, Bruno Senjean, Vincent Robert, and Emmanuel Fromager, *Ensemble density functional theory of neutral and charged excitations*, Topics in Current Chemistry **380** (2022).

[7] Ivar Ekeland and Roger Temam, *Convex Analysis and Variational Problems (Ivar Ekeland and Roger Temam)*, SIAM Review **20** (1978), no. 1, 192–194.

[8] Michael Filatov, Seunghoon Lee, and Chong Ju Choi, *Description of sudden polarization in the excited electronic states with an ensemble density functional theory method*, Journal of chemical theory and computation **17** (2021).

[9] Tim Gould, Leeor Kronik, and Stefano Pittalis, *Double excitations in molecules from ensemble density functionals: Theory and approximations*, (2021).

[10] E. K. U. Gross, L. N. Oliveira, and W. Kohn, *Density-functional theory for ensembles of fractionally occupied states. i. basic formalism*, Phys. Rev. A **37** (1988), 2809–2820.

[11] M. J. P. Hodgson, J. Wetherell, and Emmanuel Fromager, *Exact exchange-correlation potentials for calculating the fundamental gap with a fixed number of electrons*, Phys. Rev. A **103** (2021), 012806.

[12] P. Hohenberg and W. Kohn, *Inhomogeneous electron gas*, Phys. Rev. **136** (1964), B864–B871.

[13] Simen Kvaal, Ulf Ekström, Andrew M. Teale, and Trygve Helgaker, *Dif-ferentiable but exact formulation of density-functional theory*, Journal of Chemical Physics **140** (2014), no. 18, 1–13.

[14] Mathieu Lewin, *Geometric methods for nonlinear many-body quantum systems*, J. Funct. Anal **260** (2011), 3535–3595.

[15] ————, *Semi-classical limit of the Levy Lieb functional in Density Functional Theory*, Comptes Rendus Mathematique **356** (2018), no. 4, 449–455.
[16] Elliott H Lieb, *Density functionals for coulomb systems*, International Journal of Quantum Chemistry **24** (1983), no. 3, 243–277.

[17] Robert G. Parr, Mel Levy, and Jose L. Balduz, *Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy, Physical Review Letters* **49** (1982).

[18] Ralph Tyrell Rockafellar, *Convex analysis*, Princeton University Press (1970).

[19] Bruno Senjean and Emmanuel Fromager, *Unified formulation of fundamental and optical gap problems in density-functional theory for ensembles*, Phys. Rev. A **98** (2018), 022513.

[20] ————, *N-centered ensemble density-functional theory for open systems*, International Journal of Quantum Chemistry **120** (2020), no. 21.

[21] L. H. Thomas, *The calculation of atomic fields*, Mathematical Proceedings of the Cambridge Philosophical Society **23** (1927), no. 5, 542–548.