Density driven correlations in ensemble density functional theory: insights from simple excitations in atoms

Gross-Oliveira-Kohn (GOK) derived a series of proofsGross, Oliveira, and Kohn (1988a, b); Oliveira, Gross, and Kohn (1988) that form the foundation of ensemble DFT (EDFT). A consequence of GOK’s work is that we can write an energy density functional that finds ${\cal E}$ through a variational principle based on the ensemble density:

	$\displaystyle{\cal E}^{\boldsymbol{w}}=$	$\displaystyle\min_{n}\big{\{}{\cal F}^{\boldsymbol{w},1}[n]+\int n(\boldsymbol{r})v_{{\text{Ext}}}(\boldsymbol{r})d\boldsymbol{r}\big{\}}$		(5)
	$\displaystyle=$	$\displaystyle\sum_{\kappa}w_{\kappa}E_{\kappa}\;.$		(6)

The set of ensemble weights $\boldsymbol{w}=\{w_{\kappa}\}$ form an additional input, on top of the external potential $v_{{\text{Ext}}}$. The energy is thus a function of $\boldsymbol{w}$ and a functional of $v_{{\text{Ext}}}$.

Eq. (5) involves energies $E_{\kappa}=\langle\kappa|\hat{H}|\kappa\rangle$ which are eigenvalues of $\hat{H}=\hat{T}+\hat{W}+\hat{v}_{{\text{Ext}}}$ for state $|\kappa\rangle$. The density $n^{\boldsymbol{w}}=\text{Tr}[\hat{\Gamma}\hat{n}]=\sum_{\kappa}w_{\kappa}n_{\kappa}$ that minimizes (5) is the density of $\hat{\Gamma}=\sum_{\kappa}w_{\kappa}|\kappa\rangle\langle\kappa|$. Note that $\hat{\Gamma}$ is passive, so that the smallest eigenvalues are always associated with the largest weights.

The universal functional ${\cal F}^{\boldsymbol{w},1}[n]$ represents the special $\lambda=1$ case of:

	$\displaystyle{\cal F}^{\boldsymbol{w},\lambda}[n]=$	$\displaystyle\min_{\hat{\Gamma}\to n,\boldsymbol{w}}{{\text{Tr}}\left[\hat{\Gamma}(\hat{T}+\lambda\hat{W})\right]}$
	$\displaystyle\equiv$	$\displaystyle{{\text{Tr}}\left[\hat{\Gamma}^{n^{\boldsymbol{w}},\lambda}(\hat{T}+\lambda\hat{W})\right]},$		(7)

where $\hat{\Gamma}\to n,\boldsymbol{w}$ means $\text{Tr}[\hat{\Gamma}\hat{n}(\boldsymbol{r})]=n(\boldsymbol{r})$ and the weights of $\hat{\Gamma}$ are given by $\boldsymbol{w}$. Here, $\text{Tr}[\hat{\Gamma}^{n^{\boldsymbol{w}},\lambda}\hat{n}]\equiv n^{\boldsymbol{w}}$. Except when clarification is required, we shall now proceed to drop explicit mention of $\boldsymbol{w}$ on functionals and EDMs.

Furthermore, the GOK theorems guarantee that there is a unique mapping $n\to v^{\lambda}[n]$ for any given set of ensemble weights $\boldsymbol{w}$ – at least for “well-behaved” densities which we shall assume throughout. Thus, assuming that interacting and non-interacting ensemble $v$-representability is not a problem, we can define

$\displaystyle\hat{\Gamma}^{\lambda}\equiv$

$\displaystyle\sum_{\kappa}w_{\kappa}|\kappa^{\lambda}\rangle\langle\kappa^{\lambda}|,$

(8)

where $|\kappa^{\lambda}\rangle$ is the $\kappa$th eigenfunction of $\hat{H}^{\lambda}=\hat{T}+\lambda\hat{W}+\hat{v}^{\lambda}$. In the case of equi-ensembles (see next section), these EDMs are unique, but this is not necessarily so in general.

Taking $\lambda\rightarrow 0$, and working with equi-ensembles to preserve symmetries, we can identify $\hat{H}^{0}=\hat{T}+\hat{v}^{0}$ and the corresponding non-interacting states $|\kappa_{s}\rangle\equiv|\kappa^{0}\rangle$. These states are of the form of Slater-determinants or, when necessary, a superposition thereof (more below). They are formed on the same set of single-particle orbitals, obeying self-consistent equations,

$\displaystyle\big{\{}\hat{t}+v^{0}[n](\boldsymbol{r})\big{\}}\phi_{i}[n](\boldsymbol{r})=$

$\displaystyle\epsilon_{i}[n]\phi_{i}[n](\boldsymbol{r})\;,$

(9)

Here $\hat{t}=-\frac{1}{2}\nabla^{2}$ is the one-body kinetic energy operator, and $v^{0}$ is a spin-unpolarised potential. The orbitals, $\phi_{i}$, thus have the same spatial form for $\mathord{\uparrow}$ and $\mathord{\downarrow}$ spin-channels, from which it follows that the ensemble density is given by,

