eDFT_FUEG/Manuscript/eDFT.tex

\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{mathrsfs}

\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{\manu}[1]{\textcolor{blue}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashEF}[1]{\textcolor{blue}{\sout{#1}}}

%useful stuff
\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{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}

% operators
\newcommand{\hH}{\Hat{H}}
\newcommand{\hh}{\Hat{h}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\vne}{v_\text{ne}}
\newcommand{\hWee}{\Hat{W}_\text{ee}}
\newcommand{\WHF}{W_\text{HF}}

% functionals, potentials, densities, etc
\newcommand{\eps}{\epsilon}
\newcommand{\e}[2]{\eps_\text{#1}^{#2}}
\newcommand{\E}[2]{E_\text{#1}^{#2}}
\newcommand{\bE}[2]{\overline{E}_\text{#1}^{#2}}
\newcommand{\be}[2]{\overline{\eps}_\text{#1}^{#2}}
\newcommand{\bv}[2]{\overline{f}_\text{#1}^{#2}}
\newcommand{\n}[2]{n_{#1}^{#2}}
\newcommand{\DD}[2]{\Delta_\text{#1}^{#2}}
\newcommand{\LZ}[2]{\Xi_\text{#1}^{#2}}

% energies
\newcommand{\EHF}{E_\text{HF}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\Ecat}{E_\text{cat}}
\newcommand{\Eneu}{E_\text{neu}}
\newcommand{\Eani}{E_\text{ani}}
\newcommand{\EPT}{E_\text{PT2}}
\newcommand{\EFCI}{E_\text{FCI}}

% matrices/operator
\newcommand{\br}[1]{\boldsymbol{r}_{#1}}
\newcommand{\bx}[1]{\boldsymbol{x}_{#1}}
\newcommand{\bw}{{\boldsymbol{w}}}
\newcommand{\bG}{\boldsymbol{G}}
\newcommand{\bS}{\boldsymbol{S}}
\newcommand{\bGam}[1]{\boldsymbol{\Gamma}^{#1}}
\newcommand{\bgam}[1]{\boldsymbol{\gamma}^{#1}}
\newcommand{\opGam}[1]{\hat{\Gamma}^{#1}}
\newcommand{\bh}{\boldsymbol{h}}
\newcommand{\bF}[1]{\boldsymbol{F}^{#1}}
\newcommand{\Ex}[1]{\Omega^{#1}}


% elements
\newcommand{\ew}[1]{w_{#1}}
\newcommand{\eG}[1]{G_{#1}}
\newcommand{\eS}[1]{S_{#1}}
\newcommand{\eGam}[2]{\Gamma_{#1}^{#2}}
\newcommand{\hGam}[2]{\Hat{\Gamma}_{#1}^{#2}}
\newcommand{\eh}[2]{h_{#1}^{#2}}
\newcommand{\eF}[2]{F_{#1}^{#2}}
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\dbERI}[2]{(#1||#2)}

% Numbers
\newcommand{\nEl}{N}
\newcommand{\nBas}{K}

% AO and MO basis
\newcommand{\Det}[1]{\Phi^{#1}}
\newcommand{\MO}[2]{\phi_{#1}^{#2}}
\newcommand{\SO}[2]{\varphi_{#1}^{#2}}
\newcommand{\cMO}[2]{c_{#1}^{#2}}
\newcommand{\AO}[1]{\chi_{#1}}

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

\newcommand{\SI}{\textcolor{blue}{supplementary material}}

\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LCQ}{Laboratoire de Chimie Quantique, Institut de Chimie, CNRS, Universit\'e de Strasbourg, Strasbourg, France}

%%% added by Manu %%%
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\bmk}{\boldsymbol{\kappa}} % orbital rotation vector
\newcommand{\bmg}{\boldsymbol{\Gamma}} % orbital rotation vector
\newcommand{\bxi}{\boldsymbol{\xi}} 
\newcommand{\bfx}{{\bf{x}}}
\newcommand{\bfr}{{\bf{r}}}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
%%%% 

\begin{document}	

\title{Weight-dependent local density-functional approximations for ensembles}

\author{Pierre-Fran\c{c}ois Loos}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Emmanuel Fromager}
\email{fromagere@unistra.fr}
\affiliation{\LCQ}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
We report a first generation of local, weight-dependent correlation density-functional approximations (DFAs) that incorporate information about both ground and excited states in the context of density-functional theory for ensembles (eDFT). 
These density-functional approximations for ensembles (eDFAs) are specially designed for the computation of single and double excitations within eDFT, and can be seen as a natural extension of the ubiquitous local-density approximation for ensemble (eLDA).
The resulting eDFAs, based on both finite and infinite uniform electron gas models, automatically incorporate the infamous derivative discontinuity contributions to the excitation energies through their explicit ensemble weight dependence. 
Their accuracy is illustrated by computing single and double excitations in one-dimensional many-electron systems in the weak, intermediate and strong correlation regimes.
\titou{Although the present weight-dependent functional has been specifically designed for one-dimensional systems, the methodology proposed here is directly applicable to the construction of weight-dependent functionals for realistic three-dimensional systems, such as molecules and solids.}
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Over the last two decades, density-functional theory (DFT)
\cite{Hohenberg_1964,Kohn_1965} has become the method of choice for
modeling the electronic structure of large molecular systems and
materials. \cite{ParrBook}\manu{why this ref?}
The main reason is that, within DFT, the quantum contributions to the
electronic repulsion energy --- the so-called exchange-correlation (xc)
energy --- is rewritten as a functional of the electron density $n\equiv\n{}{}(\br{})$, the latter being a much simpler quantity than the many-electron wave function. 
The complexity of the many-body problem is then transferred to the xc
density functional. 
Despite its success, the standard Kohn-Sham (KS) formulation of DFT \cite{Kohn_1965} (KS-DFT) suffers, in practice, from various deficiencies. \cite{Woodcock_2002, Tozer_2003,Tozer_1999,Dreuw_2003,Sobolewski_2003,Dreuw_2004,Tozer_1998,Tozer_2000,Casida_1998,Casida_2000,Tapavicza_2008,Levine_2006}
The description of strongly multiconfigurational ground states (often
referred to as ``strong correlation problem'') still remains a
challenge. \cite{Gori-Giorgi_2010,Fromager_2015,Gagliardi_2017}
Another issue, which is partly connected to the previous one, is the description of electronically-excited states. 

The standard approach for modeling excited states in DFT is
linear-response time-dependent DFT (TDDFT). \cite{Runge_1984,Casida,Casida_2012}
In this case, the electronic spectrum relies on the (unperturbed) ground-state KS picture, which may break down when electron correlation is strong. 
Moreover, in exact TDDFT, the xc energy is in fact an xc action \cite{Vignale_2008} which is a
functional of the time-dependent density $n\equiv n(\br,t)$ and, as
such, has memory. Standard implementations of TDDFT rely on 
the adiabatic approximation where memory effects are neglected. In other
words, the xc functional is assumed to be local in time. \cite{Casida,Casida_2012}
As a result, double electronic excitations are completely absent from the TDDFT spectrum, thus reducing further the applicability of TDDFT. \cite{Maitra_2004,Cave_2004,Mazur_2009,Romaniello_2009a,Sangalli_2011,Mazur_2011,Huix-Rotllant_2011,Elliott_2011,Maitra_2012,Sundstrom_2014,Loos_2019}

When affordable (\ie, for relatively small molecules), time-independent
state-averaged wave function methods
\cite{Roos,Andersson_1990,Angeli_2001a,Angeli_2001b,Angeli_2002,Helgakerbook} can be employed to fix the various issues mentioned above. 
The basic idea is to describe a finite (canonical) ensemble of ground
and excited states altogether, \ie, with the same set of orbitals. 
Interestingly, a similar approach exists in DFT. Referred to as
Gross--Oliveira--Kohn (GOK) DFT\cite{Gross_1988a,Gross_1988b,Oliveira_1988}, it was proposed at the end of the 80's as a generalization
of Theophilou's DFT for equiensembles. \cite{Theophilou_1979}
In exact GOK-DFT, the ensemble xc energy not only a functional of the
density but also a  
function of the ensemble weights. Note that, unlike in conventional
Boltzmann ensembles~\cite{Pastorczak_2013}, the weights (each state in the ensemble
is assigned a given and fixed weight) are allowed to vary
independently in a GOK ensemble.   
The weight dependence of the xc functional plays a crucial role in the
calculation of excitation energies.
\cite{Gross_1988b,Yang_2014,Deur_2017,Deur_2019,Senjean_2018,Senjean_2020} 
It actually accounts for the infamous derivative discontinuity contribution to energy gaps. \cite{Levy_1995, Perdew_1983}
%\titou{Shall we further discuss the derivative discontinuity? Why is it important and where is it coming from?}

Despite its formal beauty and the fact that GOK-DFT can in principle tackle near-degenerate situations and multiple excitations, it has not been given much attention until recently. \cite{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,Senjean_2018,Smith_2016}
The main reason is simply the absence of density-functional approximations (DFAs) for ensembles in the literature. 
Recent works on this topic are still fundamental and exploratory, as they rely either on simple (but nontrivial) models like the Hubbard dimer \cite{Carrascal_2015,Deur_2017,Deur_2018,Deur_2019,Senjean_2015,Senjean_2016,Senjean_2018,Sagredo_2018} or on atoms for which highly accurate or exact-exchange-only calculations have been performed. \cite{Yang_2014,Yang_2017}
In both cases, the key problem, namely the design of weight-dependent DFAs for ensembles (eDFAs), remains open. 
A first step towards this goal is presented in the present manuscript with the ambition to turn, in the forthcoming future, GOK-DFT into a practical computational method for modeling excited states in molecules and extended systems.
The present eDFA is specially designed for the computation of single and double excitations within GOK-DFT, and can be seen as a natural extension of the ubiquitous local-density approximation (LDA) for ensemble. 
Consequently, we will refer to this eDFA as eLDA in the remaining of this paper.

In the following, the present methodology is illustrated on \emph{strict} one-dimensional (1D), spin-polarized electronic systems. \cite{Loos_2012, Loos_2013a, Loos_2014a, Loos_2014b}
In other words, the Coulomb interaction used in this work describes particles which are \emph{strictly} restricted to move within a 1D sub-space of three-dimensional space.
Despite their simplicity, 1D models are scrutinized as paradigms for quasi-1D materials \cite{Schulz_1993, Fogler_2005a} such as carbon nanotubes \cite{Bockrath_1999, Ishii_2003, Deshpande_2008} or nanowires. \cite{Meyer_2009, Deshpande_2010}
%Early models of 1D atoms using this interaction have been used to study the effects of external fields upon Rydberg atoms \cite{Burnett_1993, Mayle_2007} and the dynamics of surface-state electrons in liquid helium. \cite{Nieto_2000, Patil_2001}
This description of 1D systems also has interesting connections with the exotic chemistry of ultra-high magnetic fields (such as those in white dwarf stars), where the electronic cloud is dramatically compressed perpendicular to the magnetic field. \cite{Schmelcher_1990, Lange_2012, Schmelcher_2012}
In these extreme conditions, where magnetic effects compete with Coulombic forces, entirely new bonding paradigms emerge. \cite{Schmelcher_1990, Schmelcher_1997, Tellgren_2008, Tellgren_2009, Lange_2012, Schmelcher_2012, Boblest_2014, Stopkowicz_2015}

The paper is organized as follows. 
Section \ref{sec:eDFT} introduces the equations behind GOK-DFT.
In Sec.~\ref{sec:eDFA}, we detail the construction of the weight-dependent local correlation functional specially designed for the computation of single and double excitations within eDFT.
Computational details needed to reproduce the results of the present work are reported in Sec.~\ref{sec:comp_details}.
In Sec.~\ref{sec:res}, we illustrate the accuracy of the present eDFA by computing single and double excitations in one-dimensional many-electron systems in the weak, intermediate and strong correlation regimes.
Finally, we draw our conclusion in Sec.~\ref{sec:conclusion}.
Atomic units are used throughout.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:eDFT}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{GOK-DFT}\label{subsec:gokdft}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The GOK ensemble energy~\cite{Gross_1988a,Oliveira_1988,Gross_1988b} is defined as
\beq\label{eq:exact_GOK_ens_ener}
	\E{}{\bw}=\sum_{K \geq 0} \ew{K} \E{}{(K)},
\eeq
where the $K$th energy level $\E{}{(K)}$ [$K=0$ refers to the ground state] is the eigenvalue of the electronic Hamiltonian $\hH = \hh + \hWee$, where 
\beq
	\hh = \sum_{i=1}^\nEl \qty[ -\frac{1}{2} \nabla_{i}^2 + \vne(\br{i}) ]
\eeq
is the one-electron operator describing kinetic and nuclear attraction energies, and $\hat{W}_{\rm ee}$ is the electron repulsion operator. 
The (positive) ensemble weights $\ew{K}$ decrease with increasing index $K$. 
They are normalized, \ie,
\beq\label{eq:weight_norm_cond}
	\ew{0} = 1 - \sum_{K>0} \ew{K},
\eeq
so that only the weights $\bw \equiv \qty( \ew{1}, \ew{2}, \ldots, \ew{K}, \ldots )$ assigned to the excited states can vary independently.  
For simplicity we will assume in the following that the energies are not degenerate. 
Note that the theory can be extended to multiplets simply by assigning the same ensemble weight to all degenerate states~\cite{Gross_1988b}. 
In the KS formulation of GOK-DFT, \manu{which is simply referred to as
KS ensemble DFT (KS-eDFT) in the following}, the ensemble energy is determined variationally as follows~\cite{Gross_1988b}:
\beq\label{eq:var_ener_gokdft}
	\E{}{\bw}
 	= \min_{\opGam{\bw}}
	\qty{ 
	\Tr[\opGam{\bw} \hh] + \E{Hx}{\bw} \qty[\n{\opGam{\bw}}{}] + \E{c}{\bw} \qty[\n{\opGam{\bw}}{}]
	},
\eeq
where $\Tr$ denotes the trace and the trial ensemble density matrix operator reads
\beq
	\opGam{\bw}=\sum_{K \geq 0} \ew{K} \dyad*{\Det{(K)}}.
\eeq
The KS determinants [or configuration state functions~\cite{Gould_2017}]
$\Det{(K)}$ are all constructed from the same set of ensemble KS
orbitals that are variationally optimized.
The trial ensemble density in Eq.~(\ref{eq:var_ener_gokdft}) is simply
the weighted sum of the individual KS densities, \ie,
\beq\label{eq:KS_ens_density}
	\n{\opGam{\bw}}{}(\br{}) = \sum_{K\geq 0} \ew{K} \n{\Det{(K)}}{}(\br{}).
\eeq
As readily seen from Eq.~\eqref{eq:var_ener_gokdft}, both Hartree-exchange (Hx) and correlation (c) energies are described with density functionals that are \textit{weight dependent}. 
We focus in the following on the (exact) Hx part, which is defined as~\cite{Gould_2017}
\beq\label{eq:exact_ens_Hx}
	\E{Hx}{\bw}[\n{}{}]=\sum_{K \geq 0} \ew{K} \mel*{\Det{(K),\bw}[\n{}{}]}{\hWee}{\Det{(K),\bw}[\n{}{}]},
\eeq
where the KS wavefunctions fulfill the ensemble density constraint
\beq
	\sum_{K\geq 0} \ew{K} \n{\Det{(K),\bw}[\n{}{}]}{}(\br{}) = \n{}{}(\br{}).
\eeq 
The (approximate) description of the correlation part is discussed in
Sec.~\ref{sec:eDFA}.\\

In practice, the ensemble energy is not the most interesting quantity, and one is more concerned with excitation energies or individual energy levels (for geometry optimizations, for example). 
As pointed out recently in Ref.~\cite{Deur_2019}, the latter can be extracted
exactly from a single ensemble calculation as follows:
\beq\label{eq:indiv_ener_from_ens}
	\E{}{(I)} = \E{}{\bw} + \sum_{K>0} \qty(\delta_{IK} - \ew{K} )
\pdv{\E{}{\bw}}{\ew{K}},
\eeq
where, according to the normalization condition of Eq.~(\ref{eq:weight_norm_cond}), 
\beq
\pdv{\E{}{\bw}}{\ew{K}}= \E{}{(K)} -
\E{}{(0)}\equiv\Ex{(K)}
\eeq
corresponds to the $K$th excitation energy.
According to the {\it variational} ensemble energy expression of
Eq.~\eqref{eq:var_ener_gokdft}, the derivative with respect to $\ew{K}$
can be evaluated from the minimizing weight-dependent KS wavefunctions
$\Det{(K)} \equiv \Det{(K),\bw}$ as follows:
\beq\label{eq:deriv_Ew_wk}
\begin{split}
	\pdv{\E{}{\bw}}{\ew{K}} 
	& = \mel*{\Det{(K)}}{\hh}{\Det{(K)}}-\mel*{\Det{(0)}}{\hh}{\Det{(0)}}
	\\
	& + \Bigg\{\int \fdv{\E{Hx}{\bw}[\n{}{}]}{\n{}{}(\br{})} \qty[ \n{\Det{(K)}}{}(\br{}) - \n{\Det{(0)}}{}(\br{}) ] d\br{}
	+ \pdv{\E{Hx}{\bw} [\n{}{}]}{\ew{K}}
	\\
	& + \int \fdv{\E{c}{\bw}[n]}{\n{}{}(\br{})} \qty[ \n{\Det{(K)}}{}(\br{}) - \n{\Det{(0)}}{}(\br{}) ] d\br{} 
	+ \pdv{\E{c}{\bw}[n]}{\ew{K}}
	\Bigg\}_{\n{}{} = \n{\opGam{\bw}}{}}.
\end{split}
\eeq
The Hx contribution from Eq.~\eqref{eq:deriv_Ew_wk} can be recast as
\beq\label{eq:_deriv_wk_Hx}
	\left.
	\pdv{}{\xi_K} \qty(\E{Hx}{\bxi} [\n{}{\bxi,\bxi}]
	- \E{Hx}{\bw}[\n{}{\bw,\bxi}] )
	\right|_{\bxi=\bw},
\eeq
where $\bxi \equiv (\xi_1,\xi_2,\ldots,\xi_K,\ldots)$ and the
auxiliary double-weight ensemble density reads
\beq
	\n{}{\bw,\bxi}(\br{}) = \sum_{K\geq 0} \ew{K} \n{\Det{(K),\bxi}}{}(\br{}).
\eeq 
Since, for given ensemble weights $\bw$ and $\bxi$, the ensemble densities $\n{}{\bxi,\bxi}$ and $\n{}{\bw,\bxi}$ are generated from the \textit{same} KS potential (which is unique up to a constant), it comes
from the exact expression in Eq.~(\ref{eq:exact_ens_Hx}) that 
\beq
	\E{Hx}{\bxi}[\n{}{\bxi,\bxi}] = \sum_{K \geq 0} \xi_K \mel*{\Det{(K),\bxi}}{\hWee}{\Det{(K),\bxi}}
\eeq  
and
\beq
	\E{Hx}{\bw}[\n{}{\bw,\bxi}] = \sum_{K \geq 0} \ew{K} \mel*{\Det{(K),\bxi}}{\hWee}{\Det{(K),\bxi}}.
\eeq  
This yields, according to Eqs.~\eqref{eq:deriv_Ew_wk} and \eqref{eq:_deriv_wk_Hx}, the simplified expression
\beq\label{eq:deriv_Ew_wk_simplified}
\begin{split}
	\pdv{\E{}{\bw}}{\ew{K}} 
	& = \mel*{\Det{(K)}}{\hH}{\Det{(K)}} 
	- \mel*{\Det{(0)}}{\hH}{\Det{(0)}}
	\\
	& + \qty{
	\int \fdv{\E{c}{\bw}[\n{}{}]}{\n{}{}({\br{}})} 
	\qty[ \n{\Det{(K)}}{}(\br{}) - \n{\Det{(0)}}{}(\br{}) ]
	+
	\pdv{\E{c}{\bw} [\n{}{}]}{\ew{K}} 
	}_{\n{}{} = \n{\opGam{\bw}}{}} d\br{}.
\end{split}
\eeq
Since, according to Eqs.~(\ref{eq:var_ener_gokdft}) and (\ref{eq:exact_ens_Hx}), the ensemble energy can be evaluated as 
\beq
	\E{}{\bw} = \sum_{K \geq 0} \ew{K} \mel*{\Det{(K)}}{\hH}{\Det{(K)}} + \E{c}{\bw}[\n{\opGam{\bw}}{}],
\eeq 
with $\Det{(K)} = \Det{(K),\bw}$ [note that, when the minimum is reached in Eq.~\eqref{eq:var_ener_gokdft}, $\n{\opGam{\bw}}{} = \n{}{\bw,\bw}$], 
we finally recover from Eqs.~\eqref{eq:KS_ens_density} and
\eqref{eq:indiv_ener_from_ens} the {\it exact} expression of Ref.~\cite{Fromager_2020} for the $I$th energy level:
\beq\label{eq:exact_ener_level_dets}
\begin{split}
	\E{}{(I)} 
	& = \mel*{\Det{(I)}}{\hH}{\Det{(I)}} + \E{c}{{\bw}}[\n{\opGam{\bw}}{}]
	\\
	& + \int \fdv{\E{c}{\bw}[\n{\opGam{\bw}}{}]}{\n{}{}(\br{})}
	\qty[ \n{\Det{(I)}}{}(\br{}) - \n{\opGam{\bw}}{}(\br{}) ] d\br{}
	\\
	&+
	\sum_{K>0} \qty(\delta_{IK} - \ew{K} )
	\left.
	\pdv{\E{c}{\bw}[\n{}{}]}{\ew{K}}
	\right|_{\n{}{} = \n{\opGam{\bw}}{}}.
\end{split}
\eeq
Note that, when $\bw=0$, the ensemble correlation functional reduces to the
conventional (ground-state) correlation functional $E_{\rm c}[n]$. As a
result, the regular KS-DFT expression is recovered from
Eq.~(\ref{eq:exact_ener_level_dets}) for the ground-state energy:
\beq
\E{}{(0)}=\mel*{\Det{(0)}}{\hH}{\Det{(0)}} +
\E{c}{}[\n{\Det{(0)}}{}],
\eeq
or, equivalently,
\beq\label{eq:gs_ener_level_gs_lim}
\E{}{(0)}=\mel*{\Det{(0)}}{\hat{H}[\n{\Det{(0)}}{}]}{\Det{(0)}} 
,
\eeq
where the density-functional Hamiltonian reads
\beq\label{eq:dens_func_Hamilt}
\hat{H}[n]=\hH+
\sum^N_{i=1}\left(\fdv{\E{c}{}[n]}{\n{}{}(\br{i})}
+C_{\rm c}[n]
\right),
\eeq
and 
\beq\label{eq:corr_LZ_shift}
C_{\rm c}[n]=\dfrac{\E{c}{}[n]
        -\int
\fdv{\E{c}{}[n]}{\n{}{}(\br{})}n(\br{})d\br{}}{\int n(\br{})d\br{}}
\eeq
is the correlation component of
Levy--Zahariev's constant shift in potential~\cite{Levy_2014}. 
Similarly, the excited-state ($I>0$) energy level expressions
can be recast as follows: 
\beq\label{eq:excited_ener_level_gs_lim}
	\E{}{(I)} 
	 = \mel*{\Det{(I)}}{\hat{H}[\n{\Det{(0)}}{}]}{\Det{(I)}} 
      +
      \left. 
	\pdv{\E{c}{\bw}[\n{\Det{(0)}}{}]}{\ew{I}}
	\right|_{\bw=0}.
\eeq 
As readily seen from Eqs.~(\ref{eq:dens_func_Hamilt}) and
(\ref{eq:corr_LZ_shift}), introducing any constant shift $\delta
\E{c}{}[\n{\Det{(0)}}{}]/\delta n({\bf r})\rightarrow \delta
\E{c}{}[\n{\Det{(0)}}{}]/\delta n({\bf r})+C$ into the correlation
potential leaves the density-functional Hamiltonian $\hat{H}[n]$ (and
therefore the individual energy levels) unchanged. As a result, in
this context,
the correlation derivative discontinuities induced by the 
excitation process~\cite{Levy_1995} will be fully described by the ensemble
correlation derivatives [second term on the right-hand side of
Eq.~(\ref{eq:excited_ener_level_gs_lim})].    

%%%%%%%%%%%%%%%%
\subsection{One-electron reduced density matrix formulation}
%%%%%%%%%%%%%%%%
For implementation purposes, we will use in the rest of this work 
(one-electron reduced) density matrices 
as basic variables, rather than Slater determinants. If we expand the
ensemble KS (spin) orbitals [from which the latter determinants are constructed] in an atomic orbital (AO) basis,
\titou{\beq
	\SO{p}{}(\bx{}) = s(\omega) \sum_{\mu} \cMO{\mu p}{} \AO{\mu}(\br{}),
\eeq
where $\bx{}=(\omega,\br{})$ is a composite coordinate gathering spin and spatial degrees of freedom, and
\beq
	s(\omega)
	= 
	\begin{cases}
		\alpha(\omega),	&	\text{for spin-up electrons,}	\\
		\text{or}	\\
		\beta(\omega),	&	\text{for spin-down electrons,}	
	\end{cases}
\eeq
}then the density matrix of the   
determinant $\Det{(K)}$ can be expressed as follows in the AO basis:
\beq
	\bGam{(K)} \equiv \eGam{\mu\nu}{(K)} = \sum_{\SO{p}{} \in (K)} \cMO{\mu p}{} \cMO{\nu p}{},
\eeq
where the summation runs over the spinorbitals that are occupied in $\Det{(K)}$. 
\trashPFL{Note that, as the theory is applied later on to spin-polarized
systems, we drop spin indices in the density matrices, for convenience.}
\manu{Is the latter sentence ok with you?}
\titou{I don't think we need it anymore. What do you think?}
The electron density of the $K$th KS determinant can then be evaluated
as follows:
\beq
	\n{\bGam{(K)}}{}(\br{}) = \sum_{\mu\nu} \AO{\mu}(\br{}) \eGam{\mu\nu}{(K)} \AO{\nu}(\br{}),
\eeq 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manu's derivation %%%
\iffalse%%
\blue{
\beq
n_{\bmg^{(K)}}(\br{})&=&\sum_\sigma\left\langle\hat{\Psi}^\dagger(\br{}\sigma)\hat{\Psi}(\br{}\sigma)\right\rangle^{(K)}
\nonumber\\
&=&\sum_\sigma\sum_{pq}\varphi^\sigma_p(\br{})\varphi^\sigma_q(\br{})\left\langle\hat{a}_{p^\sigma,\sigma}^\dagger\hat{a}_{q^\sigma,\sigma}\right\rangle^{(K)}
\nonumber\\
&=&\sum_\sigma\sum_{\varphi^\sigma_p\in(K)}\left(\varphi^\sigma_p(\br{})\right)^2
\nonumber\\
&=&\sum_\sigma\sum_{\varphi^\sigma_p\in(K)}\sum_{\mu\nu}c^\sigma_{{\mu
p}}c^\sigma_{{\nu p}}\AO{\mu}(\br{})\AO{\nu}(\br{})
\nonumber\\
&=&\sum_{\mu\nu}\AO{\mu}(\br{})\AO{\nu}(\br{})\sum_\sigma\sum_{\varphi^\sigma_p\in(K)}c^\sigma_{{\mu
p}}c^\sigma_{{\nu p}}
\eeq
}
\fi%%%
%%%% end Manu
while the ensemble density matrix
and ensemble density read
\beq
	\bGam{\bw} 
	= \sum_{K\geq 0} \ew{K} \bGam{(K)}
	\equiv \eGam{\mu\nu}{\bw}
	= \sum_{K\geq 0} \ew{K} \eGam{\mu\nu}{(K)},
\eeq
and
\beq
	\n{\bGam{\bw}}{}(\br{}) = \sum_{\mu\nu} \AO{\mu}(\br{}) \eGam{\mu\nu}{\bw} \AO{\nu}(\br{}),
\eeq
respectively. 
The exact individual energy expression in Eq.~\eqref{eq:exact_ener_level_dets} can then be rewritten as
\beq\label{eq:exact_ind_ener_rdm}
\begin{split}
	\E{}{(I)} 
	& =\Tr[\bGam{(I)} \bh]
	+ \frac{1}{2} \Tr[\bGam{(I)} \bG \bGam{(I)}]
	+ \E{c}{{\bw}}[\n{\bGam{\bw}}{}]
	\\
	& + \int \fdv{\E{c}{\bw}[\n{\bGam{\bw}}{}]}{\n{}{}(\br{})} 
	\qty[ \n{\bGam{(I)}}{}(\br{}) - \n{\bGam{\bw}}{}(\br{}) ] d\br{} 
	\\
	& + \sum_{K>0} \qty(\delta_{IK} - \ew{K})
	\left. \pdv{\E{c}{\bw}[\n{}{}]}{\ew{K}}\right|_{\n{}{} = \n{\bGam{\bw}}{}}
	,
\end{split}
\eeq
where 
\beq
	\bh \equiv h_{\mu\nu} = \mel*{\AO{\mu}}{\hh}{\AO{\nu}}	
\eeq
denotes the one-electron integrals matrix.
The exact individual Hx energies are obtained from the following trace formula
\beq
	\Tr[\bGam{(K)} \bG \bGam{(L)}]
	= \sum_{\mu\nu\la\si} \eGam{\mu\nu}{(K)} \eG{\mu\nu\la\si} \eGam{\la\si}{(L)},
\eeq
where the antisymmetrized two-electron integrals read
\beq
	\bG 
	\equiv G_{\mu\nu\la\si} 
	= \dbERI{\mu\nu}{\la\si} 
	= \ERI{\mu\nu}{\la\si} - \ERI{\mu\si}{\la\nu},
\eeq
with
\beq
	\ERI{\mu\nu}{\la\si} = \iint \frac{\AO{\mu}(\br{1}) \AO{\nu}(\br{1}) \AO{\la}(\br{2}) \AO{\si}(\br{2})}{\abs{\br{1} - \br{2}}} d\br{1} d\br{2}.
\eeq
%Note that, in Sec.~\ref{sec:results}, the theory is applied to (1D) spin
%polarized systems in which $\eGam{\mu\nu}{(K)\beta}=0$ and
%$G_{\mu\nu\lambda\omega}^{\alpha\alpha}\equiv G_{\mu\nu\lambda\omega}=({\mu}{\nu}\vert{\lambda}{\omega})
%-(\mu\omega\vert\lambda\nu)$. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Hx energy ...
%%% Manu's derivation
\iffalse%%%%
\blue{
\beq
&&\dfrac{1}{2}\sum_{PQRS}\langle PQ\vert\vert
RS\rangle\eGam{PR}^{(K)}\eGam{QS}^{(L)}
\nonumber\\
&&
=\dfrac{1}{2}\sum_{\sigma,\tau}\sum_{p^{\sigma} q^{\tau}RS}
\nonumber\\
&&\Big(\langle p^\sigma\sigma q^\tau\tau\vert RS\rangle -\langle
p^\sigma\sigma q^\tau\tau
\vert SR\rangle
\Big)\Gamma^{(K)}_{p^\sigma\sigma,R}\Gamma^{(L)}_{q^\tau\tau, S}
\nonumber\\
&&
=\dfrac{1}{2}\sum_{\sigma,\tau}\sum_{p^{\sigma} q^{\tau}}
\nonumber\\
&&\Big(\sum_{r^\sigma s^\tau}\langle p^\sigma q^\tau\vert r^\sigma s^\tau\rangle
\Gamma^{(K)\sigma}_{p^\sigma r^\sigma}\Gamma^{(L)\tau}_{q^\tau s^\tau}
\nonumber\\
&& -\sum_{s^\sigma r^\tau}\langle
p^\sigma q^\tau
\vert s^\sigma r^\tau\rangle
\delta_{\sigma\tau}\Gamma^{(K)\sigma}_{p^\sigma
r^\sigma}\Gamma^{(L)\sigma}_{q^\sigma s^\sigma}\Big)
\nonumber\\
&&=\dfrac{1}{2}\sum_{\sigma,\tau}\sum_{p^{\sigma} q^{\tau}}
\nonumber\\
&&\left(\langle p^\sigma q^\tau\vert p^\sigma q^\tau\rangle
n_{p^\sigma}^{(K)\sigma}n_{q^\tau}^{(L)\tau}
-\delta_{\sigma\tau}\langle p^\sigma q^\sigma\vert q^\sigma p^\sigma \rangle
n_{p^\sigma}^{(K)\sigma}n_{q^\sigma}^{(L)\sigma}\right) 
\nonumber\\
&&=\dfrac{1}{2}\sum_{\mu\nu\lambda\omega}\sum_{\sigma,\tau}\Big(\langle{\mu}{\lambda}\vert{\nu}{\omega}\rangle
\Gamma_{\mu\nu}^{(K)\sigma}\Gamma_{\lambda\omega}^{(L)\tau}
-\delta_{\sigma\tau}\langle\mu\lambda\vert\omega\nu\rangle\Gamma_{\mu\nu}^{(K)\sigma}\Gamma_{\lambda\omega}^{(L)\sigma}
\Big)
\nonumber\\
&&=\dfrac{1}{2}\sum_{\mu\nu\lambda\omega}\sum_{\sigma,\tau}\Big(\langle{\mu}{\lambda}\vert{\nu}{\omega}\rangle
-\delta_{\sigma\tau}\langle\mu\lambda\vert\omega\nu\rangle
\Big)
\Gamma_{\mu\nu}^{(K)\sigma}\Gamma_{\lambda\omega}^{(L)\tau}
\nonumber\\
&&=\dfrac{1}{2}\sum_{\mu\nu\lambda\omega}\sum_{\sigma,\tau}\Big[({\mu}{\nu}\vert{\lambda}{\omega})
-\delta_{\sigma\tau}(\mu\omega\vert\lambda\nu)
\Big]
\Gamma_{\mu\nu}^{(K)\sigma}\Gamma_{\lambda\omega}^{(L)\tau}
\eeq
}
\fi%%%%%%%
%%%%
%%%%%%%%%%%%%%%%%%%%%
\iffalse%%%% Manu's derivation ...
\blue{
\beq
n^{\bw}({\br{}})&=&\sum_{K\geq 0}\sum_{\sigma=\alpha,\beta}{\tt
w}_Kn^{(K)}({\bfx})
\nonumber\\
&=&
\sum_{K\geq 0}\sum_{\sigma=\alpha,\beta}{\tt
w}_K\sum_{pq}\varphi_p({\bfx})\varphi_q({\bfx})\Gamma_{pq}^{(K)}
\nonumber\\
&=&
\sum_{\sigma=\alpha,\beta}
\sum_{K\geq 0}
{\tt
w}_K\sum_{p\in (K)}\varphi^2_p({\bfx})
\nonumber\\
&=&
\sum_{\sigma=\alpha,\beta}
\sum_{K\geq 0}
{\tt
w}_K
\sum_{\mu\nu}
\sum_{p\in (K)}c_{\mu p}c_{\nu p}\AO{\mu}({\bfx})\AO{\nu}({\bfx})
\nonumber\\
&=&\sum_{\sigma=\alpha,\beta}\sum_{\mu\nu}\AO{\mu}({\bfx})\AO{\nu}({\bfx}){\Gamma}^{\bw}_{\mu\nu}
\eeq
}
\fi%%%%%%%% end
%%%%%%%%%%%%%%%
%\subsection{Hybrid GOK-DFT}
%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%
\subsection{Approximations}\label{subsec:approx}
%%%%%%%%%%%%%%%


In the following, GOK-DFT will be applied
to one-dimension
spin-polarized systems where 
Hartree and exchange energies cannot be separated. 
For that reason, we will substitute the Hartree--Fock (HF) density-matrix-functional interaction energy,
\beq\label{eq:eHF-dens_mat_func}
	\WHF[\bGam{}] = \frac{1}{2} \Tr[\bGam{} \bG \bGam{}], 
\eeq
for the Hx density-functional energy in the variational energy
expression of Eq.~\eqref{eq:var_ener_gokdft}, thus leading to the
following approximation:  
\beq\label{eq:min_with_HF_ener_fun}
	\bGam{\bw} 
	\rightarrow \argmin_{\bgam{\bw}} 
	\qty{ 
	\Tr[\bgam{\bw} \bh ] + \WHF[ \bgam{\bw}] + \E{c}{\bw}[\n{\bgam{\bw}}{}]
	}.
\eeq
The minimizing ensemble density matrix in Eq.~(\ref{eq:min_with_HF_ener_fun}) fulfills the following
stationarity condition:
\beq\label{eq:commut_F_AO}
	\bF{\bw} \bGam{\bw} \bS = \bS \bGam{\bw} \bF{\bw},
\eeq
where $\bS \equiv \eS{\mu\nu} = \braket*{\AO{\mu}}{\AO{\nu}}$ is the
overlap matrix and the ensemble Fock-like matrix reads
\beq
	\bF{\bw} \equiv \eF{\mu\nu}{\bw} = \eh{\mu\nu}{\bw} +
\sum_{\la\si} \eG{\mu\nu\la\si} \eGam{\la\si}{\bw},
\eeq 
with
\beq
	\eh{\mu\nu}{\bw} 
	= \eh{\mu\nu}{} + \int \AO{\mu}(\br{}) \fdv{\E{c}{\bw}[\n{\bGam{\bw}}{}]}{\n{}{}(\br{})} \AO{\nu}(\br{}) d\br{}.
\eeq
%%%%%%%%%%%%%%%
\iffalse%%%%%%
% Manu's derivation %%%%
\color{blue}
I am teaching myself ...\\
Stationarity condition
\beq
&&0=\sum_{K\geq 0}w_K\sum_{t^\sigma}\Big(f_{p^\sigma\sigma,t^\sigma\sigma}\Gamma^{(K)\sigma}_{t^\sigma
q^\sigma}-\Gamma^{(K)\sigma}_{p^\sigma
t^\sigma}f_{t^\sigma\sigma,q^\sigma\sigma}\Big)
\nonumber\\
&&=\sum_{K\geq 0}w_K
\Big(f_{p^\sigma\sigma,q^\sigma\sigma}n^{(K)\sigma}_{q^\sigma}-n^{(K)\sigma}_{p^\sigma}f_{p^\sigma\sigma,q^\sigma\sigma}\Big)
\nonumber\\
&&
=\sum_{\mu\nu}\sum_{K\geq 0}w_KF_{\mu\nu}^\sigma c^\sigma_{\mu
p}c^\sigma_{\nu q}\left(n^{(K)\sigma}_{q^\sigma}-n^{(K)\sigma}_{p^\sigma}\right)
\eeq
thus leading to
\beq
&&0=\sum_{p^\sigma q^\sigma}c^\sigma_{\lambda
p}c^\sigma_{\omega q}\left(\sum_{\mu\nu}\sum_{K\geq 0}w_KF_{\mu\nu}^\sigma c^\sigma_{\mu
p}c^\sigma_{\nu q}\left(n^{(K)\sigma}_{q^\sigma}-n^{(K)\sigma}_{p^\sigma}\right)\right)
\nonumber\\
&&=\sum_{\mu\nu}\sum_{K\geq 0}w_K
F_{\mu\nu}^\sigma\left(\Gamma^{(K)\sigma}_{\nu\omega}\sum_{p^\sigma}c^\sigma_{\lambda
p}c^\sigma_{\mu
p}-\Gamma^{(K)\sigma}_{\mu\lambda}\sum_{q^\sigma}c^\sigma_{\omega q}c^\sigma_{\nu q}\right)
\nonumber\\
\eeq
If we denote $M^\sigma_{\lambda\mu}=\sum_{p^\sigma}c^\sigma_{\lambda
p}c^\sigma_{\mu
p}$ it comes 
\beq
S_{\mu\nu}=\sum_{\lambda\omega}S_{\mu\lambda}M^\sigma_{\lambda\omega}S_{\omega\nu}
\eeq
which simply means that
\beq
{\bm S}={\bm S}{\bm M}{\bm S}
\eeq
or, equivalently,
\beq
{\bm M}={\bm S}^{-1}.
\eeq
The stationarity condition simply reads
\beq
\sum_{\mu\nu}F_{\mu\nu}^\sigma\left(\Gamma^{\bw\sigma}_{\nu\omega}
\left[{\bm S}^{-1}\right]_{\lambda\mu}
-\Gamma^{\bw\sigma}_{\mu\lambda}\left[{\bm S}^{-1}\right]_{\omega\nu}\right)
=0
\eeq
thus leading to
\beq
{\bm S}^{-1}{{\bm F}^\sigma}{\bm \Gamma}^{\bw\sigma}={\bm \Gamma}^{\bw\sigma}{{\bm F}^\sigma}{\bm S}^{-1}
\eeq
or, equivalently,
\beq
{{\bm F}^\sigma}{\bm \Gamma}^{\bw\sigma}{\bm S}={\bm S}{\bm
\Gamma}^{\bw\sigma}{{\bm F}^\sigma}.
\eeq
%%%%%

Fock operator:\\
\beq
&&f_{p^\sigma\sigma,q^\sigma\sigma}-\langle\varphi_p^\sigma\vert\hat{h}\vert\varphi_q^\sigma\rangle
\nonumber\\
&&=\sum_{L\geq 0}w_L\sum_{\tau}\sum_{r^\tau s^\tau}
\nonumber\\
&&
\Big(\langle p^\sigma r^\tau\vert
q^\sigma s^\tau\rangle
-\delta_{\sigma\tau}\langle p^\sigma r^\sigma\vert
s^\sigma q^\sigma\rangle
\Big)
\Gamma^{(L)\tau}_{r^\tau
s^\tau}
\nonumber\\
&&
=\sum_{L\geq 0}w_L\sum_{\tau}\sum_{r^\tau}\Big(\langle p^\sigma r^\tau\vert
q^\sigma r^\tau\rangle
-\delta_{\sigma\tau}\langle p^\sigma r^\tau\vert
r^\tau q^\sigma\rangle
\Big)
n^{(L)\tau}_{r^\tau}
\nonumber\\
&&=\sum_{L\geq 0}w_L
\sum_{\lambda\omega}\sum_{\tau}\Big[\langle
p^\sigma\lambda\vert q^\sigma\omega\rangle
-\delta_{\sigma\tau}
\langle
p^\sigma\lambda\vert \omega q^\sigma\rangle\Big]
\Gamma^{(L)\tau}_{\lambda\omega}
\nonumber\\
&&=
\sum_{\lambda\omega}\sum_{\tau}\Big[\langle
p^\sigma\lambda\vert q^\sigma\omega\rangle
-\delta_{\sigma\tau}
\langle
p^\sigma\lambda\vert \omega q^\sigma\rangle\Big]
\Gamma^{\bw\tau}_{\lambda\omega}
\nonumber\\
&&=\sum_{\mu\nu\lambda\omega}\sum_{\tau}
\Big(\langle{\mu}{\lambda}\vert{\nu}{\omega}\rangle
-\delta_{\sigma\tau}\langle\mu\lambda\vert\omega\nu\rangle
\Big)\Gamma^{\bw\tau}_{\lambda\omega}c^\sigma_{\mu p}c^\sigma_{\nu q}
\nonumber\\
\eeq
or, equivalently,
\beq
f_{p^\sigma\sigma,q^\sigma\sigma}=\sum_{\mu\nu}F_{\mu\nu}^\sigma c^\sigma_{\mu p}c^\sigma_{\nu q}
\eeq
where
\beq
F_{\mu\nu}^\sigma=h_{\mu\nu}+\sum_{\lambda\omega}\sum_\tau
G_{\mu\nu\lambda\omega}^{\sigma\tau}\Gamma^{\bw\tau}_{\lambda\omega}
\eeq
and
\color{black}
\\
\fi%%%%%%%%%%%
%%%%% end Manu 
%%%%%%%%%%%%%%%%%%%%
Note that, within the approximation of Eq.~(\ref{eq:min_with_HF_ener_fun}), the ensemble density matrix is
optimized with a non-local exchange potential rather than a
density-functional local one, as expected from
Eq.~\eqref{eq:var_ener_gokdft}. This procedure is actually general, \ie,
applicable to not-necessarily spin polarized and real (higher-dimension) systems. 
As readily seen from Eq.~\eqref{eq:eHF-dens_mat_func}, inserting the
ensemble density matrix into the HF interaction energy functional
introduces unphysical \textit{ghost interaction} errors~\cite{Gidopoulos_2002, Pastorczak_2014, Alam_2016, Alam_2017, Gould_2017}
as well as {\it curvature}~\cite{Alam_2016,Alam_2017}:
\beq\label{eq:WHF}
\begin{split}
	\WHF[\bGam{\bw}] 
	& = \frac{1}{2} \sum_{K\geq 0} \ew{K}^2 \Tr[\bGam{(K)} \bG \bGam{(K)}]
	\\
	& + \sum_{L>K\geq 0} \ew{K} \ew{L}\Tr[\bGam{(K)} \bG \bGam{(L)}].
\end{split}
\eeq
The ensemble energy is of course expected to vary linearly with the ensemble
weights [see Eq.~(\ref{eq:exact_GOK_ens_ener})].
These errors are essentially removed when evaluating the individual energy
levels on the basis of Eq.~\eqref{eq:exact_ind_ener_rdm}.\\    

Turning to the density-functional ensemble correlation energy, the
following ensemble local density approximation (eLDA) will be employed:  
\beq\label{eq:eLDA_corr_fun}
	\E{c}{\bw}[\n{}{}]\approx \int \n{}{}(\br{}) \e{c}{\bw}(\n{}{}(\br{})) d\br{},
\eeq
where the correlation energy per particle $\e{c}{\bw}(\n{}{})$ is \textit{weight dependent}. 
As shown in Sec.~\ref{sec:eDFA}, the latter can be constructed, for
example, from a finite uniform electron gas model. 
\titou{Manu, I think we should clearly define here what the expression of the ensemble energy with and without GOC.
What do you think?}

Combining Eq.~\eqref{eq:exact_ind_ener_rdm} with Eq.~\eqref{eq:eLDA_corr_fun} leads to our final energy level expression within KS-eLDA: 
\beq\label{eq:EI-eLDA}
\begin{split}
	\E{{eLDA}}{(I)} 
	& = 
        \E{HF}{(I)}
	\\
%\Tr[\bGam{(I)} \bh] + \frac{1}{2} \Tr[\bGam{(I)} \bG \bGam{(I)}]
	& + \int \e{c}{\bw}(\n{\bGam{\bw}}{}(\br{})) \n{\bGam{(I)}}{}(\br{}) d\br{}
	\\
	&
	+ \int \n{\bGam{\bw}}{}(\br{}) \qty[ \n{\bGam{(I)}}{}(\br{}) - \n{\bGam{\bw}}{}(\br{}) ]
	\left. \pdv{\e{c}{{\bw}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{\bGam{\bw}}{}(\br{})} d\br{}
	\\
	& + \int \sum_{K>0} \qty(\delta_{IK} - \ew{K} ) \n{\bGam{\bw}}{}(\br{})
	\left. \pdv{\e{c}{\bw}(\n{}{})}{\ew{K}} \right|_{\n{}{}=\n{\bGam{\bw}}{}(\br{})} d\br{},
\end{split}
\eeq
where
\beq\label{eq:ind_HF-like_ener}
\E{HF}{(I)}=\Tr[\bGam{(I)} \bh] + \frac{1}{2} \Tr[\bGam{(I)} \bG \bGam{(I)}]
\eeq   
is the analog for ground and excited states (within an ensemble) of the HF energy.  
If, for analysis purposes, we Taylor expand the density-functional
correlation contributions
around the $I$th KS state density
$\n{\bGam{(I)}}{}(\br{})$, the sum of 
the second and third terms on the right-hand side
of Eq.~\eqref{eq:EI-eLDA} can be simplified as follows through first order in
$\n{\bGam{\bw}}{}(\br{})-\n{\bGam{(I)}}{}(\br{})$: 
\beq
\int \e{c}{\bw}(\n{\bGam{(I)}}{}(\br{})) \n{\bGam{(I)}}{}(\br{}) d\br{}
+\mathcal{O}\left([\n{\bGam{\bw}}{}-\n{\bGam{(I)}}{}]^2\right),
\eeq
and it can therefore be identified as 
an individual-density-functional correlation energy where the density-functional
correlation energy per particle is approximated by the ensemble one for
all the states within the ensemble. 
Let us finally stress that, to the best of our knowledge, eLDA is the first
density-functional approximation that incorporates ensemble weight
dependencies explicitly, thus allowing for the description of derivative
discontinuities [see Eq.~\eqref{eq:excited_ener_level_gs_lim} and the
comment that follows] {\it via} the last term on the right-hand side
of Eq.~\eqref{eq:EI-eLDA}.\\
\titou{In order to test the influence of the derivative discontinuity on the excitation energies, it is useful to perform ensemble HF (labeled as eHF) calculations in which the correlation effects are removed.
In this case, the individual energies are simply defined as
\beq\label{eq:EI-eHF}
	\E{eHF}{(I)} \approx \Tr[\bGam{(I)} \bh] + \frac{1}{2} \Tr[\bGam{(I)} \bG \bGam{(I)}].
\eeq
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Density-functional approximations for ensembles}
\label{sec:eDFA}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Paradigm}
\label{sec:paradigm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Most of the standard local and semi-local DFAs rely on the infinite uniform electron gas (IUEG) model (also known as jellium). \cite{ParrBook, Loos_2016}
One major drawback of the jellium paradigm, when it comes to develop eDFAs, is that the ground and excited states are not easily accessible like in a molecule. \cite{Gill_2012, Loos_2012, Loos_2014a, Loos_2014b, Agboola_2015, Loos_2017a}
Moreover, because the IUEG model is a metal, it is gapless, which means that both the fundamental and optical gaps are zero.
From this point of view, using finite UEGs (FUEGs), \cite{Loos_2011b,
Gill_2012} which have, like an atom, discrete energy levels and non-zero
gaps, can be seen as more relevant in this context. \cite{Loos_2014a, Loos_2014b, Loos_2017a}
However, an obvious drawback of using FUEGs is that the resulting eDFA
will inexorably depend on the number of electrons in the FUEG (see below).
Here, we propose to construct a weight-dependent eLDA for the
calculations of excited states in 1D systems by combining FUEGs with the
usual IUEG. 

As a FUEG, we consider the ringium model in which electrons move on a perfect ring (\ie, a circle) but interact \textit{through} the ring. \cite{Loos_2012, Loos_2013a, Loos_2014b}
The most appealing feature of ringium regarding the development of
functionals in the context of eDFT is the fact that both ground- and
excited-state densities are uniform, and therefore {\it equal}. 
As a result, the ensemble density will remain constant (and uniform) as the ensemble weights vary. 
This is a necessary condition for being able to model the ensemble
correlation derivatives with respect to the weights [last term
on the right-hand side of Eq.~(\ref{eq:exact_ener_level_dets})].
Moreover, it has been shown that, in the thermodynamic limit, the ringium model is equivalent to the ubiquitous IUEG paradigm. \cite{Loos_2013,Loos_2013a}
Let us stress that, in a FUEG like ringium, the interacting and
noninteracting densities match individually for all the states within the
ensemble
(these densities are all equal to the uniform density), which means that
so-called density-driven correlation
effects~\cite{Gould_2019,Gould_2019_insights,Senjean_2020,Fromager_2020} are absent from the model.
Here, we will consider the most simple ringium system featuring electronic correlation effects, \ie, the two-electron ringium model.

The present weight-dependent eDFA is specifically designed for the
calculation of excited-state energies within GOK-DFT.
In order to take into account both single and double excitations simultaneously, we consider a three-state ensemble including: 
(i) the ground state ($I=0$), (ii) the first singly-excited state ($I=1$), and (iii) the first doubly-excited state ($I=2$) of the (spin-polarized) two-electron ringium system.
All these states have the same (uniform) density $\n{}{} = 2/(2\pi R)$, where $R$ is the radius of the ring where the electrons are confined.
We refer the interested reader to Refs.~\onlinecite{Loos_2012, Loos_2013a, Loos_2014b} for more details about this paradigm.
Generalization to a larger number of states is straightforward and is left for future work.
To ensure the GOK variational principle, \cite{Gross_1988a} the
triensemble weights must fulfil the following conditions: \cite{Deur_2019}
\titou{$0 \le \ew{2} \le 1/3$ and $\ew{2} \le \ew{1} \le (1-\ew{2})/2$}. 
%The constraint in \titou{red} is wrong. If $\ew{2}=0$, you should be allowed
%to consider an equi-bi-ensemble
%for which $\ew{1}=1/2$. This possibility is excluded with your
%inequalities. The correct constraints are given in Ref.~\cite{Deur_2019}
%and are the ones you also mentioned, \ie, $0 \le \ew{2} \le 1/3$ and
%$\ew{2} \le \ew{1} \le (1-\ew{2})/2$.} 
%\manu{
%Just in case, starting from
%\beq
%\begin{split}
%0\leq \ew{2}\leq \ew{1}\leq (1-\ew{1}-\ew{2})
%\\
%\end{split}
%\eeq
%we obtain
%\beq
%0\leq \ew{2}\leq \ew{1}\leq (1-\ew{2})/2 
%\eeq
%which implies $\ew{2}\leq(1-\ew{2})/2$ or, equivalently, $\ew{2}\leq
%1/3$.
%}
%%% TABLE 1 %%%
\begin{table*}
	\caption{
	\label{tab:OG_func}
	Parameters of the weight-dependent correlation DFAs defined in Eq.~\eqref{eq:ec}.}
%	\begin{ruledtabular}
		\begin{tabular}{lcddd}
			\hline\hline
			State					&	$I$		&	\tabc{$a_1^{(I)}$}	&	\tabc{$a_2^{(I)}$}	&	\tabc{$a_3^{(I)}$}	\\
			\hline
			Ground state			&	$0$		&	-0.0137078			&	0.0538982			&	0.0751740			\\
			Singly-excited state	&	$1$		&	-0.0238184			&	0.00413142			&	0.0568648			\\
			Doubly-excited state	&	$2$		&	-0.00935749			&	-0.0261936			&	0.0336645			\\
			\hline\hline
		\end{tabular}
%	\end{ruledtabular}
\end{table*}
%%% %%% %%% %%%

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

Based on highly-accurate calculations (see {\SI} for additional details), one can write down, for each state, an accurate analytical expression of the reduced (\ie, per electron) correlation energy \cite{Loos_2013a, Loos_2014a} via the following Pad\'e approximant
\begin{equation}
\label{eq:ec}
	\e{c}{(I)}(\n{}{}) = \frac{a_1^{(I)}\,\n{}{}}{\n{}{} + a_2^{(I)} \sqrt{\n{}{}} + a_3^{(I)}},
\end{equation}
where the $a_k^{(I)}$'s are state-specific fitting parameters 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_2013a, Loos_2014a}
Equation \eqref{eq:ec} provides three state-specific correlation DFAs based on a two-electron system.
Combining these, one can build the following three-state weight-dependent correlation eDFA:
\begin{equation}
\label{eq:ecw}
	%\e{c}{\bw}(\n{}{})
 \tilde{\epsilon}_{\rm c}^\bw(n)=  (1-\ew{1}-\ew{2}) \e{c}{(0)}(\n{}{}) + \ew{1} \e{c}{(1)}(\n{}{}) + \ew{2} \e{c}{(2)}(\n{}{}).
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{LDA-centered functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
One of the main driving force behind the popularity of DFT is its ``universal'' nature, as xc density functionals can be applied to any electronic system.
Obviously, the two-electron-based eDFA defined in Eq.~\eqref{eq:ecw} does not have this feature as it does depend on the number of electrons constituting the FUEG.
However, one can partially cure this dependency by applying a simple embedding scheme in which the two-electron FUEG (the impurity) is embedded in the IUEG (the bath).
The weight-dependence of the correlation functional is then carried exclusively by the impurity [\ie, the functional defined in Eq.~\eqref{eq:ecw}], while the remaining correlation effects are provided by the bath (\ie, the usual LDA correlation functional).
Following this simple strategy, which can be further theoretically justified by the generalized adiabatic connection formalism for ensembles (GACE) originally derived by Franck and Fromager, \cite{Franck_2014} we propose to \emph{shift} the two-electron-based eDFA defined in Eq.~\eqref{eq:ecw} as follows:
\begin{equation}
\label{eq:becw}
	\tilde{\epsilon}_{\rm c}^\bw(n)\rightarrow{\e{c}{\bw}(\n{}{})} =  (1-\ew{1}-\ew{2}) \be{c}{(0)}(\n{}{}) + \ew{1} \be{c}{(1)}(\n{}{}) + \ew{2} \be{c}{(2)}(\n{}{}),
\end{equation}
where
\begin{equation}
	\be{c}{(I)}(\n{}{}) = \e{c}{(I)}(\n{}{}) + \e{c}{\text{LDA}}(\n{}{}) - \e{c}{(0)}(\n{}{}).
\end{equation}
In the following, we will use the LDA correlation functional that has been specifically designed for 1D systems in
Ref.~\onlinecite{Loos_2013}:
\begin{equation}
\label{eq:LDA}
	\e{c}{\text{LDA}}(\n{}{}) 
	= a_1^\text{LDA} F\qty[1,\frac{3}{2},a_3^\text{LDA}, \frac{a_1^\text{LDA}(1-a_3^\text{LDA})}{a_2^\text{LDA}} {\n{}{}}^{-1}],
\end{equation}
where $F(a,b,c,x)$ is the Gauss hypergeometric function, \cite{NISTbook} and
\begin{subequations}
\begin{align}
	a_1^\text{LDA} & = - \frac{\pi^2}{360},
	\\
	a_2^\text{LDA} & = \frac{3}{4} - \frac{\ln{2\pi}}{2},
	\\ 
	a_3^\text{LDA} & = 2.408779.
\end{align}
\end{subequations}
Note that the strategy described in Eq.~(\ref{eq:becw}) is general and
can be applied to real (higher-dimensional) systems. In order to make the
connection with the GACE formalism \cite{Franck_2014,Deur_2017} more explicit, one may  
recast Eq.~\eqref{eq:becw} as
\begin{equation}
\label{eq:eLDA}
\begin{split}
	{\e{c}{\bw}(\n{}{})}
	& =  \e{c}{\text{LDA}}(\n{}{}) 
	\\
	& + \ew{1} \qty[\e{c}{(1)}(\n{}{})-\e{c}{(0)}(\n{}{})] + \ew{2} \qty[\e{c}{(2)}(\n{}{})-\e{c}{(0)}(\n{}{})],
\end{split}
\end{equation}
or, equivalently,
\begin{equation}
\label{eq:eLDA_gace}
	{\e{c}{\bw}(\n{}{})}
	=  \e{c}{\text{LDA}}(\n{}{}) 
	+ \sum_{K>0}\int_0^{\ew{K}}
\qty[\e{c}{(K)}(\n{}{})-\e{c}{(0)}(\n{}{})]d\xi_K,
\end{equation}
where the $K$th correlation excitation energy (per electron) is integrated over the
ensemble weight $\xi_K$ at fixed (uniform) density $\n{}{}$. 
Equation \eqref{eq:eLDA_gace} nicely highlights the centrality of the LDA in the present eDFA. 
In particular, ${\e{c}{(0,0)}(\n{}{})} = \e{c}{\text{LDA}}(\n{}{})$.
Consequently, in the following, we name this correlation functional ``eLDA'' as it is a natural extension of the LDA for ensembles.
Finally, we note that, by construction,
\begin{equation}
	{\pdv{\e{c}{\bw}(\n{}{})}{\ew{J}} = \e{c}{(J)}(\n{}{}) - \e{c}{(0)}(\n{}{}).}
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:comp_details}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Having defined the eLDA functional in the previous section [see Eq.~\eqref{eq:eLDA}], we now turn to its validation.
Our testing playground for the validation of the eLDA functional is the ubiquitous ``electrons in a box'' model where $\nEl$ electrons are confined in a 1D box of length $L$, a family of systems that we call $\nEl$-boxium in the following.
In particular, we investigate systems where $L$ ranges from $\pi/8$ to $8\pi$ and $2 \le \nEl \le 7$. 
These inhomogeneous systems have non-trivial electronic structure properties which can be tuned by varying the box length.
For small $L$, the system is weakly correlated, while strong correlation effects dominate in the large-$L$ regime. \cite{Rogers_2017,Rogers_2016}

We use as basis functions the (orthonormal) orbitals of the one-electron system, \ie,
\begin{equation}
	\AO{\mu}(x) = 
	\begin{cases}
		\sqrt{2/L} \cos(\mu \pi x/L),	&	\mu \text{ is odd,}
		\\
		\sqrt{2/L} \sin(\mu \pi x/L),	&	\mu \text{ is even,}
	\end{cases}
\end{equation} 
with $ \mu = 1,\ldots,\nBas$ and $\nBas = 30$ for all calculations.
\manu{The convergence threshold $\tau = \max{ \abs{ \bF{\bw} \bGam{\bw}
\bS - \bS \bGam{\bw} \bF{\bw}}}$ [see Eq.~(\ref{eq:commut_F_AO})] is set
to $10^{-5}$. For comparison, regular HF and KS-DFT calculations 
are performed with the same threshold.
In order to compute the various density-functional
integrals that cannot be performed in closed form, 
a 51-point Gauss-Legendre quadrature is employed.}

In order to test the present eLDA functional we perform various sets of calculations.
To get reference excitation energies for both the single and double excitations, we compute full configuration interaction (FCI) energies with the Knowles-Handy FCI program described in Ref.~\onlinecite{Knowles_1989}.
For the single excitations, we also perform time-dependent LDA (TDLDA)
calculations [\ie, TDDFT with the LDA functional defined in
Eq.~\eqref{eq:LDA}], and the effect of the Tamm-Dancoff approximation
(TDA) has been also investigated. \cite{Dreuw_2005}\manu{Manu: has been
studied previously (if so why do you mention this?) or will be discussed
in the present work?}
Concerning the \manu{ensemble}
%KS-eDFT and eHF
calculations, two sets of weight are tested: the zero-weight
\manu{(ground-state)} limit where $\bw = (0,0)$ and the
equi\manu{-tri}-ensemble (or \manu{equal-weight} state-averaged) limit where $\bw = (1/3,1/3)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}
\label{sec:res}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% FIG 1 %%%
\begin{figure*}
	\includegraphics[width=\linewidth]{EvsW_n5}
	\caption{
	\label{fig:EvsW}
	Weight dependence of the KS-eLDA ensemble energy $\E{\titou{eLDA}}{(\ew{1},\ew{2})}$ with (dashed lines) and without (solid lines) ghost interaction correction (GIC) for 5-boxium (\ie, $\nEl = 5$) with a box of length $L = \pi/8$ (left), $L = \pi$ (center), and $L = 8\pi$ (right).
	}
\end{figure*}
%%% %%% %%%

First, we discuss the linearity of the ensemble energy.
To do so, we consider 5-boxium with box lengths of $L = \pi/8$, $L = \pi$, and $L = 8\pi$, which correspond (qualitatively at least) to the weak, intermediate, and strong correlation regimes, respectively.
The three-state ensemble energy $\E{}{(\ew{1},\ew{2})}$ is represented
in Fig.~\ref{fig:EvsW} as a function of both $\ew{1}$ and $\ew{2}$ while
fulfilling the restrictions on the ensemble weights to ensure the GOK
variational principle [\ie, $0 \le \manu{\ew{2}} \le 1/3$ and \manu{$\ew{2} \le \ew{1} \le (1-\ew{2})/2$}].
To illustrate the magnitude of the ghost interaction error (GIE), we report the KS-eLDA ensemble energy with and without ghost interaction correction (GIC) as explained above [see Eqs.~\eqref{eq:WHF} and \eqref{eq:EI-eLDA}].
\manu{Manu: Just to be sure. What you refer to as the GIC ensemble
energy is
\beq
\E{GIC-eLDA}{\bw}=\sum_{I\geq0}\ew{I}\E{{eLDA}}{(I)},
\eeq
right? (I will move this to the theory section later on). The ensemble
energy with GIE is the one computed in
Eq.~\eqref{eq:min_with_HF_ener_fun},
\beq
\E{HF-eLDA}{\bw}=\E{GIC-eLDA}{\bw}+\WHF[
\bGam{\bw}]-\sum_{I\geq0}\ew{I}\WHF[ \bGam{(I)}]
\eeq
\underline{Some suggestions for Fig. 1}: In order to "see" the curvature
it might be convenient to
plot $E^{(w_1,0)}-E^{(0,0)}$ and $E^{(1/3,w_2)}-E^{(1/3,0)}$ rather than $E^\bw$. Adding the exact curves
would be nice (we could see that the slope is also substantially
improved when introducing the GIC, at least in the strongly correlated
regime). Showing the linearly-interpolated energies also helps in
"seeing" the curvature.\\}
As one can see in Fig.~\ref{fig:EvsW}, \manu{without GIC}, the
\trashEF{GOC-free} ensemble energy becomes less and less linear as $L$
gets larger, while the \manu{GIC} makes the ensemble energy almost
perfectly linear. \manu{Manu: well, after all, it is not that stricking
for the bi-ensemble (black curves), as you point out in the following.
"Perfectly linear" is maybe too strong.}
In other words, the GIE increases as the correlation gets stronger.
\manu{Manu: discussing GIE while focusing exclusively on the linearity
is not completely relevant. The GIE is about interactions between two
different states. Individual interaction terms also have
quadratic-in-weight factors in front, which contribute to the curvature
of course. Our GIC removes not only the GIE (I guess we should see the
improvement by looking at the slope) but also the wrong factors in front
of individual interactions.}
Because the GIE can be easily computed via Eq.~\eqref{eq:WHF} even for
real, three-dimensional systems, this provides a cheap way of
quantifying strong correlation in a given electronic system.\manu{This
is a strong statement I am not sure about. The nature of the excitation
should also be invoked I guess (charge transfer or not, etc ...). If we look at the GIE:
\beq
\WHF[
\bGam{\bw}]-\sum_{I\geq0}\ew{I}\WHF[ \bGam{(I)}]
\eeq
For a bi-ensemble ($w_1=w$) it can be written as
\beq
\dfrac{1}{2}\left[(w^2-1)W_0+w(w-2)W_1\right]+w(1-w)W_{01} 
\eeq
If, for some reason, $W_0\approx W_1\approx W_{01}=W$, then the error
reduces to $-W/2$, which is weight-independent (it fits for example with
what you see in the weakly correlated regime). Such an assumption depends on the nature of the
excitation, not only on the correlation strength, right? Neverthless,
when looking at your curves, this assumption cannot be made when the
correlation is strong. It is not clear to me which integral ($W_{01}?$)
drives the all thing.\\} 
It is important to note that, even though the GIC removes the explicit
quadratic terms from the ensemble energy, a weak \manu{Manu: is it that weak
when correlation is strong? Look at the bi-ensemble case} non-linearity
remains in the GIC ensemble energy due to the optimization of the
ensemble KS orbitals in the presence of GIE [see Eq.~\eqref{eq:min_with_HF_ener_fun}].
However, this \manu{orbital-driven} error is small \manu{Manu: again, can we
really say "small" when looking at the strongly correlated case. It
seems to me that there is some residual curvature which is a signature
of the error in the orbitals} (in our case at
least) \trashEF{as the correlation part of the ensemble KS potential $\delta
\E{c}{\bw}[\n{}{}] /\delta \n{}{}(\br{})$ is relatively small compared
to the Hx contribution}.\manu{Manu: well, I guess that the problem arises
from the density matrices (or orbitals) that are used to compute
individual Coulomb-exchange energies (I would not expect the DFT
correlation part to have such an impact, as you say). The best way to check is to plot the
ensemble energy without the correlation functional.}\\
\\
\manu{Manu: another idea. As far as I can see we do
not show any individual energies (excitation energies are plotted in the
following). Plotting individual energies (to be compared with the FCI
ones) would immediately show if there is some curvature (in the ensemble
energy). The latter would
be induced by any deviation from the expected horizontal straight lines.}

%%% FIG 2 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{EvsL_5}
	\caption{
	\label{fig:EvsL}
	Excitation energies (multiplied by $L^2$) associated with the single excitation $\Ex{(1)}$ (bottom) and double excitation $\Ex{(2)}$ (top) of 5-boxium for various methods and box length $L$.
	Graphs for additional values of $\nEl$ can be found as {\SI}.
	}
\end{figure}
%%% %%% %%%

Figure \ref{fig:EvsL} reports the excitation energies (multiplied by $L^2$) for various methods and box sizes in the case of 5-boxium (\ie, $\nEl = 5$).
Similar graphs are obtained for the other $\nEl$ values and they can be found in the {\SI} alongside the numerical data associated with each method.
For small $L$, the single and double excitations can be labeled as ``pure''.
In other words, each excitation is dominated by a sole, well-defined reference Slater determinant.
However, when the box gets larger (\ie, $L$ increases), there is a strong mixing between the different excitation degrees.
In particular, the single and double excitations strongly mix, which makes their assignment as single or double excitations more discutable. \cite{Loos_2019}
This can be clearly evidenced by the weights of the different configurations in the FCI wave function.
Therefore, it is paramount to construct a two-weight \manu{correlation} functional
(\manu{\ie, a triensemble functional}, as we have done here) which
allows the mixing of \trashEF{single and double} \manu{singly- and doubly-excited} configurations.
Using a single-weight (\ie, a biensemble) functional where only the ground state and the lowest singly-excited states are taken into account, one would observe a neat deterioration of the excitation energies (as compared to FCI) when the box gets larger.
\titou{Shall we add results for $\ew{2} = 0$ to illustrate
this?}\manu{Well, neglecting the second excited state is not the same as
considering the $w_2=0$ limit. I thought you were referring to an
approximation where the triensemble calculation is performed with
the biensemble functional. This is not the same as taking $w_2=0$
because, in this limit, you may still have a derivative discontinuity
correction. The latter is absent if you truly neglect the second excited
state in your ensemble functional. This should be clarified.}\\
\manu{Are the results in the supp mat? We could just add "[not
shown]" if not. This is fine as long as you checked that, indeed, the
results deteriorate ;-)}
\manu{Should we add that, in the bi-ensemble case, the ensemble
correlation derivative $\partial \epsilon^\bw_{\rm c}(n)/\partial w_2$
is neglected (if this is really what you mean (?)). I guess that this is the reason why
the second excitation energy would not be well described (?)}

As shown in Fig.~\ref{fig:EvsL}, all methods provide accurate estimates of the excitation energies in the weak correlation regime (\ie, small $L$).
When the box gets larger, they start to deviate.
For the single excitation, TDLDA is extremely accurate up to $L = 2\pi$, but yields more significant errors at larger $L$ by underestimating the excitation energies.
TDA-TDLDA slightly corrects this trend thanks to error compensation. 
Concerning the eLDA functional, our results clearly evidence that the equiweight [\ie, $\bw = (1/3,1/3)$] excitation energies are much more accurate than the ones obtained in the zero-weight limit [\ie, $\bw = (0,0)$].
This is especially true for the single excitation\manu{Manu: in the
light of your comments about the mixed singly-excited/doubly-excited
character of the first and second excited states when correlation is
strong, I would refer to the
"first excitation" rather than the "single excitation" (to be corrected 
everywhere in the discussion if adopted)} which is significantly
improved by using state-averaged weights\manu{Manu: you mean equal-weight?
State-averaged does not mean equal-weight, don't you think? In the state-averaged CASSCF
you do not have to use equal weights, even though most people do}.
The effect on the \trashEF{double} \manu{second?} excitation is less pronounced.
Overall, one clearly sees that, with \trashEF{state-averaged}
\manu{equal} weights, KS-eLDA yields accurate excitation energies for both single and double excitations.
This conclusion is verified for smaller and larger numbers of electrons
(see {\SI}).\\
\manu{Manu: now comes the question that is, I believe, central in this
work. How important are the
ensemble correlation derivatives $\partial \epsilon^\bw_{\rm
c}(n)/\partial w_I$  that, unlike any functional
in the literature, the eLDA functional contains. We have to discuss this
point... I now see, after reading what follows that this question is
addressed later on. We should say something here and then refer to the
end of the section, or something like that ...}


%%% FIG 3 %%%
\begin{figure*}
	\includegraphics[width=\linewidth]{EvsN}
	\caption{
	\label{fig:EvsN}
	Error with respect to FCI in single and double excitation energies for $\nEl$-boxium for various methods and electron numbers $\nEl$ at $L=\pi/8$ (left), $L=\pi$ (center), and $L=8\pi$ (right).
	}
\end{figure*}
%%% %%% %%%

For the same set of methods, Fig.~\ref{fig:EvsN} reports the error (in \%) in excitation energies (as compared to FCI) as a function of $\nEl$ for three values of $L$ ($\pi/8$, $\pi$, and $8\pi$).
We draw similar conclusions as above: irrespectively of the number of
electrons, the eLDA functional with \trashEF{state-averaged} equal
weights is able to accurately model single and double excitations, with
a very significant improvement brought by the \trashEF{state-averaged}
\manu{equiensemble} KS-eLDA orbitals as compared to their zero-weight
\manu{(\ie, conventional ground-state)} analogs.
\manu{As a rule of thumb, in the weak and intermediate correlation regimes, we
see that the \trashEF{single} \manu{first
excitation} obtained from \manu{equiensemble} KS-eLDA is of
the same quality as the one obtained in the linear response formalism
(such as TDLDA). On the other hand, the \trashEF{double} second
excitation energy only deviates
from the FCI value by a few tenth of percent} \trashEF{for these two box
lengths}.
Moreover, we note that, in the strong correlation regime (left graph of
Fig.~\ref{fig:EvsN}), the \trashEF{single} \manu{first} excitation
energy obtained at the equiensemble KS-eLDA level remains in good
agreement with FCI and is much more accurate than the TDLDA and TDA-TDLDA excitation energies which can deviate by up to $60 \%$.
This also applies to \trashEF{double} \manu{the second} excitation
\manu{(which has a strong doubly-excited character)}, the discrepancy
between FCI and \manu{equiensemble} KS-eLDA remaining of the order of a few percents in the strong correlation regime.
These observations nicely illustrate the robustness of the
\trashEF{present state-averaged} GOK-DFT scheme in any correlation regime for both single and double excitations.
This is definitely a very pleasing outcome, which additionally shows
that, even though we have designed the eLDA functional based on a
two-electron model system, the present methodology is applicable to any
1D electronic system, \manu{\ie, a system that has more than two
electrons}. 

%%% FIG 4 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{EvsL_5_HF}
	\caption{
	\label{fig:EvsLHF}
	Error with respect to FCI (in \%) associated with the single excitation $\Ex{(1)}$ (bottom) and double excitation $\Ex{(2)}$ (top) as a function of the box length $L$ for 5-boxium at the KS-eLDA (solid lines) and eHF (dashed lines) levels. 
	Zero-weight (\ie, $\ew{1} = \ew{2} = 0$, red lines) and state-averaged (\ie, $\ew{1} = \ew{2} = 1/3$, blue lines) calculations are reported. 
	}
\end{figure}
%%% %%% %%%

\titou{T2: there is a micmac with the derivative discontinuity as it is
only defined at zero weight. We should clean up this.}\manu{I will!}
It is also interesting to investigate the influence of the derivative discontinuity on both the single and double excitations.
To do so, we have reported in Fig.~\ref{fig:EvsLHF} the error percentage
(with respect to FCI) on the excitation energies obtained at the KS-eLDA
and HF\manu{-like} levels [see Eqs.~\eqref{eq:EI-eLDA} and
\eqref{eq:ind_HF-like_ener}, respectively] as a function of the box
length $L$ in the case of 5-boxium.\\
\manu{Manu: there is something I do not understand. If you want to
evaluate the importance of the ensemble correlation derivatives you
should only remove the following contribution from the $K$th KS-eLDA
excitation energy:  
\beq\label{eq:DD_term_to_compute}
\int  \n{\bGam{\bw}}{}(\br{})
      \left. \pdv{\e{c}{\bw}(\n{}{})}{\ew{K}} \right|_{\n{}{}=\n{\bGam{\bw}}{}(\br{})} d\br{}
\eeq
%rather than $E^{(I)}_{\rm HF}$
}
The influence of the ensemble correlation derivative is clearly more important in the strong correlation regime.
Its contribution is also significantly larger in the case of the single
excitation; the ensemble correlation derivative hardly influences the double excitation.
Importantly, one realizes that the magnitude of the ensemble correlation
derivative is much smaller in the case of equal-weight calculations (as compared to the zero-weight calculations).
This could explain why equiensemble calculations are clearly more
accurate as it reduces the influence of the ensemble correlation derivative:
for a given method, equiensemble orbitals partially remove the burden
of modeling properly the ensemble correlation derivative.\manu{Manu: well, we
would need the exact derivative value to draw such a conclusion. We can
only speculate. Let us first see how important the contribution in
Eq.~\eqref{eq:DD_term_to_compute} is. What follows should also be
updated in the light of the new results.}

%%% FIG 5 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{EvsN_HF}
	\caption{
	\label{fig:EvsN_HF}
	Error with respect to FCI in single and double excitation energies for $\nEl$-boxium (with a box length of $L=8\pi$) as a function of the number of electrons $\nEl$ at the KS-eLDA (solid lines) and eHF (dashed lines) levels. 
	Zero-weight (\ie, $\ew{1} = \ew{2} = 0$, black and red lines) and
equal-weight (\ie, $\ew{1} = \ew{2} = 1/3$, blue and green lines) calculations are reported. 
	}
\end{figure}
%%% %%% %%%

Finally, in Fig.~\ref{fig:EvsN_HF}, we report the same quantities as a function of the electron number for a box of length $8\pi$ (\ie, in the strong correlation regime). 
The difference between the eHF and KS-eLDA excitation energies
undoubtedly show that, even in the strong correlation regime, the
ensemble correlation derivative has a small impact on the double
excitations with a slight tendency of worsening the excitation energies
in the case of equal weights, and a rather large influence on the single
excitation energies obtained in the zero-weight limit, showing once
again that the usage of equal weights has the benefit of significantly reducing the magnitude of the ensemble correlation derivative.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Concluding remarks}
\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A local and ensemble-weight-dependent correlation density-functional approximation
(eLDA) has been constructed in the context of GOK-DFT for spin-polarized
triensembles in
1D. The approach is actually general and can be extended to real
(three-dimensional)
systems~\cite{Loos_2009,Loos_2009c,Loos_2010,Loos_2010d,Loos_2017a} 
and larger ensembles in order to  
model excited states in molecules and solids. Work is currently in
progress in this direction.

Unlike any standard functional, eLDA incorporates derivative
discontinuities through its weight dependence. The latter originates
from the finite uniform electron gas eLDA is
(partially) based on. The KS-eLDA scheme, where exact exchange is
combined with eLDA, delivers accurate excitation energies for both
single and double excitations, especially when an equiensemble is used.
In the latter case, the same weights are assigned to each state belonging to the ensemble.
{\it We have observed that, although the derivative discontinuity has a
non-negligible effect on the excitation energies (especially for the
single excitations), its magnitude can be significantly reduced by
performing state-averaged calculations instead of zero-weight
calculations.}\manu{to be updated ...}

Let us finally stress that the present methodology can be extended
straightforwardly to other types of ensembles like, for example, the
$N$-centered ones\cite{Senjean_2018,Senjean_2020}, thus allowing for the design an LDA-type functional for the
calculation of ionization potentials, electron affinities, and
fundamental gaps. 
Like in the present
eLDA, such a functional would incorporate the infamous derivative
discontinuity contribution to the gap through its explicit weight
dependence. We hope to report on this in the near future.


\trashEF{This can be done by constructing a functional for the one- and
three-electron ground-state systems, and combining them with the
two-electron DFA in complete analogy with Eqs.~\eqref{eq:ec} and
\eqref{eq:ecw}.}\manu{I find the sentence too technical for a
conclusion.} 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supplementary material}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
See {\SI} for the additional details about the construction of the functionals, raw data and additional graphs.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
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''.} 
EF thanks the \textit{Agence Nationale de la Recherche} (MCFUNEX project, Grant No.~ANR-14-CE06-0014-01) for funding.
\end{acknowledgements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




\bibliography{eDFT}

\end{document}