text:
\section{Introduction}

Time-dependent density-functional theory (TD-DFT)~\cite{runge1984density} has become over the last two
decades the method of choice for modeling excited-state
properties~\cite{Casida_tddft_review_2012}.
Despite this success, it still suffers, in its standard (adiabatic)
formulation, from various limitations. The absence of
multiple-electron-excitation energies in the spectrum is one well-known
example~\cite{Casida_tddft_review_2012}. Moreover, as it relies on a
ground-state DFT calculation, linear response
TD-DFT does not provide a balanced description of
low-lying excited states. Such a description is of primary importance in
photochemistry when approching,
for example, an avoided crossing or a conical intersection, but also for
modeling the electronic structure of open $d$- or $f$-shell systems.\\

One way to overcome these limitations is to extend DFT to (canonical)
ensembles of ground and excited
states~\cite{PRA_GOK_EKSDFT,GOK3}.
Ensemble DFT relies on the Gross--Oliveira--Kohn (GOK)
variational principle~\cite{PRA_GOK_RRprinc}, which is a generalization
of Theophilou's variational principle for
equi-ensembles~\cite{JPC79_Theophilou_equi-ensembles,theophilou_book},
hence the name GOK-DFT.
Even though it is rarely mentioned, theses principles
provide a rigorous justification for the state-averaging
procedure that is routinely used in
complete active space
self-consistent field (CASSCF) calculations~\cite{mcscf_pinkbook}.
GOK-DFT has been formulated thirty years ago and, despite important conceptual
progress~\cite{Nagy_ensAC,ensemble_ghost_interaction}, it did not attract as much attention as TD-DFT until now.
Quite recently, numerous important contributions (both formal and
practical) appeared in the
literature~\cite{PRA13_Pernal_srEDFT,franck2014generalised,yang2014exact,JCP14_Burke_ensemble,pernal2015excitation,yang2017direct,gould2017hartree,gould2018charge,deur2017exact,deur2018exploring,arXiv18_Gould_ens_corr_ener_not_just_sum,sagredo2018can,senjean2018unified,JCP14_Filatov_conical_inter_REKS,Filatov-2015-Wiley,filatov2015ensemble,filatov2016self,filatov2017description}, thus making GOK-DFT an active
field of research and a promising
time-independent alternative to TD-DFT.\\

Modeling the correlation energy of an ensemble with a density functional
is a complicated task since it is not, in general,  a simple sum
of individual correlation
energies~\cite{arXiv18_Gould_ens_corr_ener_not_just_sum}. Extracting
individual energy levels is therefore not straightforward in
GOK-DFT~\cite{senjean2015linear,yang2017direct}. In the state-averaged CASSCF method the situation is different since the
(wavefunction-based) energy of each state is always computed, thus
giving access to excited-state properties (like energy gradients). From
that point of view, a state-specific DFT~\cite{PRA12_Nagy_TinD-DFT_ES,JCP15_Ayers_KS-DFT_excit-states_Coulomb,Ayers2018_Article_Time-independentDensityFunctio}
might be more appropriate.
Nevertheless, as mentioned
previously, it is often important, for example in photochemistry, to
have a balanced description (in terms of orbitals) of
ground and lower excited states. In such cases, using the ensemble formalism is
clearly relevant. Surprizingly, the flexibility of the theory regarding
the choice of the ensemble weights~\cite{PRA_GOK_RRprinc} has not been
fully explored yet. In standard GOK-DFT-based methods,
excitation energies are usually extracted from a sequence of ensemble
calculations (each of them involving a {\it single} ensemble
weight)~\cite{PRA_GOK_EKSDFT,senjean2015linear,yang2017direct}.
The Kohn--Sham DFT
limit (where all the excited-state weights become zero)
has been explored in this context, thus leading to the {\it direct ensemble
correction} (DEC) scheme of Yang
{\it et al.}~\cite{yang2017direct,sagredo2018can}. In this paper, we explore an
alternative formulation of GOK-DFT where a single {\it many}-weight-dependent ensemble exchange-correlation (xc)
functional is employed. In this formalism, all the weights can vary
independently. We show that, with such a flexibility, all individual
energy levels can be extracted, in principle exactly, from a single
GOK-DFT calculation where the ensemble weights can be freely chosen.
In contrast to TD-DFT, which gives access to excitation
energies only, this {many}-weight-dependent
formulation of GOK-DFT provides {\it total} excited-state energies.
Therefore, it should allow for a direct calculation of excited-state
properties by differentiation of the latter energies with respect to any perturbation
strength (like the nuclear displacements for the optimization of
equilibrium structures, for example).
We show that our {many}-weight-dependent approach is nothing but a generalization of DEC to
non-zero weights. As a result, it allows for a balanced description of the
states within the ensemble through the adjustment of the
weights, exactly like
in a state-averaged CASSCF calculation.\\

The paper is organized as follows. After a brief review of the GOK
principle and the various extraction procedures of individual energy levels from
an ensemble calculation (Sec.~\ref{subsec:Ens_ener_deriv_ind_ener}), we
derive in Sec.~\ref{subsec:ind_ener_gok-dft} a {\it
many}-weight-dependent version of GOK-DFT where all the energy levels
can be determined from a single calculation. The connection with
existing ensemble DFT methods is made in Sec.~\ref{subsec:connection-standard}. The theory is then applied to
the asymmetric Hubbard dimer in Sec.~\ref{sec:Hubb_dimer_theory}. The
results are discussed in Sec.~\ref{sec:results-perspec}. Comparison is
then made, at the ensemble exact exchange level of
approximation,
with the more standard extraction technique, where a
sequence of ensemble calculations is performed (see
Sec.~\ref{sec:compar_gokII}). Conclusions and perspectives are given in Sec.~\ref{conclusions}.

\section{Theory}

\subsection{Extracting individual energy levels from an ensemble
energy}\label{subsec:Ens_ener_deriv_ind_ener}

Let us consider a canonical ensemble consisting of the ground and $M$
first
excited states of the electronic Hamiltonian $\hat{H}=\hat{T}
+\hat{W}_{\rm ee}+\hat{V}_{\rm ext}$. The operators $\hat{T}$ and $\hat{W}_{\rm
ee}$ describe the electronic kinetic and repulsion
energies, respectively. The local external potential operator reads $\hat{V}_{\rm ext}=\int{d}\mathbf{r}\,v_{\rm
ext}(\mathbf{r})\hat{n}(\mathbf{r})$ where $\hat{n}(\mathbf{r})$ is the
density operator and $v_{\rm
ext}(\mathbf{r})$ will simply be the nuclear Coulomb potential in this
work.
For the sake of clarity, we will assume in the following that none of
these states are degenerate. The formalism can be easily extended to
degenerate ensembles by assigning the same weight to degenerate
states~\cite{PRA_GOK_EKSDFT,yang2017direct}. In the most general
formulation of the GOK variational principle~\cite{PRA_GOK_RRprinc}, the exact ensemble energy reads
\begin{eqnarray}\label{eq:ens_ener}
E^{\mathbf{w}}=\left(1-\sum^M_{I=1}{\tt w}_I\right)E_0+\sum^M_{I=1}{\tt
w}_IE_I,
\end{eqnarray}
where $E_0$ is the ground-state energy, $\left\{E_I\right\}_{1\leq I\leq
M}$ are
the $M$ first excited-state energies, and $\mathbf{w}\equiv ({\tt
w}_1,{\tt
w}_2,\ldots,{\tt
w}_M)$ denotes the collection of weights that are assigned to each
individual {\it excited} state. In their seminal
paper~\cite{PRA_GOK_EKSDFT}, Gross {{\it et al.}} considered a
{\it sequence} of ensemble DFT calculations in order to extract excitation
energies. In their approach, each (non-degenerate here) ensemble is a linear interpolation (controlled by a {\it single}
ensemble weight $w$) between equi-ensembles:
\begin{eqnarray}\label{eq:gokI_weights}
{\tt w}_{1\leq I<M}=\dfrac{1-w}{M},\hspace{0.2cm}{\tt w}_M=w.
\end{eqnarray}
More recently, Yang {\it et al.}~\cite{yang2017direct} used another set
of ensembles (the approach was referred
to as GOKII) which are
also characterized by a {single} weight $w$:
\begin{eqnarray}\label{eq:gokII_weights}
{\tt w}_{1\leq I\leq M}=w.
\end{eqnarray}
The practical advantage of Eq.~(\ref{eq:gokII_weights}) over
Eq.~(\ref{eq:gokI_weights}) is that two ensemble calculations are
sufficient for extracting any excitation energy~\cite{yang2017direct}.
In Ref.~\cite{yang2017direct}, the authors implemented
Eq.~(\ref{eq:gokII_weights}) in the $w\rightarrow 0$ limit, thus
providing a {\it direct ensemble
correction} (DEC) to Kohn--Sham (KS) excitation energies.\\

One practical drawback of both DEC and linear response TD-DFT
is that, in contrast to state-averaged CASSCF~\cite{mcscf_pinkbook}, it
is not straightforward to study, within their formalisms, the potential
energy
curve of one or more excited states, simply because a sequence of different
calculations is
needed. Moreover (and perhaps, more importantly) none of them provides a balanced description (in
terms of orbitals) of the
ground and lower excited states. This can become problematic, for example, in the
vicinity of a conical intersection.\\

In order to address these deficiencies, we explore in this paper a more general formulation of GOK-DFT
where the ensemble weights can all vary {\it independently}. Note that the ensemble energy can be obtained
variationally if the weights decrease with increasing
index~\cite{PRA_GOK_RRprinc}, i.e. if, for
$1\leq J\leq (M-1)$,
\begin{eqnarray}\label{eq:convex_weights_1}
{\tt w}_J\geq {\tt w}_{J+1}\geq0,
\end{eqnarray}
and
\begin{eqnarray}\label{eq:convex_weights_0}
\left(1-\sum^M_{I=1}{\tt w}_I\right)\geq{\tt w}_1.
\end{eqnarray}
Before introducing our alternative extraction procedure, we would like
to stress that, unlike state-averaged wavefunction-based methods,
GOK-DFT gives a direct access to the ensemble energy
$E^{\mathbf{w}}$ only, and not to its individual-state components (i.e. the
energy levels). The reason
is that, in GOK-DFT, a single density functional is used for describing the
xc energy of the ensemble. In the latter are mixed, in a non-trivial
way, the
individual correlation energies of all the states that belong to the
ensemble~\cite{arXiv18_Gould_ens_corr_ener_not_just_sum}.\\

Even though
excitation (or individual) energies cannot be extracted from a single
ensemble energy value $E^{\mathbf{w}}$, infinitesimal variations in the ensemble
weights will immediately give access to its individual components.
Indeed, starting from the fact that the derivative of the ensemble
energy with respect to ${\tt w}_I$ is equal to the $I$th excitation energy,
\begin{eqnarray}\label{eq:Xener}
\dfrac{\partial E^{\mathbf{w}}}{\partial {\tt w}_I}=E_I-E_0,
\end{eqnarray}
and keeping in mind that the ensemble energy varies {\it linearly} with the
ensemble weights (see Eq.~(\ref{eq:ens_ener})),
\begin{eqnarray}
E^{\mathbf{w}}=E_0+\sum^M_{I=1}{\tt w}_I\dfrac{\partial E^{\mathbf{w}}}{\partial
{\tt w}_I},
\end{eqnarray}
or, equivalently,
\begin{eqnarray}\label{eq:E0_wrt_ens_ener}
E_0=
E^{\mathbf{w}}
-\sum^M_{I=1}{\tt w}_I\dfrac{\partial E^{\mathbf{w}}}{\partial
{\tt w}_I},
\end{eqnarray}
we can rewrite any individual (ground- or excited-state) energy
as
\begin{eqnarray}\label{eq:ind_ener_exp_deriv_wft}
E_K&=&E_0+\sum^M_{I=1}\delta_{IK}(E_I-E_0)
\nonumber\\
&=&E^{\mathbf{w}}+\sum^M_{I=1}\left(\delta_{IK}-{\tt
w}_I\right)\dfrac{\partial E^{\mathbf{w}}}{\partial {\tt w}_I},
\end{eqnarray}
where $0\leq K\leq M$.
The derivation of
Eq.~(\ref{eq:ind_ener_exp_deriv_wft}) is trivial. Nevertheless, to the best
of our knowledge, it has never been
used in the context of GOK-DFT. As shown in the following, the
expression in Eq.~(\ref{eq:ind_ener_exp_deriv_wft}) is convenient for connecting the
exact individual energy levels to the KS orbital
energies. Most importantly, it will
enable us to show that a single GOK-DFT calculation (where the
weights can be freely chosen) is in principle sufficient for extracting
all the energy
levels.