$\displaystyle n(\boldsymbol{r})=\sum_{\kappa}w_{\kappa}n_{s,\kappa}(\boldsymbol{r})\equiv\sum_{i}f_{i}|\phi_{i}(\boldsymbol{r})|^{2}\;,$

(10)

where $0\leq f_{i}\leq 2$ is an average occupation factor. The factors, $f_{i}=\sum_{\kappa}w_{\kappa}\theta_{i}^{\kappa}$, are derived from the orbital expansion of the non-interacting densities,

$\displaystyle n_{s,\kappa}(\boldsymbol{r})=\sum_{i}\theta_{i}^{\kappa}|\phi_{i}(\boldsymbol{r})|^{2}\;,$

(11)

where $\theta_{i}^{\kappa}\in\{0,1,2\}$ is the occupancy of orbital $i$ in $|\kappa_{s}\rangle$. The potential

$\displaystyle v_{s}[n](\boldsymbol{r})\equiv v^{0}[n](\boldsymbol{r})$

(12)

is known as the (ensemble) Kohn-Sham potential. For the equi-ensembles considered in this work, $v_{s}[n]$ is unique.

Through plausible hypotheses on the ordering of the $\lambda$-dependent eigenvalues for $\lambda\rightarrow 0$, we can define the ensemble versions of the Kohn-Sham kinetic energy, Hartree-exchange energy, and correlation energy. Respectively, these areGould and Pittalis (2017):

	$\displaystyle{\cal T}_{s}[n]\equiv$	$\displaystyle{\cal F}^{0}[n]=\sum_{\kappa}w_{\kappa}T[n_{s,\kappa}]\;,$		(13)
	$\displaystyle{\cal E}_{\text{Hx}}[n]\equiv$	$\displaystyle\lim_{\eta\to 0^{+}}\frac{{\cal F}^{\eta}[n]-{\cal T}_{s}[n]}{\eta}\;=\sum_{\kappa}w_{\kappa}\Lambda_{\text{Hx},\kappa}[n_{s,\kappa}]\;.$		(14)
	$\displaystyle{\cal E}_{\text{c}}[n]=$	$\displaystyle{\cal F}^{1}[n]-{\cal T}_{s}[n]-{\cal E}_{\text{Hx}}[n]\;.$		(15)

Here, $T_{s,\kappa}$ and $\Lambda_{\text{Hx},\kappa}$ are orbital-dependent kinetic and Hartree-exchange-like energy termsGould and Pittalis (2017). As long as we are concerned with the kinetic energy and density, a restriction to single Slater-determinant may be harmless and practical. But when evaluating the Hartree-exchange energy, linear combination of Slater-determinants must be admitted.Gould and Pittalis (2017) Such a formalism maximally avoids “ghost-interactions”,Brandi, Matos, and Ferreira (1980); Gidopoulos, Papaconstantinou, and Gross (2002a); Pastorczak and Pernal (2014); Yang et al. (2014) which can cause problems in EDFT.

Next, turning to the correlation energy functional, we recognise that it is more complex in EDFT than in DFT because – besides the usual reasons – it has an additional term that is zero in pure states.Gould and Pittalis (2019) The additional term is required to account for the fact that a state $\kappa$ is associated with an interacting density $n_{\kappa}\equiv n_{\kappa}^{1}$ that differs from the corresponding non-interacting density $n_{s,\kappa}\equiv n_{\kappa}^{0}$. Only in pure states or in certain ensembles are these two densities guaranteed to be the same.

As a result, we arrive at

$\displaystyle{\cal E}_{\text{c}}[n]=$

$\displaystyle\sum_{\kappa}w_{\kappa}\{E_{\text{c},\kappa}^{{\text{SD}}}[n_{\kappa},n]+E_{\text{c},\kappa}^{{\text{DD}}}[n_{\kappa},n]\}$

(16)

where $E_{\text{c},\kappa}^{{\text{SD}}}[n_{\kappa},n]$ is a state-driven (SD) contribution – a bifunctional – which resembles the correlation term in pure state DFT; and $E_{\text{c},\kappa}^{{\text{DD}}}[n_{\kappa},n]$ is a density-driven (DD) contribution – another bifunctional – that is unique to ensembles.Gould and Pittalis (2019) Thus, it is the difference in the densities $n_{\kappa}$ and $n_{s,\kappa}$ that gives rise to the DD term in ensembles. Sometimes, it also makes the SD term difficult to pin down. In the cases studied in the following sections the SD terms are easily identified.

In order to formally describe the terms in Eq. (16) let us introduce, for each interacting state $\kappa$ in the ensemble, single-particle orbitals $\psi_{i}^{\kappa}$ which obey a KS-like equation (9) but with a state-dependent potential $v^{\kappa}_{s}$ that ensures,

$\displaystyle n_{\kappa}(\boldsymbol{r})=\sum_{i}\theta_{i}^{\kappa}|\psi_{i}^{\kappa}(\boldsymbol{r})|^{2}\;$

