\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.
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}
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{}{}(\br{})$, the latter being a much simpler quantity than the many-electron wave function.
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,Gagliardi_2017}
Moreover, in exact TDDFT, the xc functional is time dependent.
The simplest and most widespread approximation in state-of-the-art electronic structure programs where TDDFT is implemented consists in neglecting memory effects. \cite{Casida}
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} can be employed to fix the various issues mentioned above.
The basic idea is to describe a finite ensemble of states (ground and excited) altogether, \ie, with the same set of orbitals.
Ensemble DFT (eDFT) was proposed at the end of the 80's by Gross, Oliveira and Kohn (GOK), \cite{Gross_1988a, Oliveira_1988, Gross_1988b} and is a generalization of Theophilou's variational principle for equiensembles. \cite{Theophilou_1979}
In GOK-DFT (\ie, eDFT for excited states), the (time-independent) xc functional depends explicitly on the weights assigned to the states that belong to the ensemble of interest.
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}
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}
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}
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}.
The determinants (or configuration state functions) $\Det{(K)}$ are all constructed from the same set of ensemble KS orbitals that is variationally optimized.
The trial ensemble density is simply the weighted sum of the individual densities, \ie,
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 here on the (exact) Hx part which is defined as
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).
where, 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 KS wavefunctions $\Det{(K)}=\Det{(K),\bw}$, \ie,
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
As Hartree and exchange energies cannot be separated in the one-dimensional systems considered in the rest of this work, we will substitute the Hartree--Fock (HF) density-matrix-functional interaction energy,
Note that this approximation, where the ensemble density matrix is optimized from a non-local exchange potential [rather than a local one, as expected from Eq.~\eqref{eq:var_ener_gokdft}] is applicable to real (three-dimensional) systems.
As readily seen from Eq.~\eqref{eq:eHF-dens_mat_func}, \textit{ghost interactions}~\cite{Gidopoulos_2002, Pastorczak_2014, Alam_2016, Alam_2017, Gould_2017}
\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
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, to construct eDFAs can be seen as more relevant. \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 that composed this FUEG (see below).
Here, we propose to construct a weight-dependent eDFA for the calculations of excited states in 1D systems by combining these FUEGs with the usual IUEG to construct a weigh-dependent LDA functional for ensembles (eLDA).
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.
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 derivative discontinuities.
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}
Here, we will consider the most simple ringium system featuring electronic correlation effects, \ie, the two-electron ringium model.
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.
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
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 is 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:
specifically designed for 1D systems in Ref.~\onlinecite{Loos_2013} as been used, where $F(a,b,c,x)$ is the Gauss hypergeometric function, \cite{NISTbook} and
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 particular, we investigate systems where $L$ ranges from $\pi/8$ to $8\pi$ and $2\le\nEl\le7$.
For the self-consistent calculations (such as HF, KS-DFT or KS-eDFT), the convergence threshold $\tau=\max{\abs{\bF{\bw}\bGam{\bw}\bS-\bS\bGam{\bw}\bF{\bw}}}$ been set to $10^{-5}$.
For KS-DFT, KS-eDFT and TDDFT calculations, a Gauss-Legendre quadrature is employed to compute the various integrals that cannot be performed in closed form.
In order to test the present eLDA functional we have performed various sets of calculations.
To get reference excitation energies for both the single and double excitations, we have performed full configuration interaction (FCI) calculations with the Knowles-Handy FCI program described in Ref.~\onlinecite{Knowles_1989}.
For the single excitations, we have also performed 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}
Concerning the KS-eDFT and eHF calculations, two sets of weight have been tested: the zero-weight limit where $\bw=(0,0)$ and the equi-ensemble (or state-averaged) limit where $\bw=(1/3,1/3)$.
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).
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\ew{1}\le1/3$ and $0\le\ew{2}\le\ew{1}$].
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}].
As one can see in Fig.~\ref{fig:EvsW}, the GOC-free ensemble energy becomes less and less linear as $L$ gets larger, while the GOC makes the ensemble energy almost perfectly linear.
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.
It is important to note that, even though the GIC removes the explicit quadratic terms from the ensemble energy, a weak non-linearity remains in the GIC ensemble energy due to the optimization of the ensemble KS orbitals in the presence of GIE.
However, this ``density-driven'' type of error is small (in our case at least) 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.
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}.
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}
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.
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.
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)$].
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).
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 state-averaged weights is able to accurately model single and double excitations, with a very significant improvement brought by the state-averaged KS-eLDA orbitals as compared to their zero-weight analogs.
As a rule of thumb, in the weak and intermediate correlation regimes, we see that KS-eLDA single excitations are of the same quality as the ones obtained in the linear response formalism (such as TDLDA), while double excitations only deviates from the FCI values by a few tenth of percent for these two box lengths.
Moreover, we note that, in the strong correlation regime (left graph of Fig.~\ref{fig:EvsN}), the single excitation energies obtained at the state-averaged KS-eLDA level remain in good agreement with FCI and are much more accurate than the TDLDA and TDA-TDLDA excitation energies which can deviate by up to $60\%$.
This also applies to double excitations, the discrepancy between FCI and KS-eLDA remaining of the order of a few percents in the strong correlation regime.
These observations nicely illustrate the robustness of the present state-averaged GOK-DFT scheme in any correlation regime for both single and double excitations.
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.
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 eHF levels [see Eqs.~\eqref{eq:EI-eLDA} and \eqref{eq:EI-eHF}, respectively] as a function of the box length $L$ in the case of 5-boxium.
The influence of the derivative discontinuity is clearly more important in the strong correlation regime.
Its contribution is also significantly larger in the case of the single excitation; the derivative discontinuity hardly influences the double excitation.
Importantly, one realizes that the magnitude of the derivative discontinuity is much smaller in the case of state-averaged 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 derivative discontinuity: for a given method, state-averaged orbitals partially remove the burden of modeling properly the derivative discontinuity.
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 state-averaged (\ie, $\ew{1}=\ew{2}=1/3$, blue and green lines) calculations are reported.
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 derivative discontinuity has a small impact on the double excitations with a slight tendency of worsening the excitation energies in the case of state-averaged weights, and a rather large influence on the single excitation energies obtained in the zero-weight limit, showing once again that the usage of state-averaged weights has the benefit of significantly reducing the magnitude of the derivative discontinuity.
In the present article, we have constructed a local, weight-dependent three-state DFA in the context of ensemble DFT.
The KS-eLDA scheme delivers accurate excitation energies for both single and double excitations, especially within its state-averaged version where the same weights are assigned to each state belonging to the ensemble.
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.
Using similar ideas, a three-dimensional version \cite{Loos_2009,Loos_2009c,Loos_2010,Loos_2010d,Loos_2017a} of the present eDFA is currently under development to model excited states in molecules and solids.
Similar to the present excited-state methodology for ensembles, one can easily design a local eDFA for the calculations of the ionization potential, electron affinity, and fundamental gap. \cite{Senjean_2018}
This can be done by constructing DFAs 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}.
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.
