clean up BSE theory

This commit is contained in:
Pierre-Francois Loos 2020-05-22 08:49:43 +02:00
parent a54ad63474
commit f4708431b1

View File

@ -80,14 +80,14 @@
\newcommand{\EcBSE}{E_\text{c}^\text{BSE}}
% 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{\e}[1]{\eps_{#1}}
\newcommand{\eHF}[1]{\eps^\text{HF}_{#1}}
\newcommand{\eKS}[1]{\eps^\text{KS}_{#1}}
\newcommand{\eQP}[1]{\eps^\text{QP}_{#1}}
\newcommand{\eGOWO}[1]{\eps^\text{\GOWO}_{#1}}
\newcommand{\eGW}[1]{\eps^{GW}_{#1}}
\newcommand{\eevGW}[1]{\eps^\text{\evGW}_{#1}}
\newcommand{\eGnWn}[2]{\eps^\text{\GnWn{#2}}_{#1}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
\newcommand{\tOm}[2]{\Tilde{\Omega}_{#1}^{#2}}
@ -180,14 +180,18 @@
\newcommand{\pis}{\pi^*}
\newcommand{\ra}{\rightarrow}
\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\vari{{\eps}_i}
\newcommand\vara{{\eps}_a}
\newcommand\varj{{\eps}_j}
\newcommand\varb{{\eps}_b}
\newcommand\varn{{\eps}_n}
\newcommand\varm{{\eps}_m}
\newcommand\Oms{{\Omega}_s}
\newcommand\hOms{\frac{{\Omega}_s}{2}}
\newcommand{\hpsi}{\Hat{\psi}}
\newcommand{\ha}{\Hat{a}}
\newcommand{\tchi}{\Tilde{\chi}}
\newcommand{\bx}{\mathbf{x}}
\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}
@ -227,11 +231,11 @@ In recent years, it has been shown to be a valuable tool for computational theor
Taking the optical gap (\ie, the lowest optical excitation energy) as an example, BSE builds on top of a $GW$ calculation by adding up excitonic effects $\EB$ to the $GW$ HOMO-LUMO gap
\begin{equation}
\Eg^{\GW} = \eps_{\LUMO}^{\GW} - \varepsilon_{\HOMO}^{\GW},
\Eg^{\GW} = \eps_{\LUMO}^{\GW} - \eps_{\HOMO}^{\GW},
\end{equation}
which is itself a corrected version of the Kohn-Sham (KS) gap
\begin{equation}
\Eg^{\KS} = \eps_{\LUMO}^{\KS} - \varepsilon_{\HOMO}^{\KS} \ll \Eg^{\GW} \approx \EgFun,
\Eg^{\KS} = \eps_{\LUMO}^{\KS} - \eps_{\HOMO}^{\KS} \ll \Eg^{\GW} \approx \EgFun,
\end{equation}
in order to approximate the optical gap
\begin{equation}
@ -272,100 +276,123 @@ Unless otherwise stated, atomic units are used.
\label{sec:theory}
%%%%%%%%%%%%%%%%%%%%%%%%
In this Section, we describe the theoretical foundations leading to the dynamical Bethe-Salpeter equation, following the seminal work by Strinati, \cite{Strinati_1988} presenting in a second step the perturbative implementation \cite{Rohlfing_2000,Ma_2009a,Ma_2009b} of the dynamical correction as compared to the standard adiabatic calculations. More details of the derivation are provided in ...
In this Section, following the seminal work by Strinati, \cite{Strinati_1988} we describe, first, the theoretical foundations leading to the dynamical Bethe-Salpeter equation.
We present, in a second step, the perturbative implementation \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b} of the dynamical correction as compared to the standard static approximation.
More details of the derivation are provided as {\SI}.
%================================
\subsection{General dynamical BSE theory}
%=================================
The resolution \cite{Strinati_1988} of the Bethe-Salpeter equation:
\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 system excited states, 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^+$, with e.g. $t_2^+ = t_2 + 0^+$ where $0^+$ is a small positive infinitesimal, 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 in $L(1,2; 1',2')$ and $L(6,2;5,2')$, simplifying further by ${\tilde \chi}_s(x_2,x_{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 ) = \int \frac{ d \omega }{ 2\pi } \; 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} \times (\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 }
\left[ \theta( \tau ) e^{i ( \vari + \hOms) \tau }
+ \theta( - \tau ) e^{i (\vara - \hOms \tau } \right]
\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}
where $\; \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
The resolution \cite{Strinati_1988} of the Bethe-Salpeter equation
\begin{multline}
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{multline}
with
\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}
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')
& = \mel{N}{T \hpsi(1) \hpsi(2) \hpsi^{\dagger}(2') \hpsi^{\dagger}(1')}{N},
\end{align}
where, \eg, $1 = (\bx_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 $\ket{N,s}$ of the $N$-electron system excited states, with $\ket{N} = \ket{N,0}$ the ground-state.
In the optical limit of instantaneous electron-hole creation and destruction, imposing $t_{2'} = t_2^+$ and $t_{1'} = t_1^+$, with, \eg, $t_2^+ = t_2 + 0^+$ where $0^+$ is a small positive infinitesimal, one gets
\begin{equation}
\begin{split}
iL(1,2;1',2')
& = \theta(+\tau_{12}) \sum_{s > 0} \chi_s(x_1,x_{1'}) \tchi_s(x_2,x_{2'}) e^{ - i \Oms \tau_{12} }
\\
& - \theta(-\tau_{12}) \sum_{s > 0} \chi_s(x_2,x_{2'}) \tchi_s(x_1,x_{1'}) e^{ + i \Oms \tau_{12} },
\end{split}
\end{equation}
with $\tau_{12} = t_1 - t_2$ and
\begin{align}
\chi_s(x_1,x_{1'}) & = \mel{N}{T \hpsi(x_1) \hpsi^{\dagger}(x_{1'})}{N,s}
\\
\tchi_s(x_2,x_{2'}) & = \mel{N,s}{T \hpsi(x_2) \hpsi^{\dagger}(x_{2'})}{N}
\end{align}
The $\Oms$'s are the neutral excitation energies of interest.
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(x_2,x_{2'})$ on both side of the BSE, 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 ionization 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 ) = \int \frac{d\omega}{2\pi} \; G\qty(\omega - \frac{\omega_1}{2} ) G\qty( {\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, \eg,
\begin{equation}
G(x_1,x_2 ; \omega) = \sum_n \frac{ \phi_n(x_1) \phi_n^*(x_2) } { \omega - \eps_n + i \eta \text{sgn} \times (\eps_n - \mu) }
\end {equation}
and projecting on $\phi_a^*(x_1) \phi_i(x_{1'})$, one obtains the $\omega_1= \Oms$ component
\begin{multline}
\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} }
\frac{ \phi_a^*(x_3) \phi_i(x_4) } { \Oms - ( \vara - \vari ) + i \eta }
\times \qty[ \theta( \tau ) e^{i ( \vari + \hOms) \tau } + \theta( - \tau ) e^{i (\vara - \hOms \tau } ]
\end{multline}
with $\tau = \tau_{34}$. Adopting now the $GW$ approximation for the exchange-correlation self-energy leads to a simplification of the BSE kernel:
\begin{equation}
\Xi(3,5;4,6) = v(3,6) \delta(34) \delta(56) - W(3^+,4) \delta(36) \delta(45)
\end{equation}
We further obtain the needed spectral representation of $\mel{N}{T \hpsi(3) \hpsi^{\dagger}(4)}{N,s}$ expanding the field operators over a complete orbital basis creation/destruction operators:
\begin{multline}
\mel{N}{T \hpsi(3) \hpsi^{\dagger}(4)}{N,s}
\\
= - \qty( e^{ -i \Omega_s t^{34} } ) \sum_{mn} \phi_m(x_3) \phi_n^*(x_4)
\mel{N}{\ha_n^{\dagger} \ha_m}{N,s}
\\
\times \qty[ \theta( \tau ) e^{- i ( \eps_m - \hOms ) \tau } + \theta( -\tau ) e^{ - i ( \eps_n + \hOms) \tau } ]
\end{multline}
with $\tau = \tau_{34}$ and where the $ \lbrace \eps_{n/m} \rbrace$ are proper addition/removal energies such that
\begin{equation}
e^{i H \tau} \ha_m^{\dagger} \ket{N} = e^{ i (E_0^N + \eps_m ) \tau } \ha _m^{\dagger} \ket{N}
\end{equation}
The $GW$ quasiparticle energies $\eps{n/m}^{GW}$ are good approximations to such removal/addition energies.
Selecting $(n,m)=(j,b)$ yields the largest components
$A_{jb}^{s} = \mel{N}{\ha_j^{\dagger} \ha_b}{N,s}$, while $(n,m)=(b,j)$ yields much weaker
$B_{jb}^{s} = \mel{N}{\ha_b^{\dagger} \ha_j}{N,s}$ 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 $\mel{N}{T \hpsi(5) \hpsi^{\dagger}(5)}{N,s}$ and $\mel{N}{T \hpsi(x_1) \hpsi^{\dagger}(x_{1'})}{N,s}$, and projecting onto $\phi_a^*(x_1) \phi_i(x_{1'})$, one obtains after a few tedious manipulations (see {\SI}) the dynamical Bethe-Salpeter equation (dBSE) :
\begin{equation}
\begin{split}
( \eps_a - \eps_i - \Omega_s ) A_{ia}^{s}
& + \sum_{jb} \qty[ v_{ai,bj} - \widetilde{W}_{ij,ab}(\Oms) ] A_{jb}^{s} \\
& + \sum_{jb} \qty[ v_{ai,jb} - \widetilde{W}_{ib,aj}(\Oms) ] B_{jb}^{s}
= 0
\end{split}
\end{equation}
with an effective dynamically screened Coulomb potential (see Pina eq. 24):
\begin{multline}
\widetilde{W}_{ij,ab}(\Oms)
= \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 } ]
\end{multline}
where $\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{multline}
W_{ij,ab}(\omega)
= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m]
\\
\times \qty( \frac{1}{ \omega-\Omega_m^{RPA} + i\eta } - \frac{1}{ \omega + \Omega_m^{RPA} - i\eta } )
\end{multline}
($\Omega_m^{RPA} > 0 $) so that
\begin{multline}
\widetilde{W}_{ij,ab}( \Oms )
= (ij|ab) + 2 \sum_m^{OV} [ij|m] [ab|m]
\\
\times \qty( \frac{1}{\Omega_{ib}^{s} - \Omega_m^{RPA} + i\eta} + \frac{1}{\Omega_{ja}^{s} - \Omega_m^{RPA} + i\eta} )
\end{multline}
\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 $\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.
$$
\left| \frac{ 1 }{ \Omega_{ib}^{s} - \Omega_m^{RPA} } \right|
< \frac{1}{ \Omega_m^{RPA} }
\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. }
@ -738,6 +765,7 @@ All the BSE calculations have been performed with our locally developed $GW$ sof
& 5.81 & 5.73 & 6.30 & 5.72 & 6.26
& 5.63 & 5.85 & 6.22 & 5.94 & 6.33 \\
\\
%T2: check state ordering in BSE calculation
\ce{C2H4} & $^1B_{3u}(\pi \ra 3s)$ & 11.49 & 7.64 & 7.62 & -0.03 & 1.004
& 7.35 & 7.34 & 7.42 & 7.29 & 7.35
& 6.63 & 6.88 & 6.94 & 6.93 & 7.57 \\