(17)

analogous to Eq. (11). In this step, restriction to single Slater-determinants is admissible. More details are given in Refs Ayers, Levy, and Nagy, 2015; Gould and Pittalis, 2019.

Hence, we are ready to write

	$\displaystyle E_{c,\kappa}^{{\text{SD}}}[n_{\kappa},n]\equiv$	$\displaystyle\langle\kappa|\hat{T}+\hat{W}|\kappa\rangle-F_{{\text{EXX}},\kappa}[\{\psi_{i}^{\kappa}\}],$		(18)
	$\displaystyle E_{c,\kappa}^{{\text{DD}}}[n_{\kappa},n]\equiv$	$\displaystyle F_{{\text{EXX}},\kappa}[\{\psi_{i}^{\kappa}\}]-F_{{\text{EXX}},\kappa}[\{\phi_{i}\}],$		(19)

where $F_{{\text{EXX}},\kappa}=T_{s,\kappa}+\Lambda_{\text{Hx},\kappa}$ is the contribution of the state $\kappa$ to the exact exchange (EXX) component of ${\cal F}^{1}$. Finally, we can also define

$\displaystyle{\cal F}^{1}=$

$\displaystyle\sum_{\kappa}w_{\kappa}\{F_{{\text{EXX}},\kappa}+E_{\text{c},\kappa}^{{\text{SD}}}+E_{\text{c},\kappa}^{{\text{DD}}}\}\;,$

(20)

as an alternative expression for the interacting universal functional in eq. (7) at $\lambda=1$.

As complicated as the bifunctionals of Eq. (18) and Eq. (19) might appear, in Ref. Gould and Pittalis, 2019 we have shown an approximate way to deal with them practically. In Sec. III, we shall introduce another way which allows us to compute approximate DD correlations without having to work with bifunctional explicitly.

Although ensemble DFT may be formulated for general weights, exploitation and fulfillment of symmetries imply some restrictions. In order to explain this aspect and also to introduce the reader to the procedure implemented in our numerical results, let us start by briefly reviewing a few important facts.

Firstly, the non-relativistic many-electron Hamiltonian does not depend on the spin operators. Thus, trivially, it is invariant under spin rotations and the many-electron states may be chosen as eigenstates of $\hat{S}^{2}$ where $\hat{S}$ is the the total spin operator.

Secondly, some systems can also exhibit additional symmetry in real space. The Hamiltonian of an isolated atom, for example, is invariant under rotations in real space. This is the reason why atomic eigenstates may be chosen also as eigenstates of $\hat{L}^{2}$ where $\hat{L}$ is the total angular momentum operator. Relative to a given nuclear configuration, molecules can be invariant under the symmetry operations of point groups. Crystal requires consideration of space groups, which add (discrete) translational invariance to the point groups. Therefore the spatial symmetry of a state is also specified by an irreducible representation, $\alpha$, of a space symmetry group (which can, trivially, be the identity).

To further specify the terminology, let $\hat{H}$ be a Hamiltonian that is invariant under the transformations $\hat{G}$ that forms a group $\mathcal{G}$. Provided that there are no “accidental” degeneracy, every group of degenerate eigenstates provides an irreducible representation (Irrep) of $\mathcal{G}$ with the dimension of the degeneracy. Such a representation is “irreducible” because the degenerate manifold is an invariant under the action of $\mathcal{G}$ and contains no other invariant subspace.

Thus, strictly, the corresponding eigenstates $\Psi^{\alpha}_{l}$ of $\hat{H}_{\mathcal{G}}$ are not invariant. They transforms as

$\displaystyle\hat{G}\left|\Psi^{\alpha}_{l}\right\rangle=$

$\displaystyle\sum_{m}D^{\alpha}_{ml}\left|\Psi^{\alpha}_{m}\right\rangle$

(21)

where the upper index $\alpha$ labels the Irrep and the lower index denotes the partners in the basis for the invariant subspace. $D^{\alpha}_{ml}$ is a matrix that provides the actual representation of $\hat{G}$. States that transforms as the basis of a Irrep are said to be symmetry-adapted, whether or not they may be the exact states of the system Hamiltonian.

The Hamiltonian we work with is also invariant under exchange of particles. The fermionic nature of the corresponding many-body states can be taken into account by means of Slater determinants (Sldet) built from orthonormal single-particle states. These determinants can be combined linearly to get spin symmetry-adapted many-particle states – i.e., configuration state functions – to provide basis functions through which we may represent either approximate or exact many-particle states Helgaker, Jørgensen, and Olsen (2002). Both for computational and interpretational convenience, especially when dealing with excited multiplets, it is best to work with single-particle states that are symmetry adapted as well.