\subsection{Density-functional theory for
ensembles}\label{subsec:ind_ener_gok-dft}

In GOK-DFT, the ensemble energy is determined variationally as
follows~\cite{PRA_GOK_EKSDFT},
\begin{eqnarray}\label{eq:min_gamma_Ew}
 E^{\mathbf{w}} &=& \min_{{\hat{\gamma}}^{\mathbf{w}}}
\Big\{
{\rm Tr}\left[\hat{\gamma}^{\mathbf{w}}
\left(\hat{T}
+\hat{V}_{\rm ext}\right)
\right]+ E_{\rm
Hxc}^{\mathbf{w}}[n_{\hat{\gamma}^{\mathbf{w}}}]
\Big\}
\nonumber\\
&=&
{\rm Tr}\left[\hat{\gamma}_{\rm s}^{\mathbf{w}}
\left(\hat{T}
+\hat{V}_{\rm ext}\right)
\right]+ E_{\rm Hxc}^{\mathbf{w}}[n_{\hat{\gamma}_{\rm s}^{\mathbf{w}}}]
,
\end{eqnarray}
where $n_{\hat{\gamma}^{\mathbf{w}}}(\mathbf{r})={\rm Tr}\left[\hat{\gamma}^{\mathbf{w}}
\hat{n}(\mathbf{r})
\right]$ is a trial ensemble density and
\begin{eqnarray}\label{eq:Hxc_ens_fun}
E^{\mathbf{w}}_{\rm Hxc}[n]
=
 \frac{1}{2} \iint
{d}{\bf r} {d}{\bf r'}
~ \dfrac{n({\bf r})n({\bf r'})}{\mid{{\bf r}-{\bf r'}}\mid}
+
E^{\mathbf{w}}_{\rm xc}[n]
\end{eqnarray}
is the ensemble Hartree xc (Hxc) functional. We use here the original
in-principle-exact decomposition of the Hxc functional~\cite{PRA_GOK_EKSDFT} where, for a given and fixed
density $n$, the xc part only varies with $\mathbf{w}$. In practical
(approximate) calculations, it might be worth using another
decomposition~\cite{arXiv18_Gould_ens_corr_ener_not_just_sum} which is
ghost-interaction-free~\cite{ensemble_ghost_interaction}. In this work, we will
always use exact Hxc (or Hx) functionals. Returning to Eq.~(\ref{eq:min_gamma_Ew}), the ground and excited KS determinants in the minimizing
non-interacting density matrix operator $\hat{\gamma}_{\rm
s}^{\mathbf{w}}=\left(1-\sum^M_{I=1}{\tt w}_I\right)
\left\vert\Phi^{\mathbf{w}}_0\right\rangle\left\langle\Phi^{\mathbf{w}}_0\right\vert
+
\sum^M_{I=1}{\tt
w}_I\left\vert\Phi^{\mathbf{w}}_I\right\rangle\left\langle\Phi^{\mathbf{w}}_I\right\vert
$
are determined by solving the ensemble KS equations self-consistently,
\begin{eqnarray}\label{eq:sc_ks-edft}
\left(\hat{T}+\int {d}\mathbf{r}\,v_{\rm
s}^{\mathbf{w}}(\mathbf{r})\hat{n}(\mathbf{r})\right)\left\vert\Phi^{\mathbf{w}}_K\right\rangle
=\mathcal{E}^{\mathbf{w}}_K\left\vert\Phi^{\mathbf{w}}_K\right\rangle,
\end{eqnarray}
where the ensemble KS potential reads $v_{\rm
s}^{\mathbf{w}}(\mathbf{r})=v_{\rm
ext}(\mathbf{r})+\delta E_{\rm Hxc}^{\mathbf{w}}[n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}]/\delta n(\mathbf{r})$ and $0\leq K\leq M$. Note that
the (weight-dependent) KS energy $\mathcal{E}^{\mathbf{w}}_K$ is simply obtained by
summing up the energies of the spin-orbitals that are occupied in
$\Phi^{\mathbf{w}}_K$.\\

From the GOK-DFT ensemble energy expression in
Eq.~(\ref{eq:min_gamma_Ew}) and the expression for the individual
energies in Eq.~(\ref{eq:ind_ener_exp_deriv_wft}), we can now derive
exact density-functional expressions for all the energy levels included
into the ensemble. Indeed, according to the Hellmann--Feynman theorem
and Eq.~(\ref{eq:Hxc_ens_fun}),
we can first express the ensemble energy derivative as follows,
\begin{eqnarray}
\dfrac{\partial E^{\mathbf{w}}}{\partial
{\tt w}_I}&=&{\rm Tr}\left[{\Delta}\hat{\gamma}_{{\rm s},I}^{\mathbf{w}}
\left(\hat{T}
+\hat{V}_{\rm ext}\right)
\right]
+
\left.\dfrac{\partial E_{\rm xc}^{\mathbf{w}}[n]}{\partial
{\tt w}_I}\right|_{n=n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}}
\nonumber\\
&&+\int {d}\mathbf{r}\,\dfrac{\delta E_{\rm Hxc}^{\mathbf{w}}[n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}]}{\delta n(\mathbf{r})}{\rm
Tr}\left[{{\Delta}\hat{\gamma}_{{\rm
s},I}^{\mathbf{w}}}\hat{n}(\mathbf{r})\right],
\end{eqnarray}
where ${\Delta}\hat{\gamma}_{{\rm
s},I}^{\mathbf{w}}=\left\vert\Phi^{\mathbf{w}}_I\right\rangle\left\langle\Phi^{\mathbf{w}}_I\right\vert-\left\vert\Phi^{\mathbf{w}}_0\right\rangle\left\langle\Phi^{\mathbf{w}}_0\right\vert$,
thus leading to the following exact expression for the $I$th excitation
energy (see Eq.~(\ref{eq:sc_ks-edft})),
\begin{eqnarray}\label{eq:final_exp_dEw_over_dwI}
\dfrac{\partial E^{\mathbf{w}}}{\partial
{\tt w}_I}=E_I-E_0=\mathcal{E}^{\mathbf{w}}_I-\mathcal{E}^{\mathbf{w}}_0+
\left.\dfrac{\partial E_{\rm xc}^{\mathbf{w}}[n]}{\partial
{\tt w}_I}\right|_{n=n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}},
\end{eqnarray}
which generalizes the GOK-DFT expression for the optical gap~\cite{PRA_GOK_EKSDFT} to
higher excitations. Note that, in the original formulation of GOK-DFT~\cite{PRA_GOK_EKSDFT},
higher excitation energies were obtained from a sequence of
{\it single}-weight-dependent ensemble calculations instead. This is not
necessary anymore here as we use a {\it many}-weight-dependent xc
functional.\\

For formal convenience, we now propose to extend the Levy--Zahariev (LZ)
shift-in-potential procedure~\cite{levy2014ground} to canonical ensembles, in
complete analogy with
Ref.~\cite{senjean2018unified},
\begin{eqnarray}\label{eq:LZ_shift}
\dfrac{\delta E_{\rm
Hxc}^{\mathbf{w}}[n]}{\delta n({\bf r})}
&\rightarrow&\overline{v}^{\mathbf{w}}_{\rm Hxc}[n]({\bf r})=
\dfrac{\delta E_{\rm
Hxc}^{\mathbf{w}}[n]}{\delta n({\bf r})}
\nonumber\\
&&
+
\dfrac{E^{\mathbf{w}}_{\rm Hxc}[n]-\int{d}{\bf r}\;
\dfrac{\delta E_{\rm
Hxc}^{\mathbf{w}}[n]}{\delta n({\bf r})}
n({\bf r})}{\int{d}{\bf r}\;n({\bf r})}
.
\end{eqnarray}
Thus we obtain the following shifted KS energy expressions,
\begin{eqnarray}\label{eq:ens_LZ_shifted_KS_ener}
\mathcal{E}^{\mathbf{w}}_K\rightarrow
\overline{\mathcal{E}}^{\mathbf{w}}_K
&=&\mathcal{E}^{\mathbf{w}}_K+
E_{\rm Hxc}^{\mathbf{w}}[n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}]
\nonumber\\
&&-
\int{d}{\bf r}\;
\dfrac{\delta E_{\rm
Hxc}^{\mathbf{w}}[n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}]}{\delta n({\bf r})}
n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}({\bf r}).
\end{eqnarray}
As a result [see Eqs.~(\ref{eq:min_gamma_Ew}) and
(\ref{eq:sc_ks-edft})], the exact ensemble energy can be written as a
weighted sum of shifted KS energies,
\begin{eqnarray}\label{eq:sum_shifted_KS_ener}
{E}^{\mathbf{w}}=\left(1-\sum^M_{I=1}{\tt
w}_I\right)\overline{\mathcal{E}}^{\mathbf{w}}_0+\sum^M_{I=1}{\tt
w}_I\overline{\mathcal{E}}^{\mathbf{w}}_I.
\end{eqnarray}
Let us stress that, as readily seen from
Eq.~(\ref{eq:sum_shifted_KS_ener}), the LZ shifting procedure is a way to truly fix (i.e. not anymore up to a
constant) the KS (orbital) energies and, consequently, the ensemble KS
potential. Indeed, as shown in Eq.~(\ref{eq:LZ_shift}), any constant
added to the ensemble Hxc potential will be automatically removed by the
LZ shift. Note also that, by construction, the ensemble Hxc
density-functional energy reads
\begin{eqnarray}
E^{\mathbf{w}}_{\rm Hxc}[n]=\int{d}{\bf
r}\;\overline{v}^{\mathbf{w}}_{\rm Hxc}[n]({\bf r})n({\bf r}).
\end{eqnarray}
\\
As a result, we could think of modeling the shifted Hxc
ensemble potential $\overline{v}^{\mathbf{w}}_{\rm Hxc}[n]({\bf r})$
directly rather than the Hxc ensemble energy, in complete
analogy with Ref.~\cite{levy2014ground}. This is where, in this context, the LZ shift
becomes (much) more than a convenient formal trick. This path will not
be explored further in the rest of the paper and is left for future
work.\\

Turning finally to the extraction of individual energies,
we should keep in mind that the (global) LZ shift
does not affect KS energy differences,
\begin{eqnarray}\label{eq:LZ_shift_noeffect_XE}
\overline{\mathcal{E}}^{\mathbf{w}}_I-\overline{\mathcal{E}}^{\mathbf{w}}_0=\mathcal{E}^{\mathbf{w}}_I-\mathcal{E}^{\mathbf{w}}_0,
\end{eqnarray}
and therefore, as readily seen from
Eq.~(\ref{eq:final_exp_dEw_over_dwI}), it leaves the true excitation
energies unchanged. It only plays a role in the calculation of exact
energy levels. Indeed, if we combine
Eq.~(\ref{eq:LZ_shift_noeffect_XE}) with Eqs.~(\ref{eq:ind_ener_exp_deriv_wft}),
(\ref{eq:final_exp_dEw_over_dwI}), and (\ref{eq:sum_shifted_KS_ener}),
we obtain the following compact expressions,
\begin{eqnarray}\label{eq:ind_ener_with_DD}
E_K=\overline{\mathcal{E}}^{\mathbf{w}}_K+\sum^M_{I=1}\left(\delta_{IK}-{\tt
w}_I\right)\left.\dfrac{\partial E_{\rm xc}^{\mathbf{w}}[n]}{\partial
{\tt w}_I}\right|_{n=n_{\hat{\gamma}_{\rm
s}^{\mathbf{w}}}}.
\end{eqnarray}
Once the ensemble xc derivative corrections (second term on the right-hand
side of Eq.~(\ref{eq:ind_ener_with_DD})) have been added to the unshifted KS
energies, applying the LZ shift gives immediately access to {\it any} energy level
in the ensemble, and therefore to any ground- or excited-state molecular property. Unlike in the standard DFT+TD-DFT procedure, a
single calculation is in principle sufficient.\\

