resubmission

Manuscript/BSE_JPCL.bib

@@ -6,12 +6,6 @@
 
 %% Saved with string encoding Unicode (UTF-8) 
 
-
-
-@article{cite-key,
-	Date-Added = {2020-06-16 09:59:45 +0200},
-	Date-Modified = {2020-06-16 09:59:45 +0200}}
-
 @article{Packer_1996,
 	Author = {Packer, M. K. and Dalskov, E. K. and Enevoldsen, T. and Jensen, H. J. and Oddershede, J.},
 	Date-Added = {2020-06-08 21:57:16 +0200},
@@ -2657,18 +2651,6 @@
 	Year = {2012},
 	Bdsk-Url-1 = {https://doi.org/10.1002/cphc.201100200}}
 
-@article{Holzer_2018,
-	Author = {Christof Holzer and Xin Gui and Michael E. Harding and Georg Kresse and Trygve Helgaker and Wim Klopper},
-	Date-Added = {2020-01-04 20:49:55 +0100},
-	Date-Modified = {2020-02-05 20:58:26 +0100},
-	Doi = {10.1063/1.5047030},
-	Journal = {J. Chem. Phys.},
-	Pages = {144106},
-	Title = {Bethe--Salpeter Correlation Energies of Atoms and Molecules},
-	Volume = {149},
-	Year = {2018},
-	Bdsk-Url-1 = {https://doi.org/10.1063/1.5047030}}
-
 @article{Olsen_2014,
 	Author = {T. Olsen and K. S. Thygesen},
 	Date-Added = {2020-01-04 20:48:50 +0100},

Manuscript/BSE_JPCL.tex

@@ -200,12 +200,12 @@
 %%%%%%%%%%%%%%%%
 
 \begin{abstract}
-The many-body Green's function Bethe-Salpeter equation (BSE) formalism is steadily asserting itself as a new efficient and accurate tool in the ensemble of computational methods available to chemists in order to predict optical excitations in molecular systems.
+\titou{The Bethe-Salpeter equation (BSE) formalism is steadily asserting itself as a new efficient and accurate tool in the ensemble of computational methods available to chemists in order to predict optical excitations in molecular systems.
-In particular, the combination of the so-called $GW$ approximation of many-body perturbation theory, giving access to reliable ionization energies and electron affinities, and the BSE formalism, able to model UV/Vis spectra by catching excitonic effects, has shown to provide accurate singlet excitation energies in many chemical scenarios with a typical error of $0.1$--$0.3$ eV. 
+In particular, the combination of the so-called $GW$ approximation, giving access to reliable ionization energies and electron affinities, and the BSE formalism, able to model UV/Vis spectra, has shown to provide accurate singlet excitation energies with a typical error of $0.1$--$0.3$ eV. 
-With a similar computational cost as time-dependent density-functional theory (TD-DFT), the BSE formalism is able to provide an accuracy on par with the most accurate global and range-separated hybrid functionals without the unsettling choice of the exchange-correlation functional, resolving further known issues (\textit{e.g.}, charge-transfer excitations) and offering a well-defined path to dynamical kernels.
+With a similar computational cost as time-dependent density-functional theory (TD-DFT), BSE is able to provide an accuracy on par with the most accurate global and range-separated hybrid functionals without the unsettling choice of the exchange-correlation functional, resolving further known issues (\textit{e.g.}, charge-transfer excitations).
-In this \textit{Perspective} article, we provide a historical overview of the BSE formalism, with a particular focus on its condensed-matter roots.
+In this \textit{Perspective} article, we provide a historical overview of BSE, with a particular focus on its condensed-matter roots.
-We also propose a critical review of its strengths and weaknesses in different chemical situations.
+We also propose a critical review of its strengths and weaknesses in different chemical situations.}
-Future directions of developments and improvements are also discussed.
+%\trashPFL{Future directions of developments and improvements are also discussed.}
 \end{abstract}
 
 
@@ -223,9 +223,8 @@ Future directions of developments and improvements are also discussed.
 \noindent
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Introduction}
+\paragraph{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} 
 Martin Karplus' Nobel lecture moderated this statement, introducing his presentation by a 1929 quote from Dirac emphasizing that laws of quantum mechanics are \textit{``much too complicated to be soluble''}, urging scientists to develop \textit{``approximate practical methods''}. This is where the electronic structure community stands, attempting to develop robust approximations to study with increasing accuracy the properties of ever more complex systems. 
@@ -234,9 +233,8 @@ The present \textit{Perspective} aims at describing the current status and upcom
 \\
 
 %%%%%%%%%%%%%%%%%%%%%%
-\section{Theory}
+\paragraph{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,Onida_2002,ReiningBook} together with, for example, the algebraic-diagrammatic construction (ADC) techniques \cite{Dreuw_2015} or the polarization propagator approaches (like SOPPA\cite{Packer_1996})  in quantum chemistry.
 While the one-body density stands as the basic variable in density-functional theory (DFT), \cite{Hohenberg_1964,Kohn_1965} the pillar of Green's function MBPT is the (time-ordered) one-body Green's function
 \begin{equation}
@@ -247,16 +245,16 @@ The operators $\Hat{\psi}(\bx t)$ and $\Hat{\psi}^{\dagger}(\bx't')$ remove and
 For $t > t'$, $G$ provides the amplitude of probability of finding, on top of the ground-state Fermi sea (\ie, higher in energy than the highest-occupied energy level, also known as Fermi level), an electron in ($\bx t$) that was previously introduced in ($\bx't'$), while for $t < t'$ the propagation of an electron hole (often simply called a hole) is monitored. 
 
 This definition indicates that the one-body Green's function is well suited to obtain ``charged excitations", more commonly labeled as electronic energy levels, as obtained, \eg, in a direct or inverse photo-emission experiments where an electron is ejected or added to the $N$-electron system. 
-In particular, and as opposed to Kohn-Sham (KS) DFT, the Green's function formalism offers a more rigorous and systematically improvable path for the obtention of the ionization potential $I^\Nel = E_0^{\Nel-1} - E_0^\Nel$, the electronic affinity $A^\Nel = E_0^{\Nel+1} - E_0^\Nel$, and the experimental (photoemission) fundamental gap
+In particular, and as opposed to Kohn-Sham (KS) DFT, the Green's function formalism offers a more rigorous and systematically improvable path for the obtention of the ionization potential $I^\Nel = E_0^{\Nel-1} - E_0^\Nel$, the electronic affinity $A^\Nel = E_0^{\Nel} - E_0^{\Nel+1}$, and the experimental (photoemission) fundamental gap
 \begin{equation}\label{eq:IPAEgap}
-	\EgFun = I^\Nel - A^\Nel,
+	\EgFun = I^\Nel - A^\Nel
 \end{equation}
 of the $\Nel$-electron system, where $E_0^\Nel$ corresponds to its ground-state energy.
 Since these energy levels are key input quantities for the subsequent BSE calculation, we start by discussing these in some details.
 \\
 
 %===================================
-\subsection{Charged excitations}
+\paragraph{Charged excitations.}
 %===================================
 A central property of the one-body Green's function is that its frequency-dependent (\ie, dynamical) spectral representation has poles at the charged excitation energies (\ie, the ionization potentials and electron affinities) of the system
 \begin{equation}\label{eq:spectralG}
@@ -278,16 +276,15 @@ Here, $\delta$ is Dirac's delta function, $h$ is the one-body Hartree Hamiltonia
 Using the spectral representation of $G$ [see Eq.~\eqref{eq:spectralG}],
 dropping spin variables for simplicity, one gets the familiar eigenvalue equation, \ie,
 \begin{equation}
-	h(\br) f_s(\br) + \int d\br' \, \Sigma(\br,\br'; \varepsilon_s ) f_s(\br) = \varepsilon_s f_s(\br),
+	\titou{h(\br) f_s(\br) + \int d\br' \, \Sigma(\br,\br'; \varepsilon_s ) f_s(\br') = \varepsilon_s f_s(\br),}
 \end{equation}
 which formally resembles the KS equation \cite{Kohn_1965} with the difference that the self-energy $\Sigma$ is non-local, energy-dependent and non-hermitian. 
 The knowledge of $\Sigma$ allows to access 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. 
 \\
 
 %===================================
-\subsection{The $GW$ self-energy}
+\paragraph{The $GW$ self-energy.}
 %===================================
-
 While the equations reported above are formally exact, it remains to provide an expression for the xc 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
@@ -296,20 +293,28 @@ The resulting equation, when compared with the equation for the time-evolution o
 \begin{equation}\label{eq:Sig}
 	\Sigma(1,2) = i \int d34 \, G(1,4) W(3,1^{+}) \Gamma(42,3),
 \end{equation}
-where $W$ is the dynamically-screened Coulomb potential and $\Gamma$ is the so-called ``vertex" function.
+where $W$ is the dynamically-screened Coulomb potential and $\Gamma$ is the so-called ``vertex" function. \titou{The notation $1^+$ means that the time $t_1$ is taken at $t_1^{+} = t_1 + 0^+$ for sake of causality, where $0^+$ is a positive infinitesimal.}
 %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$. 
 The neglect of the vertex, \ie, $\Gamma(42,3) = \delta(23) \delta(24)$, leads to the so-called $GW$ approximation of the self-energy 
 \begin{equation}\label{eq:SigGW}
 	\Sigma^{\GW}(1,2) = i \, G(1,2) W(2,1^{+}),
 \end{equation}
 that can be regarded as the lowest-order perturbation in terms of the screened Coulomb potential $W$ with
+\begin{subequations}
 \begin{gather}
-	W(1,2) = v(1,2) + \int d34 \, v(1,2) \chi_0(3,4) W(4,2),
+	W(1,2) = v(1,2) + \int d34 \, v(1,\titou{3}) \chi_0(3,4) W(4,2),
 	\label{eq:defW}
 	\\
-	\chi_0(1,2) = -i \int d34 \, G(2,3) G(4,2),
+  \titou{ 	\chi_0(1,2) = -i   G(1,2^{+}) G(2,1^{+}), }
 \end{gather}
+\end{subequations}
 where $\chi_0$ is the independent electron susceptibility and $v$ the bare Coulomb potential. 
+\titou{Equation~\eqref{eq:defW} can be recast as
+\begin{gather}
+	W(1,2) = v(1,2) + \int d34 \, v(1,\titou{3}) \chi(3,4) v(4,2),
+	\label{eq:defW2}
+\end{gather}
+where $\chi$ is the interacting susceptibility. In this latter expression, $(v{\chi}v)$ represents the field created in $(2)$ by the charge rearrangement of the $N$-electron system generated by a (unit) charge added in $(1)$. As such, this term contains the effect of dielectric screening (or polarization in a quantum chemist language). As in a standard $\Delta$SCF calculation, the $GW$ formalism contains the response of the $N$-electron system to an electron added (removed) to any virtual (occupied) molecular orbital, but without the restriction that only frontier orbitals can be tackled. This explains that the $GW$ one-electron energies are proper addition/removal energies.}
 
 %%% FIG 1 %%%
 \begin{figure}[ht]
@@ -354,7 +359,7 @@ The existence of a well defined quasiparticle energy requires a solution with a
 
 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} and
-surfaces, \cite{Northrup_1991,Blase_1994,Rohlfing_1995} allowing to dramatically reduced the errors associated with KS eigenvalues in conjunction with common local or gradient-corrected approximations to the xc potential.
+surfaces, \cite{Northrup_1991,Blase_1994,Rohlfing_1995} allowing to dramatically \titou{reduce} 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 local-density approximation (LDA) KS level, was dramatically reduced, bringing the agreement with experiment to within a few tenths of an eV with a computational cost scaling quartically with the system size (see below). A compilation of data for $G_0W_0$ applied to extended inorganic semiconductors can be found in Ref.~\citenum{Shishkin_2007}.
 
 Although $G_0W_0$ provides accurate results (at least for weakly/moderately correlated systems), it is strongly starting-point dependent due to its perturbative nature. 
@@ -390,10 +395,10 @@ However, remaining a low-order perturbative approach starting with a single-dete
 %%% %%% %%%
 
 %===================================
-\subsection{Neutral excitations}
+\paragraph{Neutral excitations.}
 %===================================
 Like TD-DFT, BSE deals with the calculations of optical (or neutral) excitations, as measured by optical (\eg, absorption) spectroscopy, 
-However, while TD-DFT starts with the variation of the charge density $\rho(1)$ with respect to an external local perturbation $U(1)$, the BSE formalism considers a generalized 4-points susceptibility, or two-particle correlation function, that monitors the variation of the one-body Green's function $G(1,1')$ with respect to a non-local external perturbation $U(2,2')$: \cite{Strinati_1988}
+However, while TD-DFT starts with the variation of the charge density $\rho(1)$ with respect to an external local perturbation $U(1)$, the BSE formalism considers a generalized \titou{four-point} susceptibility, or two-particle correlation function, that monitors the variation of the one-body Green's function $G(1,1')$ with respect to a non-local external perturbation $U(2,2')$: \cite{Strinati_1988}
 \begin{equation}
 	\chi(1,2) \stackrel{\DFT}{=} \pdv{\rho(1)}{U(2)}
 	\quad  \rightarrow \quad
@@ -492,9 +497,8 @@ This defines the standard (static) BSE@$GW$ scheme that we discuss in this \text
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Historical overview}
+\paragraph{Historical overview.}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 Originally developed in the framework of nuclear physics, \cite{Salpeter_1951} the BSE formalism has emerged in condensed-matter physics around the 1960's at the tight-binding level with the study of the optical properties of simple semiconductors. \cite{Sham_1966,Strinati_1984,Delerue_2000}
 Three decades later, the first \textit{ab initio} implementations, starting with small clusters, \cite{Onida_1995,Rohlfing_1998} extended semiconductors, and wide-gap insulators \cite{Albrecht_1997,Benedict_1998,Rohlfing_1999b} paved the way to the popularization in the solid-state physics community of the BSE formalism. 
 
@@ -515,6 +519,7 @@ but still too small as compared to the experimental value, \ie,
 \begin{equation}
 	\Eg^{\KS} \ll \Eg^{G_0W_0} < \EgFun.
 \end{equation}
+\titou{Such a residual discrepancy has been attributed by several authors to ``overscreening", namely the effect associated with building the susceptibility $\chi$ based on a grossly underestimated (KS) band gap. This leads to a spurious enhancement of the screening or polarization and, consequently, to an underestimated $G_0W_0$ gap as compared to the (exact) fundamental gap. More prosaically, the $G_0W_0$ approach is constructed as a first-order perturbation theory, so by correcting a very ``bad" zeroth-order KS system one cannot expect to obtain an accurate corrected gap.}
 Such an underestimation of the fundamental gap leads to a similar underestimation of the optical gap $\EgOpt$, \ie, the lowest optical excitation energy:
 \begin{equation}
 	\EgOpt = E_1^{\Nel} - E_0^{\Nel} = \EgFun + \EB,
@@ -529,12 +534,12 @@ For sake of illustration, an average error of $0.2$ eV was found for the well-kn
 This is equivalent to the best TD-DFT results obtained by scanning a large variety of hybrid functionals with various amounts of exact exchange. 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Successes \& Challenges} 
+%\section{Successes \& Challenges} 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 %==========================================
-\subsection{Charge-transfer excited states} 
+\paragraph{Charge-transfer excited states.} 
 %==========================================
 A very remarkable success of the BSE formalism lies in the description of charge-transfer (CT) excitations, a notoriously difficult problem for TD-DFT adopting standard (semi-)local functionals. \cite{Dreuw_2004} 
 Similar difficulties emerge in solid-state physics for semiconductors where extended Wannier excitons, characterized by weakly overlapping electrons and holes (Fig.~\ref{fig:CTfig}), cause a dramatic deficit of spectral weight at low energy. \cite{Botti_2004} 
@@ -555,9 +560,8 @@ The success of the BSE formalism to treat CT excitations has been demonstrated i
 %%% %%% %%%
 
 %==========================================
-\subsection{Combining BSE with PCM and QM/MM models} 
+\paragraph{Combining BSE with PCM and QM/MM models.} 
 %==========================================
-
 The ability to account for the effect on the excitation energies of an electrostatic and dielectric environment (an electrode, a solvent, a molecular interface\ldots) is an important step towards the description of realistic systems. 
 Pioneering BSE studies demonstrated, for example, the large renormalization of charged and neutral excitations in molecular systems and nanotubes close to a metallic electrode or in bundles. \cite{Lastra_2011,Rohlfing_2012,Spataru_2013}
 Recent attempts to merge the $GW$ and BSE formalisms with model polarizable environments at the PCM or QM/MM levels
@@ -582,7 +586,7 @@ We now leave the description of successes to discuss difficulties and future dir
 \\
 
 %==========================================
-\subsection{The computational challenge} 
+\paragraph{The computational challenge.} 
 %==========================================
 As emphasized above, the BSE eigenvalue equation in the single-excitation space [see Eq.~\eqref{eq:BSE-eigen}] is formally equivalent to that of TD-DFT or TD-HF. \cite{Dreuw_2005}
 Searching iteratively for the lowest eigenstates exhibits the same $\order*{\Norb^4}$ matrix-vector multiplication computational cost within BSE and TD-DFT. 
@@ -600,38 +604,36 @@ These ongoing developments pave the way to applying the $GW$@BSE formalism to sy
 \\
 
 %==========================================
-\subsection{The triplet instability challenge}
+\paragraph{The triplet instability challenge.}
 %==========================================
 The analysis of the singlet-triplet splitting is central to numerous applications such as singlet fission or thermally activated delayed fluorescence (TADF).
-From a more theoretical point of view, triplet instabilities that often plagues the applicability of TD-DFT are intimately linked to the stability analysis of restricted closed-shell solutions at the HF \cite{Seeger_1977} and KS \cite{Bauernschmitt_1996} levels.
+From a more theoretical point of view, triplet instabilities that often \titou{plague} the applicability of TD-DFT are intimately linked to the stability analysis of restricted closed-shell solutions at the HF \cite{Seeger_1977} and KS \cite{Bauernschmitt_1996} levels.
 While TD-DFT with range-separated hybrids can benefit from tuning the range-separation parameter(s) as a mean to act on the triplet instability, \cite{Sears_2011} BSE calculations do not offer this pragmatic way-out since the screened Coulomb potential that builds the kernel does not offer any parameter to tune.
 
 Benchmark calculations \cite{Jacquemin_2017b,Rangel_2017} clearly concluded that triplets are notably too low in energy within BSE and that the use of the TDA was able to partly reduce this error. 
 However, as it stands, the BSE accuracy for triplets remains rather unsatisfactory for reliable applications. 
-An alternative cure was offered by hybridizing TD-DFT and BSE, that is, by 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}
+An alternative cure was offered by hybridizing TD-DFT and BSE, that is, by 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_2018a}
 \\
 
 
 %==========================================