The need for introducing symmetry-adaptation when accessing excited states through ensemble density functional methods has been emphasized by Theophilou (and various collaborators) Theophilou and Gidopoulos (1995); Theophilou (1997); Theophilou and Papaconstantinou (2000). In these works, weights were assigned equally to all the states. Thus, in order to access the individual energies of the levels of interest, several computations would be required which each involve different ensemble energy functionals Theophilou (1979); Hadjisavvas and Theophilou (1985). GOKGross, Oliveira, and Kohn (1988b); Oliveira, Gross, and Kohn (1988) showed that the choice of the weights could be made quite flexible.¹¹1Recently, Fromager and coworkersDeur and Fromager (2019); Senjean and Fromager (2018) have pointed out that the flexibility of GOK can be exploited even further, in such a way to access single excited levels from calculations that need to refer only to a single ensemble-density-functional Equal weights, however, were still employed for degenerate states in the same manifold. Use of equal weights is important because the particle densities of $|\Psi^{\alpha}_{m}\rangle$ does not necessarily transform according the same Irrep to which the pure state belong. Therefore, invariant Kohn-Sham potentials can be obtained by forming ensembles that are totally invariant by construction.

To be more explicit, let us consider the equi-ensemble for a degenerate set of ground-states,

$\displaystyle\hat{\Gamma}=\frac{1}{d_{\alpha}}$

$\displaystyle\sum_{l=1}^{d_{\alpha}}|\Psi^{\alpha}_{l}\rangle\langle\Psi^{\alpha}_{l}|$

(22)

where $d_{\alpha}$ is the dimension (degeneracy) of $\alpha$. Through use of the orthogonality theorem, it becomes apparent that this ensemble is invariant for any $\hat{G}\in\mathcal{G}$. The corresponding ensemble density, $n=\text{Tr}[\hat{n}\hat{\Gamma}]$ inherits this invariance. Thus, the corresponding KS potential has the symmetry of the external potential by construction. Ensembles for excited states can be formed along similar lines. Explicit examples are given in the next sections for atoms.

Furthermore, because the Hamiltonian is an invariant w.r.t. $\mathcal{G}$, it transforms states from invariant subspaces to states of the same subspaces. From this it readily follows that both the regular variational principle and ensemble-variational principles may be applied equally well to states which are symmetry adapted. Equivalently, we may work with the projected Hamiltonian

$\displaystyle\hat{H}_{\mathcal{P}}=\hat{P}\hat{H}\hat{P}$

(23)

where $\hat{P}$ is provided, essentially, by the same procedure used to extract a symmetry-adapted basis from a set of generic states.

Among the advantages of this approach are:

1. (a)

   When we are concerned with the lowest states of some symmetry species, such as the lowest triplet or the lowest ${}^{1}D$ level in atoms, we do not need to consider any exited-states formulation of DFTGunnarsson and Lundqvist (1976);

2. (b)

   When we are concerned with the lowest two (or more) excited states of some symmetry species, such as the optical gap between singlet states, we can ignore all the states in betweenTheophilou and Papaconstantinou (2000); Gidopoulos, Papaconstantinou, and Gross (2002b); Pribram-Jones et al. (2014);

3. (c)

   We can find passive EDMs that are stationary on $\hat{H}_{\mathcal{P}}$ and thus satisfy the weight ordering required by the GOK formulation despite other “orthogonal” states having lower energies. We use “$\mathcal{P}$-passive” to denote such EDMs. The states described above in (a) are the simplest example of $\mathcal{P}$-passive states.

Finally, we note that the above discussion allows linear combinations of Sldets (configuration state functions) at the non-interacting KS level, so long as the terms involved all have the same energy and density. This gives rise, through the formalism outlined previously by the authorsGould and Pittalis (2017), to non-vanishing, spontaneously emergent singlet-triplet splitting already at the Hartree-exchange level (mentioned briefly in the previous section). For the purpose of the present work, however, we average over such spin splittings and, thus, can consistently reduce to single Sldets, as detailed later.

Specifically, our cases involve an excitation between $2s$–$2p$ (Li/Be) or $3s$–$3p$ (Na/Mg) orbitals. Such transitions are very interesting, per the results of the previous section, as they allow calculations to be carried out across all weights, not just passive states. This is possible because the limit $W=1$ is $\mathcal{P}$-passive and thus amenable to an EDFT treatment. Thus we can connect (via the ensemble) states that belong to different Irreps. Importantly, the states at $W=0$ and $W=1$ may be regarded, effectively, as being ground states.

In order to illustrate the EDMs we work with, let us provide the details for the non-interacting ones,

$\displaystyle\hat{\Gamma}_{s}^{W}=$

$\displaystyle(1-W)\hat{\Gamma}_{s,S_{0}}+W\hat{\Gamma}_{s,P_{0}}\;,$

(24)

Here, $S_{0}$ refers to the equi-ensemble of states that have angular momentum $L=0$, which can be singlet or a degenerate doublet. $P_{0}$ refers to the equi-ensemble of states that have have angular momentum $L=1$, which can be singlet or triplets. The index “0” refers to the fact that we are restricting to the lowest states in energy for each $L$. As a result, we can reduce to work with ensemble of single Sldets.

In Eq. (24), the first term involves an average over a (degenerate) doublet for Li:

