PTEROSOR/BSEdyn
10
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
9be0f759f7
BSEdyn/BSEdyn.tex
Pierre-Francois Loos 9be0f759f7 Xav push
2020-05-18 20:00:24 +02:00

565 lines
22 KiB
TeX
Raw Blame History

\documentclass[aps,prb,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
\usepackage[version=4]{mhchem}
%\usepackage{natbib}
%\usepackage[extra]{tipa}
%\bibliographystyle{achemso}
%\AtBeginDocument{\nocite{achemso-control}}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}
\usepackage[
colorlinks=true,
citecolor=blue,
breaklinks=true
]{hyperref}
\urlstyle{same}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{HTML}{009900}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\xavier}[1]{\textcolor{darkgreen}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\XB}[1]{\xavier{(\underline{\bf XB}: #1)}}
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\SI}{\textcolor{blue}{supporting information}}
\newcommand{\QP}{\textsc{quantum package}}
\newcommand{\T}[1]{#1^{\intercal}}
% coordinates
\newcommand{\br}[1]{\mathbf{r}_{#1}}
\newcommand{\dbr}[1]{d\br{#1}}
% methods
\newcommand{\evGW}{ev$GW$}
\newcommand{\qsGW}{qs$GW$}
\newcommand{\GOWO}{$G_0W_0$}
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\xc}{\text{xc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\co}{\text{x}}
%
\newcommand{\Norb}{N}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}
\newcommand{\IS}{\lambda}
% operators
\newcommand{\hH}{\Hat{H}}
% methods
\newcommand{\RPA}{\text{RPA}}
\newcommand{\BSE}{\text{BSE}}
% 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}[2]{[#1|#2]}
%% bold in Table
\newcommand{\bb}[1]{\textbf{#1}}
\newcommand{\rb}[1]{\textbf{\textcolor{red}{#1}}}
\newcommand{\gb}[1]{\textbf{\textcolor{darkgreen}{#1}}}
% excitation energies
\newcommand{\OmRPA}[1]{\Omega_{#1}^{\text{RPA}}}
\newcommand{\OmRPAx}[1]{\Omega_{#1}^{\text{RPAx}}}
\newcommand{\OmBSE}[1]{\Omega_{#1}^{\text{BSE}}}
\newcommand{\spinup}{\downarrow}
\newcommand{\spindw}{\uparrow}
\newcommand{\singlet}{\uparrow\downarrow}
\newcommand{\triplet}{\uparrow\uparrow}
% Matrices
\newcommand{\bO}{\mathbf{0}}
\newcommand{\bI}{\mathbf{1}}
\newcommand{\bvc}{\mathbf{v}}
\newcommand{\bSig}{\mathbf{\Sigma}}
\newcommand{\bSigX}{\mathbf{\Sigma}^\text{x}}
\newcommand{\bSigC}{\mathbf{\Sigma}^\text{c}}
\newcommand{\bSigGW}{\mathbf{\Sigma}^{GW}}
\newcommand{\be}{\mathbf{\epsilon}}
\newcommand{\beGW}{\mathbf{\epsilon}^{GW}}
\newcommand{\beGnWn}[1]{\mathbf{\epsilon}^\text{\GnWn{#1}}}
\newcommand{\bde}{\mathbf{\Delta\epsilon}}
\newcommand{\bdeHF}{\mathbf{\Delta\epsilon}^\text{HF}}
\newcommand{\bdeGW}{\mathbf{\Delta\epsilon}^{GW}}
\newcommand{\bOm}[1]{\mathbf{\Omega}^{#1}}
\newcommand{\bA}[1]{\mathbf{A}^{#1}}
\newcommand{\btA}[1]{\Tilde{\mathbf{A}}^{#1}}
\newcommand{\bB}[1]{\mathbf{B}^{#1}}
\newcommand{\bX}[2]{\mathbf{X}_{#1}^{#2}}
\newcommand{\bY}[2]{\mathbf{Y}_{#1}^{#2}}
\newcommand{\bZ}[2]{\mathbf{Z}_{#1}^{#2}}
\newcommand{\bK}{\mathbf{K}}
\newcommand{\bP}[1]{\mathbf{P}^{#1}}
% units
\newcommand{\IneV}[1]{#1 eV}
\newcommand{\InAU}[1]{#1 a.u.}
\newcommand{\InAA}[1]{#1 \AA}
\newcommand{\kcal}{kcal/mol}
\DeclareMathOperator*{\argmax}{argmax}
\DeclareMathOperator*{\argmin}{argmin}
\newcommand\vari{{\varepsilon}_i}
\newcommand\vara{{\varepsilon}_a}
\newcommand\varj{{\varepsilon}_j}
\newcommand\varb{{\varepsilon}_b}
\newcommand\varn{{\varepsilon}_n}
\newcommand\varm{{\varepsilon}_m}
\newcommand\Oms{{\Omega}_s}
\newcommand\hOms{\frac{{\Omega}_s}{2}}
\newcommand{\NEEL}{Universit\'e Grenoble Alpes, CNRS, Institut NEEL, F-38042 Grenoble, France}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\begin{document}
\title{Dynamical Correction to the Bethe-Salpeter Equation}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Xavier \surname{Blase}}
\email{xavier.blase@neel.cnrs.fr }
\affiliation{\NEEL}
\begin{abstract}
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}
%%%%%%%%%%%%%%%%%%%%%%%%
%================================
\subsection{Theory for physics}
%=================================
The resolution of the dynamical Bethe-Salpeter equation (dBSE) [Strinati]
\begin{align*}
L(1,2; & 1',2') = L_0(1,2;1',2') + \\
&+ \int d3456 \;
L_0(1,4;1',3) \Xi(3,5;4,6) L(6,2;5,2')
\end{align*}
with:
\begin{align*}
iL(1,2; 1',2') &= -G_2(1,2;1',2') + G(1,1')G(2,2') \\
i^2 G_2(1,2;1',2') &= \langle N | T {\hat \psi}(1) {\hat \psi}(2) {\hat \psi}^{\dagger}(2') {\hat \psi}^{\dagger}(1') | N \rangle
\end{align*}
where e.g. $1 = (x_1,t_1)$ a space-spin plus time variable, starts with the expansion of the 2-body Green's function $G_2$ and response function $L$ over the complete orthonormalized set $ |N,s \rangle $ of the N-electron excited state with $| N \rangle = | N,0 \rangle$ the ground-state. In the optical limit of instantaneous electron-hole creation and destruction, imposing
$t_{2'} = t_2^+$ and $t_{1'} = t_1^+$, one obtains:
\begin{align*}
iL(1,2;1',2') &= \theta(\tau_{12}) \sum_{s > 0} \chi_s(x_1,x_{1'}) {\tilde \chi}_s(x_2,x_{2'})
e^{ +i \Oms \tau_{12} } \\
&- \theta(-\tau_{12}) \sum_{s > 0} \chi_s(x_2,x_{2'}) {\tilde \chi}_s(x_1,x_{1'})
e^{ - i \Oms \tau_{12} }
\end{align*}
with $\tau_{12} = t_1 - t_2$ and
\begin{align*}
\chi_s(x_1,x_{1'}) = \langle N | T {\hat \psi}(x_1) {\hat \psi}^{\dagger}(x_{1'}) | N,s \rangle \\
{\tilde \chi}_s(x_2,x_{2'}) = \langle N,s | T {\hat \psi}(x_2) {\hat \psi}^{\dagger}(x_{2'}) | N \rangle
\end{align*}
The $\Oms$ are the neutral excitation energies of interest. Picking up the $e^{+i \Oms t_2 }$ component and simplifying by ${\tilde \chi}_s(x_2,x_{2'})
e^{ i \Oms t_{2} }$ on both side of the Bethe-Salpeter equation, we are left with the search of the $e^{-i \Oms t_1 }$ Fourier component associated with the right-hand side of the BSE. For the lowest $\Oms$ excitation energies falling in the quasiparticle gap of the system, the $L_0(1,2;1',2')$ term cannot contribute since its lowest excitation energy is precisely the quasiparticle gap, namely the difference between the electronic affinity and the ionisation potential.
The Fourier components with respect to time $t_1$ of $iL_0(1, 4; 1', 3) = G(1, 3)G(4, 1')$ reads, dropping the (space/spin)-variables:
\begin{align*}
[iL_0]( \omega_1 ) = \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}$. Adopting now the $GW$ approximation for the exchange-correlation self-energy leads to a simplification of the BSE kernel:
$$
\Xi(3,5;4,6) = v(3,6) \delta(34) \delta(56) - W(3^+,4) \delta(36) \delta(45)
$$
We further obtain the needed spectral representation of
$\langle N | T {\hat \psi}(3) {\hat \psi}^{\dagger}(4) | N,s \rangle$
expanding the field operators over a complete orbital basis creation/destruction operators:
\begin{align*}
\langle N | T {\hat \psi}(3) & {\hat \psi}^{\dagger}(4) | N,s \rangle = - \Big( e^{ -i \Omega_s t^{34} } \Big) \sum_{mn} \phi_m(x_3) \phi_n^*(x_4) \times \nonumber \\
\times & \langle N | {\hat a}_n^{\dagger} {\hat a}_m | N,s \rangle \;\Big[ \theta( \tau ) e^{- i ( \varepsilon_m - \hOms ) \tau } + \theta( -\tau ) e^{ - i ( \varepsilon_n + \hOms) \tau } \Big]
\end{align*}
with $\tau = \tau_{34}$ and where the $ \lbrace \varepsilon_{n/m} \rbrace$ are proper addition/removal energies such that e.g.
$$
e^{ i H \tau } {\hat a}_m^{\dagger} | N \rangle = e^{ i (E_0^N + \varepsilon_m ) \tau } {\hat a} _m^{\dagger} | N \rangle
$$
The $GW$ quasiparticle energies $ \varepsilon_{n/m}^{GW}$ are good approximations to such removal/addition energies. Selecting (n,m)=(j,b) yields the largest components
$A_{jb}^{s} = \langle N | {\hat a}_j^{\dagger} {\hat a}_b | N,s \rangle $, while (n,m)=(b,j) yields much weaker
$B_{jb}^{s} = \langle N | {\hat a}_b^{\dagger} {\hat a}_j | N,s \rangle $ contributions. We used chemist notations with (i,j) indexing occupied orbitals and (a,b) virtual ones. Neglecting the $B_{jb}^{s}$ weights leads to the Tamm Dancoff approximation (TDA). Working out the same expansion for $ \langle N | T {\hat \psi}(5) {\hat \psi}^{\dagger}(5) | N,s \rangle$ and $ \langle N | T {\hat \psi}(x_1) {\hat \psi}^{\dagger}(x_{1'}) | N,s \rangle$, and projecting onto $\phi_a^*(x_1) \phi_i(x_{1'})$,
one obtains after a few tedious manipulations (see Supplemental Information) the dynamical Bethe-Salpeter equation (dBSE) :
\begin{align}
( \varepsilon_a - \varepsilon_i - \Omega_s ) A_{ia}^{s}
&+ \sum_{jb} \Big( v_{ai,bj} - \widetilde{W}_{ij,ab}(\Oms) \Big) A_{jb}^{s} \\
&+ \sum_{bj} \Big( v_{ai,jb} - \widetilde{W}_{ib,aj}(\Oms) \Big) B_{jb}^{s}
= 0
\end{align}
with an effective dynamically screened Coulomb potential (see Pina eq. 24):
\begin{align}
\widetilde{W}_{ij,ab}(\Oms) &= \frac{ i }{ 2 \pi} \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega) \times \\
\hskip 1cm &\times \left[ \frac{1}{ \Omega_{ib}^s - \omega +i \eta } + \frac{1}{ \Omega_{ja}^{s} + \omega + i\eta } \right] \nonumber
\end{align}
with $\; \Omega_{ib}^s = \Oms - ( \varb - \vari )$ and $\; \Omega_{ja}^s = \Oms - ( \vara - \varj ).$
In the present study, we use the exact spectral representation of $W(\omega)$ at the RPA level:
\begin{align*}
W_{ij,ab}(\omega) &= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m] \times \\
& \times \Big( \frac{1}{ \omega-\Omega_m^{RPA} + i\eta } - \frac{1}{ \omega + \Omega_m^{RPA} - i\eta } \Big)
\end{align*}
($\Omega_m^{RPA} > 0 $) so that
\begin{align}
\widetilde{W}_{ij,ab}( \Oms ) &= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m] \times \\
& \times \left[ \frac{ 1 }{ \Omega_{ib}^{s} - \Omega_m^{RPA} + i\eta } + \frac{ 1}{ \Omega_{ja}^{s} - \Omega_m^{RPA} + i\eta }
\right] \nonumber
\end{align}
\textcolor{red}{Due to excitonic effects, the lowest BSE ${\Omega}_1$ excitation energy stands lower than the lowest $\Omega_m^{RPA}$ excitation energy, so that
e.g. $( \Omega_{ib}^{s} - \Omega_m^{RPA} )$ is strictly negative and cannot diverge. Further, $\Omega_{ib}^{s}$ and $\Omega_{ja}^{s}$ are necessarily negative for in-gap low lying BSE excitations, such that e.g.
$$
\left| \frac{ 1 }{ \Omega_{ib}^{s} - \Omega_m^{RPA} } \right|
< \frac{1}{ \Omega_m^{RPA} }
$$
This leads to reduced electron-hole screening, namely larger electron-hole stabilising binding energy, as compared to the standard adiabatic BSE, leading to smaller (blue-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{Theory for chemists}
%=================================
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
\begin{equation}
\label{eq:LR-dyn}
\begin{pmatrix}
\bA{}(\omega) & \bB{}(\omega) \\
-\bB{}(\omega) & -\bA{}(\omega) \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{m}{}(\omega) \\
\bY{m}{}(\omega) \\
\end{pmatrix}
=
\omega
\begin{pmatrix}
\bX{m}{}(\omega) \\
\bY{m}{}(\omega) \\
\end{pmatrix},
\end{equation}
where the dynamical matrices $\bA{}(\omega)$, $\bB{}(\omega)$, $\bX{}{}(\omega)$, and $\bY{}{}(\omega)$ 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.
The BSE matrix elements read
\begin{subequations}
\begin{align}
\label{eq:BSE-Adyn}
\A{ia,jb}{}(\omega) & = \delta_{ij} \delta_{ab} \eGW{ia} + 2 \ERI{ia}{jb} - \W{ij,ab}{}(\omega),
\\
\label{eq:BSE-Bdyn}
\B{ia,jb}{}(\omega) & = 2 \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,
\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 (linear) static response problem
\begin{equation}
\label{eq:LR-stat}
\begin{pmatrix}
\bA{\RPA} & \bB{\RPA} \\
-\bB{\RPA} & -\bA{\RPA} \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{m}{\RPA} \\
\bY{m}{\RPA} \\
\end{pmatrix}
=
\OmRPA{m}
\begin{pmatrix}
\bX{m}{\RPA} \\
\bY{m}{\RPA} \\
\end{pmatrix},
\end{equation}
with
\begin{subequations}
\begin{align}
\label{eq:LR_RPA-A}
\A{ia,jb}{\RPA} & = \delta_{ij} \delta_{ab} (\e{a} - \e{i}) + 2 \ERI{ia}{jb},
\\
\label{eq:LR_RPA-B}
\B{ia,jb}{\RPA} & = 2 \ERI{ia}{bj},
\end{align}
\end{subequations}
where the $\e{p}$'s are taken as the Hartree-Fock (HF) orbital energies in the case of $G_0W_0$ or as the $GW$ quasiparticle energies in the case of self-consistent scheme such as ev$GW$.
Now, let us decompose, using basis perturbation theory, the eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static part and a first-order dynamic part, such that
\begin{equation}
\label{eq:LR-dyn}
\begin{pmatrix}
\bA{}(\omega) & \bB{}(\omega) \\
-\bB{}(\omega) & -\bA{}(\omega) \\
\end{pmatrix}
=
\begin{pmatrix}
\bA{(0)} & \bB{(0)} \\
-\bB{(0)} & -\bA{(0)} \\
\end{pmatrix}
+
\begin{pmatrix}
\bA{(1)}(\omega) & \bB{(1)}(\omega) \\
-\bB{(1)}(\omega) & -\bA{(1)}(\omega) \\
\end{pmatrix}
\end{equation}
where
\begin{subequations}
\begin{align}
\label{eq:BSE-0}
\A{ia,jb}{(0)} & = \delta_{ij} \delta_{ab} \eGW{ia} + 2 \ERI{ia}{jb} - \W{ij,ab}{\text{stat}},
\\
\label{eq:BSE-0}
\B{ia,jb}{(0)} & = 2 \ERI{ia}{bj} - \W{ib,aj}{\text{stat}},
\end{align}
\end{subequations}
and
\begin{subequations}
\begin{align}
\label{eq:BSE-1}
\A{ia,jb}{(1)}(\omega) & = - \W{ij,ab}{}(\omega) + \W{ij,ab}{\text{stat}},
\\
\label{eq:BSE-1}
\B{ia,jb}{(1)}(\omega) & = - \W{ib,aj}{}(\omega) + \W{ib,aj}{\text{stat}},
\end{align}
\end{subequations}
The static version of the screened Coulomb potential reads
\begin{equation}
\label{eq:Wstat}
\W{ij,ab}{\text{stat}} = \ERI{ij}{ab} - 4 \sum_m^{\Nocc \Nvir} \frac{\sERI{ij}{m} \sERI{ab}{m}}{\OmRPA{m}{} - i \eta}.
\end{equation}
The $m$th BSE excitation energy and its corresponding eigenvector can then decomposed as
\begin{subequations}
\begin{gather}
\Om{m}{} = \Om{m}{(0)} + \Om{m}{(1)} + \ldots
\\
\begin{pmatrix}
\bX{m}{} \\
\bY{m}{} \\
\end{pmatrix}
=
\begin{pmatrix}
\bX{m}{(0)} \\
\bY{m}{(0)} \\
\end{pmatrix}
+
\begin{pmatrix}
\bX{m}{(1)} \\
\bY{m}{(1)} \\
\end{pmatrix}
+ \ldots
\end{gather}
\end{subequations}
Solving the zeroth-order static problem yields
\begin{equation}
\begin{pmatrix}
\bA{(0)} & \bB{(0)} \\
-\bB{(0)} & -\bA{(0)} \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{m}{(0)} \\
\bY{m}{(0)} \\
\end{pmatrix}
=
\Om{m}{(0)}
\begin{pmatrix}
\bX{m}{(0)} \\
\bY{m}{(0)} \\
\end{pmatrix},
\end{equation}
Thanks to first-order perturbation theory, the first-order correction to the $m$th excitation energy is
\begin{equation}
\Om{m}{(1)} =
\T{\begin{pmatrix}
\bX{m}{(0)} \\
\bY{m}{(0)} \\
\end{pmatrix}}
\cdot
\begin{pmatrix}
\bA{(1)}(\Om{m}{(0)}) & \bB{(1)}(\Om{m}{(0)}) \\
-\bB{(1)}(\Om{m}{(0)}) & -\bA{(1)}(\Om{m}{(0)}) \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{m}{(0)} \\
\bY{m}{(0)} \\
\end{pmatrix}.
\end{equation}
From a practical point of view, if one enforces the Tamm-Dancoff approximation (TDA), we obtain the very simple expression
\begin{equation}
\Om{m}{(1)} = \T{(\bX{m}{(0)})} \cdot \bA{(1)}(\Om{m}{(0)}) \cdot \bX{m}{(0)}.
\end{equation}
This correction can be renormalized by computing, at basically no extra cost, the renormalization factor
\begin{equation}
Z_{m} = \qty[ \T{(\bX{m}{(0)})} \cdot \left. \pdv{\bA{(1)}(\omega)}{\omega} \right|_{\omega = \Om{m}{(0)}} \cdot \bX{m}{(0)} ]^{-1}.
\end{equation}
which finally yields
\begin{equation}
\Om{m}{} \approx \Om{m}{(0)} + Z_{m} \Om{m}{(1)}.
\end{equation}
This is our final expression.
%%% 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}
Reference in New Issue View Git Blame Copy Permalink