# Exact ensemble density functional theory for excited states in a model system: investigating the weight dependence of the correlation energy

###### Abstract

Ensemble density functional theory (eDFT) is an exact time-independent alternative to time-dependent DFT (TD-DFT) for the calculation of excitation energies. Despite its formal simplicity and advantages in contrast to TD-DFT (multiple excitations, for example, can be easily taken into account in an ensemble), eDFT is not standard which is essentially due to the lack of reliable approximate exchange-correlation () functionals for ensembles. Following Smith et al. [Phys. Rev. B 93, 245131 (2016)], we propose in this work to construct an exact eDFT for the nontrivial asymmetric Hubbard dimer, thus providing more insight into the weight dependence of the ensemble energy in various correlation regimes. For that purpose, an exact analytical expression for the weight-dependent ensemble exchange energy has been derived. The complementary exact ensemble correlation energy has been computed by means of Legendre–Fenchel transforms. Interesting features like discontinuities in the ensemble potential in the strongly correlated limit have been rationalized by means of a generalized adiabatic connection formalism. Finally, functional-driven errors induced by ground-state density-functional approximations have been studied. In the strictly symmetric case or in the weakly correlated regime, combining ensemble exact exchange with ground-state correlation functionals gives relatively accurate ensemble energies. However, when approaching the equiensemble in the strongly correlated regime, this approximation leads to highly curved ensemble energies with negative slope which is unphysical. Using both ground-state exchange and correlation functionals gives much better results in that case. In fact, exact ensemble energies are almost recovered in some density domains. The analysis of density-driven errors is left for future work.

^{†}

^{†}preprint:

## I Introduction

Despite its success, time-dependent density functional theory
(TD-DFT) Runge and Gross (1984)
within the adiabatic local or semi-local approximation
still suffers from various deficiencies like the underestimation of
charge transfer excitation energies or the absence of multiple electron
excitations in the spectrum Casida and Huix-Rotllant (2012). In order to describe excited states in the
framework of DFT, it is in principle not necessary to work within the
time-dependent regime. Various time-independent DFT approaches have been
investigated over the years, mostly at the formal
level Görling (1999); Levy and Nagy (1999); Gaudoin and Burke (2004); Ayers and Levy (2009); Ziegler et al. (2009); Ayers et al. (2012); Krykunov and Ziegler (2013).
In this paper, we will focus on ensemble DFT (eDFT) for excited
states Theophilou (1979); Gross et al. (1988a). The latter relies on the extension
of the variational principle to an ensemble of ground and excited states
which is characterized by a set of ensemble
weights Gross et al. (1988b). Note that Boltzmann
weights can be used Pastorczak et al. (2013) but it is not compulsory.
In fact, any set of ordered weights can
be considered Gross et al. (1988b). Since the ensemble
energy (i.e. the weighted sum of ground- and excited-state
energies) is a functional of the ensemble density, which is the
weighted sum of ground- and excited-state densities, a mapping between
the physical interacting and Kohn–Sham (KS) non-interacting ensembles can
be established. Consequently, a weight-dependent ensemble
exchange-correlation () functional must be introduced in order to obtain
the exact ensemble energy and, consequently, exact excitation
energies. Despite its formal simplicity (exact optical and KS gaps are easily
related in this context Gross et al. (1988a)) and advantages in contrast to
TD-DFT (it is straightforward to describe multiple excitations with an
ensemble), eDFT is not standard essentially because, so far, not much
effort has been
put in the development of
approximate functionals for ensembles. In particular, designing density-functional
approximations that remove the so-called ”ghost interaction”
error Gidopoulos et al. (2002),
which is induced by the ensemble Hartree energy, is
still challenging Pastorczak and Pernal (2014). Employing an ensemble exact exchange energy is of course
possible but then optimized effective potentials should in principle be
used, which is computationally demanding. Recently, accurate eDFT
calculations have been performed for the helium
atom Yang et al. (2014), the hydrogen
molecule Borgoo et al. (2015), and for two electrons in boxes or in a three-dimensional harmonic
well (Hooke’s atom) Pribram-Jones et al. (2014),
thus providing more insight into the ensemble energy and
potential. The key feature of the density functional in eDFT is
that it varies with the ensemble weight, even if the electron density is
fixed. This weight dependence plays a crucial role in the calculation of
the excitation energies Gross et al. (1988a). Developing
weight-dependent functionals is a complicated task that has not drawn
much attention so far. This explains why
eDFT is not a standard approach. There is clearly a need for models that can be solved
exactly in eDFT and, consequently, that can provide more insight into
the weight dependence of ensemble energies.

