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.
%When one is looking for high accuracy, the workhorse of molecular electronic structure is certainly coupled cluster theory.
%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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Bethe-Salpeter equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 RPAx eigensystem
\begin{equation}
\label{eq:RPA}
\begin{pmatrix}
\bA{}{}&\bB{}{}\\
-\bB{}{}& -\bA{}{}\\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{}{}\\
\bY{}{}\\
\end{pmatrix}
=
\begin{pmatrix}
\bX{}{}\\
\bY{}{}\\
\end{pmatrix}
\cdot
\bOm{}{}
\end{equation}
from which one gets, by introducing $\bT{}{}=\bY{}{}\cdot\bX{}{-1}$,
\begin{equation}
\begin{pmatrix}
\bA{}{}&\bB{}{}\\
-\bB{}{}& -\bA{}{}\\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bI\\
\bT{}{}\\
\end{pmatrix}
=
\begin{pmatrix}
\bI\\
\bT{}{}\\
\end{pmatrix}
\cdot
\bR{}{}
\end{equation}
where $\bR{}{}=\bX{}{}\cdot\bOm{}{}\cdot\bX{}{-1}$, or equivalently, the two following equations
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.
In the following, $O$ and $V$ are the number of occupied and virtual spinorbitals, respectively, and $N = O + V$.
There are various ways of computing the RPAx correlation energy, \cite{Angyan_2011} but the usual plasmon (or trace) formula (or, equivalently, the adiabatic connection fluctuation dissipation theorem) yields
because $\Tr(\bOm{}{}-\bA{}{})=\Tr(\bR{}{}-\bA{}{})=\Tr(\bB{}{}\cdot\bT{}{})$, as evidenced by Eq.~\eqref{eq:RPA_1}.
This proof was subsequently extended to excitation energies by Berkerbach, \cite{Berkelbach_2018} who showed that similitudes between EOM-rCCD and RPAx exist when the EOM space is restricted to 1h1p configurations and only the two-body terms are dressed by the rCCD correlation.
To be more specific, by restricting ourselves to CCD, \ie, $\hT=\hT_2$, the elements of the 1h1p block of the EOM Hamiltonian read
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
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 can be be also extended to the pp and hh sectors as shown by Peng \textit{et al.}\cite{Peng_2013} and Scuseria \textit{et al.}\cite{Scuseria_2013} (see also Ref.~\onlinecite{Berkelbach_2018} for an extension to excitation energies).
At the BSE level, and within the static approximation, one must solve a very similar linear eigenvalue problem
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.
Therefore, following a similar procedure, one can show that the BSE correlation energy obtained using the trace formula
can be obtained via a set of rCCD equations \eqref{eq:rCCD}, where one substitutes 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.
\titou{Comments on cost and extension to gradients.}
Following Berkelbach's analysis, one can extend the connection to excitation energies.
However, there is a significant difference as the BSE equations 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) HF orbital energies
It is therefore interesting to investigate if it possible also the recast the $GW$ equations as a set of CC-like equations.
Connections between EOM-CC and the $GW$ approximation have been already studied by Lange and Berkelbach, \cite{Lange_2018} but we believe that the present work proposes a different perspective on this particular subject.
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
\begin{equation}
\begin{pmatrix}
\be{}{}&\bV{}{\text{2h1p}}&\bV{}{\text{2p1h}}\\
\T{(\bV{}{\text{2h1p}})}&\bC{}{\text{2h1p}}&\bO\\
\T{(\bV{}{\text{2p1h}})}&\bO&\bC{}{\text{2p1h}}\\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{}{}\\
\bY{}{\text{2h1p}}\\
\bY{}{\text{2p1h}}\\
\end{pmatrix}
=
\begin{pmatrix}
\bX{}{}\\
\bY{}{\text{2h1p}}\\
\bY{}{\text{2p1h}}\\
\end{pmatrix}
\cdot
\bOm{}{}
\end{equation}
where $\be{}{}$ is a diagonal matrix gathered the orbital energies, the 2h1p and 2p1h matrix elements are
In the CC language, in order to determine the 2h1p and 2p1h amplitudes, $t_{ij}^{a}$ and $t_{i}^{ab}$, it means that one must solve the following amplitude equations
%In terms of normal-ordered (with respect to the reference wave function $\Psi_0$) quantities (which are indicated with curly braces), the Schr\"odinger equation reads
%\begin{equation}
% \hHN \Psi = \Delta E \Psi
%\end{equation}
%where $\Delta E = E - E_0 = E - \mel{\Psi_0}{\hH}{\Psi_0}$ is the correlation energy and the normal-ordered Hamiltonian $\hHN = \hH - E_0$ is
%Within coupled-cluster theory, thanks to the following exponential ansatz of the wave function
%\begin{equation}
% \Psi = e^{\hT} \Psi_0
%\end{equation}
%and the introducing the following similarity-transformed Hamiltonian
%\begin{equation}
% \bH = e^{-\hT} \hHN e^{\hT}
%\end{equation}
%one recasts the Schr\"odinger equation as
%\begin{equation}
% \bH \Psi_0 = \Delta E \Psi_0
%\end{equation}
%where $e^{\hT}$ is the wave operator and $\hT$ is the excitation operator.
%The Schr\"odinger equation is recasted as
%
%Matrix elements involving $\bH$ can be efficiently and elegantly evaluated thanks to Wick's theorem and the so-called Campbell-Baker-Hausdorff (BCH) formula which leads to a linear combination of nested communtators
%creates, when applied to $\ket*{\Psi_0}$, $k$-fold excited $\Ne$-electron determinants of the form $\ket*{\Psi_{ij\cdots}^{ab\cdots}}$.
%It is worth mentioning here that the excitation operators do commute with each other, \ie, $\comm{\hT_{k_1}}{\hT_{k_2}} = 0$ but do not commute with $\hHN$, \ie, $\comm{\hT_k}{\hHN} \neq 0$.
%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.