$\displaystyle\hat{\Gamma}_{s,S_{0}}^{\text{Li}}=\frac{1}{2}\sum_{\sigma\in\mathord{\uparrow},\mathord{\downarrow}}|2s^{\sigma}\rangle\langle 2s^{\sigma}|\;;$

(25)

and is a singlet for Be:

$\displaystyle\hat{\Gamma}_{s,S_{0}}^{\text{Be}}=|2s^{2}\rangle\langle 2s^{2}|\;.$

(26)

The second term involves averaging over three real- (or, equivalently, complex-) valued $p$-orbitals. It also involves averaging over a spin-doublet for Li:

$\displaystyle\hat{\Gamma}_{s,P_{0}}^{\text{Li}}=\frac{1}{6}\sum_{\mu=x,y,z}\sum_{\sigma\in\mathord{\uparrow},\mathord{\downarrow}}|2p_{\mu}^{\sigma}\rangle\langle 2p_{\mu}^{\sigma}|\;;$

(27)

and equally averaging over the singlet plus triplet for Be:

$\displaystyle\hat{\Gamma}_{s,P_{0}}^{\text{Be}}=\frac{1}{12}\sum_{\mu=x,y,z}\sum_{\sigma,\sigma^{\prime}\in\mathord{\uparrow},\mathord{\downarrow}}|2s^{\sigma}2p_{\mu}^{\sigma^{\prime}}\rangle\langle 2s^{\sigma}2p_{\mu}^{\sigma^{\prime}}|\;.$

(28)

In the expressions above, the core configuration $1s^{2}$ is implied. The expressions for Na and Mg can be obtained by replacing $1s^{2}\to 1s^{2}2s^{2}$, $2s\to 3s$, and $2p\to 3p$. For brevity the notation does not emphasize that the non-interacting states also depend on $W$ (interacting quantities, of course, do not depend on $W$).

Finally note an important point: in the expressions above, we have also used the fact that given a singlet-state (ss) and three triplet states (ts_-1, ts₀ and ts₁) involving any pair of single-particle orbitals, we can write (in short-hand) $\frac{1}{4}\big{[}|\text{ss}\rangle\langle\text{ss}|+\sum_{s_{z}}|\text{ts}_{s_{z}}\rangle\langle\text{ts}_{s_{z}}|\big{]}=\frac{1}{4}\sum_{\sigma,\sigma^{\prime}\in\mathord{\uparrow},\mathord{\downarrow}}|\sigma\sigma^{\prime}\rangle\langle\sigma\sigma^{\prime}|$, where the latter are the four Sldets formed on the same set of spin-restricted spatial orbitals. Working with such averaged states is acceptable here, because we are interested in excitations behaviours between different spatial Irreps, rather than excitation behaviours between different spin states.

It then follows from the properties of ensembles that the interacting energy,

$\displaystyle{\cal E}^{W}\equiv{\cal E}^{W}[n^{W}]=$

$\displaystyle(1-W){\cal E}_{S_{0}}+W{\cal E}_{P_{0}},$

(29)

and the density,

$\displaystyle n^{W}(r)=$

$\displaystyle(1-W)n_{S_{0}}(r)+Wn_{P_{0}}(r)\;,$

(30)

are piece-wise linear. Note, all densities depend only on $r=|\boldsymbol{r}|$ due to preservation of fundamental symmetries by equi-ensembles. Here, we reintroduce a weight superscript, $W$, to make explicit that we are referring to an ensemble like eq. (24) taken with excitation weights $\boldsymbol{w}=\{1-W,W\}$ and equi-ensembles, per previous paragraphs.

Next, we consider the exact exchange (EXX) contribution to the universal functional

	$\displaystyle{\cal F}_{{\text{EXX}}}^{W}[n^{W}]=$	$\displaystyle{\cal F}^{1,W}[n^{W}]-{\cal E}_{\text{c}}^{W}[n^{W}]$		(31)
	$\displaystyle=$	$\displaystyle{\cal T}_{s}^{W}[n^{W}]+{\cal E}_{\text{Hx}}^{W}[n^{W}]$		(32)
	$\displaystyle\equiv$	$\displaystyle(1-W)F_{{\text{EXX}},S_{0}}[\{\phi_{i}[n^{W}]\}]$
		$\displaystyle+WF_{{\text{EXX}},P_{0}}[\{\phi_{i}[n^{W}]\}]\;.$		(33)

Here, the orbital functional ${\cal F}_{{\text{EXX}},S_{0}}$ is as defined in previous workGould and Dobson (2013) and detailed in Section IV. Extension to $P_{0}$ is described in Section IV. Note, in the final expression we use $F$, rather than calligraphic ${\cal F}$, to indicate that the ground state ensembles are treated like “pure” states. We will further clarify this choice later.

The ensemble exact exchange (EEXX) approximation then involves finding,

	$\displaystyle{\cal E}_{{\text{EXX}}}^{W}=$	$\displaystyle\min_{n}\big{\{}{\cal F}_{{\text{EXX}}}[n]+{\cal E}_{{\text{Ext}}}[n]\big{\}}$
	$\displaystyle\equiv$	$\displaystyle{\cal F}_{{\text{EXX}}}^{W}[n^{W}_{{\text{EXX}}}]+{\cal E}_{{\text{Ext}}}[n^{W}_{{\text{EXX}}}]\;,$		(34)