It was shown very recently Carrascal et al. (2015); Smith et al. (2016) that the nontrivial asymmetric Hubbard dimer can be used for understanding the limitations of standard approximate DFT in the strongly correlated regime and also for developing functionals in thermal DFT Smith et al. (2016). In the same spirit, we propose in this work to construct an exact eDFT for this model system. The paper is organized as follows. After a brief introduction to eDFT (Sec. II.1), a generalization of the adiabatic connection formalism to ensembles will be presented in Sec. II.2. The formulation of eDFT for the Hubbard dimer is discussed in Sec. III and exact results are given and analyzed in Sec. IV. Ground-state density-functional approximations are finally proposed and tested in Sec. V. Conclusions are given in Sec. VI.

## Ii Theory

### ii.1 Ensemble density functional theory for excited states

According to the Gross–Oliveira–Kohn (GOK) variational principle Gross et al. (1988b), that generalizes the seminal work of Theophilou Theophilou (1979) on equiensembles, the following inequality

(1) |

is fulfilled for any ensemble characterized by an arbitrary set (i.e. not necessarily a Boltzmann one) of weights with and a set of orthonormal trial -electron (with fixed) wavefunctions . The lower bound in Eq. (1) is the exact ensemble energy, i.e. the weighted sum of ground- and excited-state energies,

(2) |

where is the exact th eigenfunction of the Hamiltonian operator with energy and . A consequence of the GOK principle is that the ensemble energy is a functional of the ensemble density Gross et al. (1988a), i.e. the weighted sum of ground- and excited-state densities,

(3) |

Note that, in the standard formulation of eDFT Gross et al. (1988a), the additional condition is used so that the ensemble density integrates to the number of electrons. In the rest of this work, we will focus on non-degenerate two-state ensembles. In the latter case, a single weight parameter in the range can be used, since and , so that Eq. (1) becomes

(4) |

For convenience, the trial density matrix operator

(5) |

where and are orthonormal, has been introduced. denotes the trace and the ensemble energy equals

(6) |

For any electronic system, the Hamiltonian can be decomposed as where is the kinetic energy operator, denotes the two-electron repulsion operator, is the nuclear potential and is the density operator. Like in conventional (ground-state) DFT, the exact ensemble energy can be expressed variationally as follows Gross et al. (1988a),

(7) |

where

(8) | |||||

is the analog of the Levy–Lieb (LL) functional for ensembles. The minimization in Eq. (8) is performed over all ensemble density matrix operators with density ,

(9) |

Note that, according to the GOK variational principle, the following inequality is fulfilled for any local potential ,

(10) |

where is the ensemble energy of , so that the ensemble LL functional can be rewritten as a Legendre–Fenchel transform Eschrig (2003); Kutzelnigg (2006); van Leeuwen (2003); Lieb (1983); Franck and Fromager (2014); Borgoo et al. (2015),

(11) |

Note also that, in Eq. (7), the minimizing density is the exact physical ensemble density

(12) |

Like in standard ground-state DFT, the KS decomposition,

(13) |

is usually considered, where

(14) | |||||

is the non-interacting ensemble kinetic energy and is the (-dependent) ensemble Hartree-exchange-correlation functional. Applying the GOK principle to non-interacting systems leads to the following Legendre–Fenchel transform,

(15) |

where is the ensemble energy of . Combining Eq. (7) with Eq. (13) leads to the following KS expression for the exact ensemble energy,

(16) | |||||

The minimizing non-interacting ensemble density matrix in Eq. (16),

(17) |

reproduces the exact physical ensemble density,

(18) |

It is obtained by solving the self-consistent equations Gross et al. (1988a)

(19) |

As readily seen in Eq. (6), the exact (neutral) excitation energy is simply the first derivative of the ensemble energy with respect to the ensemble weight ,

(20) |

Using Eq. (16) and the Hellmann–Feynman theorem leads to

(21) | |||||

where By using Eq. (II.1), we finally obtain

(22) |

If the ground and first-excited states differ by a single electron excitation then the KS excitation energy (first term on the right-hand side of Eq. (22)) becomes the weight-dependent KS HOMO-LUMO gap . If, in addition, we use the decomposition

