done 1st iteration with theory

Manuscript/BSE_JPCL.tex

@@ -70,6 +70,7 @@
 \newcommand{\bI}{\boldsymbol{I}}
 \newcommand{\bA}{\boldsymbol{A}}
 \newcommand{\bH}{\boldsymbol{H}}
+\newcommand{\br}{\boldsymbol{r}}
 \newcommand{\bx}{\boldsymbol{x}}
 \newcommand{\bb}{\boldsymbol{b}}
 \newcommand{\bc}{\boldsymbol{c}}
@@ -206,9 +207,9 @@ Future directions of developments and improvements are also discussed.
 
 \noindent
 
-%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%% INTRODUCTION %%%
+\section{Introduction}
-%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 In its press release announcing the attribution of the 2013 Nobel prize in Chemistry to  Karplus,  Levitt and  Warshel, the Royal Swedish Academy of Sciences concluded by stating \textit{``Today the computer is just as important a tool for chemists as the test tube. 
 Simulations are so realistic that they predict the outcome of traditional experiments.''} \cite{Nobel_2003} 
@@ -217,18 +218,40 @@ Martin Karplus Nobel lecture moderated this bold statement, introducing his pres
 The study of neutral electronic excitations in condensed matter systems, from molecules to extended solids, has witnessed the development of a large number of such approximate methods with numerous applications to a large variety of fields, from the  prediction of the colour of precious metals and stones for jewellery, to the understanding, \eg, of the basic principles behind photovoltaics, photocatalysis or DNA damage under irradiation in the context of biology. 
 The present Perspective aims at describing the current status and upcoming challenges for the Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} that, while sharing many features with time-dependent density functional theory (TD-DFT), \cite{Runge_1984,Casida_1995,Dreuw_2005} including computational cost scaling with system size, relies on a different formalism, with specific difficulties but also potential solutions to known issues. \cite{Blase_2018}
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{History}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Originally developed in the framework of nuclear physics, \cite{Salpeter_1951}
+the use of the BSE formalism in condensed-matter physics emerged in the 1960's at the semi-empirical tight-binding level with the study of the optical properties of simple semiconductors.  \cite{Sham_1966,Strinati_1984,Delerue_2000}
+Three decades latter, the first \textit{ab initio} implementations, starting with small clusters \cite{Onida_1995,Rohlfing_1998} and extended  semiconductors and wide-gap insulators,  \cite{Albrecht_1997,Benedict_1998,Rohlfing_1999} 
+paved the way to the popularization in the solid-state physics community of the BSE formalism. 
+
+Following early applications to periodic polymers and molecules, [REFS] the   BSE formalism  gained much momentum in the quantum chemistry community with in particular several benchmarks \cite{Korbel_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017,Krause_2017,Gui_2018} on large molecular systems performed with the very same running parameters (geometries, basis sets) than the available reference higher-level calculations such as CC3. Such comparisons were grounded in the development of codes replacing the planewave solid-state physics paradigm by   well documented correlation-consistent Gaussian basis sets, together with adequate auxiliary bases when resolution-of-the-identity techniques were used. [REFS] 
+
+An important conclusion drawn from these calculations was that the quality of the BSE excitation energies are strongly correlated to the deviation of the preceding $GW$ HOMO-LUMO gap with the experimental (IP-AE) photoemission gap. Standard $G_0W_0$ calculations starting with Kohn-Sham (KS) eigenstates generated with (semi)local functionals yield   much larger HOMO-LUMO gaps than the input KS one, but still too small as compared to the experimental (AE-IP) value. Such an underestimation of the (IP-AE) gap leads to a similar underestimation of the lowest optical excitation energies.
+
+Such a residual HOMO-LUMO gap problem can be significantly improved by adopting exchange-correlation (xc) functionals with a tuned amount of exact exchange that yield a much improved KS HOMO-LUMO gap as a starting point for the $GW$ correction. \cite{Bruneval_2013,Rangel_2016,Knight_2016} Alternatively, self-consistent schemes, where corrected eigenvalues, and possibly orbitals, are reinjected in the construction of $G$ and $W$, have been shown to lead to a significant improvement of the quasiparticle energies in the case of molecular systems, with the advantage of significantly removing the dependence on the starting point functional. \cite{Rangel_2016,Kaplan_2016,Caruso_2016} As a result, BSE excitation singlet energies starting from such improved quasiparticle energies were found to be in much better agreement with reference calculations such as CC3. For sake of illustration,  an average  0.2 eV error was found for the well-known Thiel set comprising more than a hundred representative singlet excitations from a large variety of representative molecules. 
+\cite{Jacquemin_2015a,Bruneval_2015,Gui_2018,Krause_2017} This is equivalent to the best TD-DFT results obtained by scanning a large variety of global hybrid functionals with varying fraction of exact exchange. 
+
+A very remarkable success of the BSE formalism lies in the description of the charge-transfer (CT) excitations, a notoriously difficult problem for TD-DFT calculations adopting standard functionals. \cite{Dreuw_2004} Similar difficulties emerge as well in solid-state physics for semiconductors where extended Wannier  excitons are characterized by weakly overlapping electrons and holes, causing a dramatic deficit of spectral weight at low energy. \cite{Botti_2004} These difficulties can be ascribed to the lack of long-range electron-hole interaction with local XC functionals that can be cured through an exact exchange contribution, a solution that explains in particular the success of range-separated hybrids for the description of CT excitations. \cite{Stein_2009} The analysis of the screened Coulomb potential matrix elements in the BSE kernel (see Eqn.~\ref{Wmatel}) reveals that such long-range (non-local) electron-hole interactions are properly described, including in environments (solvents, molecular solid, etc.) where screening reduces the long-range electron-hole interactions. The success of the BSE formalism to treat CT excitations has been demonstrated in several studies, \cite{Blase_2011b,Baumeier_2012,Duchemin_2012,Sharifzadeh_2013,Cudazzo_2010,Cudazzo_2013} opening the way to important applications such as doping, photovoltaics or photocatalysis in organic systems. We now leave the description of successes to discuss difficulties and Perspectives.\\
+
+%%%%%%%%%%%%%%%%%%%%%%
+\section{Theory}
+%%%%%%%%%%%%%%%%%%%%%%
+
 The BSE formalism \cite{Salpeter_1951,Strinati_1988,Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999} belongs to the family of  Green's function many-body perturbation theories (MBPT) \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Reining_2017,ReiningBook} together with, for example, the algebraic-diagrammatic construction (ADC) techniques in quantum chemistry. \cite{Dreuw_2015} 