by variational principles, using ${\cal E}_{{\text{Ext}}}[n]=\int n(\boldsymbol{r})v_{{\text{Ext}}}(\boldsymbol{r})d\boldsymbol{r}$. Our second expression introduce the density obtained self-consistently within the EEXX approximation, $n_{{\text{EXX}}}^{W}$, which is not the same as $n^{W}$.

The cases $W=0$ and $W=1$ are both $\mathcal{P}$-passive, and thus may be treated as ground states. In these cases, EXX gives densities,

$\displaystyle n_{{\text{EXX}}}^{W=0}(r)\approx$

$\displaystyle n_{S_{0}}(r),$

$\displaystyle n_{{\text{EXX}}}^{W=1}(r)\approx$

$\displaystyle n_{P_{0}}(r),$

(35)

that are close to their exact counterparts in the atomic systems we consider. Consequently, we can approximately treat $n_{{\text{EXX}}}^{W=0,1}$ as the true density $n^{W=0,1}$. It then follows from ${\cal E}={\cal E}_{{\text{EXX}}}+{\cal E}_{\text{c}}$ that the energies,

$\displaystyle{\cal E}_{{\text{EXX}}}^{W=0}\approx$

$\displaystyle E_{S_{0}}-E_{\text{c},S_{0}}^{{\text{SD}}},$

$\displaystyle{\cal E}_{{\text{EXX}}}^{W=1}\approx$

$\displaystyle E_{P_{0}}-E_{\text{c},P_{0}}^{{\text{SD}}},$

(36)

are correct up to SD correlation terms $E_{c}^{{\text{SD}}}$ and small errors from the densities. Here we neglect DD correlations at $W=0$ and $1$. In fact, there may be some DD correlations associated with (actual or effective) ground state degeneracies – we do not concern ourselves with these to focus on only the DD correlations that affect excitations directly.

For the cases $0<W<1$ the situation becomes more complicated. This is because the density $n_{{\text{EXX}}}^{W}$ that minimises Eq. (34) is not the right ensemble density, i.e.

$\displaystyle n_{{\text{EXX}}}^{W}\neq{(1-W)n_{{\text{EXX}}}^{W=0}+Wn_{{\text{EXX}}}^{W=1}}\;.$

(37)

Thus, when ${\cal E}_{{\text{EXX}}}^{W}[n_{{\text{EXX}}}^{W}]$ is employed as an estimation of the exact energy it has two sources of error:

1. 1.

   It misses SD and DD correlation terms as these are not included in the EXX approximation;

2. 2.

   It is not obtained at the correct density which introduces an additional source of error.

We will discuss these points further below.

First, however, we will consider how this approach affects energy terms. Figure 1 shows the deviation from piece-wise linear,

$\displaystyle\Delta{\cal E}^{W}=$

$\displaystyle{\cal E}^{W}-[(1-W){\cal E}^{W=0}+W{\cal E}^{W=1}].$

(38)

of various energy quantities. Here, energies and densities are obtained approximately self-consistently, by solving the Krieger-Li-Iafrate (KLI) approximation for the self-consistent EXX potentialKrieger, Li, and Iafrate (1992). Further details of calculations are provided in Methodology.

As the goal of these self-consistent calculations is to minimise the total EEXX energy [Eq. (34)] neither energies nor densities are expected to be piece-wise linear. It is thus notable that, despite no requirement to be piece-wise linear, the EXX energy ${\cal E}_{{\text{EXX}}}$ barely deviates from piece-wise linear. Adapting arguments laid out by Gould and DobsonGould and Dobson (2013) shows that $\Delta{\cal E}_{{\text{EXX}}}^{W}$ should be concave, but says nothing about how concave it should be. In fact, a very small amount of concavity can be seen for Be and Na. In Li and Mg the concavity is invisible to the eye. In all cases it is tiny compared to deviations of other quantities.

Other energy deviations can be convex or concave, as they are only constrained by the fact that ${\cal E}_{{\text{EXX}}}$ must be concave. Indeed, we see substantial variations in all cases except Li. In Mg, the maximum deviation of ${\cal E}_{\text{Hx}}$ and ${\cal E}_{{\text{Ext}}}$ are more than half the vertical gap ($\text{Gap}=E_{{\text{EXX}},P_{0}}-E_{{\text{EXX}},S_{0}}$ – equivalent to the optical gap in atomic systems) in magnitude, but have equal and opposite errors which thus cancel. Figure 2 shows that the deviation of the density (Mg is shown, but other atoms are qualitatively similar) is almost exclusively found around the nucleus.

Above, we identified two sources of error found when one obtains an energy and density using (34). We can directly remove the second source (errors caused by deviation from piece-wise-linearity of densities) by using “density inversion”. That is, by seeking a potential,

