\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{\denis}[1]{\textcolor{purple}{#1}}
\newcommand{\xavier}[1]{\textcolor{darkgreen}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashDJ}[1]{\textcolor{purple}{\sout{#1}}}
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\renewcommand{\DJ}[1]{\denis{(\underline{\bf DJ}: #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{\Norb}{N}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}
\newcommand{\IS}{\lambda}

% operators
\newcommand{\hH}{\Hat{H}}

% 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{\EcRPAx}{E_\text{c}^\text{RPAx}}
\newcommand{\EcBSE}{E_\text{c}^\text{BSE}}
\newcommand{\IP}{\text{IP}}
\newcommand{\EA}{\text{EA}}
\newcommand{\Req}{R_\text{eq}}

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

% 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}[3]{[#1|#2]^{#3}}

%% 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}[2]{\Omega_{#1}^{#2,\text{RPA}}}
\newcommand{\OmRPAx}[2]{\Omega_{#1}^{#2,\text{RPAx}}}
\newcommand{\OmBSE}[2]{\Omega_{#1}^{#2,\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}[1]{\mathbf{X}^{#1}}
\newcommand{\bY}[1]{\mathbf{Y}^{#1}}
\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}


\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{\CEISAM}{Laboratoire CEISAM - UMR CNRS 6230, Universit\'e de Nantes, 2 Rue de la Houssini\`ere, BP 92208, 44322 Nantes Cedex 3, France}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\CEA}{Universit\'e Grenoble Alpes, CEA, IRIG-MEM-L Sim, 38054 Grenoble, 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}
This is the abstract
%\\
%\bigskip
%\begin{center}
%	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
%\end{center}
%\bigskip
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theory}
%%%%%%%%%%%%%%%%%%%%%%%%

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 )  = \frac{ 1 }{  2\pi } \int d \omega   \; 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}(\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 }   
       \Big( \theta( \tau ) e^{i ( \vari + \hOms) \tau }   
     +   \theta( - \tau  )   e^{i (\vara -  \hOms \tau }  \Big)  
\end{align*}
with $\tau = \tau_{34}$.
We further obtain the 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)  \langle N | {\hat a}_n^{\dagger} {\hat a}_m  | N,s \rangle    \times \nonumber \\
\times & \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
$$
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}$ leads to the Tamm Dancoff approximation (TDA). Obtaining similarly the spectral representation of $ \langle N  | T {\hat \psi}(1) {\hat \psi}^{\dagger}(1') | N,s \rangle$ ($t_{1'} = t_1^{+}$) projected 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)  &= { i \over 2  \pi}  \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega) \times \\
    \hskip 1cm &\times \left[ \frac{1}{ (\Oms - \omega) - ( \varb - \vari ) +i \eta } + \frac{1}{ (\Oms + \omega) - ( \vara - \varj ) + i\eta } \right] \nonumber
\end{align}
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*}
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} 
with e.g.  $  \Omega_{ib}^{s} = \Oms - ( \varepsilon_b -  \varepsilon_i) $.   \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 cannot diverge. Further, $\Omega_{ib}^{s}$ and $\Omega_{ja}^{s}$ are necessarily negative for in-gap   low lying BSE excitations, such that
$$
\left[ \frac{  1 }{    \Omega_{ib}^{s} - \Omega_m^{RPA}    + i\eta   }   +   \frac{ 1}{   \Omega_{ja}^{s}   - \Omega_m^{RPA}   + i\eta }
     \right] 
     < 
  \Big(  \frac{1}{ \omega-\Omega_m^{RPA} + i\eta }  -   \frac{1}{ \omega + \Omega_m^{RPA} - i\eta } \Big) < 0   
$$
in the limit $(\omega \rightarrow 0)$ of the standard adiabatic BSE . WELL, do we know the sign of 
$[ij|m] [ab|m]$ ?? }

\titou{This is the theory section from the previous paper.}

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)$.

For a closed-shell system in a finite basis, to compute the singlet BSE excitation energies (within the static approximation) of the physical system (\ie, $\IS = 1$), one must solve the following linear response problem \cite{Casida,Dreuw_2005,Martin_2016} 
\begin{equation}
\label{eq:LR}
	\begin{pmatrix}
		\bA{\IS}	&	\bB{\IS}	\\
		-\bB{\IS}	&	-\bA{\IS}	\\
	\end{pmatrix}
	\begin{pmatrix}
		\bX{\IS}_m	\\
		\bY{\IS}_m	\\
	\end{pmatrix}
	=
	\Om{m}{\IS}
	\begin{pmatrix}
		\bX{\IS}_m	\\
		\bY{\IS}_m	\\
	\end{pmatrix},
\end{equation}
where $\Om{m}{\IS}$ is the $m$th excitation energy with eigenvector $\T{(\bX{\IS}_m \, \bY{\IS}_m)}$ at interaction strength $\IS$, $\T{}$ is the matrix transpose, and we assume real-valued spatial orbitals $\{\MO{p}(\br{})\}_{1 \le p \le \Norb}$.
The matrices $\bA{\IS}$, $\bB{\IS}$, $\bX{\IS}$, and $\bY{\IS}$ are all 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.
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.

In the absence of instabilities (\ie, when $\bA{\IS} - \bB{\IS}$ is positive-definite), \cite{Dreuw_2005} Eq.~\eqref{eq:LR} is usually transformed into an Hermitian eigenvalue problem of smaller dimension
\begin{equation}
\label{eq:small-LR}
	(\bA{\IS} - \bB{\IS})^{1/2} (\bA{\IS} + \bB{\IS}) (\bA{\IS} - \bB{\IS})^{1/2} \bZ{m}{\IS} = (\Om{m}{\IS})^2 \bZ{m}{\IS},
\end{equation}
where the excitation amplitudes are
\begin{subequations}
\begin{align}
	(\bX{\IS} + \bY{\IS})_m = (\Om{m}{\IS})^{-1/2} (\bA{\IS} - \bB{\IS})^{+1/2} \bZ{m}{\IS},
	\\
	(\bX{\IS} - \bY{\IS})_m = (\Om{m}{\IS})^{+1/2} (\bA{\IS} - \bB{\IS})^{-1/2} \bZ{m}{\IS}.
\end{align}
\end{subequations}
Introducing the so-called Mulliken notation for the bare two-electron integrals
\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}
and the corresponding (static) screened Coulomb potential matrix elements at coupling strength $\IS$ 
\begin{equation}
	\W{pq,rs}{\IS} = \iint  \MO{p}(\br{}) \MO{q}(\br{}) \W{}{\IS}(\br{},\br{}') \MO{r}(\br{}') \MO{s}(\br{}')  \dbr{} \dbr{}',
\end{equation}
the BSE matrix elements read
\begin{subequations}
\begin{align}
	\label{eq:LR_BSE-A}
	\ABSE{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \IS \qty[ 2 \ERI{ia}{jb}   -  \W{ij,ab}{\IS} ],
	\\ 
	\label{eq:LR_BSE-B}
	\BBSE{ia,jb}{\IS} & =  \IS \qty[ 2 \ERI{ia}{bj} -   \W{ib,aj}{\IS} ],
\end{align}
\end{subequations}
where $\eGW{p}$ are the $GW$ quasiparticle energies. 
In the standard BSE approach, $\W{}{\IS}$ is built within the direct RPA scheme, \ie,
\begin{subequations}
\label{eq:wrpa}
\begin{align}
	\W{}{\IS}(\br{},\br{}') 
	& = \int \frac{\epsilon_{\IS}^{-1}(\br{},\br{}''; \omega=0)}{\abs*{\br{}' - \br{}''}} \dbr{}'' , 
	\\ 
	\epsilon_{\IS}(\br{},\br{}'; \omega) 
	& = \delta(\br{}-\br{}') - \IS  \int \frac{\chi_{0}(\br{},\br{}''; \omega)}{\abs*{\br{}' - \br{}''}} \dbr{}'' ,
\end{align}
\end{subequations}
with $\epsilon_{\IS}$ the dielectric function at coupling constant $\IS$ and $\chi_{0}$ the non-interacting polarizability.  In the occupied-to-virtual orbital product basis, the spectral representation of $\W{}{\IS}$ can be written as follows in the case of real spatial orbitals
\begin{multline}
\label{eq:W}
	\W{ij,ab}{\IS}(\omega) = \ERI{ij}{ab} +  2 \sum_m^{\Nocc \Nvir} \sERI{ij}{m}{\IS} \sERI{ab}{m}{\IS}
	\\
	\times  \qty(\frac{1}{\omega - \OmRPA{m}{\IS} + i \eta} - \frac{1}{\omega + \OmRPA{m}{\IS} - i \eta}),
\end{multline}
where the spectral weights at coupling strength $\IS$ read
\begin{equation}
	\sERI{pq}{m}{\IS} =  \sum_i^{\Nocc} \sum_a^{\Nvir} \ERI{pq}{ia} (\bX{\IS}_m + \bY{\IS}_m)_{ia}.
\end{equation}
In the case of complex orbitals, we refer the reader to Ref.~\onlinecite{Holzer_2019} for a correct use of complex conjugation in the spectral representation of $\W{}{}$.
Note that, in the case of {\GOWO}, the RPA neutral excitations in Eq.~\eqref{eq:W} are computed using the HF orbital energies.

In Eq.~\eqref{eq:W}, $\eta$ is a positive infinitesimal, and $\OmRPA{m}{\IS}$ are the direct (\ie, without exchange) RPA neutral excitation energies computed by solving the linear eigenvalue problem \eqref{eq:LR} with the following matrix elements
\begin{subequations}
\begin{align}
	\label{eq:LR_RPA-A}
	\ARPA{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eHF{a} - \eHF{i}) +   2 \IS \ERI{ia}{jb},
	\\ 
	\label{eq:LR_RPA-B}
	\BRPA{ia,jb}{\IS} & =   2 \IS \ERI{ia}{bj},
\end{align}
\end{subequations}
where $\eHF{p}$ are the Hartree-Fock (HF) orbital energies.

The relationship between the BSE formalism and the well-known RPAx (\ie, RPA with exchange) approach can be obtained by switching off the screening so that $\W{}{\IS}$ reduces to the bare Coulomb potential $\vc{}$. 
In this limit, the $GW$ quasiparticle energies reduce to the HF eigenvalues, and Eqs.~\eqref{eq:LR_BSE-A} and \eqref{eq:LR_BSE-B} to the RPAx equations: 
\begin{subequations}
\begin{align}
	\label{eq:LR_RPAx-A}
	\ARPAx{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eHF{a} - \eHF{i}) + \IS \qty[ 2 \ERI{ia}{jb} -    \ERI{ij}{ab} ],
	\\ 
	\label{eq:LR_RPAx-B}
	\BRPAx{ia,jb}{\IS} & =  \IS \qty[ 2 \ERI{ia}{bj} - \ERI{ib}{aj} ].
\end{align}
\end{subequations}

%%% FIG 1 %%%
%\begin{figure}
%	\includegraphics[width=\linewidth]{}
%\caption{
%\label{fig:}
%}
%\end{figure}
%%% %%% %%%


%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%
This is the conclusion

%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
%%%%%%%%%%%%%%%%%%%%%%%%
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}