Let us
stress that Eq.~(\ref{eq:ind_ener_with_DD}), which is the key result of
this paper, holds for {\it any} set of ordered ensemble
weights [see Eqs.~(\ref{eq:convex_weights_1}) and
(\ref{eq:convex_weights_0})], including both ground-state ${\tt
w}_{1\leq I\leq M}\rightarrow 0$ and
equi-ensemble ${\tt w}_{1\leq I\leq M}\rightarrow 1/(M+1)$ limits. In this
respect, it
generalizes the original formulation of GOK-DFT~\cite{PRA_GOK_EKSDFT} (where
single-weight-dependent xc functionals only were introduced) as well as the
more recent DEC method~\cite{yang2017direct} which, as shown in the following, is recovered
from
the ground-state limit of Eq.~(\ref{eq:ind_ener_with_DD}). Note finally
that the latter equation extends the recent work of Senjean and
Fromager on charged excitations~\cite{senjean2018unified}
to neutral excitation processes.\\

\subsection{Connection with existing ensemble
DFT approaches}\label{subsec:connection-standard}

We should point out that our formalism may be connected to the very
recent work of Gould
and Pittalis~\cite{arXiv18_Gould_ens_corr_ener_not_just_sum} on the
expression of density-functional ensemble xc energies in terms of
individual-state contributions. Indeed, starting from
Eq.~(\ref{eq:ind_ener_with_DD}), we could derive, for each state, an individual xc functional that is a bi-functional
of the individual KS density (through the unshifted KS energy) and the
ensemble one. By taking the weighted sum of these bi-functionals we recover a decomposition for the
ensemble xc energy which resembles the one of Gould
and Pittalis~\cite{arXiv18_Gould_ens_corr_ener_not_just_sum}. The connection between the two approaches should
clearly be
explored further. This is left for future work.\\

We also note from Eq.~(\ref{eq:ind_ener_with_DD}) that, even though both
terms on the right-hand side are in principle
weight-dependent, their sum should of course be
weight-independent. As shown in the following, this will not be the case
anymore when approximate
xc density functionals are used. Note also that,
in the $\mathbf{w}=
0$ limit, which has been used in previous
works~\cite{PRA_GOK_EKSDFT,PRA_Levy_XE-N-N-1,yang2017direct,sagredo2018can},
the LZ ground-state energy expression
$E_0=\overline{\mathcal{E}}^{\mathbf{w}=0}_0$~\cite{levy2014ground} is recovered and, most importantly, the excited-state
energy expressions can be simplified further as follows,
\begin{eqnarray}\label{eq:w_zero_limit}
E_J=\overline{\mathcal{E}}^{\mathbf{w}=0}_J+\left.\dfrac{\partial E_{\rm
xc}^{\mathbf{w}}[n_0]}{\partial
{\tt w}_J}\right|_{{\mathbf{w}}=0}
\end{eqnarray}
where $n_0$ denotes the ground-state density.
As shown in the seminal work of Levy~\cite{PRA_Levy_XE-N-N-1} and readily seen from Eq.~(\ref{eq:w_zero_limit}),
both ground and $J$th excited states cannot be described with the same
KS potential. The latter should indeed exhibit a jump [see the second
term on the right-hand side of Eq.~(\ref{eq:w_zero_limit})], which is
known as the derivative
discontinuity (DD), as the
(neutral) excitation process occurs, exactly like in charged excitation
processes~\cite{senjean2018unified}. If we are able to
model the {\it many}-weight-dependence of the ensemble xc functional,
then we have access to {\it all} ensemble xc derivatives $\partial E_{\rm
xc}^{\mathbf{w}}[n]/\partial{\tt w}_I$ and therefore, by considering
the ground-state $\mathbf{w}\rightarrow 0$ limit, we obtain {\it all}
the DDs.\\

Note finally that, if
we use Eq.~(\ref{eq:w_zero_limit}) to compute the $J$th excitation energy, we
recover the bare KS excitation energy (i.e. the sum of KS orbital energy
differences) to which an ensemble xc derivative correction is applied. When rewritten as follows,
\begin{eqnarray}\label{eq:DEC_recovered_w=0}
\left.\dfrac{\partial E_{\rm
xc}^{\mathbf{w}}[n_0]}{\partial
{\tt w}_J}\right|_{{\mathbf{w}}=0}=
\left[\dfrac{{d}
E^{\mathbf{w}=\underline{w}^J}_{\rm
xc}[n_0]}{{d} w}-
\dfrac{{d}
E^{\mathbf{w}=\underline{w}^{J-1}}_{\rm
xc}[n_0]}{{d} w}
\right]_{w=0},
\end{eqnarray}
where $\underline{w}^J$ is the ensemble weight vector defined by ${\tt
w}_I=w$ for $1\leq I\leq J$ and ${\tt
w}_I=0$ for $J < I\leq M$, it becomes clear that, in the ground-state
limit, our approach reduces to the DEC one~\cite{yang2017direct}. By
considering a {\it many}-weight-dependent xc functional, we simply extend the applicability of DEC to any kind of ensemble (including
equi-ensembles). We also obtain all the energies from a single ensemble calculation.
\section{Application to the Hubbard dimer}\label{sec:Hubb_dimer_theory}

We present in the following an implementation of
Eq.~(\ref{eq:ind_ener_with_DD}) for a three-state singlet ensemble. In
the latter case, the convexity conditions in Eqs.~(\ref{eq:convex_weights_1}) and
(\ref{eq:convex_weights_0}) become
\begin{eqnarray}
0\leq {\tt w}_2\leq {1}/{3}
\end{eqnarray}
and
\begin{eqnarray}
{\tt
w}_2\leq{\tt w}_1\leq (1-{\tt
w}_2)/2.
\end{eqnarray}
The theory is applied to the (not necessarily symmetric) Hubbard
dimer~\cite{carrascal2015hubbard,carrascal2016corrigendum}. It is a simple but non-trivial toy system that is nowadays routinely
used for
exploring new concepts in
DFT~\cite{carrascal2015hubbard,carrascal2016corrigendum,li2018density,carrascal2018linear,sagredo2018can,PRB18_Carsten_dft_non-coll_spin_Hubbard_dimer,senjean2018unified,deur2017exact,deur2018exploring}.
Within this model, the Hamiltonian is simplified as follows (we write operators in second
quantization):
\begin{eqnarray}\label{eq:Hamil_Hubbard_dimer_model}
\hat{T} &\rightarrow& -t
\sum_{\sigma=\uparrow\downarrow}(\hat{c}^\dagger_{0\sigma}\hat{c}_{1\sigma} +
\hat{c}^\dagger_{1\sigma}\hat{c}_{0\sigma}),\hspace{0.2cm} \hat{W}_{\rm ee}\rightarrow
U\sum^1_{i=0}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},
\nonumber\\
\hat{V}_{\rm ext}&\rightarrow&\Delta v_{\rm ext}(\hat{n}_1 -
\hat{n}_0)/2,\hspace{0.3cm}\hat{n}_{i\sigma}=\hat{c}^\dagger_{i\sigma}\hat{c}_{i\sigma},
\end{eqnarray}
where $\hat{n}_i=\sum_{\sigma=\uparrow\downarrow}\hat{n}_{i\sigma}$ is
the density operator on site $i$ ($i=0,1$). Note that the external potential reduces to a      single number $\Delta v_{\rm
ext}$ which controls the asymmetry of the model.
The density also reduces
to a single number $n=n_0$ which is the occupation of site 0, given that
$n_1 = 2-n$ (we consider 2-electron canonical ensembles only in this
work).\\

The bi-ensemble consisting of the ground and first singlet excited
states has been extensively studied in
Refs.~\cite{deur2017exact,deur2018exploring}. Very recently, Sagredo and
Burke~\cite{sagredo2018can} added one more (doubly excited) singlet state
to the ensemble. As proven in Appendix~\ref{appendix:reduction_tri-to_bi}, the
tri-ensemble analog of the Hohenberg--Kohn functional can be expressed in
terms of the bi-ensemble one.
As a result, both the ensemble non-interacting kinetic energy $T_{\rm
s}^{\mathbf{w}}(n)$ and the ensemble exact exchange (EEXX) one
$E^{\mathbf{w}}_{\rm x}(n)$ [here $\mathbf{w}\equiv
({\tt w}_1,{\tt w}_2)$] can be determined from their bi-ensemble analogs (see Eqs.~(57) and (62) in Ref.~\cite{deur2017exact}), thus leading to the
simple expressions
\begin{eqnarray}\label{eq:Ts_ens_3states}
T_{\rm s}^{\mathbf{w}}(n)=-2t\sqrt{(1-{\tt
w}_1-2{\tt w}_2)^2-(1-n)^2},
\end{eqnarray}
and
\begin{eqnarray}\label{eq:final_expression_3-state_EEXX}
E^{\mathbf{w}}_{\rm x}(n)
&=&
\dfrac{U}{2}\left[1+{\tt w}_1-\dfrac{(3{\tt w}_1-1)(1-n)^2}{(1-{\tt
w}_1-2{\tt w}_2)^2}\right]
\nonumber
\\
&&-E_{\rm
H}(n),
\end{eqnarray}
where the Hartree energy reads $E_{\rm H}(n)=
U\left(1 + (n - 1)^2\right)$~\cite{deur2017exact}.
The tri-ensemble density-functional correlation energy is then obtained as follows (
see Appendix~\ref{appendix:reduction_tri-to_bi}),
\begin{eqnarray}\label{eq:scaling_relation_tri-ens-corr}
E^{\mathbf{w}}_{\rm c}(n)&=&(1-3{\tt w}_2)E_{\rm c}^{w
}(\nu),
\end{eqnarray}
where $E_{\rm c}^{w}(\nu)$ a bi-ensemble correlation energy, with {\it
effective} weight
$w=({{\tt w}_1-{\tt
w}_2})/({1-3{\tt w}_2})$ and density $\nu=(n-3{\tt w}_2)/(1-3{\tt
w}_2)$, which can be computed to arbitrary accuracy by Lieb
maximization~\cite{deur2017exact}. The tri- to bi-ensemble
reduction in Eq.~(\ref{eq:scaling_relation_tri-ens-corr}) is of course
{\it not} a general result. It only applies to the Hubbard dimer and originates
from the fact that, in this system, the three singlet energies
sum up to $2U$ (see Eq.~\ref{eq:simplification_tri-to_bi}).\\

Turning to the non-interacting KS system
with potential $\Delta v_{\rm s}^{{\mathbf{w}}}$, the
(unshifted) energies of the ground-, singly- and doubly-excited states
read $\mathcal{E}^{\mathbf{w}}_0=2\varepsilon_{\rm H}\left(\Delta
v_{\rm s}^{{\mathbf{w}}}\right)$, $\mathcal{E}^{\mathbf{w}}_1=0$, and
$\mathcal{E}^{\mathbf{w}}_2=-2\varepsilon_{\rm H}\left(\Delta
v_{\rm s}^{{\mathbf{w}}}\right)$, respectively, where $\varepsilon_{\rm H}(\Delta v)=
- \sqrt{t^2 + (\Delta v^2 / 4)}$~\cite{deur2017exact}. Note that the density-functional KS potential
can be simply calculated as $\Delta v_{\rm
s}^{{\mathbf{w}}}(n)=\partial T_{\rm
s}^{\mathbf{w}}(n)/\partial n$~\cite{deur2017exact}. The Hxc potential, which is needed
in the LZ shift-in-potential procedure (see
Eq.~(\ref{eq:LZ_shift})), is then determined as follows, $\Delta v_{\rm
Hxc}^{{\mathbf{w}}}(n)=\Delta v_{\rm s}^{{\mathbf{w}}}(n)-\Delta v_{\rm
ext}$, where $n$ is
the physical tri-ensemble density obtained from the Hamiltonian in
Eq.~(\ref{eq:Hamil_Hubbard_dimer_model}). As shown in Appendix~\ref{appendix:symmetric_dimer}, in the symmetric
case ($\Delta v_{\rm
ext}=0$), the full problem can be solved analytically.\\

