PTEROSOR/BSEdyn
10
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

saving work

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-05-26 22:47:22 +02:00
parent fb06abc6d6
commit 5365979f32
1 changed files with 84 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

202
BSEdyn.tex
View File

@ -106,6 +106,7 @@
\newcommand{\XiBSE}[1]{\Xi_{#1}} \newcommand{\XiBSE}[1]{\Xi_{#1}}
\newcommand{\Po}[1]{P_{#1}} \newcommand{\Po}[1]{P_{#1}}
\newcommand{\W}[2]{W_{#1}^{#2}} \newcommand{\W}[2]{W_{#1}^{#2}}
\newcommand{\tW}[2]{\widetilde{W}_{#1}^{#2}}
\newcommand{\Wc}[1]{W^\text{c}_{#1}} \newcommand{\Wc}[1]{W^\text{c}_{#1}}
\newcommand{\vc}[1]{v_{#1}} \newcommand{\vc}[1]{v_{#1}}
\newcommand{\Sig}[1]{\Sigma_{#1}} \newcommand{\Sig}[1]{\Sigma_{#1}}
@ -368,10 +369,10 @@ Dropping the (space/spin) variables, the Fourier components with respect to $t_1
with $\tau_{34} = t_3 - t_4$ and $t^{34} = (t_3 + t_4)/2$. with $\tau_{34} = t_3 - t_4$ and $t^{34} = (t_3 + t_4)/2$.
We now adopt the Lehman representation of the one-body Green's function in the quasiparticle approximation, \ie, We now adopt the Lehman representation of the one-body Green's function in the quasiparticle approximation, \ie,
\begin{equation} \label{eq:G-Lehman} \begin{equation} \label{eq:G-Lehman}
G(\bx_1,\bx_2 ; \omega) = \sum_p \frac{ \phi_p(\bx_1) \phi_p^*(\bx_2) } { \omega - \e{p} + i \eta \times \text{sgn} (\e{p} - \mu) }, G(\bx_1,\bx_2 ; \omega) = \sum_p \frac{ \MO{p}(\bx_1) \MO{p}^*(\bx_2) } { \omega - \e{p} + i \eta \times \text{sgn} (\e{p} - \mu) },
\end{equation} \end{equation}
where $\mu$ is the chemical potential. where $\mu$ is the chemical potential.
The $\e{p}$'s in Eq.~\eqref{eq:G-Lehman} are quasiparticle energies (\ie, proper addition/removal energies) and the $\phi_p$'s are their associated one-body (spin)orbitals. The $\e{p}$'s in Eq.~\eqref{eq:G-Lehman} are quasiparticle energies (\ie, proper addition/removal energies) and the $\MO{p}$'s are their associated one-body (spin)orbitals.
%where the $\eps_{p}$'s are proper addition/removal energies such that %where the $\eps_{p}$'s are proper addition/removal energies such that
%\begin{equation} %\begin{equation}
% e^{i \hH \tau} \ha_p^{\dagger} \ket{N} = e^{ i (E_0^N + \e{p} ) \tau } \ha_p^{\dagger} \ket{N}, % e^{i \hH \tau} \ha_p^{\dagger} \ket{N} = e^{ i (E_0^N + \e{p} ) \tau } \ha_p^{\dagger} \ket{N},
@ -379,12 +380,12 @@ The $\e{p}$'s in Eq.~\eqref{eq:G-Lehman} are quasiparticle energies (\ie, proper
%$\hH$ being the exact many-body Hamiltonian. %$\hH$ being the exact many-body Hamiltonian.
In the following, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals. In the following, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals.
%\titou{namely $GW$ quasiparticle energies and input Hartree-Fock molecular orbitals in the present study. (T2: shall we really mention this here?)} %\titou{namely $GW$ quasiparticle energies and input Hartree-Fock molecular orbitals in the present study. (T2: shall we really mention this here?)}
Projecting the Fourier component $L_0(\bx_1,4;\bx_{1'},3; \omega_1 = \Oms )$ onto $\phi_a^*(\bx_1) \phi_i(\bx_{1'})$ yields Projecting the Fourier component $L_0(\bx_1,4;\bx_{1'},3; \omega_1 = \Oms )$ onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$ yields
\begin{multline} \begin{multline}
\iint d\bx_1 d\bx_{1'} \; \phi_a^*(\bx_1) \phi_i(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Oms) \iint d\bx_1 d\bx_{1'} \, \MO{a}^*(\bx_1) \MO{i}(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Oms)
\\ \\
= =
\frac{ \phi_a^*(\bx_3) \phi_i(\bx_4) e^{i \Oms t^{34} }} { \Oms - ( \e{a} - \e{i} ) + i \eta } \frac{ \MO{a}^*(\bx_3) \MO{i}(\bx_4) e^{i \Oms t^{34} }} { \Oms - ( \e{a} - \e{i} ) + i \eta }
\qty[ \theta( \tau_{34} ) e^{i ( \e{i} + \hOms) \tau_{34} } + \theta( - \tau_{34} ) e^{i (\e{a} - \hOms \tau_{34}) } ]. \qty[ \theta( \tau_{34} ) e^{i ( \e{i} + \hOms) \tau_{34} } + \theta( - \tau_{34} ) e^{i (\e{a} - \hOms \tau_{34}) } ].
\end{multline} \end{multline}
% and $(i,j)$/$(a,b)$ index occupied/virtual orbitals, respectively. % and $(i,j)$/$(a,b)$ index occupied/virtual orbitals, respectively.
@ -395,7 +396,7 @@ For example, we have
\begin{multline} \begin{multline}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s} \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
\\ \\
= - \qty( e^{ -i \Omega_s \tau^{56} } ) \sum_{pq} \phi_p(\bx_6) \phi_q^*(\bx_5) = - \qty( e^{ -i \Omega_s \tau^{56} } ) \sum_{pq} \MO{p}(\bx_6) \MO{q}^*(\bx_5)
\mel{N}{\ha_q^{\dagger} \ha_p}{N,s} \mel{N}{\ha_q^{\dagger} \ha_p}{N,s}
\\ \\
\times \qty[ \theta( \tau_{56} ) e^{- i ( \e{p} - \hOms ) \tau_{56} } + \theta( - \tau_{56} ) e^{ - i ( \e{q} + \hOms) \tau_{56} } ] \times \qty[ \theta( \tau_{56} ) e^{- i ( \e{p} - \hOms ) \tau_{56} } + \theta( - \tau_{56} ) e^{ - i ( \e{q} + \hOms) \tau_{56} } ]
@ -419,12 +420,12 @@ The $GW$ quasiparticle energies $\eGW{p}$ are good approximations to the removal
%Neglecting the $Y_{jb}^{s}$ weights leads to the Tamm-Dancoff approximation (TDA). %Neglecting the $Y_{jb}^{s}$ weights leads to the Tamm-Dancoff approximation (TDA).
%Working out similar expressions for $\mel{N}{T [\hpsi(5) \hpsi^{\dagger}(5)] }{N,s}$ and $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1'})] }{N,s}$, %Working out similar expressions for $\mel{N}{T [\hpsi(5) \hpsi^{\dagger}(5)] }{N,s}$ and $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1'})] }{N,s}$,
Substituting Eq.~\eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, working out similar expressions for the remaining terms, and projecting onto $\phi_a^*(\bx_1) \phi_i(\bx_{1'})$, one gets after a few tedious manipulations (see {\SI}) the dynamical BSE (dBSE): Substituting Eq.~\eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, working out similar expressions for the remaining terms, and projecting onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$, one gets after a few tedious manipulations (see {\SI}) the dynamical BSE (dBSE):
\begin{equation} \label{eq:BSE-final} \begin{equation} \label{eq:BSE-final}
\begin{split} \begin{split}
( \eGW{a} - \eGW{i} - \Oms ) X_{ia}^{s} ( \eGW{a} - \eGW{i} - \Oms ) X_{ia}^{s}
& + \sum_{jb} \qty[ (ia|jb) - \widetilde{W}_{ij,ab}(\Oms) ] X_{jb}^{s} \\ & + \sum_{jb} \qty[ \ERI{ia}{jb} - \widetilde{W}_{ij,ab}(\Oms) ] X_{jb}^{s} \\
& + \sum_{jb} \qty[ (ia|bj) - \widetilde{W}_{ib,aj}(\Oms) ] Y_{jb}^{s} & + \sum_{jb} \qty[ \ERI{ia}{bj} - \widetilde{W}_{ib,aj}(\Oms) ] Y_{jb}^{s}
= 0, = 0,
\end{split} \end{split}
\end{equation} \end{equation}
@ -432,9 +433,9 @@ with $X_{jb}^{s} = \mel{N}{\ha_j^{\dagger} \ha_b}{N,s}$ and $Y_{jb}^{s} = \mel{N
Neglecting the term $Y_{jb}^{s}$ in the dBSE, which is much smaller than $X_{jb}^{s}$, leads to the well-known Tamm-Dancoff approximation (TDA). Neglecting the term $Y_{jb}^{s}$ in the dBSE, which is much smaller than $X_{jb}^{s}$, leads to the well-known Tamm-Dancoff approximation (TDA).
In Eq.~\eqref{eq:BSE-final}, In Eq.~\eqref{eq:BSE-final},
\begin{equation} \begin{equation}
(pq|rs) = \iint d\br d\br' \, \phi_p^*(\br) \phi_q(\br) v(\br -\br') \phi_r^*(\br') \phi_s(\br'), \ERI{pq}{rs} = \iint d\br d\br' \, \MO{p}^*(\br) \MO{q}(\br) v(\br -\br') \MO{r}^*(\br') \MO{s}(\br'),
\end{equation} \end{equation}
are the bare two-electron integrals, $\phi_p(\br)$ is a spatial orbital, and are the bare two-electron integrals in the spatial orbital basis $\lbrace \MO{p}(\br{}) \rbrace$, and
\begin{multline} \begin{multline}
\widetilde{W}_{ij,ab}(\Oms) \widetilde{W}_{ij,ab}(\Oms)
= \frac{ i }{ 2 \pi} \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega) = \frac{ i }{ 2 \pi} \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega)
@ -444,124 +445,26 @@ are the bare two-electron integrals, $\phi_p(\br)$ is a spatial orbital, and
is an effective dynamically-screened Coulomb potential, \cite{Romaniello_2009b} where $\Om{pq}{s} = \Oms - ( \eGW{q} - \eGW{p} )$ and is an effective dynamically-screened Coulomb potential, \cite{Romaniello_2009b} where $\Om{pq}{s} = \Oms - ( \eGW{q} - \eGW{p} )$ and
\begin{equation} \begin{equation}
W_{pq,rs}({\omega}) W_{pq,rs}({\omega})
= \iint d\br d\br' \, \phi_p(\br) \phi_q^*(\br) W(\br ,\br'; \omega) \phi_r^*(\br') \phi_s(\br'). = \iint d\br d\br' \, \MO{p}(\br) \MO{q}^*(\br) W(\br ,\br'; \omega) \MO{r}^*(\br') \MO{s}(\br').
\end{equation} \end{equation}
\xavier{A second coupled equation for the $(X_{ia}^{s}, Y_{ia}^{s} )$ vector can be obtained by projecting now onto the $\mel{N}{T \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1'})}{N,s}$ left-hand side and right-hand-side of the BSE, leading to : } \xavier{A second coupled equation for the $(X_{ia}^{s}, Y_{ia}^{s} )$ vector can be obtained by projecting now onto the $\mel{N}{T \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1'})}{N,s}$ left-hand side and right-hand-side of the BSE, leading to : }
\titou{T2: I propose to move what's below in the next section. What do you think Xavier?} In the present study, we consider the exact spectral representation of $W(\omega)$ at the random-phase approximation (RPA) level, which reads
In the present study, we consider the exact spectral representation of $W(\omega)$ at the RPA level, which reads
\begin{multline} \begin{multline}
\label{eq:W}
W_{ij,ab}(\omega) W_{ij,ab}(\omega)
= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m] = \ERI{ij}{ab} + 2 \sum_m \sERI{ij}{m} \sERI{ab}{m}
\\ \\
\times \qty[ \frac{1}{ \omega-\Om{m}{\RPA} + i\eta } - \frac{1}{ \omega + \Om{m}{\RPA} - i\eta } ] \times \qty[ \frac{1}{ \omega-\Om{m}{\RPA} + i\eta } - \frac{1}{ \omega + \Om{m}{\RPA} - i\eta } ]
\end{multline} \end{multline}
where where
\begin{equation} \begin{equation}
\label{eq:sERI} \label{eq:sERI}
\sERI{pq}{m} = \sum_i^{\Nocc} \sum_a^{\Nvir} \ERI{pq}{ia} (\bX{m}{\RPA} + \bY{m}{\RPA})_{ia} \sERI{pq}{m} = \sum_{ia} \ERI{pq}{ia} (\bX{m}{\RPA} + \bY{m}{\RPA})_{ia}
\end{equation} \end{equation}
are the spectral weights. are the spectral weights.
Because $\Om{m}{\RPA} > 0$, we have In Eqs.~\eqref{eq:W} and \eqref{eq:sERI}, $\OmRPA{m}{}$ and $(\bX{m}{\RPA} + \bY{m}{\RPA})$ are RPA neutral excitations and their corresponding transition vectors computed by solving the (static) linear response problem
\begin{multline}
\widetilde{W}_{ij,ab}( \Oms )
= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m]
\\
\times \qty( \frac{1}{\Om{ib}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{s} - \Om{m}{\RPA} + i\eta} )
\end{multline}
\titou{Due to excitonic effects, the lowest BSE excitation energy, ${\Omega}_1$, stands lower than the lowest RPA excitation energy, $\Omega_m^{RPA}$, 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.
$$
\abs{ \frac{1}{\Omega_{ib}^{s} - \Omega_m^{RPA}} } < \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{}(\titou{-}\omega) & -\bA{}(\titou{-}\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} \begin{equation}
\label{eq:LR-RPA} \label{eq:LR-RPA}
\begin{pmatrix} \begin{pmatrix}
@ -591,6 +494,69 @@ with
\end{align} \end{align}
\end{subequations} \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$. 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$.
Because $\Om{m}{\RPA} > 0$, we have
\begin{multline}
\widetilde{W}_{ij,ab}( \Oms )
= \ERI{ij}{ab} + 2 \sum_m^{OV} \sERI{ij}{m} \sERI{ab}{m}
\\
\times \qty( \frac{1}{\Om{ib}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{s} - \Om{m}{\RPA} + i\eta} )
\end{multline}
\titou{Due to excitonic effects, the lowest BSE excitation energy, ${\Omega}_1$, stands lower than the lowest RPA excitation energy, $\Omega_m^{RPA}$, 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.
$$
\abs{ \frac{1}{\Omega_{ib}^{s} - \Omega_m^{RPA}} } < \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. }
%================================
\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{}(\titou{-}\omega) & -\bA{}(\titou{-}\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} - \tW{ij,ab}{}(\omega),
\\
\label{eq:BSE-Bdyn}
\B{ia,jb}{}(\omega) & = 2 \sigma \ERI{ia}{bj} - \tW{ib,aj}{}(\omega),
\end{align}
\end{subequations}
where $\eGW{ia} = \eGW{a} - \eGW{i}$ are occupied-to-virtual differences of $GW$ quasiparticle energies, and $\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}
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 Now, let us decompose, using basic perturbation theory, the non-linear eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static (hence linear) reference and a first-order dynamic (hence non-linear) perturbation, such that
\begin{multline} \begin{multline}
@ -625,16 +591,16 @@ and
\begin{subequations} \begin{subequations}
\begin{align} \begin{align}
\label{eq:BSE-A1} \label{eq:BSE-A1}
\A{ia,jb}{(1)}(\omega) & = - \W{ij,ab}{}(\omega) + \W{ij,ab}{\text{stat}}, \A{ia,jb}{(1)}(\omega) & = - \tW{ij,ab}{}(\omega) + \W{ij,ab}{\text{stat}},
\\ \\
\label{eq:BSE-B1} \label{eq:BSE-B1}
\B{ia,jb}{(1)}(\omega) & = - \W{ib,aj}{}(\omega) + \W{ib,aj}{\text{stat}}, \B{ia,jb}{(1)}(\omega) & = - \tW{ib,aj}{}(\omega) + \W{ib,aj}{\text{stat}},
\end{align} \end{align}
\end{subequations} \end{subequations}
where we have defined the static version of the screened Coulomb potential where we have defined the static version of the screened Coulomb potential
\begin{equation} \begin{equation}
\label{eq:Wstat} \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}. \W{ij,ab}{\text{stat}} = W_{ij,ab}(\omega = 0) = \ERI{ij}{ab} - 4 \sum_m^{\Nocc \Nvir} \frac{\sERI{ij}{m} \sERI{ab}{m}}{\OmRPA{m}{} - i \eta}.
\end{equation} \end{equation}
According to perturbation theory, the $m$th BSE excitation energy and its corresponding eigenvector can then decomposed as According to perturbation theory, the $m$th BSE excitation energy and its corresponding eigenvector can then decomposed as
\begin{subequations} \begin{subequations}