-\subsection{The challenge of the ground-state energy}
+\paragraph{The challenge of the ground-state energy.}
 %==========================================
-
 In contrast to TD-DFT which relies on KS-DFT as its ground-state analog, the ground-state BSE energy is not a well-defined quantity, and no clear consensus has been found regarding its formal definition.
-Consequently, the BSE ground-state formalism remains in its infancy with very few available studies for atomic and molecular systems. \cite{Olsen_2014,Holzer_2018,Li_2019,Li_2020,Loos_2020}
+Consequently, the BSE ground-state formalism remains in its infancy with very few available studies for atomic and molecular systems. \cite{Olsen_2014,Holzer_2018b,Li_2019,Li_2020,Loos_2020}
 
 A promising route, which closely follows RPA-type formalisms, \cite{Angyan_2011} is to calculate the ground-state BSE energy within the adiabatic-connection fluctuation-dissipation theorem (ACFDT) framework. \cite{Furche_2005}
 Thanks to comparisons with both similar and state-of-art computational approaches, it was recently shown that the ACFDT@BSE@$GW$ approach yields extremely accurate PES around equilibrium, and can even compete with high-order coupled cluster methods in terms of absolute ground-state energies and equilibrium distances. \cite{Loos_2020}
-However, their accuracy near the dissociation limit remains an open question. \cite{Caruso_2013,Olsen_2014,Colonna_2014,Hellgren_2015,Holzer_2018}
+However, their accuracy near the dissociation limit remains an open question. \cite{Caruso_2013,Olsen_2014,Colonna_2014,Hellgren_2015,Holzer_2018b}
-Indeed, in the largest available benchmark study \cite{Holzer_2018} encompassing the total energies of the atoms \ce{H}--\ce{Ne}, the atomization energies of the 26 small molecules forming the HEAT test set, and the bond lengths and harmonic vibrational frequencies of $3d$ transition-metal monoxides, the BSE correlation energy, as evaluated within the ACFDT framework, \cite{Furche_2005} was mostly discarded from the set of tested techniques due to instabilities (negative frequency modes in the BSE polarization propagator) and replaced by an approximate (RPAsX) approach where the screened-Coulomb potential matrix elements was removed from the resonant electron-hole contribution. \cite{Maggio_2016,Holzer_2018} 
+Indeed, in the largest available benchmark study \cite{Holzer_2018b} encompassing the total energies of the atoms \ce{H}--\ce{Ne}, the atomization energies of the 26 small molecules forming the HEAT test set, and the bond lengths and harmonic vibrational frequencies of $3d$ transition-metal monoxides, the BSE correlation energy, as evaluated within the ACFDT framework, \cite{Furche_2005} was mostly discarded from the set of tested techniques due to instabilities (negative frequency modes in the BSE polarization propagator) and replaced by an approximate (RPAsX) approach where the screened-Coulomb potential matrix elements was removed from the resonant electron-hole contribution. \cite{Maggio_2016,Holzer_2018} 
 Moreover, it was also observed in Ref.~\citenum{Loos_2020} that, in some cases, unphysical irregularities on the ground-state PES appear  due to the appearance of discontinuities as a function of the bond length for some of the $GW$ quasiparticle energies. 
 Such an unphysical behavior stems from defining the quasiparticle energy as the solution of the quasiparticle equation with the largest spectral weight in cases where several solutions can be found [see Eq.~\eqref{eq:QP-eq}].
 We refer the interested reader to Refs.~\citenum{vanSetten_2015,Maggio_2017,Loos_2018,Veril_2018,Duchemin_2020} for detailed discussions.
 \\
 
 %==========================================