\section{Results and discussion}\label{sec:results-perspec}

We have shown in Sec.~\ref{subsec:ind_ener_gok-dft} that individual energy levels can be
extracted, in principle exactly, from a single many-weight-dependent ensemble GOK-DFT calculation by
adding to each (ground- and excited-state) KS energy
a global LZ-type shift {\it and} an ensemble-based state-specific xc
derivative
correction (see Eqs.~(\ref{eq:ens_LZ_shifted_KS_ener}) and (\ref{eq:ind_ener_with_DD})). In order to assess
the importance of both corrections, we first investigate the deviation
of the KS energies from the exact physical ones. The former are
simply obtained
by summing up (unshifted) KS orbital energies. Note that, in
contrast to the LZ-shifted ones, these energies are not uniquely
defined because the KS potential is unique {up to an arbitrary
constant}. In the Hubbard dimer model, the latter is chosen such that the
potential sums to zero over the two sites (see
Eq.~(\ref{eq:Hamil_Hubbard_dimer_model})).
\begin{figure}
\resizebox{0.45\textwidth}{!}{
\includegraphics[scale=1.0]{newfig1_final.pdf}
}
\caption{
Unshifted $\mathcal{E}^{\mathbf{w}}_{K=0,1,2}$ and LZ-shifted
$\overline{\mathcal{E}}^{\mathbf{w}}_{K=0,1,2}$ KS energies obtained for
the asymmetric Hubbard dimer by
varying the first ensemble weight while fixing the second one to
zero. Results are shown for $U/t=2$ (top panel) and $U/t=10$ (bottom
panel). Comparison is made with the
exact energies $E_{K=0,1,2}$. First and second excited-state energies
are shown in green and blue, respectively. EEXX-only results (including both LZ shift
and ensemble xc derivative corrections) are also plotted in the top panel (with squares) for
analysis
purposes.
}
\label{fig:unshift-LZ-shift_ener_asym_U1U5}
\end{figure}
\begin{figure}
\centering
\resizebox{0.45\textwidth}{!}{
\includegraphics[scale=1.0]{newfig2_final.pdf}
}
\caption{
Exact ensemble xc derivatives obtained from Eqs.~(\ref{eq:Xener}) and (\ref{eq:final_exp_dEw_over_dwI}) for the asymmetric ($\Delta
v_{\rm ext}/t=2$) Hubbard dimer
with $U/t=2$. First- ($I=1$) and
second-excited-state ($I=2$) derivatives are calculated for the true
(weight-dependent) ensemble density generated from $\Delta v_{\rm
ext}$ and are plotted (in green and blue,
respectively) as functions of the first ensemble weight ${\tt
w}_1$ for ${\tt
w}_2=0$ (top panel) and ${\tt
w}_2=1/4$ (bottom panel).}
\label{fig:each_exact_DDs_U1_deltav1}
\end{figure}
As illustrated in
Fig.~\ref{fig:unshift-LZ-shift_ener_asym_U1U5}, the unshifted KS energies
are found to be substantially lower than the exact energies. It is particularly
striking for the first excited state whose unshifted KS energy
equals zero, by construction (see Sec.~\ref{sec:Hubb_dimer_theory}). In the
symmetric case (see the Appendix~\ref{appendix:symmetric_dimer} and the supplementary material), the
ground- and second-excited-state unshifted KS energies are equal to $-2t$ and $2t$, respectively.
As a result, varying the ensemble weights has no impact. The situation
is different in the asymmetric case since the unshifted energies can vary with
the weights through the density-functional KS potential. The second (doubly-) excited-state
energy can for example be substantially improved when increasing the
weights. However, the ground-state energy
deteriorates in that case.\\

If we now apply the (weight-dependent) LZ shift, more accurate energies
are obtained, as shown in Fig.~\ref{fig:unshift-LZ-shift_ener_asym_U1U5}. Note that, by
construction, the LZ-shifted KS ground-state energy is exact
when ${\tt w}_1={\tt w}_2=0$. It is important to notice that,
unlike the exact energies, the LZ-shifted ones are (sometimes strongly)
weight-dependent, thus illustrating the importance of modeling
ensemble xc derivative corrections. The latter are plotted in
Fig.~\ref{fig:each_exact_DDs_U1_deltav1}.
Interestingly, both first- and
second-excited-state derivatives are non-negligible and will therefore contribute to the exact
ground-state energy away from the $\mathbf{w}=0$ limit (see the
second term on the right-hand side of Eq.~(\ref{eq:ind_ener_with_DD})).
Note also that these derivatives are strongly state-dependent. In the
asymmetry and correlation regimes considered in
Fig.~\ref{fig:each_exact_DDs_U1_deltav1}, each derivative vanish for
particular (state-dependent) weight values. In this
case, the corresponding excitation energy is exactly equal to the KS
one (see Eq.~(\ref{eq:final_exp_dEw_over_dwI})).\\

Returning to the LZ-shifted KS energies, their weight dependence becomes even
more important in stronger correlation regimes, as shown in the bottom
panel of
Fig.~\ref{fig:unshift-LZ-shift_ener_asym_U1U5}.
Note that, in this case, the first and
second excited
states are single- and double-charge transfer states, respectively~\cite{deur2017exact}.
Note also that the first-weight-dependence of the shifted energies is sensitive to
the value of the second weight, as shown in the
supplementary material.
Interestingly, increasing the ensemble weights can
provide more accurate excited-state LZ-shifted energies, often at the
expense of deteriorated ground-state energies. As shown in
Fig.~\ref{fig:unshift-LZ-shift_ener_asym_U1U5} (see also the
supplementary material), this is a general trend that can be seen in all correlation regimes.\\

Finally, in order to assess the importance of correlation effects in the
calculation of individual energy levels,
we computed EEXX-only
LZ shift and ensemble derivative corrections to the unshifted KS energies. Since we used
exact densities (and therefore exact KS potentials), the LZ shift has
been computed with the full
(exact) Hxc
potential in conjunction with the EEXX energy, for the sake of
consistency. In the moderately correlated $U/t=2$ regime (see
the top panel of Fig.~\ref{fig:unshift-LZ-shift_ener_asym_U1U5}),
relatively good total energies
are obtained, which is in agreement with the DEC/EEXX results of Ref.~\cite{sagredo2018can}.
Interestingly, the doubly-excited state energy is the one that exhibits the weakest
weight dependence. As shown in Fig.~\ref{fig:EEXX_U5}, in the $\Delta v_{\rm
ext}/t=2$ asymmetry regime, EEXX fails
dramatically for the larger $U/t=10$ value. Total energies become strongly weight-dependent and
their ordering is wrong for a wide range of weight values. The latter
observation was actually expected for small weight values on the basis
of
Ref.~\cite{deur2017exact} (where we see in Fig. 1 that, for $2t=1$, $U=5$ and
$\Delta v_{\rm ext}=1$, the ground-state density is close to 1, which
corresponds to the symmetric case) and
Appendix~\ref{appendix:symmetric_dimer}, where the EEXX
energies are derived for the symmetric Hubbard dimer (see
Eq.~(\ref{eq:symmetric_EEXX_ener})).   \\
\begin{figure}
\centering
\resizebox{0.45\textwidth}{!}{
\includegraphics[scale=1.0]{newfig3_final.pdf}
}
\caption{Ground- (red), first- (green) and second- (blue) excited-state
energies (exact and EEXX-only) plotted as functions of the first
ensemble weight (with ${\tt w}_2=0$) for the asymmetric ($\Delta
v_{\rm ext}/t=2$) Hubbard
dimer with $U/t=10$.}
\label{fig:EEXX_U5}
\end{figure}

\section{Single versus sequence of ensemble calculations
}\label{sec:compar_gokII}

While, in conventional GOK-DFT approaches, excitation (or individual) energies are
extracted from a sequence of ensemble calculations (where ensemble
weights are controlled by a single one $w$), we have shown in this work
that a single ensemble calculation is sufficient provided, of course,
that the {\it many}-weight-dependence of the ensemble xc functional is
known. The two approaches are equivalent in the exact theory
but they may give different results when density-functional approximations
are used. This is analyzed further in the rest of this section at the
EEXX level of approximation.\\

Let us first rewrite the exact individual energy expressions
within the GOKII approach~\cite{yang2017direct} (see
Eq.~(\ref{eq:gokII_weights})) where both bi- and
tri-ensemble calculations are needed for extracting the three lowest
energies. From the bi-ensemble energy
\begin{eqnarray}
E^{(w,0)}=(1-w)E_0+wE_1,
\end{eqnarray}
we can extract both ground- and first-excited-state energies as follows,
\begin{eqnarray}\label{eq:bi-ensemble-E0E1}
E_0&=&E^{(w,0)}-w\dfrac{dE^{(w,0)}}{dw},
\nonumber\\
E_1&=&E^{(w,0)}+(1-w)\dfrac{dE^{(w,0)}}{dw},
\end{eqnarray}
which is equivalent (for these two states) to a tri-ensemble
calculation where ${\tt w}_1=w$ and ${\tt
w}_2=0$. On the other hand, we have the tri-ensemble energy (with ${\tt w}_1={\tt
w}_2=w$),
\begin{eqnarray}
E^{(w,w)}=(1-2w)E_0+wE_1+wE_2,
\end{eqnarray}
from which we can extract, when combined with the bi-ensemble one, the
second-excited-state energy:
\begin{eqnarray}\label{eq:2nd_ex_state_gokII}
E_2&=&E_0-(E_1-E_0)+\dfrac{dE^{(w,w)}}{dw}
\nonumber\\
&=&E^{(w,0)}-(1+w)\dfrac{dE^{(w,0)}}{dw}
+\dfrac{dE^{(w,w)}}{dw}.
\end{eqnarray}
If, like in Sec.~\ref{sec:Hubb_dimer_theory}, we use a single
ensemble calculation instead (with ${\tt w}_1={\tt
w}_2=w$ for ease of comparison), then
individual energies will be determined as follows (see Eq.~(\ref{eq:ind_ener_exp_deriv_wft})),
\begin{eqnarray}\label{eq:E0-E2_tri-ens}
E_0&=&E^{(w,w)}-w\left.\dfrac{\partial E^{(w,{\tt w}_2)}}{\partial w}\right|_{{\tt w}_2=w}
-w\left.\dfrac{\partial E^{({\tt w}_1,w)}}{\partial w}\right|_{{\tt w}_1=w}
\nonumber\\
&=&E^{(w,w)}-w\dfrac{dE^{(w,w)}}{dw},
\end{eqnarray}
and
\begin{eqnarray}\label{eq:E0-E2_tri-ens_excited_states}
E_1&=&E_0+\left.\dfrac{\partial E^{(w,{\tt w}_2)}}{\partial w}\right|_{{\tt w}_2=w},
\nonumber\\
E_2&=&E_0+\left.\dfrac{\partial E^{({\tt w}_1,w)}}{\partial
w}\right|_{{\tt w}_1=w}.
\end{eqnarray}
Note that we use the latter expressions rather than the (equivalent)
ones in Eq.~(\ref{eq:ind_ener_with_DD}) for ease of comparison.\\

As readily seen from
Eqs.~(\ref{eq:bi-ensemble-E0E1})-(\ref{eq:E0-E2_tri-ens}), the two
approaches become identical (and equivalent to DEC~\cite{yang2017direct,sagredo2018can}) in the $w\rightarrow 0$ limit, even when
approximate ensemble energies are used. For larger $w$ values (in the
range $0<w\leq 1/3$), the two methods will
give substantially different results for the excited states
when the EEXX-only approximation is used, as illustrated
in Fig.~\ref{fig:compar_with_gok2}. While, in our (single-calculation-based)
approach, individual
energies are increasingly insensitive to the value of the tri-ensemble
weight as $U/t$ increases, the GOKII excited-state energies exhibit an
important weight dependence. Interestingly, as $w$ increases, they
become closer to the exact energies. In the large $U/t$ regime (see
the bottom panel of Fig.~\ref{fig:compar_with_gok2}), increasing $w$
restores the correct ordering of the excited states.\\

