We discuss the physical properties and accuracy of three distinct dynamical (\ie, frequency-dependent) kernels for the computation of optical excitations within linear response theory:
i) an \textit{a priori} built kernel inspired by the dressed time-dependent density-functional theory (TDDFT) kernel proposed by Maitra and coworkers [\href{https://doi.org/10.1063/1.1651060}{J.~Chem.~Phys.~120, 5932 (2004)}],
ii) the dynamical kernel stemming from the Bethe-Salpeter equation (BSE) formalism derived originally by Strinati [\href{https://doi.org/10.1007/BF02725962}{Riv.~Nuovo Cimento 11, 1--86 (1988)}], and
iii) the second-order BSE kernel derived by Yang and coworkers [\href{https://doi.org/10.1063/1.4824907}{J.~Chem.~Phys.~139, 154109 (2013)}].
In particular, using a simple two-level model, we analyze, for each kernel, the appearance of spurious excitations, as first evidenced by Romaniello and collaborators [\href{https://doi.org/10.1063/1.3065669}{J.~Chem.~Phys.~130, 044108 (2009)}], due to the approximate nature of the kernels.
Linear response theory is a powerful approach that allows to directly access the optical excitations $\omega_s$ of a given electronic system (such as a molecule) and their corresponding oscillator strengths $f_s$ [extracted from their eigenvectors $\T{(\bX_s \bY_s)}$] via the response of the system to a weak electromagnetic field. \cite{Casida_1995}
From a practical point of view, these quantities are obtained by solving non-linear, frequency-dependent Casida-like equations in the space of single excitations and de-excitations \cite{Casida_1995}
\begin{equation}\label{eq:LR}
\begin{pmatrix}
\bR(\omega_s) &\bC(\omega_s)
\\
-\bC(-\omega_s) & -\bR(-\omega_s)
\end{pmatrix}
\begin{pmatrix}
\bX_s
\\
\bY_s
\end{pmatrix}
=
\omega_s
\begin{pmatrix}
\bX_s
\\
\bY_s
\end{pmatrix}
\end{equation}
where the explicit expressions of the resonant and coupling blocks, $\bR(\omega)$ and $\bC(\omega)$, depend on the level of approximation that one employs.
Neglecting the coupling block between the resonant and anti-resonants parts, $\bR(\omega)$ and $-\bR(-\omega)$, is known as the Tamm-Dancoff approximation (TDA).
The non-linear eigenvalue problem defined in Eq.~\eqref{eq:LR} has particle-hole symmetry which means that it is invariant via the transformation $\omega\to-\omega$, and, thanks to its non-linear nature stemming from its frequency dependence, it potentially generates more than just single excitations.
In a wave function context, introducing a spatial orbital basis $\lbrace\MO{p}\rbrace$, we assume here that the elements of the matrices defined in Eq.~\eqref{eq:LR} have the following generic form:
Here, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, and $f^{\sigma}(\omega)$ is the correlation part of the spin-resolved kernel.
As a chemist, it is maybe difficult to understand the concept of dynamical properties, the motivation behind their introduction, and their actual usefulness.
Here, we will try to give a pedagogical example showing the importance of dynamical quantities and their main purposes. \cite{Romaniello_2009b,Sangalli_2011,ReiningBook}
To do so, let us consider the usual chemical scenario where one wants to get the optical excitations of a given system.
If we assume that the operator $\bA$ has a matrix representation of size $N \times N$, this \textit{linear} set of equations yields $N$ excitation energies.
However, in practice, $N$ might be (very) large (\eg, equal to the total number of single and double excitations generated from a reference Slater determinant), and it might therefore be practically useful to recast this system as two smaller coupled systems, such that
where the blocks $\bA_1$ and $\bA_2$, of sizes $N_1\times N_1$ and $N_2\times N_2$ (with $N_1+ N_2= N$), can be associated with, for example, the single and double excitations of the system.
For example, an operator $\Tilde{\bA}_1(\omega)$ built in the single-excitation basis can potentially provide excitation energies for double excitations thanks to its frequency-dependent nature, the information from the double excitations being ``folded'' into $\Tilde{\bA}_1(\omega)$ via Eq.~\eqref{eq:row2}. \cite{ReiningBook}
This procedure converting degrees of freedom into frequency or energy dependence is very general and can be applied in various contexts. \cite{Garniron_2018,QP2}
However, because there is no free lunch, this non-linear system is obviously harder to solve than its corresponding linear analog given by Eq.~\eqref{eq:lin_sys}.
For example, assuming that $\bA_2$ is a diagonal matrix is of common practice (see, for example, Ref.~\onlinecite{Garniron_2018} and references therein).
Another of these approximations is the so-called \textit{static} approximation, where one sets the frequency to a particular value.
For example, as commonly done within the Bethe-Salpeter equation (BSE) formalism of many-body perturbation theory (MBPT), \cite{Strinati_1988}$\Tilde{\bA}_1(\omega)=\Tilde{\bA}_1\equiv\Tilde{\bA}_1(\omega=0)$.
A similar example in the context of time-dependent density-functional theory (TDDFT) \cite{Runge_1984} is provided by the ubiquitous adiabatic approximation, \cite{Tozer_2000} which neglects all memory effects by making static the exchange-correlation (xc) kernel (\ie, frequency-independent). \cite{Maitra_2016}
These approximations come with a heavy price as the number of solutions provided by the system of equations \eqref{eq:non_lin_sys} has now been reduced from $N$ to $N_1$.
Coming back to our example, in the static (or adiabatic) approximation, the operator $\Tilde{\bA}_1$ built in the single-excitation basis cannot provide double excitations anymore, and the $N_1$ excitation energies are associated with single excitations.
In the next section, we illustrate these concepts and the various tricks that can be used to recover some of these dynamical effects starting from the static eigenproblem.
Let us consider a two-level quantum system made of two orbitals \cite{Romaniello_2009b} in its singlet ground state (\ie, the lowest orbital is doubly occupied).
We will label these two orbitals, $\MO{v}$ and $\MO{c}$, as valence ($v$) and conduction ($c$) orbitals with respective one-electron Hartree-Fock (HF) energies $\e{v}$ and $\e{c}$.
The ground state $\ket{0}$ has a one-electron configuration $\ket{v\bar{v}}$, while the doubly-excited state $\ket{D}$ has a configuration $\ket{c\bar{c}}$.
There is then only one single excitation possible which corresponds to the transition $v \to c$ with different spin-flip configurations.
As usual, this can produce a singlet singly-excited state $\ket{S}=(\ket{v\bar{c}}+\ket{c\bar{v}})/\sqrt{2}$, and a triplet singly-excited state $\ket{T}=(\ket{v\bar{c}}-\ket{c\bar{v}})/\sqrt{2}$. \cite{SzaboBook}
For the sake of illustration, we will use the same numerical example throughout this study, and consider the singlet ground state of the \ce{He} atom in Pople's 6-31G basis set.
This system contains two orbitals and the numerical values of the various quantities defined above are
\begin{subequations}
\begin{align}
\e{v}& = -0.914\,127
&
\e{c}& = + 1.399\,859
\\
\ERI{vv}{vv}& = 1.026\,907
&
\ERI{cc}{cc}& = 0.766\,363
\\
\ERI{vv}{cc}& = 0.858\,133
&
\ERI{vc}{cv}& = 0.227\,670
\\
\ERI{vv}{vc}& = 0.316\,490
&
\ERI{vc}{cc}& = 0.255\,554
\end{align}
\end{subequations}
This yields the following exact singlet and triplet excitation energies
\begin{align}\label{sec:exact}
\omega_{1}^{\updw}& = 1.92145
&
\omega_{3}^{\updw}& = 3.47880
&
\omega_{1}^{\upup}& = 1.47085
\end{align}
that we are going to use a reference for the remaining of this study.
The kernel proposed by Maitra and coworkers \cite{Maitra_2004,Cave_2004} in the context of dressed TDDFT (D-TDDFT) corresponds to an \textit{ad hoc} many-body theory correction to TDDFT.
More specifically, D-TDDFT adds manually to the static kernel a frequency-dependent part by reverse-engineering the exact Hamiltonian \eqref{eq:H-exact}.
The very same idea was taking further by Huix-Rotllant, Casida and coworkers. \cite{Huix-Rotllant_2011}
The expression \eqref{eq:f-Maitra} can be easily obtained by folding the double excitation onto the single excitation starting from the Hamiltonian \eqref{eq:H-exact}, as explained in Sec.~\ref{sec:dyn}.
It is clear that one must know \textit{a priori} the structure of the Hamiltonian to construct such dynamical kernel, and this obviously hampers its applicability to realistic photochemical systems where it is sometimes hard to get a clear picture of the interplay between excited states. \cite{Boggio-Pasqua_2007}
For the two-level model, the non-linear equations defined in Eq.~\eqref{eq:LR} provides the following effective Hamiltonian
Assuming that the dynamically-screened Coulomb potential has been calculated at the random-phase approximation (RPA) level and within the Tamm-Dancoff approximation (TDA), the expression of the $\GW$ quasiparticle energy is
In Eq.~\eqref{eq:SigC}, $\Omega=\Delta\eGW{}+2\ERI{vc}{vc}$ is the sole (singlet) RPA excitation energy of the system, with $\Delta\eGW{}=\eGW{c}-\eGW{v}$.
are the elements of the dynamically-screened Coulomb potential for the resonant and coupling blocks of the dBSE Hamiltonian.
It can be easily shown that solving the equation
\begin{equation}
\det[\bH^{\dBSE}(\omega) - \omega \bI] = 0
\end{equation}
yields 6 solutions (per spin manifold): 3 pairs of frequencies opposite in sign, which corresponds to the 3 resonant states and the 3 anti-resonant states.
As mentioned in Ref.~\cite{Romaniello_2009b}, spurious solutions appears due to the approximate nature of the dBSE kernel.
Indeed, diagonalizing the exact Hamiltonian would produce two singlet solutions corresponding to the singly- and doubly-excited states, while there is only one triplet state (see discussion earlier in the section).
Therefore, there is one spurious solution for the singlet manifold and two spurious solution for the triplet manifold.
Within the static approximation, the BSE Hamiltonian is
In the static approximation, only one pair of solutions (per spin manifold) is obtained by diagonalizing $\bH^{\BSE}$.
There are, like in the dynamical case, opposite in sign.
Therefore, the static BSE Hamiltonian does not produce spurious excitations but misses the (singlet) double excitation.
Enforcing the TDA, which corresponds to neglecting the coupling term between the resonant and anti-resonant part of the BSE Hamiltonian, \ie, $C(\omega)=0$, allows to remove some of these spurious excitations.
In this case, the excitation energies are obtained by solving the simple equation $R(\omega)-\omega=0$, which yields two solutions for each spin manifold.
There is thus only one spurious excitation in the triplet manifold, the two solutions of the singlet manifold corresponding to the single and double excitations.
Another way to access dynamical effects while staying in the static framework is to use perturbation theory.
To do so, one must decompose the dBSE Hamiltonian into a (zeroth-order) static part and a dynamical perturbation, such that
This corresponds to a dynamical correction to the static excitations, and the TDA can be applied to the dynamical correction, a scheme we label as dTDA in the following.
$\det[\bH^{\dBSE}(\omega)-\omega\bI]$ as a function of $\omega$ for both the singlet (red) and triplet (blue) manifolds.
\label{fig:dBSE}
}
\end{figure}
%%% %%% %%% %%%
Figure \ref{fig:dBSE} shows the three resonant solutions (for the singlet and triplet spin manifold) of the dynamical BSE Hamiltonian $\bH(\omega)$ defined in Eq.~\eqref{eq:HBSE}, the curve being invariant with respect to the transformation $\omega\to-\omega$ (electron-hole symmetry).
Numerically, we find
\begin{align}
\omega_{1,\updw}^{\dBSE}& = 1.90527
&
\omega_{2,\updw}^{\dBSE}& = 2.78377
&
\omega_{3,\updw}^{\dBSE}& = 4.90134
\end{align}
for the singlet states, and
\begin{align}
\omega_{1,\upup}^{\dBSE}& = 1.46636
&
\omega_{2,\upup}^{\dBSE}& = 2.76178
&
\omega_{3,\upup}^{\dBSE}& = 4.91545
\end{align}
for the triplet states.
it is interesting to mention that, around $\omega=\omega_1^{\sigma}$ ($\sigma=$$\updw$ or $\upup$), the slope of the curves depicted in Fig.~\ref{fig:dBSE} is small, while the two other solutions, $\omega_2^{\sigma}$ and $\omega_3^{\sigma}$, stem from poles and consequently the slope is very large around these frequency values.
Diagonalizing the static BSE Hamiltonian yields the following singlet and triplet excitation energies:
\begin{align}
\omega_{1,\updw}^{\BSE}& = 1.92778
&
\omega_{1,\upup}^{\BSE}& = 1.48821
\end{align}
which shows that the physical single excitation stemming from the dynamical BSE Hamiltonian is the lowest one for each spin manifold, \ie, $\omega_1^{\updw}$ and $\omega_1^{\upup}$.
%%% FIGURE 2 %%%
\begin{figure}
\includegraphics[width=\linewidth]{dBSE-TDA}
\caption{
$\det[\bH^{\TDAdBSE}(\omega)-\omega\bI]$ as a function of $\omega$ for both the singlet (red) and triplet (blue) manifolds within the TDA.
\label{fig:dBSE-TDA}
}
\end{figure}
%%% %%% %%% %%%
Figure \ref{fig:dBSE-TDA} shows the same curves as Fig.~\ref{fig:dBSE} but in the TDA.
As one can see, the spurious solution $\omega_2^{\sigma}$ has disappeared, and two pairs of solutions remain for each spin manifold.
This evidences that BSE reproduces qualitatively well the singlet and triplet single excitations, but quite badly the double excitation which is off by more than 1 hartree.
Here, we follow a different strategy and compute the dynamical second-order BSE kernel as illustrated by Yang and collaborators \cite{Zhang_2013}, and Rebolini and Toulouse \cite{Rebolini_2016}.
First, let us compute the second-order quasiparticle energies, which reads
This expression can be easily obtained in the present case by the substitution $\Omega\to\e{c}-\e{v}$ which transforms the $GW$ self-energy into its GF2 analog.
The static Hamiltonian for this theory is just the usual RPAx (or TDHF) Hamiltonian, \ie,
The dynamical part of the kernel for BSE2 (that we will call dRPAx for notational consistency) is a bit ugly but it simplifies greatly in the case of the present model to yield
In Ref.~\cite{Sangalli_2011}, Sangalli proposed a dynamical kernel (based on the second RPA) without (he claims) spurious excitations thanks to the design of a number-conserving approach which correctly describes particle indistinguishability and Pauli exclusion principle.
We will first start by writing down explicitly this kernel as it is given in obscure physicist notations.
The dynamical BSE Hamiltonian with Sangalli's kernel is
The author thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481) for financial support.