BSEdyn/BSEdyn.tex

962 lines
50 KiB
TeX

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
\usepackage[version=4]{mhchem}
%\usepackage{natbib}
%\usepackage[extra]{tipa}
%\bibliographystyle{achemso}
%\AtBeginDocument{\nocite{achemso-control}}
\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{\xavier}[1]{\textcolor{darkgreen}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\XB}[1]{\xavier{(\underline{\bf XB}: #1)}}
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\SI}{\textcolor{blue}{supporting information}}
\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{\BSE}{\text{BSE}}
\newcommand{\GW}{GW}
\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{\eGOWO}[1]{\eps^\text{\GOWO}_{#1}}
\newcommand{\eGW}[1]{\eps^{GW}_{#1}}
\newcommand{\eevGW}[1]{\eps^\text{\evGW}_{#1}}
\newcommand{\eGnWn}[2]{\eps^\text{\GnWn{#2}}_{#1}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
\newcommand{\tOm}[2]{\Tilde{\Omega}_{#1}^{#2}}
% Matrix elements
\newcommand{\A}[2]{A_{#1}^{#2}}
\newcommand{\tA}[2]{\Tilde{A}_{#1}^{#2}}
\newcommand{\B}[2]{B_{#1}^{#2}}
\renewcommand{\S}[1]{S_{#1}}
\newcommand{\ABSE}[2]{A_{#1}^{#2,\text{BSE}}}
\newcommand{\BBSE}[2]{B_{#1}^{#2,\text{BSE}}}
\newcommand{\ARPA}[2]{A_{#1}^{#2,\text{RPA}}}
\newcommand{\BRPA}[2]{B_{#1}^{#2,\text{RPA}}}
\newcommand{\ARPAx}[2]{A_{#1}^{#2,\text{RPAx}}}
\newcommand{\BRPAx}[2]{B_{#1}^{#2,\text{RPAx}}}
\newcommand{\G}[1]{G_{#1}}
\newcommand{\LBSE}[1]{L_{#1}}
\newcommand{\XiBSE}[1]{\Xi_{#1}}
\newcommand{\Po}[1]{P_{#1}}
\newcommand{\W}[2]{W_{#1}^{#2}}
\newcommand{\tW}[2]{\widetilde{W}_{#1}^{#2}}
\newcommand{\Wc}[1]{W^\text{c}_{#1}}
\newcommand{\vc}[1]{v_{#1}}
\newcommand{\Sig}[1]{\Sigma_{#1}}
\newcommand{\SigGW}[1]{\Sigma^{GW}_{#1}}
\newcommand{\Z}[1]{Z_{#1}}
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\sERI}[2]{[#1|#2]}
%% bold in Table
\newcommand{\bb}[1]{\textbf{#1}}
\newcommand{\rb}[1]{\textbf{\textcolor{red}{#1}}}
\newcommand{\gb}[1]{\textbf{\textcolor{darkgreen}{#1}}}
% 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{\bI}{\mathbf{1}}
\newcommand{\bvc}{\mathbf{v}}
\newcommand{\bSig}{\mathbf{\Sigma}}
\newcommand{\bSigX}{\mathbf{\Sigma}^\text{x}}
\newcommand{\bSigC}{\mathbf{\Sigma}^\text{c}}
\newcommand{\bSigGW}{\mathbf{\Sigma}^{GW}}
\newcommand{\be}{\mathbf{\epsilon}}
\newcommand{\beGW}{\mathbf{\epsilon}^{GW}}
\newcommand{\beGnWn}[1]{\mathbf{\epsilon}^\text{\GnWn{#1}}}
\newcommand{\bde}{\mathbf{\Delta\epsilon}}
\newcommand{\bdeHF}{\mathbf{\Delta\epsilon}^\text{HF}}
\newcommand{\bdeGW}{\mathbf{\Delta\epsilon}^{GW}}
\newcommand{\bOm}[1]{\mathbf{\Omega}^{#1}}
\newcommand{\bA}[1]{\mathbf{A}^{#1}}
\newcommand{\btA}[1]{\Tilde{\mathbf{A}}^{#1}}
\newcommand{\bB}[1]{\mathbf{B}^{#1}}
\newcommand{\bX}[2]{\mathbf{X}_{#1}^{#2}}
\newcommand{\bY}[2]{\mathbf{Y}_{#1}^{#2}}
\newcommand{\bZ}[2]{\mathbf{Z}_{#1}^{#2}}
\newcommand{\bK}{\mathbf{K}}
\newcommand{\bP}[1]{\mathbf{P}^{#1}}
% 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{\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{\si}{\sigma}
\newcommand{\sis}{\sigma^*}
\newcommand{\pis}{\pi^*}
\newcommand{\ra}{\rightarrow}
\newcommand{\hpsi}{\Hat{\psi}}
\newcommand{\ha}{\Hat{a}}
\newcommand{\tchi}{\Tilde{\chi}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\NEEL}{Universit\'e Grenoble Alpes, CNRS, Institut NEEL, F-38042 Grenoble, France}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\begin{document}
\title{ \textcolor{red}{Assessing} Dynamical Corrections to the Bethe-Salpeter Equation}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Xavier \surname{Blase}}
\email{xavier.blase@neel.cnrs.fr }
\affiliation{\NEEL}
\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}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%
The Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} is to the $GW$ approximation \cite{Hedin_1965,Golze_2019} of many-body perturbation theory (MBPT) \cite{Martin_2016} what time-dependent density-functional theory (TD-DFT) \cite{Runge_1984,Casida_1995} is to Kohn-Sham density-functional theory (KS-DFT), \cite{Hohenberg_1964,Kohn_1965} an affordable way of computing the neutral excitations of a given electronic system.
In recent years, it has been shown to be a valuable tool for computational chemists with a large number of systematic benchmark studies on large molecular systems appearing in the literature \cite{Korbel_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Krause_2017,Gui_2018} (see Ref.~\onlinecite{Blase_2018} for a recent review).
Taking the optical gap (\ie, the lowest optical excitation energy) as an example, BSE builds on top of a $GW$ calculation by adding up excitonic effects $\EB$ to the $GW$ HOMO-LUMO gap
\begin{equation}
\Eg^{\GW} = \eps_{\LUMO}^{\GW} - \eps_{\HOMO}^{\GW},
\end{equation}
which is itself a corrected version of the Kohn-Sham (KS) gap
\begin{equation}
\Eg^{\KS} = \eps_{\LUMO}^{\KS} - \eps_{\HOMO}^{\KS} \ll \Eg^{\GW} \approx \EgFun,
\end{equation}
in order to approximate the optical gap
\begin{equation}
\EgOpt = E_1^{N} - E_0^{N} = \EgFun + \EB,
\end{equation}
where
\begin{equation} \label{eq:Egfun}
\EgFun = I^N - A^N
\end{equation}
is the the fundamental gap, \cite{Bredas_2014} $I^N = E_0^{N-1} - E_0^N$ and $A^N = E_0^{N+1} - E_0^N$ being the ionization potential and the electron affinity of the $N$-electron system.
Here, $E_s^{N}$ is the total energy of the $s$th excited state of the $N$-electron system, and $E_0^N$ corresponds to its ground-state energy.
Because the excitonic effect corresponds physically to the stabilization implied by the attraction of the excited electron and its hole left behind, we have $\EgOpt < \EgFun$.
Most of BSE implementations rely on the so-called static approximation, which approximates the dynamical (\ie, frequency-dependent) BSE kernel by its static limit.
In complete analogy with the ubiquitous adiabatic approximation in TD-DFT where the exchange-correlation (xc) kernel is made static, one key consequence of the static approximation within BSE is that double (and higher) excitations are completely absent from the BSE spectrum.
Indeed, a frequency-dependent kernel has the ability to create additional poles in the response function, which describe states with a multiple-excitation character, and, in particular, double excitations.
Although these double excitations are usually experimentally dark (which means that they usually cannot be observed in photo-absorption spectroscopy), these states play, indirectly, a key role in many photochemistry mechanisms, \cite{Boggio-Pasqua_2007} and are particularly important in the faithful description of the ground state of open-shell molecules. \cite{Casida_2005,Romaniello_2009a,Huix-Rotllant_2011,Loos_2020c}
They are, moreover, a real challenge for high-level computational methods. \cite{Loos_2018a,Loos_2019,Loos_2020b,Loos_2020c}
Double excitations play also a significant role in the correct location of the excited states of polyenes that are closely related to rhodopsin which is involved in the visual transduction. \cite{Olivucci_2010,Robb_2007,Manathunga_2016}
In butadiene, for example, while the bright $1 ^1B_u$ state has a clear ($\HOMO \ra \LUMO$) single-excitation character, the dark $2 ^1A_g$ state includes a substantial fraction of doubly-excited character from the $\HOMO^2 \ra \LUMO^2$ double excitation (roughly $30\%$), yet dominant contributions from the $\HOMO-1 \ra \LUMO$ and $\HOMO \ra \LUMO+1$ single excitations. \cite{Maitra_2004,Cave_2004,Saha_2006,Watson_2012,Shu_2017,Barca_2018a,Barca_2018b,Loos_2019}
Going beyond the static approximation is tricky and very few groups have dared to take the plunge. \cite{Strinati_1988,Rohlfing_2000,Sottile_2003,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Huix-Rotllant_2011,Zhang_2013,Rebolini_2016,Olevano_2019,Lettmann_2019}
Nonetheless, it is worth mentioning the seminal work of Strinati, \cite{Strinati_1988} who \titou{bla bla bla.}
Following Strinati's footsteps, Rohlfing and coworkers have developed an efficient way of taking into account, thanks to first-order perturbation theory, the dynamical effects via a plasmon-pole approximation combined with the Tamm-Dancoff approximation (TDA). \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b}
With such a scheme, they have been able to compute the excited states of biological chromophores, showing that taking into account the electron-hole dynamical screening is important for an accurate description of the lowest $n \ra \pi^*$ excitations. \cite{Ma_2009a,Ma_2009b,Baumeier_2012b}
Indeed, studying PYP, retinal and GFP chromophore models, Ma \textit{et al.}~found that \textit{``the influence of dynamical screening on the excitation energies is about $0.1$ eV for the lowest $\pi \ra \pis$ transitions, but for the lowest $n \ra \pis$ transitions the influence is larger, up to $0.25$ eV.''} \cite{Ma_2009b}
A similar conclusion was reached in Ref.~\onlinecite{Ma_2009a}.
Zhang \textit{et al.}~have studied the frequency-dependent second-order Bethe-Salpeter kernel and they have observed an appreciable improvement over configuration interaction with singles (CIS), time-dependent Hartree-Fock (TDHF), and adiabatic TD-DFT results. \cite{Zhang_2013}
Rebolini and Toulouse have performed a similar investigation in a range-separated context, and they have reported a modest improvement over its static counterpart. \cite{Rebolini_2016,Rebolini_PhD}
In these two latter studies, they also followed a (non-self-consistent) perturbative approach within the TDA with a renormalization of the first-order perturbative correction.
It is important to note that, although all the studies mentioned above are clearly going beyond the static approximation of BSE, they are not able to recover additional excitations as the perturbative treatment makes ultimately the BSE kernel static.
However, it does permit to recover, for transitions with a dominant single-excitation character, additional relaxation effects coming from higher excitations (and, in particular, non-interacting double excitations).
These higher excitations would be explicitly present in the BSE Hamiltonian by ``unfolding'' the dynamical BSE kernel, and one would recover a linear eigenvalue problem with, nonetheless, a much larger dimension.
Based on a rather simple model (the Hubbard dimer) which permits to analytically solve the dynamical equations, Romaniello and coworkers \cite{Romaniello_2009b,Sangalli_2011} evidenced that one can genuinely access additional excitations by solving the non-linear, frequency-dependent eigenvalue problem.
For this particular system, it was shown that a BSE kernel based on the random-phase approximation (RPA) produces indeed double excitations, but also unphysical excitations. \cite{Romaniello_2009b}
The appearance of these spurious excitations was attributed to the self-screening problem. \cite{Romaniello_2009a}
This was fixed in a follow-up paper by Sangalli \textit{et al.} \cite{Sangalli_2011} thanks to the design of a number-conserving approach based on the second RPA.
Finally, let us mention efforts to borrow ingredients from BSE in order to go beyond the adiabatic approximation of TD-DFT.
For example, Huix-Rotllant and Casida \cite{Casida_2005,Huix-Rotllant_2011} proposed a nonadiabatic correction to the xc kernel by using the formalism of superoperators, which includes as a special case the dressed TD-DFT method of Maitra and coworkers. \cite{Maitra_2004,Cave_2004,Elliott_2011,Maitra_2012}
Following a similar strategy, Romaniello \textit{et al.} \cite{Romaniello_2009b} took advantages of the dynamically-screened Coulomb potential from BSE to obtain a dynamic TD-DFT kernel.
In this regard, MBPT provides key insights about what is missing in adiabatic TD-DFT, as discussed at length by Casida and Huix-Rotllant in Ref.~\onlinecite{Casida_2016}.
In the present study, we extend the work of Rohlfing and coworkers \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b} by proposing a renormalized first-order perturbative correction to the static BSE excitation energies.
Importantly, our 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 TDA by computing both the resonant and anti-resonant dynamical corrections to the BSE excitation energies.
Unless otherwise stated, atomic units are used.
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theory}
%%%%%%%%%%%%%%%%%%%%%%%%
In this Section, following Strinati's seminal work, \cite{Strinati_1988} we first derive in some details the theoretical foundations leading to the dynamical BSE.
Additional details about this derivation are provided as {\SI}.
We present, in a second step, the perturbative implementation of the dynamical correction as compared to the standard static approximation.
%================================
\subsection{General dynamical BSE}
%=================================
The two-particle correlation function $L(1,2; 1',2')$ --- a central quantity in the BSE formalism --- relates the variation of the one-body Green's function $G(1,1')$ with respect to an external non-local perturbation $U(2',2)$, \ie,
\begin{equation}
iL(1,2; 1',2') = \pdv{G(1,1')}{U(2',2)},
\end{equation}
where, \eg, $1 \equiv (\bx_1 t_1)$ is a space-spin plus time composite variable.
The relation between $G$ and the one-body charge density $\rho(1) = -i G(1,1^+)$ provides a direct connection with the density-density susceptibility $\chi(1,2) = L(1,2;1^+,2^+)$ at the core of TD-DFT.
(The notation $1^+$ means that the time $t_1$ is taken at $t_1^{+} = t_1 + 0^+$, where $0^+$ is a small positive infinitesimal.)
The two-body correlation function $L$ satisfies the self-consistent BSE \cite{Strinati_1988}
\begin{multline} \label{eq:BSE}
L(1,2; 1',2') = L_0(1,2;1',2')
\\
+ \int d3456 \, L_0(1,4;1',3) \Xi(3,5;4,6) L(6,2;5,2'),
\end{multline}
where
\begin{subequations}
\begin{align}
\label{eq:L0}
iL_0(1, 4; 1', 3) & = G(1, 3)G(4, 1'),
\\
\label{eq:L}
iL(1,2; 1',2') & = - G_2(1,2;1',2') + G(1,1') G(2,2'),
\end{align}
\end{subequations}
can be expressed as a function of the one- and two-body Green's functions
\begin{subequations}
\begin{align}
\label{eq:G1}
G(1,2) & = - i \mel{N}{T [ \hpsi(1) \hpsi^{\dagger}(2) ] }{N},
\\
\label{eq:G2}
G_2(1,2;1',2') & = - \mel{N}{T [ \hpsi(1) \hpsi(2) \hpsi^{\dagger}(2') \hpsi^{\dagger}(1') ]}{N},
\end{align}
\end{subequations}
and
\begin{equation}
\Xi(3,5;4,6) = i \fdv{[v_\text{H}(3) \delta(3,4) + \Sigma_\text{xc}(3,4)]}{G(6,5)}
\end{equation}
is the BSE kernel that takes into account the self-consistent variation of the Hartree potential
\begin{equation}
v_\text{H}(1) = - i \int d2 \, v(1,2) G(2,2^+),
\end{equation}
[where $\delta$ is Dirac's delta function and $v$ is the bare Coulomb operator] and the xc self-energy $ \Sigma_\text{xc}$ with respect to the variation of $G$.
In Eqs.~\eqref{eq:G1} and \eqref{eq:G2}, the field operators $\Hat{\psi}(\bx t)$ and $\Hat{\psi}^{\dagger}(\bx't')$ remove and add (respectively) an electron to the $N$-electron ground state $\ket{N}$ in space-spin-time positions ($\bx t$) and ($\bx't'$), while $T$ is the time-ordering operator.
The resolution of the dynamical BSE equation \cite{Strinati_1988} starts with the expansion of $L_0$ and $L$ [see Eqs.~\eqref{eq:L0} and \eqref{eq:L}] over the complete orthonormalized set of $N$-electron excited states $\ket{N,s}$ (with $\ket{N,0} \equiv \ket{N}$).
In the optical limit of instantaneous electron-hole creation and destruction, imposing $t_{2'} = t_2^+$ and $t_{1'} = t_1^+$, and using the relation between the field operators in their time-dependent (Heisenberg) and time-independent (Schr\"{o}dinger) representations, \eg,
\begin{equation}
\hpsi(1) = e^{ i \hH t_1 } \hpsi(\bx_1) e^{-i \hH t_1 },
\end{equation}
($\hH$ being the exact many-body Hamiltonian), one gets
\begin{equation}
\begin{split}
iL(1,2; 1',2')
& = \theta(+\tau_{12}) \sum_{s > 0} \chi_s(\bx_1,\bx_{1'}) \tchi_s(\bx_2,\bx_{2'}) e^{ - i \Om{s}{} \tau_{12} }
\\
& - \theta(-\tau_{12}) \sum_{s > 0} \chi_s(\bx_2,\bx_{2'}) \tchi_s(\bx_1,\bx_{1'}) e^{ + i \Om{s}{} \tau_{12} },
\end{split}
\end{equation}
where $\tau_{12} = t_1 - t_2$, $\theta$ is the Heaviside step function, and
\begin{subequations}
\begin{align}
\chi_s(\bx_1,\bx_{2}) & = \mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N,s},
\\
\tchi_s(\bx_1,\bx_{2}) & = \mel{N,s}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N}.
\end{align}
\end{subequations}
The $\Om{s}{}$'s are the neutral excitation energies of interest.
Picking up the $e^{+i \Om{s}{} t_2 }$ component in $L(1,2; 1',2')$ and $L(6,2;5,2')$, simplifying further by $\tchi_s(\bx_2,\bx_{2'})$ on both side of the BSE [see Eq.~\eqref{eq:BSE}], we are left with seeking the $e^{-i \Om{s}{} t_1 }$ Fourier component associated with the right-hand side of a modified dynamical BSE, which reads
\begin{multline} \label{eq:BSE_2}
\mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,s} e^{ - i \Om{s}{} t_1 }
\theta ( \tau_{12} )
\\
= \int d3456 \, L_0(1,4;1',3) \Xi(3,5;4,6)
\\
\times \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
\theta [\min(t_5,t_6) - t_2].
\end{multline}
For the lowest neutral excitation energies falling in the fundamental gap of the system (\ie, $\Om{s}{} < \EgFun$ due to excitonic effects), $L_0(1,2;1',2')$ cannot contribute to the $e^{-i \Om{s}{} t_1 }$ response since its lowest excitation energy is precisely the fundamental gap [see Eq.~\eqref{eq:Egfun}].
Consequently, special care has to be taken for high-lying excited states (like core or Rydberg excitations) where additional terms have to be taken into account (see Refs.~\onlinecite{Strinati_1982,Strinati_1984}).
Dropping the (space/spin) variables, the Fourier components with respect to $t_1$ of $L_0(1,4;1',3)$ reads
\begin{align} \label{eq:iL0}
[iL_0]( \omega_1 )
= \int \frac{d\omega}{2\pi} \; G\qty(\omega - \frac{\omega_1}{2} ) G\qty( {\omega} + \frac{\omega_1}{2} )
e^{ i \omega \tau_{34} } e^{ i \omega_1 t^{34} },
\end{align}
with $\tau_{34} = t_3 - t_4$ and $t^{34} = (t_3 + t_4)/2$.
We now adopt the Lehman representation of the one-body Green's function in the quasiparticle approximation, \ie,
\begin{equation} \label{eq:G-Lehman}
G(\bx_1,\bx_2 ; \omega) = \sum_p \frac{ \MO{p}(\bx_1) \MO{p}^*(\bx_2) } { \omega - \e{p} + i \eta \times \text{sgn} (\e{p} - \mu) },
\end{equation}
where $\eta$ is a positive infinitesimal and $\mu$ is the chemical potential.
The $\e{p}$'s in Eq.~\eqref{eq:G-Lehman} are quasiparticle energies (\ie, proper addition/removal energies) and the $\MO{p}(\bx)$'s are their associated one-body (spin)orbitals.
In the following, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$ and $q$ indicate arbitrary orbitals.
Projecting the Fourier component $L_0(\bx_1,4;\bx_{1'},3; \omega_1 = \Om{s}{} )$ onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$ yields
\begin{multline} \label{eq:iL0bis}
\iint d\bx_1 d\bx_{1'} \, \MO{a}^*(\bx_1) \MO{i}(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Om{s}{})
\\
=
\frac{ \MO{a}^*(\bx_3) \MO{i}(\bx_4) e^{i \Om{s}{} t^{34} }} { \Om{s}{} - ( \e{a} - \e{i} ) + i \eta }
\qty[ \theta( \tau_{34} ) e^{i ( \e{i} + \frac{\Om{s}{}}{2}) \tau_{34} } + \theta( - \tau_{34} ) e^{i (\e{a} - \frac{\Om{s}{}}{2}) \tau_{34} } ].
\end{multline}
% and $(i,j)$/$(a,b)$ index occupied/virtual orbitals, respectively.
As a final step, we express the terms $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}')] }{N,s}$ and $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}$ from Eq.~\eqref{eq:BSE_2} in the standard electron-hole product (or single-excitation) space.
% with $(6,5) \rightarrow (5,5) \; \text{or} \; (3,4)$ when multiplied by $\delta(5,6)$ or $\delta(3,6) \delta(4,5)$, respectively.
This is done by expanding the field operators over a complete orbital basis of creation/destruction operators.
For example, we have
\begin{multline} \label{eq:spectral65}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
\\
= - \qty( e^{ -i \Omega_s t^{65} } ) \sum_{pq} \MO{p}(\bx_6) \MO{q}^*(\bx_5)
\mel{N}{\ha_q^{\dagger} \ha_p}{N,s}
\\
\times \qty[ \theta( \tau_{65} ) e^{- i ( \e{p} - \frac{\Om{s}{}}{2} ) \tau_{65} } + \theta( - \tau_{65} ) e^{ - i ( \e{q} + \frac{\Om{s}{}}{2}) \tau_{65} } ],
\end{multline}
with $t^{65} = (t_5 + t_6)/2$ and $\tau_{65} = t_6 -t_5$.
%with a similar expression for $\mel{N}{T [\hpsi(\bx_3) \hpsi^{\dagger}(\bx_4)] }{N,s}$.
%================================
\subsection{Dynamical BSE within the $GW$ approximation}
%=================================
Adopting now the $GW$ approximation for the xc self-energy, \ie,
\begin{equation}
\Sigma_\text{xc}^{\GW}(1,2) = i G(1,2) W(1^+,2),
\end{equation}
leads to the following simplified BSE kernel
\begin{equation} \label{eq:Xi_GW}
\Xi(3,5;4,6) = v(3,6) \delta(3,4) \delta(5,6) - W(3^+,4) \delta(3,6) \delta(4,5),
\end{equation}
where $W$ is the dynamically-screened Coulomb operator.
The $GW$ quasiparticle energies $\eGW{p}$ are usually good approximations to the removal/addition energies $\e{p}$ introduced in Eq.~\eqref{eq:G-Lehman}.
%Selecting $(p,q)=(j,b)$ yields the largest components
%$X_{jb}^{s} = \mel{N}{\ha_j^{\dagger} \ha_b}{N,s}$, while $(p,q)=(b,j)$ yields much weaker
%$Y_{jb}^{s} = \mel{N}{\ha_b^{\dagger} \ha_j}{N,s}$ contributions.
%Neglecting the $Y_{jb}^{s}$ weights leads to the Tamm-Dancoff approximation (TDA).
%Working out similar expressions for $\mel{N}{T [\hpsi(5) \hpsi^{\dagger}(5)] }{N,s}$ and $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1'})] }{N,s}$,
Substituting Eqs.~\eqref{eq:iL0bis}, \eqref{eq:spectral65}, and \eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, and projecting onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$, one gets after a few tedious manipulations (see {\SI}) the dynamical BSE (dBSE):
\begin{equation} \label{eq:BSE-final}
\begin{split}
( \eGW{a} - \eGW{i} - \Om{s}{} ) X_{ia,s}
& + \sum_{jb} \qty[ \kappa \ERI{ia}{jb} - \widetilde{W}_{ij,ab}(\Om{s}{}) ] X_{jb,s} \\
& + \sum_{jb} \qty[ \kappa \ERI{ia}{bj} - \widetilde{W}_{ib,aj}(\Om{s}{}) ] Y_{jb,s}
= 0,
\end{split}
\end{equation}
with $X_{jb,s} = \mel{N}{\ha_j^{\dagger} \ha_b}{N,s}$ and $Y_{jb,s} = \mel{N}{\ha_b^{\dagger} \ha_j}{N,s}$, and where $\kappa = 2 $ or $0$ for singlet and triplet excited states (respectively).
Neglecting the anti-resonant terms, $Y_{jb,s}$, in the dBSE, which are much smaller than their resonant counterparts, $X_{jb,s}$, leads to the well-known Tamm-Dancoff approximation (TDA).
In Eq.~\eqref{eq:BSE-final},
\begin{equation}
\ERI{ia}{jb} = \iint d\br d\br' \, \MO{i}^*(\br) \MO{a}(\br) v(\br -\br') \MO{j}^*(\br') \MO{b}(\br'),
\end{equation}
are the bare two-electron integrals in the spatial orbital basis $\lbrace \MO{p}(\br{}) \rbrace$, and
\begin{multline} \label{eq:wtilde}
\widetilde{W}_{ij,ab}(\Om{s}{})
= \frac{ i }{ 2 \pi} \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega)
\\
\times \qty[ \frac{1}{ \Omega_{ib}^s - \omega + i \eta } + \frac{1}{ \Omega_{ja}^{s} + \omega + i\eta } ],
\end{multline}
is an effective dynamically-screened Coulomb potential, \cite{Romaniello_2009b} where $\Om{pq}{s} = \Om{s}{} - ( \eGW{q} - \eGW{p} )$ and
\begin{equation}
W_{ij,ab}({\omega})
= \iint d\br d\br' \, \MO{i}(\br) \MO{j}^*(\br) W(\br ,\br'; \omega) \MO{a}^*(\br') \MO{b}(\br').
\end{equation}
\xavier{A second coupled equation for the $(X_{ia}^{s}, Y_{ia}^{s} )$ vector can be obtained by projecting now onto the $\mel{N}{T \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1'})}{N,s}$ left-hand side and right-hand-side of the BSE, leading to : }
%================================
\subsection{Dynamical screening}
%=================================
In the present study, we consider the exact spectral representation of $W$ at the RPA level:
\begin{multline}
\label{eq:W}
W_{ij,ab}(\omega)
= \ERI{ij}{ab} + 2 \sum_m \sERI{ij}{m} \sERI{ab}{m}
\\
\times \qty[ \frac{1}{ \omega-\Om{m}{\RPA} + i\eta } - \frac{1}{ \omega + \Om{m}{\RPA} - i\eta } ],
\end{multline}
where $m$ labels single excitations, and
\begin{equation}
\label{eq:sERI}
\sERI{pq}{m} = \sum_{ia} \ERI{pq}{ia} (\bX{m}{\RPA} + \bY{m}{\RPA})_{ia}
\end{equation}
are the spectral weights.
In Eqs.~\eqref{eq:W} and \eqref{eq:sERI}, $\OmRPA{m}{}$ and $(\bX{m}{\RPA} + \bY{m}{\RPA})$ are RPA neutral excitations and their corresponding transition vectors computed by solving the (static) linear response problem
\begin{equation}
\label{eq:LR-RPA}
\begin{pmatrix}
\bA{\RPA} & \bB{\RPA} \\
-\bB{\RPA} & -\bA{\RPA} \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{m}{\RPA} \\
\bY{m}{\RPA} \\
\end{pmatrix}
=
\OmRPA{m}
\begin{pmatrix}
\bX{m}{\RPA} \\
\bY{m}{\RPA} \\
\end{pmatrix},
\end{equation}
with
\begin{subequations}
\begin{align}
\label{eq:LR_RPA-A}
\A{ia,jb}{\RPA} & = \delta_{ij} \delta_{ab} (\e{a} - \e{i}) + 2 \ERI{ia}{jb},
\\
\label{eq:LR_RPA-B}
\B{ia,jb}{\RPA} & = 2 \ERI{ia}{bj},
\end{align}
\end{subequations}
where the $\e{p}$'s are taken as the HF orbital energies in the case of $G_0W_0$ \cite{Hybertsen_1985a, Hybertsen_1986} or as the $GW$ quasiparticle energies in the case of self-consistent schemes such as ev$GW$. \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011}
The RPA matrices $\bA{\RPA}$ and $\bB{\RPA}$ in Eq.~\eqref{eq:LR-RPA} are of size $\Nocc \Nvir \times \Nocc \Nvir$, where $\Nocc$ and $\Nvir$ are the number of occupied and virtual orbitals (\ie, $\Norb = \Nocc + \Nvir$ is the total number of spatial orbitals), respectively, and $\bX{m}{\RPA}$, and $\bY{m}{\RPA}$ are (eigen)vectors of length $\Nocc \Nvir$.
The analysis of the poles of the integrand in Eq.~\eqref{eq:wtilde} yields
\begin{multline}
\widetilde{W}_{ij,ab}( \Om{s}{} )
= \ERI{ij}{ab} + 2 \sum_m \sERI{ij}{m} \sERI{ab}{m}
\\
\times \qty[ \frac{1}{\Om{ib}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{s} - \Om{m}{\RPA} + i\eta} ].
\end{multline}
Due to excitonic effects, the lowest BSE excitation energy, $\Om{1}{}$, stands lower than the lowest RPA excitation energy, $\Om{1}{\RPA}$, so that, $\Om{ib}{s} - \Om{m}{\RPA} < 0 $ and $\widetilde{W}_{ij,ab}(\Om{s}{})$ has no resonances.
Furthermore, $\Om{ib}{s}$ and $\Om{ja}{s}$ are necessarily negative quantities for in-gap low-lying BSE excitations.
Thus, we have $\abs*{\Omega_{ib}^{s} - \Om{m}{\RPA}} > \Omega_m^{\RPA}$.
As a consequence, we observe a reduction of the electron-hole screening, \ie, an enhancement of electron-hole binding energy, as compared to the standard static BSE, and yields smaller (red-shifted) excitation energies.
%================================
\subsection{Perturbative dynamical correction}
%=================================
From a more practical point of view, to compute the BSE excitation energies of a closed-shell system, one must solve the following (non-linear) dynamical (\ie, frequency-dependent) response problem \cite{Strinati_1988}
\begin{equation}
\label{eq:LR-dyn}
\begin{pmatrix}
\bA{}(\Om{s}{}) & \bB{}(\Om{s}{})
\\
-\bB{}(\titou{-}\Om{s}{}) & -\bA{}(\titou{-}\Om{s}{})
\\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{s}{} \\
\bY{s}{} \\
\end{pmatrix}
=
\Om{s}{}
\begin{pmatrix}
\bX{s}{} \\
\bY{s}{} \\
\end{pmatrix},
\end{equation}
where the dynamical matrices $\bA{}$ and $\bB{}$ have the same $\Nocc \Nvir \times \Nocc \Nvir$ size than their RPA counterparts.
Same comment applies to the eigenvectors $\bX{s}{}$, and $\bY{s}{}$ of length $\Nocc \Nvir$.
Note that, due to its non-linear nature, Eq.~\eqref{eq:LR-dyn} may provide more than one solution for each value of $s$. \cite{Romaniello_2009b,Sangalli_2011,Martin_2016}
Accordingly to Eq.~\eqref{eq:BSE-final}, the BSE matrix elements in Eq.~\eqref{eq:LR-dyn} read
\begin{subequations}
\begin{align}
\label{eq:BSE-Adyn}
\A{ia,jb}{}(\Om{s}{}) & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \kappa \ERI{ia}{jb} - \tW{ij,ab}{}(\Om{s}{}),
\\
\label{eq:BSE-Bdyn}
\B{ia,jb}{}(\Om{s}{}) & = \kappa \ERI{ia}{bj} - \tW{ib,aj}{}(\Om{s}{}).
\end{align}
\end{subequations}
%\begin{equation}
% \ERI{pq}{rs} = \iint \frac{\MO{p}(\br{}) \MO{q}(\br{}) \MO{r}(\br{}') \MO{s}(\br{}')}{\abs*{\br{} - \br{}'}} \dbr{} \dbr{}'
%\end{equation}
%are the bare two-electron integrals in the molecular orbital basis $\lbrace \MO{p}(\br{}) \rbrace_{1 \le p \le \Norb}$, and the dynamically-screened Coulomb potential reads
%\begin{multline}
%\label{eq:W}
% \W{ij,ab}{}(\omega) = \ERI{ij}{ab} + 2 \sum_m^{\Nocc \Nvir} \sERI{ij}{m} \sERI{ab}{m}
% \\
% \times \qty(\frac{1}{\omega - \OmRPA{m}{} - \eGW{ib} + i \eta} + \frac{1}{\omega - \OmRPA{m}{} - \eGW{ja} + i \eta}).
%\end{multline}
Now, let us decompose, using basic perturbation theory, the non-linear eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static (hence linear) reference and a first-order dynamic (hence non-linear) perturbation, such that
\begin{multline}
\label{eq:LR-PT}
\begin{pmatrix}
\bA{}(\Om{s}{}) & \bB{}(\Om{s}{}) \\
-\bB{}(\titou{-}\Om{s}{}) & -\bA{}(\titou{-}\Om{s}{}) \\
\end{pmatrix}
\\
=
\begin{pmatrix}
\bA{(0)} & \bB{(0)}
\\
-\bB{(0)} & -\bA{(0)}
\\
\end{pmatrix}
+
\begin{pmatrix}
\bA{(1)}(\Om{s}{}) & \bB{(1)}(\Om{s}{}) \\
-\bB{(1)}(\titou{-}\Om{s}{}) & -\bA{(1)}(\titou{-}\Om{s}{}) \\
\end{pmatrix},
\end{multline}
with
\begin{subequations}
\begin{align}
\label{eq:BSE-A0}
\A{ia,jb}{(0)} & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \kappa \ERI{ia}{jb} - \W{ij,ab}{\text{stat}},
\\
\label{eq:BSE-B0}
\B{ia,jb}{(0)} & = \kappa \ERI{ia}{bj} - \W{ib,aj}{\text{stat}}.
\end{align}
\end{subequations}
and
\begin{subequations}
\begin{align}
\label{eq:BSE-A1}
\A{ia,jb}{(1)}(\Om{s}{}) & = - \tW{ij,ab}{}(\Om{s}{}) + \W{ij,ab}{\text{stat}},
\\
\label{eq:BSE-B1}
\B{ia,jb}{(1)}(\Om{s}{}) & = - \tW{ib,aj}{}(\Om{s}{}) + \W{ib,aj}{\text{stat}},
\end{align}
\end{subequations}
where we have defined the static version of the screened Coulomb potential
\begin{equation}
\label{eq:Wstat}
\W{ij,ab}{\text{stat}} = W_{ij,ab}(\omega = 0) = \ERI{ij}{ab} - 4 \sum_m \frac{\sERI{ij}{m} \sERI{ab}{m}}{\OmRPA{m}{} - i \eta}.
\end{equation}
According to perturbation theory, the $s$th BSE excitation energy and its corresponding eigenvector can then expanded as
\begin{subequations}
\begin{gather}
\Om{s}{} = \Om{s}{(0)} + \Om{s}{(1)} + \ldots,
\\
\begin{pmatrix}
\bX{s}{} \\
\bY{s}{} \\
\end{pmatrix}
=
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\end{pmatrix}
+
\begin{pmatrix}
\bX{s}{(1)} \\
\bY{s}{(1)} \\
\end{pmatrix}
+ \ldots.
\end{gather}
\end{subequations}
Solving the zeroth-order static problem
\begin{equation}
\label{eq:LR-BSE-stat}
\begin{pmatrix}
\bA{(0)} & \bB{(0)} \\
-\bB{(0)} & -\bA{(0)} \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\end{pmatrix}
=
\Om{s}{(0)}
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\end{pmatrix},
\end{equation}
yields the zeroth-order (static) $\Om{s}{(0)}$ excitation energies and their corresponding eigenvectors $\bX{s}{(0)}$ and $\bY{s}{(0)}$.
Thanks to first-order perturbation theory, the first-order correction to the $s$th excitation energy is
\begin{equation}
\label{eq:Om1}
\Om{s}{(1)} =
\T{\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\end{pmatrix}}
\cdot
\begin{pmatrix}
\bA{(1)}(\Om{s}{(0)}) & \bB{(1)}(\Om{s}{(0)}) \\
-\bB{(1)}(\titou{-}\Om{s}{(0)}) & -\bA{(1)}(\titou{-}\Om{s}{(0)}) \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\end{pmatrix}.
\end{equation}
From a practical point of view, if one enforces the TDA for the dynamical correction (which we label dTDA in the following), we obtain the very simple expression
\begin{equation}
\label{eq:Om1-TDA}
\Om{s}{(1)} = \T{(\bX{s}{(0)})} \cdot \bA{(1)}(\Om{s}{(0)}) \cdot \bX{s}{(0)}.
\end{equation}
This correction can be renormalized by computing, at basically no extra cost, the renormalization factor which reads, in the dTDA,
\begin{equation}
\label{eq:Z}
Z_{s} = \qty[ 1 - \T{(\bX{s}{(0)})} \cdot \left. \pdv{\bA{(1)}(\Om{s}{})}{\Om{s}{}} \right|_{\Om{s}{} = \Om{s}{(0)}} \cdot \bX{s}{(0)} ]^{-1}.
\end{equation}
This finally yields
\begin{equation}
\Om{s}{\text{dyn}} = \Om{s}{\text{stat}} + \Delta\Om{s}{\text{dyn}} = \Om{s}{(0)} + Z_{s} \Om{s}{(1)}.
\end{equation}
with $\Om{s}{\text{stat}} \equiv \Om{s}{(0)}$ and $\Delta\Om{s}{\text{dyn}} = Z_{s} \Om{s}{(1)}$.
This is our final expression.
%%% FIG 1 %%%
%\begin{figure}
% \includegraphics[width=\linewidth]{}
%\caption{
%\label{fig:}
%}
%\end{figure}
%%% %%% %%%
In terms of computational cost, if one decides to compute the dynamical correction of the $M$ lowest excitation energies, one must perform, first, a conventional (static) BSE calculation and extract the $M$ lowest eigenvalues and their corresponding eigenvectors [see Eq.~\eqref{eq:LR-BSE-stat}].
These are then used to compute the first-order correction from Eq.~\eqref{eq:Om1} or Eq.~\eqref{eq:Om1-TDA}, which also require to construct and evaluate the dynamical part of the BSE Hamiltonian for each excitation one wants to dynamically correct.
The static BSE Hamiltonian is computed once during the static BSE calculation and does not dependent on the targeted excitation.
Searching iteratively for the lowest eigenstates, via Davidson's algorithm for instance, can be performed in $\order*{\Norb^4}$ computational cost.
Constructing the static and dynamic BSE Hamiltonians is much more expensive as it requires the complete diagonalization of the $(\Nocc \Nvir \times \Nocc \Nvir)$ RPA linear response matrix [see Eq.~\eqref{eq:LR-RPA}], which corresponds to a $\order*{\Nocc^3 \Nvir^3} = \order*{\Norb^6}$ computational cost.
Although it might be reduced to $\order*{\Norb^4}$ operations with standard resolution-of-the-identity techniques, \cite{Duchemin_2019,Duchemin_2020} this step is the computational bottleneck in our current implementation.
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:compdet}
%%%%%%%%%%%%%%%%%%%%%%%%
All systems under investigation have close-shell singlet ground states.
We then adopt a restricted formalism throughout this work.
The $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies are done using a (restricted) HF starting point.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986} quasiparticle energies are employed as starting points to compute the BSE neutral excitations.
These quasiparticle energies are obtained by linearizing the frequency-dependent quasiparticle equation, and the entire set of orbitals is corrected.
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018b,Veril_2018}.
As one-electron basis sets, we employ the augmented Dunning family (aug-cc-pVXZ) defined with cartesian Gaussian functions.
Finally, the infinitesimal $\eta$ is set to $100$ meV for all calculations.
For comparison purposes, we employ the theoretical best estimates and geometries of Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} from which coupled cluster (CC) excitation energies, namely, CC2 \cite{Christiansen_1995}, CCSD, \cite{Purvis_1982} and CC3, \cite{Christiansen_1995b} are also extracted.
All the BSE calculations have been performed with our locally developed $GW$ software, \texttt{QuAcK}, \cite{QuAcK} freely available on \texttt{github}, where the present perturbative correction has been implemented.
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and Discussion}
\label{sec:resdis}
%%%%%%%%%%%%%%%%%%%%%%%%
%%% TABLE I %%%
\begin{table*}
\caption{
Singlet and triplet excitation energies (in eV) of \ce{N2} at the BSE@{\GOWO}@HF level for various basis sets.
\label{tab:N2}
}
\begin{ruledtabular}
\begin{tabular}{lddddddddd}
& \mc{3}{c}{aug-cc-pVDZ ($\Eg^{\GW} = 19.49$ eV)}
& \mc{3}{c}{aug-cc-pVTZ ($\Eg^{\GW} = 19.20$ eV)}
& \mc{3}{c}{aug-cc-pVQZ ($\Eg^{\GW} = 19.00$ eV)} \\
\cline{2-4} \cline{5-7} \cline{8-10}
State & \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$}
& \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$}
& \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$} \\
\hline
$^1\Pi_g(n \ra \pis)$ & 10.18 & -0.41 & -0.43 & 10.42 & -0.42 & -0.40 & 10.52 & -0.43 & -0.40 \\
$^1\Sigma_u^-(\pi \ra \pis)$ & 9.95 & -0.44 & -0.44 & 10.11 & -0.45 & -0.45 & 10.20 & -0.45 & -0.45 \\
$^1\Delta_u(\pi \ra \pis)$ & 10.57 & -0.41 & -0.40 & 10.75 & -0.42 & -0.41 & 10.85 & -0.42 & -0.42 \\
$^1\Sigma_g^+$(R) & 13.72 & -0.04 & -0.04 & 13.60 & -0.03 & -0.03 & 13.55 & -0.02 & -0.02 \\
$^1\Pi_u$(R) & 14.07 & -0.05 & -0.05 & 13.98 & -0.04 & -0.04 & 13.96 & -0.03 & -0.04 \\
$^1\Sigma_u^+$(R) & 13.80 & -0.08 & -0.08 & 13.98 & -0.07 & -0.08 & 14.08 & -0.06 & -0.06 \\
$^1\Pi_u$(R) & 14.22 & -0.04 & -0.03 & 14.24 & -0.03 & -0.03 & 14.26 & -0.03 & -0.02 \\
\\
$^3\Sigma_u^+(\pi \ra \pis)$ & 9.21 & -1.01 & -0.71 & 9.50 & -1.04 & -0.73 & 9.61 & -1.04 & -0.71 \\
$^3\Pi_g(n \ra \pis)$ & 9.58 & -0.57 & -0.34 & 9.85 & -0.58 & -0.33 & 9.96 & -0.57 & -0.23 \\
$^3\Delta_u(\pi \ra \pis)$ & 9.97 & -0.92 & -0.58 & 10.19 & -0.95 & -0.36 & 10.29 & -0.95 & -0.70 \\
$^3\Sigma_u^-(\pi \ra \pis)$ & 10.71 & -0.81 & -0.68 & 10.89 & -0.82 & -0.30 & 11.00 & -0.83 & -0.53 \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%%% TABLE I %%%
%\begin{squeezetable}
\begin{table*}
\caption{
Singlet excitation energies (in eV) for various molecules obtained with the aug-cc-pVTZ basis set at various levels of theory.
The dynamical correction is computed in the TDA.
CT and R stand respectively for charge transfer and Rydberg.
\label{tab:BigTabSi}
}
\begin{ruledtabular}
\begin{tabular}{lldddddddddd}
& & \mc{5}{c}{BSE@{\GOWO}@HF} & \mc{5}{c}{Wave function-based methods} \\ %& \mc{5}{c}{Density-based methods} \\
\cline{3-7} \cline{8-12} %\cline{13-17}
Mol. & State & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{s}{\stat}$} & \tabc{$\Om{m}{\dyn}$} & \tabc{$\Delta\Om{s}{\dyn}$} & \tabc{$Z_{s}$}
& \tabc{CIS(D)} & \tabc{ADC(2)} & \tabc{CCSD} & \tabc{CC2} & \tabc{CC3} \\
% & \tabc{B3LYP} & \tabc{PBE0} & \tabc{M06-2X} & \tabc{CAM-B3LYP} & \tabc{LC-$\omega$HPBE} \\
\hline
\ce{HCl} & $^1\Pi$(CT) & 13.43 & 8.30 & 8.19 & -0.11 & 1.009
& 6.07 & 7.97 & 7.91 & 7.96 & 7.84 \\
% & 7.33 & 7.59 & 7.56 & 7.52 & 7.96 \\
\\
\ce{H2O} & $^1B_1(n \ra 3s)$ & 13.58 & 8.09 & 8.00 & -0.09 & 1.007
& 7.62 & 7.18 & 7.60 & 7.23 & 7.65 \\
% & 6.92 & 7.18 & 7.46 & 7.13 & 7.50 \\
& $^1A_2(n \ra 3p)$ & & 9.79 & 9.72 & -0.07 & 1.005
& 9.41 & 8.84 & 9.36 & 8.89 & 9.43 \\
% & 8.33 & 8.61 & 8.93 & 8.69 & 9.11 \\
& $^1A_1(n \ra 3s)$ & & 10.42 & 10.35 & -0.07 & 1.006
& 9.99 & 9.52 & 9.96 & 9.58 & 10.00 \\
% & 9.08 & 9.37 & 9.64 & 9.28 & 9.65 \\
\\
\ce{N2} & $^1\Pi_g(n \ra \pis)$ & 19.20 & 10.42 & 9.99 & -0.42 & 1.031
& 9.66 & 9.48 & 9.41 & 9.44 & 9.34 \\
% & 9.23 & \\
& $^1\Sigma_u^-(\pi \ra \pis)$ & & 10.11 & 9.66 & -0.45 & 1.029
& 10.31 & 10.26 & 10.00 & 10.32 & 9.88 \\
& $^1\Delta_u(\pi \ra \pis)$ & & 10.75 & 10.33 & -0.42 & 1.030
& 10.85 & 10.79 & 10.44 & 10.86 & 10.29 \\
& $^1\Sigma_g^+$(R) & & 13.60 & 13.57 & -0.03 & 1.003
& 13.67 & 12.99 & 13.15 & 12.83 & 13.01 \\
& $^1\Pi_u$(R) & & 13.98 & 13.94 & -0.04 & 1.004
& 13.64 & 13.32 & 13.43 & 13.15 & 13.22 \\
& $^1\Sigma_u^+$(R) & & 13.98 & 13.91 & -0.07 & 1.008
& 13.75 & 13.07 & 13.26 & 12.89 & 13.12 \\
& $^1\Pi_u$(R) & & 14.24 & 14.21 & -0.03 & 1.002
& 14.52 & 14.00 & 13.67 & 13.96 & 13.49 \\
\\
\ce{CO} & $^1\Pi(n \ra \pis)$ & 16.46 & 9.54 & 9.19 & -0.34 & 1.029 & 8.78 & 8.69 & 8.59 & 8.64 & 8.49 \\
& $^1\Sigma^-(\pi \ra \pis)$ & & 10.25 & 9.90 & -0.35 & 1.023 & 10.13 & 10.03 & 9.99 & 10.30 & 9.99 \\
& $^1\Delta(\pi \ra \pis)$ & & 10.71 & 10.39 & -0.32 & 1.023 & 10.41 & 10.30 & 10.12 & 10.60 & 10.12 \\
& $^1\Sigma^+$(R) & & 11.88 & 11.85 & -0.03 & 1.005 & 11.48 & 11.32 & 11.22 & 11.11 & 10.94 \\
& $^1\Sigma^+$(R) & & 12.39 & 12.37 & -0.02 & 1.003 & 11.71 & 11.83 & 11.75 & 11.63 & 11.49 \\
& $^1\Pi$(R) & & 12.37 & 12.32 & -0.05 & 1.004 & 12.06 & 12.03 & 11.96 & 11.83 & 11.69 \\
\\
\ce{HNO} & $^1A''(n \ra \pis)$ & 11.71 & 2.46 & 1.98 & -0.48 & 1.035
& 1.80 & 1.68 & 1.76 & 1.74 & 1.75 \\
% & 1.55 & 1.51 & 0.99 & 1.51 & 1.46 \\
& $^1A'$(R) & & 7.05 & 7.01 & -0.04 & 1.003
& 5.81 & 5.73 & 6.30 & 5.72 & 6.26 \\
% & 5.63 & 5.85 & 6.22 & 5.94 & 6.33 \\
\\
%T2: check state ordering in BSE calculation
\ce{C2H4} & $^1B_{3u}(\pi \ra 3s)$ & 11.49 & 7.64 & 7.62 & -0.03 & 1.004
& 7.35 & 7.34 & 7.42 & 7.29 & 7.35 \\
% & 6.63 & 6.88 & 6.94 & 6.93 & 7.57 \\
& $^1B_{1u}(\pi \ra \pis)$ & & 8.18 & 8.03 & -0.15 & 1.022
& 7.95 & 7.91 & 8.02 & 7.92 & 7.91 \\
% & 8.06 & 7.51 & 7.50 & 7.46 & 7.64 \\
& $^1B_{1g}(\pi \ra 3p)$ & & 8.29 & 8.26 & -0.03 & 1.003
& 8.01 & 7.99 & 8.08 & 7.95 & 8.03 \\
% & 7.18 & 7.45 & 7.47 & 7.54 & 8.15 \\
\\
\ce{CH2O} & $^1A_2(n \ra \pis)$ & 12.00 & 5.03 & 4.68 & -0.35 & 1.027 & 4.04 & 3.92 & 4.01 & 4.07 & 3.97 \\
& $^1B_2(n \ra 3s)$ & & 7.87 & 7.85 & -0.02 & 1.001 & 6.64 & 6.50 & 7.23 & 6.56 & 7.18 \\
& $^1B_2(n \ra 3p)$ & & 8.76 & 8.72 & -0.04 & 1.003 & 7.56 & 7.53 & 8.12 & 7.57 & 8.07 \\
& $^1A_1(n \ra 3p)$ & & 8.85 & 8.84 & -0.01 & 1.000 & 8.16 & 7.47 & 8.21 & 7.52 & 8.18 \\
& $^1A_2(n \ra 3p)$ & & 8.87 & 8.85 & -0.02 & 1.002 & 8.04 & 7.99 & 8.65 & 8.04 & 8.64 \\
& $^1B_1(\si \ra \pis)$ & & 10.18 & 9.77 & -0.42 & 1.032 & 9.38 & 9.17 & 9.28 & 9.32 & 9.19 \\
& $^1A_1(\pi \ra \pis)$ & & 10.05 & 9.81 & -0.24 & 1.026 & 9.08 & 9.46 & 9.67 & 9.54 & 9.48 \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%\end{squeezetable}
%%% TABLE II %%%
%\begin{squeezetable}
\begin{table*}
\caption{
Triplet excitation energies (in eV) for various molecules obtained with the aug-cc-pVTZ basis set at various levels of theory.
The dynamical correction is computed in the TDA.
\label{tab:BigTabTr}
}
\begin{ruledtabular}
\begin{tabular}{lldddddddddd}
& & \mc{5}{c}{BSE@{\GOWO}@HF} & \mc{5}{c}{Wave function-based methods} \\%& \mc{5}{c}{Density-based methods} \\
\cline{3-7} \cline{8-12} %\cline{13-17}
Mol. & State & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{m}{\stat}$} & \tabc{$\Om{m}{\dyn}$} & \tabc{$\Delta\Om{m}{\dyn}$} & \tabc{$Z_{m}$}
& \tabc{CIS(D)} & \tabc{ADC(2)} & \tabc{CCSD} & \tabc{CC2} & \tabc{CC3} \\
% & \tabc{B3LYP} & \tabc{PBE0} & \tabc{M06-2X} & \tabc{CAM-B3LYP} & \tabc{LC-$\omega$HPBE} \\
\hline
\ce{H2O} & $^3B_1(n \ra 3s)$ & 13.58 & 8.14 & 7.98 & -0.15 & 1.014
& 7.25 & 6.86 & 7.20 & 6.91 & 7.28 \\
% & 6.55 & 6.75 & 7.12 & 6.72 & 7.04 \\
& $^3A_2(n \ra 3p)$ & & 9.97 & 9.89 & -0.07 & 1.008
& 9.24 & 8.72 & 9.20 & 8.77 & 9.26 \\
% & 8.22 & 8.45 & 8.77 & 8.54 & 8.92 \\
& $^3A_1(n \ra 3s)$ & & 10.28 & 10.13 & -0.15 & 1.012
& 9.54 & 9.15 & 9.49 & 9.20 & 9.56 \\
% & 8.60 & 8.82 & 9.24 & 8.79 & 9.11 \\
\\
\ce{N2} & $^3\Sigma_u^+(\pi \ra \pis)$ & 19.20 & 9.50 & 8.46 & -1.04 & 1.060 & 8.20 & 8.15 & 7.66 & 8.19 & 7.68 \\
& $^3\Pi_g(n \ra \pis)$ & & 9.85 & 9.27 & -0.58 & 1.050 & 8.33 & 8.20 & 8.09 & 8.19 & 8.04 \\
& $^3\Delta_u(\pi \ra \pis)$ & & 10.19 & 9.24 & -0.95 & 1.060 & 9.30 & 9.25 & 8.91 & 9.30 & 8.87 \\
& $^3\Sigma_u^-(\pi \ra \pis)$ & & 10.89 & 10.06 & -0.82 & 1.058 & 10.29 & 10.23 & 9.83 & 10.29 & 9.68 \\
\\
\ce{CO} & $^3\Pi(n \ra \pis)$ & 16.46 & 8.10 & 7.33 & -0.77 & 1.055 & 6.51 & 6.45 & 6.36 & 6.42 & 6.30 \\
& $^3\Sigma^+(\pi \ra \pis)$ & & 9.61 & 9.04 & -0.57 & 1.037 & 8.63 & 8.54 & 8.34 & 8.72 & 8.45 \\
& $^3\Delta(\pi \ra \pis)$ & & 10.20 & 9.69 & -0.50 & 1.036 & 9.44 & 9.33 & 9.23 & 9.56 & 9.30 \\
& $^3\Sigma_u^-(\pi \ra \pis)$ & & 10.79 & 10.38 & -0.42 & 1.034 & 10.10 & 10.01 & 9.81 & 10.27 & 9.82 \\
& $^3\Sigma_u^+$(R) & & 11.48 & 11.38 & -0.10 & 1.010 & 10.98 & 10.83 & 10.71 & 10.60 & 10.45 \\
\\
\ce{HNO} & $^3A''(n \ra \pis)$ & 11.71 & 3.05 & 2.35 & -0.71 & 1.069
& 0.91 & 0.78 & 0.85 & 0.84 & 0.88 \\
% & -0.47 & -0.61 & 0.36 & -0.49 & -0.58 \\
& $^3A'(\pi \ra \pis)$ & & 6.69 & 6.70 & 0.01 & 1.000
& 5.72 & 5.46 & 5.49 & 5.44 & 5.59 \\
% & 4.73 & 4.46 & 5.27 & 4.55 & 4.57 \\
\\
\ce{C2H4} & $^3B_{1u}(\pi \ra \pis)$ & 11.49 & 6.54 & 5.85 & -0.69 & 1.065
& 4.62 & 4.59 & 4.46 & 4.59 & 4.53 \\
% & 4.07 & 3.84 & 4.54 & 3.92 & 3.55 \\
& $^3B_{3u}(\pi \ra 3s)$ & & 7.61 & 7.55 & -0.06 & 1.008
& 7.26 & 7.23 & 7.29 & 7.19 & 7.24 \\
% & 6.54 & 6.74 & 6.90 & 6.83 & 7.41 \\
& $^3B_{1g}(\pi \ra 3p)$ & & 8.34 & 8.31 & -0.03 & 1.003
& 7.97 & 7.95 & 8.03 & 7.91 & 7.98 \\
% & 7.14 & 7.34 & 7.46 & 7.45 & 7.53 \\
\\
\ce{CH2O} & $^3A_2(n \ra \pis)$ & 12.00 & 5.53 & 5.05 & -0.47 & 1.049 & 3.58 & 3.46 & 3.56 & 3.59 & 3.57 \\
& $^3A_1(\pi \ra \pis)$ & & 8.15 & 7.32 & -0.83 & 1.067 & 6.27 & 6.20 & 5.97 & 6.30 & 6.05 \\
& $^3B_2(n \ra 3s)$ & & 7.51 & 7.54 & 0.03 & 0.994 & 6.66 & 6.39 & 7.08 & 6.44 & 7.03 \\
& $^3B_2(n \ra 3p)$ & & 8.62 & 8.61 & -0.00 & 0.998 & 7.52 & 7.41 & 7.94 & 7.45 & 7.92 \\
& $^3A_1(n \ra 3p)$ & & 8.75 & 8.69 & -0.06 & 1.007 & 7.68 & 7.40 & 8.09 & 7.44 & 8.08 \\
& $^3B_1(n \ra 3d)$ & & 8.82 & 8.82 & -0.01 & 1.000 & 8.57 & 8.39 & 8.43 & 8.52 & 8.41 \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%\end{squeezetable}
%%% TABLE III %%%
%\begin{table}
% \caption{
% Excitation energies (in eV) of CN3 obtained with the aug-cc-pVDZ basis set at various levels of theory.
% %$\Eg^{\GW} = 13.79$ eV.
% \label{tab:CN3}
% }
% \begin{ruledtabular}
% \begin{tabular}{lcc}
% & \mc{2}{c}{Excitation} \\
% Method & $^1B_2(\pi \ra \pis)$ & $^3B_2(\pi \ra \pis)$ \\
% \hline
% BSE@{\GOWO}@HF & 7.66 & 6.52 \\
% dBSE(TDA)@{\GOWO}@HF & 7.51 & 6.11 \\
% FCI & 7.14 & 5.47 \\
% \end{tabular}
% \end{ruledtabular}
%\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%
This is the conclusion
%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
%%%%%%%%%%%%%%%%%%%%%%%%
PFL thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481) for financial support.
This work was performed using HPC resources from GENCI-TGCC (Grant No.~2019-A0060801738) and CALMIP (Toulouse) under allocation 2020-18005.
Funding from the \textit{``Centre National de la Recherche Scientifique''} is acknowledged.
This work has also been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''.}}
%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Supporting Information}
%%%%%%%%%%%%%%%%%%%%%%%%
%See {\SI} for plenty of stuff
\bibliography{BSEdyn}
\end{document}