Here, we build on the works of Scuseria \textit{et al.} [\href{http://dx.doi.org/10.1063/1.3043729}{J.~Chem.~Phys.~\textbf{129}, 231101 (2008)}] and Berkelbach [\href{https://doi.org/10.1063/1.5032314}{J.~Chem.~Phys.~\textbf{149}, 041103 (2018)}] to show intimate connections between the Bethe-Salpeter equation (BSE) formalism and $GW$ approximation from many-body perturbation theory and coupled-cluster (CC) theory at the ground- and excited-state levels.
Similitudes between BSE@$GW$ and the similarity-transformed equation-of-motion CC method introduced by Nooijen are put forward.
The present work allows to easily transfer key developments and general knowledge gathered in CC theory to many-body perturbation theory.
In particular, it provides a clear path for the computation of ground- and excited-state properties (such as nuclear gradients) within the BSE framework.
The random-phase approximation (RPA), introduced by Bohm and Pines in the context of the uniform electron gas, \cite{Loos_2016} is a quasibosonic approximation where one treats fermion products as bosons.
In the particle-hole (ph) channel, which is quite popular in the electronic structure community, particle-hole fermionic excitations and deexcitations are assumed to be bosons.
Because RPA corresponds to a resummation of all ring diagrams, it is adequate in the high-density (or weakly correlated) regime and catch effectively long-range correlation effects (such as dispersion). \cite{Gell-Mann_1957,Nozieres_1958}
Roughly speaking, the Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988,Blase_2018,Blase_2020} of many-body perturbation theory \cite{Martin_2016} can be seen as a cheap and effective way of introducing correlation in order to go \textit{beyond} RPA physics.
In the ph channel, BSE is commonly performed on top of a $GW$ calculation \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019,Bruneval_2021} from which one extracts the quasiparticle energies as well as the dynamically-screened Coulomb potential $W$.
This scheme has been shown to be highly successful to compute low-lying excited states of various natures (charge transfer, Rydberg, valence, etc) in molecular systems with a very attractive accuracy/cost ratio.\cite{Rohlfing_1999a,Horst_1999,Puschnig_2002,Tiago_2003,Rocca_2010,Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018,Blase_2018,Liu_2020,Blase_2020,Holzer_2018a,Holzer_2018b,Loos_2020e,Loos_2021}
Interestingly, RPA has strong connections with coupled-cluster (CC) theory, the workhorse of molecular electronic structure is when one is looking for high accuracy. \cite{Freeman_1977,Scuseria_2008,Jansen_2010,Scuseria_2013,Peng_2013,Berkelbach_2018,Rishi_2020}
%At low orders, \eg, coupled-cluster with singles and doubles, CC can treat weakly correlated systems, while when one can afford to ramp up the excitation degree to triplets or quadruples, even strongly correlated systems can be treated, hence eschewing to resort to multireference methods.
%Link between BSE and STEOM-CC.
%A route towards the obtention of BSE gradients for ground and excited states.
In a landmark paper, Scuseria \textit{et al.}\cite{Scuseria_2008} have proven that rCCD is equivalent to RPAx for the computation of the correlation energy.
Assuming the existence of $\bX{}{-1}$ (which can be proven as long as the RPAx problem is stable), this proof can be quickly summarised starting from the ubiquitous RPAx eigensystem
The composite variable $\bx$ gathers spin and spatial ($\br$) variables.
The indices $i$, $j$, $k$, and $l$ are occupied (hole) orbitals; $a$, $b$, $c$, and $d$ are unoccupied (particle) orbitals; $p$, $q$, $r$, and $s$ indicate arbitrary orbitals; and $m$ labels single excitations or deexcitations.
This simple and elegant proof was subsequently extended to excitation energies by Berkerbach, \cite{Berkelbach_2018} who showed that similitudes between EE-EOM-rCCD and RPAx exist when the EOM space is restricted to the 1h1p configurations and only the two-body terms are dressed by rCCD correlation.
where $\bHN= e^{-\hT}\hH e^{\hT}- E_\text{CCD}$ is the normal-ordered similarity-transformed Hamiltonian, $\Psi_{i}^{a}$ are singly-excited determinants, the one-body terms are
Neglecting the effect of $\hT_2$ on the one-body terms [see Eqs.~\eqref{eq:cFab} and \eqref{eq:cFij}] and relying on the rCCD amplitudes in the two-body terms, Eq.~\eqref{eq:cWibaj}, yields
Although the excitation energies of this approximate EE-EOM-rCCD scheme are equal to the RPAx ones, it has been shown that the transition amplitudes (or residues) are distinct and only agrees at the lowest order in the Coulomb interaction. \cite{Emrich_1981,Berkelbach_2018}
This does not hold only for RPAx as it is actually quite general and can be applied to any ph problem with the structure of Eq.~\eqref{eq:RPA}, such as TD-DFT, BSE, and many others.
This has been also extended to the pp and hh sectors by Peng \textit{et al.}\cite{Peng_2013} and Scuseria \textit{et al.}\cite{Scuseria_2013} (see also Ref.~\onlinecite{Berkelbach_2018} for the extension to excitation energies for the pp and hh channels).
are the elements of the correlation part of the dynamically-screened Coulomb potential which is set to its static limit \ie, $W_{pq,rs}^\text{stat}= W_{pq,rs}^\text{c}(\omega=0)$.
In Eq.~\eqref{eq:W}, $\eta$ is a positive infinitesimal, the screened two-electron integrals are
and $\Om{m}{\dRPA}$ is the $m$th (positive) eigenvalues and $\bX{m}{\dRPA}+\bY{m}{\dRPA}$ is constructed from the corresponding eigenvectors of the direct (\ie, without exchange) RPA (dRPA) problem defined as
As readily seen in Eqs.~\eqref{eq:A_RPAx}, \eqref{eq:B_RPAx}, \eqref{eq:A_BSE} and \eqref{eq:B_BSE}, the only difference between RPAx and BSE lies in the definition of the matrix elements, where one includes, via the screening of the electron-electron interaction [see Eq.~\eqref{eq:W}], correlation effects at the BSE level.
Note also that, within BSE, the orbital energies in Eqs.~\eqref{eq:A_BSE} and \eqref{eq:B_BSE} are usually computed at the $GW$ level to be consistent with the introduction of $W$ in the BSE kernel.
For example, at the one-shot $GW$ level, \cite{Strinati_1980,Hybertsen_1985a,Hybertsen_1986,Godby_1988,Linden_1988,Northrup_1991,Blase_1994,Rohlfing_1995,Shishkin_2007} we have
can be obtained via a set of rCCD equations, where one substitutes in Eq.~\eqref{eq:rCCD} all the two-electron integrals $\dbERI{pq}{rs}$ by $\dbERI{pq}{rs}- W_{ps,qr}^\text{stat}$.
Similarly to the diagonalization of the eigensystem \eqref{eq:BSE}, these rCCD-based equations can be solved in $\order*{N^6}$ cost, and provides a clear path for the computation of ground- and excited-state properties (such as nuclear gradients) within the BSE framework.
However, there is a significant difference as the BSE involves $GW$ quasiparticle energies [see Eq.~\eqref{eq:A_BSE}], where some of the correlation has been already dressed, while the RPAx equations only involves (undressed) one-electron orbital energies
Connections between approximate IP- and EA-EOM-CC schemes and the $GW$ approximation have been already studied in details by Lange and Berkelbach, \cite{Lange_2018} but we believe that the present work proposes a different perspective on this particular subject, as we derive genuine CC equations and we do not decouple the 2h1p and 2p1h sectors.
As shown by Bintrim and Berkelbach, \cite{Bintrim_2021} the non-linear $GW$ equations can be recast as a set of linear equations, which reads in the Tamm-Dancoff approximation
In the CC language, in order to determine the 2h1p and 2p1h amplitudes, $t_{ija,p}^{\text{2h1p}}$ and $t_{iab,p}^{\text{2p1h}}$, it means that one must solve the following amplitude equations
The $G_0W_0$ quasiparticle energies can be easily obtained via the procedure described in Ref.~\onlinecite{Monino_2022} by solving the previous equation for each value of $p$ separately.
%It is interesting to note that, in the case where $t_{ij}^{ab} = 0$, the 1h1p block reduces to the well-known random-phase approximation (with exchange).
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Data availability statement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The data that supports the findings of this study are available within the article and its supplementary material.