In the Hubbard dimer, the EEXX-only
individual energies can be expressed as explicit functionals of the
ensemble density (see
Appendix~\ref{appendix:dens_fun_ener_compar_gok2}), thus allowing for a
better understanding of these results. The key difference between GOKII and
the single tri-ensemble calculation approach is the ensemble density
itself. As $U/\Delta
v_{\rm ext}$ and $\Delta
v_{\rm ext}/t$ increase, the tri-ensemble density becomes closer to 1
(see Appendix~\ref{appendix:dens_fun_ener_compar_gok2}), thus explaining why tri-ensemble-based-only energies
are essentially the (weight-independent) ones obtained at the symmetric
EEXX level. Note that, in the latter case, the excited states are wrongly ordered (see
Appendix~\ref{appendix:symmetric_dimer}). On the other hand, the
bi-ensemble density varies as $1+w$ in the same asymmetry and
correlation regime. As a result, analytical expressions can be derived
for the variation in $U$ and $w$ of the
GOKII/EEXX energies (see Eqs.~(\ref{eq:bi-ens_contr_E2_largeU}) and
(\ref{eq:bi-ens_contr_E0E1_largeU})), thus providing a rationale for the
results shown in Fig.~\ref{fig:compar_with_gok2}. As proven in Appendix~\ref{appendix:dens_fun_ener_compar_gok2}, the improvement of the
second-excited-state energy as $w$ increases is exclusively due to the bi-ensemble
contribution (two first terms on the right-hand side of
Eq.~(\ref{eq:2nd_ex_state_gokII})). The good performance of GOKII/EEXX (in terms of
total excited-state energies)
may be
specific to the Hubbard dimer. Nevertheless, it clearly shows that
the choice of ensemble and extraction procedure is crucial when using
density-functional approximations.

\begin{figure}
\centering
\resizebox{0.45\textwidth}{!}{
\includegraphics[scale=1.0]{newfig4_final.pdf}
}
\caption{
Comparison of single tri-ensemble (with ${\tt w}_1={\tt w}_2=w$) and GOKII approaches at the
EEXX-only level of approximation in the
asymmetric Hubbard dimer for $U/t=2$ (top panel) and $U/t=10$ (bottom panel). Red, green and blue
colors are used for the ground, first and second excited states,
respectively. See text for further details.}
\label{fig:compar_with_gok2}
\end{figure}

\section{Conclusions and perspectives}\label{conclusions}

A generalized many-weight-dependent formulation of GOK-DFT has been
explored, thus leading to an in-principle-exact energy level extraction
procedure that applies to any (ground or excited) state in the ensemble
and relies on a single ensemble DFT
calculation. The latter consists, like a conventional DFT
calculation, in solving a single set of self-consistent KS equations
where the
orbitals are fractionally occupied (the occupation numbers are determined from the ensemble
weights). The theory has been applied to the Hubbard dimer.
The two corrections that should in principle be added to the bare KS energies (namely the global LZ shift and
a state-specific (ensemble-based) xc
derivative correction) were both shown to be important in the calculation of accurate and
weight-independent energy levels. In order to turn the method
into a practical computational tool, {\it ab initio}
many-weight-dependent xc density-functional approximations should be
developed. This can be achieved, for example, by applying GOK-DFT to finite uniform electron
gases~\cite{JCP17_T2_GLDA_EXX}. A nice feature of such model systems is that
both ground and excited states share the same density which is the
ensemble density itself. Consequently, in this particular case, the
density-functional ensemble xc energy is simply the weighted sum of the
individual-state xc energies. Work is currently in progress in this
direction. Finally, regarding the application of the theory to
photochemical processes, we would like to explore the possibility of
extracting non-adiabatic couplings from a GOK-DFT calculation. It may be
useful, for that
purpose, to extend the theory to the time-dependent linear response
regime. This is left for future work.\\

\section{Supplementary Material}

We provide complementary curves showing the variation in the first
ensemble weight ${\tt w}_1$ of individual energy levels (before and after the LZ shift)
for ${\tt w}_2=0$ or ${\tt
w}_2=1/4$ in various correlation and asymmetry regimes of the Hubbard
dimer.

\begin{acknowledgments}
The authors thank the ANR (MCFUNEX project, Grant No. ANR-14-CE06-
0014-01) for funding. E.~F. would also like to thank P.~F.~Loos for
stimulating discussions.

\end{acknowledgments}


text:
\section{Introduction}
Quasicrystals were discovered by \citet{Shechtman1984} opening up a new  and promising,
interdisciplinary research field. A comprehensive presentation of the current state of the art
 focusing on the mathematical theory of elasticity can be found in the book of \citet{Fan Book}.
 Within the scientific community of quasicrystals it is well-known that
there has not yet been an agreement, concerning the description of the
dynamic behavior of quasicrystals. Different points of view date
back to \citet{Bak1985a, Bak1985b} and to \citet{Lubensky85}. \citet{Bak1985a, Bak1985b}
argued that phason modes represent structural disorder or
structural fluctuations. \citet{Ding1993} generalized first the classical theory of linear elastodynamics towards quasicrystals. They extended the Newton's law of motion to quasicrystals describing both, phonons and phasons, with wave
propagations. This model can be found in the literature \citep{Fan and Mai, Shi2005, Fan Book} as
{\it elastodynamics of quasicrystals}. To avoid confusion, we propose to use for this model the name {\it elastodynamics of wave type}. On the other hand in analogy to hydrodynamic theory, \citep{Lubensky85, Lubensky86} proposed {\it the hydrodynamics of quasicrystals}
where phonons
and phasons play very different roles in the dynamics, since phasons are insensitive to spatial translations. In general, phasons are diffusive modes rather than propagating. \citet{Rochal2002} proposed a combination of
the two previous models, the so-called {\it minimal model of
the phonon-phason elastodynamics}, according to which the
equations of motion for the phonons are of wave type and for the
phasons are of diffusion type. During the last few years, this
model has been adopted by many researchers (e.g. \citet{RadiMariano2011, Fanetal2012}) and can usually be found under the name {\it
elasto-hydrodynamics of quasicrystals} \citep{Fan2009}.
 \par
 The aforementioned models investigate the phasons in a very different manner. Pure diffusion is very different from wave propagation without attenuation. Our primary motivation has been to formulate a model, which unifies the two approaches (the elastodynamic model of wave type and the elasto-hydrodynamic model) as asymptotic cases of it and moreover provides a theory for the phason dynamics valid in the whole range of possible wavelengths.
\par
We introduce herein {\it the elastodynamic model of wave-telegraph type} for the description of dynamics of quasicrystals. Phonons are represented by waves, and phasons by waves damped in time and propagating with finite velocity; that means the equations of motion for the phonons are partial differential equations of wave type, and for the phasons partial differential equations of telegraph type. In the telegraph type equation we encounter both the first and the second time derivative of the  phason displacement field. In fact the phason velocity field is multiplied by a crucial parameter, which is generally a tensor, the so-called {\it tensor of phason friction coefficients} which is representative for the damping character of phason dynamics and is directly related to the characteristic time of damping of the process. It is interesting that when the phason friction coefficient becomes very large (in comparison with the average mass density of the material per unit time) the telegraph type equation degenerates to a diffusion type one, while when this coefficient becomes small enough, then the telegraph type equation gives birth to the wave type phason dynamics. In other words, the elasto-hydrodynamic model and the elastodynamic model of wave type can be recovered as asymptotic cases of the elastodynamic model of wave-telegraph type by letting the phason friction coefficient to obtain extreme values. This advantage of the telegraph equation has been exploited herein to formulate our approach in the form of a model unifying the existing models and simultaneously offering the possibility to describe phason dynamics for all intermediate values of the phason friction coefficient.
\par
Based on a proper mathematical formulation, that is a variational principal formulation allowing for dissipative processes for the phason fields, we derive the equations of motion for the displacement fields in the theory of incompatible as well as of compatible
 elastodynamics of wave-telegraph type of quasicrystals. The incompatible case is important, for instance, for the
 study of dislocations and plasticity while
the compatible case is important for the study of dynamics of cracks, elastic waves, body forces, as well
as dispersion relations in quasicrystals.

\section{Phason dynamics}
This section is devoted to the description of phason dynamics. We propose that phasons should be represented by waves damped in time and propagating with finite velocity, which corresponds to a telegraph type partial differential equation.
\par
It has been written by \citet{Francoual2003, Francoual2006} in an experimental study of the dynamics of long-wavelength
phason fluctuations in quasicrystals (i-AlPdMn), that phasons are collective diffusive modes. \citet{Francoual2003, Francoual2006}  found that phason waves should
decay exponentially
 in time, with ``characteristic time" (relaxation time) proportional to the square of the phason
 wavelength. In their approach which is based on the hydrodynamics, they use the
 dispersion relation of the diffusion equation producing actually damped standing waves (see eq.~(\ref{standing})) and finding the relaxation time $\tau$ (see eq.~(\ref{theta})). However, for a standing wave,
 there is no oscillation of the examined field, no propagation and
 every term dies away, the terms of shorter wavelength dying away most rapidly \citep{Webster}. Moreover, as it has been well written in \citet{deBoissieu2008} and
\citet{de Boissieu et al 2008} ``hydrodynamics deals with excitations in the long-wavelength limit"  and
``phason modes as a result of the hydrodynamic theory are always diffusive modes, provided the wavelength is large enough".
The question that comes up at this point is what happens in the full
regime of all possible wavelengths and what has to be implemented in order for
the model to incorporate the realistic case of a finite propagation velocity.

\par
The elastodynamic model of wave type can describe phason waves with finite velocity but without attenuation. On the other hand, the hydrodynamic and the elasto-hydrodynamic model support a diffusion type equation for the dynamics of the phason fields. However, as it is well-known the diffusion equation\footnote{Diffusion equation: $\frac{\partial \psi}{\partial t}=\D_{\text {dif}}\, \Delta \psi$,
where $\Delta$ is the Laplacian and $\D_{\text {dif}}>0$ the diffusion coefficient.} supports exclusively an infinite velocity of propagation. For example, when an initial disturbance
 with compact support is applied in space then the signal is transmitted spontaneously with infinite velocity in the whole space.  In fact, using the dispersion relation, it is possible to construct standing waves damped exponentially in time or propagating waves damped exponentially in space but the initial
 disturbance (in the second case) expands at once with infinite velocity to the whole spatial structure of the problem \citep{Kneubuhl}.
 In order to model waves damped exponentially in time and propagating with finite velocity, a telegraph type equation is needed.
 A typical form of {\it the telegraph equation} is given by \citep{MF, Iwanenko, Barton, Kneubuhl}
 \begin{align}    \label{tel}
 \frac{\partial^2 \psi}{\partial t^2} +\frac{2}{\tau_ {\text {tel}}} \frac{\partial\psi}{\partial t}=c^2\Delta \psi, \quad  {\tau_ {\text {tel}}}>0
 \end{align}
 where $\Delta$ is {\it the Laplacian}, $c$ indicates a {\it wave velocity} and $\tau_{\text {tel}}$ a {\it characteristic time of damping}.  The quantity
 $\bf{{\tau_ {\text {tel}}}}$ is also called {\it the modulus of decay of the wave}. Eq.~(\ref{tel}) describes the situation where a wave propagated by the wave equation is subjected to a dissipative effect caused by the term
 $({2}/{\tau_ {\text {tel}}}){\partial\psi}/{\partial t}$ \citep{Jeffrey2003}.
