598 lines
32 KiB
TeX
598 lines
32 KiB
TeX
\documentclass[aip,jcp,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[utf8]{inputenc}
|
|
\usepackage[T1]{fontenc}
|
|
\usepackage{txfonts}
|
|
|
|
\usepackage[
|
|
colorlinks=true,
|
|
citecolor=blue,
|
|
breaklinks=true
|
|
]{hyperref}
|
|
\urlstyle{same}
|
|
|
|
\begin{document}
|
|
|
|
\title{Spin-Conserved and Spin-Flip Optical Excitations From the Bethe-Salpeter Equation Formalism}
|
|
|
|
\author{Enzo \surname{Monino}}
|
|
\affiliation{\LCPQ}
|
|
\author{Pierre-Fran\c{c}ois \surname{Loos}}
|
|
\email{loos@irsamc.ups-tlse.fr}
|
|
\affiliation{\LCPQ}
|
|
|
|
\begin{abstract}
|
|
\alert{Here comes the abstract.}
|
|
%\bigskip
|
|
%\begin{center}
|
|
% \boxed{\includegraphics[width=0.5\linewidth]{TOC}}
|
|
%\end{center}
|
|
%\bigskip
|
|
\end{abstract}
|
|
|
|
\maketitle
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Introduction}
|
|
\label{sec:intro}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\alert{Here comes the introduction.}
|
|
Unless otherwise stated, atomic units are used, and we assume real quantities throughout this manuscript.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Unrestricted $GW$ formalism}
|
|
\label{sec:UGW}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
Let us consider an electronic system consisting of $n = n_\up + n_\dw$ electrons (where $n_\up$ and $n_\dw$ are the number of spin-up and spin-down electrons, respectively) and $N$ one-electron basis functions.
|
|
The number of spin-up and spin-down occupied orbitals are $O_\up = n_\up$ and $O_\dw = n_\dw$, respectively, and, assuming no linear dependencies in the one-electron basis set, there is $V_\up = N - O_\up$ and $V_\dw = N - O_\dw$ spin-up and spin-down virtual (\ie, unoccupied) orbitals.
|
|
The number of spin-conserved (sc) single excitations is then $S^\spc = S_{\up\up}^\spc + S_{\dw\dw}^\spc = O_\up V_\up + O_\dw V_\dw$, while the number of spin-flip excitations is $S^\spf = S_{\up\dw}^\spf + S_{\dw\up}^\spf = O_\up V_\dw + O_\dw V_\up$.
|
|
Let us denote as $\MO{p_\sig}(\br)$ the $p$th (spin)orbital of spin $\sig$ (where $\sig =$ $\up$ or $\dw$) and $\e{p_\sig}{}$ its one-electron energy.
|
|
It is important to understand that, in a spin-conserved excitation the hole orbital $\MO{i_\sig}$ and particle orbital $\MO{a_\sig}$ have the same spin $\sig$.
|
|
In a spin-flip (sf) excitation, the hole and particle states, $\MO{i_\sig}$ and $\MO{a_\bsig}$, have opposite spins, $\sig$ and $\bsig$.
|
|
In the following, we assume real quantities throughout this manuscript, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, $p$, $q$, $r$, and $s$ indicate arbitrary orbitals, and $m$ labels single excitations.
|
|
Moreover, we consider systems with collinear spins and a spin-independent hamiltonian without contributions such as spin-orbit interaction.
|
|
|
|
%================================
|
|
\subsection{The dynamical screening}
|
|
%================================
|
|
The pillar of Green's function many-body perturbation theory is the (time-ordered) one-body Green's function, which has poles at the charged excitations (i.e., ionization potentials and electron affinities) of the system. \cite{ReiningBook}
|
|
The spin-$\sig$ component of the one-body Green's function reads \cite{ReiningBook,Bruneval_2016a}
|
|
\begin{equation}
|
|
\label{eq:G}
|
|
G^{\sig}(\br_1,\br_2;\omega)
|
|
= \sum_i \frac{\MO{i_\sig}(\br_1) \MO{i_\sig}(\br_2)}{\omega - \e{i_\sig}{} - i\eta}
|
|
+ \sum_a \frac{\MO{a_\sig}(\br_1) \MO{a_\sig}(\br_2)}{\omega - \e{a_\sig}{} + i\eta}
|
|
\end{equation}
|
|
where $\eta$ is a positive infinitesimal.
|
|
As readily seen in Eq.~\eqref{eq:G}, the Green's function can be evaluated at different levels of theory depending on the choice of orbitals and energies, $\MO{p_\sig}$ and $\e{p_\sig}{}$.
|
|
For example, $G_{\KS}^{\sig}$ is the independent-particle Green's function built with KS orbitals $\MO{p_\sig}^{\KS}(\br)$ and one-electron energies $\e{p_\sig}^{\KS}$.
|
|
Within self-consistent schemes, these quantities can be replaced by quasiparticle energies and orbitals evaluated within the $GW$ approximation (see below). \cite{Hedin_1965,Golze_2019}
|
|
|
|
Based on the spin-up and spin-down components of $G$ defined in Eq.~\eqref{eq:G}, one can easily compute the non-interacting polarizability (which is a sum over spins)
|
|
\begin{equation}
|
|
\label{eq:chi0}
|
|
\chi_0(\br_1,\br_2;\omega) = - \frac{i}{2\pi} \sum_\sig \int G^{\sig}(\br_1,\br_2;\omega+\omega') G^{\sig}(\br_1,\br_2;\omega') d\omega'
|
|
\end{equation}
|
|
and subsequently the dielectric function
|
|
\begin{equation}
|
|
\label{eq:eps}
|
|
\epsilon(\br_1,\br_2;\omega) = \delta(\br_1 - \br_2) - \int \frac{\chi_0(\br_1,\br_3;\omega) }{\abs{\br_2 - \br_3}} d\br_3
|
|
\end{equation}
|
|
where $\delta(\br)$ is the Dirac delta function.
|
|
Based on this latter ingredient, one can access the dynamically-screened Coulomb potential
|
|
\begin{equation}
|
|
\label{eq:W}
|
|
W(\br_1,\br_2;\omega) = \int \frac{\epsilon^{-1}(\br_1,\br_3;\omega) }{\abs{\br_2 - \br_3}} d\br_3
|
|
\end{equation}
|
|
which is naturally spin independent as the bare Coulomb interaction $\abs{\br_1 - \br_2}^{-1}$ does not depend on spin coordinates.
|
|
|
|
Within the $GW$ formalism, \cite{Hedin_1965,Onida_2002,Golze_2019} the dynamical screening is computed at the random-phase approximation (RPA) level by considering only the manifold of the spin-conserved neutral excitations.
|
|
In the orbital basis, the spectral representation of $W$ is
|
|
\begin{multline}
|
|
\label{eq:W_spectral}
|
|
W_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega) = \ERI{p_\sig q_\sig}{r_\sigp s_\sigp}
|
|
+ \sum_m \ERI{p_\sig q_\sig}{m}\ERI{r_\sigp s_\sigp}{m}
|
|
\\
|
|
\times \qty[ \frac{1}{\omega - \Om{m}{\spc,\RPA} + i \eta} - \frac{1}{\omega + \Om{m}{\spc,\RPA} - i \eta} ]
|
|
\end{multline}
|
|
where the bare two-electron integrals are \cite{Gill_1994}
|
|
\begin{equation}
|
|
\ERI{p_\sig q_\tau}{r_\sigp s_\taup} = \int \frac{\MO{p_\sig}(\br_1) \MO{q_\tau}(\br_1) \MO{r_\sigp}(\br_2) \MO{s_\taup}(\br_2)}{\abs{\br_1 - \br_2}} d\br_1 d\br_2
|
|
\end{equation}
|
|
and the screened two-electron integrals (or spectral weights) are explicitly given by
|
|
\begin{equation}
|
|
\label{eq:sERI}
|
|
\ERI{p_\sig q_\sig}{m} = \sum_{ia\sigp} \ERI{p_\sig q_\sig}{r_\sigp s_\sigp} (\bX{m}{\spc,\RPA}+\bY{m}{\spc,\RPA})_{i_\sigp a_\sigp}
|
|
\end{equation}
|
|
In Eqs.~\eqref{eq:W_spectral} and \eqref{eq:sERI}, the spin-conserved RPA neutral excitations $\Om{m}{\spc,\RPA}$ and their corresponding eigenvectors, $\bX{m}{\spc,\RPA}$ and $\bY{m}{\spc,\RPA}$, are obtained by solving a linear response system of the form
|
|
\begin{equation}
|
|
\label{eq:LR-RPA}
|
|
\begin{pmatrix}
|
|
\bA{}{} & \bB{}{} \\
|
|
-\bB{}{} & -\bA{}{} \\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX{m}{} \\
|
|
\bY{m}{} \\
|
|
\end{pmatrix}
|
|
=
|
|
\Om{m}{}
|
|
\begin{pmatrix}
|
|
\bX{m}{} \\
|
|
\bY{m}{} \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
where the expressions of the matrix elements of $\bA{}{}$ and $\bB{}{}$ are specific of the method and of the spin manifold.
|
|
The spin structure of these matrices, though, is general
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR-RPA-AB-sc}
|
|
\bA{}{\spc} & = \begin{pmatrix}
|
|
\bA{}{\upup,\upup} & \bA{}{\upup,\dwdw} \\
|
|
\bA{}{\dwdw,\upup} & \bA{}{\dwdw,\dwdw} \\
|
|
\end{pmatrix}
|
|
&
|
|
\bB{}{\spc} & = \begin{pmatrix}
|
|
\bB{}{\upup,\upup} & \bB{}{\upup,\dwdw} \\
|
|
\bB{}{\dwdw,\upup} & \bB{}{\dwdw,\dwdw} \\
|
|
\end{pmatrix}
|
|
\\
|
|
\label{eq:LR-RPA-AB-sf}
|
|
\bA{}{\spf} & = \begin{pmatrix}
|
|
\bA{}{\updw,\updw} & \bO \\
|
|
\bO & \bA{}{\dwup,\dwup} \\
|
|
\end{pmatrix}
|
|
&
|
|
\bB{}{\spf} & = \begin{pmatrix}
|
|
\bO & \bB{}{\updw,\dwup} \\
|
|
\bB{}{\dwup,\updw} & \bO \\
|
|
\end{pmatrix}
|
|
\end{align}
|
|
\end{subequations}
|
|
In the absence of instabilities, the linear eigenvalue problem \eqref{eq:LR-RPA} has particle-hole symmetry which means that the eigenvalues are obtained by pairs $\pm \Om{m}{}$.
|
|
In such a case, $(\bA{}{}-\bB{}{})^{1/2}$ is positive definite, and Eq.~\eqref{eq:LR-RPA} can be recast as a Hermitian problem of half its original dimension
|
|
\begin{equation}
|
|
\label{eq:small-LR}
|
|
(\bA{}{} - \bB{}{})^{1/2} \cdot (\bA{}{} + \bB{}{}) \cdot (\bA{}{} - \bB{}{})^{1/2} \cdot \bZ{}{} = \bOm{2} \cdot \bZ{}{}
|
|
\end{equation}
|
|
where the excitation amplitudes are
|
|
\begin{equation}
|
|
\bX{}{} + \bY{}{} = \bOm{-1/2} \cdot (\bA{}{} - \bB{}{})^{1/2} \cdot \bZ{}{}
|
|
\end{equation}
|
|
Within the Tamm-Dancoff approximation (TDA), the coupling terms between the resonant and anti-resonant parts, $\bA{}{}$ and $-\bA{}{}$, are neglected, which consists in setting $\bB{}{} = \bO$.
|
|
In such a case, Eq.~\eqref{eq:LR-RPA} reduces to a straightforward Hermitian problem of the form:
|
|
\begin{equation}
|
|
\bA{}{} \cdot \bX{m}{} = \Om{m}{} \bX{m}{}
|
|
\end{equation}
|
|
Note that, for spin-flip excitations, it is quite common to enforce the TDA especially when one considers a triplet reference as the first ``excited-state'' is usually the ground state of the closed-shell system (hence, corresponding to a negative excitation energy).
|
|
|
|
At the RPA level, the matrix elements of $\bA{}{}$ and $\bB{}{}$ are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-A}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} & = \delta_{ij} \delta_{ab} \delta_{\sig \sigp} \delta_{\tau \taup} (\e{a_\tau} - \e{i_\sig}) + \ERI{i_\sig a_\tau}{b_\sigp j_\taup}
|
|
\\
|
|
\label{eq:LR_RPA-B}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} & = \ERI{i_\sig a_\tau}{j_\sigp b_\taup}
|
|
\end{align}
|
|
\end{subequations}
|
|
from which we obtain the following expressions
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-Asc}
|
|
\A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} & = \delta_{ij} \delta_{ab} \delta_{\sig \sigp} (\e{a_\sig} - \e{i_\sig}) + \ERI{i_\sig a_\sig}{b_\sigp j_\sigp}
|
|
\\
|
|
\label{eq:LR_RPA-Bsc}
|
|
\B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} & = \ERI{i_\sig a_\sig}{j_\sigp b_\sigp}
|
|
\end{align}
|
|
\end{subequations}
|
|
for the spin-conserved excitations and
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-Asf}
|
|
\A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\RPA} & = \delta_{ij} \delta_{ab} (\e{a_\bsig} - \e{i_\sig})
|
|
\\
|
|
\label{eq:LR_RPA-Bsf}
|
|
\B{i_\sig a_\bsig,j_\bsig b_\sig}{\spf,\RPA} & = 0
|
|
\end{align}
|
|
\end{subequations}
|
|
for the spin-flip excitations.
|
|
|
|
%================================
|
|
\subsection{The $GW$ self-energy}
|
|
%================================
|
|
Within the acclaimed $GW$ approximation, \cite{Hedin_1965,Golze_2019} the exchange-correlation (xc) part of the self-energy
|
|
\begin{equation}
|
|
\label{eq:Sig}
|
|
\begin{split}
|
|
\Sig{}{\xc,\sig}(\br_1,\br_2;\omega)
|
|
& = \Sig{}{\x,\sig}(\br_1,\br_2) + \Sig{}{\co,\sig}(\br_1,\br_2;\omega)
|
|
\\
|
|
& = \frac{i}{2\pi} \int G^{\sig}(\br_1,\br_2;\omega+\omega') W(\br_1,\br_2;\omega') e^{i \eta \omega'} d\omega'
|
|
\end{split}
|
|
\end{equation}
|
|
is, like the one-body Green's function, spin-diagonal, and its spectral representation reads
|
|
\begin{subequations}
|
|
\begin{gather}
|
|
\Sig{p_\sig q_\sig}{\x}
|
|
= - \sum_{i} \ERI{p_\sig i_\sig}{i_\sig q_\sig}
|
|
\\
|
|
\begin{split}
|
|
\Sig{p_\sig q_\sig}{\co}(\omega)
|
|
& = \sum_{im} \frac{\ERI{p_\sig i_\sig}{m} \ERI{q_\sig i_\sig}{m}}{\omega - \e{i_\sig} + \Om{m}{\spc,\RPA} - i \eta}
|
|
\\
|
|
& + \sum_{am} \frac{\ERI{p_\sig a_\sig}{m} \ERI{q_\sig a_\sig}{m}}{\omega - \e{a_\sig} - \Om{m}{\spc,\RPA} + i \eta}
|
|
\end{split}
|
|
\end{gather}
|
|
\end{subequations}
|
|
where the self-energy has been split in its exchange (x) and correlation (c) contributions.
|
|
The Dyson equation linking the Green's function and the self-energy holds separately for each spin component
|
|
\begin{equation}
|
|
\label{eq:Dyson_G}
|
|
\begin{split}
|
|
\qty[ G^{\sig}(\br_1,\br_2;\omega) ]^{-1}
|
|
& = \qty[ G_{\KS}^{\sig}(\br_1,\br_2;\omega) ]^{-1}
|
|
\\
|
|
& + \Sig{}{\xc,\sig}(\br_1,\br_2;\omega) - v^{\xc}(\br_1) \delta(\br_1 - \br_2)
|
|
\end{split}
|
|
\end{equation}
|
|
where $v^{\xc}(\br)$ is the KS (local) exchange-correlation potential.
|
|
The target quantities here are the quasiparticle energies $\eGW{p_\sig}$, \ie, the poles of $G$ [see Eq.~\eqref{eq:G}], which correspond to well-defined addition/removal energies (unlike the KS orbital energies).
|
|
Because the exchange-correlation part of the self-energy is, itself, constructed with the Green's function [see Eq.~\eqref{eq:Sig}], the present process is, by nature, self-consistent.
|
|
The same comment applies to the dynamically-screened Coulomb potential $W$ entering the definition of $\Sig{}{\xc}$ [see Eq.~\eqref{eq:Sig}] which is also constructed from $G$ [see Eqs.~\eqref{eq:chi0}, \eqref{eq:eps}, and \eqref{eq:W}].
|
|
|
|
%================================
|
|
\subsection{Level of self-consistency}
|
|
%================================
|
|
This is where $GW$ schemes differ.
|
|
In its simplest perturbative (\ie, one-shot) version, known as {\GOWO}, \cite{Strinati_1980,Hybertsen_1985a,Hybertsen_1986,Godby_1988,Linden_1988,Northrup_1991,Blase_1994,Rohlfing_1995,Shishkin_2007} a single iteration is performed, and the quasiparticle energies $\eGW{p_\sig}$ are obtained by solving the frequency-dependent quasiparticle equation
|
|
\begin{equation}
|
|
\label{eq:QP-eq}
|
|
\omega = \e{p_\sig}{} + \Sig{p_\sig}{\xc}(\omega) - V_{p_\sig}^{\xc}
|
|
\end{equation}
|
|
where $\Sig{p_\sig}{\xc}(\omega) \equiv \Sig{p_\sig p_\sig}{\xc}(\omega)$ and its offspring quantities have been constructed at the KS level, and
|
|
\begin{equation}
|
|
V_{p_\sigma}^{\xc} = \int \MO{p_\sig}(\br) v^{\xc}(\br) \MO{p_\sig}(\br) d\br
|
|
\end{equation}
|
|
Because, from a practical point of view, one is usually interested by the so-called quasiparticle solution (or peak), the quasiparticle equation \eqref{eq:QP-eq} is often linearized around $\omega = \e{p_\sig}^{\KS}$, yielding
|
|
\begin{equation}
|
|
\eGW{p_\sig}
|
|
= \e{p_\sig}^{\KS} + Z_{p_\sig} [\Sig{p_\sig}{\xc}(\e{p_\sig}^{\KS}) - V_{p_\sig}^{\xc} ]
|
|
\end{equation}
|
|
where
|
|
\begin{equation}
|
|
\label{eq:Z_GW}
|
|
Z_{p_\sig} = \qty[ 1 - \left. \pdv{\Sig{p_\sig}{\xc}(\omega)}{\omega} \right|_{\omega = \e{p_\sig}^{\KS}} ]^{-1}
|
|
\end{equation}
|
|
is a renormalization factor (with $0 \le Z_{p_\sig} \le 1$) which also represents the spectral weight of the quasiparticle solution.
|
|
In addition to the principal quasiparticle peak which, in a well-behaved case, contains most of the spectral weight, the frequency-dependent quasiparticle equation \eqref{eq:QP-eq} generates a finite number of satellite resonances with smaller weights.
|
|
|
|
Within the ``eigenvalue'' self-consistent $GW$ scheme (known as ev$GW$), \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016,Kaplan_2016,Gui_2018} several iterations are performed during which only the one-electron energies entering the definition of the Green's function [see Eq.~\eqref{eq:G}] are updated by the quasiparticle energies obtained at the previous iteration (the corresponding orbitals remain evaluated at the KS level).
|
|
|
|
Finally, within the quasiparticle self-consistent $GW$ (qs$GW$) scheme, both the one-electron energies and the orbitals are updated until convergence is reached.
|
|
These are obtained via the diagonalization of an effective Fock matrix which includes explicitly a frequency-independent and hermitian self-energy defined as
|
|
\begin{equation}
|
|
\Tilde{\Sigma}_{p_\sig q_\sig}^{\xc} = \frac{1}{2} \qty[ \Sig{p_\sig q_\sig}{\xc}(\e{p_\sig}{}) + \Sig{q_\sig p_\sig}{\xc}(\e{p_\sig}{}) ]
|
|
\end{equation}
|
|
|
|
%================================
|
|
\section{Unrestricted Bethe-Salpeter equation formalism}
|
|
%================================
|
|
Like its TD-DFT cousin, the BSE formalism \cite{Salpeter_1951,Strinati_1988,Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999} deals with the calculation of (neutral) optical excitations as measured by absorption spectroscopy. \cite{Rohlfing_1999a,Horst_1999,Puschnig_2002,Tiago_2003,Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018}
|
|
Using the BSE formalism, one can access the spin-conserved and spin-flip excitations.
|
|
In a nutshell, BSE builds on top of a $GW$ calculation by adding up excitonic effects (\ie, the electron-hole binding energy) to the $GW$ fundamental gap which is itself a corrected version of the KS gap.
|
|
The purpose of the underlying $GW$ calculation is to provide quasiparticle energies and a dynamically-screened Coulomb potential that are used to build the BSE Hamiltonian from which the vertical excitations of the system are extracted.
|
|
|
|
%================================
|
|
\subsection{Static approximation}
|
|
%================================
|
|
|
|
|
|
The Dyson equation that links the generalized four-point susceptibility $L^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)$ and the BSE kernel $\Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6;\omega)$ is
|
|
\begin{multline}
|
|
L^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
|
|
= L_{0}^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
|
|
\\
|
|
+ \int L_{0}^{\sig\sigp}(\br_1,\br_4;\br_1',\br_3;\omega)
|
|
\Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6;\omega)
|
|
\\
|
|
\times L^{\sig\sigp}(\br_6,\br_2;\br_5,\br_2';\omega)
|
|
d\br_3 d\br_4 d\br_5 d\br_6
|
|
\end{multline}
|
|
where
|
|
\begin{multline}
|
|
L_{0}^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
|
|
\\
|
|
= \frac{1}{2\pi} \int G^{\sig}(\br_1,\br_2';\omega+\omega') G^{\sig}(\br_1',\br_2;\omega') d\omega'
|
|
\end{multline}
|
|
is the non-interacting analog of the two-particle correlation function $L$.
|
|
|
|
Within the $GW$ approximation, the BSE kernel is
|
|
\begin{multline}
|
|
i \Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6;\omega)
|
|
= \frac{\delta(\br_3 - \br_4) \delta(\br_5 - \br_6) }{\abs{\br_3-\br_6}}
|
|
\\
|
|
- \delta_{\sig\sigp} W(\br_3,\br_4;\omega) \delta(\br_3 - \br_6) \delta(\br_4 - \br_6)
|
|
\end{multline}
|
|
where, as usual, we have not considered the higher-order terms in $W$ by neglecting the derivative $\partial W/\partial G$. \cite{Hanke_1980, Strinati_1982, Strinati_1984, Strinati_1988}
|
|
Within the static approximation which consists in neglecting the frequency dependence of the dynamically-screened Coulomb potential, the spin-conserved and spin-flip BSE optical excitations are obtained by solving the usual Casida-like linear response (eigen)problem:
|
|
\begin{equation}
|
|
\label{eq:LR-BSE}
|
|
\begin{pmatrix}
|
|
\bA{}{\BSE} & \bB{}{\BSE} \\
|
|
-\bB{}{\BSE} & -\bA{}{\BSE} \\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX{m}{\BSE} \\
|
|
\bY{m}{\BSE} \\
|
|
\end{pmatrix}
|
|
=
|
|
\Om{m}{\BSE}
|
|
\begin{pmatrix}
|
|
\bX{m}{\BSE} \\
|
|
\bY{m}{\BSE} \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
Defining the elements of the static screening as $W^{\stat}_{p_\sig q_\sig,r_\sigp s_\sigp} = W_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega = 0)$, the general expressions of the BSE matrix elements are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_BSE-A}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{\BSE} & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\taup a_\tau}
|
|
\\
|
|
\label{eq:LR_BSE-B}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{\BSE} & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\taup,j_\sigp a_\tau}
|
|
\end{align}
|
|
\end{subequations}
|
|
from which we obtain the following expressions for the spin-conserved and spin-flip BSE excitations:
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_BSE-Asc}
|
|
\A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\BSE} & = \A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\sigp a_\sig}
|
|
\\
|
|
\label{eq:LR_BSE-Bsc}
|
|
\B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\BSE} & = \B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\sigp,j_\sigp a_\sig}
|
|
\\
|
|
\label{eq:LR_BSE-Asf}
|
|
\A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\BSE} & = \A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\RPA} - W^{\stat}_{i_\sig j_\sig,b_\bsig a_\bsig}
|
|
\\
|
|
\label{eq:LR_BSE-Bsf}
|
|
\B{i_\sig a_\bsig,j_\bsig b_\sig}{\spf,\BSE} & = - W^{\stat}_{i_\sig b_\sig,j_\bsig a_\bsig}
|
|
\end{align}
|
|
\end{subequations}
|
|
|
|
At this stage, it is of particular interest to discuss the form of the spin-flip matrix elements defined in Eqs.~\eqref{eq:LR_BSE-Asf} and \eqref{eq:LR_BSE-Bsf}.
|
|
As readily seen from Eq.~\eqref{eq:LR_RPA-Asf}, at the RPA level, the spin-flip excitations are given by the difference of one-electron energies, hence missing out on key exchange and correlation effects.
|
|
This is also the case at the TD-DFT level when one relies on (semi-)local functionals.
|
|
This explains why most of spin-flip TD-DFT calculations are performed with hybrid functionals containing a substantial amount of Hartree-Fock exchange as only the exact exchange integral of the form $\ERI{i_\sig j_\sig}{b_\bsig a_\bsig}$ survive spin-symmetry requirements.
|
|
At the BSE level, these matrix elements are, of course, also present thanks to the contribution of $W^{\stat}_{i_\sig j_\sig,b_\bsig a_\bsig}$ as evidenced in Eq.~\eqref{eq:W_spectral} but it also includes correlation effects.
|
|
|
|
%================================
|
|
\subsection{Dynamical correction}
|
|
%================================
|
|
In order to go beyond the ubiquitous static approximation of BSE \cite{Strinati_1988,Rohlfing_2000,Sottile_2003,Myohanen_2008,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Huix-Rotllant_2011,Sakkinen_2012,Zhang_2013,Rebolini_2016,Olevano_2019,Lettmann_2019} (which is somehow similar to the adiabatic approximation of TD-DFT \cite{Casida_2005,Huix-Rotllant_2011,Casida_2016,Maitra_2004,Cave_2004,Elliott_2011,Maitra_2012}), we have recently implemented, following Strinati's seminal work \cite{Strinati_1982,Strinati_1984,Strinati_1988} (see also the work of Romaniello \textit{et al.} \cite{Romaniello_2009b} and Sangalli \textit{et al.} \cite{Sangalli_2011}), a renormalized first-order perturbative correction in order to take into consideration the dynamical nature of the screened Coulomb potential $W$. \cite{Loos_2020e,Authier_2020}
|
|
This dynamical correction to the static BSE kernel (dubbed as dBSE in the following) does permit to recover additional relaxation effects coming from higher excitations.
|
|
|
|
Our implementation follows closely the work of Rohlfing and co-workers \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b} in which they computed the dynamical correction in the TDA and plasmon-pole approximation.
|
|
However, our scheme goes beyond the plasmon-pole approximation as the spectral representation of the dynamically-screened Coulomb potential is computed exactly at the RPA level consistently with the underlying $GW$ calculation:
|
|
\begin{multline}
|
|
\widetilde{W}_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega) = \ERI{p_\sig q_\sig}{r_\sigp s_\sigp}
|
|
+ \sum_m \ERI{p_\sig q_\sig}{m}\ERI{r_\sigp s_\sigp}{m}
|
|
\\
|
|
\times \Bigg[ \frac{1}{\omega - (\e{s_\sigp}{} - \e{q_\sig}{}) - \Om{m}{\spc,\RPA} + i \eta}
|
|
\\
|
|
+ \frac{1}{\omega - (\e{r_\sigp}{} - \e{p_\sig}{}) - \Om{m}{\spc,\RPA} + i \eta} \Bigg]
|
|
\end{multline}
|
|
The dBSE non-linear response problem is
|
|
\begin{multline}
|
|
\label{eq:LR-dyn}
|
|
\begin{pmatrix}
|
|
\bA{}{\dBSE}(\Om{m}{\dBSE}) & \bB{}{\dBSE}(\Om{m}{\dBSE})
|
|
\\
|
|
-\bB{}{\dBSE}(-\Om{m}{\dBSE}) & -\bA{}{\dBSE}(-\Om{m}{\dBSE})
|
|
\\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX{m}{\dBSE} \\
|
|
\bY{m}{\dBSE} \\
|
|
\end{pmatrix}
|
|
\\
|
|
=
|
|
\Om{m}{\dBSE}
|
|
\begin{pmatrix}
|
|
\bX{m}{\dBSE} \\
|
|
\bY{m}{\dBSE} \\
|
|
\end{pmatrix}
|
|
\end{multline}
|
|
where the dynamical matrices are generally defined as
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_dBSE-A}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{\dBSE}(\omega) & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} \widetilde{W}_{i_\sig j_\sigp,b_\taup a_\tau}(\omega)
|
|
\\
|
|
\label{eq:LR_dBSE-B}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{\dBSE}(\omega) & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} \widetilde{W}_{i_\sig b_\taup,j_\sigp a_\tau}(\omega)
|
|
\end{align}
|
|
\end{subequations}
|
|
from which one can easily obtained the matrix elements for the spin-conserved and spin-flip manifolds similarly to Eqs.~\eqref{eq:LR_BSE-Asc}, \eqref{eq:LR_BSE-Bsc}, \eqref{eq:LR_BSE-Asf}, and \eqref{eq:LR_BSE-Bsf}.
|
|
Following Rayleigh-Schr\"odinger perturbation theory, we then decompose the non-linear eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static (\ie, linear) reference and a first-order dynamic (\ie, non-linear) perturbation such that
|
|
\begin{multline}
|
|
\label{eq:LR-PT}
|
|
\begin{pmatrix}
|
|
\bA{}{\dBSE}(\omega) & \bB{}{\dBSE}(\omega) \\
|
|
-\bB{}{\dBSE}(-\omega) & -\bA{}{\dBSE}(-\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{multline}
|
|
with
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:BSE-A0}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{(0)} & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\BSE}
|
|
\\
|
|
\label{eq:BSE-B0}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{(0)} & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\BSE}
|
|
\end{align}
|
|
\end{subequations}
|
|
and
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:BSE-A1}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{(1)}(\omega) & = - \delta_{\sig \sigp} \widetilde{W}_{i_\sig j_\sigp,b_\taup a_\tau}(\omega) + \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\taup a_\tau}
|
|
\\
|
|
\label{eq:BSE-B1}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{(1)}(\omega) & = - \delta_{\sig \sigp} \widetilde{W}_{i_\sig b_\taup,j_\sigp a_\tau}(\omega) + \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\taup,j_\sigp a_\tau}
|
|
\end{align}
|
|
\end{subequations}
|
|
The dBSE excitation energies are then obtained via
|
|
\begin{equation}
|
|
\Om{m}{\dBSE} = \Om{m}{\BSE} + \zeta_{m} \Om{m}{(1)}
|
|
\end{equation}
|
|
where $\Om{m}{\BSE} \equiv \Om{m}{(0)}$ are the static (zeroth-order) BSE excitation energies obtained by solving Eq.~\eqref{eq:LR-BSE}, and
|
|
\begin{equation}
|
|
\label{eq:Om1-TDA}
|
|
\Om{m}{(1)} = \T{(\bX{m}{\BSE})} \cdot \bA{}{(1)}(\Om{m}{\BSE}) \cdot \bX{m}{\BSE}
|
|
\end{equation}
|
|
are first-order corrections (with $\bX{m}{\BSE} \equiv \bX{m}{(0)}$) obtained within the dynamical TDA (dTDA) with the renormalization factor
|
|
\begin{equation}
|
|
\label{eq:Z}
|
|
\zeta_{m} = \qty[ 1 - \T{(\bX{m}{\BSE})} \cdot \left. \pdv{\bA{}{(1)}(\omega)}{\omega} \right|_{\omega = \Om{m}{\BSE}} \cdot \bX{m}{\BSE} ]^{-1}
|
|
\end{equation}
|
|
which, unlike the $GW$ case [see Eq.~\eqref{eq:Z_GW}], is not restricted to be between $0$ and $1$.
|
|
In most cases, the value of $\zeta_{m}$ is close to unity which indicates that the perturbative expansion behaves nicely.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\subsection{Oscillator strengths}
|
|
\label{sec:os}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
For the spin-conserved transitions, the $x$ component of the transition dipole moment is
|
|
\begin{equation}
|
|
\mu_{x,m}^{\spc} = \sum_{ia\sig} (i_\sig|x|a_\sig)(\bX{m}{\spc}+\bY{m}{\spc})_{i_\sig a_\sig}
|
|
\end{equation}
|
|
where
|
|
\begin{equation}
|
|
(p_\sig|x|q_\sigp) = \int \MO{p_\sig}(\br) \, x \, \MO{q_\sigp}(\br) d\br
|
|
\end{equation}
|
|
are one-electron integrals in the orbital basis.
|
|
The total oscillator strength is given by
|
|
\begin{equation}
|
|
f_{m}^{\spc} = \frac{2}{3} \Om{m}{\spc} \qty[ \qty(\mu_{x,m}^{\spc})^2 + \qty(\mu_{x,m}^{\spc})^2 + \qty(\mu_{x,m}^{\spc})^2 ]
|
|
\end{equation}
|
|
For spin-flip transitions, we have $f_{m}^{\spf} = 0$ as the transition matrix elements $(i_\sig|x|a_\bsig)$ vanish via integration over the spin coordinate.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\subsection{Spin contamination}
|
|
\label{sec:spin}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
One of the key issues of linear response formalism based on unrestricted references is spin contamination.
|
|
As nicely explained in Ref.~\onlinecite{Casanova_2020}, there are two sources of spin contamination: i) spin contamination of the reference configuration for which, for example, $\expval{\hS^2} > 2$ for a high-spin triplets, and ii) spin-contamination of the excited states due to spin incompleteness of the configuration interaction expansion.
|
|
The latter issue is an important source of spin contamination in the present context as BSE is limited to single excitations with respect to the reference configuration.
|
|
Specific schemes have been developed to palliate these shortcomings and we refer the interested reader to Ref.~\onlinecite{Casanova_2020} for a detailed discussion on this matter.
|
|
|
|
In order to monitor closely how contaminated are these states, we compute
|
|
\begin{equation}
|
|
\expval{\hS^2}_m = \expval{\hS^2}_0 + \Delta \expval{\hS^2}_m
|
|
\end{equation}
|
|
where
|
|
\begin{equation}
|
|
\expval{\hS^2}_{0}
|
|
= \frac{n_{\up} - n_{\dw}}{2} \qty( \frac{n_{\up} - n_{\dw}}{2} + 1 )
|
|
+ n_{\dw} - \sum_p (p_{\up}|p_{\dw})^2
|
|
\end{equation}
|
|
is the expectation value of $\hS^2$ for the reference configuration, the first term correspoding to the exact value of $\expval{\hS^2}$, and
|
|
\begin{equation}
|
|
\label{eq:OV}
|
|
(p_\sig|q_\sigp) = \int \MO{p_\sig}(\br) \MO{q_\sigp}(\br) d\br
|
|
\end{equation}
|
|
are overlap integrals between spin-up and spin-down orbitals.
|
|
|
|
For a given single excitation $m$, the explicit expressions of $\Delta \expval{\hS^2}_m^{\spc}$ and $\Delta \expval{\hS^2}_m^{\spf}$ can be found in the Appendix of Ref.~\onlinecite{Li_2011a} for spin-conserved and spin-flip excitations, and are functions of the $\bX{m}{}$ and $\bY{m}{}$ vectors and the orbital overlaps defined in Eq.~\eqref{eq:OV}.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Computational details}
|
|
\label{sec:compdet}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
For the systems under investigation here, we consider either an open-shell doublet or triplet reference state.
|
|
We then adopt the unrestricted formalism throughout this work.
|
|
The $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies are done using a (unrestricted) UHF starting point.
|
|
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a,Hybertsen_1986,vanSetten_2013} quasiparticle energies are employed as starting points to compute the BSE neutral excitations.
|
|
These quasiparticle energies are obtained by linearizing the frequency-dependent quasiparticle equation, and the entire set of orbitals is corrected.
|
|
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018b,Veril_2018,Loos_2020,Loos_2020e}.
|
|
%Note that, for the present (small) molecular systems, {\GOWO}@UHF and ev$GW$@UHF yield similar quasiparticle energies and fundamental gap.
|
|
%Moreover, {\GOWO} allows to avoid rather laborious iterations as well as the significant additional computational effort of ev$GW$.
|
|
%In the present study, the zeroth-order Hamiltonian [see Eq.~\eqref{eq:LR-PT}] is always the ``full'' BSE static Hamiltonian, \ie, without TDA.
|
|
The dynamical correction is computed in the TDA throughout.
|
|
As one-electron basis sets, we employ the Dunning families cc-pVXZ and aug-cc-pVXZ (X = D, T, and Q) defined with cartesian Gaussian functions.
|
|
Finally, the infinitesimal $\eta$ is set to $100$ meV for all calculations.
|
|
%It is important to mention that the small molecular systems considered here are particularly challenging for the BSE formalism, \cite{Hirose_2015,Loos_2018b} which is known to work best for larger systems where the amount of screening is more important. \cite{Jacquemin_2017b,Rangel_2017}
|
|
|
|
%For comparison purposes, we employ the theoretical best estimates (TBEs) and geometries of Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} from which CIS(D), \cite{Head-Gordon_1994,Head-Gordon_1995} ADC(2), \cite{Trofimov_1997,Dreuw_2015} CC2, \cite{Christiansen_1995a} CCSD, \cite{Purvis_1982} and CC3 \cite{Christiansen_1995b} excitation energies are also extracted.
|
|
%Various statistical quantities are reported in the following: the mean signed error (MSE), mean absolute error (MAE), root-mean-square error (RMSE), and the maximum positive [Max($+$)] and maximum negative [Max($-$)] errors.
|
|
All the static and dynamic BSE calculations have been performed with the software \texttt{QuAcK}, \cite{QuAcK} freely available on \texttt{github}.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Results}
|
|
\label{sec:res}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%% TABLE I %%%
|
|
\begin{squeezetable}
|
|
\begin{table*}
|
|
\caption{
|
|
Spin-flip excitations of \ce{Be} obtained for various methods with the 6-31G basis.
|
|
The $GW$ calculations are performed with a HF starting point.
|
|
\label{tab:Be}}
|
|
\begin{ruledtabular}
|
|
\begin{tabular}{lcccccccccccc}
|
|
State & TD-BLYP & TD-B3LYP & TD-BHHLYP & CIS
|
|
& BSE@{\GOWO} & BSE@{\evGW} & dBSE@{\GOWO} & dBSE@{\evGW}
|
|
& ADC(2) & ADC(2)-x & ADC(3) & FCI \\
|
|
\hline
|
|
$^3P(1s^22s2p)$ & 3.210 & 3.332 & 2.874 & 2.111 & 2.399& 2.407& 2.363& 2.369 & 2.433 & 2.866 & 2.863 & 2.862 \\
|
|
$^1P(1s^22s2p)$ & 3.210 & 4.275& 4.922& 6.036 & 6.191& 6.199 & 6.263 & 6.273 & 6.255 &6.581 & 6.579 & 6.577 \\
|
|
$^3P(1s^22p^2)$ & 6.691 & 6.864& 7.112 & 7.480 & 7.792 & 7.788 & 7.824 & 7.820 & 7.745 & 7.664 & 7.658 & 7.669 \\
|
|
$^1P(1s^22p^2)$ & 7.598 & 7.762& 8.188 & 8.945 &9.373 & 9.388 & 9.424 & 9.441 & 9.047 & 8.612 & 8.618 & 8.624 \\
|
|
\end{tabular}
|
|
\end{ruledtabular}
|
|
\end{table*}
|
|
\end{squeezetable}
|
|
%%% %%% %%% %%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Conclusion}
|
|
\label{sec:ccl}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\acknowledgements{
|
|
We would like to thank Xavier Blase and Denis Jacquemin for insightful discussions.
|
|
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section*{Data availability statement}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
The data that supports the findings of this study are available within the article and its supplementary material.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\bibliography{sfBSE}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
\end{document}
|