$\displaystyle v_{s}[n^{W}]\to n^{W}=(1-W)n^{W=0}+Wn^{W=1},$

(39)

that yields a set of orbials $\phi_{i}[n^{W}]$ that give our target piece-wise linear density $n=\sum_{i}f_{i}|\phi_{i}|^{2}$. This lets us obtain ${\cal F}_{{\text{EXX}}}^{W}[n^{W}]$ and ${\cal E}_{{\text{EXX}}}^{W}[n^{W}]={\cal F}_{{\text{EXX}}}^{W}[n^{W}]+{\cal E}_{{\text{Ext}}}[n^{W}]>{\cal E}_{{\text{EXX}}}^{W}[n_{{\text{EXX}}}^{W}]$, and thus remove any error caused by deviations from piece-wise linearity of the density. Note, the external energy ${\cal E}_{{\text{Ext}}}[n^{W}]=\int d\boldsymbol{r}n^{W}(\boldsymbol{r})v_{{\text{Ext}}}$ is piece-wise linear, unlike ${\cal E}_{{\text{Ext}}}[n_{{\text{EXX}}}^{W}]$. Remember, in practice we approximate $n^{W=0,1}\approx n_{{\text{EXX}}}^{W=0,1}$.

In fact, we can go further than just removing contributions from densities. We can also use inversion to obtain a good approximation to the DD correlation energy contribution to the excitation energy, and thus shed light on the first identified issue. The key elements of this approach are illustrated in Figure 3

Consider the following relationships obeyed by the exact universal functional of our systems:

	$\displaystyle{\cal F}^{W}[n^{W}]=$	$\displaystyle{\cal F}_{{\text{EXX}}}^{W}[n^{W}]+{\cal E}_{\text{c}}^{W}[n^{W}]$		(40)
	$\displaystyle\equiv$	$\displaystyle(1-W)F_{S_{0}}[n_{S_{0}}]+WF_{P_{0}}[n_{P_{0}}]\;.$		(41)

Next, define the linearised (lin) EXX:

	$\displaystyle{\cal F}^{\text{lin},W}_{{\text{EXX}}}[n^{W}]=$	$\displaystyle(1-W){\cal F}_{{\text{EXX}}}^{0}[n^{0}]+W{\cal F}_{{\text{EXX}}}^{1}[n^{1}]$		(42)
	$\displaystyle\equiv$	$\displaystyle(1-W)F_{{\text{EXX}},S_{0}}+WF_{{\text{EXX}},P_{0}}\;,$		(43)

where $F_{{\text{EXX}},\kappa}=E_{{\text{EXX}},\kappa}-E_{{\text{Ext}},\kappa}$ for $\kappa\in S_{0},P_{0}$. We remind the reader we use $F$ rather than ${\cal F}$ to indicate that the ensembles specified by the index are treated as “pure” states. The convention can now be understood in the sense that DD correlations within (real or effective) degenerate ground-states are treated as fixed throughout the excitation. Their contributions may therefore be ignored. We introduce such an assumption for the purpose of directly focusing on the more important DD correlations in excitations.

By definition, the SD correlation terms are, $E^{{\text{SD}}}_{c,\kappa}=F^{1}[n_{\kappa}]-F_{{\text{EXX}}}[n_{\kappa}]$. It thus follows that,

	$\displaystyle{\cal E}_{\text{c}}^{{\text{SD}},W}[n^{W}]=$	$\displaystyle{\cal F}^{W}[n^{W}]-{\cal F}^{\text{lin},W}_{{{\text{EXX}}}}[n^{W}]\;,$		(44)
	$\displaystyle{\approx}$	$\displaystyle{{\cal F}^{W}[n_{{\text{EXX}}}^{W}]-{\cal F}^{\text{lin},W}_{{\text{EXX}}}[n_{{\text{EXX}}}^{W}]\;,}$		(45)

at least as far as excitations are concerned. This relationship is exact if we use the interacting density $n^{W}$ for $W=0,1$ or approximate if we use $n^{W}_{{\text{EXX}}}$. For arbitrary $W$ we can also write,

	$\displaystyle{\cal F}_{{\text{EXX}}}^{W}[n^{W}]{\approx}$	$\displaystyle(1-W)F_{{\text{EXX}},S_{0}}[n_{s,S_{0}}^{W}]$
		$\displaystyle+WF_{{\text{EXX}},P_{0}}[n_{s,P_{0}}^{W}],$		(46)

where $n^{W}=(1-W)n_{s,S_{0}}^{W}+Wn_{s,P_{0}}^{W}=(1-W)n_{S_{0}}+Wn_{P_{0}}$. Here, $n_{s,S_{0}}$ and $n_{s,P_{0}}$ are formed on the same set of orbitals from the KS potential $v_{s}[n^{W}]$ found by inversion of $n^{W}$. Here and henceforth, approximately equals signs indicate errors caused by treating the EXX densities for $W=0$ and $W=1$ as exact. For brevity we also drop the EXX subscript on densities $n^{W=0}$ and $n^{W=1}$.