\par
For the description of several physical phenomena, for example wave propagation in cables and wires, heat waves, and oceanic diffusion,
the diffusion equation has been widely used as a first approximation.
However, the diffusion equation, which is a parabolic partial differential equation,
possesses the unphysical property and disadvantage
that admits infinite velocity of propagation \citep{Strauss}. For instance,
the infinite propagation of a thermal signal; that means that the mechanism of heat conduction is
established instantaneously under all conditions, which is a paradox. A proper improvement is obtained by using the telegraph equation,
which is the simplest mathematical model combining waves and diffusion, and it has been proven that it
can describe the real physics more appropriately.  From the mathematical point of view, the telegraph equation
is a hyperbolic partial differential equation overcoming, in this way, the physically inadmissible infinite
velocity of ``information'' propagation. It is worth noting that the first model in the theory of electric signalling, that is for the description of wave
propagation in telegraphs, telephones and cables, in the 19th century was based on the diffusion equation.
It was \citet{Heaviside1876}
(see also~\citet{Heaviside, Griffiths}) who
first understood that this leads to unphysical results and then he proposed his famous equation, the so-called {\it telegraph equation} or
{\it telegrapher's equation}, to improve the theory of telegraphy. For that reason, it was written by Whittaker in the preface of \cite{Heaviside}, that it should be also properly called {\it Heaviside's equation}. Later,
the same ``story'' has been repeated in problems with heat waves~\citep{MF, Vernotte, Cattaneo, Joseph} and
oceanic diffusion~\citep{Okubo1971}.
In the case of heat waves, the telegraph equation is usually called {\it hyperbolic heat conduction equation} (HHC). Sometimes, the name {\it hyperbolic diffusion equation} is also used for the telegraph equation.
\par
For the modeling of damped phason waves the concept of a phason frictional force is adopted. This force is given in terms of an anisotropic tensor, the tensor of the phason friction coefficients, in order to capture anisotropic effects of the damped waves. The anisotropic structure of this tensor is inspired by \citet{Francoual2003}, since in experiments strong anisotropies in the phason diffuse scattering intensity have been observed. Therefore, the proposed model accounts also for anisotropy; in the sense that the phason waves can be damped anisotropically.

\section{Theoretical framework}
In this section the mathematical modeling of the proposed model is presented. We start with the necessary preliminaries concerning the incompatible elasticity theory of quasicrystals in subsection 3.1. For a more detailed exposition on this matter one can also see in \cite{Agiasofitouetal2010}. We continue in subsection 3.2 with the formulation of our problem based on a variational principle allowing for dissipation of the phason fields.
\subsection{Preliminaries}
The displacement field $\UU$ in a quasicrystal is constituted by two components, the usual {\it phonon displacement field} $u_i^\|(\Bx,t) \in E_\|$, which belongs to the physical or parallel space $E_\|$,
and the {\it phason displacement field} $u_i^\bot (\Bx, t) \in E_\bot$ that
belongs to the perpendicular space $E_\bot$, that is
\begin{equation}\label{U}
\UU=(\Bu^\|, \Bu^\bot) \ \in E_\| \oplus E_\bot.
\end{equation}
All quantities depend on the physical space coordinates $\Bx \in E_\|$ and the time~$t$. The differentiation of the displacement fields with respect to the physical coordinates gives the total distortion tensors for the phonon and phason fields, ${\beta}^{\|\, \T}_{ij}$ and ${\beta}^{\bot \T}_{ij}$, which can be decomposed into {\it the elastic distortion tensors} $\beta^\|_{ij}, \ \beta^\bot_{ij} $
and {\it the plastic distortion tensors} ${\beta}^{\|\, \P}_{ij}, \ {\beta}^{\bot \P}_{ij}$,
respectively,
\begin{alignat}{2}\label{ugrad}
u^\|_{i,j}:={\beta}^{\|\, \T}_{ij}=\beta^\|_{ij}+{\beta}^{\|\, \P}_{ij}, \qquad
u^\bot_{i,j}:={\beta}^{\bot \T}_{ij}=\beta^\bot_{ij}+{\beta}^{\bot \P}_{ij}.
\end{alignat}
A comma denotes differentiation with respect to the spatial physical coordinates. The differentiation of the displacement fields with respect to the time $t$ gives the total velocities, $v_i^{\|\, \T}$ and $v_i^{\bot\, \T}$, which can be decomposed into {\it the elastic velocities}, $v^\|_i$ and $v^\bot_i$, and {\it the plastic} or { \it initial velocities} \citep{Kossecka1975, Kossecka-deWit1977}, $v_i^{\|\, \P}$ and $v_i^{\bot\, \P}$, for the phonon and phason fields, respectively,
\begin{alignat}{2}\label{udot}
\dot{u}_i^\|:=v_i^{\|\, \T}=v_i^\|+v_i^{\|\, \P},\qquad
\dot{u}_i^\bot:=v_i^{\bot\, \T}=v_i^\bot+v_i^{\bot\, \P}.
\end{alignat}
A superimposed dot denotes differentiation with respect to time. The elastic distortion tensors and the elastic velocity vectors are the physical state quantities in the theory of quasicrystals.
\par
It should be emphasized that in the dislocation theory of
quasicrystals the incorporation of the plastic fields is essential, since a dislocation is the elementary carrier of plasticity. The plastic fields, especially the plastic distortion tensors, have to be included (see also~\citet{Agiasofitouetal2010}). For defects, the plastic distortion tensors give the so-called eigenstrains~\citep{Mura}. In the theory of linear elasticity of dislocations in quasicrystals the eigenstrains have been used by~\cite{Ding94, Ding95} and~\citet{Hu2000}. Another important aspect in the study of defects in quasicrystals that the reader should always take into consideration is that of the temperature. It is reported for example in \cite{Francoual2003} that phason fluctuations are being frozen at room temperature and plastic deformation
is possible through dislocations movements at enough high temperatures.

\subsection{Variational principal formulation with dissipation}
As we have discussed in Section 2, phasons should be modeled as damped waves. This means that dissipative processes are behind the damping. The presence of dissipative processes leads to the appearance of {\it phason frictional forces} (dissipative forces) which are linear functions of the phason velocities
\begin{align}\label{fr}
f_i^{\text{fr}}=-\D_{ij}v^\bot_j.
\end{align}
It is proven by methods of statistical physics \citep{Landau Statistical Physics} that the tensor $D_{ij}$ is symmetric, that is
\begin{align}\label{Dij}
\D_{ij}=\D_{ji}.
\end{align}
These forces can be written as the derivatives, with respect to the corresponding velocities, of the quadratic function
\begin{align}\label{R}
\RR=\frac{1}{2}\, \D_{ij}v^\bot_i v^\bot_j,
\end{align}
which is called {\it the dissipative function} \citep{ Landau Mechanics, Landau Statistical Physics, Landau Elasticity, Goldstein}.
Therefore,
\begin{align}\label{fr-R}
f_i^{\text{fr}}=-\frac{\partial \RR}{\partial v^\bot_i}.
\end{align}
The dissipative function gives the rate of dissipation of the energy in the system, making its physical significance very important.
It has been proven (see, e.g., \cite{Landau Mechanics, Landau Statistical Physics}) that $\RR$ should be positive definite ($\RR>0$), therefore the tensor $\D_{ij}$ should also be positive definite. We may call here $\D_{ij}$ the {\it tensor of phason friction coefficients}.

In the framework of incompatible elasticity theory, we consider
the Lagrangian density $\LL$ as a smooth function of the physical state quantities, that is $\LL=\LL(\Bv^\|, \Bv^\bot, \Bbeta^\|, \Bbeta^\bot)$ as well as the dissipation function density $\RR=\RR(\Bv^\bot)$ as it is given by the relation~(\ref{R}).  The corresponding Euler-Lagrange equations allowing for dissipative processes for the phason fields (\cite{Landau Mechanics, Landau Statistical Physics, Goldstein}) are
\begin{align}
&\label{EL1}\frac{\partial \LL}{\partial u^\|_i}-\frac{\partial}{\partial t}
\bigg(\frac{\partial \LL}{\partial \dot{u}^\|_i}\bigg)-
\frac{\partial}{\partial x_j}\bigg(\frac{\partial \LL}{\partial u^\|_{i,j}}\bigg)=0,\\
&\label{EL2} \frac{\partial \LL}{\partial u^\bot_i}
-\frac{\partial}{\partial t}\bigg(\frac{\partial \LL}{\partial \dot{u}^\bot_i}\bigg)-
\frac{\partial}{\partial x_j}\bigg(\frac{\partial \LL}{\partial u^\bot_{i,j}}\bigg)=\frac{\partial \RR}{\partial  \dot{u}^\bot_i}.
\end{align}
It can be seen that the Euler-Lagrange equation for the phason fields (eq.~(\ref{EL2})) differs from the usual form of the Euler-Lagrange equations by the presence of the derivative of the dissipative function density on the right-hand side. \\
Explicitly, the Lagrangian density is given in terms of the kinetic energy density $T$, the elastic energy density $W$ and the potential $V$ of the external forces
\begin{equation}
\LL=T-W-V.
\end{equation} The elastic energy density (for the unlocked state) is considered as a
quadratic function of the phonon and phason strains \citep{Ding1993, Agiasofitouetal2010}
\begin{equation}\label{w}
W=\frac{1}{2}\beta^\|_{ij}C_{ijkl}\beta^\|_{kl}+\beta^\|_{ij}D_{ijkl}\beta^\bot_{kl}+\frac{1}{2}\beta^\bot_{ij}E_{ijkl}\beta^\bot_{kl}
\end{equation}
with the tensors of elastic constants to possess the following symmetries
\begin{align}\label{Coef}
C_{ijkl}=C_{klij}=C_{ijlk}=C_{jikl},\quad D_{ijkl}=D_{jikl}, \quad E_{ijkl}=E_{klij}.
\end{align}
The kinetic energy density is considered to admit phonon and phason contributions
\begin{equation}\label{T3}
T=\frac{1}{2}\rho(v^\|_i)^2+\frac{1}{2}\rho(v^\bot_i)^2,
\end{equation}
where $\rho$ is the average mass density of the material.  \citet{Rochal2000} proposed that the phasons should possess
a generalized phason mass density, the  effective phason density  $\rho_{\text {eff}}$, which
may be different from the ordinary density $\rho$. However, later in their suggested minimal model \citep{Rochal2002}, they supported that the phason momentum in quasicrystals is not ``conserved", implying $\rho_{\text {eff}}=0$. In any case, the meaning of this quantity is not clarified and its measurement
is difficult (see, e.g., \cite{Fan Book}). In the case that the phasons indeed possess a different mass density, $\rho_{\text {eff}}$, the form of the kinetic energy density~(\ref{T3}) as well as of the corresponding derived equations bellow can be easily modified.\\
The potential $V$ can be expressed in terms of the {\it conventional (phonon) body force density} $f_i^\|$ and a  {\it generalized (phason) body force density} $f_i^\bot$ as follows
\begin{equation}\label{V3}
V=-u^\|_i f_i^\|-{u}_i^\bot f_i^\bot.
\end{equation}
The constitutive relations for the considered Lagrangian density are
\begin{align}
&p_i^\|=\frac{\partial T}{\partial v^{\|}_i}=\rho v^\|_i, \label{CR1}\\
&p^\bot_i=\frac{\partial T}{\partial v^\bot_i}=\rho v^\bot_i, \label{CR2}\\
&\sigma^\|_{ij}=\frac{\partial W}{\partial \beta^\|_{ij}}=C_{ijkl}\beta^\|_{kl}+D_{ijkl}\beta^\bot_{kl}, \label{CR3}\\
&\sigma^\bot_{ij}=\frac{\partial W}{\partial \beta^\bot_{ij}}=D_{klij}\beta^\|_{kl}+E_{ijkl}\beta^\bot_{kl} \label{CR4}.
\end{align}
In the above relations, $p_i^\|$ and $p^\bot_i$ are {\it the phonon} and {\it phason momentum vectors}
and $\sigma^\|_{ij}$ and $\sigma^\bot_{ij}$ are {\it the phonon} and {\it phason stress tensors}, respectively. We note that the phonon stress is symmetric,
$\sigma^\|_{ij}=\sigma^\|_{ji}$, but the phason stress
is asymmetric, $\sigma^\bot_{ij}\neq\sigma^\bot_{ji}$, (see also, e.g., \cite{Ding1993}).
\par