(23) |

where is the conventional (weight-independent) ground-state Hartree functional,

(24) |

we then recover the KS-eDFT expression for the excitation energy Gross et al. (1988a),

(25) |

where . Interestingly, in the limit, the excitation energy can be expressed exactly in terms of the usual ground-state KS HOMO-LUMO gap as

(26) |

As shown analytically by Levy Levy (1995) and illustrated numerically by Yang et al. Yang et al. (2014), corresponds to the jump in the potential when moving from (-electron ground state) to (ensemble of -electron ground and excited states). It is therefore a derivative discontinuity (DD) contribution to the optical gap that should not be confused with the conventional ground-state DD Stein et al. (2012); Kraisler and Kronik (2013, 2014); Gould and Toulouse (2014),

(27) |

where the fundamental gap is expressed in terms of , and ground-state energies as follows,

(28) |

For simplicity, we will also refer to the weight-dependent
quantity (see Eq. (25)) as DD.

Returning to the decomposition in Eq. (23), the contribution is usually split as follows,

(29) |

where

(30) |

is the exact ensemble exchange energy functional and is the non-interacting ensemble density matrix operator with density (see Eq. (14)). Consequently, according to Eqs. (8), (13) and (14), the ensemble correlation energy equals

(31) | |||||

### ii.2 Generalized adiabatic connection for ensembles

In order to construct the ensemble functional from the ground-state one (), Franck and Fromager Franck and Fromager (2014) have derived a generalized adiabatic connection for ensembles (GACE) where an integration over both the interaction strength parameter () and an ensemble weight in the range is performed. The major difference between conventional ACs Langreth and Perdew (1975); Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977); Savin et al. (2003); Nagy (1995) and the GACE is that, along a GACE path, the ensemble density is held constant and equal to when both and vary. Consequently, the integration over can be performed in the ground state while the deviation of the ensemble energy from the ground-state one is obtained when varying only. Formally, the GACE can be summarized as follows. Let us consider the Schrödinger,

(32) |

and KS

(33) |

equations where . The potentials and are adjusted so that the GACE density constraint is fulfilled,

(34) |

where

(35) |

and

(36) | |||||

According to Eqs. (13) and (23), the ensemble energy can be expressed as

(37) | |||||

where is the ground-state functional. Since and are the maximizing (and therefore stationary) potentials in the Legendre–Fenchel transforms of Eqs. (11) and (15) when , respectively, we finally obtain

(38) |

where the GACE integrand is simply equal to the difference in excitation energy between the interacting and non-interacting electronic systems whose ensemble density with weight is equal to :

(39) |

Note that, when the density equals the physical ensemble density
(see Eq. (12)) and , the GACE integrand
equals the
DD introduced in
Eq. (25).

An open and critical question is whether the GACE can actually be constructed for all weights in and densities of interest. In other words, does the GACE density constraint lead to interacting and/or non-interacting -representability problems ? So far, the GACE has been constructed only for the simple hydrogen molecule in a minimal basis and near the dissociation limit Franck and Fromager (2014), which basically corresponds to the strongly correlated symmetric Hubbard dimer. In the following, we extend this work to the nontrivial asymmetric Hubbard dimer. An important feature of such a model is that, in contrast to the symmetric case, the density (which is simply a collection of two site occupations) can vary, thus allowing for the construction of density functionals Carrascal et al. (2015); Smith et al. (2016).

## Iii Asymmetric Hubbard dimer

In the spirit of recent works by Carrascal et al. Carrascal et al. (2015) as well as Senjean et al. Senjean et al. (2016), we propose to apply eDFT to the asymmetric two-electron Hubbard dimer. The corresponding model Hamiltonian is decomposed as follows,

(40) |

where the two sites are labelled as 0 and 1, and is the hopping operator () which plays the role of the kinetic energy operator. The two-electron repulsion becomes an on-site repulsion,

(41) |

where is the spin-occupation operator. The last two contributions on the right-hand side of Eq. (40) play the role of the local nuclear potential. In this context, the density operator is . For convenience, we will assume that

(42) |

Note that the latter condition is fulfilled by any potential once it has been shifted by . Therefore, the final expression for the Hamiltonian is

(43) |

where

(44) |

