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

\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}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{\red}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\SupMat}{\textcolor{blue}{supplementary material}}

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

% coordinates
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bx}{\boldsymbol{x}}


% methods
\newcommand{\NO}[1]{\{#1\}}
\newcommand{\HF}{\text{HF}}
\newcommand{\GW}{GW}
\newcommand{\dRPA}{\text{dRPA}}
\newcommand{\RPAx}{\text{RPAx}}
\newcommand{\BSE}{\text{BSE}}
\newcommand{\CC}{\text{CC}}

% 
\newcommand{\Ne}{N}
\newcommand{\Norb}{K}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}

% operators
\newcommand{\hH}{\Hat{H}}
\newcommand{\hHN}{\Hat{H}_{\text{N}}}
\newcommand{\hFN}{\Hat{F}_{\text{N}}}
\newcommand{\hVN}{\Hat{V}_{\text{N}}}
\newcommand{\hh}{\Hat{h}}
\newcommand{\hf}{\Hat{f}}
\newcommand{\bH}{\Bar{H}}
\newcommand{\bHN}{\Bar{H}_{\text{N}}}
\newcommand{\tHN}{\Tilde{H}_{\text{N}}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hS}{\Hat{S}}
\newcommand{\cre}[1]{\Hat{a}_{#1}^{\dag}}
\newcommand{\ani}[1]{\Hat{a}_{#1}}
\newcommand{\ca}[2]{\Hat{a}^{#1}_{#2}}

\newcommand{\cF}{\mathcal{F}}
\newcommand{\cW}{\mathcal{W}}

% energies
\newcommand{\Enuc}{E^\text{nuc}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\EHF}{E_\text{HF}}
\newcommand{\ECC}{E_\text{CC}}

% orbital energies
\newcommand{\e}[2]{\epsilon_{#1}^{#2}}
\newcommand{\eHF}[1]{\epsilon_{#1}^\text{HF}}
\newcommand{\eKS}[1]{\epsilon_{#1}^\text{KS}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
\newcommand{\om}[2]{\omega_{#1}^{#2}}


% Matrix elements
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\SO}[1]{\psi_{#1}}
\newcommand{\ERI}[2]{\braket{#1}{#2}}
\newcommand{\sERI}[2]{(#1|#2)}
\newcommand{\dbERI}[2]{\mel{#1}{}{#2}}
\newcommand{\SigC}[1]{\Sigma^\text{c}_{#1}}

% Matrices
\newcommand{\bO}{\boldsymbol{0}}
\newcommand{\bI}{\boldsymbol{1}}
\newcommand{\bbH}{\Bar{\boldsymbol{H}}}
\newcommand{\bOm}[2]{\boldsymbol{\Omega}_{#1}^{#2}}
\newcommand{\bom}[2]{\boldsymbol{\omega}_{#1}^{#2}}
\newcommand{\bA}[2]{\boldsymbol{A}_{#1}^{#2}}
\newcommand{\bB}[2]{\boldsymbol{B}_{#1}^{#2}}
\newcommand{\bC}[2]{\boldsymbol{C}_{#1}^{#2}}
\newcommand{\bD}[2]{\boldsymbol{D}_{#1}^{#2}}
\newcommand{\bF}[2]{\boldsymbol{F}_{#1}^{#2}}
\newcommand{\bR}[2]{\boldsymbol{R}_{#1}^{#2}}
\newcommand{\bT}[2]{\boldsymbol{T}_{#1}^{#2}}
\newcommand{\bU}[2]{\boldsymbol{U}_{#1}^{#2}}
\newcommand{\bV}[2]{\boldsymbol{V}_{#1}^{#2}}
\newcommand{\bX}[2]{\boldsymbol{X}_{#1}^{#2}}
\newcommand{\bY}[2]{\boldsymbol{Y}_{#1}^{#2}}
\newcommand{\bZ}[2]{\boldsymbol{Z}_{#1}^{#2}}
\newcommand{\be}[2]{\boldsymbol{\epsilon}_{#1}^{#2}}
\newcommand{\bSig}[2]{\boldsymbol{\Sigma}_{#1}^{#2}}

% orbitals, gaps, etc
\newcommand{\IP}{I}
\newcommand{\EA}{A}
\newcommand{\HOMO}{\text{HOMO}}
\newcommand{\LUMO}{\text{LUMO}}
\newcommand{\Eg}{E_\text{g}}
\newcommand{\EgFun}{\Eg^\text{fund}}
\newcommand{\EgOpt}{\Eg^\text{opt}}
\newcommand{\EB}{E_B}

\newcommand{\ii}{\mathrm{i}}

% addresses
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LPT}{Laboratoire de Physique Th\'eorique, Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\ETSF}{European Theoretical Spectroscopy Facility (ETSF)}
\newcommand{\NEEL}{Universit\'e Grenoble Alpes, CNRS, Institut N\'EEL, F-38042 Grenoble, France}

\begin{document}	

\title{Connections between many-body perturbation and coupled-cluster theories}

\author{Ra\'ul \surname{Quintero-Monsebaiz}}
	\affiliation{\LCPQ}
\author{Enzo \surname{Monino}}
	\affiliation{\LCPQ}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
	\email{loos@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}

\begin{abstract}
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.
%\bigskip
%\begin{center}
%	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
%\end{center}
%\bigskip
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%
%\section{RPA Physics and Beyond}
%%%%%%%%%%%%%%%%%%%%%%
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 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}


%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 ubiquitous 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
\begin{subequations}
\begin{align}
	\label{eq:RPA_1}
	\bA{}{}	+ \bB{}{} \cdot \bT{}{} & = \bR{}{}
	\\
	\label{eq:RPA_2}
	-\bB{}{} - \bA{}{} \cdot \bT{}{} & = \bT{}{} \cdot \bR{}{}
\end{align}
\end{subequations}

Substituting Eq.~\eqref{eq:RPA_1} into Eq.~\eqref{eq:RPA_2} yields the following Riccati equation
\begin{equation}
	\bB{}{} + \bA{}{} \cdot \bT{}{} + \bT{}{} \cdot \bA{}{} + \bT{}{} \cdot \bB{}{} \cdot \bT{}{} = \bO
\end{equation}
that matches the well-known rCCD equations
\begin{multline}
\label{eq:rCCD}
	\dbERI{ij}{ab} 
	+ \qty( \e{a}{} + \e{b}{} - \e{i}{} - \e{j}{} ) t_{ij}^{ab}
	+ \sum_{kc} \dbERI{ic}{ak} t_{kj}^{cb} 
	\\
	+ \sum_{kc} \dbERI{kb}{cj} t_{ik}^{ac} 
	+ \sum_{klcd} \dbERI{kl}{cd} t_{ik}^{ac} t_{lj}^{db} = 0
\end{multline}
knowing that
\begin{subequations}
\begin{align}
	\label{eq:A_RPAx}
	A_{ia,jb} & = (\e{a}{} - \e{i}{}) \delta_{ij} \delta_{ab} + \dbERI{ib}{aj}
	\\ 
	\label{eq:B_RPAx}
	B_{ia,jb} & = \dbERI{ij}{ab}
\end{align}
\end{subequations}
We assume real quantities throughout this paper, $\e{p}{}$ is the one-electron energy associated with the Hartree-Fock spinorbital $\SO{p}(\bx)$ and 
\begin{equation}
	\label{eq:ERI}
	\braket{pq}{rs} = \iint \SO{p}(\bx_1) \SO{q}(\bx_2) \frac{1}{\abs{\br_1 - \br_2}} \SO{r}(\bx_1) \SO{s}(\bx_2)  d\bx_1 d\bx_2
\end{equation}
are two-electron repulsion integrals, while 
\begin{equation}
	\dbERI{pq}{rs} = \braket{pq}{rs} - \braket{pq}{sr}
\end{equation}
are their anti-symmetrized versions.
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$ is the total number.


There are various ways of computing the RPAx correlation energy, \cite{Angyan_2011} but the usual plasmon (or trace) formula yields
\begin{equation}
\label{eq:EcRPAx}
	\Ec^\text{RPAx} = \frac{1}{4} \Tr(\bOm{}{} - \bA{}{})
\end{equation}
and matches the rCCD correlation energy
\begin{equation}
	\Ec^\text{rCCD} = \frac{1}{4} \sum_{ijab} \dbERI{ij}{ab} t_{ij}^{ab} = \frac{1}{4} \Tr(\bB{}{} \cdot \bT{}{})
\end{equation}
because $\Tr(\bOm{}{} - \bA{}{}) = \Tr(\bR{}{} - \bA{}{}) = \Tr(\bB{}{} \cdot \bT{}{})$, as evidenced by Eq.~\eqref{eq:RPA_1}.
Note that the same expression as Eq.~\eqref{eq:EcRPAx} can be derived from the adiabatic connection fluctuation dissipation theorem.\cite{Angyan_2011}

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 1h1p configurations and only the two-body terms are dressed by rCCD correlation.

To be more specific, restricting ourselves to CCD, \ie, $\hT = \hT_2$, the elements of the 1h1p block of the EOM Hamiltonian read
\begin{equation}
	\mel*{ \Psi_{i}^{a} }{ \bHN }{ \Psi_{j}^{b} } = \cF_{ab} \delta_{ij} - \cF_{ij} \delta_{ab} + \cW_{jabi}
\end{equation}
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
\begin{subequations}
\begin{align}
	\cF_{ab} & = \e{a}{} \delta_{ab} - \frac{1}{2} \sum_{klc} \ERI{kl}{bc} t_{kl}^{ac}
	\\
	\cF_{ij} & = \e{i}{} \delta_{ij} + \frac{1}{2} \sum_{kcd} \ERI{ik}{cd} t_{jk}^{cd}
\end{align}
\end{subequations}
and the two-body term is
\begin{equation}
	\cW_{ibaj} = \ERI{ib}{aj} + \sum_{kc} \ERI{ik}{ac} t_{kj}^{cb}
\end{equation}
Neglecting the effect of $\hT_2$ on the one-body terms and relying on the rCCD ampltiudes in the two-body terms yields 
\begin{equation}
\begin{split}
	\mel*{ \Psi_{i}^{a} }{ \tHN }{ \Psi_{j}^{b} } 
	& = (\e{a}{} - \e{i}{}) \delta_{ij} \delta_{ab} + \dbERI{ib}{aj} + \sum_{kc} \dbERI{ik}{ac} t_{kj}^{cb}
	\\
	& = (\bA{}{} + \bB{}{} \cdot \bT{}{})_{ia,jb}
\end{split}
\end{equation}
which exactly matches Eq.~\eqref{eq:RPA_1}.
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 distincts 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 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
\begin{equation}
	\label{eq:BSE}
	\begin{pmatrix}
		\bA{}{\BSE}		&	\bB{}{\BSE}	\\
		-\bB{}{\BSE}	&	-\bA{}{\BSE}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{}{\BSE}	\\
		\bY{}{\BSE}	\\
	\end{pmatrix}
	=
	\begin{pmatrix}
		\bX{}{\BSE}	\\
		\bY{}{\BSE}	\\
	\end{pmatrix}
	\cdot
	\bOm{}{\BSE}
\end{equation}
where
\begin{subequations}
\begin{align}
	\label{eq:A_BSE}
	A_{ia,jb}^{\BSE} & = \delta_{ij} \delta_{ab} (\e{a}{\GW} - \e{i}{\GW}) + \dbERI{ib}{aj} - W_{ij,ba}^\text{stat}
	\\ 
	\label{eq:B_BSE}
	B_{ia,jb}^{\BSE} & = \dbERI{ij}{ab} - W_{ib,ja}^\text{stat}
\end{align}
\end{subequations}
and 
\begin{multline}
	\label{eq:W}
	W_{pq,rs}^\text{c}(\omega) = \sum_{m} \sERI{pq}{m} \sERI{rs}{m} 
	\\
	\times \qty[ \frac{1}{\omega - \Om{m}{\dRPA} + \ii \eta} - \frac{1}{\omega + \Om{m}{\dRPA} - \ii \eta} ]
\end{multline}
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
\begin{equation}
	\label{eq:sERI}
	\sERI{pq}{m} = \sum_{ia} \ERI{pi}{qa} \qty( \bX{m}{\dRPA} + \bY{m}{\dRPA} )_{ia}
\end{equation}
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
\begin{equation}
	\label{eq:dRPA}
	\begin{pmatrix}
		\bA{}{\dRPA}		&	\bB{}{\dRPA}	\\
		-\bB{}{\dRPA}	&	-\bA{}{\dRPA}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{}{\dRPA}	\\
		\bY{}{\dRPA}	\\
	\end{pmatrix}
	=
	\begin{pmatrix}
		\bX{}{\dRPA}	\\
		\bY{}{\dRPA}	\\
	\end{pmatrix}
	\cdot
	\bOm{}{\dRPA}
\end{equation}
with
\begin{subequations}
\begin{align}
	\label{eq:A_dRPA}
	A_{ia,jb}^{\dRPA} & = \delta_{ij} \delta_{ab} (\e{a}{} - \e{i}{}) + \ERI{ib}{aj}
	\\ 
	\label{eq:B_dRPA}
	B_{ia,jb}^{\dRPA} & = \ERI{ij}{ab}
\end{align}
\end{subequations}
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
\begin{equation}
	\label{eq:G0W0}
	\e{p}{\GW} = \e{p}{\HF} + Z_p \SigC{p}(\omega = \e{p}{\HF})
\end{equation}
where the elements of the correlation part of the dynamical self-energy are
\begin{equation}
\begin{split}
	\SigC{pq}(\omega) 
	& = \sum_{im} \frac{\sERI{pi}{m} \sERI{qi}{m}}{\omega - \e{i}{\HF} + \Om{m}{\dRPA} - \ii \eta}
	\\
	& + \sum_{am} \frac{\sERI{pa}{m} \sERI{qa}{m}}{\omega - \e{a}{\HF} - \Om{m}{\dRPA} + \ii \eta}
\end{split}
\end{equation}
and the renormalization factor is
\begin{equation}
	Z_p = \qty[ 1 - \eval{\pdv{\SigC{p}(\omega)}{\omega}}_{\omega = \e{p}{\HF}} ]^{-1}
\end{equation}


Therefore, following a similar procedure, one can show that the BSE correlation energy obtained using the trace formula 
\begin{equation}
	\Ec^\text{BSE} = \frac{1}{4} \Tr(\bOm{}{\BSE} - \bA{}{\BSE})
\end{equation}
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. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{$GW$ approximation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
\begin{subequations}
\begin{align}
	C^\text{2h1p}_{ija,klc} & = \qty[ \qty( \e{i}{} + \e{j}{} - \e{a}{}) \delta_{jl} \delta_{ac} - \ERI{jc}{al} ] \delta_{ik} 
	\\
	C^\text{2p1h}_{iab,kcd} & = \qty[ \qty( \e{a}{} + \e{b}{} - \e{i}{}) \delta_{ik} \delta_{ac} + \ERI{ak}{ic} ] \delta_{bd} 
\end{align}
\end{subequations}
and the corresponding coupling blocks read
\begin{align}
	V^\text{2h1p}_{p,klc} & = \ERI{pc}{kl}
	&
	V^\text{2p1h}_{p,kcd} & = \ERI{pk}{dc}
\end{align}
Introducing $\bT{}{\text{2h1p}} = \bY{}{\text{2h1p}} \cdot \bX{}{-1}$ and $\bT{}{\text{2p1h}} = \bY{}{\text{2p1h}} \cdot \bX{}{-1}$, we have
\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}
		\bI	\\
		\bT{}{\text{2h1p}}	\\
		\bT{}{\text{2p1h}}	\\
	\end{pmatrix}
	=
	\begin{pmatrix}
		\bI	\\
		\bT{}{\text{2h1p}}	\\
		\bT{}{\text{2p1h}}	\\
	\end{pmatrix}
	\cdot
	\bR{}{}
\end{equation}
with $\bR{}{} = \bX{}{} \cdot \bOm{}{} \cdot \bX{}{-1}$.
The matrix $\bX{}{-1}$ can be proven to exist when the system is stable.

This yields the three following equations
\begin{subequations}
\begin{align}
	\be{}{} + \bV{}{\text{2h1p}} \cdot \bT{}{\text{2h1p}} + \bV{}{\text{2p1h}} \cdot \bT{}{\text{2p1h}} & = \bR{}{}
	\label{eq:R}
	\\
	\T{(\bV{}{\text{2h1p}})} + \bC{}{\text{2h1p}} \cdot \bT{}{\text{2h1p}} & = \bT{}{\text{2h1p}} \cdot \bR{}{}
	\label{eq:T1R}
	\\
	\T{(\bV{}{\text{2p1h}})} + \bC{}{\text{2p1h}} \cdot \bT{}{\text{2p1h}} & = \bT{}{\text{2p1h}} \cdot \bR{}{}
	\label{eq:T2R}
\end{align}
\end{subequations}
Substituting Eq.~\eqref{eq:R} into Eqs.~\eqref{eq:T1R} and \eqref{eq:T2R} yields the two following coupled Riccati equations 
\begin{subequations}
\begin{align}
	\begin{split}
		\T{(\bV{}{\text{2h1p}})} 
		+ \bC{}{\text{2h1p}} \cdot \bT{}{\text{2h1p}} 
		- \bT{}{\text{2h1p}} \cdot \be{}{} 
		- \bT{}{\text{2h1p}} \cdot \bV{}{\text{2h1p}} \cdot \bT{}{\text{2h1p}}
		\\
		- \bT{}{\text{2h1p}} \cdot \bV{}{\text{2p1h}} \cdot \bT{}{\text{2p1h}}
		= \bO
	\end{split}
	\\
	\begin{split}
		\T{(\bV{}{\text{2p1h}})} 
		+ \bC{}{\text{2p1h}} \cdot \bT{}{\text{2p1h}} 
		- \bT{}{\text{2p1h}} \cdot \be{}{} 
		- \bT{}{\text{2p1h}} \cdot \bV{}{\text{2h1p}} \cdot \bT{}{\text{2h1p}}
		\\
		- \bT{}{\text{2p1h}} \cdot \bV{}{\text{2p1h}} \cdot \bT{}{\text{2p1h}}
		= \bO
	\end{split}
\end{align}
\end{subequations}
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
\begin{subequations}
\begin{align}
	\begin{split}
		r_{ija,p}^{\text{2h1p}}
		& = \ERI{pa}{ij}
		+ \Delta_{ija,p}^{\text{2h1p}} t_{ija,p}^{\text{2h1p}} 
		- \sum_{kc} \ERI{jc}{ak} t_{ikc,p}^{\text{2h1p}} 
		\\
		& - \sum_{klcq} t_{ija,q}^{\text{2h1p}} \ERI{qc}{kl} t_{klc,p}^{\text{2h1p}} 
		- \sum_{kcdq} t_{ija,q}^{\text{2h1p}} \ERI{qk}{dc} t_{kcd,p}^{\text{2p1h}} 
	\end{split}
	\\
	\begin{split}
		r_{iab,p}^{\text{2p1h}}
		& = \ERI{pi}{ba}
		+ \Delta_{iab,p}^{\text{2p1h}} t_{ija,p}^{\text{2p1h}} 
		+ \sum_{kc} \ERI{ak}{ic} t_{kcb,p}^{\text{2p1h}} 
		\\
		& - \sum_{klcq} t_{iab,q}^{\text{2p1h}} \ERI{qc}{kl} t_{klc,p}^{\text{2h1p}} 
		- \sum_{kcdq} t_{iab,q}^{\text{2p1h}} \ERI{qk}{dc} t_{kcd,p}^{\text{2p1h}} 
	\end{split}
\end{align}
\end{subequations}
with 
\begin{subequations}
\begin{align}
	\Delta_{ija,p}^{\text{2h1p}} & = \e{i}{} + \e{j}{} - \e{a}{} - \e{p}{}
	\\
	\Delta_{iab,p}^{\text{2p1h}} & = \e{a}{} + \e{b}{} - \e{i}{} - \e{p}{} 
\end{align}
\end{subequations}
\begin{subequations}
\begin{align}
	t_{ija,p}^{\text{2h1p}} & \leftarrow t_{ija,p}^{\text{2h1p}} - \frac{r_{ija,p}^{\text{2h1p}}}{\Delta_{ij}^{ap}}
	\\
	t_{iab,p}^{\text{2h1p}} & \leftarrow t_{iab,p}^{\text{2p1h}} - \frac{r_{iab,p}^{\text{2p1h}}}{\Delta_{ab}^{ip}}
\end{align}
\end{subequations}

The quasiparticle energies $\e{p}{GW}$ are thus provided by the eigenvalues of $\be{}{} + \bSig{}{\GW}$, where
\begin{equation}
	\bSig{}{\GW} = \bV{}{\text{2h1p}} \cdot \bT{}{\text{2h1p}} + \bV{}{\text{2p1h}} \cdot \bT{}{\text{2p1h}} 
\end{equation}
is the correlation part of the $GW$ self-energy.
The $G_0W_0$ quasiparticle energies can be easily obtained via the procedure described in Ref.~\onlinecite{Monino_2022}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Coupled-cluster theory}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The Schr\"odinger equation is
%\begin{equation}
%	\hH \Psi = E \Psi
%\end{equation}
%where the bare Hamiltonian reads
%\begin{equation}
%	\hH = \sum_{pq} h_{p}^{q} \cre{p} \ani{q} + \frac{1}{4} \sum_{pqrs} v_{pq}^{rs} \cre{p} \cre{q} \ani{s} \ani{r}
%\end{equation}
%where $h_{p}^{q} = \mel{p}{\hh}{q}$ are the core Hamiltonian elements and $v_{pq}^{rs} = \ERI{pq}{rs}$ are antisymmetrized two-electron integrals. 
%The creation and annihilation operators fulfils the following anticommutation relationships:
%\begin{align}
%	\cre{p} \cre{q} + \cre{q} \cre{p} & = 0
%	&
%	\ani{p} \ani{q} + \ani{q} \ani{p} & = 0
%	&
%	\cre{p} \ani{q} + \ani{p} \cre{q} & = \delta_{pq}
%\end{align}
%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 
%\begin{equation}
%	\hHN = \hFN + \hVN = \sum_{pq} f_{p}^{q} \NO{\cre{p}\ani{q}} + \frac{1}{4} \sum_{pqrs} v_{pq}^{rs} \NO{\cre{p}\cre{q}\ani{s}\ani{r}} 
%\end{equation}
%with $f_{p}^{q} = \mel{p}{\hf}{q}$ the elements of the Fock operator.
% and where the creation and annihilation operators are written compactly as
%\begin{equation}
%	\ca{pq\cdots}{rs\cdots} = \cre{p} \cre{q} \cdots \ani{s} \ani{r}
%\end{equation}

%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
%\begin{equation}
%\begin{split}
%	\bH
%	& = \hHN + \comm{\hHN}{\hT{}{}} 
%	+ \frac{1}{2!} \comm{\comm{\hHN}{\hT{}{}}}{\hT{}{}}
%	\\
%	& + \frac{1}{3!} \comm{\comm{\comm{\hHN}{\hT{}{}}}{\hT{}{}}}{\hT{}{}}
%	+ \frac{1}{4!} \comm{\comm{\comm{\hHN}{\hT{}{}}}{\hT{}{}}}{\hT{}{}}
%\end{split}
%\end{equation}
%which naturally truncate at the fourth term due to the two-body nature of the bare Hamiltonian.
%
%In coupled-cluster theory, it is usual to truncate the excitation operator at the $n$th order, \ie,
%\begin{equation}
%	\hT = \sum_{k=1}^{n} \hT_k
%\end{equation}
%where the normal-ordered $k$-fold ``neutral'' excitation operator
%\begin{equation}
%	\hT_{k} = \frac{1}{(k!)^2} \sum_{ij\cdots} \sum_{ab\cdots} t_{ij\cdots}^{ab\cdots} 
%	\NO{\cre{a}\cre{b}\cdots}{\ani{i}\ani{j}\cdots}
%\end{equation}
%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$.

%======================
%\subsection{EE-EOM-CCD}
%======================

%Up to 2h2p, the EE-EOM-CC matrix has the simple form
%\begin{equation}
%	\bdH_\text{2h2p}^\text{EE} = 
%	\begin{pmatrix}
%		\mel*{ \Phi_{i}^{a} }{ \bHN }{ \Phi_{k}^{c} }		&	\mel*{ \Phi_{i}^{a} }{ \bHN }{ \Phi_{kl}^{cd} } 
%		\\
%		\mel*{ \Phi_{ij}^{ab} }{ \bHN }{ \Phi_{kl}^{cd} }	&	\mel*{ \Phi_{ij}^{ab} }{ \bHN }{ \Phi_{kl}^{cd} } 
%		\\
%	\end{pmatrix}
%\end{equation}
%Restricting ourselves to CCD, the elements of the 1h1p block read
%\begin{equation}
%	\mel*{ \Phi_{i}^{a} }{ \bHN }{ \Phi_{j}^{b} } = \cF_{ab} \delta_{ij} - \cF_{ij} \delta_{ab} + \cW_{jabi}
%\end{equation}
%where the one-body terms are
%\begin{subequations}
%\begin{align}
%	\cF_{ab} & = \e{a}{} \delta_{ab} - \frac{1}{2} \sum_{klc} \ERI{kl}{bc} t_{kl}^{ac}
%	\\
%	\cF_{ij} & = \e{i}{} \delta_{ij} + \frac{1}{2} \sum_{kcd} \ERI{ik}{cd} t_{jk}^{cd}
%\end{align}
%\end{subequations}
%and the two-body term is
%\begin{equation}
%	\cW_{ibaj} = \ERI{ib}{aj} + \sum_{kc} \ERI{ik}{ac} t_{kj}^{cb}
%\end{equation}
%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).

%================================
%\subsection{IP- and EA-EOM-CCD}
%================================
%Up to 2h1p and 2p1h, the IP-EOM-CC and EA-EOM-CC matrices, respectively, have the simple form
%\begin{subequations}
%\begin{align}
%	\bdH_\text{2h1p}^\text{IP} & = 
%	\begin{pmatrix}
%		\mel*{ \Phi_{i} }{ \bHN }{ \Phi_{k} } &	\mel{ \Phi_{i} }{ \bHN }{ \Phi_{kl}^{c} } 
%		\\
%		\mel*{ \Phi_{ij}^{a} }{ \bHN }{ \Phi_{k} } &	\mel{ \Phi_{ij}^{a} }{ \bHN }{ \Phi_{kl}^{c} } 
%		\\
%	\end{pmatrix}
%	\\
%	\bdH_\text{2p1h}^\text{EA} & = 
%	\begin{pmatrix}
%		\mel*{ \Phi^{a} }{ \bHN }{ \Phi^{c} }		&	\mel*{ \Phi^{a} }{ \bHN }{ \Phi_{k}^{cd} } 
%		\\
%		\mel*{ \Phi_{i}^{ab} }{ \bHN }{ \Phi_{k}^{cd} }	&	\mel*{ \Phi_{i}^{ab} }{ \bHN }{ \Phi_{k}^{cd} } 
%		\\
%	\end{pmatrix}
%\end{align}
%\end{subequations}
%
%\begin{align}
%	\mel*{ \Phi_{i} }{ \bHN }{ \Phi_{k} } & = \titou{??}
%	\\
%	\mel*{ \Phi_{ij}^{a} }{ \bHN }{ \Phi_{kl}^{c} } & = \titou{??}
%	\\
%	\mel*{ \Phi_{i} }{ \bHN }{ \Phi_{kl}^{c} } & = \cW_{ickl}
%	\\
%	\mel*{ \Phi_{ij}^{a} }{ \bHN }{ \Phi_{k} } & = \ERI{ij}{ka}
%\end{align}
%and
%\begin{multline}
%	\cW_{iakl} = \ERI{ia}{kl} + \sum_{me} \ERI{im}{ke} t_{lm}^{ae} 
%	\\
%	- \sum_{me} \ERI{im}{le} t_{km}^{ae} 
%	+ \frac{1}{2} \sum_{ef} \ERI{ia}{ef} t_{kl}^{ef}
%\end{multline}
%================================
%\subsection{DIP- and DEA-EOM-CCD}
%================================
%Up to 3h1p and 3p1h, the DIP-EOM-CC and DEA-EOM-CC matrices, respectively, have the simple form
%\begin{subequations}
%\begin{align}
%	\bdH_\text{3h2p}^\text{DIP} & = 
%	\begin{pmatrix}
%		\mel*{ \Phi_{ij} }{ \bHN }{ \Phi_{lm} }			&	\mel*{ \Phi_{ij} }{ \bHN }{ \Phi_{lmn}^{d} } 
%		\\
%		\mel*{ \Phi_{ijk}^{a} }{ \bHN }{ \Phi_{lm} }	&	\mel*{ \Phi_{ijk}^{a} }{ \bHN }{ \Phi_{lmn}^{d} } 
%		\\
%	\end{pmatrix}
%	\\
%	\bdH_\text{3p2h}^\text{DEA} & = 
%	\begin{pmatrix}
%		\mel*{ \Phi_{i}^{abc} }{ \bHN }{ \Phi^{de} }	&	\mel*{ \Phi_{i}^{abc} }{ \bHN }{ \Phi_{l}^{def} } 
%		\\
%		\mel*{ \Phi^{ab} }{ \bHN }{ \Phi^{de} }			&	\mel*{ \Phi^{ab} }{ \bHN }{ \Phi_{l}^{def} } 
%	\end{pmatrix}
%\end{align}
%\end{subequations}
%
%\begin{subequations}
%\begin{align}
%	\mel*{ \Phi_{ij} }{ \bHN }{ \Phi_{kl} } & = - \cF_{ij} \delta_{kl} -  \delta_{ij} \cF_{kl}+ \cW_{ijkl}
%		\\
%	\mel*{ \Phi^{ab} }{ \bHN }{ \Phi^{cd} } & = + \cF_{ab} \delta_{cd} + \delta_{ab} \cF_{cd} + \cW_{abcd}
%\end{align}
%\end{subequations}
%with
%\begin{subequations}
%\begin{align}
%	\cW_{ijkl} & = \ERI{ij}{kl} + \sum_{a<b} \ERI{ij}{ab} t_{ij}^{ab}
%	\\
%	\cW_{abcd} & = \ERI{ab}{cd} + \sum_{i<j} \ERI{ab}{ij} t_{ij}^{ab}
%\end{align}
%\end{subequations}

%%%%%%%%%%%%%%%%%%%%%
%\section{Conclusion}
%%%%%%%%%%%%%%%%%%%%%
Here comes the conclusion.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Supplementary Material}
%\label{sec:supmat}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%
\acknowledgments{
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.

%%%%%%%%%
%\appendix
%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{CCvsMBPT}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