-\subsection{The challenge of analytical nuclear gradients}
+\paragraph{The challenge of analytical nuclear gradients.}
 %==========================================
-
 The features of ground- and excited-state potential energy surfaces (PES) are critical for the faithful description and a deeper understanding of photochemical and photophysical processes. \cite{Olivucci_2010}
 For example, chemoluminescence and fluorescence are associated with geometric relaxation of excited states, and structural changes upon electronic excitation. \cite{Navizet_2011}
 Reliable predictions of these mechanisms, which have attracted much experimental and theoretical interest lately, require exploring the ground- and excited-state PES. 
@@ -648,12 +650,12 @@ In this seminal work devoted to small molecules (\ce{CO} and \ce{NH3}), only the
 \\
 
 %==========================================
-\subsection{Beyond the static approximation}
+\paragraph{Beyond the static approximation.}
 %==========================================
 Going beyond the static approximation is a difficult challenge which has been, nonetheless, embraced by several groups.\cite{Strinati_1988,Rohlfing_2000,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Huix-Rotllant_2011,Zhang_2013,Rebolini_2016,Olevano_2019}
 As mentioned earlier in this \textit{Perspective}, most of BSE calculations are performed within the so-called static approximation, which substitutes the dynamically-screened (\ie, frequency-dependent) Coulomb potential $W(\omega)$ by its static limit $W(\omega = 0)$ [see Eq.~\eqref{eq:Wmatel}].
 It is important to mention that diagonalizing the BSE Hamiltonian in the static approximation corresponds to solving a \textit{linear} eigenvalue problem in the space of single excitations, while it is, in its dynamical form, a non-linear eigenvalue problem (in the same space) which is much harder to solve from a numerical point of view.