Finally, using ${\cal E}_{\text{c}}={\cal E}_{\text{c}}^{{\text{SD}}}+{\cal E}_{\text{c}}^{{\text{DD}}}={\cal F}-{\cal F}_{{\text{EXX}}}$ and ${\cal E}_{\text{c}}^{{\text{SD}}}={\cal F}-{\cal F}^{\text{lin}}_{{\text{EXX}}}$, we obtain

$\displaystyle{\cal E}_{\text{c}}^{{\text{DD}},W}[n^{W}]{\approx}$

$\displaystyle{\cal F}^{\text{lin},W}_{{\text{EXX}}}[n^{W}]-{\cal F}^{W}_{{\text{EXX}}}[n^{W}]\;,$

(47)

for the DD correlation energy associated with excitations. Thus, if we calculate EXX energies at $W=0$ and $W=1$ and interpolate, this is approximately equivalent to removing SD energies at all $W$. Furthermore, deviations

	$\displaystyle\Delta{\cal T}_{s}^{W}{\approx}$	$\displaystyle{\cal T}_{s}[n^{W}]-\{(1-W){\cal T}_{s}[n^{W=0}]+W{\cal T}_{s}[n^{W=1}]\},$		(48)
	$\displaystyle\Delta{\cal E}_{\text{Hx}}^{W}{\approx}$	$\displaystyle{\cal E}_{\text{Hx}}[n^{W}]-\{(1-W){\cal E}_{\text{Hx}}[n^{W=0}]+W{\cal E}_{\text{Hx}}[n^{W=1}]\},$		(49)

have an equal and opposite component in the DD correlation terms. We are consequently able to estimate ${\cal E}_{\text{c}}^{{\text{DD}}}$ and its components without needing to know ${\cal E}_{\text{c}}^{{\text{SD}}}$ – albeit, under the assumption that $n_{{\text{EXX}}}^{W=0,1}\approx n^{W=0,1}$.

Figure 4 reveals that the density-driven correlation energy is very small in our cases, because ${\cal E}$ barely budges from piece-wise linearity, and is well within the margin of errors ($\pm 1$ kcal/mol) from density inversion (here, estimated by considering the deviation of ${\cal E}_{{\text{Ext}}}$ which is theoretically constrained to be zero). In the case of $S$–$P$ transitions within atoms it seems that density-driven correlation energies are likely to be small. This bodes well for calculations of excitations obtained from simple EDFT approximations.

Finally, we recognise that the $S$–$P$ transitions so far considered are not the only ones we have available by using ensembles composed of two “members”. The results of Section II.2 mean we can also study the lowest lying transition between $S$ states, as such ensembles are $\mathcal{P}$-passive.

Focusing on non-interacting states for illustration, let us detail the form of our EDMs. The full EDM is,

$\displaystyle\hat{\Gamma}_{s}^{W}\equiv$

$\displaystyle(1-W)\hat{\Gamma}_{s,S_{0}}+W\hat{\Gamma}_{s,S_{1}},$

(50)

where $S_{0}$ is as before and $S_{1}$ now refers to the fact that we average over first excited states with $L=0$. Specifically, for Be, we study the excitation involving promotion of the $2s$ orbital to $3s$. $\hat{\Gamma}_{s,S_{0}}$ is defined in Eq. (26). For $S_{1}$, averaging over singlet-triplet splitting, we get,

$\displaystyle\hat{\Gamma}^{\text{Be}}_{s,S_{1}}=$

$\displaystyle\frac{1}{4}\sum_{\sigma,\sigma^{\prime}\in\mathord{\uparrow},\mathord{\downarrow}}|2s^{\sigma}3s^{\sigma^{\prime}}\rangle\langle 2s^{\sigma}3s^{\sigma^{\prime}}|\;.$

(51)

For Mg we set $1s^{2}\to 1s^{2}2s^{2}$, $2s\to 3s$ and $3s\to 4s$.

Formally, the excitations we consider can only be studied for $0\leq W\leq\frac{1}{2}$. However, we found that self-consistently evaluating the fully excited state ($W=1$ – also $\frac{1}{2}<W\leq 1$) worked, despite not being formally guaranteed.

Figure 5 shows results obtained self-consistently for the $2s$–$3s$ transition in Be and the $3s$–$4s$ transition in Mg. Energies have deviations that are, compared to the total transition energy, very similar to the $S$–$P$ transition shown in Figure 1. However, the total EXX energy is noticeably more concave.

Figure 6 further highlights the concavity, and thus the importance of DD correlations. Unlike the cases shown in Figure 4, which show little deviation, the energy ${\cal E}^{W}$ deviates from piece-wise linear by up to $3$ kcal/mol for Be and $2$ kcal/mol for Mg. As discussed above, this deviation is a decent approximation to ${\cal E}_{\text{c}}^{{\text{DD}}}$, which is consequently around 2.5% of the total excitation energy of Be, and 2% of Mg. This is not a result of poor inversion either, as ${\cal E}_{{\text{Ext}}}$ shows minimal deviations from piece-wise linear.