\documentclass[10pt]{book}
\usepackage{graphicx}
\usepackage{amsmath}    %% pour mettre du texte en mode math
\usepackage{tcolorbox}

\newcommand\uom{\underline{\omega}}
\newcommand\Oms{{\Omega}_s}
\newcommand\hOms{\frac{{\Omega}_s}{2}}
\begin{document}

\chapter{Strinati for the big kids}

We wish to solve the self-consistent Bethe-Salpeter equation:
\begin{equation*}
    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{equation*}
where e.g. $1=(x_1,t_1)$ is a space/spin and time position and:
\begin{align*}
    iL(1,2;1',2') &= -G_2(1,2 ; 1',2') + G(1,1')G(2,2')   \\
    iL_0(1,2;1',2') &=  G(1,2')G(2,1') \\
    iG(1,2') &= \langle N | T   {\hat \psi}(1)  {\hat \psi}^{\dagger}(2')  | N \rangle \\
    i^2 G_2(1,2 ; 1',2') &= \langle N | T \Big[ {\hat \psi}(1)   {\hat \psi}(2) {\hat \psi}^{\dagger}(2') {\hat \psi}^{\dagger}(1') \Big] | N \rangle
\end{align*}
For optical properties, with instantaneous creation/destruction of excitations, we are interested in the collapse of the time variables starting with imposing $t_2'=t_2^{+}$ in the right-hand side $L(1,2;1',2')$ and left hand-side $L(6,2;5,2')$ response functions. Since $t_2'=t_2^{+}=(t_2 + 0^+)$  indicates that   the (2,2') operators cannot be separated, we see that the ordering operator selects two possible channels depending on the times ordering, namely:
\begin{align*}
i^2 G_2^{I}(1,2 ; 1',2') &= \theta(  t_{m}^{11'} -t_2 )\langle N | T [{\hat \psi}(1) {\hat \psi}^{\dagger}(1') ] T [ {\hat \psi}(2) {\hat \psi}^{\dagger}(2') ] | N \rangle \\
i^2 G_2^{II}(1,2 ; 1',2') &= \theta( t_2 - t_{M}^{11'} )\langle N | T [ {\hat \psi}(2) {\hat \psi}^{\dagger}(2') ]
T [{\hat \psi}(1) {\hat \psi}^{\dagger}(1') ] | N \rangle   
\end{align*}
with $t_{M}^{11'} = \max(t_1,t_{1'}) = | \tau_{11'} |/2 + t^{11'}$ and $t_{m}^{11'} = \min(t_1,t_{1'}) = - | \tau_{11'} |/2 + t^{11'}$ where $\tau_{11'} = t_1 - t_{1'}$ and $t^{11'} = (t_1 + t_{1'})/2$.
Introducing the complete set of eigenstates of the N-electron systems $\lbrace | N,s \rangle \rbrace$, where s index the eigenstates and $| N,s=0 \rangle = | N \rangle$ the ground-state, one obtains for $G^{I}$ :
\begin{align*}
i^2 G_2^{I}(1,2 ; 1',2') &= \theta( t_{m}^{11'} -t_2  )  \sum_{s>0} \langle N | T   {\hat \psi}(1) {\hat \psi}^{\dagger}(1')  | N,s \rangle 
 \langle N,s | T   {\hat \psi}(2) {\hat \psi}^{\dagger}(2')  | N \rangle  \\
  & + \theta(  t_{m}^{11'} -t_2  ) \frac{  G(1,1') }{ i } \frac{ G(2,2') }{ i }   \\
\end{align*}
where in the right-hand side of the first line the sum $(s>0)$ avoids the (s=0) ground-state contribution that builds the 1-body Green's functions product in the 2nd line. 
This yields a first channel for $L(1,2;1',2')$ namely:
$$
iL^{I}(1,2;1',2') = \theta( t_{m}^{11'} -t_2  )  \sum_{s>0} \langle N | T   {\hat \psi}(1) {\hat \psi}^{\dagger}(1')  | N,s \rangle
 \langle N,s | T   {\hat \psi}(2) {\hat \psi}^{\dagger}(2')  | N \rangle 
$$
Switching from the Heisenberg representation (time-dependent operators) to the Schr{\"o}dinger one for the field operators, with e.g.
$$
{\hat \psi}(2)  = {\hat \psi}(x_2 t_2) = e^{ i H t_2} {\hat \psi}(x_2) e^{ -i H t_2}
$$
one obtains
\begin{align*}
   \langle N,s | T  {\hat \psi}(2) {\hat \psi}^{\dagger}(2')   | N  \rangle  &=
 - \langle N,s |  e^{i H t_2^+} {\hat \psi}^{\dagger}(x_2') e^{i H (t_2 - t_2^+) } {\hat \psi}(x_2)   e^{-i H t_2} | N  \rangle  \\
  &= - \langle N,s |   {\hat \psi}^{\dagger}(x_{2'})   {\hat \psi}(x_2)    | N  \rangle \times e^{i  \Omega_s t_2 } \\
  & = {\tilde \chi}_s(x_2,x_{2'}) e^{i  \Omega_s t_2 }
\end{align*}
with $\;{\tilde \chi}_s(x_2,x_{2'}) = \langle N,s |  {\hat \psi}(x_2) {\hat \psi}^{\dagger}(x_{2'})       | N  \rangle$.
 The energy $\Omega_s = E_s^N - E_0^N$ is the s-th neutral excitation energy of the system, namely the energy we are looking for.
The (-) sign in the right-hand side comes from the commutation of the creation and destruction operators with $t_2^+ > t_2$. 

A similar calculation for the second channel [$t_2 > \max(t_1,t_{1'})$] leads to a  frequency dependence in $e^{ -i  \Omega_s t_2 }$.
Putting everything together, one obtains :
\begin{align*}
iL (1,2;1',2')  = & \;   \theta( t_{m}^{11'} -t_2  )  \sum_{s>0} \langle N | T   {\hat \psi}(1) {\hat \psi}^{\dagger}(1')  | N,s \rangle
     {\tilde \chi}_s(x_2,x_{2'}) \times e^{i  \Omega_s t_2 } \\
                      -  & \; \theta( t_2 - t_{M}^{11'}    )  \sum_{s>0} \langle N,s | T   {\hat \psi}(1) {\hat \psi}^{\dagger}(1')  | N  \rangle
 {\chi}_s(x_2,x_{2'})   \times e^{ -i  \Omega_s t_2 }
 \end{align*}
with $\;{\chi}_s(x_2,x_{2'}) = \langle N |  {\hat \psi}(x_2) {\hat \psi}^{\dagger}(x_{2'})  | N,s  \rangle$. \\


Performing the very same analysis for $L(6,2;5,2')$ and picking up a specific   $e^{ +i  \Omega_s t_2 }$  contribution, dividing by 
$ \langle N,s | T  {\hat \psi}(2) {\hat \psi}^{\dagger}(2')   | N  \rangle $   both sides of the Bethe-Salpeter equation leads to: \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
\begin{align}
\langle N & |  T   {\hat \psi}(1)   {\hat \psi}^{\dagger}(1')  | N,s \rangle  \theta( t_{m}^{11'} -t_2  )   = \\
&=  \int d3456 \; L_0(1,4;1',3) \Xi[3,5;4,6]   \langle N | T   {\hat \psi}(6) {\hat \psi}^{\dagger}(5)  | N,s \rangle   \theta( t_{m}^{56} -t_2  )   \nonumber
\end{align}
\end{tcolorbox}  
\noindent where the $L_0(1,2;1',2')$ does not contribute if \textbf{ we select some $\Omega_s$ frequency below the quasiparticle gap of the system}. This is the common situation in molecular systems where due to excitonic effects, namely electron-hole interaction, the lowest optical excitation energies are lower than the photoemission quasiparticle gap.  $ L_0(1,2;1',2')$   has indeed a lowest excitation energy at the energy gap $E_g$ and cannot participate  to the $e^{i \Omega_s t_2}$ response of the system if $\Omega_s < E_g$. \\
  This is equation (11.3) of Strinati with e.g. $$ \langle N  | T   {\hat \psi}(1) {\hat \psi}^{\dagger}(1')   | N,s \rangle = \chi_s(x_1,x_{1'}; \tau_1) e^{-i \Omega_s t^1}$$
(see equation G.3a of Appendix G) but keeping the Heaviside fonctions for a better definition of time integrals.     \\

\noindent \textbf{Spectral representation of $G(t_1-t_3)G(t_4-t_1^+)$ } \\

We now consider the $t_1$ time. Using the same algebra as above, one obtains
$$
\langle N | T  {\hat \psi}(1)    {\hat \psi}^{\dagger}(1')   | N,s \rangle  =  
  \langle N  | T   {\hat \psi}(x_1)  {\hat \psi}^{\dagger}(x_{1'})    | N,s \rangle  \times e^{-  i  \Oms t_1 }  
$$
imposing now $t_1' = t_1^+$. To extract the $e^{-i \Omega_s t_1 }$ behavior in the right hand-side of the BSE equation, one must find the Fourier components with respect to the $t_1$ variable of  $L_0(1,4;1',3) = -i G(1,3) G(4,1')$ by multiplying by $ \int dt_1 e^{i \omega_1 t_1}$   both sides of the BSE equation, \textcolor{red}{ with $( \omega_1 \rightarrow \Oms + i\delta )$ where we introduce the small infinitesimal $\delta = 0^+$ for defining properly the $t_1$   time integral on the  left-hand side: }  
$$
\int dt_1  e^{i \omega_1 t_1} e^{ -  i \Oms t_1 } \theta( t_1 -t_2  )  =    e^{-\delta t_2} \int d \tau  e^{- \delta \tau } \theta( \tau  ) = ( \delta^{-1} )
$$
This yields an annoying but useful factor that will canceled out in the following. 
Expressing again the   field operators in their Schr\"{o}dinger representation, one see that the  1-body Green's functions depend on the time-differences $\tau_{13} = t_1 - t_3$ and $\tau_{41} = t_4 - t_1$, respectively, with e.g.
$$
\langle N | \psi(1) \psi^{\dagger}(3) | N \rangle = \langle N | \psi(x_1) e^{-i H \tau_{13}} \psi^{\dagger}(x_3) | N \rangle e^{+i E_N^0 \tau_{13}} 
$$
Taking by convention the time difference as that of the destruction operator minus that of the creation analog and introducing our notations for the  Fourier representation  :
\begin{equation*}
  G( \tau_{13} ) = {1 \over 2\pi} \int d \omega \; G(\omega) e^{ - i \omega  \tau_{13}  }  \;\;\;\;\;
  G(\tau_{41} ) = {1 \over 2\pi} \int d \omega' \; G(\omega') e^{ - i {\omega}'  \tau_{41}  }  
\end{equation*}
one obtains (keeping only time/frequency variables for clarity)
$$
i L_0(1,4;1',3)  =  \int { d\omega \over 2\pi } { d{\omega}' \over 2\pi } \; G(\omega) G({\omega}')  e^{ -i \omega  \tau_{13}  } e^{ - i {\omega}'  \tau_{41} }   
$$
The Fourier components of ($iL_0$) with respect to the ($t_1$) time are: 
\begin{align*}
[iL_0]( \omega_1 )  &=   \int dt_1 e^{   i \omega_1 t_1 }   \int  { d\omega \over 2\pi } { d{\omega}' \over 2\pi }  \; G(\omega ) G( {\omega}' )     
    e^{- i \omega  (t_1 - t_3)  } e^{  -i {\omega}'  (t_4 - t_1) }   \\
    &=   {1 \over (2\pi)^2 }   \int d \omega d{\omega}'  \; G(\omega ) G( {\omega}' )      2\pi \delta(  -\omega + \omega_1 +  {\omega}' )
    e^{  i \omega    t_3   } e^{ - i {\omega}'   t_4   } \\
   &=  {1 \over  2\pi   } \int d \omega   \; G(\omega ) G( {\omega} - \omega_1 )   e^{   i \omega    t_3   } e^{  -i ( \omega - \omega_1)   t_4   } \\
   &= { 1 \over 2\pi } \int d \omega   \; G(\omega - {\omega_1 \over 2} ) G( {\omega} + {\omega_1 \over 2}  )   e^{  i  \omega  \tau_{34}  } e^{  i  \omega_1   t^{34}   }   
   \end{align*}
   where in the last line we made the change of variable $( \omega \rightarrow \omega +   {\omega_1} / 2 ) $ and introduced the average time
   $t^{34} = (t_3 + t_4) /2$.
   We now make the identification $( \omega_1 = {\Oms} + i\delta )$ as done to properly pick-up the left-hand side $e^{-i \Oms t_1}$ contribution and obtain   :
\begin{align*}
[iL_0 ](  \Omega_s)  
    = { 1 \over 2\pi } \int d \omega   \; G(\omega - \hOms ) G( {\omega} + \hOms )   e^{  i  \omega \tau_{34}  } e^{  i  \Oms^{+}   t^{34}   }   
\end{align*}   
with $( \Oms^{+} = \Oms + i\delta )$ keeping the infinitesimal $\delta$ where it will be needed.   
Using the   Lehman  representation of the 1-body Green's functions, 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) }
$$
where we reintroduced the (space, spin) variables, one obtains by picking the residues associated with the poles of each Green's function:
\begin{align*}
 \int d \omega  & \; G(x_1,x_3; \omega -  \hOms  )    G(x_4,x_{1'}; \omega + \hOms  )  e^{  i    \omega  \tau_{34}    } =   \\
  &  2 {i} \pi   \theta( \tau_{34} ) \sum_{nj}  \frac{ \phi_n(x_1) \phi_n^*(x_3) \phi_j(x_4) \phi_j^*(x_{1'}) } { \varepsilon_j + \Omega_s - \varepsilon_n + i \eta \times \text{sgn}(\varepsilon_n - \mu)  } e^{i (\varepsilon_j + {\Omega_s \over 2}  ) \tau_{34}  }  \\
  + &    2 {i} \pi   \theta( \tau_{34} ) \sum_{jn}  \frac{ \phi_j(x_1) \phi_j^*(x_3) \phi_n(x_4) \phi_n^*(x_{1'}) } { \varepsilon_j - \Omega_s - \varepsilon_n + i \eta \times \text{sgn}(\varepsilon_n - \mu)  } e^{i (\varepsilon_j - {\Omega_s \over 2}  ) \tau_{34}  }  \\
 - &   2 {i} \pi   \theta(- \tau_{34} ) \sum_{nb}  \frac{ \phi_n(x_1) \phi_n^*(x_3) \phi_b(x_4) \phi_b^*(x_{1'}) } { \varepsilon_b + \Omega_s - \varepsilon_n + i \eta \times \text{sgn}(\varepsilon_n - \mu)  } e^{i (\varepsilon_b + {\Omega_s \over 2}  ) \tau_{34}  }  \\
-  &   2 {i} \pi   \theta(- \tau_{34} ) \sum_{bn}  \frac{ \phi_b(x_1) \phi_b^*(x_3) \phi_n(x_4) \phi_n^*(x_{1'}) } { \varepsilon_b - \Omega_s - \varepsilon_n + i \eta \times \text{sgn}(\varepsilon_n - \mu)  } e^{i (\varepsilon_b - {\Omega_s \over 2}  ) \tau_{34}  }  
\end{align*}
where for $\tau_{34} > 0$ (resp. $\tau_{34} < 0$) the contour is close in the upper (resp. lower) half complex plane.  
We adopted chemist notations with (i,j) and (a,b) indexing occupied and virtual states, respectively. Projecting now on
$  \phi_a^*(x_1) \phi_i(x_{1'}) \; $  one obtains by orthonormalization :
\begin{align*}
\int dx_1 dx_{1'} \; \phi_a^*(x_1) \phi_i(x_{1'})   \int d \omega  & \; G(x_1,x_3; \omega + {\Omega_s \over 2} )    G(x_4,x_{1'}; \omega - {\Omega_s \over 2} )  e^{  i     \omega  \tau_{34}    } =   \\
     2  & {i} \pi   \theta( \tau_{34} )   \frac{   \phi_a^*(x_3) \phi_i(x_4)  } { \varepsilon_i + \Omega_s - \varepsilon_a + i \eta    } e^{i (\varepsilon_i + {\Omega_s \over 2}  ) \tau_{34}  }  \\
     -2 & {i} \pi   \theta( - \tau_{34} )   \frac{   \phi_a^*(x_3) \phi_i(x_4)  } { \varepsilon_a - \Omega_s - \varepsilon_i - i \eta    } e^{i (\varepsilon_a - {\Omega_s \over 2}  ) \tau_{34}  }
\end{align*}
After multiplication by $ (     \varepsilon_a - \varepsilon_i   -  \Omega_s -i \eta  )$ (and setting $\eta \rightarrow 0$), the Bethe-Salpeter equation becomes: \\

\begin{tcolorbox}[ width=\textwidth, colback={white}]
\begin{align*}
     ( \varepsilon_a - \varepsilon_i    & - \Omega_s    )   \int dx_1 \phi_a^*(x_1) \phi_i(x_{1'})      \langle N | T  {\hat \psi}(x_1)  {\hat \psi}^{\dagger}(x_{1'})     | N,s \rangle  (  \delta^{-1} )
    =   \nonumber \\
 &   - \int d34  \;  \phi_a^*(x_3) \phi_i(x_4) 
   \left[ \theta( \tau_{34} )    e^{i (\varepsilon_i + {\Omega_s \over 2}  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - {\Omega_s \over 2}  ) \tau_{34}  } \right] \\
&  \times 
 e^{  i     \Omega_s^{+}     t^{34} } \int d56 \;        \Xi[3,6;4,5]   \langle N  | T  {\hat \psi}(5) {\hat \psi}^{\dagger}(6)   | N,s \rangle     \times \theta( t^{56}_m - t_2 )
\end{align*}
\end{tcolorbox}
\vskip .3cm \noindent This is equation  (11.5) by Strinati with the same
$$
  \theta( \tau  )    e^{i (\varepsilon_i + {\Omega_s \over 2}  ) \tau   }  + \theta( - \tau  ) e^{ i (\varepsilon_a - {\Omega_s \over 2} ) \tau  } = 
 e^{i   {\Omega_s \over 2}  \tau   } \Big(  \theta( \tau )    e^{i \varepsilon_i   \tau   }  + \theta( - \tau  ) e^{i (\varepsilon_a -  \Omega_s )\tau  }   \Big)
 $$ 
 time-structure where $\tau = \tau_{34} = t_3 - t_4 $. \\

\noindent{\textbf{The $GW$ approximation. }   Adopting the $GW$ approximation, and neglecting the $(\partial W / \partial G)$ of higher order in $W$, one obtains:
\begin{equation}
    \Xi[3,6;4,5] = v(3,5) \delta(3,4) \delta(5,6) - W(3^+,4) \delta(3,5) \delta(4,6)
\end{equation}
so that the Bethe-Salpeter equation becomes:
\begin{align*}
     ( \varepsilon_a - \varepsilon_i    & - \Omega_s    )   \int dx_1 \phi_a^*(x_1) \phi_i(x_{1'})   \langle N | T {\hat \psi}(x_1)   {\hat \psi}^{\dagger}(x_{1'})    | N,s \rangle  (  \delta^{-1} )
       =   \nonumber \\
 &   - \int d35 \;  \phi_a^*(x_3) \phi_i(x_3) e^{  i     \Omega_s^{+}     t_3 }   v(3,5)    \langle N  | T [{\hat \psi}(5) {\hat \psi}^{\dagger}(5) ] | N,s \rangle \times \theta( t_5 - t_2 ) \\
 & + \int d34 \;  \phi_a^*(x_3) \phi_i(x_4)
   \left[  \theta( \tau_{34} )    e^{i (\varepsilon_i + {\Omega_s \over 2}  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - {\Omega_s \over 2}  ) \tau_{34}  } \right] \times     \\
& \hskip 1cm  \times 
 e^{  i     \Omega_s^{+}   t^{34} } W(3^+,4)   \langle N  | T [{\hat \psi}(3) {\hat \psi}^{\dagger}(4) ] | N,s \rangle   \times \theta( t^{34}_m - t_2 )
\end{align*}
\\

\noindent\textbf{Spectral representation of  $\langle N  | T {\hat \psi}(3) {\hat \psi}^{\dagger}(4) | N,s \rangle$ } \\

\noindent Starting with:
\begin{align*}
\langle N  | T  {\hat \psi}(3) {\hat \psi}^{\dagger}(4)   | N,s \rangle &= \theta( \tau_{34} )  \langle N  |  {\hat \psi}(3) {\hat \psi}^{\dagger}(4)   | N,s \rangle \\
&- \theta( -\tau_{34} )  \langle N  |   {\hat \psi}^{\dagger}(4) {\hat \psi}(3)  | N,s \rangle
\end{align*}
and decomposing the static (Schr\"{o}dinger) field operators over the Hartree-Fock   molecular orbitals creation/destruction operators:
$$
 {\hat \psi}^{\dagger}(x_4) =    \sum_n \phi_n^*(x_4) {\hat a}_n^{\dagger}    \;\;\; \text{ and} \;\;\; {\hat \psi}(x_3) = \sum_m \phi_m(x_3) {\hat a}_m $$
one obtains
\begin{align*}
\langle N  | T {\hat \psi}(3) {\hat \psi}^{\dagger}(4)  | N,s \rangle &= \sum_{mn}  \phi_m(x_3) \phi_n^*(x_4)  \times \\
\times & \Big( \theta( \tau_{34} ) \langle N  |  {\hat a}_m e^{- i H \tau_{34} }  {\hat a}^{\dagger}_n   | N,s \rangle e^{-i E_s^N t_4} e^{+ i E_0^N t_3}  \\
&- \theta( -\tau_{34} )  \langle N  |   {\hat a_n}^{\dagger} e^{  i H \tau_{34} } {\hat a}_m    | N,s \rangle  e^{ - i E_s^N t_3} e^{  + i E_0^N t_4}  \Big)
\end{align*}
Taking the conjugate of the $  {\hat a}_m e^{- i H \tau_{34} } $ and
$ {\hat a_n}^{\dagger} e^{  i H \tau_{34} } $ operators to act on the ground-state $| N \rangle$ state, one obtains:
$$
e^{  i H \tau  } {\hat a}_m^{\dagger}   | N  \rangle = e^{ i (E_0^N + \varepsilon_m ) \tau  }  {\hat a} _m^{\dagger}   | N  \rangle 
\;\;\;\;\; \text{and} \;\;\;\;
e^{ - i H \tau  } {\hat a}_n   | N \rangle = e^{ - i ( E_0^N - \varepsilon_n )  \tau  } {\hat a}_n   | N  \rangle
  $$
This properties hold in the case of Hartree-Fock formalism (Koopman's theorem) but holds more generally if $\lbrace \varepsilon_m , \varepsilon_n \rbrace$ are quasiparticle energies, 
namely true electronic removal or addition energies such as in the $GW$ formalism.
Forming the corresponding bra, one obtains :  
\begin{align*}
\langle N  | T   {\hat \psi}(3)  & {\hat \psi}^{\dagger}(4)  | N,s \rangle  = - \sum_{mn}  \phi_m(x_3) \phi_n^*(x_4)   \langle N | {\hat a}_n^{\dagger} {\hat a}_m  | N,s \rangle \times \\
\times & \Big( \theta( \tau_{34} ) e^{- i   \varepsilon_m   \tau_{34}  }  e^{ -i    \Oms   t_4    }    
 + \theta( -\tau_{34} )    e^{ - i   \varepsilon_n    \tau_{34}  }  e^{ - i    \Oms  t_3    }    \Big)
\end{align*}
 With $t_3 =  \tau_{34}/2 + t^{34}$ and $t_4 = -\tau_{34}/2 + t^{34}$, one finally obtains :
\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_{34} )   e^{- i   ( \varepsilon_m  -  \hOms ) \tau_{34}  }  
 + \theta( -\tau_{34} )   e^{ - i   ( \varepsilon_n  + \hOms)  \tau_{34}  }      \Big)
\end{align*}
We must now discuss the (m,n) indexes. The most natural channel for ${\hat a}_m^{\dagger} {\hat a}_n | N \rangle$ to project on an excited state occurs
for (n,m)   indexing (occupied,virtual) orbitals, respectively, in order to promote an electron from the occupied to the empty orbitals. 
We will label $A_{jb}^{s} = \langle N | {\hat a}_j^{\dagger} {\hat a}_b  | N,s \rangle$ these large amplitude projections. However, in the case of a generic correlated system,
the ground-state does not project exclusively on the  mono-determinant built from the occupied 1-body orbitals $\lbrace \phi_i \rbrace$, but contains also contributions from
singly, doubly, etc. excited determinants in the full $\lbrace \phi_n \rbrace$ space. As such, the $\langle N | {\hat a}_b^{\dagger} {\hat a}_j  | N,s \rangle$ weights are small in general but not strictly zero.
We will label $B_{jb}^{s}$ these small amplitudes. Neglecting these $B_{jb}^{s}$ coefficients leads to the Tamm-Dancoff approximation (TDA). We start with the TDA in the following for simplicity, adding the "small B" contribution in a second step.  \\

\noindent In the TDA, namely with $B_{jb}^{s} = 0$, one obtains :
\begin{align*}
\langle N  | T   {\hat \psi}(3)   & {\hat \psi}^{\dagger}(4)    | N,s \rangle   \cong - \Big( e^{  -i \Omega_s t^{34} }  \Big)    \sum_{jb} \phi_b(x_3) \phi_j^*(x_4)  A_{jb}^{s}    \times \nonumber \\
\times & \Big( \theta( \tau_{34} )   e^{- i   ( \varepsilon_b  -  \hOms ) \tau_{34}  }  
 + \theta( -\tau_{34} )   e^{ - i   ( \varepsilon_j  + \hOms)  \tau_{34}  }      \Big)
\end{align*}
As a result, the integral  associated with the screened Coulomb potential reads:
\begin{align*}
  \int d34 \;  \phi_a^*(x_3) & \phi_i(x_4)  
   \left[  \theta( \tau_{34} )    e^{i (\varepsilon_i + {\Omega_s \over 2}  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - {\Omega_s \over 2}  ) \tau_{34}  } \right] \times     \\
&    \times 
 e^{  i     \Omega_s^{+}   t^{34} } W(3^+,4)   \langle N  | T [{\hat \psi}(3) {\hat \psi}^{\dagger}(4) ] | N,s \rangle   \times \theta( t^{34}_m - t_2 )   \\
 = -   \int d34 \;  \phi_a^*(x_3) & \phi_i(x_4)  \theta( t^{34}_m   - t_2 )  W(3^+,4) \sum_{jb} \phi_b(x_3) \phi_j^*(x_4)  A_{jb}^{s} e^{ - \delta   t^{34} }  \times \\
  &  \times \left[  \theta( \tau_{34} )    e^{i ( \varepsilon_i - \varepsilon_b + \Oms  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - \varepsilon_j - \Oms  ) \tau_{34}  } \right]      \\
 = -   \sum_{jb} A_{jb}^{s}   &    \int d\tau_{34} dt^{34} e^{ - \delta   t^{34} } \theta( t^{34}  - | \tau_{34} |/2 - t_2 ) W_{ij,ab}( \tau_{34}^{+} )   \times \\
 &  \times \left[  \theta( \tau_{34} )    e^{i ( \varepsilon_i - \varepsilon_b + \Oms  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - \varepsilon_j - \Oms  ) \tau_{34}  } \right] 
\end{align*}
where we have explicited $t^{34}_m$ and made a change of variable   $(t_3,t_4) \rightarrow (\tau_{34}, t^{34})$ with unit wronskian, defining the screened Coulomb matrix elements
\begin{align*}
W_{ij,ab}( \tau_{34}^{+} )&= \int d{\bf r}_3 d{\bf r}_4 \;  \phi_i(x_4)  \phi_j^*(x_4)   W(3^+,4) \phi_a^*(x_3)  \phi_b(x_3)  \;\;\;\; 
\end{align*}
with $(\tau_{34}^{+} = t_3^{+} - t_4)$ and adopting Mulliken notations where we group together the index associated with orbitals taken at the same space position (and putting complexe conjugate orbitals as inner indexes).
The $t^{34}$ time integral  yields again $( \delta^{-1} )$ as done previously for the  $t_1$-time integral, leaving only the dependence on the time difference $\tau_{34} = t_3 - t_4$.
We can understand the origin of this problematic $( \delta^{-1} )$ factor : the screened- Coulomb potential contribution only depends on the time difference $\tau_{34}$ so that the time variable $t^{34}$ becomes problematic.  

To deal now with the bare Coulomb contribution, we need 
 $\langle N  | T   {\hat \psi}(5)   {\hat \psi}^{\dagger}(5)  | N,s \rangle$ that in the TDA reads  :
$$
\langle N  | T   {\hat \psi}(5)   {\hat \psi}^{\dagger}(5)  | N,s \rangle  \cong - \Big( e^{  -i \Omega_s t_5 }  \Big)    \sum_{jb} \phi_b(x_5) \phi_j^*(x_5)  A_{jb}^{s} 
$$
leading to the following contribution :
\begin{align*}
    - \int d35 \;  \phi_a^*(x_3) & \phi_i(x_3) e^{  i     \Omega_s^{+}     t_3 }   v(3,5)    \langle N  | T [{\hat \psi}(5) {\hat \psi}^{\dagger}(5) ] | N,s \rangle \times \theta( t_5 - t_2 ) \\
 & =  \sum_{jb}   A_{jb}^{s} v_{ia,jb} \int dt_3 dt_5 \delta( t_3 - t_5) e^{  i     \Omega_s^{+}     t_3 }   e^{  -i \Omega_s t_5 }\theta( t_5 - t_2 )
\end{align*}
with $v(3,5) = v( {\bf r}_3 - {\bf r}_5 ) \delta( t_3 - t_5)$ a static interaction and :
$$
v_{ia,jb}    = \int d{\bf r}_3 d{\bf r}_5 \;  \phi_i(x_3) \phi_a^*(x_3) v( {\bf r}_3 - {\bf r}_5 ) \phi_j^*(x_5) \phi_b(x_5) 
$$
The time integral reads now $\int dt_3 e^{- \delta t_3} \theta(t_3 - t_2 )$ yielding $( \delta^{-1} )$ again. 
This factor in front of all terms can now be factorized out.  \\

\noindent Finally, for the right-hand side of the Bethe-Salpeter equation, we develop: 
\begin{equation*}
\langle N  |   {\hat \psi}^{\dagger}(x_{1'})   {\hat \psi}(x_1)   | N,s \rangle   =    \sum_{jb} \phi_b(x_1) \phi_j^*(x_1)  A_{jb}^{s}
\end{equation*}
so that by orthonormalization:
\begin{equation*}
   \int dx_1 dx_{1'} \; \phi_a^*(x_1) \phi_i(x_{1'})    \langle N |   {\hat \psi}^{\dagger}(x_{1'})   {\hat \psi}(x_1)     | N,s \rangle   =  
      \sum_{jb}   A_{jb}^{s} \delta_{ij} \delta_{ab}   =    A_{ia}^{s}
\end{equation*}
Note that the projection on $\phi_a^*(x_1) \phi_i(x_{1'})$ will always cancel the small $B_{jb}^{s}$ contribution so that this last expression is always valid.
We obtain  thus a much simplified reformulation of the dynamical TDA Bethe-Salpeter equation: \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
\begin{align*}
      ( \varepsilon_a - \varepsilon_i    - \Omega_s    )   &   A_{ia}^{s}         
       + \sum_{jb} A_{jb}^{s} \times v_{ai,bj}   \stackrel{\text{TDA}}{=}         \sum_{jb} A_{jb}^{s}     \times \int d\tau   \;   W_{ij,ab}( {\tau}^{+} )     \times  \\         
    &   \times \left[  \theta( \tau )    e^{i (\varepsilon_i -  \varepsilon_b+  \Oms  ) \tau  }  + \theta( - \tau ) e^{i (\varepsilon_a - \varepsilon_j  - \Oms  ) \tau  } \right]  
\end{align*}
\end{tcolorbox} 
\noindent where we dropped the index of ${\tau}_{34}$. Taking finally the Fourier representation of $W$, the time integration reads:
\begin{align*}
  \widetilde{W}_{ij,ab}(\Oms) =&     \int    d\tau   \;    W_{ij,ab}( \tau^{+} )      
   \times \;  \left[  \theta( \tau  )    e^{i (\varepsilon_i -  \varepsilon_b+  \Oms  ) \tau   }  + \theta( - \tau  ) e^{i (\varepsilon_a - \varepsilon_j  - \Oms  ) \tau   } \right]     \\
  =& \int { d\omega \over 2\pi } W_{ij,ab}( \omega  ) \int d\tau \;  e^{-i \omega ( \tau + 0^+) }  
     \left[  \theta( \tau  )    e^{i (\varepsilon_i -  \varepsilon_b+  \Oms  ) \tau   }  + \theta( - \tau  ) e^{i (\varepsilon_a - \varepsilon_j  - \Oms  ) \tau   } \right]  \\
   =&  {  i \over 2\pi} \int   d\omega    e^{-i \omega 0^{+}} \left[   \frac{  W_{ij,ab}( \omega  ) }{   ( \Oms - \omega ) - ( \varepsilon_b -  \varepsilon_i) + i\eta   }   
        +  \frac{   W_{ij,ab}( \omega  ) }{   ( \Oms + \omega) - ( \varepsilon_a -  \varepsilon_j ) + i\eta   }  \right]
\end{align*}
where we have introduced again a small damping $( \Oms \rightarrow \Oms + i\eta)$ of the excitation for convergency. This leads to
equation (11.8) by Strinati: \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
\begin{align*}
      ( \varepsilon_a - \varepsilon_i    - \Omega_s    )   &   A_{ia}^{s}          
       + \sum_{jb} A_{jb}^{s} \times v_{ai,bj}  \stackrel{\text{TDA}}{=}    { i \over 2 \pi}   \sum_{jb} A_{jb}^{s}    \int d\omega e^{-i \omega 0^+ } W_{ij,ab} \times \\
    \hskip 1cm &\times \left[ \frac{1}{( \Oms - \omega) - (\varepsilon_b - \varepsilon_i) +i \eta } + \frac{1}{ (\Oms + \omega) - (\varepsilon_a - \varepsilon_j) + i\eta } \right]   
\end{align*}
\end{tcolorbox}
\vskip .3cm  \noindent or to a form similar to the static  limit : \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
 \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} \stackrel{\text{TDA}}{=} 0
\end{align}
\end{tcolorbox}
\vskip .3cm  \noindent but with an effective dynamical  screened Coulomb potential (see Pina  eq. 24): \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
 \begin{align}
 \widetilde{W}_{ij,ab}(\Oms)  &= { i \over 2  \pi}  \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab} \times \\
    \hskip 1cm &\times \left[ \frac{1}{ (\Oms - \omega) - (\varepsilon_b - \varepsilon_i) +i \eta } + \frac{1}{ (\Oms + \omega) - (\varepsilon_a - \varepsilon_j) + i\eta } \right] \nonumber
\end{align}
\end{tcolorbox}
\vskip .3cm  \noindent Taking now the  spectral representation for the (time-ordered) $W_{ij,ab}( \omega )$, namely:
\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*}
the frequency integral becomes, closing in the lower half plane thanks to the $e^{-i\omega 0^+}$ factor.
 \textcolor{red}{(This is the only place where the (+) in $3^+$ makes a difference as for the $GW$ operator ... check carefully origin ...)}:
\begin{align*}
  \int  d\omega    e^{-i \omega 0^{+}} & \Big(   \frac{  W_{ij,ab}( \omega  ) }{   ( \Oms - \omega ) - ( \varepsilon_b -  \varepsilon_i) + i\eta   }   
        +  \frac{   W_{ij,ab}( \omega  ) }{   ( \Oms + \omega) - ( \varepsilon_a -  \varepsilon_j ) + i\eta   }  \Big)= \\
        = &  (-2i \pi) W_{ij,ab}( - \Oms + (\varepsilon_a - \varepsilon_j)  ) +  (-2i \pi) \times 2 \sum_m^{OV} [ij|m] [ab|m] \times \\
        & \times \Big(   \frac{  1 }{   ( \Oms - \Omega_m^{RPA}  ) - ( \varepsilon_b -  \varepsilon_i) + i\eta   }   
        +  \frac{  1 }{   ( \Oms + \Omega_m^{RPA} ) - ( \varepsilon_a -  \varepsilon_j ) + i\eta   }  \Big)
\end{align*}
that is :
\begin{align*}
 \widetilde{W}_{ij,ab}( \Oms ) = &  (ij|ab) + 2 \sum_m^{OV}    [ij|m] [ab|m] \times \\ 
   \times  &   \Big(  \frac{ 1}{ -\Oms + (\varepsilon_a - \varepsilon_j) -\Omega_m^{RPA} + i\eta }  -   \frac{ 1}{ -\Oms + (\varepsilon_a - \varepsilon_j) + \Omega_m^{RPA} - i\eta }   \\
 +     &     \frac{  1 }{   ( \Oms - \Omega_m^{RPA}  ) - ( \varepsilon_b -  \varepsilon_i) + i\eta   }   
        +  \frac{  1  }{   ( \Oms + \Omega_m^{RPA} ) - ( \varepsilon_a -  \varepsilon_j ) + i\eta   }  \Big)
\end{align*} \\
or
\begin{align*}
 \widetilde{W}_{ij,ab}( \Oms ) = &  (ij|ab) + 2 \sum_m^{OV}    [ij|m] [ab|m] \times \\ 
   \times  &   \Big( \frac{  1 }{   ( \Oms - \Omega_m^{RPA}  ) - ( \varepsilon_b -  \varepsilon_i) + i\eta   }   +   \frac{ 1}{ ( \Oms   - \Omega_m^{RPA} ) - (\varepsilon_a - \varepsilon_j)  + i\eta }   \\
 +     &         \frac{  1  }{   ( \Oms + \Omega_m^{RPA} ) - ( \varepsilon_a -  \varepsilon_j ) + i\eta   } - \frac{ 1}{ ( \Oms +\Omega_m^{RPA} ) - (\varepsilon_a - \varepsilon_j)  - i\eta } \Big)
\end{align*} \\
Keeping only \textbf{ the two resonant contributions : } (ACTUALLY the other ones in the last line cancel  when $\eta \rightarrow 0$ !! check ...) one obtains Pina's equation (27) \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
\begin{align}
 \widetilde{W}_{ij,ab}&( \Oms ) =   (ij|ab) + 2 \sum_m^{OV}    [ij|m] [ab|m] \times \\ 
   \times  &   \left[ \frac{  1 }{   ( \Oms - \Omega_m^{RPA}  ) - ( \varepsilon_b -  \varepsilon_i) + i\eta   }   +   \frac{ 1}{ ( \Oms   - \Omega_m^{RPA} ) - (\varepsilon_a - \varepsilon_j)  + i\eta }
     \right]   \nonumber
\end{align} 
\end{tcolorbox} 
\noindent This dynamical $\widetilde{W}( \Oms )$ kernel can be rewritten:
\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) $. This formulation is similar to the spectral representation of $W_{ij,ab}(\omega)$ but with 
$( \omega \rightarrow \Omega_{ib}^{s} )$ and $( \Omega \rightarrow \Omega_{ja}^{s} )$.

\vskip 1cm
\noindent {\textbf {Non-resonant contributions.}} We now add the $B_{jb}^{s}$ small contributions to the $\langle N  | T   {\hat \psi}(3)   {\hat \psi}^{\dagger}(4)  | N,s \rangle $ weight:
\begin{align*}
\langle N  | T   {\hat \psi}(3)   {\hat \psi}^{\dagger}(4)  | N,s \rangle^{\text{corr}} &= - \Big( e^{  -i \Omega_s t^{34} }  \Big)    \sum_{bj} \phi_j(x_3) \phi_b^*(x_4) B_{jb}^{s}     \times \nonumber \\
\times & \Big( \theta( \tau_{34} )   e^{- i   ( \varepsilon_j  -  \hOms ) \tau_{34}  }  
 + \theta( -\tau_{34} )   e^{ - i   ( \varepsilon_b  + \hOms)  \tau_{34}  }      \Big)
\end{align*}
where ``corr" stands for correction. As a result, the integral  associated with the screened Coulomb potential must be corrected by:
\begin{align*}
  \int d34 \;  \phi_a^*(x_3) & \phi_i(x_4)  
   \left[  \theta( \tau_{34} )    e^{i (\varepsilon_i + {\Omega_s \over 2}  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - {\Omega_s \over 2}  ) \tau_{34}  } \right] \times     \\
&    \times 
 e^{  i     \Omega_s^{+}   t^{34} } W(3^+,4)   \langle N  | T [{\hat \psi}(3) {\hat \psi}^{\dagger}(4) ] | N,s \rangle   \times \theta( t^{34}_m - t_2 )   \\
 = -   \int d34 \;  \phi_a^*(x_3) & \phi_i(x_4)  \theta( t^{34}_m   - t_2 )  W(3^+,4) \sum_{bj} \phi_j(x_3) \phi_b^*(x_4)  B_{jb}^{s} e^{ - \delta   t^{34} }  \times \\
  &  \times \left[  \theta( \tau_{34} )    e^{i ( \varepsilon_i - \varepsilon_j + \Oms  ) \tau_{34}  }  + \theta( - \tau_{34} ) e^{i (\varepsilon_a - \varepsilon_b - \Oms  ) \tau_{34}  } \right]      \\
 = - ({ \delta}^{-1} )   \sum_{jb} B_{jb}^{s}  &   \int d\tau  \;  W_{ib,aj}( {\tau}^{+} )  \times \\
   & \times \left[  \theta( \tau  )    e^{i ( \varepsilon_i - \varepsilon_j + \Oms  ) \tau   }  + \theta( - \tau  ) e^{i (\varepsilon_a - \varepsilon_b - \Oms  ) \tau   } \right]
 \end{align*}
The time integration yields now:
\begin{align*}
  \widetilde{W}_{ib,aj}(\Oms) =&     \int    d\tau   \;    W_{ib,aj}( \tau^{+} )      
   \times \;  \left[  \theta( \tau  )    e^{i (\varepsilon_i -  \varepsilon_j+  \Oms  ) \tau   }  + \theta( - \tau  ) e^{i (\varepsilon_a - \varepsilon_b  - \Oms  ) \tau   } \right]     \\
  =& \int { d\omega \over 2\pi } W_{ib,aj}( \omega  ) \int d\tau \;  e^{-i \omega ( \tau + 0^+) }  
     \left[  \theta( \tau  )    e^{i (\varepsilon_i -  \varepsilon_j+  \Oms  ) \tau   }  + \theta( - \tau  ) e^{i (\varepsilon_a - \varepsilon_b  - \Oms  ) \tau   } \right]  \\
   =&  {  i \over 2\pi} \int   d\omega    e^{-i \omega 0^{+}} \left[   \frac{  W_{ib,aj}( \omega  ) }{   ( \Oms - \omega ) - ( \varepsilon_j -  \varepsilon_i) + i\eta   }   
        +  \frac{   W_{ib,aj}( \omega  ) }{   ( \Oms + \omega) - ( \varepsilon_a -  \varepsilon_b ) + i\eta   }  \right]
\end{align*} \\
Plugging now the spectral representation of $W_{ib,aj}( \omega  )$:
\begin{align*}
 \widetilde{W}_{ib,aj}( \Oms ) = &  (ib|aj) + 2 \sum_m^{OV}    [ib|m] [aj|m] \times \\ 
   \times  &   \Big( \frac{  1 }{   ( \Oms - \Omega_m^{RPA}  ) - ( \varepsilon_j -  \varepsilon_i) + i\eta   }   +   \frac{ 1}{ ( \Oms   - \Omega_m^{RPA} ) - (\varepsilon_a - \varepsilon_b)  + i\eta }   \\
 +     &         \frac{  1  }{   ( \Oms + \Omega_m^{RPA} ) - ( \varepsilon_a -  \varepsilon_b ) + i\eta   } - \frac{ 1}{ ( \Oms +\Omega_m^{RPA} ) - (\varepsilon_a - \varepsilon_b)  - i\eta } \Big)
\end{align*} 
where again the last line contribution cancel. \\

   
\noindent The correction to the bare Coulomb contribution stems from the correction to  
 $\langle N  | T   {\hat \psi}(5)   {\hat \psi}^{\dagger}(5)  | N,s \rangle$, namely   :
$$
\langle N  | T   {\hat \psi}(5)   {\hat \psi}^{\dagger}(5)  | N,s \rangle^{\text{corr}}  = - \Big( e^{  -i \Omega_s t_5 }  \Big)    \sum_{bj} \phi_j(x_5) \phi_b^*(x_5)  B_{jb}^{s} 
$$
leading to the following correction :
\begin{align*}
    - \int d35 \;  \phi_a^*(x_3) & \phi_i(x_3) e^{  i     \Omega_s^{+}     t_3 }   v(3,5)    \langle N  | T [{\hat \psi}(5) {\hat \psi}^{\dagger}(5) ] | N,s \rangle \times \theta( t_5 - t_2 ) \\
 & =  \sum_{bj}   B_{bj}^{s} v_{ia,bj} ( \delta^{-1} )
 \end{align*}

As discussed above, there is no correction to the diagonal part of the BSE Hamiltonian so that one obtains beyond TDA so that: \\
\begin{tcolorbox}[ width=\textwidth, colback={white}]
 \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} \nonumber \\
        + & \sum_{bj}\Big(  v_{ai,jb}  - \widetilde{W}_{ib,aj}(\Oms)    \Big)  B_{jb}^{s} = 0
\end{align}
\end{tcolorbox}
\noindent with:
\begin{align*}
 \widetilde{W}_{ib,aj}( \Oms ) = &  (ib|aj) + 2 \sum_m^{OV}    [ib|m] [aj|m] \times \\ 
   \times  &   \left[ \frac{  1 }{   ( \Oms - \Omega_m^{RPA}  ) - ( \varepsilon_j -  \varepsilon_i) + i\eta   }   +   \frac{ 1}{ ( \Oms   - \Omega_m^{RPA} ) - (\varepsilon_a - \varepsilon_b)  + i\eta }    \right]
\end{align*}


%%%%%%%%%%
\end{document}