-While the density \titou{and density matrix} stand as the basic variables in DFT \titou{and Hartree-Fock}, Green's function MBPT takes the one-body  Green's function as the central quantity. The (time-ordered) one-body Green's function reads
+While the density stands as the basic variable in DFT, Green's function MBPT takes the one-body  Green's function as the central quantity. The (time-ordered) one-body Green's function reads
 \begin{equation}
     G(\bx t,\bx't') = -i \mel{N}{T \qty[ \Hat{\psi}(\bx t) \Hat{\psi}^{\dagger}(\bx't') ]}{N},
 \end{equation}
 where $\ket{N}$ is the $N$-electron ground-state wave function. 
 The operators $\Hat{\psi}(\bx t)$ and $\Hat{\psi}^{\dagger}(\bx't')$ remove and add an electron (respectively) in space-spin-time positions ($\bx t$) and ($\bx't'$), while \titou{$T$ is the time-ordering operator}. 
-For ($t > t'$) the one-body Green's function provides the amplitude of probability of finding, on top of the ground-state Fermi sea, an electron in ($\bx t$) that was previously introduced in ($\bx't'$), while for ($t < t'$) the propagation of a hole is monitored. 
+For ($t > t'$) the one-body Green's function provides the amplitude of probability of finding, on top of the ground-state Fermi sea, an electron in ($\bx t$) that was previously introduced in ($\bx't'$), while for ($t < t'$) the propagation of a hole is monitored.\\
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%===================================
-\paragraph{Charged excitations.}
+\subsection{Charged excitations}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%===================================
 A central property of the one-body Green's function is that its spectral representation presents poles at the charged excitation energies of the system
 \begin{equation}\label{eq:spectralG}
     G(\bx,\bx'; \omega ) = \sum_s \frac{ f_s(\bx) f^*_s(\bx') }{ \omega - \varepsilon_s + i \eta \, \text{sgn}(\varepsilon_s - \mu ) },  
@@ -236,145 +259,153 @@ A central property of the one-body Green's function is that its spectral represe
 where $\mu$ is the chemical potential, $\eta$ is a positive infinitesimal, $\varepsilon_s = E_s(N+1) - E_0(N)$ for $\varepsilon_s > \mu$, and $\varepsilon_s = E_0(N) - E_s(N-1)$ for $\varepsilon_s < \mu$. 
 The quantities $E_s(N+1)$ and $E_s(N-1)$ are the total energy of the $s$th excited state of the $(N+1)$ and $(N-1)$-electron systems, while $E_0(N)$ is the $N$-electron ground-state energy. 
 \titou{The $f_s$'s  are the so-called Lehmann amplitudes that reduce to one-body orbitals in the case of single-determinant many-body wave functions [more ??].}
-Unlike Kohn-Sham eigenvalues, the Green's function poles $\lbrace \varepsilon_s \rbrace$ are thus the proper \titou{charging} energies of the $N$-electron system, leading to well-defined ionization potentials and electronic affinities. Contrary to standard $\Delta$SCF techniques, the knowledge of $G$ provides the full ionization spectrum, as measured by direct and inverse photoemission, not only that associated with frontier orbitals. 
+Unlike KS eigenvalues, the Green's function poles $\lbrace \varepsilon_s \rbrace$ are thus the proper \titou{charging} energies of the $N$-electron system, leading to well-defined ionization potentials and electronic affinities. Contrary to standard $\Delta$SCF techniques, the knowledge of $G$ provides the full ionization spectrum, as measured by direct and inverse photoemission, not only that associated with frontier orbitals. 
 
 Using the equation-of-motion formalism for the creation/destruction operators, it can be shown formally that $G$ verifies
 \begin{equation}\label{eq:Gmotion}
-    \qty[  \pdv{}{t_1} - h({\bf r}_1) ] G(1,2) - \int d3 \; \Sigma(1,3) G(3,2),
+    \qty[  \pdv{}{t_1} - h(\br_1) ] G(1,2) - \int d3 \, \Sigma(1,3) G(3,2),
     = \delta(1,2)
 \end{equation}
 where we introduce the usual composite index, \eg, $1 \equiv (\bx_1,t_1)$. 
-Here, $h$ is the \titou{one-body Hartree Hamiltonian} and $\Sigma$ the so-called exchange-correlation self-energy operator.
+Here, $h$ is the \titou{one-body Hartree Hamiltonian} and $\Sigma$ the so-called xc self-energy operator.
 Using the spectral representation of $G$ [see Eq.~\eqref{eq:spectralG}], one obtains the familiar eigenvalue equation, \ie,
 \begin{equation}
-	h({\bf r}) f_s({\bf r}) + \int d{\bf r}' \; \Sigma({\bf r},{\bf r}'; \varepsilon_s ) f_s({\bf r}) = \varepsilon_s f_s({\bf r})
+	h(\br) f_s(\br) + \int d\br' \, \Sigma(\br,\br'; \varepsilon_s ) f_s(\br) = \varepsilon_s f_s(\br)
 \end{equation}
-which resembles formally the Kohn-Sham equation with the difference that the self-energy $\Sigma$ is non-local, energy dependent and non-hermitian.  
+which resembles formally the KS equation with the difference that the self-energy $\Sigma$ is non-local, energy dependent and non-hermitian.  
 The knowledge of $\Sigma$ allows thus to obtain the true addition/removal energies, namely the entire spectrum of occupied and virtual electronic energy levels, at the cost of solving a generalized one-body eigenvalue equation.  [INTRODUCE QUASIPARTICLES and OTHER solutions ??] \\
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%===================================
-\paragraph{The $GW$ self-energy.}
+\subsection{The $GW$ self-energy}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%===================================
-While the presented equations are formally exact, it remains to provide an expression for the exchange-correlation self-energy operator $\Sigma$.  This is where Green's function practical theories differ. Developed by Lars Hedin in 1965 with application  to the interacting homogeneous electron gas, \cite{Hedin_1965} the $GW$ approximation
+
-\cite{Aryasetiawan_1998,Farid_1999,Onida_2002,Ping_2013,Leng_2016,Golze_2019rev} follows the path of linear response by considering the variation   of $G$ with respect to an external perturbation. The obtained equation, when compared with the equation for the time-evolution of $G$ [Eqn.~\eqref{eq:Gmotion}], leads to a formal expression for the self-energy :
+While the equations reported above are formally exact, it remains to provide an expression for the xc self-energy operator $\Sigma$.  
-\begin{equation}
+This is where Green's function practical theories differ. 
-    \Sigma(1,2) = i \int d34 \; G(1,4) W(3,1^{+}) \Gamma(42,3)
+Developed by Lars Hedin in 1965 with application to the interacting homogeneous electron gas, \cite{Hedin_1965} the $GW$ approximation
-\end{equation}
+\cite{Aryasetiawan_1998,Farid_1999,Onida_2002,Ping_2013,Leng_2016,Golze_2019rev} follows the path of linear response by considering the variation of $G$ with respect to an external perturbation. 
-where $W$ is the dynamically screened Coulomb potential and $\Gamma$ a ``vertex" function that can be written as $\Gamma(12,3) = \delta(12)\delta(13) + \mathcal{O}(W)$ where   $\mathcal{O}(W)$ means a corrective term with leading linear order in terms of $W$. The neglect of the vertex leads to the so-called $GW$ approximation for $\Sigma$ that can be regarded as the lowest-order perturbation in terms of the screened Coulomb potential $W$ with :
+The obtained equation, when compared with the equation for the time-evolution of $G$ [Eq.~\eqref{eq:Gmotion}], leads to a formal expression for the self-energy
-\begin{align}
+\begin{equation}
-    W(1,2) &= v(1,2) + \int d34 \; v(1,2)   \chi_0(3,4) W(4,2)  \\
+    \Sigma(1,2) = i \int d34 \, G(1,4) W(3,1^{+}) \Gamma(42,3),
-      \chi_0(1,2) &= -i \int d34 \;  G(2,3) G(4,2) 
+\end{equation}
-\end{align}
+where $W$ is the dynamically-screened Coulomb potential and $\Gamma$ is a ``vertex" function that can be written as $\Gamma(12,3) = \delta(12)\delta(13) + \order{W}$, where $\order{W}$ means a corrective term with leading linear order in terms of $W$. 
-with   $\chi_0$ the well-known independent electron susceptibility and $v$ the bare Coulomb potential. In practice, the input $G$ and $\chi_0$ needed to buld $\Sigma$ are taken to be the best Green's function and susceptibility that can be easily calculated, namely the DFT (or HF) ones where the $\lbrace \varepsilon_s, f_s \rbrace$ of equation~\ref{spectralG} are taken to be DFT  Kohn-Sham (or HF) eigenstates. Taking then $( \Sigma^{\GW}-V^{\XC} )$ as a correction to the   $V^{\XC}$   DFT exchange correlation potential, a first-order correction to the input Kohn-Sham $\lbrace \varepsilon_n^{KS} \rbrace$ energies is obtained as follows:
+The neglect of the vertex leads to the so-called $GW$ approximation for $\Sigma$ that can be regarded as the lowest-order perturbation in terms of the screened Coulomb potential $W$ with
-\begin{equation}
+\begin{align}
-    \varepsilon_n^{\GW} = \varepsilon_n^{\KS} +
+	W(1,2) & = v(1,2) + \int d34 \, v(1,2)   \chi_0(3,4) W(4,2),
-    \langle \phi_n^{\KS} | \Sigma^{\GW}(\varepsilon_n^{\GW}) -V^{\XC}  | \phi_n^{\KS} \rangle
+	\\
-\end{equation}
+	\chi_0(1,2) & = -i \int d34 \,  G(2,3) G(4,2),
-Such an approach, where input Kohn-Sham energies are corrected to yield better electronic energy levels, is labeled the single-shot, or perturbative, $G_0W_0$ technique. This simple scheme was used in the early $GW$ studies of extended semiconductors and insulators, \cite{Strinati_1980,Hybertsen_1986,Godby_1988,Linden_1988} 
+\end{align}
-surfaces, \cite{Northrup_1991,Blase_1994,Rohlfing_1995} and 2D systems, \cite{Blase_1995} allowing to dramatically reduced the errors associated with Kohn-Sham eigenvalues in conjunction with  common local or gradient-corrected approximations to the exchange-correlation potential.
+where $\chi_0$ is the independent electron susceptibility and $v$ the bare Coulomb potential. 
-In particular, the well-known ``band gap" problem, \cite{Perdew_1983,Sham_1983} namely the underestimation of the occupied to unoccupied bands energy gap at the LDA Kohn-Sham level, was dramatically reduced, bringing the agreement with experiment to within a very few tenths of an eV  [REFS] with an $\mathcal{O}(N^4)$ computational cost scaling (see below). As another important feature   compared to other perturbative techniques, the $GW$ formalism can tackle finite size and periodic systems, and does not present any divergence in the limit of zero gap (metallic) systems. \cite{Campillo_1999} However, remaining a low order perturbative approach starting with mono-determinental mean-field solutions, it is not intended to explore strongly correlated systems. \cite{Verdozzi_1995} \\
+In practice, the input $G$ and $\chi_0$ needed to build $\Sigma$ are taken to be the ``best'' Green's function and susceptibility that can be easily calculated, namely the DFT or Hartree-Fock (HF) ones where the $\lbrace \varepsilon_p, f_p \rbrace$ of Eq.~\eqref{eq:spectralG} are taken to be KS (or HF) eigenstates. 
+Taking then $( \Sigma^{\GW}-V^{\XC} )$ as a correction to the DFT xc potential $V^{\XC}$, a first-order correction to the input KS energies $\lbrace \varepsilon_p^{\KS} \rbrace$ is obtained as follows:
+\begin{equation}
+    \varepsilon_p^{\GW} = \varepsilon_p^{\KS} +
+    \mel{\phi_p^{\KS}}{\Sigma^{\GW}(\varepsilon_p^{\GW}) - V^{\XC}}{\phi_p^{\KS}}.
+\end{equation}
+Such an approach, where input KS energies are corrected to yield better electronic energy levels, is labeled as the single-shot, or perturbative, $G_0W_0$ technique. This simple scheme was used in the early $GW$ studies of extended semiconductors and insulators, \cite{Strinati_1980,Hybertsen_1986,Godby_1988,Linden_1988} 
+surfaces, \cite{Northrup_1991,Blase_1994,Rohlfing_1995} and 2D systems, \cite{Blase_1995} allowing to dramatically reduced the errors associated with KS eigenvalues in conjunction with common local or gradient-corrected approximations to the xc potential.
+In particular, the well-known ``band gap" problem, \cite{Perdew_1983,Sham_1983} namely the underestimation of the occupied to unoccupied bands energy gap at the LDA KS level, was dramatically reduced, bringing the agreement with experiment to within a few tenths of an eV  [REFS] with an $\mathcal{O}(N^4)$ computational scaling (see below). 
+Another important feature compared to other perturbative techniques, the $GW$ formalism can tackle finite size and periodic systems, and does not present any divergence in the limit of zero gap (metallic) systems. \cite{Campillo_1999} 
+However, remaining a low order perturbative approach starting with single-determinant mean-field solutions, it is not intended to explore strongly correlated systems. \cite{Verdozzi_1995} \\
+
+%===================================
+\subsection{Neutral excitations}
+%===================================
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\paragraph{Neutral excitations.}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 While TD-DFT   starts with the  variation of the charge density $\rho$ with respect to an external local perturbation, the BSE formalism considers a generalized 4-points susceptibility that monitors the variation of the Green's function with respect to a non-local external perturbation: 
 \begin{equation}
-    \chi(1,2) \stackrel{\DFT}{=} \frac{ \partial \rho(1) }{\partial U(2) }
+    \chi(1,2) \stackrel{\DFT}{=} \pdv{\rho(1)}{U(2)}
-   \;\;  \rightarrow \;\; 
+   \quad  \rightarrow \quad
-    L(12,34) \stackrel{\BSE}{=} -i \frac{ \partial G(1,2) } { \partial U(3,4) }
+    L(12,34) \stackrel{\BSE}{=} -i \pdv{G(1,2)}{U(3,4)}.
 \end{equation}
 %with the relation $\chi(1,2) = L(11,22)$ since $\rho(1) = -iG(1,1^{+})$, as a first bridge between the TD-DFT and BSE worlds. 
-The equation of motion for $G$ (Eqn.~\ref{Gmotion}) can be reformulated in the form of a Dyson equation:
+The equation of motion for $G$ [Eq.~\ref{eq:Gmotion}] can be reformulated in the form of a Dyson equation
 \begin{equation}
- G = G_0 + G_0 \Sigma G
+	G = G_0 + G_0 \Sigma G,
 \end{equation}
-that relates the full (interacting) $G$ to  the Hartree $G_0$  obtained by replacing the $\lbrace \varepsilon_s, f_s \rbrace$ by the Hartree eigenvalues and eigenfunctions. 
+that relates the full (interacting) Green's function, $G$, to its Hartree version, $G_0$, obtained by replacing the $\lbrace \varepsilon_s, f_s \rbrace$ by the Hartree eigenvalues and eigenfunctions. 
-The derivation by $U$ of the Dyson equation yields :
+The derivative with respect to $U$ of the Dyson equation yields
-\begin{align}
+\begin{multline}\label{eq:DysonL}
-    L(12,34) &= L^0(12,34)  + \nonumber  \\
+    L(12,34) = L_0(12,34) 
-     & \int d5678 \; L^0(12,34) \Xi(5,6,7,8) L(78,34)
+    \\
-     \label{DysonL}
+     \int d5678 \, L_0(12,34) \Xi^{\BSE}(5,6,7,8) L(78,34),
-\end{align}
+\end{multline}
-with $L_0 = \partial G_0 / \partial U$  the Hartree 4-point susceptibility and:
+where $L_0 = \partial G_0 / \partial U$ is the Hartree 4-point susceptibility and
+\begin{equation}
+	\Xi^{\BSE}(5,6,7,8) = v(5,7) \delta(56) \delta(78)  + \pdv{\Sigma(5,6)}{G(7,8)}
+\end{equation}  
+is the so-called BSE kernel.
+
+This equation can be compared to its TD-DFT analog
+\begin{equation}
+    \chi(1,2) = \chi_0(1,2) + \int d34 \, \chi_0(1,3) \Xi^{\DFT}(3,4) \chi(4,2),
+\end{equation}
+where
+\begin{equation}
+    \Xi^{\DFT}(3,4) = v(3,4) + \pdv{ V^{\XC}(3)}{\rho(4)}
+\end{equation} 
+is the TD-DFT kernel.
+Plugging now the $GW$ self-energy, in a scheme that we label the BSE@$GW$ approach, leads to an approximation to the BSE kernel
 \begin{align*}
- \Xi(5,6,7,8) = v(5,7) \delta(56) \delta(78)  + \frac{ \partial \Sigma(5,6) }{ \partial G(7,8) }
+	\Xi^{\BSE}(5,6,7,8) = v(5,7) \delta(56) \delta(78)  -W(5,6) \delta(57) \delta(68 ),
 \end{align*}
-This equation can be compared to its TD-DFT analog: 
+where it is  traditional to neglect the derivative $( \partial W / \partial G)$ that introduces again higher orders in $W$. 
+Taking the static limit, \ie, $W(\omega=0)$, for the screened Coulomb potential, that replaces thus the static DFT xc kernel, and expressing Eq.~\eqref{eq:DysonL} in the standard product space $\lbrace \phi_i(\br) \phi_a(\br') \rbrace$ [where $(i,j)$ are occupied orbitals and $(a,b)$ are unoccupied orbitals), leads to an eigenvalue problem similar to the so-called Casida's equations in TD-DFT:
 \begin{equation}
-    \chi(1,2) = \chi_0(1,2) + \int d34 \; \chi_0(1,3) \Xi^{\DFT}(3,4) \chi(4,2)
+    \begin{pmatrix}
-\end{equation}
+		R & C  
-with $\Xi^{\DFT}$ the TD-DFT kernel :
+		\\
-\begin{equation}
-    \Xi^{\DFT}(3,4) = v(3,4) + \frac{ \partial V^{\XC}(3) }{ \partial \rho(4) }
-\end{equation} 
-Plugging now the $GW$ self-energy, in a scheme that we  label the BSE/$GW$ approach, leads to an approximation to the BSE kernel:
-\begin{align*}
- \Xi(5,6,7,8) = v(5,7) \delta(56) \delta(78)  -W(5,6) \delta(57) \delta(68 ) 
-\end{align*}
-where it is  traditional to neglect the derivative $( \partial W / \partial G)$ that introduces again higher orders in $W$. Taking the static limit $W(\omega=0)$ for the screened Coulomb potential, that replaces thus the static DFT exchange-correlation kernel, and expressing equation~\ref{DysonL} in the standard product space $\lbrace \phi_i({\bf r}) \phi_a({\bf r}') \rbrace$  where (i,j) and (a,b) index  occupied and virtual orbitals,  leads to an eigenvalue problem similar to the so-called Casida's equations in TD-DFT :
-\begin{equation}
-\left( 
-    \begin{matrix}
-R & C  \\
 		-C^*  & R^{*}
-\end{matrix}
+	\end{pmatrix}
-\right) \cdot 
+    \begin{pmatrix}
-\left( 
+		X^m
-    \begin{matrix}
+		\\
-X^{\lambda}  \\
+		Y^m
-Y^{\lambda}
+	\end{pmatrix}
-\end{matrix}
+	=
-\right) =
+	\Omega_m
-\Omega_{\lambda}
+    \begin{pmatrix}
-\left( 
+		X^m  
-    \begin{matrix}
+		\\
-X^{\lambda}  \\
+		Y^m
-Y^{\lambda}
+	\end{pmatrix},
-\end{matrix}
-\right)
 \end{equation}
-with electron-hole (e-h) eigenstates written:
+with electron-hole (e-h) eigenstates written as
 \begin{equation}
-    \psi_{\lambda}^{eh}(r_e,r_h) = \sum_{ia} \left( X_{ia}^{ \lambda}
+    \psi_{m}^{eh}(\br_e,\br_h) 
-    \phi_i({r}_h) \phi_a({r}_e) + Y_{ia}^{ \lambda}
+    = \sum_{ia} \qty[ X_{ia}^{m} \phi_i(\br_h) \phi_a(\br_e)
-    \phi_i({r}_e) \phi_a({r}_h) \right)
+	+ Y_{ia}^{m} \phi_i(\br_e) \phi_a(\br_h) ],
 \end{equation}
-where $\lambda$ index the electronic excitations. The $\lbrace \phi_{i/a} \rbrace$ are the input (Kohn-Sham) eigenstates used to build the $GW$ self-energy. The resonant part of the BSE Hamiltonian reads:
+where $\lambda$ index the electronic excitations. 
-\begin{align*}
+The $\lbrace \phi_{i/a} \rbrace$ are the input (KS) eigenstates used to build the $GW$ self-energy. 
-    R_{ai,bj} = \left( \varepsilon_a^{\GW} - \varepsilon_i^{\GW} \right) \delta_{ij} \delta_{ab} + \eta (ai|bj) - W_{ai,bj}
+The resonant and anti-resonant parts of the BSE Hamiltonian read
-\end{align*}
+\begin{gather}
-with $\eta=2,0$ for singlets/triplets and:
+    R_{ai,bj} = \qty( \varepsilon_a^{\GW} - \varepsilon_i^{\GW} ) \delta_{ij} \delta_{ab} + \eta (ai|bj) - W_{ai,bj},
-\begin{equation}
+    \\
-W_{ai,bj} = \int d{\bf r} d{\bf r}' 
+    C_{ai,bj} = \eta (ai|jb) - W_{ai,jb},
-  \phi_i({\bf r}) \phi_j({\bf r}) W({\bf r},{\bf r}'; \omega=0)
+\end{gather}
-  \phi_a({\bf r}') \phi_b({\bf r}')
+with $\eta=2,0$ for singlets/triplets and
-  \label{Wmatel}
+\begin{equation}\label{eq:Wmatel}
+	W_{ai,bj} = \iint d\br d\br' 
+	\phi_i(\br) \phi_j(\br) W(\br,\br'; \omega=0)
+	\phi_a(\br') \phi_b(\br'),
 \end{equation}
-where we notice that the 2 occupied (virtual) eigenstates are taken at the same space position, in contrast with the 
+where we notice that the two occupied (virtual) eigenstates are taken at the same position of space, in contrast with the 
-$(ai|bj)$ bare Coulomb term.  As compared to TD-DFT :
+$(ai|bj)$ bare Coulomb term.  
+As compared to TD-DFT,
 \begin{itemize}
-    \item the $GW$ quasiparticle energies $\lbrace \varepsilon_{i/a}^{\GW} \rbrace$ replace the Kohn-Sham eigenvalues
+    \item the $GW$ quasiparticle energies $\lbrace \varepsilon_{i/a}^{\GW} \rbrace$ replace the KS eigenvalues
-    \item the non-local screened Coulomb matrix elements replaces the DFT exchange-correlation kernel.
+    \item the non-local screened Coulomb matrix elements replaces the DFT xc kernel.
 \end{itemize}    
-We emphasise that these equations can be solved at exactly the same cost as the standard TD-DFT equations once the quasiparticle energies and screened Coulomb potential $W$ are inherited from preceding $GW$ calculations. This defines the standard (static) BSE/$GW$scheme that we discuss in this Perspective, emphasizing its pros and cons.  \\
+We emphasise that these equations can be solved at exactly the same cost as the standard TD-DFT equations once the quasiparticle energies and screened Coulomb potential $W$ are inherited from preceding $GW$ calculations. This defines the standard (static) BSE@$GW$ scheme that we discuss in this \textit{Perspective}, emphasizing its pros and cons.  \\
 
-%% BSE historical
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Challenges} 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-Originally developed in the framework of nuclear physics, \cite{Salpeter_1951}
+%==========================================
-the use of the BSE formalism in condensed-matter physics emerged in the 60s at the semi-empirical tight-binding level with the study of the optical properties of simple semiconductors.  \cite{Sham_1966,Strinati_1984,Delerue_2000}
+\subsection{The computational challenge} 
-Three decades latter, the first \textit{ab initio} implementations, starting with small clusters \cite{Onida_1995,Rohlfing_1998} and extended  semiconductors and wide-gap insulators,  \cite{Albrecht_1997,Benedict_1998,Rohlfing_1999} 
+%==========================================
-paved the way to the popularization in the solid-state physics community of the BSE formalism. 
+As emphasized above, the BSE eigenvalue equation in the occupied-to-virtual product space is formally equivalent to that of TD-DFT or TD-HF. Searching iteratively for the lowest eigenstates presents the same $\mathcal{O}(N^4)$ matrix-vector multiplication  computational cost within BSE and TD-DFT. Concerning the construction of the BSE Hamiltonian, it is no more expensive than building the TD-DFT one with hybrid functionals, reducing again to $\mathcal{O}(N^4)$ operations with standard resolution-of-identity techniques. At the price of sacrifying the knowledge of the eigenvectors, the BSE absorption spectrum can be known with $\mathcal{O}(N^3)$ operations using iterative techniques. \cite{Ljungberg_2015} With the same restriction on the eigenvectors, a time-propagation approach, similar to that implemented for TD-DFT, \cite{Yabana_1996} combined with stochastic techniques to reduce the cost of building the BSE Hamiltonian matrix elements, allows quadratic scaling with systems size. \cite{Rabani_2015}
-
-Following early applications to periodic polymers and molecules, [REFS] the   BSE formalism  gained much momentum in the quantum chemistry community with in particular several benchmarks \cite{Korbel_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017,Krause_2017,Gui_2018} on large molecular systems performed with the very same running parameters (geometries, basis sets) than the available reference higher-level calculations such as CC3. Such comparisons were grounded in the development of codes replacing the planewave solid-state physics paradigm by   well documented correlation-consistent Gaussian basis sets, together with adequate auxiliary bases when resolution-of-the-identity techniques were used. [REFS] 
-
-An important conclusion drawn from these calculations was that the quality of the BSE excitation energies are strongly correlated to the deviation of the preceding $GW$ HOMO-LUMO gap with the experimental (IP-AE) photoemission gap. Standard $G_0W_0$ calculations starting with Kohn-Sham eigenstates generated with (semi)local functionals yield   much larger HOMO-LUMO gaps than the input Kohn-Sham one, but still too small as compared to the experimental (AE-IP) value. Such an underestimation of the (IP-AE) gap leads to a similar underestimation of the lowest optical excitation energies.
-
-Such a residual HOMO-LUMO gap problem can be significantly improved by adopting exchange-correlation functionals with a tuned amount of exact exchange that yield a much improved Kohn-Sham HOMO-LUMO gap as a starting point for the $GW$ correction. \cite{Bruneval_2013,Rangel_2016,Knight_2016} Alternatively, self-consistent schemes, where corrected eigenvalues, and possibly orbitals, are reinjected in the construction of $G$ and $W$, have been shown to lead to a significant improvement of the quasiparticle energies in the case of molecular systems, with the advantage of significantly removing the dependence on the starting point functional. \cite{Rangel_2016,Kaplan_2016,Caruso_2016} As a result, BSE excitation singlet energies starting from such improved quasiparticle energies were found to be in much better agreement with reference calculations such as CC3. For sake of illustration,  an average  0.2 eV error was found for the well-known Thiel set comprising more than a hundred representative singlet excitations from a large variety of representative molecules. 
-\cite{Jacquemin_2015a,Bruneval_2015,Gui_2018,Krause_2017} This is equivalent to the best TD-DFT results obtained by scanning a large variety of global hybrid functionals with varying fraction of exact exchange. 
-
-A very remarkable success of the BSE formalism lies in the description of the charge-transfer (CT) excitations, a notoriously difficult problem for TD-DFT calculations adopting standard functionals. \cite{Dreuw_2004} Similar difficulties emerge as well in solid-state physics for semiconductors where extended Wannier  excitons are characterized by weakly overlapping electrons and holes, causing a dramatic deficit of spectral weight at low energy. \cite{Botti_2004} These difficulties can be ascribed to the lack of long-range electron-hole interaction with local XC functionals that can be cured through an exact exchange contribution, a solution that explains in particular the success of range-separated hybrids for the description of CT excitations. \cite{Stein_2009} The analysis of the screened Coulomb potential matrix elements in the BSE kernel (see Eqn.~\ref{Wmatel}) reveals that such long-range (non-local) electron-hole interactions are properly described, including in environments (solvents, molecular solid, etc.) where screening reduces the long-range electron-hole interactions. The success of the BSE formalism to treat CT excitations has been demonstrated in several studies, \cite{Blase_2011b,Baumeier_2012,Duchemin_2012,Sharifzadeh_2013,Cudazzo_2010,Cudazzo_2013} opening the way to important applications such as doping, photovoltaics or photocatalysis in organic systems. We now leave the description of successes to discuss difficulties and Perspectives.\\
-
-\noindent {\textbf{The computational challenge.}} As emphasized above, the BSE eigenvalue equation in the occupied-to-virtual product space is formally equivalent to that of TD-DFT or TD-Hartree-Fock. Searching iteratively for the lowest eigenstates presents the same $\mathcal{O}(N^4)$ matrix-vector multiplication  computational cost within BSE and TD-DFT. Concerning the construction of the BSE Hamiltonian, it is no more expensive than building the TD-DFT one with hybrid functionals, reducing again to $\mathcal{O}(N^4)$ operations with standard resolution-of-identity techniques. At the price of sacrifying the knowledge of the eigenvectors, the BSE absorption spectrum can be known with $\mathcal{O}(N^3)$ operations using iterative techniques. \cite{Ljungberg_2015} With the same restriction on the eigenvectors, a time-propagation approach, similar to that implemented for TD-DFT, \cite{Yabana_1996} combined with stochastic techniques to reduce the cost of building the BSE Hamiltonian matrix elements, allows quadratic scaling with systems size. \cite{Rabani_2015}
 
 In practice, the main bottleneck for standard BSE calculations as compared to TD-DFT resides in the preceding $GW$ calculations that scale  as   $\mathcal{O}(N^4)$ with system size using  PWs  or RI techniques, but with a rather large prefactor. 
 %%Such a cost is mainly associated with calculating the free-electron susceptibility with its entangled summations over occupied and virtual states. 
@@ -382,9 +413,9 @@ In practice, the main bottleneck for standard BSE calculations as compared to TD
 The field of low-scaling $GW$ calculations is however witnessing significant advances. While the sparcity of ..., \cite{Foerster_2011,Wilhelm_2018} efficient real-space-grid  and time  techniques are blooming, \cite{Rojas_1995,Liu_2016} borrowing in particular the well-known Laplace transform approach used in quantum chemistry. \cite{Haser_1992}
 Together with a stochastic  sampling of virtual states, this family of techniques allow to set up linear scaling $GW$ calculations. \cite{Vlcek_2017} The separability of occupied and virtual states summations lying at the heart of these approaches are now blooming in quantum chemistry withing the Interpolative Separable Density Fitting (ISDF) approach applied to calculating with cubic scaling the susceptibility needed in RPA and $GW$ calculations. \cite{Lu_2017,Duchemin_2019,Gao_2020}   These ongoing developments pave the way to applying the $GW$/BSE formalism  to systems comprising several hundred atoms on standard laboratory clusters.  \\
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%==========================================
-\paragraph{The Triplet Instability Challenge.}
+\subsection{The triplet instability challenge}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%==========================================
 The analysis of the singlet-triplet splitting is central to numerous applications such as singlet fission, thermally activated delayed fluorescence (TADF) or 
 stability analysis of restricted closed-shell solutions at the HF \cite{Seeger_1977} and TD-DFT \cite{Bauernschmitt_1996} levels. 
 contaminating as well TD-DFT calculations with popular range-separated hybrids (RSH) that generally contains a large fraction of exact exchange in the long-range. \cite{Sears_2011}
@@ -394,33 +425,18 @@ benchmarks \cite{Jacquemin_2017b,Rangel_2017}
 
 a first cure was offered by hybridizing TD-DFT and BSE, namely adding to the BSE kernel the correlation part of the underlying DFT functional used to build the susceptibility and resulting screened Coulomb potential $W$.  \cite{Holzer_2018b}
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%==========================================
-\paragraph{The challenge of analytic gradients.}
+\subsection{The challenge of analytic gradients}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%==========================================
 An additional issue concerns the formalism taken to calculate the ground-state energy for a given atomic configuration. Since the BSE formalism presented so far concerns the calculation of the electronic excitations, namely the difference of energy between the GS and the ES, gradients of the ES absolute energy require 
 
 This points to another direction for the BSE foramlism, namely the calculation of GS total energy with the correlation energy calculated at the BSE level. Such a task was performed by several groups using in particular the adiabatic connection fluctuation-dissipation theorem (ACFDT), focusing in particular on small dimers. \cite{Olsen_2014,Holzer_2018b,Li_2020,Loos_2020}
 
-
+%==========================================
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{The double excitation challenge}
-\paragraph{Dynamical kernels and multiple excitations.}
+%==========================================
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\cite{Zhang_2013}
-
-
-\noindent {\textbf{Core-level spectroscopy.}}. \\
-XANES, 
-\cite{Olovsson_2009,Vinson_2011}
-
-
-diabatization and conical intersections \cite{Kaczmarski_2010}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\paragraph{The Concept of dynamical properties.}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 As a chemist, it is maybe difficult to understand the concept of dynamical properties, the motivation behind their introduction, and their actual usefulness.
-Here, we will try to give a pedagogical example showing the importance of dynamical quantities and their main purposes. \cite{Romaniello_2009,Sangalli_2011,ReiningBook}
+Here, we will try to give a pedagogical example showing the importance of dynamical quantities and their main purposes. \cite{Romaniello_2009,Sangalli_2011,Zhang_2013,ReiningBook}
 To do so, let us consider the usual chemical scenario where one wants to get the neutral excitations of a given system.
 In most cases, this can be done by solving a set of linear equations of the form
 \begin{equation}
@@ -485,8 +501,17 @@ In such a way, the operator $\Tilde{\bA}_1$ is made linear again by removing its
 This approximation comes with a heavy price as the number of solutions provided by the system of equations \eqref{eq:non_lin_sys} has now been reduced from $K$ to $K_1$.
 Coming back to our example, in the static approximation, the operator $\Tilde{\bA}_1$ built in the basis of single excitations cannot provide double excitations anymore, and the only $K_1$ excitation energies are associated with single excitations.
 
+\noindent {\textbf{Core-level spectroscopy.}}. \\
+XANES, 
+\cite{Olovsson_2009,Vinson_2011}
+
+
+diabatization and conical intersections \cite{Kaczmarski_2010}
+
+
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\paragraph{Conclusion.}
+\section{Conclusion.}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Here goes the conclusion.
 

Manuscript/xbbib.bib

@@ -1,7 +1,7 @@
 %% This BibTeX bibliography file was created using BibDesk.
 %% http://bibdesk.sourceforge.net/
 
-%% Created for Pierre-Francois Loos at 2020-05-13 10:12:24 +0200 
+%% Created for Pierre-Francois Loos at 2020-05-13 11:30:11 +0200 
 
 
 %% Saved with string encoding Unicode (UTF-8) 

TOC/TOC.tex

@@ -0,0 +1,77 @@
+\documentclass{standalone}
+\usepackage{amsmath,amssymb,amsfonts,pgfpages,graphicx,subfigure,xcolor,bm,multirow,microtype,wasysym,multimedia,hyperref,tabularx,amscd,pgfgantt,mhchem,physics}
+\usetikzlibrary{shapes.gates.logic.US,trees,positioning,arrows}
+
+\usepackage{tgchorus}
+\usepackage[T1]{fontenc}
+
+\begin{document}
+
+\begin{tikzpicture}
+	\begin{scope}[very thick,
+		node distance=5cm,on grid,>=stealth',
+		block1/.style={rectangle,draw,fill=violet!20},
+		block2/.style={rectangle,draw,fill=blue!20},
+		comp1/.style={rectangle,draw,fill=purple!40},
+		comp2/.style={circle,draw,fill=green!40},
+		comp3/.style={rectangle,fill=white!40}]
+		
+		\node [block1, text width=7cm, align=center]	(KS)		
+		{\textbf{\LARGE Kohn-Sham DFT}
+			$$
+				\qty[ -\frac{\nabla^2}{2} + v_\text{ext} + V^{\text{Hxc}} ] \phi_p^{\text{KS}} = \varepsilon^{\text{KS}}_p \phi_p^{\text{KS}}
+			$$	
+		};									
+
+		\node [block2, text width=7cm, align=center]	(GW)	[below=of KS, yshift=2cm]				
+		{\textbf{\LARGE $GW$ approximation} 
+			$$   
+				\varepsilon_p^{GW} = \varepsilon_p^{\text{KS}} +
+				\mel{\phi_p^{\text{KS}}}{\Sigma^{GW}(\varepsilon_p^{GW}) - V^{\text{xc}}}{\phi_p^{\text{KS}}}
+			$$
+
+		};												
+
+		\node [comp1, text width=7cm, align=center]	(BSE)	[below=of GW, yshift=2cm]									
+		{\textbf{\LARGE Bethe-Salpeter equation}
+			$$
+				\begin{pmatrix}
+					R & C  \\
+					-C^*  & R^{*}
+				\end{pmatrix}
+				\begin{matrix}
+					X_m  \\
+					Y_m
+				\end{matrix}
+				=
+				\Omega_{m}
+				\begin{pmatrix}
+					X_m  \\
+					Y_m
+				\end{pmatrix}
+			$$
+		};
+								
+
+		\node [comp2, align=center] 					(phys)	[below=of BSE, xshift=-4cm]										
+		{\LARGE Molecules};
+
+		\node [comp2, align=center] 					(chem)	[below=of BSE, xshift=0cm]			
+		{\LARGE Materials};
+		
+		\node [comp2, align=center] 					(bio)	[below=of BSE, xshift=4cm]		
+		{\LARGE Clusters};
+		
+		\path 
+		(KS) 			edge [->,color=black] 			node [right,black] {Fundamental gap}		(GW)
+		(GW) 			edge [->,color=black] 			node [right,black] {Excitonic effect}	(BSE)
+		(BSE) 			edge [->,color=black] 			node [above,black] {}	(phys)
+		(BSE) 			edge [->,color=black] 			node [above,black] {}	(chem)
+		(BSE) 			edge [->,color=black] 			node [above,black] {}	(bio)
+		;
+
+		
+	\end{scope}
+\end{tikzpicture}
+
+\end{document}