The many-body Green's function Bethe-Salpeter equation (BSE) formalism has shown to be a promising alternative to time-dependent density-functional theory (TD-DFT) in order to compute vertical transition energies for medium and large molecular systems.
In particular, BSE@\textit{GW} faithfully treats charge-transfer states and avoids the delicate choice of the exchange-correlation functional from an ever-growing zoo of functionals.
However, contrary to TD-DFT, there is no clear definition of the ground- and excited-state energies within BSE, which makes the development of analytical nuclear gradients a particularly tricky task.
Here, we provide an unambiguous definition of the ground-state BSE energy, and we calculate the excited-state BSE energy of a given state by adding the corresponding BSE excitation energy to the ground-state BSE energy.
Embracing this definition which treats on equal footing ground and excited states at the BSE level, we study the topological features of the ground- and (singlet and triplet) excited-state potential energy surfaces (PES) for several diatomic molecules.
Our aim is to know whether or not the BSE formalism is able to reproduce faithfully the main features of these PES near equilibrium, and, in particular, the location of the minima on the ground- and excited-state PES.
\titou{Thanks to comparison with both similar and state-of-art computational approaches, we show that ...}
\end{abstract}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%
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{Bernardi_1996,Olivucci_2010,Robb_2007}
For example, chemoluminescence, fluorescence and other related processes 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.
From a theoretical point of view, the accurate prediction of excited electronic states remains a challenge, \cite{Gonzales_2012, Loos_2020a} especially for large systems where state-of-the-art computational techniques (such as multiconfigurational methods \cite{Andersson_1990,Andersson_1992,Roos_1996,Angeli_2001}) cannot be afforded.
For such systems, one has to rely on more approximate, yet computationally cheaper approaches. \cite{Grimme_2004a,Ghosh_2018}
For the last two decades, time-dependent density-functional theory (TD-DFT) \cite{Casida} has been the go-to method to compute absorption and emission spectra in large molecular systems.
At a relatively low computational cost, TD-DFT can provide accurate vertical and adiabatic transition energies for low-lying excited states of organic molecules with a typical error of $0.2$--$0.4$ eV. \cite{Loos_2019b}
Similar to density-functional theory (DFT), \cite{Hohenberg_1964,Kohn_1965,ParrBook} TD-DFT is an in-principle exact theory which formal foundation relies on the Runge-Gross theorem. \cite{Runge_1984}
The Kohn-Sham (KS) formalism of TD-DFT transfers the complexity of the many-body problem to the exchange-correlation (xc) functional thanks to a judicious mapping between a time-dependent non-interacting reference system and its interacting analog which have both the exact same electronic density.
In TD-DFT, the excitation energies are obtained as the poles of the ground-state frequency-dependent linear response function, which is obtained by solving a Dyson equation in which the exchange and correlation effects are cast in the xc kernel.
This xc kernel plays for the excited states the same role as the xc functional for the ground state.
Hence, the PES for the excited states can be easily and efficiently obtained as a function of the molecular geometry by simply adding the ground-state DFT energy to the excitation energy of the selected state.
One of the strongest assets of TD-DFT is the availability of first- and second-order analytic nuclear gradients (\ie, the first derivatives of the excited-state energy with respect to the nuclear displacements), which enables the exploration of excited-state PES.\cite{Furche_2002}
Although the exact spectrum of the interacting system can be obtained in principle, both the ground-state xc potential and the xc kernel (which enters the definition of the response function) have to be approximated in practice.
This has major consequences, and it is well documented that, TD-DFT fails badly to model, for examples, charge transfer,\cite{Tozer_1999,Dreuw_2003,Sobolewski_2003,Dreuw_2004,Maitra_2017} Rydberg \cite{Tozer_1998,Tozer_2000,Casida_1998,Casida_2000,Tozer_2003} and doubly excited states. \cite{Levine_2006,Tozer_2000,Elliott_2011}
An important point from a practical point of view is that the TD-DFT xc kernel is usually considered as static instead of being frequency dependent. \cite{Romaniello_2009a,Sangalli_2011}
One key consequence of this so-called adiabatic approximation (based on the assumption that the density varies slowly with time) is that double excitations are completely absent from the TD-DFT spectra. \cite{Levine_2006,Tozer_2000,Huix-Rotllant_2010,Elliott_2011}
Yet another problem is the choice of the xc functionals as the quality of excitation energies are substantially dependent on this choice.
With a similar computational cost, the many-body Green's function Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} is a very valuable alternative to TD-DFT with early \textit{ab initio} calculations in condensed matter physics appearing at the end of the 90's. \cite{Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999}
In the past few years, BSE has gained momentum for the study of molecular systems \cite{Ma_2009,Pushchnig_2002,Tiago_2003,Tiago_2008,Sai_2008,Palumno_2009,Rocca_2010,Sharifzadeh_2012,Cudazzo_2012,Boulanger_2014,Ljungberg_2015,Hirose_2015,Cocchi_2015,Ziaei_2017,Abramson_2017} and is now a serious candidate as a computationally inexpensive method that can effectively model excited states with a typical error of $0.1$--$0.3$ eV according to large and systematic benchmark calculations. \cite{Jacquemin_2015,Bruneval_2015,Blase_2016,Jacquemin_2016,Hung_2016,Hung_2017,Krause_2017,Jacquemin_2017,Blase_2018}
One of the main advantage of BSE compared to TD-DFT is that it allows a faithful description of charge-transfer states. \cite{Lastra_2011,Blase_2011b,Baumeier_2012,Duchemin_2012,Cudazzo_2013,Ziaei_2016}
Moreover, when performed on top of a (partially) self-consistently {\evGW} calculation, \cite{Hybertsen_1986, Shishkin_2007, Blase_2011, Faber_2011,Rangel16,Kaplan_2018,Gui_2018} BSE@{\evGW} has been shown to be weakly dependent on its starting point (\ie, on the xc functional selected for the underlying DFT calculation). \cite{Jacquemin_2016,Gui_2018}
However, similar to TD-DFT, the static version of BSE cannot describe multiple excitations. \cite{Romaniello_2009a,Sangalli_2011}
A significant limitation of the BSE formalism as compared to TD-DFT lies in the lack of analytic forces for the excited states, preventing efficient applications to the study of photoluminescence, photocatalysis, etc. While calculations of the {\GW} quasiparticle energies ionic gradients is becoming very popular,
\cite{Lazzeri_2008,Faber_2011b,Yin_2013,Faber_2015,Montserrat_2016,Zhenglu_2019} only one pioneering study of the excited-state BSE gradients has been published so far. \cite{Beigi_2003} In this study devoted to small molecules (\ce{CO} and \ce{NH3}), only the BSE excitation energy gradients were calculated, while computing the KS-DFT (LDA) forces as ground-state analogue.
Contrary to TD-DFT, the ground-state correlation energy calculated at the BSE level remains in its infancy with very few available studies for atomic and molecular systems. \cite{Olsen_2014,Holzer_2018,Li_2019,Li_2020}
As a matter of fact, in the largest recent available benchmark study of the 26 small molecules forming the HEAT test set, \cite{Holzer_2018} the BSE correlation energy, as evaluated within the adiabatic connection formulation (AC-BSE), 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}
Such a modified BSE polarization propagator was inspired by a previous study on the homogeneous electron gas. \cite{Maggio_2016}
With such an approximation, amounting to neglect excitonic effects in the electron-hole propagator, the question of using either KS-DFT or {\GW} eigenvalues in the construction of the propagator becomes further relevant, increasing accordingly the number of possible definitions for the ground-state correlation energy.
Finally, renormalizing or not the Coulomb interaction by the coupling parameter $\lambda$ in the Dyson equation for the interacting polarizability leads to two different versions of the AC-BSE correlation energy.
Here, in analogy to the random-phase approximation (RPA) formalism, \cite{Furche_2008} the ground-state BSE energy is calculated via the ``trace'' formula (see below). The excited-state BSE energy is then computed by adding the BSE excitation energy of the selected state to the ground-state BSE energy.
This definition of the energy has the advantage of treating at the same level of theory the ground state and the excited states.
Embracing this definition, the purpose of the present study is to investigate the quality of ground- and excited-state PES near equilibrium obtained within the BSE approach.
The location of the minima on the ground- and (singlet and triplet) excited-state PES is of particular interest.
In Sec.~\ref{sec:theo}, we introduce the equations behind the BSE formalism.
In particular, the general equations are reported in Subsec.~\ref{sec:BSE}, and the corresponding equations obtained in a finite basis in Subsec.~\ref{sec:BSE_basis}.
Subsection \ref{sec:BSE_energy} defines the BSE total energy, and various other quantities of interest for this study.
In Sec.~\ref{sec:comp_details}, computational details are reported.
Section \ref{sec:PES} reports PES of the ground- and excited-states for various diatomic molecules.
Finally, we draw our conclusion in Sec.~\ref{sec:conclusion}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theo}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The Bethe-Salpeter equation}
\label{sec:BSE}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to compute the neutral (optical) excitations of the system and their associated oscillator strengths, the BSE expresses the two-body propagator \cite{Strinati_1988, Martin_2016}
as the linear response of the one-body Green's function $\G{}$ with respect to a general non-local external potential
\begin{equation}
\XiBSE{}(3,5,4,6) = i \fdv{[\vc{\Ha}(3) \delta(3,4) + \Sig{\xc}(3,4)]}{\G{}(6,5)},
\end{equation}
which takes into account the self-consistent variation of the Hartree potential
\begin{equation}
\vc{\Ha}(1) = - i \int d2 \vc{}(2) \G{}(2,2^+),
\end{equation}
(where $\vc{}$ is the bare Coulomb operator) and the exchange-correlation self-energy $\Sig{\xc}$.
In Eq.~\eqref{eq:BSE}, $\LBSE{0}(1,2,1',2')=- i \G{}(1,2')\G{}(2,1')$, and $(1)=(\br{1},\sigma_1,t_1)$ is a composite index gathering space, spin and time variables.
In the {\GW} approximation, \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Martin_2016,Reining_2017} we have
\begin{equation}
\SigGW{\xc}(1,2) = i \G{}(1,2) \W{}(1^+,2),
\end{equation}
where $\W{}$ is the screened Coulomb operator, and hence the BSE reduces to
where, as commonly done, we have neglected the term $\delta\W{}/\delta\G{}$ in the functional derivative of the self-energy. \cite{Hanke_1980,Strinati_1984,Strinati_1982}
Finally, the static approximation is enforced, \ie, $\W{}(1,2)=\W{}(\{\br{1}, \sigma_1, t_1\},\{\br{2}, \sigma_2, t_2\})\delta(t_1-t_2)$, which corresponds to restricting $\W{}$ to its static limit, \ie, $\W{}(1,2)=\W{}(\{\br{1},\sigma_1\},\{\br{2},\sigma_2\}; \omega=0)$.
To compute the singlet and triplet BSE excitation energies in a finite basis within the static approximation, one must solve the following linear response problem \cite{Casida,Dreuw_2005,Martin_2016}
where $\OmBSE{m}$ is the $m$th excitation energy with eigenvector $\T{(\bX_m \,\bY_m)}$, and we have assumed real-valued orbitals $\{\MO{p}(\br{})\}_{1\le p \le\Norb}$.
The matrices $\bA$, $\bB$, $\bX$, and $\bY$ are all of size $\Nocc\Nvir\times\Nocc\Nvir$ where $\Nocc$ and $\Nvir$ are the number of occupied and virtual orbitals (\ie, $\Norb=\Nocc+\Nvir$), respectively.
In the following, the index $m$ labels the $\Nocc\Nvir$ single excitations, $i$, $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals.
In the absence of triplet instabilities, \cite{Dreuw_2005} Eq.~\eqref{eq:LR} is usually transformed into an Hermitian eigenvalue problem of smaller dimension
In Eq.~\eqref{eq:W}, $\OmRPA{m}$ are the neutral (direct) RPA excitation energies computed during the {\GW} calculation, and $\eta$ is a positive infimitesimal.
is the ground-state BSE correlation energy computed with the so-called trace formula, \cite{Schuck_Book, Rowe_1968, Sawada_1957b} and $\OmBSE{m}$ is the $m$th BSE excitation energy with the convention that, for the ground state ($m=0$), $\OmBSE{0}=0$.
Equation \eqref{eq:EtotBSE} defines unambiguously the total BSE energy of the system for both ground and (singlet and triplet) excited states.
Note that the prime in the summation of Eq.~\eqref{eq:EcBSE} means that we only consider the \textit{positive} excitation energies from the singlet and triplet manifolds.
Indeed, in case of triplet instabilities, some of the triplet excitation energies are negative, and must be discarded as the resonant-only part of the BSE excitonic Hamiltonian has to be considered.
Note also that an alternative formulation does exist which consists in integrating along the adiabatic connection path. \cite{Gell-Mann_1957}
These two formulation have been found to be equivalent in practice for both the uniform electron gas \cite{Sawada_1957b, Fukuta_1964, Furche_2008} and in molecules. \cite{Li_2020}\xavier{WELL, THIS IS UNCLEAR TO ME : this is true for direct RPA, but is this true for BSE ... I am really not sure ...}
From a practical point of view, it is also convenient to split $\EcBSE=\EcsBSE+\EctBSE$ in its singlet ($\sigma\neq\sigma^{\prime}$) and triplet ($\sigma=\sigma^{\prime}$) contributions, \ie,
As a final remark, we point out that Eq.~\eqref{eq:EtotBSE} can be easily generalized to other theories (such as CIS, RPA, or TDHF) by computing $\Ec$ and $\Om{m}$ at the different level of theory.
All the preliminary {\GW} calculations performed to obtain the screened Coulomb operator and the quasiparticle energies have been done using a Hartree-Fock (HF) starting point, which is a very adequate choice in the case of the (small) systems that we consider here.
Both perturbative {\GW} (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986} and partially self-consistent {\evGW}\cite{Hybertsen_1986, Shishkin_2007, Blase_2011, Faber_2011} calculations are performed here.
These will be labeled as BSE@{\GOWO} and BSE@{\evGW}, respectively.
In the case of {\GOWO}, the quasiparticle energies have been obtained by linearizing the non-linear, frequency-dependent quasiparticle equation.
For {\evGW}, the quasiparticle energies are obtained self-consistently and we have used the DIIS convergence accelerator technique proposed by Pulay \cite{Pulay_1980,Pulay_1982} to avoid convergence issues.
Further details about our implementation of {\GOWO} and {\evGW} can be found in Refs.~\onlinecite{Loos_2018,Veril_2018}.
Finally, the infinitesimal $\eta$ has been set to $10^{-3}$ for all calculations.
Because Eq.~\eqref{eq:EcBSE} requires the entire BSE excitation spectrum (both singlet and triplet), we perform a complete diagonalization of the $OV \times OV$ BSE linear response matrix [see Eq.~\eqref{eq:small-LR}], which corresponds to a $\order{O^3V^3}$ computational cost.
This step is, by far, the computational bottleneck in our current implementation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Potential energy surfaces}
\label{sec:PES}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hydrogen molecule}
\label{sec:H2}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Lithium hydride and lithium fluoride}
\label{sec:LiX}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The isoelectronic sequence: \ce{N2}, \ce{CO}, and \ce{BF}}
\label{sec:isoN2}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%
\titou{A nice conclusion here saying that what we have done is pretty awesome.}
%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supporting Information}
%%%%%%%%%%%%%%%%%%%%%%%%
See {\SI} for \titou{bla bla bla.}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
This work was performed using HPC resources from GENCI-TGCC (Grant No.~2018-A0040801738) and CALMIP (Toulouse) under allocation 2019-18005.
Funding from the \textit{``Centre National de la Recherche Scientifique''} is acknowledged.
This work has been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''}.
