\documentclass[aps,prb,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
\usepackage[version=4]{mhchem}

\usepackage{natbib}
\usepackage[extra]{tipa}
\bibliographystyle{achemso}
\AtBeginDocument{\nocite{achemso-control}}


\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}

\usepackage[
	colorlinks=true,
    citecolor=blue,
    breaklinks=true
	]{hyperref}
\urlstyle{same}

\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{HTML}{009900}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\denis}[1]{\textcolor{purple}{#1}}
\newcommand{\xavier}[1]{\textcolor{darkgreen}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashDJ}[1]{\textcolor{purple}{\sout{#1}}}
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\renewcommand{\DJ}[1]{\denis{(\underline{\bf DJ}: #1)}}
\newcommand{\XB}[1]{\xavier{(\underline{\bf XB}: #1)}}

\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\SI}{\textcolor{blue}{supporting information}}
\newcommand{\QP}{\textsc{quantum package}}
\newcommand{\T}[1]{#1^{\intercal}}

% coordinates
\newcommand{\br}[1]{\mathbf{r}_{#1}}
\newcommand{\dbr}[1]{d\br{#1}}

% methods
\newcommand{\evGW}{ev$GW$}	
\newcommand{\qsGW}{qs$GW$}	
\newcommand{\GOWO}{$G_0W_0$}	
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\xc}{\text{xc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\co}{\text{x}}

% 
\newcommand{\Norb}{N}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}
\newcommand{\IS}{\lambda}

% operators
\newcommand{\hH}{\Hat{H}}

% energies
\newcommand{\Enuc}{E^\text{nuc}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\EHF}{E^\text{HF}}
\newcommand{\EBSE}{E^\text{BSE}}
\newcommand{\EcRPA}{E_\text{c}^\text{RPA}}
\newcommand{\EcRPAx}{E_\text{c}^\text{RPAx}}
\newcommand{\EcBSE}{E_\text{c}^\text{BSE}}
\newcommand{\IP}{\text{IP}}
\newcommand{\EA}{\text{EA}}

% orbital energies
\newcommand{\e}[1]{\epsilon_{#1}}
\newcommand{\eHF}[1]{\epsilon^\text{HF}_{#1}}
\newcommand{\eKS}[1]{\epsilon^\text{KS}_{#1}}
\newcommand{\eQP}[1]{\epsilon^\text{QP}_{#1}}
\newcommand{\eGOWO}[1]{\epsilon^\text{\GOWO}_{#1}}
\newcommand{\eGW}[1]{\epsilon^{GW}_{#1}}
\newcommand{\eevGW}[1]{\epsilon^\text{\evGW}_{#1}}
\newcommand{\eGnWn}[2]{\epsilon^\text{\GnWn{#2}}_{#1}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}

% Matrix elements
\newcommand{\A}[2]{A_{#1}^{#2}}
\newcommand{\tA}[2]{\Tilde{A}_{#1}^{#2}}
\newcommand{\B}[2]{B_{#1}^{#2}}
\renewcommand{\S}[1]{S_{#1}}
\newcommand{\ABSE}[2]{A_{#1}^{#2,\text{BSE}}}
\newcommand{\BBSE}[2]{B_{#1}^{#2,\text{BSE}}}
\newcommand{\ARPA}[2]{A_{#1}^{#2,\text{RPA}}}
\newcommand{\BRPA}[2]{B_{#1}^{#2,\text{RPA}}}
\newcommand{\ARPAx}[2]{A_{#1}^{#2,\text{RPAx}}}
\newcommand{\BRPAx}[2]{B_{#1}^{#2,\text{RPAx}}}
\newcommand{\G}[1]{G_{#1}}
\newcommand{\LBSE}[1]{L_{#1}}
\newcommand{\XiBSE}[1]{\Xi_{#1}}
\newcommand{\Po}[1]{P_{#1}}
\newcommand{\W}[2]{W_{#1}^{#2}}
\newcommand{\Wc}[1]{W^\text{c}_{#1}}
\newcommand{\vc}[1]{v_{#1}}
\newcommand{\Sig}[1]{\Sigma_{#1}}
\newcommand{\SigGW}[1]{\Sigma^{GW}_{#1}}
\newcommand{\Z}[1]{Z_{#1}}
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\sERI}[2]{[#1|#2]}

%% bold in Table
\newcommand{\bb}[1]{\textbf{#1}}
\newcommand{\rb}[1]{\textbf{\textcolor{red}{#1}}}
\newcommand{\gb}[1]{\textbf{\textcolor{darkgreen}{#1}}}

% excitation energies
\newcommand{\OmRPA}[2]{\Omega_{#1}^{#2,\text{RPA}}}
\newcommand{\OmRPAx}[2]{\Omega_{#1}^{#2,\text{RPAx}}}
\newcommand{\OmBSE}[2]{\Omega_{#1}^{#2,\text{BSE}}}

\newcommand{\spinup}{\downarrow}
\newcommand{\spindw}{\uparrow}
\newcommand{\singlet}{\uparrow\downarrow}
\newcommand{\triplet}{\uparrow\uparrow}

% Matrices
\newcommand{\bO}{\mathbf{0}}
\newcommand{\bI}{\mathbf{1}}
\newcommand{\bvc}{\mathbf{v}}
\newcommand{\bSig}{\mathbf{\Sigma}}
\newcommand{\bSigX}{\mathbf{\Sigma}^\text{x}}
\newcommand{\bSigC}{\mathbf{\Sigma}^\text{c}}
\newcommand{\bSigGW}{\mathbf{\Sigma}^{GW}}
\newcommand{\be}{\mathbf{\epsilon}}
\newcommand{\beGW}{\mathbf{\epsilon}^{GW}}
\newcommand{\beGnWn}[1]{\mathbf{\epsilon}^\text{\GnWn{#1}}}
\newcommand{\bde}{\mathbf{\Delta\epsilon}}
\newcommand{\bdeHF}{\mathbf{\Delta\epsilon}^\text{HF}}
\newcommand{\bdeGW}{\mathbf{\Delta\epsilon}^{GW}}
\newcommand{\bOm}[1]{\mathbf{\Omega}^{#1}}
\newcommand{\bA}[1]{\mathbf{A}^{#1}}
\newcommand{\btA}[1]{\Tilde{\mathbf{A}}^{#1}}
\newcommand{\bB}[1]{\mathbf{B}^{#1}}
\newcommand{\bX}[1]{\mathbf{X}^{#1}}
\newcommand{\bY}[1]{\mathbf{Y}^{#1}}
\newcommand{\bZ}[1]{\mathbf{Z}^{#1}}
\newcommand{\bK}{\mathbf{K}}
\newcommand{\bP}[1]{\mathbf{P}^{#1}}

% units
\newcommand{\IneV}[1]{#1 eV}
\newcommand{\InAU}[1]{#1 a.u.}
\newcommand{\InAA}[1]{#1 \AA}
\newcommand{\kcal}{kcal/mol}

\newcommand{\NEEL}{Universit\'e Grenoble Alpes, CNRS, Institut NEEL, F-38042 Grenoble, France}
\newcommand{\CEISAM}{Laboratoire CEISAM - UMR CNRS 6230, Universit\'e de Nantes, 2 Rue de la Houssini\`ere, BP 92208, 44322 Nantes Cedex 3, France}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\CEA}{Universit\'e Grenoble Alpes, CEA, IRIG-MEM-L Sim, 38054 Grenoble, France}

\begin{document}	

\title{Ground-State Potential Energy Surfaces Within the Bethe-Salpeter Formalism: Pros and Cons}

\author{Pierre-Fran\c{c}ois \surname{Loos}}
	\email{loos@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}
\author{Anthony \surname{Scemama}}
	\email{scemama@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}
\author{Ivan \surname{Duchemin}}
	\email{ivan.duchemin@cea.fr}
	\affiliation{\CEA}
\author{Denis \surname{Jacquemin}}
	\email{denis.jacquemin@univ-nantes.fr}
	\affiliation{\CEISAM}
\author{Xavier \surname{Blase}}
	\email{xavier.blase@neel.cnrs.fr }
	\affiliation{\NEEL}

\begin{abstract}
%\begin{wrapfigure}[12]{o}[-1.25cm]{0.4\linewidth}
%	\centering
%  \includegraphics[width=\linewidth]{TOC}
%\end{wrapfigure}
The combined many-body Green's function $GW$ and 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 of molecular systems. 
The BSE formalism can also be employed to compute ground-state correlation energies thanks to the adiabatic-connection fluctuation-dissipation theorem (ACFDT).
Here, we study the topological features of the ground-state potential energy surfaces (PES) of several diatomic molecules near their equilibrium bond length. 
Thanks to comparisons with state-of-art computational approaches, we show that ACFDT@BSE is surprisingly accurate, and can even compete with coupled cluster methods in terms of total energies and equilibrium bond distances.  
However, we sometimes observe unphysical irregularities on the ground-state PES in relation with the appearance of satellite resonances with a weight similar to that of the $GW$ quasiparticle peak.
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Introduction}
%\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%

With a similar computational cost to time-dependent density-functional theory (TD-DFT), \cite{Runge_1984,Casida} the many-body Green's function Bethe-Salpeter equation (BSE) formalism 
\cite{Salpeter_1951,Strinati_1988,Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999} is a valuable alternative that has gained momentum in the past few years for studying molecular systems.\cite{Ma_2009,Pushchnig_2002,Tiago_2003,Palumno_2009,Rocca_2010,Sharifzadeh_2012,Cudazzo_2012,Boulanger_2014,Ljungberg_2015,Hirose_2015,Cocchi_2015,Ziaei_2017,Abramson_2017} 
It now stands as  a computationally inexpensive method that can effectively model excited states \cite{Gonzales_2012,Loos_2020a} 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 advantages 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,Rangel_2016,Kaplan_2016,Gui_2018} BSE@{\evGW} has been shown to be weakly dependent on its starting point (\ie, on the exchange-correlation functional selected for the underlying DFT calculation). \cite{Jacquemin_2016,Gui_2018}
However, similar to adiabatic TD-DFT, \cite{Levine_2006,Tozer_2000,Huix-Rotllant_2010,Elliott_2011} the static version of BSE cannot describe multiple excitations. \cite{Romaniello_2009a,Sangalli_2011,Loos_2019}

A significant limitation of the BSE formalism, as compared to TD-DFT, lies in the lack of analytical nuclear gradients (\ie, the first derivatives of the energy with respect to the nuclear displacements) for both the ground and excited states, \cite{Furche_2002} preventing efficient applications to the study of chemoluminescence, fluorescence and other related processes \cite{Navizet_2011} associated with geometric relaxation of ground and excited states, and structural changes upon electronic excitation. \cite{Bernardi_1996,Olivucci_2010,Robb_2007}
While calculations of the $GW$ quasiparticle energy ionic gradients is becoming 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 Kohn-Sham (KS) LDA forces as its ground-state analog.

Contrary to TD-DFT which relies on KS-DFT \cite{Hohenberg_1964,Kohn_1965,ParrBook} as its ground-state counterpart, the BSE ground-state energy is not a well-defined quantity, and no clear consensus has been found regarding its formal definition.
It then 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 \cite{Holzer_2018} of the total energies of the atoms \ce{H}--\ce{Ne}, the atomization energies of the 26 small molecules forming the HEAT test set, \cite{Harding_2008} and the bond lengths and harmonic vibrational frequencies of $3d$ transition-metal monoxides, the BSE correlation energy, as evaluated within the adiabatic-connection fluctuation-dissipation theorem (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} 
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 interaction strength $\IS$ in the Dyson equation for the interacting polarizability (see below) leads to two different versions of the BSE correlation energy. \cite{Holzer_2018}

Here, in analogy to the random-phase approximation (RPA)-type formalisms \cite{Furche_2008,Toulouse_2009,Toulouse_2010,Angyan_2011,Ren_2012} and similarly to Refs.~\onlinecite{Olsen_2014,Maggio_2016,Holzer_2018}, the ground-state BSE energy is calculated in the adiabatic connection framework. 
Embracing this definition, the purpose of the present study is to investigate the quality of ground-state PES near equilibrium obtained within the BSE approach for several diatomic molecules.
The location of the minima on the ground-state PES is of particular interest.
This study is a first preliminary step towards the development of analytical nuclear gradients within the BSE@$GW$ formalism.
Thanks to comparisons with both similar and state-of-art computational approaches, we show that the ACFDT@BSE@$GW$ approach is surprisingly accurate, and can even compete with high-order coupled cluster (CC) methods in terms of absolute energies and equilibrium distances.
However, we also observe that, in some cases, unphysical irregularities on the ground-state PES, which are due to the appearance of a satellite resonance with a weight similar to that of the $GW$ quasiparticle peak. \cite{vanSetten_2015,Maggio_2017,Loos_2018,Veril_2018,Duchemin_2020}

%The paper is organized as follows. 
%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.
%Computational details needed to reproduce the results of the present work are reported in Sec.~\ref{sec:comp_details}.
%Section \ref{sec:PES} reports  ground-state PES 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} 
\begin{multline}
\label{eq:BSE}
	\LBSE{}(1,2,1',2') = \LBSE{0}(1,2,1',2') 
	\\
	+ \int d3 d4 d5  d6 \LBSE{0}(1,4,1',3) \XiBSE{}(3,5,4,6) \LBSE{}(6,2,5,2')
\end{multline}
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 
\begin{equation}
	\XiBSE{}(3,5,4,6) = \delta(3,4) \delta(5,6) \vc{}(3,6) - \delta(3,6) \delta(4,5) \W{}{}(3,4),
\end{equation}
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)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{BSE in a finite basis}
%\label{sec:BSE_basis}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

For a closed-shell system in a finite basis, to compute the singlet BSE excitation energies (within the static approximation) of the physical system (\ie, $\IS = 1$), one must solve the following linear response problem \cite{Casida,Dreuw_2005,Martin_2016} 
\begin{equation}
\label{eq:LR}
	\begin{pmatrix}
		\bA{\IS}	&	\bB{\IS}	\\
		-\bB{\IS}	&	-\bA{\IS}	\\
	\end{pmatrix}
	\begin{pmatrix}
		\bX{\IS}_m	\\
		\bY{\IS}_m	\\
	\end{pmatrix}
	=
	\Om{m}{\IS}
	\begin{pmatrix}
		\bX{\IS}_m	\\
		\bY{\IS}_m	\\
	\end{pmatrix},
\end{equation}
where $\Om{m}{\IS}$ is the $m$th excitation energy with eigenvector $\T{(\bX{\IS}_m \, \bY{\IS}_m)}$ at interaction strength $\IS$, $\T{}$ is the matrix transpose, and we assume real-valued spatial orbitals $\{\MO{p}(\br{})\}_{1 \le p \le \Norb}$.
The matrices $\bA{\IS}$, $\bB{\IS}$, $\bX{\IS}$, and $\bY{\IS}$ 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$ is the total number of spatial orbitals), respectively.
In the following, the index $m$ labels the $\Nocc \Nvir$ single excitations, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals.

In the absence of instabilities (\ie, $\bA{\IS} - \bB{\IS}$ is positive-definite), \cite{Dreuw_2005} Eq.~\eqref{eq:LR} is usually transformed into an Hermitian eigenvalue problem of smaller dimension
\begin{equation}
\label{eq:small-LR}
	(\bA{\IS} - \bB{\IS})^{1/2} (\bA{\IS} + \bB{\IS}) (\bA{\IS} - \bB{\IS})^{1/2} \bZ{\IS} = (\bOm{\IS})^2 \bZ{\IS},
\end{equation}
where the excitation amplitudes are
\begin{subequations}
\begin{align}
	\bX{\IS} + \bY{\IS} = (\bOm{\IS})^{-1/2} (\bA{\IS} - \bB{\IS})^{+1/2} \bZ{\IS},
	\\
	\bX{\IS} - \bY{\IS} = (\bOm{\IS})^{+1/2} (\bA{\IS} - \bB{\IS})^{-1/2} \bZ{\IS}.
\end{align}
\end{subequations}
Introducing the so-called Mulliken notation for the bare two-electron integrals
\begin{equation}
	\ERI{pq}{rs} = \iint \frac{\MO{p}(\br{}) \MO{q}(\br{}) \MO{r}(\br{}') \MO{s}(\br{}')}{\abs*{\br{} - \br{}'}} \dbr{} \dbr{}',
\end{equation}
and the corresponding (static) screened Coulomb potential matrix elements at coupling strength $\IS$ 
\begin{equation}
	\W{pq,rs}{\IS} = \iint  \MO{p}(\br{}) \MO{q}(\br{}) \W{}{\IS}(\br{},\br{}') \MO{r}(\br{}') \MO{s}(\br{}')  \dbr{} \dbr{}',
\end{equation}
the BSE matrix elements read
\begin{subequations}
\begin{align}
	\label{eq:LR_BSE-A}
	\ABSE{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \IS \qty[ 2 \ERI{ia}{jb}   -  \W{ij,ab}{\IS} ],
	\\ 
	\label{eq:LR_BSE-B}
	\BBSE{ia,jb}{\IS} & =  \IS \qty[ 2 \ERI{ia}{bj} -   \W{ib,aj}{\IS} ],
\end{align}
\end{subequations}
where $\eGW{p}$ are the $GW$ quasiparticle energies. 
%In the standard BSE implementation,  the screened Coulomb potential $\W{}{\IS}$ is taken to be static $(\omega \rightarrow 0)$. 
In the standard BSE approach, $\W{}{\IS}$ is built within the direct RPA scheme, \ie,
\begin{subequations}
\label{eq:wrpa}
\begin{align}
	\W{}{\IS}(\br{},\br{}') 
	& = \int \frac{\epsilon_{\IS}^{-1}(\br{},\br{}''; \omega=0)}{\abs*{\br{}' - \br{}''}} \dbr{}'' , 
	\\ 
	\epsilon_{\IS}(\br{},\br{}'; \omega) 
	& = \delta(\br{}-\br{}') - \IS  \int \frac{\chi_{0}(\br{},\br{}''; \omega)}{\abs*{\br{}' - \br{}''}} \dbr{}'' ,
\end{align}
\end{subequations}
with $\epsilon_{\IS}$ the dielectric function at coupling constant $\IS$ and $\chi_{0}$ the non-interacting polarizability.  In the occupied-to-virtual orbital product basis, the spectral representation of $\W{}{\IS}$ can be written as follows in the case of real spatial orbitals
\begin{multline}
\label{eq:W}
	\W{ij,ab}{\IS}(\omega) = \ERI{ij}{ab} +  2 \sum_m^{\Nocc \Nvir} \sERI{ij}{m} \sERI{ab}{m}
	\\
	\times  \qty(\frac{1}{\omega - \OmRPA{m}{\IS} + i \eta} - \frac{1}{\omega + \OmRPA{m}{\IS} - i \eta}),
\end{multline}
where the spectral weights at coupling strength $\IS$ read
\begin{equation}
	\sERI{pq}{m} =  \sum_i^{\Nocc} \sum_a^{\Nvir} \ERI{pq}{ia} (\bX{\IS}_m + \bY{\IS}_m)_{ia}.
\end{equation}
In the case of complex molecular orbitals, see Ref.~\onlinecite{Holzer_2019} for a correct use of complex conjugation in the spectral representation of $\W{}{}$.

In Eq.~\eqref{eq:W},  $\eta$ is a positive infinitesimal, and $\OmRPA{m}{\IS}$ are the direct (\ie, without exchange) RPA neutral excitation energies computed by solving the linear eigenvalue problem \eqref{eq:LR} with the following matrix elements
\begin{subequations}
\begin{align}
	\label{eq:LR_RPA-A}
	\ARPA{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eHF{a} - \eHF{i}) +   2 \IS \ERI{ia}{jb},
	\\ 
	\label{eq:LR_RPA-B}
	\BRPA{ia,jb}{\IS} & =   2 \IS \ERI{ia}{bj},
\end{align}
\end{subequations}
where $\eHF{p}$ are the HF orbital energies.

The relationship between the BSE formalism and the well-known RPAx (\ie, RPA with exchange) approach can be obtained by switching off the screening  
%namely setting  $\epsilon_{\IS}({\bf r},{\bf r}'; \omega) = \delta({\bf r}-{\bf r}')$ 
so that $\W{}{\IS}$ reduces to the bare Coulomb potential $\vc{}$. 
In this limit,  the $GW$ quasiparticle energies reduce to the Hartree-Fock (HF) eigenvalues, and Eqs.~\eqref{eq:LR_BSE-A} and \eqref{eq:LR_BSE-B} to the RPAx equations: 
\begin{subequations}
\begin{align}
	\label{eq:LR_RPAx-A}
	\ARPAx{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eHF{a} - \eHF{i}) + \IS \qty[ 2 \ERI{ia}{jb} -    \ERI{ij}{ab} ],
	\\ 
	\label{eq:LR_RPAx-B}
	\BRPAx{ia,jb}{\IS} & =  \IS \qty[ 2 \ERI{ia}{bj} - \ERI{ib}{aj} ].
\end{align}
\end{subequations}
%This allows to understand that the strength parameter $\IS$ enters twice in the $\IS W^{\IS}$ contribution, one time to renormalize the screening efficiency, and a second  time to renormalize the direct electron-hole interaction.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Ground-state BSE energy}
%\label{sec:BSE_energy}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The key quantity to define in the present context is the total BSE ground-state energy $\EBSE$.
Although this choice is not unique, \cite{Holzer_2018} we propose here to define it as
\begin{equation}
\label{eq:EtotBSE}
	\EBSE = \Enuc + \EHF + \EcBSE,
\end{equation}
where $\Enuc$ and $\EHF$ are the nuclear repulsion energy and ground-state HF energy (respectively), and
\begin{equation}
\label{eq:EcBSE}
	\EcBSE = \frac{1}{2} \int_0^1 \Tr(\bK \bP{\IS}) d\IS
\end{equation}
is the ground-state BSE correlation energy computed in the adiabatic connection framework, where 
\begin{equation}
\label{eq:K}
	\bK = 
	\begin{pmatrix}
		\btA{\IS=1}	&	\bB{\IS=1}	\\
		\bB{\IS=1}	&	\btA{\IS=1}	\\
	\end{pmatrix}
\end{equation}
is the interaction kernel \cite{Angyan_2011, Holzer_2018} [with $\tA{ia,jb}{\IS} = \IS \ERI{ia}{bj}$],
\begin{equation}
\label{eq:2DM}
	\bP{\IS} = 
	\begin{pmatrix}
		\bY{\IS} \T{(\bY{\IS})}		&	\bY{\IS} \T{(\bX{\IS})}	\\
		\bX{\IS} \T{(\bY{\IS})}		&	\bX{\IS} \T{(\bX{\IS})}	\\
	\end{pmatrix}
	-
	\begin{pmatrix}
		\bO		&	\bO	\\
		\bO		&	\bI	\\
	\end{pmatrix}
\end{equation} 
is the correlation part of the two-electron density matrix at interaction strength $\IS$, and $\Tr$ denotes the matrix trace.
Note that the present definition of the BSE correlation energy [see Eq.~\eqref{eq:EcBSE}], which we refer to as BSE@$GW$@HF here, has been named ``XBS'' for ``extended Bethe Salpeter'' by Holzer \textit{et al.} \cite{Holzer_2018}
Contrary to DFT where the electron density is fixed along the adiabatic path, in the present formalism, the density is not maintained as $\IS$ varies.
Therefore, an additional contribution to Eq.~\eqref{eq:EcBSE} originating from the variation of the Green's function along the adiabatic connection should be added.
However, as commonly done within RPA and RPAx, \cite{Toulouse_2009, Toulouse_2010, Holzer_2018} we shall neglect it in the present study. 

Equation \eqref{eq:EcBSE} can also be straightforwardly applied to RPA and RPAx, the only difference being the expressions of $\bA{\IS}$ and $\bB{\IS}$ used to obtain the eigenvectors $\bX{\IS}$ and $\bY{\IS}$ entering the definition of $\bP{\IS}$ [see Eq.~\eqref{eq:2DM}].
For RPA, these expressions have been provided in Eqs.~\eqref{eq:LR_RPA-A} and \eqref{eq:LR_RPA-B}, and their RPAx analogs in Eqs.~\eqref{eq:LR_RPAx-A} and \eqref{eq:LR_RPAx-B}.
In the following, we will refer to these two types of calculations as RPA@HF and RPAx@HF, respectively.
Finally, we will also consider the RPA@$GW$@HF scheme which consists in replacing the HF orbital energies in Eq.~\eqref{eq:LR_RPA-A} by the $GW$ quasiparticles energies.

Several important comments are in order here.
For spin-restricted closed-shell molecular systems around their equilibrium geometry (such as the ones studied here), it is rare to encounter singlet instabilities as these systems can be classified as weakly correlated.
However, singlet instabilities may appear in the presence of strong correlation, \eg, when the bond is stretched, hampering in particular the calculation of atomization energies. \cite{Holzer_2018}
Even for weakly correlated systems, triplet instabilities are much more common, but triplet excitations do not contribute to the correlation energy in the ACFDT formulation. \cite{Toulouse_2009, Toulouse_2010, Angyan_2011}

%\xavier{ IS THIS NEEDED NOW ? However, contrary to the plasmon formulation (an alternative expression of the BSE correlation energy), \cite{Schuck_Book, Sawada_1957b, Rowe_1968} the triplet excitations do not contribute in the ACFDT formulation, which is an indisputable advantage of this approach. 
%Indeed, although the plasmon and adiabatic connection formulations are equivalent for RPA, \cite{Sawada_1957b, Gell-Mann_1957, Fukuta_1964, Furche_2008} this is not the case at the BSE and RPAx levels. \cite{Toulouse_2009, Toulouse_2010, Angyan_2011, Li_2020}
%For RPAx, an alternative plasmon expression (equivalent to its ACFDT analog) can be found if exchange is also added to the interaction kernel [see Eq.~\eqref{eq:K}]. \cite{Angyan_2011}
%However, to the best of our knowledge, such alternative plasmon expression does not exist for BSE. }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Computational details}
%\label{sec:comp_details}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
All the $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies have been done using a (restricted) HF starting point, which is a very adequate choice in the case of the (small) systems that we have considered here.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986} calculations are employed as starting point to compute the BSE neutral excitations.
In the case of {\GOWO}, the quasiparticle energies have been obtained by linearizing the non-linear, frequency-dependent quasiparticle equation.
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018, Veril_2018}.
Finally, the infinitesimal $\eta$ has been set to zero for all calculations.
The numerical integration required to compute the correlation energy along the adiabatic path [see Eq.~\eqref{eq:EcBSE}] has been performed with a 21-point Gauss-Legendre quadrature. 
Comparison with the so-called plasmon (or trace) formula \cite{Furche_2008} at the RPA level has confirmed the excellent accuracy of this quadrature scheme over $\IS$. 

For comparison purposes, we have also computed the PES at the MP2, CC2 \cite{Christiansen_1995}, CCSD, \cite{Purvis_1982} and CC3 \cite{Christiansen_1995b} levels of theory using DALTON. \cite{dalton}
The computational cost of these methods, in their usual implementation, scale as $\order*{N^5}$, $\order*{N^5}$, $\order*{N^6}$, and $\order*{N^7}$, respectively.
All the other calculations have been performed with our locally developed $GW$ software. \cite{Loos_2018,Veril_2018}
As one-electron basis sets, we employ the Dunning family (cc-pVXZ) defined with cartesian gaussian functions. 
Unless, otherwise stated, the frozen-core approximation has not been enforced in our calculations in order to provide a fair comparison between methods.
However, we have found that the conclusions drawn in the present study hold within the frozen-core approximation (see {\SI}).

Because Eq.~\eqref{eq:EcBSE} requires the entire BSE singlet excitation spectrum for each quadrature point, we perform several complete diagonalization of the $\Nocc \Nvir \times \Nocc \Nvir$ BSE linear response matrix [see Eq.~\eqref{eq:small-LR}], which corresponds to a $\order{\Nocc^3 \Nvir^3} = \order{\Norb^6}$ computational cost.
This step is, by far, the computational bottleneck in our current implementation.
However, we are currently pursuing different avenues to lower this cost by computing the two-electron density matrix of Eq.~\eqref{eq:2DM} via a quadrature in frequency space. \cite{Duchemin_2019,Duchemin_2020}

%%% FIG 1 %%%
\begin{figure*}
	\includegraphics[width=0.49\linewidth]{H2_GS_VQZ}
	\includegraphics[width=0.49\linewidth]{LiH_GS_VQZ}
\caption{
Ground-state PES of \ce{H2} (left) and \ce{LiH} (right) around their respective equilibrium geometry obtained at various levels of theory with the cc-pVQZ basis set. 
%Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-H2-LiH}
}
\end{figure*}
%%% %%% %%%

%%% FIG 2 %%%
\begin{figure*}
	\includegraphics[height=0.35\linewidth]{LiF_GS_VQZ}
	\includegraphics[height=0.35\linewidth]{HCl_GS_VQZ}
\caption{
Ground-state PES of \ce{LiF} (left) and \ce{HCl} (right) around their respective equilibrium geometry obtained at various levels of theory with the cc-pVQZ basis set. 
%Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-LiF-HCl}
}
\end{figure*}
%%% %%% %%%

%%% FIG 3 %%%
\begin{figure*}
	\includegraphics[height=0.26\linewidth]{N2_GS_VQZ}
	\includegraphics[height=0.26\linewidth]{CO_GS_VQZ}
	\includegraphics[height=0.26\linewidth]{BF_GS_VQZ}
\caption{
Ground-state PES of the isoelectronic series \ce{N2} (left), \ce{CO} (center), and \ce{BF} (right) around their respective equilibrium geometry obtained at various levels of theory with the cc-pVQZ basis set. 
%Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-N2-CO-BF}
}
\end{figure*}
%%% %%% %%%

%%% FIG 4 %%%
\begin{figure}
	\includegraphics[width=\linewidth]{F2_GS_VQZ}
\caption{
Ground-state PES of \ce{F2} around its equilibrium geometry obtained at various levels of theory with the cc-pVQZ basis set. 
%Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-F2}
}
\end{figure}
%%% %%% %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Potential energy surfaces}
%\label{sec:PES}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to illustrate the performance of the BSE-based adiabatic connection formulation, we have computed the ground-state PES of several closed-shell diatomic molecules around their equilibrium geometry: \ce{H2}, \ce{LiH}, \ce{LiF}, \ce{HCl}, \ce{N2}, \ce{CO}, \ce{BF}, and \ce{F2}.
The PES of these molecules for various methods are represented in Figs.~\ref{fig:PES-H2-LiH}, \ref{fig:PES-LiF-HCl}, \ref{fig:PES-N2-CO-BF}, and \ref{fig:PES-F2}, while the computed equilibrium distances for various basis sets are gathered in Table \ref{tab:Req}.
Additional graphs for other basis sets can be found in the {\SI}.

%%% TABLE I %%%
\begin{table*}
\caption{
Equilibrium bond length (in bohr) of the ground state of diatomic molecules obtained at various levels of theory and basis sets. 
The reference CC3 and corresponding BSE@{\GOWO}@HF data are highlighted in bold black and bold red for visual convenience, respectively.
The values in parenthesis have been obtained by fitting a Morse potential to the PES.
}
\label{tab:Req}
	\begin{ruledtabular}
	\begin{tabular}{llcccccccc}
				&				&	\mc{8}{c}{Molecules}	\\
									\cline{3-10}
	Method		&	Basis		&	\ce{H2}		&	\ce{LiH}	&	\ce{LiF}	&	\ce{HCl}	&	\ce{N2}			&	\ce{CO}	&	\ce{BF}	&	\ce{F2}	\\
	\hline
	CC3				&	cc-pVDZ		&	1.438	&	3.043		&	3.012		&	2.435		&	2.114		&	2.166		&	2.444		&	2.740	\\
					&	cc-pVTZ		&	1.403	& 	3.011		&	2.961		&	2.413		&	2.079		&	2.143		&	2.392		&	2.669	\\
					&	cc-pVQZ		&\bb{1.402}	&\bb{3.019}		&\bb{2.963}		&\bb{2.403}		&\bb{2.075}		&\bb{2.136}		&\bb{2.390}		&\bb{2.663}	\\
	CCSD			&	cc-pVDZ		&	1.438	&	3.044		&	3.006		&	2.433		&	2.101		&	2.149		&	2.435		&	2.695	\\
					&	cc-pVTZ		&	1.403	&	3.012		&	2.954		&	2.409		&	2.064		&	2.126		&	2.382		&	2.629	\\
					&	cc-pVQZ		&	1.402	&	3.020		&	2.953		&	2.398		&	2.059		&	2.118		&	2.118		&	2.621	\\
	CC2				&	cc-pVDZ		&	1.426	&	3.046		&	3.026		&	2.427		&	2.146		&	2.187		&	2.444		&	2.710	\\
					&	cc-pVTZ		&	1.393	&	3.008		&	2.995		&	2.406		&	2.109		&	2.163		&	2.394		&	2.664	\\
					&	cc-pVQZ		&	1.391	&	2.989		&	2.982		&	2.396		&	2.106		&	2.156		&	2.393		&	2.665	\\
	MP2				&	cc-pVDZ		&	1.426	&	3.041		&	3.010		&	2.426		&	2.133		&	2.166		&	2.431		&	2.681	\\
					&	cc-pVTZ		&	1.393	&	3.004		&	2.968		&	2.405		&	2.095		&	2.144		&	2.383		&	2.636	\\
					&	cc-pVQZ		&	1.391	&	3.008		&	2.970		&	2.395		&	2.091		&	2.137		&	2.382		&	2.634	\\
	BSE@{\GOWO}@HF	&	cc-pVDZ		&	1.437	&	3.042		&	3.000		&	2.454		&	2.107		&	2.153		&	2.407		&	(2.698)	\\
					&	cc-pVTZ		&	1.404	&	3.023		&	(2.982)		&	2.410		&	2.068		&	2.116		&	(2.389)		&	(2.647)	\\
					&	cc-pVQZ		&\rb{1.399}	&\rb{3.017}		&\rb{(2.974)}	&\gb{(2.408)}	&\gb{(2.070)}	&\gb{(2.130)}	&\gb{(2.383)}	&\rb{(2.640)}\\
	RPA@{\GOWO}@HF	&	cc-pVDZ		&	1.426	&	3.019		&	2.994		&	2.436		&	2.083		&	2.144		&	2.403		&	(2.629)	\\
					&	cc-pVTZ		&	1.388	&	2.988		&	(2.965)		&	2.408		&	2.055		&	2.114		&	(2.370)		&	(2.584)	\\
					&	cc-pVQZ		&	1.382	&	2.997		&	(2.965)		&\gb{(2.389)}	&\gb{(2.045)}	&\gb{(2.110)}	&\gb{(2.367)}	&	(2.571)	\\
	RPAx@HF			&	cc-pVDZ		&	1.428	&	3.040		&	2.998		&	2.424		&	2.077		&	2.130		&	2.417		&	2.611	\\
					&	cc-pVTZ		&	1.395	&	3.003		&	2.943		&	2.400		&	2.046		&	2.110		&	2.368		&	2.568	\\
					&	cc-pVQZ		&	1.394	&	3.011		&	2.944		&\gb{(2.393)}	&\gb{(2.041)}	&\gb{(2.105)}	&\gb{(2.367)}	&\gb{(2.563)}	\\
	RPA@HF			&	cc-pVDZ		&	1.431	&	3.021		&	2.999		&	2.424		&	2.083		&	2.134		&	2.416		&	2.623	\\
					&	cc-pVTZ		&	1.388	&	2.978		&	2.939		&	2.396		&	2.045		&	2.110		&	2.362		&	2.579	\\
					&	cc-pVQZ		&	1.386	&	2.994		&	2.946		&\gb{(2.385)}	&\gb{(2.042)}	&\gb{(2.104)}	&\gb{(2.365)}	&\gb{(2.571)}	\\
	\end{tabular}
	\end{ruledtabular}
\end{table*}

Let us start with the two smallest molecules, \ce{H2} and \ce{LiH}, which are both held by covalent bonds.
Their corresponding PES computed with the cc-pVQZ basis are reported in Fig.~\ref{fig:PES-H2-LiH}.
For \ce{H2}, we take as reference the full configuration interaction (FCI) energies \cite{QP2} and we also report the MP2 curve and its third-order variant (MP3), which improves upon MP2 towards FCI.
RPA@HF and RPA@{\GOWO}@HF yield almost identical results, and significantly underestimate the FCI energy, while RPAx@HF and BSE@{\GOWO}@HF slightly over and undershoot the FCI energy, respectively, RPAx@HF being the best match in the case of \ce{H2}.
Interestingly though, the BSE@{\GOWO}@HF scheme yields a more accurate equilibrium bond length than any other method irrespectively of the basis set.
For example, with the cc-pVQZ basis set, BSE@{\GOWO}@HF is only off by $0.003$ bohr compared to FCI, while RPAx@HF, MP2, and CC2 underestimate the bond length by $0.008$, $0.011$, and $0.011$ bohr, respectively.
The RPA-based schemes are much less accurate, with even shorter equilibrium bond lengths.

Albeit the shallow nature of the \ce{LiH} PES, the scenario is almost identical for \ce{LiH} for which we report the CC2, CCSD and CC3 energies in addition to MP2.
In this case, RPAx@HF and BSE@{\GOWO}@HF nestle the CCSD and CC3 energy curves, and they are almost perfectly parallel.
Here again, the BSE@{\GOWO}@HF equilibrium bond length (obtained with cc-pVQZ) is extremely accurate ($3.017$ bohr) as compared to FCI ($3.019$ bohr).

The cases of \ce{LiF} and \ce{HCl} (see Fig.~\ref{fig:PES-LiF-HCl}) are interesting as they corresponds to strongly polarized bonds towards the halogen atoms which are much more electronegative than the first row elements.
For these ionic bonds, the performance of BSE@{\GOWO}@HF are terrific with an almost perfect match to the CC3 curve.
%For \ce{LiF}, the two curves starting to deviate a few tenths of bohr after the equilibrium geometry, but they remain tightly bound for much longer in the case of \ce{HCl}.
Maybe surprisingly, BSE@{\GOWO}@HF is on par with both CC2 and CCSD, and outperforms RPAx@HF by a big margin for these two molecules exhibiting charge transfer.
However, in the case of \ce{LiF}, the attentive reader would have observed a small glitch in the $GW$-based curves very close to their minimum.
As observed in Refs.~\onlinecite{vanSetten_2015,Maggio_2017,Loos_2018} and explained in details in Refs.~\onlinecite{Veril_2018,Duchemin_2020}, these irregularities, which makes particularly tricky the location of the minima, are due to ``jumps'' between distinct solutions of the $GW$ quasiparticle equation.
Including a broadening via the increasing the value of $\eta$ in the $GW$ self-energy and the screened Coulomb operator soften the problem, but does not remove it completely.
Note that these irregularities would be genuine discontinuities in the case of {\evGW}. \cite{Veril_2018}


Let us now look at the isoelectronic series \ce{N2}, \ce{CO}, and \ce{BF}, which have a decreasing bond order (from triple bond to single bond).
In that case again, the performance of BSE@{\GOWO}@HF are outstanding, as shown in Fig.~\ref{fig:PES-N2-CO-BF}, and systematically outperforms both CC2 and CCSD.

The \ce{F2} molecule is a notoriously difficult case to treat due to the weakness of its covalent bond (see Fig.~\ref{fig:PES-F2}), hence its relatively long equilibrium bond length ($2.663$ bohr at the CC3/cc-pVQZ level).
Similarly to what we have observed for \ce{LiF} and \ce{BF}, there is an irregularities near the minimum of the {\GOWO}-based curves.
However, BSE@{\GOWO}@HF is the closest to the CC3 curve

%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Conclusion}
%\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%
In this Letter, we hope to have illustrated that the ACFDT@BSE formalism is a promising methodology for the computation of accurate ground-state PES and their corresponding equilibrium structures.
To do so, we have shown that calculating the BSE correlation energy computed within the ACFDT framework yield extremely accurate PES around equilibrium.
(Their accuracy near the dissociation limit remains an open question.)
We have illustrated this for 8 diatomic molecules for which we have also computed reference ground-state energies using coupled cluster methods (CC2, CCSD, and CC3).
However, we have also observed that, in some cases, unphysical irregularities on the ground-state PES due to the appearance of a satellite resonance with a weight similar to that of the $GW$ quasiparticle peak.
This shortcoming, which is entirely due to the quasiparticle nature of the underlying $GW$ calculation, has been thoroughly described in several previous studies.\cite{vanSetten_2015,Maggio_2017,Loos_2018,Veril_2018,Duchemin_2020}
We believe that this central issue must be resolved if one wants to expand the applicability of the present method.
In the perspective of developing analytical nuclear gradients within the BSE@$GW$ formalism, we are currently investigating the accuracy of the ACFDT@BSE scheme for excited-state PES.
We hope to be able to report on this in the near future.

%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
%PFL would like to thank Julien Toulouse for enlightening discussions about RPA, and 
XB is indebted to Valerio Olevano for numerous discussions.
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''.} 
\end{acknowledgements}
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supporting Information}
%%%%%%%%%%%%%%%%%%%%%%%%
See {\SI} for additional potential energy curves with other basis sets and within the frozen-core approximation.

\bibliography{BSE-PES,BSE-PES-control}

\end{document}