In this work, we will consider the singlet two-electron ground and
first excited states for which analytical solutions exist (see
Refs. Carrascal et al. (2015); Smith et al. (2016) and the Appendix).
Note that, in order to yield the first singlet transition, the minimization in the GOK variational principle (see
Eq. (1)) can be restricted to singlet
wavefunctions, since singlet and triplet states are not coupled.
Consequently, eDFT can be formulated for singlet ensembles only. Obviously,
in He for example, singlet eDFT would not describe the lowest
transition . In the following, the first singlet excited
state (which is the excited state studied in this work) will be referred
to as ”first excited state” for simplicity.

For convenience, the occupation of site 0 is denoted and we have since the number of electrons is held constant and equal to 2. Therefore, in this simple system, the density is given by a single number that can vary from 0 to 2. Consequently, in this context, DFT becomes a site-occupation functional theory Chayes et al. (1985); Gunnarsson and Schönhammer (1986); Schönhammer et al. (1995); Capelle and Campo Jr. (2013) and the various functionals introduced previously will now be functions of . The ensemble LL functional in Eq. (8) becomes

(45) |

where the density constraint reads . By analogy with Eq. (11) and using , we obtain the following Legendre–Fenchel transform expression,

(46) | |||||

where and are the ground- and
first-excited-state energies of . Note that, even
though analytical expressions exist for the energies, has no
simple expression in terms of the density . Nevertheless, as readily
seen in Eq. (46), it can be computed exactly by performing
so-called Lieb maximizations. Note that an accurate parameterization has
been provided by Carrascal et al. Carrascal et al. (2015) for the
ground-state LL functional ().

Similarly, the ensemble non-interacting kinetic energy in Eq. (15) becomes

(47) | |||||

where and are the ground- and first-excited-state energies of the KS Hamiltonian

(48) |

From the simple analytical expressions for the HOMO and LUMO energies,

(49) |

and

(50) |

it comes that

(51) |

and

(52) |

According to the Hellmann–Feynman theorem, combining Eqs. (48) and (52) leads to

(53) |

where is the first singlet (two-electron) excited state of . Therefore, the density (i.e. the occupation of site 0) in the non-interacting first excited state is equal to 1 for any and values, as illustrated in the top left-hand panel of Fig. 1. Consequently, a density will be ensemble non-interacting representable in this context if it can be written as where the non-interacting ground-state density varies in the range (see the top left-hand panel of Fig. 1), thus leading to the non-interacting representability condition

(54) |

or, equivalently,

(55) |

For such densities, the maximizing KS potential in Eq. (47) equals

(56) |

and, consequently, the ensemble non-interacting kinetic energy functional can be expressed analytically as follows,

(57) |

The ensemble correlation energy, which is the key quantity studied in this work, is defined as follows,

(58) |

where the Hartree energy equals Smith et al. (2016)

(59) | |||||

Note that the latter expression is simply obtained from the conventional one in Eq. (24) by substituting a Dirac-delta interaction with strength for the regular two-electron repulsion,

(60) |

and by summing over sites rather than integrating over the (continuous) real space. The exact ensemble exchange energy in Eq. (30) becomes in this context

(61) | |||||

thus leading, according to the Appendix, to the analytical expression

(62) | |||||

where

(63) |

is the ground-state exchange energy for two unpolarized electrons. Note that the exchange contribution to the GACE integrand (see Eq. (38)) will therefore have a simple analytical expression,

(64) | |||||

Finally, the maximizing potential in Eq. (46) which reproduces the ensemble density fulfills, according to the inverse Legendre–Fenchel transform,

(65) |

where the minimizing density is . Therefore,

(66) |

and, since (see Eqs. (56) and (57))

(67) |

the ensemble Hartree- potential reads

(68) | |||||

As a result, the ensemble correlation potential can be calculated exactly as follows,

(69) | |||||

where all contributions but have an analytical expression. The Hartree potential equals and, according to Eq. (62), the ensemble exchange potential reads

(70) | |||||

Note the unexpected minus sign on the right-hand side of Eq. (68). It originates from the definition of the potential difference (see Eq. (44)) and the choice of (occupation of site 0) as variable, the occupation of site 1 being . Therefore, can be rewritten as and

(71) | |||||

