\documentclass[aps,prb,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}[1]{\mathbf{r}_{#1}}
\newcommand{\dbr}[1]{d\br{#1}}

% 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{\Nel}{N}
\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]{\epsilon_{#1}}
\newcommand{\eHF}[1]{\epsilon^\text{HF}_{#1}}
\newcommand{\eKS}[1]{\epsilon^\text{KS}_{#1}}
\newcommand{\eQP}[1]{\epsilon^\text{QP}_{#1}}
\newcommand{\eGOWO}[1]{\epsilon^\text{\GOWO}_{#1}}
\newcommand{\eGW}[1]{\epsilon^{GW}_{#1}}
\newcommand{\eevGW}[1]{\epsilon^\text{\evGW}_{#1}}
\newcommand{\eGnWn}[2]{\epsilon^\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{\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\vari{{\varepsilon}_i}
\newcommand\vara{{\varepsilon}_a}
\newcommand\varj{{\varepsilon}_j}
\newcommand\varb{{\varepsilon}_b}
\newcommand\varn{{\varepsilon}_n}
\newcommand\varm{{\varepsilon}_m}
\newcommand\Oms{{\Omega}_s}
\newcommand\hOms{\frac{{\Omega}_s}{2}}

\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{Dynamical Correction 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}
Similarly 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 molecular excitation energies thanks to a renormalized first-order perturbative correction to the static BSE excitation energies.
The present 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 theoretical chemists with a large number of systematic benchmark studies on large molecular systems appearing in the scientific 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} - \varepsilon_{\HOMO}^{\GW},
\end{equation}
which is itself a corrected version of the Kohn-Sham (KS) gap
\begin{equation}
	\Eg^{\KS} = \eps_{\LUMO}^{\KS} - \varepsilon_{\HOMO}^{\KS} \ll \Eg^{\GW} \approx \EgFun,
\end{equation}
in order to approximate the optical gap
\begin{equation}
	\EgOpt = E_1^{\Nel} - E_0^{\Nel} = \EgFun + \EB,
\end{equation}
where $\EgFun = I^\Nel - A^\Nel$ is the the fundamental gap, \cite{Bredas_2014} $I^\Nel = E_0^{\Nel-1} - E_0^\Nel$ and $A^\Nel = E_0^{\Nel+1} - E_0^\Nel$ being the ionization potential and the electron affinity of the $\Nel$-electron system.
Here, $E_s^{\Nel}$ is the total energy of the $s$th excited state of the $\Nel$-electron system, and $E_0^\Nel$ 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, one key consequence of the static approximation is that double (and higher) excitations are completely absent from the BSE spectrum.
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} 
They are, moreover, a real challenge for high-level computational methods. \cite{Loos_2018a,Loos_2019,Loos_2020b}

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 as 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 double 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.
Finally, let us also mentioned the work of Romaniello and coworkers, \cite{Romaniello_2009b,Sangalli_2011} in which the authors genuinely accessed additional excitations by solving the non-linear, frequency-dependent eigenvalue problem. 
However, it is based on a rather simple model (the Hubbard dimer) which permits to analytically solve the dynamical equations.

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, we describe the theoretical foundations leading to the dynamical Bethe-Salpeter equation, following the seminal work by Strinati, \cite{Strinati_1988} presenting in a second step the perturbative implementation \cite{Rohlfing_2000,Ma_2009a,Ma_2009b} of the dynamical correction as compared to the standard adiabatic calculations. More details of the derivation are provided in ...
%================================
\subsection{General dynamical BSE theory}
%=================================

The resolution \cite{Strinati_1988} of the  Bethe-Salpeter equation:
\begin{align*}
    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{align*}
 with:
\begin{align*}
 iL(1,2;   1',2') &= -G_2(1,2;1',2') + G(1,1')G(2,2')    \\
 i^2 G_2(1,2;1',2') &= \langle N | T {\hat \psi}(1) {\hat \psi}(2) {\hat \psi}^{\dagger}(2') {\hat \psi}^{\dagger}(1') | N \rangle
\end{align*}
where e.g. $1 = (x_1,t_1)$ a space-spin plus time variable, starts with the expansion of the 2-body Green's function $G_2$ and response function $L$ over the complete orthonormalized set  $  |N,s \rangle   $   of the N-electron system excited states, with $| N \rangle = | N,0 \rangle$ the ground-state. In the optical limit of instantaneous electron-hole creation and destruction, imposing 
$t_{2'} = t_2^+$ and $t_{1'} = t_1^+$, with e.g. $t_2^+ = t_2 + 0^+$ where $0^+$ is a small positive infinitesimal, one obtains:
\begin{align*}
 iL(1,2;1',2') &= \theta(\tau_{12}) \sum_{s > 0} \chi_s(x_1,x_{1'}) {\tilde \chi}_s(x_2,x_{2'})
 e^{ -i \Oms  \tau_{12}  } \\
 &- \theta(-\tau_{12}) \sum_{s > 0} \chi_s(x_2,x_{2'}) {\tilde \chi}_s(x_1,x_{1'})
 e^{ + i \Oms \tau_{12} }  
\end{align*}
     with $\tau_{12} = t_1 - t_2$ and
\begin{align*}
 \chi_s(x_1,x_{1'}) = \langle N | T {\hat \psi}(x_1) {\hat \psi}^{\dagger}(x_{1'}) | N,s \rangle \\
 {\tilde \chi}_s(x_2,x_{2'}) = \langle N,s | T {\hat \psi}(x_2) {\hat \psi}^{\dagger}(x_{2'}) | N \rangle
     \end{align*}  
The $\Oms$ are the neutral excitation energies of interest. Picking up the $e^{+i \Oms t_2 }$ component in $L(1,2;   1',2')$ and $L(6,2;5,2')$,   simplifying further by ${\tilde \chi}_s(x_2,x_{2'})$ on both side of the Bethe-Salpeter equation, we are left with the  search of the  $e^{-i \Oms t_1 }$ Fourier component  associated with the right-hand side of the  BSE.  For the lowest $\Oms$ excitation energies falling in the quasiparticle gap of the system, the  $L_0(1,2;1',2')$ term cannot contribute since its lowest excitation energy is precisely the quasiparticle gap, namely the difference between the electronic affinity and the ionisation potential.
 
 The Fourier components with respect to time $t_1$ of $iL_0(1, 4; 1', 3) = G(1, 3)G(4, 1')$   reads, dropping the (space/spin)-variables:  
\begin{align*}
[iL_0]( \omega_1 )  =  \int \frac{ d \omega }{  2\pi }   \; G(\omega - \frac{\omega_1}{2} ) G( {\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$. Plugging now the 1-body Green's function Lehman representation, e.g.
$$
G(x_1,x_3 ; \omega) = \sum_n \frac{ \phi_n(x_1) \phi_n^*(x_3) } { \omega - \varepsilon_n + i \eta \text{sgn} \times (\varepsilon_n - \mu) }
$$
and projecting on $\phi_a^*(x_1) \phi_i(x_{1'})$, one obtains the $\omega_1= \Oms$ component  
\begin{align*}
\int dx_1 dx_{1'} \; & \phi_a^*(x_1) \phi_i(x_{1'})  L_0(x_1,3;x_{1'},4; \Oms) = e^{i \Oms t^{34} } \times    \\
     &       \frac{ \phi_a^*(x_3) \phi_i(x_4)  } {   \Oms - (  \vara - \vari )  + i \eta }   
       \left[ \theta( \tau ) e^{i ( \vari + \hOms) \tau }   
     +   \theta( - \tau  )   e^{i (\vara -  \hOms \tau }  \right]  
\end{align*}
with $\tau = \tau_{34}$. Adopting now the $GW$ approximation for the exchange-correlation self-energy leads to a simplification of the BSE kernel:
$$
\Xi(3,5;4,6) = v(3,6) \delta(34) \delta(56) - W(3^+,4)   \delta(36) \delta(45)
$$
We further obtain the needed spectral representation of 
$\langle N  | T {\hat \psi}(3) {\hat \psi}^{\dagger}(4) | N,s \rangle$
expanding the field operators  over a complete orbital basis creation/destruction operators:
\begin{align*}
\langle N  | T   {\hat \psi}(3)  &   {\hat \psi}^{\dagger}(4) | N,s \rangle = - \Big( e^{  -i \Omega_s t^{34} }  \Big)    \sum_{mn} \phi_m(x_3) \phi_n^*(x_4)    \times \nonumber \\
\times & \langle N | {\hat a}_n^{\dagger} {\hat a}_m  | N,s \rangle  \;\Big[  \theta( \tau )   e^{- i (  \varepsilon_m  -  \hOms ) \tau  }   + \theta( -\tau )   e^{ - i   ( \varepsilon_n  + \hOms)  \tau  }      \Big]
\end{align*}
with $\tau = \tau_{34}$ and where the $ \lbrace \varepsilon_{n/m} \rbrace$ are proper addition/removal energies such that e.g.
$$
e^{  i H \tau  } {\hat a}_m^{\dagger}   | N  \rangle = e^{ i (E_0^N + \varepsilon_m ) \tau  }  {\hat a} _m^{\dagger}   | N  \rangle
$$
The $GW$ quasiparticle energies $    \varepsilon_{n/m}^{GW}$ are good approximations to such removal/addition energies. Selecting (n,m)=(j,b) yields the largest components
$A_{jb}^{s} = \langle N | {\hat a}_j^{\dagger} {\hat a}_b  | N,s \rangle $, while  (n,m)=(b,j) yields much weaker  
$B_{jb}^{s} = \langle N | {\hat a}_b^{\dagger} {\hat a}_j  | N,s \rangle $ contributions. We used chemist notations with (i,j) indexing occupied orbitals and (a,b) virtual ones. Neglecting the   $B_{jb}^{s}$ weights leads to the Tamm Dancoff approximation (TDA). Working out the same expansion for $ \langle N  | T {\hat \psi}(5) {\hat \psi}^{\dagger}(5) | N,s \rangle$ and $ \langle N  | T {\hat \psi}(x_1) {\hat \psi}^{\dagger}(x_{1'}) | N,s \rangle$,  and projecting onto $\phi_a^*(x_1) \phi_i(x_{1'})$, 
one obtains after a few tedious manipulations (see Supplemental Information) the dynamical  Bethe-Salpeter equation (dBSE) :
 \begin{align}
      ( \varepsilon_a - \varepsilon_i    - \Omega_s    )       A_{ia}^{s}         
       &+ \sum_{jb}  \Big(  v_{ai,bj} - \widetilde{W}_{ij,ab}(\Oms)    \Big)  A_{jb}^{s}  \\ 
       &+  \sum_{bj}  \Big(  v_{ai,jb} - \widetilde{W}_{ib,aj}(\Oms)    \Big)  B_{jb}^{s} 
       = 0
\end{align}
 with an effective dynamically  screened Coulomb potential (see Pina  eq. 24): 
 \begin{align}
 \widetilde{W}_{ij,ab}(\Oms)  &= \frac{ i }{ 2  \pi}  \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega) \times \\
    \hskip 1cm &\times \left[ \frac{1}{  \Omega_{ib}^s - \omega   +i \eta } + \frac{1}{  \Omega_{ja}^{s} + \omega   + i\eta } \right] \nonumber
\end{align}
where $\; \Omega_{ib}^s =  \Oms - ( \varb - \vari )$ and $\; \Omega_{ja}^s = \Oms - ( \vara - \varj ).$  
In the present study, we use the exact spectral representation of $W(\omega)$ at the RPA level:
\begin{align*}
W_{ij,ab}(\omega) &= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m] \times \\
   & \times \Big(  \frac{1}{ \omega-\Omega_m^{RPA} + i\eta }  -   \frac{1}{ \omega + \Omega_m^{RPA} - i\eta } \Big) 
\end{align*}
($\Omega_m^{RPA} > 0 $) so that
\begin{align}
 \widetilde{W}_{ij,ab}( \Oms ) &=   (ij|ab) + 2 \sum_m^{OV}    [ij|m] [ab|m] \times \\ 
    & \times     \left[ \frac{  1 }{    \Omega_{ib}^{s} - \Omega_m^{RPA}    + i\eta   }   +   \frac{ 1}{   \Omega_{ja}^{s}   - \Omega_m^{RPA}   + i\eta }
     \right]   \nonumber
\end{align} 
  \textcolor{red}{Due to excitonic effects, the lowest BSE ${\Omega}_1$ excitation energy stands lower than the lowest $\Omega_m^{RPA}$ excitation  energy, so that 
e.g. $( \Omega_{ib}^{s} - \Omega_m^{RPA} )$ is strictly negative and $\widetilde{W}_{ij,ab}( \Oms )$ presents no resonances. Further, $\Omega_{ib}^{s}$ and $\Omega_{ja}^{s}$ are necessarily negative for in-gap   low lying BSE excitations, such that e.g.
$$
\left| \frac{  1 }{    \Omega_{ib}^{s} - \Omega_m^{RPA}       } \right|
 <  \frac{1}{ \Omega_m^{RPA}  }
$$
This leads to reduced electron-hole screening, namely larger electron-hole stabilising binding energy, as compared to the standard adiabatic BSE, leading to smaller (red-shifted) excitation energies. }


%In order to compute the neutral (optical) excitations of the system and their associated oscillator strengths, the BSE expresses the two-body propagator \cite{Strinati_1988, Martin_2016} 
%\begin{multline}
%\label{eq:BSE}
%	\LBSE{}(1,2,1',2') = \LBSE{0}(1,2,1',2') 
%	\\
%	+ \int d3 d4 d5  d6  \LBSE{0}(1,4,1',3) \XiBSE{}(3,5,4,6) \LBSE{}(6,2,5,2')
%\end{multline}
%as the linear response of the one-body Green's function $\G{}$ with respect to a general non-local external potential
%\begin{equation}
%	\XiBSE{}(3,5,4,6) = i \fdv{[\vc{\Ha}(3) \delta(3,4) + \Sig{\xc}(3,4)]}{\G{}(6,5)},
%\end{equation}
%which takes into account the self-consistent variation of the Hartree potential 
%\begin{equation}
%	\vc{\Ha}(1) = - i \int d2 \vc{}(2) \G{}(2,2^+),
%\end{equation}
%(where $\vc{}$ is the bare Coulomb operator) and the exchange-correlation self-energy $\Sig{\xc}$.
%In Eq.~\eqref{eq:BSE}, $\LBSE{0}(1,2,1',2') = - i \G{}(1,2')\G{}(2,1')$, and $(1)=(\br{1},\sigma_1,t_1)$ is a composite index gathering space, spin and time variables.
%In the $GW$ approximation, \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Martin_2016,Reining_2017} we have
%\begin{equation}
%	\SigGW{\xc}(1,2) = i \G{}(1,2) \W{}{}(1^+,2),
%\end{equation}
%where $\W{}{}$ is the screened Coulomb operator, and hence the BSE reduces to 
%\begin{equation}
%	\XiBSE{}(3,5,4,6) = \delta(3,4) \delta(5,6) \vc{}(3,6) - \delta(3,6) \delta(4,5) \W{}{}(3,4),
%\end{equation}
%where, as commonly done, we have neglected the term $\delta \W{}{}/\delta \G{}$ in the functional derivative of the self-energy. \cite{Hanke_1980,Strinati_1984,Strinati_1982}
%Finally, the static approximation is enforced, \ie, $\W{}{}(1,2) = \W{}{}(\{\br{1}, \sigma_1, t_1\},\{\br{2}, \sigma_2, t_2\}) \delta(t_1-t_2)$, which corresponds to restricting $\W{}{}$ to its static limit, \ie, $\W{}{}(1,2) = \W{}{}(\{\br{1},\sigma_1\},\{\br{2},\sigma_2\}; \omega=0)$.

%================================
\subsection{Perturbative dynamical correction}
%=================================

For a closed-shell system in a finite basis, to compute the BSE excitation energies, 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{}(\omega)	&	\bB{}(\omega)	\\
		-\bB{}(\omega)	&	-\bA{}(\omega)	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{}(\omega)	\\
		\bY{m}{}(\omega)	\\
	\end{pmatrix}
	=
	\omega
	\begin{pmatrix}
		\bX{m}{}(\omega)	\\
		\bY{m}{}(\omega)	\\
	\end{pmatrix},
\end{equation}
where the dynamical matrices $\bA{}(\omega)$ and $\bB{}(\omega)$ 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}{}(\omega)$, and $\bY{m}{}(\omega)$ are (eigen)vectors of length $\Nocc \Nvir$. 
In the following, the index $m$ labels the $\Nocc \Nvir$ single excitations, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals.
Note that, due to its non-linear nature, Eq.~\eqref{eq:LR-dyn} may provide more than one solution for each value of $m$. \cite{Romaniello_2009b,Sangalli_2011,Martin_2016}

The BSE matrix elements read
\begin{subequations}
\begin{align}
	\label{eq:BSE-Adyn}
	\A{ia,jb}{}(\omega) & = \delta_{ij} \delta_{ab} \eGW{ia} + 2 \sigma \ERI{ia}{jb} - \W{ij,ab}{}(\omega),
	\\ 
	\label{eq:BSE-Bdyn}
	\B{ia,jb}{}(\omega) & =  2 \sigma \ERI{ia}{bj} -   \W{ib,aj}{}(\omega),
\end{align}
\end{subequations}
where $\eGW{ia} = \eGW{a} - \eGW{i}$ are occupied-to-virtual differences of $GW$ quasiparticle energies, $\sigma = 1 $ or $0$ for singlet and triplet excited states (respectively), 
\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}
where $\eta$ is a positive infinitesimal, and
\begin{equation}
\label{eq:sERI}
	\sERI{pq}{m} =  \sum_i^{\Nocc} \sum_a^{\Nvir} \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 direct (\ie, without exchange) 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$ or as the $GW$ quasiparticle energies in the case of self-consistent schemes such as ev$GW$.

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{equation}
\label{eq:LR-PT}
	\begin{pmatrix}
		\bA{}(\omega)	&	\bB{}(\omega)	\\
		-\bB{}(\omega)	&	-\bA{}(\omega)	\\
	\end{pmatrix}
	=
	\begin{pmatrix}
		\bA{(0)}	&	\bB{(0)}	\\
		-\bB{(0)}	&	-\bA{(0)}	\\
	\end{pmatrix}
	+
	\begin{pmatrix}
		\bA{(1)}(\omega)	&	\bB{(1)}(\omega)	\\
		-\bB{(1)}(\omega)	&	-\bA{(1)}(\omega)	\\
	\end{pmatrix}
\end{equation}
with
\begin{subequations}
\begin{align}
	\label{eq:BSE-A0}
	\A{ia,jb}{(0)} & = \delta_{ij} \delta_{ab} \eGW{ia} + 2 \ERI{ia}{jb} - \W{ij,ab}{\text{stat}},
	\\ 
	\label{eq:BSE-B0}
	\B{ia,jb}{(0)} & =  2 \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)}(\omega) & = - \W{ij,ab}{}(\omega) + \W{ij,ab}{\text{stat}},
	\\ 
	\label{eq:BSE-B1}
	\B{ia,jb}{(1)}(\omega) & = -  \W{ib,aj}{}(\omega) + \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}} = \ERI{ij}{ab} -  4 \sum_m^{\Nocc \Nvir} \frac{\sERI{ij}{m} \sERI{ab}{m}}{\OmRPA{m}{} - i \eta}.
\end{equation}
According to perturbation theory, the $m$th BSE excitation energy and its corresponding eigenvector can then decomposed as
\begin{subequations}
\begin{gather}
	\Om{m}{} = \Om{m}{(0)} + \Om{m}{(1)} + \ldots,
	\\
	\begin{pmatrix}
		\bX{m}{}	\\
		\bY{m}{}	\\
	\end{pmatrix}
	= 
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}
	+
	\begin{pmatrix}
		\bX{m}{(1)}	\\
		\bY{m}{(1)}	\\
	\end{pmatrix}
	+ \ldots.
\end{gather}
\end{subequations}
Solving the zeroth-order static problem yields
\begin{equation}
\label{eq:LR-BSE-stat}
	\begin{pmatrix}
		\bA{(0)}	&	\bB{(0)}	\\
		-\bB{(0)}	&	-\bA{(0)}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}
	=
	\Om{m}{(0)}
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix},
\end{equation}
and, thanks to first-order perturbation theory, the first-order correction to the $m$th excitation energy is
\begin{equation}
\label{eq:Om1}
	\Om{m}{(1)} =
	\T{\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}}
	\cdot
	\begin{pmatrix}
		\bA{(1)}(\Om{m}{(0)})	&	\bB{(1)}(\Om{m}{(0)})	\\
		-\bB{(1)}(\Om{m}{(0)})	&	-\bA{(1)}(\Om{m}{(0)})	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}.
\end{equation}
From a practical point of view, if one enforces the TDA, we obtain the very simple expression
\begin{equation}
\label{eq:Om1-TDA}
	\Om{m}{(1)} = \T{(\bX{m}{(0)})} \cdot \bA{(1)}(\Om{m}{(0)}) \cdot \bX{m}{(0)}.
\end{equation}
This correction can be renormalized by computing, at basically no extra cost, the renormalization factor which reads, in the TDA,
\begin{equation}
\label{eq:Z}
	Z_{m} = \qty[ 1 - \T{(\bX{m}{(0)})} \cdot \left. \pdv{\bA{(1)}(\omega)}{\omega} \right|_{\omega = \Om{m}{(0)}} \cdot \bX{m}{(0)} ]^{-1}.
\end{equation}
This finally yields
\begin{equation}
	\Om{m}{\text{dyn}} = \Om{m}{\text{stat}} + \Delta\Om{m}{\text{dyn}} = \Om{m}{(0)} + Z_{m} \Om{m}{(1)}.
\end{equation}
with $\Om{m}{\text{stat}} \equiv \Om{m}{(0)}$ and $\Delta\Om{m}{\text{dyn}} = Z_{m} \Om{m}{(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 the Davidson 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} and partially self-consistent ev$GW$ \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011}calculations are employed as starting points to compute the BSE neutral excitations.
For both {\GOWO} and {\evGW}, the entire set of orbitals are corrected.
In the case of {\GOWO}, the quasiparticle energies are obtained by linearizing the frequency-dependent quasiparticle equation.
For ev$GW$, the convergence criterion has been set to $10^{-5}$.
Further details about our implementation of {\GOWO} and {\evGW} 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{m}{\stat}$}	&	\tabc{$\Delta\Om{m}{\dyn}$(TDA)}	&	\tabc{$\Delta\Om{m}{\dyn}$} 	
											&	\tabc{$\Om{m}{\stat}$}	&	\tabc{$\Delta\Om{m}{\dyn}$(TDA)}	&	\tabc{$\Delta\Om{m}{\dyn}$}	
											&	\tabc{$\Om{m}{\stat}$}	&	\tabc{$\Delta\Om{m}{\dyn}$(TDA)}	&	\tabc{$\Delta\Om{m}{\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{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.
	\label{tab:BigTabSi}
	}
	\begin{ruledtabular}
		\begin{tabular}{llddddddddddddddd}
						&			&	\mc{5}{c}{BSE@{\GOWO}@HF}	&	\mc{5}{c}{BSE@{\evGW}@HF}	\\
													\cline{3-7}	\cline{8-12}
			Mol.		&	State	&	\tabc{$\Eg^{\GW}$}	&	\tabc{$\Om{m}{\stat}$}	&	\tabc{$\Om{m}{\dyn}$}	&	\tabc{$\Delta\Om{m}{\dyn}$}	&	\tabc{$Z_{m}$}	
									&	\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}	\\
			\hline
			\ce{HCl}	&	$^1\Pi$(CT)			&	13.43	&	8.30	&	8.19	&	-0.11	&	1.009	&	13.42	&	8.29	&	8.18	&	-0.11	&	1.010	&	6.07	&	7.97	&	7.91	&	7.96	&		7.84	\\
			\\
			\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	\\
						&	$^1A_2(n \ra 3p)$	&			&	9.79	&	9.72	&	-0.07	&	1.005	&	&	&	&	&	&	9.41	&	8.84	&	9.36	&	8.89	&	9.43	\\
						&	$^1A_1(n \ra 3s)$	&			&	10.42	&	10.35	&	-0.07	&	1.006	&	&	&	&	&	&	9.99	&	9.52	&	9.96	&	9.58	&	10.00	\\
			\\
			\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	\\
						&	$^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{C2H4}	&	$^1B_{3u}(\pi \ra 3s)$		&	&	&	&	&	&	&	&	&	&	&	7.35	&	7.34	&	7.42	&	7.29	&	7.35	\\
						&	$^1B_{1u}(\pi \ra \pis)$	&	&	&	&	&	&	&	&	&	&	&	7.95	&	7.91	&	8.02	&	7.92	&	7.91	\\
						&	$^1B_{1g}(\pi \ra 3p)$		&	&	&	&	&	&	&	&	&	&	&	8.01	&	7.99	&	8.08	&	7.95	&	8.03	\\
			\\
			\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*}

%%% TABLE II %%%
\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}{llddddddddddddddd}
						&			&	\mc{5}{c}{BSE@{\GOWO}@HF}	&	\mc{5}{c}{BSE@{\evGW}@HF}	\\
													\cline{3-7}	\cline{8-12}
			Mol.		&	State	&	\tabc{$\Eg^{\GW}$}	&	\tabc{$\Om{m}{\stat}$}	&	\tabc{$\Om{m}{\dyn}$}	&	\tabc{$\Delta\Om{m}{\dyn}$}	&	\tabc{$Z_{m}$}	
									&	\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}	\\
			\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	\\
						&	$^3A_2(n \ra 3p)$	&			&	9.97	&	9.89	&	-0.07	&	1.008	&	&	&	&	&	&	9.24	&	8.72	&	9.20	&	8.77	&	9.26	\\
						&	$^3A_1(n \ra 3s)$	&			&	10.28	&	10.13	&	-0.15	&	1.012	&	&	&	&	&	&	9.54	&	9.15	&	9.49	&	9.20	&	9.56	\\
			\\
			\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{C2H4}	&	$^3B_{1u}(\pi \ra \pis)$	&	&	&	&	&	&	&	&	&	&	&	4.62	&	4.59	&	4.46	&	4.59	&	4.53	\\
						&	$^3B_{3u}(\pi \ra 3s)$		&	&	&	&	&	&	&	&	&	&	&	7.26	&	7.23	&	7.29	&	7.19	&	7.24	\\
						&	$^3B_{1g}(\pi \ra 3p)$		&	&	&	&	&	&	&	&	&	&	&	7.97	&	7.95	&	8.03	&	7.91	&	7.98	\\
			\\
			\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*}

%%% 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}
