BSEdyn/Notes/BSEdyn-notes.tex

660 lines
27 KiB
TeX

\documentclass[aps,prb,reprint,noshowkeys,superscriptaddress]{revtex4-2}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
\usepackage[version=4]{mhchem}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}
\usepackage[
colorlinks=true,
citecolor=blue,
breaklinks=true
]{hyperref}
\urlstyle{same}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{HTML}{009900}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\SI}{\textcolor{blue}{supplementary material}}
\newcommand{\QP}{\textsc{quantum package}}
\newcommand{\T}[1]{#1^{\intercal}}
% coordinates
\newcommand{\br}{\mathbf{r}}
\newcommand{\dbr}{d\br}
% methods
\newcommand{\evGW}{ev$GW$}
\newcommand{\qsGW}{qs$GW$}
\newcommand{\GOWO}{$G_0W_0$}
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\xc}{\text{xc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\co}{\text{x}}
%
\newcommand{\Norb}{N_\text{orb}}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}
\newcommand{\IS}{\lambda}
% operators
\newcommand{\hH}{\Hat{H}}
% methods
\newcommand{\KS}{\text{KS}}
\newcommand{\HF}{\text{HF}}
\newcommand{\RPA}{\text{RPA}}
\newcommand{\RPAx}{\text{RPAx}}
\newcommand{\dRPAx}{\text{dRPAx}}
\newcommand{\BSE}{\text{BSE}}
\newcommand{\TDABSE}{\text{BSE(TDA)}}
\newcommand{\dBSE}{\text{dBSE}}
\newcommand{\TDAdBSE}{\text{dBSE(TDA)}}
\newcommand{\GW}{GW}
\newcommand{\GF}{\text{GF2}}
\newcommand{\stat}{\text{stat}}
\newcommand{\dyn}{\text{dyn}}
\newcommand{\TDA}{\text{TDA}}
% energies
\newcommand{\Enuc}{E^\text{nuc}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\EHF}{E^\text{HF}}
\newcommand{\EBSE}{E^\text{BSE}}
\newcommand{\EcRPA}{E_\text{c}^\text{RPA}}
\newcommand{\EcBSE}{E_\text{c}^\text{BSE}}
% orbital energies
\newcommand{\e}[1]{\eps_{#1}}
\newcommand{\eHF}[1]{\eps^\text{HF}_{#1}}
\newcommand{\eKS}[1]{\eps^\text{KS}_{#1}}
\newcommand{\eQP}[1]{\eps^\text{QP}_{#1}}
\newcommand{\eGW}[1]{\eps^{GW}_{#1}}
\newcommand{\eGF}[1]{\eps^{\text{GF2}}_{#1}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
% Matrix elements
\newcommand{\Sig}[1]{\Sigma_{#1}}
\newcommand{\SigGW}[1]{\Sigma^{\GW}_{#1}}
\newcommand{\SigGF}[1]{\Sigma^{\GF}_{#1}}
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\sERI}[2]{[#1|#2]}
% excitation energies
\newcommand{\OmRPA}[1]{\Omega_{#1}^{\text{RPA}}}
\newcommand{\OmRPAx}[1]{\Omega_{#1}^{\text{RPAx}}}
\newcommand{\OmBSE}[1]{\Omega_{#1}^{\text{BSE}}}
\newcommand{\spinup}{\downarrow}
\newcommand{\spindw}{\uparrow}
\newcommand{\singlet}{\uparrow\downarrow}
\newcommand{\triplet}{\uparrow\uparrow}
% Matrices
\newcommand{\bO}{\mathbf{0}}
\newcommand{\bH}{\mathbf{H}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bI}{\mathbf{1}}
\newcommand{\bb}{\mathbf{b}}
\newcommand{\bA}{\mathbf{A}}
\newcommand{\bB}{\mathbf{B}}
\newcommand{\bc}{\mathbf{c}}
\newcommand{\bx}{\mathbf{x}}
% units
\newcommand{\IneV}[1]{#1 eV}
\newcommand{\InAU}[1]{#1 a.u.}
\newcommand{\InAA}[1]{#1 \AA}
\newcommand{\kcal}{kcal/mol}
\DeclareMathOperator*{\argmax}{argmax}
\DeclareMathOperator*{\argmin}{argmin}
% orbitals, gaps, etc
\newcommand{\updw}{\uparrow\downarrow}
\newcommand{\upup}{\uparrow\uparrow}
\newcommand{\eps}{\varepsilon}
\newcommand{\IP}{I}
\newcommand{\EA}{A}
\newcommand{\HOMO}{\text{HOMO}}
\newcommand{\LUMO}{\text{LUMO}}
\newcommand{\Eg}{E_\text{g}}
\newcommand{\EgFun}{\Eg^\text{fund}}
\newcommand{\EgOpt}{\Eg^\text{opt}}
\newcommand{\EB}{E_B}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\begin{document}
\title{Dynamical Kernels}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\begin{abstract}
Similar to the ubiquitous adiabatic approximation in time-dependent density-functional theory, the static approximation, which substitutes a dynamical (\ie, frequency-dependent) kernel by its static limit, is usually enforced in most implementations of the Bethe-Salpeter equation (BSE) formalism.
Here, going beyond the static approximation, we compute the dynamical correction in the electron-hole screening for molecular excitation energies thanks to a renormalized first-order perturbative correction to the static BSE excitation energies.
The present dynamical correction goes beyond the plasmon-pole approximation as the dynamical screening of the Coulomb interaction is computed exactly.
Moreover, we investigate quantitatively the effect of the Tamm-Dancoff approximation by computing both the resonant and anti-resonant dynamical corrections to the BSE excitation energies.
%\\
%\bigskip
%\begin{center}
% \boxed{\includegraphics[width=0.5\linewidth]{TOC}}
%\end{center}
%\bigskip
\end{abstract}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The concept of dynamical quantities}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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{ReiningBook}.
To do so, let us consider the usual chemical scenario where one wants to get the neutral excitations of a given system.
In most cases, this can be done by solving a set of linear equations of the form
\begin{equation}
\label{eq:lin_sys}
\bA \bc = \omega \bc
\end{equation}
where $\omega$ is one of the neutral excitation energies of the system associated with the transition vector $\bc$.
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, and it might therefore be practically useful to recast this system as two smaller coupled systems, such that
\begin{equation}
\label{eq:lin_sys_split}
\begin{pmatrix}
\bA_1 & \T{\bb} \\
\bb & \bA_2 \\
\end{pmatrix}
\begin{pmatrix}
\bc_1 \\
\bc_2 \\
\end{pmatrix}
= \omega
\begin{pmatrix}
\bc_1 \\
\bc_2 \\
\end{pmatrix}
\end{equation}
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.
Note that this \textit{exact} decomposition does not alter, in any case, the values of the excitation energies, not their eigenvectors.
Solving separately each row of the system \eqref{eq:lin_sys_split} yields
\begin{subequations}
\begin{gather}
\label{eq:row1}
\bA_1 \bc_1 + \T{\bb} \bc_2 = \omega \bc_1
\\
\label{eq:row2}
\bc_2 = (\omega \bI - \bA_2)^{-1} \bb \bc_1
\end{gather}
\end{subequations}
Substituting Eq.~\eqref{eq:row2} into Eq.~\eqref{eq:row1} yields the following effective \textit{non-linear}, frequency-dependent operator
\begin{equation}
\label{eq:non_lin_sys}
\Tilde{\bA}_1(\omega) \bc_1 = \omega \bc_1
\end{equation}
with
\begin{equation}
\Tilde{\bA}_1(\omega) = \bA_1 + \T{\bb} (\omega \bI - \bA_2)^{-1} \bb
\end{equation}
which has, by construction, exactly the same solutions than the linear system \eqref{eq:lin_sys} but a smaller dimension.
For example, an operator $\Tilde{\bA}_1(\omega)$ built in the basis of single excitations 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}.
How have we been able to reduce the dimension of the problem while keeping the same number of solutions?
To do so, we have transformed a linear operator $\bA$ into a non-linear operator $\Tilde{\bA}_1(\omega)$ by making it frequency dependent.
In other words, we have sacrificed the linearity of the system in order to obtain a new, non-linear systems of equations of smaller dimension.
This procedure converting degrees of freedom into frequency or energy dependence is very general and can be applied in various contexts.
Thanks to its non-linearity, Eq.~\eqref{eq:non_lin_sys} can produce more solutions than its actual dimension.
However, because there is no free lunch, this non-linear system is obviously harder to solve than its corresponding linear analogue given by Eq.~\eqref{eq:lin_sys}.
Nonetheless, approximations can be now applied to Eq.~\eqref{eq:non_lin_sys} in order to solve it efficiently.
One of these approximations is the so-called \textit{static} approximation, which corresponds to fix the frequency to a particular value.
For example, as commonly done within the Bethe-Salpeter formalism, $\Tilde{\bA}_1(\omega) = \Tilde{\bA}_1 \equiv \Tilde{\bA}_1(\omega = 0)$.
In such a way, the operator $\Tilde{\bA}_1$ is made linear again by removing its frequency-dependent nature.
This approximation comes 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 $K$ to $K_1$.
Coming back to our example, in the static approximation, the operator $\Tilde{\bA}_1$ built in the basis of single excitations cannot provide double excitations anymore, and the only $K_1$ excitation energies are associated with single excitations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A two-level model}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let us consider a two-level quantum system made of two orbitals \cite{Romaniello_2009b}.
We will label these two orbitals as valence ($v$) and conduction ($c$) orbitals with respective one-electron energies $\e{v}$ and $\e{c}$.
In a more quantum chemical language, these correspond to the HOMO and LUMO orbitals (respectively).
The ground state has a one-electron configuration $v\bar{v}$, while the doubly-excited state has a configuration $c\bar{c}$.
There is then only one single excitation which corresponds to the transition $v \to c$.
As usual, this can produce a singlet singly-excited state of configuration $(v\bar{c} + c\bar{v})/\sqrt{2}$, and a triplet singly-excited state of configuration $(v\bar{c} - c\bar{v})/\sqrt{2}$ \cite{SzaboBook}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Dynamical BSE kernel}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Within many-body perturbation theory (MBPT), one can easily compute the quasiparticle energies associated with the valence and conduction orbitals.
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
\begin{equation}
\e{p}^{\GW} = \e{p} + Z_{p}^{\GW} \SigGW{p}(\e{p})
\end{equation}
where $p = v$ or $c$,
\begin{equation}
\label{eq:SigC}
\SigGW{p}(\omega) = \frac{2 \ERI{pv}{vc}^2}{\omega - \e{v} + \Omega} + \frac{2 \ERI{pc}{cv}^2}{\omega - \e{c} - \Omega}
\end{equation}
are the correlation parts of the self-energy associated with wither the valence of conduction orbitals,
\begin{equation}
Z_{p}^{\GW} = \qty( 1 - \left. \pdv{\SigGW{p}(\omega)}{\omega} \right|_{\omega = \e{p}} )^{-1}
\end{equation}
is the renormalization factor, and
\begin{equation}
\ERI{pq}{rs} = \iint \MO{p}(\br) \MO{q}(\br) \frac{1}{\abs{\br - \br'}} \MO{r}(\br') \MO{s}(\br') d\br d\br'
\end{equation}
are the usual (bare) two-electron integrals.
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}$.
One can now build the dynamical Bethe-Salpeter equation (dBSE) Hamiltonian, which reads
\begin{equation} \label{eq:HBSE}
\bH^{\dBSE}(\omega) =
\begin{pmatrix}
R(\omega) & C(\omega)
\\
-C(-\omega) & -R(-\omega)
\end{pmatrix}
\end{equation}
with
\begin{subequations}
\begin{gather}
R(\omega) = \Delta\eGW{} + 2 \sigma \ERI{vc}{cv} - W_R(\omega)
\\
C(\omega) = 2 \sigma \ERI{vc}{cv} - W_C(\omega)
\end{gather}
\end{subequations}
($\sigma = 1$ for singlets and $\sigma = 0$ for triplets) and
\begin{subequations}
\begin{gather}
W_R(\omega) = \ERI{vv}{cc} + \frac{4 \ERI{vv}{vc} \ERI{vc}{cc}}{\omega - \Omega - \Delta\eGW{}}
\\
W_C(\omega) = \ERI{vc}{cv} + \frac{4 \ERI{vc}{cv}^2}{\omega - \Omega}
\end{gather}
\end{subequations}
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
\begin{equation}
\bH^{\BSE} =
\begin{pmatrix}
R^{\stat} & C^{\stat}
\\
-C^{\stat} & -R^{\stat}
\end{pmatrix}
\end{equation}
with
\begin{subequations}
\begin{gather}
R^{\stat} = \Delta\eGW{} + 2 \sigma \ERI{vc}{vc} - W_R(\omega = \Delta\eGW{})
\\
C^{\stat} = 2 \sigma \ERI{vc}{vc} - W_C(\omega = 0)
\end{gather}
\end{subequations}
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
\begin{equation}
\bH^{\dBSE}(\omega) = \underbrace{\bH^{\BSE}}_{\bH^{(0)}} + \underbrace{\qty[ \bH^{\dBSE}(\omega) - \bH^{\BSE} ]}_{\bH^{(1)}}
\end{equation}
Thanks to (renormalized) first-order perturbation theory, one gets
\begin{equation}
\omega_{1,\sigma}^{\BSE1} = \omega_{1,\sigma}^{\BSE} + Z_{1} \T{\bV} \cdot \qty[ \bH^{\dBSE}(\omega = \omega_{1,\sigma}^{\BSE}) - \bH^{\BSE} ] \cdot \bV
\end{equation}
where
\begin{equation}
\bV =
\begin{pmatrix}
X \\ Y
\end{pmatrix}
\end{equation}
are the eigenvectors of $\bH^{\BSE}$, and
\begin{equation}
Z_{1} = \qty{ 1 - \T{\bV} \cdot \left. \pdv{\bH^{\dBSE}(\omega)}{\omega} \right|_{\omega = \omega_{1,\sigma}^{\BSE}} \cdot \bV }^{-1}
\end{equation}
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.
We now take a numerical example by considering the singlet ground state of the \ce{He} atom in the 6-31G basis set.
This system contains two orbitals and the numerical values of the various quantities defined above are
\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}
which yields
\begin{align}
\Omega & = 2.769\,327
&
\eGW{v} & = -0.863\,700
&
\eGW{c} & = +1.373\,640
\end{align}
%%% FIGURE 1 %%%
\begin{figure}
\includegraphics[width=\linewidth]{dBSE}
\caption{
$\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.
Numerically, we have
\begin{align}
\omega_{1,\updw}^{\TDAdBSE} & = 1.94005
&
\omega_{3,\updw}^{\TDAdBSE} & = 4.90117
\end{align}
for the singlet states, and
\begin{align}
\omega_{1,\upup}^{\TDAdBSE} & = 1.47070
&
\omega_{3,\upup}^{\TDAdBSE} & = 4.91517
\end{align}
while the static values are
\begin{align}
\omega_{1,\updw}^{\TDABSE} & = 1.95137
&
\omega_{1,\upup}^{\TDABSE} & = 1.49603
\end{align}
It is now instructive to provide the exact results, \ie, the excitation energies obtained by diagonalizing the exact Hamiltonian in the same basis set.
A quick configuration interaction with singles and doubles (CISD) calculation provide the following excitation energies:
\begin{align}
\omega_{1}^{\updw} & = 1.92145
&
\omega_{1}^{\upup} & = 1.47085
&
\omega_{3}^{\updw} & = 3.47880
\end{align}
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.
All these numerical results are gathered in Table \ref{tab:BSE}.
The perturbatively-corrected values are also reported, which shows that this scheme is very efficient at reproducing the dynamical value.
Note that, although the BSE1(dTDA) value is further from the dBSE value than BSE1, it is quite close to the exact excitation energy.
%%% TABLE I %%%
\begin{table*}
\caption{BSE singlet and triplet excitation energies at various levels of theory.
\label{tab:BSE}
}
\begin{ruledtabular}
\begin{tabular}{|c|ccccccc|c|}
Singlets & BSE & BSE1 & BSE1(dTDA) & dBSE & BSE(TDA) & BSE1(TDA) & dBSE(TDA) & Exact \\
\hline
$\omega_1$ & 1.92778 & 1.90022 & 1.91554 & 1.90527 & 1.95137 & 1.94004 & 1.94005 & 1.92145 \\
$\omega_2$ & & & & 2.78377 & & & & \\
$\omega_3$ & & & & 4.90134 & & & 4.90117 & 3.47880 \\
\hline
Triplets & BSE & BSE1 & BSE1(dTDA) & dBSE & BSE(TDA) & BSE1(TDA) & dBSE(TDA) & Exact \\
\hline
$\omega_1$ & 1.48821 & 1.46860 & 1.46260 & 1.46636 & 1.49603 & 1.47070 & 1.47070 & 1.47085 \\
$\omega_2$ & & & & 2.76178 & & & & \\
$\omega_3$ & & & & 4.91545 & & & 4.91517 & \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Second-order BSE kernel}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
\begin{equation}
\eGF{p} = \e{p} + Z_{p}^{\GF} \SigGF{p}(\e{p})
\end{equation}
where the second-order self-energy is
\begin{equation}
\label{eq:SigGF}
\SigGF{p}(\omega) = \frac{2 \ERI{pv}{vc}^2}{\omega - \e{v} + \e{c} - \e{v}} + \frac{2 \ERI{pc}{cv}^2}{\omega - \e{c} - (\e{c} - \e{v})}
\end{equation}
and
\begin{equation}
Z_{p}^{\GF} = \qty( 1 - \left. \pdv{\SigGF{p}(\omega)}{\omega} \right|_{\omega = \e{p}} )^{-1}
\end{equation}
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,
\begin{equation}
\bH^{\RPAx} =
\begin{pmatrix}
A^{\stat} & B^{\stat}
\\
-B^{\stat} & -A^{\stat}
\end{pmatrix}
\end{equation}
with
\begin{subequations}
\begin{gather}
A^{\stat} = \Delta\eGF{} + 2 \sigma \ERI{vc}{vc} - \ERI{vv}{cc}
\\
B^{\stat} = 2 \sigma \ERI{vc}{vc} - \ERI{vc}{cv}
\end{gather}
\end{subequations}
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
\begin{equation}
\bH^{\dRPAx} = \bH^{\RPAx} +
\begin{pmatrix}
A(\omega) & B
\\
-B & -A(-\omega)
\end{pmatrix}
\end{equation}
with
\begin{subequations}
\begin{gather}
A^{\updw}(\omega) = - \frac{4 \ERI{cv}{vv} \ERI{vc}{cc} - \ERI{vc}{cc}^2 - \ERI{cv}{vv}^2 }{\omega - 2 \Delta\eGF{}}
\\
B^{\updw} = - \frac{4 \ERI{vc}{cv}^2 - \ERI{cc}{cc} \ERI{vc}{cv} - \ERI{vv}{vv} \ERI{vc}{cv} }{2 \Delta\eGF{}}
\end{gather}
\end{subequations}
and
\begin{subequations}
\begin{gather}
A^{\upup}(\omega) = - \frac{ \ERI{vc}{cc}^2 + \ERI{cv}{vv}^2 }{\omega - 2 \Delta\eGF{}}
\\
B^{\upup} = - \frac{\ERI{cc}{cc} \ERI{vc}{cv} + \ERI{vv}{vv} \ERI{vc}{cv} }{2 \Delta\eGF{}}
\end{gather}
\end{subequations}
Note that the coupling blocks $B$ are frequency independent, as they should.
This has an important consequence as this lack of frequency dependence removes one of the spurious pole.
The singlet manifold has then the right number of excitations.
However, one spurious triplet excitation remains (see Fig.~\ref{fig:dBSE2}).
Numerical results for the two-level model are reported in Table \ref{tab:RPAx} with the usual approximations and perturbative treatments.
In the case of dRPAx, the perturbative partitioning is simply
\begin{equation}
\bH^{\dRPAx}(\omega) = \underbrace{\bH^{\RPAx}}_{\bH^{(0)}} + \underbrace{\qty[ \bH^{\dRPAx}(\omega) - \bH^{\RPAx} ]}_{\bH^{(1)}}
\end{equation}
This might not be the smartest way of decomposing the Hamiltonian though but it seems to work fine.
%%% TABLE II %%%
\begin{table*}
\caption{RPAx singlet and triplet excitation energies at various levels of theory.
\label{tab:RPAx}
}
\begin{ruledtabular}
\begin{tabular}{|c|ccccccc|c|}
Singlets & RPAx & RPAx1 & RPAx1(dTDA) & dRPAx & RPAx(TDA) & RPAx1(TDA) & dRPAx(TDA) & Exact \\
\hline
$\omega_1$ & 1.84903 & 1.90941 & 1.90950 & 1.90362 & 1.86299 & 1.92356 & 1.92359 & 1.92145 \\
$\omega_2$ & & & & & & & & \\
$\omega_3$ & & & & 4.47124 & & & 4.47097 & 3.47880 \\
\hline
Triplets & RPAx & RPAx1 & RPAx1(dTDA) & dRPAx & RPAx(TDA) & RPAx1(TDA) & dRPAx(TDA) & Exact \\
\hline
$\omega_1$ & 1.38912 & 1.44285 & 1.44304 & 1.42564 & 1.40765 & 1.46154 & 1.46155 & 1.47085 \\
$\omega_2$ & & & & & & & & \\
$\omega_3$ & & & & 4.47797 & & & 4.47767 & \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%%% %%% %%% %%%
%%% FIGURE 3 %%%
\begin{figure}
\includegraphics[width=\linewidth]{dBSE2}
\caption{
$\det[\bH^{\dRPAx}(\omega) - \omega \bI]$ as a function of $\omega$ for both the singlet (red) and triplet (blue) manifolds.
\label{fig:dBSE2}
}
\end{figure}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Sangalli's kernel}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
\begin{equation}
\bH^\text{NC}(\omega) =
\begin{pmatrix}
H(\omega) & K(\omega)
\\
-K(-\omega) & -H(-\omega)
\end{pmatrix}
\end{equation}
with
\begin{subequations}
\begin{gather}
H_{ia,jb}(\omega) = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \Xi_{ia,jb} (\omega)
\\
K_{ia,jb}(\omega) = \Xi_{ia,bj} (\omega)
\end{gather}
\end{subequations}
and
\begin{subequations}
\begin{gather}
\Xi_{ia,jb} (\omega) = \sum_{m \neq n} \frac{ C_{ia,mn} C_{jb,mn} }{\omega - ( \omega_{m} + \omega_{n})}
\\
C_{ia,mn} = \frac{1}{2} \sum_{jb,kc} \qty{ \qty[ \ERI{ij}{kc} \delta_{ab} + \ERI{kc}{ab} \delta_{ij} ] \qty[ R_{m,jc} R_{n,kb}
+ R_{m,kb} R_{n,jc} ] }
\end{gather}
\end{subequations}
where $R_{m,ia}$ are the elements of the RPA eigenvectors.
Here $i$, $j$, and $k$ are occupied orbitals, $a$, $b$, and $c$ are unoccupied orbitals, and $m$ and $n$ label single excitations.
For the two-level model, Sangalli's kernel reads
\begin{align}
H(\omega) & = \Delta\eGW{} + \Xi_H (\omega)
\\
K(\omega) & = \Xi_K (\omega)
\end{align}
\begin{gather}
\Xi_H (\omega) = 2 \frac{ [\ERI{vv}{vc} + \ERI{vc}{cc}]^2 }{\omega - 2\omega_1}
\\
\Xi_C (\omega) = 0
\end{gather}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Take-home messages}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
What have we learnt here?
%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
%%%%%%%%%%%%%%%%%%%%%%%%
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.
% BIBLIOGRAPHY
\bibliography{../BSEdyn}
\end{document}