\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amsmath,amssymb,amsfonts,physics,mhchem}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}

\usepackage[
	colorlinks=true,
    citecolor=blue,
    breaklinks=true
	]{hyperref}
\urlstyle{same}

\newcommand{\alert}[1]{\textcolor{red}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\cloclo}[1]{\textcolor{purple}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashCM}[1]{\textcolor{purple}{\sout{#1}}}

\newcommand{\cmark}{\color{green}{\text{\ding{51}}}}
\newcommand{\xmark}{\color{red}{\text{\ding{55}}}}

%useful stuff
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\ra}{\rightarrow}
\newcommand{\cdash}{\multicolumn{1}{c}{---}}
\newcommand{\mc}{\multicolumn}
\newcommand{\mr}{\multirow}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\la}{\lambda}
\newcommand{\si}{\sigma}
\newcommand{\cD}{\mathcal{D}}

% operators
\newcommand{\hHc}{\Hat{h}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hWee}{\Hat{W}_\text{ee}}
\newcommand{\hGam}[1]{\Hat{\Gamma}^{#1}}
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\hVext}{\Hat{V}_\text{ext}}
\newcommand{\vext}{v_\text{ext}}

% functionals, potentials, densities, etc
\newcommand{\F}[2]{F_{#1}^{#2}}
\newcommand{\Ts}[1]{T_\text{s}^{#1}}
\newcommand{\eps}[2]{\varepsilon_{#1}^{#2}}
\newcommand{\Eps}[2]{\mathcal{E}_{#1}^{#2}}
\newcommand{\e}[2]{\epsilon_{#1}^{#2}}
\newcommand{\kin}[2]{t_\text{#1}^{#2}}
\newcommand{\E}[2]{E_{#1}^{#2}}
\newcommand{\bE}[2]{\overline{E}_{#1}^{#2}}
\newcommand{\tE}[2]{\widetilde{E}_{#1}^{#2}}
\newcommand{\be}[2]{\overline{\epsilon}_{#1}^{#2}}
\newcommand{\n}[2]{n_{#1}^{#2}}
\newcommand{\Cx}[1]{C_\text{x}^{#1}}

% energies
\newcommand{\EHF}{E_\text{HF}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\HF}{\text{HF}}
\newcommand{\LDA}{\text{LDA}}
\newcommand{\eLDA}{\text{eLDA}}
\newcommand{\CID}{\text{CID}}
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\ex}{\text{x}}
\newcommand{\co}{\text{c}}
\newcommand{\xc}{\text{xc}}

% matrices
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bw}{\boldsymbol{w}}
\newcommand{\bHc}{\boldsymbol{h}}
\newcommand{\Ex}[2]{\Omega_{#1}^{#2}}
\newcommand{\tEx}[2]{\widetilde{\Omega}_{#1}^{#2}}

% elements
\newcommand{\ew}[1]{w_{#1}}
\newcommand{\eHc}[1]{h_{#1}}
\newcommand{\eJ}[1]{J_{#1}}
\newcommand{\eK}[1]{K_{#1}}
\newcommand{\eF}[2]{F_{#1}^{#2}}
\newcommand{\ON}[2]{f_{#1}^{#2}}
\newcommand{\Det}[2]{\Phi_{#1}^{#2}}

% Numbers
\newcommand{\nEns}{M}
\newcommand{\Nel}{N}
\newcommand{\Norb}{K}

% Ao and MO basis
\newcommand{\MO}[2]{\phi_{#1}^{#2}}
\newcommand{\cMO}[2]{c_{#1}^{#2}}
\newcommand{\AO}[1]{\chi_{#1}}

\newcommand{\RHH}{R_{\ce{H-H}}}

% units
\newcommand{\IneV}[1]{#1 eV}
\newcommand{\InAU}[1]{#1 a.u.}
\newcommand{\InAA}[1]{#1 \AA}
\newcommand{\kcal}{kcal/mol}

\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques, Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LCQ}{Laboratoire de Chimie Quantique, Institut de Chimie, CNRS, Universit\'e de Strasbourg, Strasbourg, France}
\newcommand{\UL}{Instituut-Lorentz, Universiteit Leiden, P.O.~Box 9506, 2300 RA Leiden, The Netherlands}
\newcommand{\VU}{Division of Theoretical Chemistry, Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands}
\begin{document}	

\title{Weight-dependent exchange-correlation functionals for molecules: I. The local-density approximation}

\author{Clotilde \surname{Marut}}
	\affiliation{\LCPQ}
\author{Bruno \surname{Senjean}}
	\affiliation{\UL}
	\affiliation{\VU}
\author{Emmanuel \surname{Fromager}}
	\affiliation{\LCQ}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
	\email{loos@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}

\begin{abstract}
Gross-Oliveira-Kohn (GOK) ensemble density-functional theory (GOK-DFT) is a time-independent formalism which allows to compute individual excitation energies via the derivative of the ensemble energy with respect to the weight of each excited state.
Contrary to the time-dependent version of density-functional theory (TD-DFT), double excitations can be easily computed within GOK-DFT.
However, to take full advantage of this formalism, one must have access to a \textit{weight-dependent} exchange-correlation functional in order to model the infamous derivative discontinuity contribution to the excitation energies.
In the present article, we report a first-rung (\ie, local) weight-dependent exchange-correlation density-functional approximation for atoms and molecules specifically designed for the computation of double excitations within GOK-DFT.
This density-functional approximation for ensembles, based on finite and infinite uniform electron gases, incorporate information about both ground and excited states. 
Its accuracy is illustrated by computing the double excitation in the prototypical H$_2$ molecule.
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%
%%% INTRODUCTION %%%
%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Time-dependent density-functional theory (TD-DFT) has been the dominant force in the calculation of excitation energies of molecular systems in the last two decades.\cite{Casida,Ulrich_2012}
At a relatively low computational cost (at least compared to the other excited-state \textit{ab initio} methods), TD-DFT can provide accurate transition energies for low-lying excited states of organic molecules (see, for example, Ref.~\onlinecite{Dreuw_2005} and references therein).
Importantly, setting up a TD-DFT calculation for a given system is an almost pain-free process from a user perspective as the only (yet essential) input variable is the choice of the so-called exchange-correlation (xc) functional.

Similar to density-functional theory (DFT), \cite{Hohenberg_1964,Kohn_1965,ParrBook} TD-DFT is an in-principle exact theory which formal foundations relie on the Runge-Gross theorem. \cite{Runge_1984}
The Kohn-Sham (KS) formalism of TD-DFT transfers the complexity of the many-body problem to the xc functional thanks to a judicious mapping between a time-dependent non-interacting reference system and its interacting analog which have both the exact same one-electron density.

However, TD-DFT is far from being perfect as, in practice, drastic approximations must be made for the xc functional.
One of its issues actually originates directly from the choice of the xc functional, and more specifically, the possible (not to say likely) substantial variations in the quality of the excitation energies for two different choices of xc functionals.
Moreover, because it is so popular, it has been studied in excruciated details by the community, and some researchers have quickly unveiled various theoretical and practical deficiencies of approximate TD-DFT.

For example, TD-DFT has problems with charge-transfer \cite{Tozer_1999,Dreuw_2003,Sobolewski_2003,Dreuw_2004,Maitra_2017} and Rydberg \cite{Tozer_1998,Tozer_2000,Casida_1998,Casida_2000,Tozer_2003} excited states (the excitation energies are usually drastically underestimated) due to the wrong asymptotic behaviour of the semi-local xc functional.
The development of range-separated hybrids provides an effective solution to this problem. \cite{Tawada_2004,Yanai_2004} 
From a practical point of view, the TD-DFT xc kernel is usually considered as static instead of being frequency dependent.
One key consequence of this so-called adiabatic approximation (based on the assumption that the density varies slowly with time) is that double excitations are completely absent from the TD-DFT spectra. \cite{Levine_2006,Tozer_2000,Elliott_2011}
Although these double excitations are usually experimentally dark (which means they usually cannot be observed in photo-absorption spectroscopy), these states play, indirectly, a key role in many photochemistry mechanisms. \cite{Boggio-Pasqua_2007,Loos_2019}

One possible solution to access double excitations within TD-DFT is provided by spin-flip TD-DFT which describes double excitations as single excitations from the lowest triplet state. \cite{Huix-Rotllant_2010,Krylov_2001,Shao_2003,Wang_2004,Wang_2006,Minezawa_2009}
However, spin contamination might be an issue. \cite{Huix-Rotllant_2010}
In order to go beyond the adiabatic approximation, a dressed TD-DFT approach has been proposed by Maitra and coworkers \cite{Maitra_2004,Cave_2004} (see also Refs.~\onlinecite{Mazur_2009,Mazur_2011,Huix-Rotllant_2011,Elliott_2011,Maitra_2012}).
In this approach the xc kernel is made frequency dependent, which allows to treat doubly-excited states. \cite{Romaniello_2009a,Sangalli_2011,Loos_2019}

Maybe surprisingly, another possible way of accessing double excitations is to resort to a time-\textit{independent} formalism. \cite{Yang_2017,Sagredo_2018,Deur_2019}
DFT for ensembles (eDFT) \cite{Theophilou_1979,Gross_1988a,Gross_1988b,Oliveira_1988} is a viable alternative following such a strategy currently under active development. \cite{Gidopoulos_2002,Franck_2014,Borgoo_2015,Kazaryan_2008,Gould_2013,Gould_2014,Filatov_2015,Filatov_2015b,Filatov_2015c,Gould_2017,Deur_2017,Gould_2018,Gould_2019,Sagredo_2018,Ayers_2018,Deur_2018,Deur_2019,Kraisler_2013,Kraisler_2014,Alam_2016,Alam_2017,Nagy_1998,Nagy_2001,Nagy_2005,Pastorczak_2013,Pastorczak_2014,Pribram-Jones_2014,Yang_2013a,Yang_2014,Yang_2017,Senjean_2015,Senjean_2016,Smith_2016,Senjean_2018}
In the assumption of monotonically decreasing weights, eDFT for excited states has the undeniable advantage to be based on a rigorous variational principle for ground and excited states,  the so-called Gross-Oliveria-Kohn (GOK) variational principle.  \cite{Gross_1988a} 
In short, GOK-DFT (\ie, eDFT for excited states) is the density-based analog of state-averaged wave function methods, and excitation energies can then be easily extracted from the total ensemble energy. \cite{Deur_2019}
Although the formal foundations of GOK-DFT have been set three decades ago, \cite{Gross_1988a,Gross_1988b,Oliveira_1988} its practical developments have been rather slow.
We believe that it is partly due to the lack of accurate approximations for GOK-DFT.
In particular, to the best of our knowledge, an explicitly weight-dependent density-functional approximation for ensemble (eDFA) has never been developed for atoms and molecules.
The present contribution is a first step towards this goal.

When one talks about constructing functionals, the local-density approximation (LDA) is never far away.
The LDA, as we know it, is based on the uniform electron gas (UEG) also known as jellium, an hypothetical infinite substance where an infinite number of electrons ``bathe'' in a (uniform) positively-charged jelly. \cite{Loos_2016}
Although the Hohenberg-Kohn theorems \cite{Hohenberg_1964} are here to provide firm theoretical grounds to DFT, modern KS-DFT rests largely on the presumed similarity between this hypothetical UEG and the electronic behaviour in a real system. \cite{Kohn_1965}
However, Loos and Gill have recently shown that there exists other UEGs which contain finite numbers of electrons (more like in a molecule), \cite{Loos_2011b,Gill_2012} and that they can be exploited to construct LDA functionals. \cite{Loos_2014a,Loos_2014b,Loos_2017a}
Electrons restricted to remain on the surface of a $\cD$-sphere (where $\cD$ is the dimensionality of the surface of the sphere) are an example of finite UEGs (FUEGs). \cite{Loos_2011b}
Here, we combine these FUEGs with the usual infinite UEG (IUEG) to construct a weigh-dependent LDA functional for ensembles (eLDA).
The present eLDA functional is specifically designed to compute double excitations within GOK-DFT, and it automatically incorporates the infamous derivative discontinuity contribution to the excitation energies through its explicit ensemble weight dependence. \cite{Levy_1995, Perdew_1983}

The paper is organised as follows.
In Sec.~\ref{sec:theo}, the theory behind GOK-DFT is presented.
Section \ref{sec:func} provides details about the construction of the weight-dependent xc LDA functional.
The results of our calculations for the prototypical \ce{H2} molecule are reported and discussed in Sec.~\ref{sec:resdis}. 
Finally, we draw our conclusions in Sec.~\ref{sec:ccl}. 
Unless otherwise stated, atomic units are used throughout.

%%%%%%%%%%%%%%%%%%%%
%%% THEORY %%%
%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theo}
As mentioned above, eDFT for excited states is based on the GOK variational principle \cite{Gross_1988a} which states that a variational principle holds for the ensemble energy
\begin{equation}
	\E{}{\bw} = \sum_{I=0}^{\nEns-1} \ew{I} \E{}{(I)}
\end{equation}
built from an ensemble of $\nEns$ electronic states with individual energies $\E{}{(0)} \le \ldots \le \E{}{(\nEns-1)}$, and (normalised) monotonically decreasing weights $\bw = (\ew{0},\ldots,\ew{M-1})$, \ie, $\sum_{I=0}^{\nEns-1} \ew{I} = 1$, and $\ew{0} \ge \ldots \ge \ew{\nEns-1}$.
Multiplet degeneracies can be easily handled by assigning the same weight to the degenerate states.

One of the key feature of GOK-DFT in the present context is that one can easily extract individual excitation energies from the ensemble energy via differentiation with respect to individual weights:
\begin{equation}
	\pdv{\E{}{\bw}}{\ew{I}} = \E{}{(I)} - \E{}{(0)} = \Ex{}{(I)},
\end{equation}
where the weights are normalised by setting $\ew{0} = 1 - \sum_{I \ne 0} \ew{I}$.

In GOK-DFT, one defines a universal (weight-dependent) ensemble functional $\F{}{\bw}[\n{}{}]$ such that 
\begin{equation}
\label{eq:Ew-GOK}
	\E{}{\bw} = \min_{\n{}{}} \qty{ \F{}{\bw}[\n{}{}] + \int \vext(\br{}) \n{}{}(\br{}) d\br{} },
\end{equation}
where $\vext(\br{})$ is the external potential.
In the KS formulation of GOK-DFT, the universal ensemble functional (the weight-dependent analog of the Hohenberg-Kohn universal functional for ensembles) is decomposed as
\begin{equation}
	\F{}{\bw}[\n{}{}] 
	= \Ts{\bw}[\n{}{}] + \E{\Hxc}{\bw}[\n{}{}]
	= \Tr[ \hGam{\bw} \hT ] + \Tr[ \hGam{\bw} \hWee ],
\end{equation}
where $\hT$ and $\hWee$ are the kinetic and electron-electron interaction potential operators, respectively, $\Ts{\bw}[\n{}{}]$ is the noninteracting ensemble kinetic energy functional,
\begin{equation}
	\hGam{\bw} = \sum_{I=0}^{\nEns} \ew{I} \dyad{\Det{I}{\bw}}
\end{equation}
is the density matrix operator, $\Det{I}{\bw}$ are single-determinant wave functions (or configuration state functions) built with KS orbitals $\MO{p}{\bw}(\br{})$, and 
\begin{equation}
\label{eq:exc_def}
\begin{split}
	\E{\Hxc}{\bw}[\n{}{}] 
	& = \E{\Ha}{}[\n{}{}] + \E{\xc}{\bw}[\n{}{}]
	\\
	& = \frac{1}{2} \iint \frac{\n{}{}(\br{}) \n{}{}(\br{}')}{\abs{\br{}-\br{}'}} d\br{} d\br{}' 
	+ \int \e{\xc}{\bw}(\n{}{}(\br{})) \n{}{}(\br{}) d\br{}.
\end{split}
\end{equation}
is the ensemble Hartree-exchange-correlation (Hxc) functional.
Note that the weight-independent Hartree functional $\E{\Ha}{}[\n{}{}]$ causes the infamous ghost-interaction error (GIE) \cite{Gidopoulos_2002, Pastorczak_2014, Alam_2016, Alam_2017, Gould_2017} in eDFT, which is supposed to be cancelled by the weight-dependent xc functional $\E{\xc}{\bw}[\n{}{}]$.

From the GOK-DFT ensemble energy expression in Eq.~\eqref{eq:Ew-GOK}, we obtain \cite{Gross_1988b,Deur_2019}
\begin{equation}
	\label{eq:dEdw}
	\pdv{\E{}{\bw}}{\ew{I}} 
	= \E{}{(I)} - \E{}{(0)} 
	= \Eps{I}{\bw} - \Eps{0}{\bw} + \left. \pdv{\E{\xc}{\bw}[\n{}{}]}{\ew{I}} \right|_{\n{}{} = \n{}{\bw}(\br{})},
\end{equation}
where 
\begin{align}
	\label{eq:nw}
	\n{}{\bw}(\br{}) & = \sum_{I=0}^{\nEns-1} \ew{I} \n{}{(I)}(\br{}),
	\\
	\label{eq:nI}
	\n{}{(I)}(\br{}) & = \sum_{p}^{\Norb} \ON{p}{(I)} [\MO{p}{\bw}(\br{})]^2
\end{align}
are the ensemble and individual one-electron densities, respectively, 
\begin{equation}
\label{eq:KS-energy}
	\Eps{I}{\bw} = \sum_{p}^{\Norb} \ON{p}{(I)} \eps{p}{\bw}
\end{equation}
is the weight-dependent KS energy, and $\eps{p}{\bw}$ is the KS orbital energy associated with $\MO{p}{\bw}(\br{})$ [$\ON{p}{(I)}$ being its occupancy for the state $I$] given by the ensemble KS equation
\begin{equation}
\label{eq:eKS}
	\qty( \hHc(\br{}) + \fdv{\E{\Hxc}{\bw}[\n{}{}]}{\n{}{}(\br{})}) \MO{p}{\bw}(\br{}) = \eps{p}{\bw} \MO{p}{\bw}(\br{}),
\end{equation}
where $\hHc(\br{}) = -\nabla^2/2  + \vext(\br{})$, and
\begin{equation}
\begin{split}
	\fdv{\E{\Hxc}{\bw}[\n{}{}]}{\n{}{}(\br{})} 
	& = \fdv{\E{\Ha}{\bw}[\n{}{}]}{\n{}{}(\br{})} + \fdv{\E{\xc}{\bw}[\n{}{}]}{\n{}{}(\br{})}
	\\
	& = \frac{1}{2} \int \frac{\n{}{}(\br{}')}{\abs{\br{}-\br{}'}} d\br{}' 
	+ \left. \pdv{\e{\xc}{\bw{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{}(\br{})} \n{}{}(\br{}) + \e{\xc}{\bw{}}(\n{}{}(\br{}))
\end{split}
\end{equation}
is the Hxc potential.
Equation \eqref{eq:dEdw} is our working equation for computing excitation energies from a practical point of view.
Note that, although we have dropped the weight-dependency in the individual densities $\n{}{(I)}(\br{})$ defined in Eq.~\eqref{eq:nI}, these do not match the \textit{exact} individual-state densities as the non-interacting KS ensemble is expected to reproduce the true interacting ensemble density $\n{}{\bw}(\br{})$ defined in Eq.~\eqref{eq:nw}, and not each individual density.

%%%%%%%%%%%%%%%%%%
%%% FUNCTIONAL %%%
%%%%%%%%%%%%%%%%%%
\section{Functional}
\label{sec:func}
The present work deals with the explicit construction of the (reduced) LDA xc functional $\e{\xc}{\bw}[\n{}{}]$ defined in Eq.~\eqref{eq:exc_def}.
Here, we restrict our study to the case of a two-state ensemble (\ie, $\nEns = 2$) where both the ground state ($I=0$) and the first doubly-excited state ($I=1$) are considered. 
Thus, we have $0 \le \ew{} \le 1/2$.
The generalisation to a larger number of states (in particular the inclusion of the first singly-excited state) is trivial and left for future work.

We adopt the usual decomposition, and write down the weight-dependent xc functional as
\begin{equation}
	\e{\xc}{\ew{}}(\n{}{}) = \e{\ex}{\ew{}}(\n{}{}) + \e{\co}{\ew{}}(\n{}{}),
\end{equation}
where $\e{\ex}{\ew{}}(\n{}{})$ and $\e{\co}{\ew{}}(\n{}{})$ are the weight-dependent exchange and correlation functionals, respectively.
The construction of these two functionals is described below.
Here, we restrict our study to spin-unpolarised systems, \ie, $\n{\uparrow}{} = \n{\downarrow}{} = \n{}{}/2$ (where $\n{\uparrow}{}$ and $\n{\downarrow}{}$ are the spin-up and spin-down electron densities).
Extension to spin-polarised systems will be reported in future work.

To build our weight-dependent xc functional, we propose to consider the singlet ground state and the first singlet doubly-excited state of a two-electron FUEG which consists of two electrons confined to the surface of a 3-sphere (also known as a glome).\cite{Loos_2009a,Loos_2009c,Loos_2010e}
Notably, these two states have the same (uniform) density $\n{}{} = 2/(2\pi^2 R^3)$, where $R$ is the radius of the 3-sphere onto which the electrons are confined.
Indeed, the orbitals for an electron on a 3-sphere of unit radius are the normalised hyperspherical harmonics $Y_{\ell\mu}$, where $\ell$ is the principal quantum number and $\mu$ is a composite index of the remaining two quantum numbers. \cite{AveryBook, Avery_1993}  
As mentioned above, we confine our attention to paramagnetic (or unpolarised) systems, and in particular to the simple two-electron system in which the orbital with $\ell = 0$ is doubly-occupied by one spin-up and one spin-down electron, thus yielding an electron density that is uniform on the 3-sphere.
Note that the present paradigm is equivalent to the IUEG model in the thermodynamic limit. \cite{Loos_2011b}
We refer the interested reader to Refs.~\onlinecite{Loos_2011b,Loos_2017a} for more details about this paradigm.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Weight-dependent exchange functional}
\label{sec:Ex}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The reduced (\ie, per electron) Hartree-Fock (HF) energies for these two states are
\begin{subequations}
\begin{align}
	\e{\HF}{(0)}(\n{}{}) & = \frac{4}{3} \qty(\frac{\n{}{}}{\pi})^{1/3},
	\\
	\e{\HF}{(1)}(\n{}{}) & = \frac{3\pi^{2}}{2} \qty(\frac{\n{}{}}{\pi})^{2/3} + \frac{176}{105} \qty(\frac{\n{}{}}{\pi})^{1/3}.
\end{align}
\end{subequations}
These two energies can be conveniently decomposed as 
\begin{equation}
	\e{\HF}{(I)}(\n{}{}) = \kin{s}{(I)}(\n{}{}) + \e{\Ha}{(I)}(\n{}{}) + \e{\ex}{(I)}(\n{}{}),
\end{equation}
with
\begin{subequations}
\begin{align}
	\kin{s}{(0)}(\n{}{}) & = 0,
	&
	\kin{s}{(1)}(\n{}{}) & = \frac{3\pi^{2}}{2} \qty(\frac{\n{}{}}{\pi})^{2/3},
	\\
	\e{\Ha}{(0)}(\n{}{}) & = \frac{8}{3} \qty(\frac{\n{}{}}{\pi})^{1/3},
	&
	\e{\Ha}{(1)}(\n{}{}) & = \frac{352}{105} \qty(\frac{\n{}{}}{\pi})^{1/3},
	\\
	\e{\ex}{(0)}(\n{}{}) & = - \frac{4}{3} \qty(\frac{\n{}{}}{\pi})^{1/3},
	&
	\e{\ex}{(1)}(\n{}{}) & = - \frac{176}{105} \qty(\frac{\n{}{}}{\pi})^{1/3}.	
\end{align}
\end{subequations}

In analogy with the conventional Dirac exchange functional, \cite{Dirac_1930} we write down the exchange functional of each individual state as
\begin{equation}
	\e{\ex}{(I)}(\n{}{}) = \Cx{(I)} \n{}{1/3},
\end{equation}
and we then obtain
\begin{align}
	 \Cx{(0)} & = - \frac{4}{3} \qty( \frac{1}{\pi} )^{1/3},
	 &
	 \Cx{(1)} & = - \frac{176}{105} \qty( \frac{1}{\pi} )^{1/3}.
\end{align}
We can now combine these two exchange functionals to create a weight-dependent exchange functional for a two-state ensemble
\begin{equation}
\label{eq:exw}
	\e{\ex}{\ew{}}(\n{}{}) 
	= (1-\ew{}) \e{\ex}{(0)}(\n{}{}) + \ew{} \e{\ex}{(1)}(\n{}{})
	= \Cx{\ew{}} \n{}{1/3}
\end{equation}
with
\begin{equation}
	\Cx{\ew{}} = (1-\ew{}) \Cx{(0)} + \ew{} \Cx{(1)}.
\end{equation}
Conveniently, the weight dependence of the exchange functional can be transferred to the $\Cx{}$ coefficient; an expected feature from a theoretical point of view, yet a nice property from a more practical aspect.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Weight-dependent correlation functional}
\label{sec:Ec}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Based on highly-accurate calculations, \cite{Loos_2009a,Loos_2009c,Loos_2010e} one can write down, for each state, an accurate analytical expression of the reduced correlation energy \cite{Loos_2013a, Loos_2014a} via the following Pad\'e approximant
\begin{equation}
\label{eq:ec}
	\e{\co}{(I)}(\n{}{}) = \frac{a_1^{(I)}}{1 + a_2^{(I)} \n{}{-1/6} + a_3^{(I)} \n{}{-1/3}},
\end{equation}
where $a_2^{(I)}$ and $a_3^{(I)}$ are state-specific fitting parameters, which are provided in Table \ref{tab:OG_func}.
The value of $a_1^{(I)}$ is obtained via the exact high-density expansion of the correlation energy. \cite{Loos_2011b}
Equation \eqref{eq:ec} is depicted in Fig.~\ref{fig:Ec} for each state alongside the data gathered in Table \ref{tab:Ref}.

Combining these, we build a two-state weight-dependent correlation functional:
\begin{equation}
\label{eq:ecw}
	\e{\co}{\ew{}}(\n{}{}) =  (1-\ew{}) \e{\co}{(0)}(\n{}{}) + \ew{} \e{\co}{(1)}(\n{}{}).
\end{equation}

%%% FIG 1 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{fig/fig1}
	\caption{
	Reduced (i.e., per electron) correlation energy $\e{\co}{(I)}$ [see Eq.~\eqref{eq:ec}] as a function of $R = 1/(\pi^2 \n{}{})^{1/3}$ for the ground state ($I=0$), and the first doubly-excited state ($I=1$) of the (spin-unpolarised) two-electron FUEG.
	The data gathered in Table \ref{tab:Ref} are also reported. 
	}
	\label{fig:Ec}
\end{figure}
%%% %%% %%%

%%% TABLE I %%%
\begin{table}
	\caption{
	\label{tab:Ref}
	$-\e{\co}{(I)}$ as a function of the radius of the glome $R = 1/(\pi^2 \n{}{})^{1/3}$ for the ground state ($I=0$), and the first doubly-excited state ($I=1$) of the (spin-unpolarised) two-electron FUEG.
	}
	\begin{ruledtabular}
	\begin{tabular}{lcc}
					&	\tabc{Ground state}	&	\tabc{Doubly-excited state}	\\
			$R$		&	\tabc{$I=0$}	&	\tabc{$I=1$}	\\
			\hline
			$0$		&	$0.023\,818$	&	$0.014\,463$	\\
			$0.1$	&	$0.023\,392$	&	$0.014\,497$	\\
			$0.2$	&	$0.022\,979$	&	$0.014\,523$	\\
			$0.5$	&	$0.021\,817$	&	$0.014\,561$	\\
			$1$		&	$0.020\,109$	&	$0.014\,512$	\\
			$2$		&	$0.017\,371$	&	$0.014\,142$	\\
			$5$		&	$0.012\,359$	&	$0.012\,334$	\\
			$10$	&	$0.008\,436$	&	$0.009\,716$	\\
			$20$	&	$0.005\,257$	&	$0.006\,744$	\\
			$50$	&	$0.002\,546$	&	$0.003\,584$	\\
			$100$	&	$0.001\,399$	&	$0.002\,059$	\\
			$150$	&	$0.000\,972$	&	$0.001\,458$	\\
	\end{tabular}
	\end{ruledtabular}
\end{table}

%%% TABLE 1 %%%
\begin{table}
	\caption{
	\label{tab:OG_func}
	Parameters of the correlation functionals for each individual state defined in Eq.~\eqref{eq:ec}.
	The values of $a_1$ are obtained to reproduce the exact high density correlation energy of each individual state, while $a_2$ and $a_3$ are fitted on the numerical values reported in Table \ref{tab:Ref}.}
	\begin{ruledtabular}
		\begin{tabular}{lll}
						&	\tabc{Ground state}	&	\tabc{Doubly-excited state}	\\
						&	\tabc{$I=0$}		&	\tabc{$I=1$}				\\
			\hline
				$a_1$	&	$-0.023\,818\,4$	&	$-0.014\,463\,3$	\\
				$a_2$	&	$+0.005\,409\,94$	&	$-0.050\,601\,9$	\\	
				$a_3$	&	$+0.083\,076\,6$	&	$+0.033\,141\,7$	\\
		\end{tabular}
	\end{ruledtabular}
\end{table}
%%% %%% %%% %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{LDA-centered functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Our intent is to incorporate into standard functionals (which are ``universal'' in the sense that they do not depend on the number of electrons) information about excited states that will be extracted from finite systems (whose properties may depend on the number of electrons).
Hence, we employ a simple embedding scheme where the two-electron FUEG (the impurity) is embedded in the IUEG (the bath).
The weight-dependence of the xc functional is then carried exclusively by the impurity [\ie, the functionals defined in Eqs.~\eqref{eq:exw} and \eqref{eq:ecw}], while the remaining effects are produced by the bath (\ie, the usual LDA xc functional).
 
Consistently with such a strategy, Eqs.~\eqref{eq:exw} and \eqref{eq:ecw} are ``centred'' on their corresponding jellium reference
\begin{equation}
\label{eq:becw}
	\be{\xc}{\ew{}}(\n{}{}) =  (1-\ew{}) \be{\xc}{(0)}(\n{}{}) + \ew{} \be{\xc}{(1)}(\n{}{})
\end{equation}
via the following shift:
\begin{equation}
	\be{\xc}{(I)}(\n{}{}) = \e{\xc}{(I)}(\n{}{}) + \e{\xc}{\LDA}(\n{}{}) - \e{\xc}{(0)}(\n{}{}).
\end{equation}
The LDA xc functional is similarly decomposed as
\begin{equation}
	\e{\xc}{\LDA}(\n{}{}) = \e{\ex}{\LDA}(\n{}{}) + \e{\co}{\LDA}(\n{}{}),
\end{equation}
where we consider here the Dirac exchange functional \cite{Dirac_1930} 
\begin{equation}
	\e{\ex}{\LDA}(\n{}{}) = \Cx{\LDA} \n{}{1/3},
\end{equation}
with 
\begin{equation}
	\Cx{\LDA} = -\frac{3}{4} \qty(\frac{3}{\pi})^{1/3},
\end{equation}
and the VWN5 correlation functional \cite{Vosko_1980}
\begin{equation}
	\e{\co}{\LDA}(\n{}{}) \equiv \e{\co}{\text{VWN5}}(\n{}{}).
\end{equation}
For the sake of clarity, the explicit expression of the VWN5 functional is not reported here but it can be found in Ref.~\onlinecite{Vosko_1980}.

Equation \eqref{eq:becw} can be recast
\begin{equation}
\label{eq:eLDA}
\begin{split}
	\be{\xc}{\ew{}}(\n{}{}) 
	& =  \e{\xc}{\LDA}(\n{}{}) + \ew{} \qty[\e{\xc}{(1)}(\n{}{}) - \e{\xc}{(0)}(\n{}{})]
	\\
	& =  \e{\xc}{\LDA}(\n{}{}) + \ew{} \pdv{\e{\xc}{\ew{}}(\n{}{})}{\ew{}},
\end{split}
\end{equation}
which nicely highlights the centrality of the LDA in the present eDFA. 
In particular, $\be{\xc}{(0)}(\n{}{}) = \e{\xc}{\LDA}(\n{}{})$.
Consequently, in the following, we name this weight-dependent xc functional ``eLDA'' as it is a natural extension of the LDA for ensembles.
Also, we note that, by construction,
\begin{equation}
\label{eq:dexcdw}
	\left. \pdv{\be{\xc}{\ew{}}[\n{}{}]}{\ew{}}\right|_{\n{}{} = \n{}{\ew{}}(\br)} 
	= \be{\xc}{(1)}(\n{}{\ew{}}(\br)) - \be{\xc}{(0)}(\n{}{\ew{}}(\br)).
\end{equation}

This embedding procedure can be theoretically justified by the generalised adiabatic connection formalism for ensembles (GACE) 
\begin{equation}
\label{eq:GACE}
	\E{\xc}{\bw}[\n{}{}] 
	= \E{\xc}{}[\n{}{}] 
	+ \sum_{I=1}^{\nEns-1} \int_0^{\ew{I}} \pdv{\E{\xc}{(0,\ldots,0,\xi,\ew{I},\ldots,\ew{\nEns-1})}[\n{}{}]}{\xi} d\xi,
\end{equation}
(where $\E{\xc}{}[\n{}{}]$ is the usual ground-state xc functional) originally derived by Franck and Fromager. \cite{Franck_2014} 
Within this in-principle-exact formalism, the (weight-dependent) xc energy of the ensemble is constructed from the (weight-independent) ground-state functional. 
In the case of a homogeneous system (or equivalently within the LDA), substituting Eq.~\eqref{eq:dexcdw} into \eqref{eq:GACE} yields, in the case of a bi-ensemble, Eq.~\eqref{eq:eLDA}.

%%% TABLE I %%%
\begin{table*}
\caption{
Total energies (in hartree) and excitation energies (in hartree) of \ce{H2} with $\RHH = 1.4$ bohr for various methods with the STO-3G minimal basis.
\label{tab:Energies}
}
\begin{ruledtabular}
\begin{tabular}{lddddddd}
	Method		&	\tabc{$\E{}{(0)}$}	&	\tabc{$\E{}{(1)}$}	&	\tabc{$\E{}{(1)} - \E{}{(0)}$}	
				&	\tabc{$\tE{}{\ew{} = 0}$}		&	\tabc{$\left. \pdv{\tE{}{\ew{}}}{\ew{}} \right|_{\ew{} = 0}$}		
				&	\tabc{$\tE{}{\ew{} = 1/2}$}		&	\tabc{$\left. \pdv{\tE{}{\ew{}}}{\ew{}} \right|_{\ew{} = 1/2}$}	\\
	\hline
	HF				&	-1.11671	&	0.460576	&	1.57729	&	-1.11671	&	2.49694			&	-0.0981563	&	1.57729	\\
	LDA				&	-1.12120	&	0.379745	&	1.50095	&	-1.12120	&	1.49536			&	-0.370725	&	1.50565	\\
	eLDA			&	-1.12120	&	0.175337	&	1.29654	&	-1.12120	&	1.31995			&	-0.462421	&	1.30839	\\
	CID				&	-1.13728	&	0.481138	&	1.61841	&					\\
	accurate\fnm[1]	&				&				&			&					\\
\end{tabular}
\end{ruledtabular}
\fnt[1]{FCI/cc-pV5Z excitation energies computed with QUANTUM PACKAGE. \cite{QP2}}
\end{table*}
%%% %%% %%% %%%

%%% FIG 1 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{fig/GSetDES_exact_HF_LDA_eLDA}
\caption{
Total energies (in hartree) of \ce{H2} as a function of $\RHH$ (in bohr) for various methods with the STO-3G minimal basis.
%\label{fig:Energies}
}
\end{figure}
%%% %%% %%% %%%

%%% FIG 2 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{fig/ExcitationEnergyExact_wHF_wLDA_weLDA_w=0etw=0.5}
\caption{
Excitation energies (in hartree) of \ce{H2} as a function of $\RHH$ (in bohr) for various methods with the STO-3G minimal basis.
%\label{fig:Energies}
}
\end{figure}
%%% %%% %%% %%%

%%% FIG 3 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{fig/EnsembleEnergy_wHF_wLDA_weLDA_wHFbarre_wLDAbarre_weLDAbarre_R=1.4}
\caption{
Ensemble energies (in hartree) of \ce{H2} with $\RHH = 1.4$ bohr as a function of the weight $\ew{}$ for various methods with the STO-3G minimal basis.
%\label{tab:Energies}
}
\end{figure}
%%% %%% %%% %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% COMPUTATIONAL DETAILS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:compdet}
The GOK-DFT calculations have been performed with the \texttt{QuAcK} software, freely available on \texttt{github}, where the present functional has been implemented. 
For all calculations, we use the aug-cc-pVXZ (X = D, T, and Q) Dunning's family of atomic basis sets.
Numerical quadratures are performed with the \texttt{numgrid} library using 194 angular points (Lebedev grid) and a radial precision of $10^{-6}$. \cite{Becke_1988,Lindh_2001}

%%%%%%%%%%%%%%%
%%% RESULTS %%%
%%%%%%%%%%%%%%%
\section{Results and discussion}
\label{sec:resdis}
\titou{Bigger basis set calculations coming soon.}
Here, we consider as testing ground the minimal-basis \ce{H2} molecule.
We select STO-3G as minimal basis, and study the behaviour of the total energy of \ce{H2} as a function of the internuclear distance $\RHH$ (in bohr).
This minimal-basis example is quite pedagogical as the molecular orbitals are fixed by symmetry. 
We have then access to the individual densities of the ground and doubly-excited states (which is not straightforward in practice).
Moreover, thanks to the spatial symmetry and the minimal basis, the individual densities extracted from the ensemble density are equal to the \textit{exact} individual densities.
In other words, there is no density-driven error and the only error that we are going to observe is the functional-driven error (and this is what we want to study).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Results}
\label{sec:res}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The bonding and antibonding orbitals of the \ce{H2} molecule are given by
\begin{subequations}
\begin{align}
	\MO{1}{}(\br{}) & = \qty[ \AO{A}(\br{}) + \AO{B}(\br{}) ]/\sqrt{2(1 + S_{AB})},
	\\
	\MO{2}{}(\br{}) & = \qty[ \AO{A}(\br{}) - \AO{B}(\br{}) ]/\sqrt{2(1 - S_{AB})},
\end{align}
\end{subequations}
where $\AO{A}$ and $\AO{B}$ are the two contracted Gaussian basis functions centred on each of the nucleus, and $S_{AB} = \braket{\AO{A}}{\AO{B}}$ is the overlap between these two basis functions.

The HF energies of the ground state and the doubly-excited states are
\begin{subequations}
\begin{align}
	\label{eq:HF0}
	\E{\HF}{(0)} & = 2 \eHc{1} + 2 \eJ{11} - \eK{11},
	\\
	\label{eq:HF1}
	\E{\HF}{(1)} & = 2 \eHc{2} + 2 \eJ{22} - \eK{22},
\end{align}
\end{subequations}
with 
\begin{subequations}
\begin{align}
	\eHc{p}  & =  \int \MO{p}{}(\br{}) \qty[-\frac{\nabla^2}{2}  + \vext(\br{})] \MO{p}{}(\br{})d\br{},
	\\
	\eJ{pq}  & = \iint \frac{\MO{p}{}(\br{})\MO{p}{}(\br{}) \MO{q}{}(\br{}')\MO{q}{}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}',
	\\
	\eK{pq}  & =  \iint \frac{\MO{p}{}(\br{})\MO{q}{}(\br{}) \MO{q}{}(\br{}')\MO{p}{}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}'.
\end{align}
\end{subequations}
Note that, in the HF case, there is no self-interaction error as $\eJ{pp} = \eK{pp}$.
We also define the HF excitation energy as 
\begin{equation}
	\Ex{\HF}{(1)} = \E{\HF}{(1)} - \E{\HF}{(0)}.
\end{equation}
%The HF orbital energies are
%\begin{subequations}
%\begin{align}
%	\eps{1}{\HF} & = \eHc{1} + 2\eJ{11} - \eK{11},
%	\\
%	\eps{2}{\HF} & = \eHc{2} + 2\eJ{12} - \eK{12}.
%\end{align}
%\end{subequations}

It is also instructive to consider the CID (configuration interaction with doubles) excitation energies computed in the same (minimal) basis set.
The CID energies of the ground state and doubly-excited states are provided by the eigenvalues of the following CID matrix:
\begin{equation}
	\bH_\CID = 
	\begin{pmatrix}
		\E{\HF}{(0)}	&	\eK{12}
		\\
		\eK{12}		&	\E{\HF}{(1)}
	\end{pmatrix}.
\end{equation}
These CID energies are explicitly given by
\begin{subequations}
\begin{align}
	\E{\CID}{(0)} & = \frac{\E{\HF}{(0)} + \E{\HF}{(1)}}{2} - \frac{1}{2} \sqrt{\qty(\E{\HF}{(1)} - \E{\HF}{(0)})^2 + 4 \eK{12}^2},
	\\
	\E{\CID}{(1)} & = \frac{\E{\HF}{(0)} + \E{\HF}{(1)}}{2} + \frac{1}{2} \sqrt{\qty(\E{\HF}{(1)} - \E{\HF}{(0)})^2 + 4 \eK{12}^2},
\end{align}
\end{subequations}
and the CID excitation energy reads
\begin{equation}
	\Ex{\CID}{(1)} = \sqrt{\qty(\E{\HF}{(1)} - \E{\HF}{(0)})^2 + 4 \eK{12}^2} \ge \Ex{\HF}{(1)}.
\end{equation}

At the (ground-state) LDA level (\ie, we only consider ground-state functionals), the individual energies reads
\begin{subequations}
\begin{align}
	\label{eq:LDA0}
	\E{\LDA}{(0)} & = 2 \eHc{1} + 2 \eJ{11} + \int \e{\xc}{\LDA}(\n{}{(0)}(\br{})) \n{}{(0)}(\br{}) d\br{},
	\\
	\label{eq:LDA1}
	\E{\LDA}{(1)} & = 2 \eHc{2} + 2 \eJ{22} + \int \e{\xc}{\LDA}(\n{}{(1)}(\br{})) \n{}{(1)}(\br{}) d\br{},
\end{align}
\end{subequations}
with
\begin{align}
	\n{}{(0)}(\br{}) & = 2 \MO{1}{2}(\br{}),
	& 
	\n{}{(1)}(\br{}) & = 2 \MO{2}{2}(\br{}),
\end{align}
Note that, contrary to the HF case, self-interaction is present in LDA. \cite{Perdew_1981}
%The KS orbital energies are given by
%\begin{subequations}
%\begin{align}
%	\eps{1}{\LDA} 
%	& = \eHc{1} + 2\eJ{11} 
%	+ \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}[\n{}{}]}{\n{}{}} \right|_{\n{}{} = \n{}{(0)}(\br{})} \n{}{(0)}(\br{}) d\br{},
%	\\
%	\eps{2}{\LDA} & = \eHc{2} + 2\eJ{12} 
%	+ \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}[\n{}{}]}{\n{}{}} \right|_{\n{}{} = \n{}{(0)}(\br{})} \n{}{(1)}(\br{}) d\br{}.
%\end{align}
%\end{subequations}

At the eLDA, we have
\begin{subequations}
\begin{align}
	\label{eq:eLDA0}
	\E{\eLDA}{(0)} & = 2 \eHc{1} + 2 \eJ{11} + \int \be{\xc}{(0)}(\n{}{(0)}(\br{})) \n{}{(0)}(\br{}) d\br{},
	\\
	\label{eq:eLDA1}
	\E{\eLDA}{(1)} & = 2 \eHc{2} + 2 \eJ{22} + \int \be{\xc}{(1)}(\n{}{(1)}(\br{})) \n{}{(1)}(\br{}) d\br{},
\end{align}
\end{subequations}
with $\be{\xc}{(0)}(\n{}{}) \equiv \e{\xc}{\LDA}(\n{}{})$ and $\be{\xc}{(1)}(\n{}{}) = \e{\xc}{\LDA}(\n{}{}) + \e{\xc}{(1)}(\n{}{}) - \e{\xc}{(0)}(\n{}{})$.

Interestingly here, there is a strong connection between the LDA and eLDA excitation energies:
\begin{equation}
\begin{split}
	\Ex{\eLDA}{(1)} 
	& = \Ex{\LDA}{(1)} + \int \qty[ \e{\xc}{(1)}(\n{}{(1)}(\br{})) - \e{\xc}{(0)}(\n{}{(1)}(\br{})) ] \n{}{(1)}(\br{}) d\br{}.
	\\
	& = \Ex{\LDA}{(1)} + \int \left. \pdv{\e{\xc}{\ew{}}(\n{}{})}{\ew{}} \right|_{\n{}{} = \n{}{(1)}(\br{})} \n{}{(1)}(\br{}) d\br{}.
\end{split}
\end{equation}

These equations can be combined to define three ensemble energies
\begin{subequations}
\begin{align}
	\label{eq:EwHF}
	\E{\HF}{\ew{}} & = (1-\ew{}) \E{\HF}{(0)} + \ew{} \E{\HF}{(1)},
	\\
	\label{eq:EwLDA}
	\E{\LDA}{\ew{}} & = (1-\ew{}) \E{\LDA}{(0)} + \ew{} \E{\LDA}{(1)},
	\\
	\label{eq:EweLDA}
	\E{\eLDA}{\ew{}} & = (1-\ew{}) \E{\eLDA}{(0)} + \ew{} \E{\eLDA}{(1)},
\end{align}
\end{subequations}
which are all, by construction, linear with respect to $\ew{}$.
Excitation energies can be easily extracted from these formulae via differentiation with respect to $\ew{}$.
Note that this is not how one would do in a ``practical'' ensemble calculation as one does not have (usually) access to the individual densities.
However, for pedagogical purposes and to study the magnitude of the ghost-interaction error, it is interesting to defined them.

Alternatively to Eqs.~\eqref{eq:EwHF}, \eqref{eq:EwLDA} and \eqref{eq:EweLDA}, one can obtain ensemble energies directly from the ensemble density 
\begin{equation}
	\n{}{\ew{}} = (1-\ew{}) \n{}{(0)} + \ew{} \n{}{(1)}.
\end{equation}
(This is what one would do in practice, \ie, by performing a KS ensemble calculation.)
We will label these energies as $\tE{}{\ew{}}$ to avoid confusion with the expressions reported in Eqs.~\eqref{eq:EwHF}, \eqref{eq:EwLDA} and \eqref{eq:EweLDA}.
For HF, we have
\begin{equation}
\label{eq:bEwHF}
\begin{split}
	\tE{\HF}{\ew{}} 
	& = \Ts{\ew{}}[\n{}{\ew{}}(\br{})] 
	+ \int \vext(\br{}) \n{}{\ew{}}(\br{}) d\br{}
	\\
	& + \frac{1}{2} \iint \frac{\n{}{\ew{}}(\br{})\n{}{\ew{}}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}'
	\\
	& = 2 (1-\ew{}) \eHc{1} + 2 \ew{} \eHc{2} + (1-\ew{})^2 (2\eJ{11}- \eK{11}) 
	\\
	& + \ew{}^2 (2\eJ{22}- \eK{22}) + 2 (1-\ew{})\ew{} (2 \eJ{12} - \eK{12}),
\end{split}
\end{equation}
which is clearly quadratic with respect to $\ew{}$ due to the ghost interaction error in the Hartree term.
In the case of the LDA, it reads
\begin{equation}
\label{eq:bEwLDA}
\begin{split}
	\tE{\LDA}{\ew{}} 
	& = \Ts{\ew{}}[\n{}{\ew{}}(\br{})]
	+ \int \vext(\br{}) \n{}{\ew{}}(\br{}) d\br{}
	\\
	& + \iint \frac{\n{}{\ew{}}(\br{})\n{}{\ew{}}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}' 
	+ \int \e{\xc}{\LDA}(\n{}{\ew{}}(\br{})) \n{}{\ew{}}(\br{}) d\br{}
	\\
	& = 2 (1-\ew{}) \eHc{1} + 2 \ew{} \eHc{2} 
	\\
	& + 2(1-\ew{})^2 \eJ{11} + 2\ew{}^2 \eJ{22} + 4 (1-\ew{})\ew{} \eJ{12}
	\\
	& + (1-\ew{}) \int \e{\xc}{\LDA}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{}
	\\
	& + \ew{} \int \e{\xc}{\LDA}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{},
\end{split}
\end{equation}
which is also clearly quadratic with respect to $\ew{}$ because the (weight-independent) LDA functional cannot compensate the ``quadraticity'' of the Hartree term.
For eLDA, the ensemble energy can be decomposed as
\begin{equation}
\label{eq:bEweLDA}
\begin{split}
	\tE{\eLDA}{\ew{}} 
	& = \Ts{\ew{}}[\n{}{\ew{}}(\br{})] 
	+ \int \vext(\br{}) \n{}{\ew{}}(\br{}) d\br{}
	\\ 
	& + \iint \frac{\n{}{\ew{}}(\br{})\n{}{\ew{}}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}' 
	+ \int \be{\xc}{\ew{}}(\n{}{\ew{}}(\br{})) \n{}{\ew{}}(\br{}) d\br{}
	\\
	& = 2 (1-\ew{}) \eHc{1} + 2 \ew{} \eHc{2} 
	\\
	& + 2(1-\ew{})^2 \eJ{11} + 2\ew{}^2 \eJ{22} + 4 (1-\ew{})\ew{} \eJ{12}
	\\
	& + \int \be{\xc}{\ew{}}(\n{}{\ew{}}(\br{})) \n{}{\ew{}}(\br{}) d\br{}
%	& + (1-\ew{})^2 \int \be{\xc}{(0)}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{}
%	\\
%	& + \ew{}^2 \int \be{\xc}{(1)}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{}
%	\\
%	& + (1-\ew{})\ew{} \int \be{\xc}{(0)}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{}
%	\\
%	& + \ew{}(1-\ew{}) \int \be{\xc}{(1)}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{}
%	\\
%	& = 2 (1-\ew{}) \eHc{1} + 2 \ew{} \eHc{2} 
%	\\
%	& 
%	+ (1-\ew{})^2 \qty[ 2\eJ{11} + \int \be{\xc}{(0)}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{} ]
%	\\
%	& + \ew{}^2 \qty[ 2\eJ{22} + \int \be{\xc}{(1)}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{} ]
%	\\
%	& + 2 (1-\ew{})\ew{} \qty[ 2\eJ{12} 
%	+ \frac{1}{2} \int \be{\xc}{(0)}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{} 
%	& + \frac{1}{2} \int \be{\xc}{(1)}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{} ],
\end{split}
\end{equation}
which \textit{could} be linear with respect to $\ew{}$ if the weight-dependent xc functional compensates exactly the quadratic terms in the Hartree term.
This would be, for example, the case with the exact xc functional.

Extracting excitation energies from Eqs.~\eqref{eq:bEwHF}, \eqref{eq:bEwLDA} and \eqref{eq:bEweLDA} is more tricky.
To do so, we will employ Eq.~\eqref{eq:dEdw}.
The two first terms are simply $\Eps{0}{\ew{}} = 2 \eps{1}{\ew{}}$, $\Eps{1}{\ew{}} = 2 \eps{2}{\ew{}}$, and the HF, LDA and eLDA weight-dependent orbital energies are, respectively,
\begin{subequations}
\begin{align}
	\eps{1}{\ew{},\HF} 
	& = \eHc{1} + (1-\ew{})(2\eJ{11} - \eK{11}) + \ew{}(2\eJ{12} - \eK{12}),
	\\
	\eps{2}{\ew{},\HF} 
	& = \eHc{2} +  (1-\ew{})(2\eJ{12} - \eK{12}) + \ew{}(2\eJ{22} - \eK{22}), 
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
\begin{split}
	\eps{1}{\ew{},\LDA} 
	& = \eHc{1} + 2(1-\ew{}) \eJ{11} + 2\ew{} \eJ{12}
	\\
	& + \frac{1}{2} \int \left. \pdv{\e{\xc}{\LDA}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{\ew{}}(\br{}) \n{}{(0)}(\br{}) d\br{},
	\\
	& + \frac{1}{2} \int \e{\xc}{\LDA}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{},
\end{split}
	\\
\begin{split}
	\eps{2}{\ew{},\LDA} 
	& = \eHc{2} +  2(1-\ew{}) \eJ{12} + 2 \ew{} \eJ{22}
	\\
	& + \frac{1}{2} \int \left. \pdv{\e{\xc}{\LDA}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{\ew{}}(\br{}) \n{}{(1)}(\br{}) d\br{},
	\\
	& + \frac{1}{2} \int \e{\xc}{\LDA}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{},
\end{split}
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
\begin{split}
	\eps{1}{\ew{},\eLDA} 
	& = \eHc{1} + 2(1-\ew{}) \eJ{11} + 2\ew{} \eJ{12} 
	\\
	& + \frac{1}{2} \int \left. \pdv{\be{\xc}{\ew{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{\ew{}}(\br{}) \n{}{(0)}(\br{}) d\br{},
	\\
	& + \frac{1}{2} \be{\xc}{\ew{}}(\n{}{\ew{}}(\br{})) \n{}{(0)}(\br{}) d\br{},
\end{split}
	\\
\begin{split}
	\eps{2}{\ew{},\eLDA} 
	& = \eHc{2} +  2(1-\ew{}) \eJ{12} + 2\ew{} \eJ{22}
	\\
	& + \frac{1}{2} \int \left. \pdv{\be{\xc}{\ew{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{\ew{}}(\br{}) \n{}{(1)}(\br{}) d\br{}.
	\\
	& + \frac{1}{2} \int \be{\xc}{\ew{}}(\n{}{\ew{}}(\br{})) \n{}{(1)}(\br{}) d\br{}.
\end{split}
\end{align}
\end{subequations}
The derivative discontinuity is modelled by the last term of the right-hand-side of Eq.~\eqref{eq:dEdw}.
Note that this contribution is only non-zero in the case of an explicitly weight-dependent functional [see Eq.~\eqref{eq:dexcdw}].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Discussion}
\label{sec:dis}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Numerical results are reported in Table \ref{tab:Energies}.


%%%%%%%%%%%%%%%%%%
%%% CONCLUSION %%%
%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:ccl}
As concluding remarks, we would like to say that what we have done, we think, is awesome.

%%%%%%%%%%%%%%%%%%%%%%%%
%%% ACKNOWLEDGEMENTS %%%
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
CM thanks the \textit{Universit\'e Paul Sabatier} (Toulouse, France) for a PhD scholarship.
%PFL thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481) for financial support.
This work has also been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''.} 
\end{acknowledgements}

%%%%%%%%%%%%%%%%%%%%
%%% BIBLIOGRAPHY %%%
%%%%%%%%%%%%%%%%%%%%
\bibliography{FarDFT}

\end{document}