We should point out that the frictional force vector (\ref{fr}) is given as the scalar product of the anisotropic tensor $\D_{ij}$ and the phason velocity $\Bv^\bot$, in order to account for anisotropic effects of the damped waves. If we consider the case of an isotropic tensor $\D_{ij}$, that is $\D_{ij}=\D\delta_{ij}, \ \D>0$, then
\begin{align}\label{force}
f_i^{\text{fr}}=-\D {v}_i^\bot.
\end{align}
 $\D$ is the {\it phason friction coefficient}, $\D>0$. Such a force (eq.~(\ref{force}) in the compatible case) was introduced  under the name phason ``bulk force" for quasicrystals in order to capture an attenuation of diffusive type of the phason modes by \citet{Rochal2000, Rochal2002}. However, this force was actually used (\citet{Rochal2000, Rochal2002}) in the place of an external force and inserted ``by hand" in the derived equations of motion.  From the physical point of view, we notice the contradiction that a force depending on the velocity (internal field) has been considered as an external force. For that reason, a straightforward procedure, clear from the mathematical (derivation from Euler-Lagrange equations with dissipation) as well as from the physical point of view has been preferred and followed here.

\section{The equations of motion}

In this section, we first derive the field equations for the considered model in terms of the momenta and the stresses
as well as in terms of the elastic fields. Next, we give the equations of motion for the phonon and phason displacement fields in the theory of incompatible and compatible elastodynamics of wave-telegraph type of quasicrystals.
\par
The Euler-Lagrange equations~(\ref{EL1}) and (\ref{EL2}),
using the relations (\ref{ugrad}) and (\ref{udot}), the eq.~(\ref{fr-R}) as well as the
constitutive relations~(\ref{CR1})--(\ref{CR4}),
give {\it the field equations for the elastodynamic model of wave-telegraph
type of quasicrystals in terms of the momenta and the stresses}
\begin{align} \label{Eq.7}
\dot{p}_i^\|-\sigma^\|_{ij,j}&=f^\|_i,\\
\label{Eq.8}
\dot{p}_i^\bot-\sigma^\bot_{ij,j}-f_i^{\text{fr}}&=f^\bot_i.
\end{align}
The above equations, using the explicit formulas of the corresponding fields (see eqs. (\ref{CR1})--(\ref{CR4}) and eq.~(\ref{fr})),
can further give {\it the field equations in terms
of the elastic fields} $\Bv^\|, \Bv^\bot, \Bbeta^\|$ and $\Bbeta^\bot$ as follows
\begin{align}\label{Eq.9}
\rho \dot{v}^\|_i-C_{ijkl}\beta^\|_{kl,j}-D_{ijkl}\beta^\bot_{kl,j}&=f^\|_i, \\
\label{Eq.10}
\rho \dot{v}^\bot_i+\D_{ij}{v}_j^\bot-
D_{klij}\beta^\|_{kl,j}-E_{ijkl}\beta^\bot_{kl,j}&=f^\bot_i.
\end{align}
Eqs.~(\ref{Eq.9}) and (\ref{Eq.10}), using eqs.~(\ref{ugrad}) and
(\ref{udot}), give {\it the equations of motion for the phonon and phason displacement fields}
\begin{align}
\label{Eq.11}
\rho \ddot{u}_i^\|-C_{ijkl}u_{k,lj}^\|-D_{ijkl}u_{k,lj}^\bot&=
\rho \dot{v}_i^{\|\, \P}-C_{ijkl}{\beta}^{\|\, \P}_{kl,j}-D_{ijkl}{\beta}^{\bot\, \P}_{kl,j}+f_i^\|,\\
\label{Eq.12}
\rho \ddot{u}_i^\bot+\D_{ij}\dot{u}_j^\bot
-D_{klij}u_{k,lj}^\|-E_{ijkl}u_{k,lj}^\bot&=
\rho \dot{v}_i^{\bot\, \P}
+\D_{ij}{v}_j^{\bot\, \P}
-D_{klij}{\beta}^{\|\, \P}_{kl,j}-E_{ijkl}{\beta}^{\bot\, \P}_{kl,j}+f^\bot_i,
\end{align}
where the plastic fields and the phonon and phason body forces are
the source terms.
Eq.~(\ref{Eq.11}) is a partial differential equation of wave
type, or in other words an elastodynamic Navier equation and
eq.~(\ref{Eq.12}) is a partial differential equation of telegraph type. Both equations are hyperbolic. The system of the partial differential equations~(\ref{Eq.11}) and (\ref{Eq.12}) is coupled due
to the tensor of phonon-phason coupling $D_{ijkl}$. Special attention should be focused on the term $\D_{ij} \dot{u}_j^\bot$ in eq.~(\ref{Eq.12}), which contains the first time derivative of the phason displacement field. This term is responsible for the damping and appears due to the existence of the phason frictional forces. To sum up, {\it eqs.~(\ref{Eq.11}) and (\ref{Eq.12}) are the equations of motion for the displacement fields
in the theory of incompatible elastodynamics of wave-telegraph type of quasicrystals}.

\par
In the theory of dislocations \citep{Lardner}, the concept of the plastic velocity, introduced by \citet{Kossecka1969} (see also \citet{Kossecka1975, Kossecka-deWit1977}), has not been widely recognized and used. For that reason, we examine also the ``usual" case of elastoplasticity, which is obtained when
the plastic velocities are zero, that is, $v_i^{\|\, \P}=0,\quad v_i^{\bot\, \P}=0$.
Then, eqs.~(\ref{Eq.11}) and (\ref{Eq.12}) (the external forces are zero) reduce to the following equations of motion
\begin{align}
\label{Eq.11-v}
&\rho \ddot{u}_i^\|-C_{ijkl}u_{k,lj}^\|-D_{ijkl}u_{k,lj}^\bot=
-C_{ijkl}{\beta}^{\|\, \P}_{kl,j}-D_{ijkl}{\beta}^{\bot\, \P}_{kl,j},\\
\label{Eq.12-v}
&\rho \ddot{u}_i^\bot+\D_{ij}\dot{u}_j^\bot
-D_{klij}u_{k,lj}^\|-E_{ijkl}u_{k,lj}^\bot=
-D_{klij}{\beta}^{\|\, \P}_{kl,j}-E_{ijkl}{\beta}^{\bot\, \P}_{kl,j}.
\end{align}
In this case, only the plastic distortion tensors,
$\beta^{\|\, \P}_{ij}$ and $\beta^{\bot \P}_{ij}$, are non-vanishing and play
the physical role of the sources for the dislocation displacement fields, $u_i^\|$ and $u_i^\bot$.
Thus, it becomes evident that the plastic distortion tensors cannot be neglected in the study of dislocations.

\par
We now proceed to investigate the compatible case. In the compatible elastodynamics all the plastic fields are zero, that is, $\beta^{\|\, \P}_{ij}=0$, $\beta^{\bot \P}_{ij}=0$, $v_i^{\|\, \P}=0$ and $v_i^{\bot\, \P}=0$.
In this case, eqs.~(\ref{Eq.7}) and (\ref{Eq.8}), and eqs.~(\ref{Eq.9}) and (\ref{Eq.10}), are valid and eqs.~(\ref{Eq.11}) and (\ref{Eq.12}) are reduced to the following ones
\begin{align}
\label{Eq.13}
&\rho \ddot{u}_i^\|-C_{ijkl}u_{k,lj}^\|-D_{ijkl}u_{k,lj}^\bot=f_i^\|,\\
\label{Eq.14}
&\rho \ddot{u}_i^\bot+\D_{ij}\dot{u}_j^\bot
-D_{klij}u_{k,lj}^\|-E_{ijkl}u_{k,lj}^\bot=f^\bot_i.
\end{align}
The mathematical character of the above equations remains the same as of the corresponding ones
in the incompatible case, that is, eq.~(\ref{Eq.13}) is a partial differential equation of wave type
and  eq.~(\ref{Eq.14}) is of telegraph type. {\it Eqs.~(\ref{Eq.13}) and (\ref{Eq.14}) are  the equations of motion for the displacement fields
in the theory of compatible elastodynamics of wave-telegraph type of quasicrystals}.
\par
The tensor telegraph differential operator (see relation~(\ref{L-T})) is the same for all three equations~(\ref{Eq.12}), (\ref{Eq.12-v}), and (\ref{Eq.14}). Comparing the scalar telegraph operator in eq.~(\ref{tel}) with the tensor telegraph operator in relation~(\ref{L-T}), one may define {\it the tensor of characteristic time of damping}, $ \tau_{ij}$, as follows
\begin{align}\label{t-ij-tel}
   ( \tau_{ij})_{\text {tel}}=2 \rho \D_{ij}^{-1},