-In complete analogy with the ubiquitous adiabatic approximation in TD-DFT, one key consequence of the static approximation is that double (and higher) excitations are completely absent from the BSE optical spectrum, which obviously hampers the applicability of BSE as double excitation may play, indirectly, a key role in photochemistry mechanisms. 
+In complete analogy with the ubiquitous adiabatic approximation in TD-DFT, one key consequence of the static approximation is that double (and higher) excitations are completely absent from the BSE optical spectrum, which obviously hampers the applicability of BSE as \titou{double excitations} may play, indirectly, a key role in photochemistry mechanisms. 
 Higher excitations would be explicitly present in the BSE Hamiltonian by ``unfolding'' the dynamical BSE kernel, and one would recover a linear eigenvalue problem with, nonetheless, a much larger dimension.
 Corrections to take into account the dynamical nature of the screening may or may not recover these multiple excitations.
 However, dynamical corrections permit, in any case, to recover, for transitions with a dominant single-excitation character, additional relaxation effects coming from higher excitations.
@@ -667,77 +669,8 @@ Rebolini and Toulouse have performed a similar investigation in a range-separate
 In these two latter studies, they also followed a (non-self-consistent) perturbative approach within TDA with a renormalization of the first-order perturbative correction.
 \\
 
-%==========================================
-%\subsection{The double excitation challenge}
-%==========================================
-%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_2009b,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}
-%	\label{eq:lin_sys}
-%	\bA \bx = \omega \bx,
-%\end{equation}
-%where $\omega$ is one of the neutral excitation energies of the system associated with the transition vector $\bx$.
-%If we assume that the operator $\bA$ has a matrix representation of size $K \times K$, this \textit{linear} set of equations yields $K$ excitation energies.
-%However, in practice, $K$ might be very large, and it might therefore be practically useful to recast this system as two smaller coupled systems, such that 
-%\begin{equation}
-%	\label{eq:lin_sys_split}
-%	\begin{pmatrix}
-%		\bA_1	&	\tr{\bb}		\\
-%		\bb			&	\bA_2		\\
-%	\end{pmatrix}
-%	\begin{pmatrix}
-%		\bx_1	\\
-%		\bx_2	\\
-%	\end{pmatrix}
-%	= \omega
-%	\begin{pmatrix}
-%		\bx_1	\\
-%		\bx_2  	\\
-%	\end{pmatrix},
-%\end{equation}
-%where the blocks $\bA_1$ and $\bA_2$, of sizes $K_1 \times K_1$ and $K_2 \times K_2$ (with $K_1 + K_2 =  K$), can be associated with, for example, the single and double excitations of the system.
-%Note that this \textit{exact} decomposition does not alter, in any case, the values of the excitation energies, not their eigenvectors.
-%
-%Solving separately each row of the system \eqref{eq:lin_sys_split} yields
-%\begin{subequations}
-%\begin{gather}
-%	\label{eq:row1}
-%	\bA_1 \bx_1  + \tr{\bb} \bx_2 = \omega \bx_1,
-%	\\
-%	\label{eq:row2}
-%	\bx_2 = (\omega \bI - \bA_2)^{-1} \bb \bx_1.
-%\end{gather}
-%\end{subequations}
-%Substituting Eq.~\eqref{eq:row2} into Eq.~\eqref{eq:row1} yields the following effective \textit{non-linear}, frequency-dependent operator
-%\begin{equation}
-%	\label{eq:non_lin_sys}
-%	\Tilde{\bA}_1(\omega) \bx_1 = \omega \bx_1
-%\end{equation}
-%with 
-%\begin{equation}
-%	\Tilde{\bA}_1(\omega) = \bA_1 +  \tr{\bb} (\omega \bI - \bA_2)^{-1} \bb 
-%\end{equation}
-%which has, by construction, exactly the same solutions than the linear system \eqref{eq:lin_sys} but a smaller dimension.
-%For example, an operator $\Tilde{\bA}_1(\omega)$ built in the basis of single excitations can potentially provide excitation energies for double excitations thanks to its frequency-dependent nature, the information from the double excitations being ``folded'' into $\Tilde{\bA}_1(\omega)$ via Eq.~\eqref{eq:row2}. \cite{Romaniello_2009b,Sangalli_2011,ReiningBook}
-%
-%How have we been able to reduce the dimension of the problem while keeping the same number of solutions?
-%To do so, we have transformed a linear operator $\bA$ into a non-linear operator $\Tilde{\bA}_1(\omega)$ by making it frequency dependent.
-%In other words, we have sacrificed the linearity of the system in order to obtain a new, non-linear systems of equations of smaller dimension. 
-%This procedure converting degrees of freedom into frequency or energy dependence is very  general and can be applied in various contexts. \cite{Gatti_2007,Garniron_2018}
-%Thanks to its non-linearity, Eq.~\eqref{eq:non_lin_sys} can produce more solutions than its actual dimension. 
-%However, because there is no free lunch, this non-linear system is obviously  harder to solve than its corresponding linear analogue given by Eq.~\eqref{eq:lin_sys}.
-%Nonetheless, approximations can be now applied to Eq.~\eqref{eq:non_lin_sys} in order to solve it efficiently.
-%
-%One of these approximations is the so-called \textit{static} approximation, which corresponds to fix the frequency to a particular value. 
-%For example, as commonly done within the Bethe-Salpeter formalism, $\Tilde{\bA}_1(\omega) = \Tilde{\bA}_1 \equiv \Tilde{\bA}_1(\omega = 0)$.
-%In such a way, the operator $\Tilde{\bA}_1$ is made linear again by removing its frequency-dependent nature.
-%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.\\
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Conclusions}
+\paragraph{Conclusion.}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Although far from being exhaustive, we hope that this \textit{Perspective} provides a concise and fair assessment of the strengths and weaknesses of the BSE formalism of many-body perturbation theory.
 To do so, we have briefly reviewed the theoretical aspects behind BSE, and its intimate link with the underlying $GW$ calculation that one must perform to compute quasiparticle energies and the dynamically-screened Coulomb potential; two of the key input ingredients associated with the BSE formalism.