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

saving work: merging parts

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-05-27 17:38:58 +02:00
parent 6bb4bd5235
commit f2e09d1532
1 changed files with 80 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

155
BSEdyn.tex
View File

@ -180,8 +180,6 @@
\newcommand{\pis}{\pi^*} \newcommand{\pis}{\pi^*}
\newcommand{\ra}{\rightarrow} \newcommand{\ra}{\rightarrow}
\newcommand\Oms{{\Omega}_s}
\newcommand\hOms{\frac{{\Omega}_s}{2}}
\newcommand{\hpsi}{\Hat{\psi}} \newcommand{\hpsi}{\Hat{\psi}}
\newcommand{\ha}{\Hat{a}} \newcommand{\ha}{\Hat{a}}
\newcommand{\tchi}{\Tilde{\chi}} \newcommand{\tchi}{\Tilde{\chi}}
@ -279,7 +277,7 @@ Additional details about this derivation are provided as {\SI}.
We present, in a second step, the perturbative implementation of the dynamical correction as compared to the standard static approximation. We present, in a second step, the perturbative implementation of the dynamical correction as compared to the standard static approximation.
%================================ %================================
\subsection{General dynamical BSE theory} \subsection{General dynamical BSE}
%================================= %=================================
The two-particle correlation function $L(1,2; 1',2')$ --- a central quantity in the BSE formalism --- relates the variation of the one-body Green's function $G(1,1')$ with respect to an external non-local perturbation $U(2',2)$, \ie, The two-particle correlation function $L(1,2; 1',2')$ --- a central quantity in the BSE formalism --- relates the variation of the one-body Green's function $G(1,1')$ with respect to an external non-local perturbation $U(2',2)$, \ie,
@ -332,9 +330,9 @@ In the optical limit of instantaneous electron-hole creation and destruction, im
\begin{equation} \begin{equation}
\begin{split} \begin{split}
iL(1,2; 1',2') iL(1,2; 1',2')
& = \theta(+\tau_{12}) \sum_{s > 0} \chi_s(\bx_1,\bx_{1'}) \tchi_s(\bx_2,\bx_{2'}) e^{ - i \Oms \tau_{12} } & = \theta(+\tau_{12}) \sum_{s > 0} \chi_s(\bx_1,\bx_{1'}) \tchi_s(\bx_2,\bx_{2'}) e^{ - i \Om{s}{} \tau_{12} }
\\ \\
& - \theta(-\tau_{12}) \sum_{s > 0} \chi_s(\bx_2,\bx_{2'}) \tchi_s(\bx_1,\bx_{1'}) e^{ + i \Oms \tau_{12} }, & - \theta(-\tau_{12}) \sum_{s > 0} \chi_s(\bx_2,\bx_{2'}) \tchi_s(\bx_1,\bx_{1'}) e^{ + i \Om{s}{} \tau_{12} },
\end{split} \end{split}
\end{equation} \end{equation}
where $\tau_{12} = t_1 - t_2$, $\theta$ is the Heaviside step function, and where $\tau_{12} = t_1 - t_2$, $\theta$ is the Heaviside step function, and
@ -345,15 +343,15 @@ where $\tau_{12} = t_1 - t_2$, $\theta$ is the Heaviside step function, and
\tchi_s(\bx_1,\bx_{2}) & = \mel{N,s}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N}. \tchi_s(\bx_1,\bx_{2}) & = \mel{N,s}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N}.
\end{align} \end{align}
\end{subequations} \end{subequations}
The $\Oms$'s are the neutral excitation energies of interest. We have used the relation between the field operators in their time-dependent (Heisenberg) and time-independent (Schr\"{o}dinger) representations, e.g. The $\Om{s}{}$'s are the neutral excitation energies of interest. We have used the relation between the field operators in their time-dependent (Heisenberg) and time-independent (Schr\"{o}dinger) representations, e.g.
$$ $$
\hpsi(1) = e^{ i {\hat H} t_1 } \hpsi(\bx_1) e^{-i {\hat H} t_1 } \hpsi(1) = e^{ i {\hat H} t_1 } \hpsi(\bx_1) e^{-i {\hat H} t_1 }
$$ $$
with $\hat H$ the exact many-body Hamiltonian. with $\hat H$ the exact many-body Hamiltonian.
Picking up the $e^{+i \Oms t_2 }$ component in $L(1,2; 1',2')$ and $L(6,2;5,2')$, simplifying further by $\tchi_s(\bx_2,\bx_{2'})$ on both side of the BSE [see Eq.~\eqref{eq:BSE}], we are left with the search of the $e^{-i \Oms t_1 }$ Fourier component associated with the right-hand side of a modified dynamical BSE, which reads Picking up the $e^{+i \Om{s}{} t_2 }$ component in $L(1,2; 1',2')$ and $L(6,2;5,2')$, simplifying further by $\tchi_s(\bx_2,\bx_{2'})$ on both side of the BSE [see Eq.~\eqref{eq:BSE}], we are left with the search of the $e^{-i \Om{s}{} t_1 }$ Fourier component associated with the right-hand side of a modified dynamical BSE, which reads
\begin{multline} \label{eq:BSE_2} \begin{multline} \label{eq:BSE_2}
\mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,s} e^{ - i \Oms t_1 } \mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,s} e^{ - i \Om{s}{} t_1 }
\theta ( \tau_{12} ) \theta ( \tau_{12} )
\\ \\
= \int d3456 \, L_0(1,4;1',3) \Xi(3,5;4,6) = \int d3456 \, L_0(1,4;1',3) \Xi(3,5;4,6)
@ -361,7 +359,7 @@ Picking up the $e^{+i \Oms t_2 }$ component in $L(1,2; 1',2')$ and $L(6,2;5,2')$
\times \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s} \times \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
\theta [\min(t_5,t_6) - t_2]. \theta [\min(t_5,t_6) - t_2].
\end{multline} \end{multline}
For the lowest neutral excitation energies falling in the fundamental gap of the system (\ie, $\Oms < \EgFun$ due to excitonic effects), $L_0(1,2;1',2')$ cannot contribute to the $e^{-i \Oms t_1 }$ response since its lowest excitation energy is precisely the fundamental gap [see Eq.~\eqref{eq:Egfun}]. For the lowest neutral excitation energies falling in the fundamental gap of the system (\ie, $\Om{s}{} < \EgFun$ due to excitonic effects), $L_0(1,2;1',2')$ cannot contribute to the $e^{-i \Om{s}{} t_1 }$ response since its lowest excitation energy is precisely the fundamental gap [see Eq.~\eqref{eq:Egfun}].
Consequently, special care has to be taken for high-lying excited states (like core or Rydberg excitations) where additional terms have to be taken into account (see Refs.~\onlinecite{Strinati_1982,Strinati_1984}). Consequently, special care has to be taken for high-lying excited states (like core or Rydberg excitations) where additional terms have to be taken into account (see Refs.~\onlinecite{Strinati_1982,Strinati_1984}).
Dropping the (space/spin) variables, the Fourier components with respect to $t_1$ of $L_0(1,4;1',3)$ reads Dropping the (space/spin) variables, the Fourier components with respect to $t_1$ of $L_0(1,4;1',3)$ reads
@ -384,13 +382,13 @@ 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 $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$ yields Projecting the Fourier component $L_0(\bx_1,4;\bx_{1'},3; \omega_1 = \Om{s}{} )$ onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$ yields
\begin{multline} \label{eq:iL0bis} \begin{multline} \label{eq:iL0bis}
\iint d\bx_1 d\bx_{1'} \, \MO{a}^*(\bx_1) \MO{i}(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Oms) \iint d\bx_1 d\bx_{1'} \, \MO{a}^*(\bx_1) \MO{i}(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Om{s}{})
\\ \\
= =
\frac{ \MO{a}^*(\bx_3) \MO{i}(\bx_4) e^{i \Oms t^{34} }} { \Oms - ( \e{a} - \e{i} ) + i \eta } \frac{ \MO{a}^*(\bx_3) \MO{i}(\bx_4) e^{i \Om{s}{} t^{34} }} { \Om{s}{} - ( \e{a} - \e{i} ) + i \eta }
\qty[ \theta( \tau_{34} ) e^{i ( \e{i} + \hOms) \tau_{34} } + \theta( - \tau_{34} ) e^{i (\e{a} - \hOms ) \tau_{34} } ]. \qty[ \theta( \tau_{34} ) e^{i ( \e{i} + \frac{\Om{s}{}}{2}) \tau_{34} } + \theta( - \tau_{34} ) e^{i (\e{a} - \frac{\Om{s}{}}{2}) \tau_{34} } ].
\end{multline} \end{multline}
% and $(i,j)$/$(a,b)$ index occupied/virtual orbitals, respectively. % and $(i,j)$/$(a,b)$ index occupied/virtual orbitals, respectively.
As a final step, we express the terms $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}')] }{N,s}$ and $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}$ from Eq.~\eqref{eq:BSE_2} in the standard electron-hole product (or single-excitation) space. As a final step, we express the terms $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}')] }{N,s}$ and $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}$ from Eq.~\eqref{eq:BSE_2} in the standard electron-hole product (or single-excitation) space.
@ -403,11 +401,15 @@ For example, we have
= - \qty( e^{ -i \Omega_s t^{65} } ) \sum_{pq} \MO{p}(\bx_6) \MO{q}^*(\bx_5) = - \qty( e^{ -i \Omega_s t^{65} } ) \sum_{pq} \MO{p}(\bx_6) \MO{q}^*(\bx_5)
\mel{N}{\ha_q^{\dagger} \ha_p}{N,s} \mel{N}{\ha_q^{\dagger} \ha_p}{N,s}
\\ \\
\times \qty[ \theta( \tau_{65} ) e^{- i ( \e{p} - \hOms ) \tau_{65} } + \theta( - \tau_{65} ) e^{ - i ( \e{q} + \hOms) \tau_{65} } ] \times \qty[ \theta( \tau_{65} ) e^{- i ( \e{p} - \frac{\Om{s}{}}{2} ) \tau_{65} } + \theta( - \tau_{65} ) e^{ - i ( \e{q} + \frac{\Om{s}{}}{2}) \tau_{65} } ]
\end{multline} \end{multline}
with $t^{65} = (t_5 + t_6)/2$ and $\tau_{65} = t_6 -t_5$. with $t^{65} = (t_5 + t_6)/2$ and $\tau_{65} = t_6 -t_5$.
%with a similar expression for $\mel{N}{T [\hpsi(\bx_3) \hpsi^{\dagger}(\bx_4)] }{N,s}$. %with a similar expression for $\mel{N}{T [\hpsi(\bx_3) \hpsi^{\dagger}(\bx_4)] }{N,s}$.
%================================
\subsection{Dynamical BSE within the $GW$ approximation}
%=================================
Adopting now the $GW$ approximation for the exchange-correlation self-energy, \ie, Adopting now the $GW$ approximation for the exchange-correlation self-energy, \ie,
\begin{equation} \begin{equation}
\Sigma_\text{xc}^{\GW}(1,2) = i G(1,2) W(1^+,2), \Sigma_\text{xc}^{\GW}(1,2) = i G(1,2) W(1^+,2),
@ -427,9 +429,9 @@ The $GW$ quasiparticle energies $\eGW{p}$ are good approximations to the removal
Substituting Eqs.~\eqref{eq:iL0bis},\eqref{eq:spectral65},\eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, and projecting onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$, one gets after a few tedious manipulations (see {\SI}) the dynamical BSE (dBSE): Substituting Eqs.~\eqref{eq:iL0bis},\eqref{eq:spectral65},\eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, and projecting onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$, one gets after a few tedious manipulations (see {\SI}) the dynamical BSE (dBSE):
\begin{equation} \label{eq:BSE-final} \begin{equation} \label{eq:BSE-final}
\begin{split} \begin{split}
( \eGW{a} - \eGW{i} - \Oms ) X_{ia}^{s} ( \eGW{a} - \eGW{i} - \Om{s}{} ) X_{ia}^{s}
& + \sum_{jb} \qty[ \ERI{ia}{jb} - \widetilde{W}_{ij,ab}(\Oms) ] X_{jb}^{s} \\ & + \sum_{jb} \qty[ \ERI{ia}{jb} - \widetilde{W}_{ij,ab}(\Om{s}{}) ] X_{jb}^{s} \\
& + \sum_{jb} \qty[ \ERI{ia}{bj} - \widetilde{W}_{ib,aj}(\Oms) ] Y_{jb}^{s} & + \sum_{jb} \qty[ \ERI{ia}{bj} - \widetilde{W}_{ib,aj}(\Om{s}{}) ] Y_{jb}^{s}
= 0, = 0,
\end{split} \end{split}
\end{equation} \end{equation}
@ -441,12 +443,12 @@ In Eq.~\eqref{eq:BSE-final},
\end{equation} \end{equation}
are the bare two-electron integrals in the spatial orbital basis $\lbrace \MO{p}(\br{}) \rbrace$, 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}(\Om{s}{})
= \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)
\\ \\
\times \qty[ \frac{1}{ \Omega_{ib}^s - \omega + i \eta } + \frac{1}{ \Omega_{ja}^{s} + \omega + i\eta } ], \times \qty[ \frac{1}{ \Omega_{ib}^s - \omega + i \eta } + \frac{1}{ \Omega_{ja}^{s} + \omega + i\eta } ],
\end{multline} \end{multline}
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} = \Om{s}{} - ( \eGW{q} - \eGW{p} )$ and
\begin{equation} \begin{equation}
W_{pq,rs}({\omega}) W_{pq,rs}({\omega})
= \iint d\br d\br' \, \MO{p}(\br) \MO{q}^*(\br) W(\br ,\br'; \omega) \MO{r}^*(\br') \MO{s}(\br'). = \iint d\br d\br' \, \MO{p}(\br) \MO{q}^*(\br) W(\br ,\br'; \omega) \MO{r}^*(\br') \MO{s}(\br').
@ -454,13 +456,18 @@ is an effective dynamically-screened Coulomb potential, \cite{Romaniello_2009b}
\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 : }
%================================
\subsection{Dynamical screening}
%=================================
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 random-phase approximation (RPA) level, which reads
\begin{multline} \begin{multline}
\label{eq:W} \label{eq:W}
W_{ij,ab}(\omega) W_{ij,ab}(\omega)
= \ERI{ij}{ab} + 2 \sum_m \sERI{ij}{m} \sERI{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}
@ -498,20 +505,19 @@ 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$.
The RPA matrices $\bA{\RPA}$ and $\bB{\RPA}$ in Eq.~\eqref{eq:LR-RPA} are of size $\Nocc \Nvir \times \Nocc \Nvir$, where $\Nocc$ and $\Nvir$ are the number of occupied and virtual orbitals (\ie, $\Norb = \Nocc + \Nvir$ is the total number of spatial orbitals), respectively, and $\bX{m}{}$, and $\bY{m}{}$ are (eigen)vectors of length $\Nocc \Nvir$.
Because $\Om{m}{\RPA} > 0$, we have Because $\Om{m}{\RPA} > 0$, we have
\begin{multline} \begin{multline}
\widetilde{W}_{ij,ab}( \Oms ) \widetilde{W}_{ij,ab}( \Om{s}{} )
= \ERI{ij}{ab} + 2 \sum_m^{OV} \sERI{ij}{m} \sERI{ab}{m} = \ERI{ij}{ab} + 2 \sum_m \sERI{ij}{m} \sERI{ab}{m}
\\ \\
\times \qty( \frac{1}{\Om{ib}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{s} - \Om{m}{\RPA} + i\eta} ) \times \qty[ \frac{1}{\Om{ib}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{s} - \Om{m}{\RPA} + i\eta} ].
\end{multline} \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 Due to excitonic effects, the lowest BSE excitation energy, $\Om{1}{}$, stands lower than the lowest RPA excitation energy, $\Om{1}{\RPA}$, so that, $\Om{ib}{s} - \Om{m}{\RPA} < 0 $ and $\widetilde{W}_{ij,ab}(\Om{s}{})$ has no resonances.
e.g. $( \Omega_{ib}^{s} - \Omega_m^{RPA} )$ is strictly negative and the $\widetilde{W}_{ij,ab}( \Oms )$ present no resonances. Furthermore, $\Om{ib}{s}$ and $\Om{ja}{s}$ are necessarily negative quantities for in-gap low-lying BSE excitations.
Further, $\Omega_{ib}^{s}$ and $\Omega_{ja}^{s}$ are necessarily negative for in-gap low lying BSE excitations, such that e.g. Thus, we have $\abs*{\Omega_{ib}^{s} - \Om{m}{\RPA}} > \Omega_m^{\RPA}$.
$$ As a consequence, we observe a reduction of the electron-hole screening, \ie, an enhancement of electron-hole stabilizing binding energy, as compared to the standard static BSE, and yields smaller (red-shifted) excitation energies.
\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} \subsection{Perturbative dynamical correction}
@ -521,36 +527,35 @@ For a closed-shell system in a finite basis, to compute the BSE excitation energ
\begin{equation} \begin{equation}
\label{eq:LR-dyn} \label{eq:LR-dyn}
\begin{pmatrix} \begin{pmatrix}
\bA{}(\omega) & \bB{}(\omega) \\ \bA{}(\Om{s}{}) & \bB{}(\Om{s}{}) \\
-\bB{}(\titou{-}\omega) & -\bA{}(\titou{-}\omega) \\ -\bB{}(\titou{-}\Om{s}{}) & -\bA{}(\titou{-}\Om{s}{}) \\
\end{pmatrix} \end{pmatrix}
\cdot \cdot
\begin{pmatrix} \begin{pmatrix}
\bX{m}{}(\omega) \\ \bX{s}{} \\
\bY{m}{}(\omega) \\ \bY{s}{} \\
\end{pmatrix} \end{pmatrix}
= =
\omega \Om{s}{}
\begin{pmatrix} \begin{pmatrix}
\bX{m}{}(\omega) \\ \bX{s}{} \\
\bY{m}{}(\omega) \\ \bY{s}{} \\
\end{pmatrix}, \end{pmatrix},
\end{equation} \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$. where the dynamical matrices $\bA{}$ and $\bB{}$, as well as $\bX{s}{}$, and $\bY{s}{}$, have the same size as their RPA counterparts.
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 $s$. \cite{Romaniello_2009b,Sangalli_2011,Martin_2016}
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 The BSE matrix elements read
\begin{subequations} \begin{subequations}
\begin{align} \begin{align}
\label{eq:BSE-Adyn} \label{eq:BSE-Adyn}
\A{ia,jb}{}(\omega) & = \delta_{ij} \delta_{ab} \eGW{ia} + 2 \sigma \ERI{ia}{jb} - \tW{ij,ab}{}(\omega), \A{ia,jb}{}(\Om{s}{}) & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + 2 \sigma \ERI{ia}{jb} - \tW{ij,ab}{}(\Om{s}{}),
\\ \\
\label{eq:BSE-Bdyn} \label{eq:BSE-Bdyn}
\B{ia,jb}{}(\omega) & = 2 \sigma \ERI{ia}{bj} - \tW{ib,aj}{}(\omega), \B{ia,jb}{}(\Om{s}{}) & = 2 \sigma \ERI{ia}{bj} - \tW{ib,aj}{}(\Om{s}{}),
\end{align} \end{align}
\end{subequations} \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). where $\sigma = 1 $ or $0$ for singlet and triplet excited states (respectively).
%\begin{equation} %\begin{equation}
% \ERI{pq}{rs} = \iint \frac{\MO{p}(\br{}) \MO{q}(\br{}) \MO{r}(\br{}') \MO{s}(\br{}')}{\abs*{\br{} - \br{}'}} \dbr{} \dbr{}' % \ERI{pq}{rs} = \iint \frac{\MO{p}(\br{}) \MO{q}(\br{}) \MO{r}(\br{}') \MO{s}(\br{}')}{\abs*{\br{} - \br{}'}} \dbr{} \dbr{}'
%\end{equation} %\end{equation}
@ -566,8 +571,8 @@ Now, let us decompose, using basic perturbation theory, the non-linear eigenprob
\begin{multline} \begin{multline}
\label{eq:LR-PT} \label{eq:LR-PT}
\begin{pmatrix} \begin{pmatrix}
\bA{}(\omega) & \bB{}(\omega) \\ \bA{}(\Om{s}{}) & \bB{}(\Om{s}{}) \\
-\bB{}(\titou{-}\omega) & -\bA{}(\titou{-}\omega) \\ -\bB{}(\titou{-}\Om{s}{}) & -\bA{}(\titou{-}\Om{s}{}) \\
\end{pmatrix} \end{pmatrix}
\\ \\
= =
@ -577,8 +582,8 @@ Now, let us decompose, using basic perturbation theory, the non-linear eigenprob
\end{pmatrix} \end{pmatrix}
+ +
\begin{pmatrix} \begin{pmatrix}
\bA{(1)}(\omega) & \bB{(1)}(\omega) \\ \bA{(1)}(\Om{s}{}) & \bB{(1)}(\Om{s}{}) \\
-\bB{(1)}(\titou{-}\omega) & -\bA{(1)}(\titou{-}\omega) \\ -\bB{(1)}(\titou{-}\Om{s}{}) & -\bA{(1)}(\titou{-}\Om{s}{}) \\
\end{pmatrix} \end{pmatrix}
\end{multline} \end{multline}
with with
@ -595,35 +600,35 @@ and
\begin{subequations} \begin{subequations}
\begin{align} \begin{align}
\label{eq:BSE-A1} \label{eq:BSE-A1}
\A{ia,jb}{(1)}(\omega) & = - \tW{ij,ab}{}(\omega) + \W{ij,ab}{\text{stat}}, \A{ia,jb}{(1)}(\Om{s}{}) & = - \tW{ij,ab}{}(\Om{s}{}) + \W{ij,ab}{\text{stat}},
\\ \\
\label{eq:BSE-B1} \label{eq:BSE-B1}
\B{ia,jb}{(1)}(\omega) & = - \tW{ib,aj}{}(\omega) + \W{ib,aj}{\text{stat}}, \B{ia,jb}{(1)}(\Om{s}{}) & = - \tW{ib,aj}{}(\Om{s}{}) + \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}} = W_{ij,ab}(\omega = 0) = \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 \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 $s$th BSE excitation energy and its corresponding eigenvector can then decomposed as
\begin{subequations} \begin{subequations}
\begin{gather} \begin{gather}
\Om{m}{} = \Om{m}{(0)} + \Om{m}{(1)} + \ldots, \Om{s}{} = \Om{s}{(0)} + \Om{s}{(1)} + \ldots,
\\ \\
\begin{pmatrix} \begin{pmatrix}
\bX{m}{} \\ \bX{s}{} \\
\bY{m}{} \\ \bY{s}{} \\
\end{pmatrix} \end{pmatrix}
= =
\begin{pmatrix} \begin{pmatrix}
\bX{m}{(0)} \\ \bX{s}{(0)} \\
\bY{m}{(0)} \\ \bY{s}{(0)} \\
\end{pmatrix} \end{pmatrix}
+ +
\begin{pmatrix} \begin{pmatrix}
\bX{m}{(1)} \\ \bX{s}{(1)} \\
\bY{m}{(1)} \\ \bY{s}{(1)} \\
\end{pmatrix} \end{pmatrix}
+ \ldots. + \ldots.
\end{gather} \end{gather}
@ -637,50 +642,50 @@ Solving the zeroth-order static problem yields
\end{pmatrix} \end{pmatrix}
\cdot \cdot
\begin{pmatrix} \begin{pmatrix}
\bX{m}{(0)} \\ \bX{s}{(0)} \\
\bY{m}{(0)} \\ \bY{s}{(0)} \\
\end{pmatrix} \end{pmatrix}
= =
\Om{m}{(0)} \Om{s}{(0)}
\begin{pmatrix} \begin{pmatrix}
\bX{m}{(0)} \\ \bX{s}{(0)} \\
\bY{m}{(0)} \\ \bY{s}{(0)} \\
\end{pmatrix}, \end{pmatrix},
\end{equation} \end{equation}
and, thanks to first-order perturbation theory, the first-order correction to the $m$th excitation energy is and, thanks to first-order perturbation theory, the first-order correction to the $s$th excitation energy is
\begin{equation} \begin{equation}
\label{eq:Om1} \label{eq:Om1}
\Om{m}{(1)} = \Om{s}{(1)} =
\T{\begin{pmatrix} \T{\begin{pmatrix}
\bX{m}{(0)} \\ \bX{s}{(0)} \\
\bY{m}{(0)} \\ \bY{s}{(0)} \\
\end{pmatrix}} \end{pmatrix}}
\cdot \cdot
\begin{pmatrix} \begin{pmatrix}
\bA{(1)}(\Om{m}{(0)}) & \bB{(1)}(\Om{m}{(0)}) \\ \bA{(1)}(\Om{s}{(0)}) & \bB{(1)}(\Om{s}{(0)}) \\
-\bB{(1)}(\titou{-}\Om{m}{(0)}) & -\bA{(1)}(\titou{-}\Om{m}{(0)}) \\ -\bB{(1)}(\titou{-}\Om{s}{(0)}) & -\bA{(1)}(\titou{-}\Om{s}{(0)}) \\
\end{pmatrix} \end{pmatrix}
\cdot \cdot
\begin{pmatrix} \begin{pmatrix}
\bX{m}{(0)} \\ \bX{s}{(0)} \\
\bY{m}{(0)} \\ \bY{s}{(0)} \\
\end{pmatrix}. \end{pmatrix}.
\end{equation} \end{equation}
From a practical point of view, if one enforces the TDA, we obtain the very simple expression From a practical point of view, if one enforces the TDA, we obtain the very simple expression
\begin{equation} \begin{equation}
\label{eq:Om1-TDA} \label{eq:Om1-TDA}
\Om{m}{(1)} = \T{(\bX{m}{(0)})} \cdot \bA{(1)}(\Om{m}{(0)}) \cdot \bX{m}{(0)}. \Om{s}{(1)} = \T{(\bX{s}{(0)})} \cdot \bA{(1)}(\Om{s}{(0)}) \cdot \bX{s}{(0)}.
\end{equation} \end{equation}
This correction can be renormalized by computing, at basically no extra cost, the renormalization factor which reads, in the TDA, This correction can be renormalized by computing, at basically no extra cost, the renormalization factor which reads, in the TDA,
\begin{equation} \begin{equation}
\label{eq:Z} \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}. Z_{s} = \qty[ 1 - \T{(\bX{s}{(0)})} \cdot \left. \pdv{\bA{(1)}(\Om{s}{})}{\Om{s}{}} \right|_{\Om{s}{} = \Om{s}{(0)}} \cdot \bX{s}{(0)} ]^{-1}.
\end{equation} \end{equation}
This finally yields This finally yields
\begin{equation} \begin{equation}
\Om{m}{\text{dyn}} = \Om{m}{\text{stat}} + \Delta\Om{m}{\text{dyn}} = \Om{m}{(0)} + Z_{m} \Om{m}{(1)}. \Om{s}{\text{dyn}} = \Om{s}{\text{stat}} + \Delta\Om{s}{\text{dyn}} = \Om{s}{(0)} + Z_{s} \Om{s}{(1)}.
\end{equation} \end{equation}
with $\Om{m}{\text{stat}} \equiv \Om{m}{(0)}$ and $\Delta\Om{m}{\text{dyn}} = Z_{m} \Om{m}{(1)}$. with $\Om{s}{\text{stat}} \equiv \Om{s}{(0)}$ and $\Delta\Om{s}{\text{dyn}} = Z_{s} \Om{s}{(1)}$.
This is our final expression. This is our final expression.
%%% FIG 1 %%% %%% FIG 1 %%%