\end{align}
where $\D_{ij}^{-1}$ is the inverse tensor of the $\D_{ij}$. It is reasonable that the tensor of characteristic time of damping depends on the average mass density of the material and the tensor of phason friction coefficients. Particularly, we see that the anisotropy of the tensor $\D_{ij}$ influences the tensor of characteristic time of damping. It is obvious that $(\tau_{ij})_{\text{tel}}=(\tau_{ji})_{\text {tel}}$.
\par
In the isotropic case of the phason friction tensor $\D_{ij}$, that is  $\D_{ij}=\D\delta_{ij}$, {\it the characteristic time of damping} simplifies to
\begin{align}\label{t-tel}
\tau_{\text {tel}}=\frac{2 \rho}{\D},
\end{align}
and the corresponding eqs.~(\ref{Eq.10}), (\ref{Eq.12}), (\ref{Eq.12-v}), and (\ref{Eq.14}) can be easily modified.  From eq.~(\ref{t-tel}), it is clear that\footnote{The phason friction coefficient $\D$  should not be confused with the ``phason diffusion constant", $D_{phason}=1/{\tau q^2}$ \citep{de Boissieu et al 2008}. It is apparent that the two quantities have different units. $D_{phason}$ is the $\D_{\text{dif}}$ in our notation (see eq.~(\ref{theta})).}  the dimension of $\D$ is $[\D]=\frac{[M]}{[L]^3[T]}$ and coincides with the dimension of  ${1}/{\Gamma_w}$, where $\Gamma_w$ is the so-called phason kinetic coefficient introduced by \citep{Lubensky85, Lubensky86}. Hence, $[\D]=1/[\Gamma_w]$.

\par
 Making a qualitative asymptotic analysis to the telegraph type equation (see, e.g., eq.~(\ref{Eq.14})) one can see that if the friction parameter $\D_{ij}$ is very large compared with the average mass density $\rho$ per unit time, then the second term of the telegraph equation is dominant in comparison with the first term. In this case, the telegraph type equation degenerates to a diffusion type equation\footnote{In this case, the phason diffusion parameter can be interrelated with the phason friction coefficient: $[\D_{\text {dif}}]\sim[E_{ijkl}]/[\D_{ij}]$.}. If $\D_{ij}$ is very small, then the wave character is dominant. For intermediate values of the $\D_{ij}$, we have waves with attenuation. This qualitative analysis of $\D_{ij}$ shows how the proposed model constitutes a unified theory of the previously mentioned models.

\par
 It is important to notice that the characteristic time of damping $\tau_{\text{tel}}$ given by eq.~(\ref{t-tel}) is qualitatively and quantitatively different from the relaxation time (eq.~(\ref{theta})) derived from the dispersion relation of the diffusion equation \citep{Francoual2003, de Boissieu et al 2008}. The characteristic time of damping $\tau_{\text{tel}}$ is appearing in the equation of motion itself and it is not based on dispersion relations. It is remarkable that for sufficiently large wavevectors it coincides with the relaxation time $\tau$ (see eq.~(\ref{t=tel})) produced by the dispersion relation of the telegraph equation as explained in the Appendix \ref{appendixA}.
\par
  A general view into the proposed model concerning the equation of motion of the phason displacement field shows that the
consideration of two quantities is important; that of the phason momentum vector and of the phason frictional forces. If the phason momentum
is zero, then eq.~(\ref{Eq.12}) (or eq.~(\ref{Eq.12-v}), eq.~(\ref{Eq.14})) will be a partial differential equation
of diffusion type (infinite velocity of propagation) and if the phason frictional forces are zero, then eq.~(\ref{Eq.12}) (or eq.~(\ref{Eq.12-v}), eq.~(\ref{Eq.14})) will be a partial differential equation
of wave type (without damping). Hence, it is evident that we need both quantities in order to obtain the desired equation
for the phason displacement field, that is the equation of telegraph type.

\section{Conclusions}
We have proposed the elastodynamic model of wave-telegraph type to describe the dynamics of quasicrystals.
The main feature of this model is that the phason modes are described by a telegraph type partial differential
equation and consequently represent waves damped in time with finite propagation velocity. The telegraph equation has drastically different solutions compared with the solutions provided by the diffusion equation. Since the two equations of motion are coupled, the choice of a new equation for the phason fields will have a very different impulse to phonon dynamics also. Therefore, the dynamic behavior of the whole system is severely influenced by the adoption of a different equation for the description of phason dynamics.
\par
Several noteworthy features characterize the proposed model. By construction, a telegraph type equation supports wave propagation with attenuated amplitude. The influence of the damping in the dynamic behavior of the phasons is expressed by the tensor of phason friction coefficients, which gives the possibility to take into account that the phason waves can be damped anisotropically. In terms of the phason friction coefficient and the average mass density of the material an important quantity, the characteristic time of damping, has been defined. Moreover, this model provides a theory valid in the whole regime of possible wavelengths for the phasons. In addition, with the new model there is no longer the drawback of the infinite propagation velocity that exists with the equation of diffusion type.
\par
 One more advantage of a telegraph type equation is that it establishes uniformity and symmetry between the phonons and phasons as far as the corresponding initial value problems are concerned. More precisely, the implication of diffusion equation creates a lack of symmetry to the whole system. Indeed, the phonons are represented by a wave type partial differential equation which leads to a solution depending on two initial data, that is,  the initial values of both phonon displacement and total velocity fields. If the phasons are represented by a diffusion type equation, then the solution depends only on the initial condition for the displacement field (and the initial total velocity is absent). While if phasons are represented by a telegraph type equation, the solution incorporates both initial data for phason displacement and total velocity fields. Therefore, a telegraph type equation restores uniformity in the settlement of initial value problems.
\par
The understanding of the phason modes is not only important for the study of quasicrystals but also of other types of aperiodic crystals. It seems that the
proposed model could describe very well the behavior of incommensurately modulated phases. Indeed, it has been observed by \citet{deBoissieu2008} that there are two regimes for the phason modes, a long-wavelength regime where the phasons are diffusive and a shorter
wavelength regime where the phasons are damped propagating modes. Therefore, the proposed model could give a unified and broad framework for the description of aperiodic crystals. Of course, its establishment depends on its
experimental validation. We believe that the introduction of the new model will bring a new perspective, promoting the
study in the field of quasicrystals and further of aperiodic crystals.

\section*{Acknowledgements}
Dr. Eleni Agiasofitou is greatly indebted to Professor Antonios Charalambopoulos
for efficient discussions and valuable remarks concerning several aspects
of wave propagation.
The authors gratefully acknowledge grants from the Deutsche Forschungsgemeinschaft
(Grant Nos. La1974/2-1, La1974/2-2, La1974/3-1).

\begin{appendix}

\section{Dispersion relations}
\label{appendixA}
\setcounter{equation}{0}
\renewcommand{\theequation}{\thesection.\arabic{equation}}

  The dispersion relations contain a lot of information about the behavior of wave propagation inside a medium. Here, we give the dispersion relations for the diffusion equation (we examine the case of standing waves damped exponentially in time for necessary comparisons) and a typical telegraph equation. For sake of simplicity, this appendix is restricted to one space dimension.
 \par
 A plane harmonic wave moving in the $z$-direction is given by
\begin{align}\label{wave}
\psi(z,t)=A\, \e^{\ii (qz-\omega t)},
\end{align}
where $A$ is the wave amplitude, $q={2\pi}/{\lambda}$ stands for the wavenumber,  $\lambda$ is the corresponding wavelength, and $\omega$ is the circular frequency.
\par
The diffusion equation is given by
\begin{align}\label{Dif}
\frac{\partial \psi}{\partial t}=\D_{\text {dif}}\, \frac{\partial^2 \psi}{\partial z^2}, \quad \D_{\text {dif}}>0.
\end{align}
Application of eq.~(\ref{wave}) to eq.~(\ref{Dif}) yields to the dispersion relation \citep{Kneubuhl}
\begin{align}\label{disp-dif}
\omega=-\ii\, \D_{\text{dif}}\, q^2.
\end{align}
Substituting eq.~(\ref{disp-dif}) into eq.~(\ref{wave}) leads to damped standing waves represented by
\begin{align}\label{standing}
\psi(z,t)=A\, \e^{\ii qz} \e^{-\theta t}=A\, \e^{\ii qz} \e^{-{t}/{\tau}},
\end{align}
 where the damping coefficient $\theta$ is given by
\begin{align}\label{theta}
\theta=\frac{1}{\tau}= \D_{\text {dif}}\, q^2
\end{align}
with $\tau$ the corresponding relaxation time.
\par
Let us consider now the telegraph equation (\ref{tel}) in one space dimension
 \begin{align}    \label{tele}
 \frac{\partial^2 \psi}{\partial t^2} +\frac{2}{\tau_ {\text {tel}}}\, \frac{\partial\psi}{\partial t}=c^2 \frac{\partial^2 \psi}{\partial z^2}.
 \end{align}
Substituting eq.~(\ref{wave}) into eq.~(\ref{tele}), we find that $\omega$ and $q$ must satisfy the dispersion relation \citep{Jeffrey2003}
\begin{align}\label{om}
\omega^2+\frac{2\ii}{\tau_ {{\text {tel}}}}\, \omega-c^2 q^2=0.
\end{align}
We can distinguish three cases:\\
i) $c\,q > {1}/\tau_ {\text {tel}}$\\
Then, from eq.~(\ref{om}) we deduce
\begin{align}\label{omega}
\omega=-\frac{\ii}{\tau_ {\text {tel}}} \pm \sqrt{c^2\, q^2\,-\frac{1}{\tau_ {\text {tel}}^2}}
\end{align}
and consequently the solution (\ref{wave}) takes the form
\begin{align}\label{finalpsi}
\psi(z,t)=A\, \e^{-\frac{1}{\tau_ {\text {tel}}}t}\, \e^{\ii \bigl(qz \mp \sqrt{c^2q^2-
\frac{1}{\tau^2_ {\text {tel}}}}\, t \bigl)}.
\end{align}
 Eq.~(\ref{finalpsi}) represents a traveling wave damped in time with relaxation time $\tau$, which is equal to the  characteristic time of damping $\tau_ {\text {tel}}$ \begin{align}\label{t=tel}
\tau=\tau_ {\text {tel}}
\end{align}
and phase velocity
\begin{align}\label{c_p}
c_p=c\, \sqrt{1-\frac{1}{c^2\, \tau^2_{\text {tel}}\, q^2}}.
\end{align}
ii) $c\,q < {1}/\tau_ {\text {tel}}$\\
Then, we deduce from eq.~(\ref{wave}) that only standing waves damped in time are supported of the form
\begin{align}
\psi(z,t)=A\, \e^{\ii qz} \e^{-\theta t}=A\, \e^{\ii qz} \e^{-{t}/{\tau}}
\end{align}
where two relaxation times appear
\begin{align}\label{eq12}
\theta=\frac{1}{\tau}=\frac{1}{\tau_{{\text {tel}}}}\pm\sqrt{\frac{1}{\tau^2_{\text {tel}}}-c^2\,q^2}.
\end{align}
 $\theta >0$ for both values.\\
 iii) $c\,q = {1}/\tau_ {\text {tel}}$\\
In this case, it is easy to see that we have standing waves damped in time with relaxation time equal to the characteristic time of damping, $\tau=\tau_ {\text {tel}}$.
\par
Consequently, the critical value $q_0={1}/{c\, \tau_ {\text {tel}}}$ can be considered as a crucial threshold for the wavenumber $q$. If $q<q_0$, then we have only exponentially damped standing waves with two relaxation times. When $q$ overpasses $q_0$, then propagating waves (damped in time) emerge with relaxation time equals the characteristic time of damping and phase velocity dependent on the wavenumber $q$ via the relation (\ref{c_p}). If $q=q_0$, then standing waves appear with relaxation time equal to the characteristic time of damping.
\par
 All these arguments can be reformulated in terms of the wavelength $\lambda$ of the process. In case that $\lambda > \lambda_0$, where $\lambda_0=2\pi c\,\tau_{\text{tel}}$ is the critical wavelength, then standing waves appear with two characteristic times depending on $q$  via eq.~(\ref{eq12}) which are different from the relaxation time (\ref{theta}). In case that $\lambda < \lambda_0$, then propagating waves damped in time of the form (\ref{finalpsi}) emerge. This has a physical interest in case that  $\lambda_0$ is sufficiently larger than the characteristic dimension $d_0$ of the interatomic distance, so that the wavelength $\lambda$ satisfies the basic principle of the continuum theory to remain larger than $d_0$.

\section{Matrix form of the coupled equations of motion of the elastodynamic model of wave-telegraph type}
\label{appendixB}
\setcounter{equation}{0}
\renewcommand{\theequation}{\thesection.\arabic{equation}}

For a numerical implementation it is useful to rewrite the equations of motion in a matrix
form. First, we define the following tensor differential operators:
\begin{itemize}
\item
tensor wave operator:\\
\begin{align}
\label{L-W}
L_{ik}^{\text{w}}=\rho\, \delta_{ik}\pd_{tt}-C_{ijkl}\pd_j\pd_l
\end{align}
\item
tensor telegraph operator:\\
\begin{align}
\label{L-T}
L_{ik}^{\text{te}}=\rho\, \delta_{ik}\pd_{tt}+\D_{ik}\, \pd_t -E_{ijkl}\pd_j\pd_l
\end{align}
\item
tensor operator of phonon-phason coupling:\\
\begin{align}
\label{L-C}
L_{ik}^{\text{c}}=-D_{ijkl}\, \pd_j \pd_l,
\end{align}
\end{itemize}
where $\pd_t$ and $\pd_j$ are the first derivatives with respect to the time and the spatial coordinates $x_j$, respectively.
\par
Using the tensor differential operators~(\ref{L-W})--(\ref{L-C}),
the coupled equations of motion (\ref{Eq.13}) and (\ref{Eq.14}) for the compatible case
read in the matrix form
\begin{align}
\label{MR}
\begin{pmatrix}
L_{ik}^{\text{w}} & L_{ik}^{\text{c}}\\
L_{ki}^{\text{c}} & L_{ik}^{\text{te}}
\end{pmatrix}
\begin{pmatrix}
u_{k}^{\|}\\
u_{k}^{\bot}
\end{pmatrix}
=\begin{pmatrix}
f_{i}^{\|}\\
f_{i}^{\bot}
\end{pmatrix},
\end{align}
and the coupled equations of motion~(\ref{Eq.11}) and (\ref{Eq.12}) for the incompatible case are given by
\begin{align}
\label{MR2}
\begin{pmatrix}
L_{ik}^{\text{w}} & L_{ik}^{\text{c}}\\
L_{ki}^{\text{c}} & L_{ik}^{\text{te}}
\end{pmatrix}
\begin{pmatrix}
u_{k}^\|\\
u_{k}^\bot
\end{pmatrix}=
\begin{pmatrix}
\rho\,\delta_{ik} \pd_t & 0\\
 0 &\rho\, \delta_{ik}\pd_t+\D_{ik}
\end{pmatrix}
\begin{pmatrix}
v_{k}^{\|\,\P }\\
v_{k}^{\bot\,\P }
\end{pmatrix}
-\begin{pmatrix}
C_{ijkl}\pd_j & D_{ijkl}\pd_j \\
D_{klij}\pd_j  &E_{ijkl}\pd_j
\end{pmatrix}
\begin{pmatrix}
\beta_{kl}^{\|\,\P }\\
\beta_{kl}^{\bot\,\P }
\end{pmatrix}
+
\begin{pmatrix}
f_{i}^{\|}\\
f_{i}^{\bot}
\end{pmatrix}\,.
\end{align}

\end{appendix}

redpajama_set_name: RedPajamaArXiv