Note finally that, as readily seen in Eq. (70), the ensemble potential can be expressed in terms of the ground-state potential () and the ensemble weight. This simple relation, which is transferable to ab initio Hamiltonians, could be used for developing ”true” approximate weight-dependent density-functional potentials.

## Iv Exact results

### iv.1 Interacting ensemble density and derivative discontinuity

In the rest of the paper, the hopping parameter is set to . For clarity, we shall refer to the local potential in the physical (fully-interacting) Hubbard Hamiltonian as . This potential is the analog of the nuclear-electron attraction potential in the ab initio Hamiltonian. The corresponding ensemble density is the weighted sum of the ground- and excited-state occupations of site 0,

(72) |

where, according to the Hellmann–Feynman theorem,

(73) |

Note that the first-order derivative of the energies with respect to
can be simply expressed in terms of the
energies (see Eq. (124)) and that, for a fixed
value, the ensemble density varies linearly with . Ground- and excited-state densities are shown in Fig. 1.
For an arbitrary potential value , in the
weakly correlated regime (), site occupations are close to 2 or
0 in the ground state and they become equal to 1
in the first excited state. Therefore, in this case, the model describes a charge
transfer excitation. On the other hand, in the strongly correlated
regime (), the ground-state density will be close to 1
(symmetric case).
When is large,
small changes in around cause large changes in the
excited-state density. As clearly seen from the Hamiltonian
expression in Eq. (43), when , site 0 ”gains” an
electron when the lowest (singlet) transition occurs if whereas, if , it ”loses” an electron. This explains why
the excited-state density curves approach a discontinuous limit at
when . Let us stress that, for large but finite
values, the latter density will vary rapidly and continuously from 0 to
2 in the vicinity of while the ground-state density remains
close to 1. This observation will enable us to interpret the
GACE integrand in the following.

Turning to the calculation of the DD (see Eq. (25)), the latter can be obtained in two ways, either by taking the difference between the physical and KS

(74) |

excitation energies, which gives

(75) |

or by differentiation,

(76) |

In the former case, we obtain from Eqs. (49), (50), and (56) the analytical expression

(77) | |||||

Regarding Eq. (76), the -dependent ensemble energy must be determined numerically by means of a Legendre–Fenchel transform calculation (see Eqs. (46) and (58)) and its derivative at is then obtained by finite difference. As illustrated in the right-hand top panel of Fig. 2, the two expressions are indeed equivalent. In the symmetric Hubbard dimer (), it is clear from Eq. (77) that the DD is weight-independent, since , and it is equal to . In this particular case, the ground and first-excited states actually belong to different symmetries. In the asymmetric case, various patterns are obtained (see Fig. 2). Interestingly, the ”fish picture” obtained by Yang et al. Yang et al. (2014) for the helium atom is qualitatively reproduced by the Hubbard dimer model when , except in the small- region where a sharp change in the DD (with positive slope) is observed for the helium atom. This feature does not occur in the two-site model. From the analytical expression,

(78) |

and Fig. 1, it becomes clear that, in the Hubbard dimer, the DD will systematically decrease with . Variations in and for various weights are shown in Figs. 3 and 4, respectively. When , is close to (according to Fig. 1) and, since the on-site repulsion becomes a perturbation, the DD can be well reproduced by the exchange-only contribution. Thus, according to Eq. (64), we obtain

(79) |

As readily seen in Eq. (79), the DD is close to zero for small weights and, when , it equals , which is in agreement with both Figs. 3 and 4. On the other hand, when , the physical energies are expanded as follows, according to Eq. (116),

(80) | |||||

thus leading to the following expansions for the derivatives,

(81) | |||||

and, according to Eqs. (72) and (73), to the following expansion for the ensemble density,

(82) |

As readily seen in Eq. (IV.1), the ensemble density is close to 1 in the small- region. Consequently, according to Eqs. (77) and (IV.1), the DD varies as , which is in agreement with the panel of Fig. 3 and the panel of Fig. 4. On the other hand, when , it comes from Eq. (IV.1),

(83) | |||||

thus leading to the following expansion for the equiensemble DD,

(84) |

The latter expansion matches the behavior observed in the and panels of Fig. 3 as well as and panels of Fig. 4, when . Note finally that, in the panel of Fig. 3, the equiensemble DD is highly sensitive to changes in around when . In the latter case, the ground-state density remains close to 1 (symmetric dimer), as shown in Fig. 1, and the DD becomes

(85) | |||||