SRGGW/Notes/Notes.tex

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts,mleftright,bbold}
\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{\etal}{\textit{et al.}}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\ant}[1]{\textcolor{green}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{r\ed}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}

\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}}
\newcommand{\Sig}[2]{\Sigma_{#1}^{#2}}
\newcommand{\dRPA}{\text{dRPA}}

% coordinates
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bx}{\boldsymbol{x}}
\newcommand{\dbr}{d\br}
\newcommand{\dbx}{d\bx}

% methods
\newcommand{\GW}{\text{$GW$}}
\newcommand{\GT}{\text{$GT$}}	
\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{c}}
\newcommand{\x}{\text{x}}
\newcommand{\KS}{\text{KS}}
\newcommand{\HF}{\text{HF}}
\newcommand{\RPA}{\text{RPA}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
\newcommand{\sERI}[2]{(#1|#2)}
\newcommand{\e}[2]{\epsilon_{#1}^{#2}}

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

% operators
\newcommand{\hH}{\Hat{H}}
\newcommand{\hS}{\Hat{S}}
\newcommand{\ani}[1]{\hat{a}_{#1}}
\newcommand{\cre}[1]{\hat{a}_{#1}^\dagger}
\newcommand{\no}[2]{\mleft\{ \hat{a}_{#1}^{#2}\mright\} }

% energies
\newcommand{\Enuc}{E^\text{nuc}}
\newcommand{\Ec}[1]{E_\text{c}^{#1}}
\newcommand{\EHF}{E^\text{HF}}

% orbital energies
\newcommand{\eps}{\epsilon}
\newcommand{\reps}{\Tilde{\epsilon}}

% Matrix elements
\newcommand{\SigC}{\Sigma^\text{c}}
\newcommand{\rSigC}{\Tilde{\Sigma}^\text{c}}
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\SO}[1]{\psi_{#1}}
\newcommand{\eri}[2]{\braket{#1}{#2}}
\newcommand{\aeri}[2]{\mel{#1}{}{#2}} 
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\rbra}[1]{(#1|}
\newcommand{\rket}[1]{|#1)}


% Matrices
\newcommand{\bO}{\boldsymbol{0}}
\newcommand{\bI}{\boldsymbol{1}}
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\bSig}{\boldsymbol{\Sigma}}
\newcommand{\bSigC}{\boldsymbol{\Sigma}^{\text{c}}}
\newcommand{\be}{\boldsymbol{\epsilon}}
\newcommand{\bOm}{\boldsymbol{\Omega}}
\newcommand{\bA}{\boldsymbol{A}}
\newcommand{\bB}{\boldsymbol{B}}
\newcommand{\bC}[2]{\boldsymbol{C}_{#1}^{#2}}
\newcommand{\bD}{\boldsymbol{D}}
\newcommand{\bF}[2]{\boldsymbol{F}_{#1}^{#2}}
\newcommand{\bR}{\boldsymbol{R}}
\newcommand{\bU}{\boldsymbol{U}}
\newcommand{\bV}[2]{\boldsymbol{V}_{#1}^{#2}}
\newcommand{\bW}{\boldsymbol{W}}
\newcommand{\bX}[2]{\boldsymbol{X}_{#1}^{#2}}
\newcommand{\bY}[2]{\boldsymbol{Y}_{#1}^{#2}}
\newcommand{\bZ}[2]{\boldsymbol{Z}_{#1}^{#2}}
\newcommand{\bc}{\boldsymbol{c}}

% 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}

% shortcuts for greek letters
\newcommand{\si}{\sigma}
\newcommand{\la}{\lambda}


\newcommand{\RHH}{R_{\ce{H-H}}}
\newcommand{\ii}{\mathrm{i}}

\newcommand{\bEta}[1]{\boldsymbol{\eta}^{(#1)}(s)} 
\newcommand{\bHd}[1]{\bH_\text{d}^{(#1)}}
\newcommand{\bHod}[1]{\bH_\text{od}^{(#1)}} 

% addresses
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}

\begin{document}	

\title{Notes on the project: Similarity Renormalization Group formalism applied to Green's function theory}

\author{Antoine \surname{Marie}}
	\email{amarie@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}

\author{Pierre-Fran\c{c}ois \surname{Loos}}
	\email{loos@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}

%\begin{abstract}
%Here comes the abstract.
%\bigskip
%\begin{center}
%	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
%\end{center}
%\bigskip
%\end{abstract}

\maketitle

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

This document is a compilation of various notes related to the similarity renormalisation group (SRG) and its application to many-body perturbation theory (MBPT).
Before tackling its application to MBPT, we give the main SRG equations in Sec.~\ref{sec:intro}.
Then, we derive in details the perturbative expressions obtained by Evangelista by applying the SRG to the electronic Hamiltonian (see Sec.~\ref{sec:srg}).
Now turning to MBPT, its various flavors are presented in Sec.~\ref{sec:folded} and the corresponding unfolded equations are given in Sec.~\ref{sec:unfolded}.
Then, a SRG perturbative analysis is performed on general matrices in Sec.~\ref{sec:matrix_srg}.
Finally, in Sec.~\ref{sec:second_quant_mbpt} we investigate the possibility to find a second quantized effective Hamiltonian corresponding to the unfolded MBPT equations.

%=================================================================%
\section{The similarity renormalisation group}
\label{sec:srg}
%=================================================================%

The similarity renormalization group aims at continuously transforming an Hamiltonian to a diagonal form, or more often to a block-diagonal form.
Therefore, the transformed Hamiltonian
\begin{equation}
  \bH(s) = \bU(s) \, \bH \, \bU^\dag(s)
\end{equation}
depends on a flow parameter $s$.
The resulting Hamiltonian possess up to $N$-body operators with $N$ the number of particle.
\begin{equation}
  \bH(s) = E_0(s) + \bF(s) + \bV{}{}(s) + \bW(s) + \dots
\end{equation}
In the following, we will truncate every contribution superior to two-body operators.
We can easily derive an evolution equation for this Hamiltonian by taking the derivative of $\bH(s)$. This gives
\begin{equation}
  \label{eq:flowEquation}
  \dv{\bH(s)}{s} = \comm{\boldsymbol{\eta}(s)}{\bH(s)}
\end{equation}
where $\boldsymbol{\eta}(s)$, the flow generator, is defined as
\begin{equation}
\boldsymbol{\eta}(s) = \dv{\bU(s)}{s}	\bU^\dag(s)	= - \boldsymbol{\eta}^\dag(s) .
\end{equation}
To solve this equation at a cost inferior to the one of diagonalizing the initial Hamiltonian, one needs to introduce approximation for $\boldsymbol{\eta}(s)$.
Before doing so, we need to define what is the blocks to suppress in order to obtain a block-diagonal Hamiltonian.
Therefore, the Hamiltonian is separated in two parts as
\begin{equation}
  \bH(s) = \underbrace{\bH^\text{d}(s)}_{\text{diagonal}} + \underbrace{\bH^\text{od}(s)}_{\text{off-diagonal}}.
\end{equation}
By definition, we have the following condition on $\bH^\text{od}$
\begin{equation}
  \bH^\text{od}(\infty) = \boldsymbol{0}.
\end{equation}

In this work, we will use Wegner's canonical generator which is defined as
\begin{equation}
  \boldsymbol{\eta}^\text{W}(s) = \comm{\bH^\text{d}(s)}{\bH(s)} = \comm{\bH^\text{d}(s)}{\bH^\text{od}(s)}.
\end{equation}
This generator has the advantage of defining a true renormalisation scheme, \ie the coupling coefficients with the highest energy determinants are removed first.
One of the flaws of this generator is that it generates a \titou{stiff} set of ODE which is difficult to solve numerically.
However, here we consider analytical perturbative expressions so we will not be affected by this problem.

%=================================================================%
\section{The SRG electronic Hamiltonian}
\label{sec:electronic_ham}
%=================================================================%

In this part, we derive the perturbative expression for the SRG applied to the non-relativistic electronic Hamiltonian
\begin{equation}
  \label{eq:hamiltonianSecondQuant}
  \hH = \sum_{pq} f_{pq} \cre{p}\ani{q} + \frac{1}{4} \sum_{pqrs} \aeri{pq}{rs} \cre{p}\cre{q}\ani{r}\ani{s}
\end{equation}
which can also be written in normal order wrt a reference determinant as
\begin{equation}
  \label{eq:hamiltonianNormalOrder}
  \hH = E_0 + \sum_{pq} + f_p^q\no{q}{p} + \frac{1}{4} \sum_{pqrs}v_{pq}^{rs}\no{rs}{pq}.
\end{equation}

In this case, we want to decouple the reference determinant from every singly and doubly excited determinants.
Hence, we define the off-diagonal Hamiltonian as
\begin{equation}
  \label{eq:hamiltonianOffDiagonal}
  \hH^{\text{od}}(s) = \sum_{ia} f_i^a(s)\no{a}{i} + \frac{1}{4} \sum_{ijab}v_{ij}^{ab}(s)\no{ab}{ij}.
\end{equation}
Note that each coefficients depend on $s$.

The perturbative parameter $\la$ is such that
\begin{equation}
  \bH(0) = E_0(0) + \bF{}{}(0) + \la \bV{}{}(0)
\end{equation}
In addition, we know the following initial conditions.
We use the HF basis set of the reference such that $\bF{}{\text{od}}(0) = 0$ and $\bF{}{\mathrm{d}}(0)_{pq}=\delta_{pq}\epsilon_p$.
Therefore, we have
\begin{align}
  \bH^\text{d}(0)&=E_0(0) + \bF{}{\mathrm{d}}(0) + \la \bV{}{\mathrm{d}}(0) & \bH^\text{od}(0)&= \la V^{\mathrm{od}}(0)
\end{align}
Now, we want to compute the terms at each order of the following development
\begin{equation}
  \label{eq:hamiltonianPTExpansion}
  \bH(s) = \bH^{(0)}(s) + \bH^{(1)}(s) + \bH^{(2)}(s) + \bH^{(3)}(s) + \dots
\end{equation}
by integrating Eq.~\eqref{eq:flowEquation}.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Zeroth order Hamiltonian}
%%%%%%%%%%%%%%%%%%%%%%

First, we start by showing that the zeroth order Hamiltonian is independant of $s$ and therefore equal to $\bH^{(0)}(0)$
\begin{align}
  \bH(\delta s) &= \bH(0) + \delta s \dv{\bH(s)}{s}\bigg|_{s=0} + O(\delta s^2) \\
  \dv{\bH(s)}{s}\bigg|_{s=0} &= [ \boldsymbol{\eta}(0), \bH(0)] \\
\end{align}
However, we have seen that $\bH^\text{od}(0)$ is of order 1.
Hence, $\boldsymbol{\eta}(0)$ does not have a zero order contribution.
Which gives us the following equality
\begin{equation}
  \bH^{(0)}(\delta s) = \bH^{(0)}(0)
\end{equation}
meaning that the zeroth order Hamiltonian is independant of $s$.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{First order Hamiltonian}
%%%%%%%%%%%%%%%%%%%%%%

The right-hand side of \eqref{eq:flowEquation} has no zeroth order, so we want to compute its first order term.
To do so, we first need to compute the first order term of $\boldsymbol{\eta}(s)$.
We have seen that $\bH^\text{od}(0)$ has no zeroth order contribution so we have
\begin{equation}
  \boldsymbol{\eta}^{(1)}(s) = \comm{\bH^{\text{d},(0)}(0)}{\bH^{\text{od},(1)}(s)}
\end{equation}
According to the appendix the one-body part of $\boldsymbol{\eta}^{(1)}(s) $ has four contributions. However, two of them involve the two-body part of $\bH^{\text{d},(0)}(0)$ which is equal to zero.
In addition, the term $\left[A_{1}, B_{2}\right]_{p}^{q}$ involves the coefficients $A_i^a = f_i^a(0) = 0$.
So finally we only have one term
\begin{align}
  \eta_a^{i,(1)}(s) &= \comm{\bH_1^{\text{d},(0)}(0)}{\bH_1^{\text{od},(1)}(s)}_a^i  \\
                    &= \sum_r \left( f_a^r(0) f_r^{i,(1)}(s)  - f_r^i(0) f_a^{r,(1)}(s) \right) \\
                    &= (\epsilon_a - \epsilon_i)f_a^{i,(1)}(s)
\end{align}
Now turning to the two-body part of $\boldsymbol{\eta}^{(1)}(s) $ and once again two terms are zero because the two-body part of $\bH^{\text{d},(0)}(0)$ is equal to zero.
So we have 
\begin{align}
  \eta_{ab}^{ij,(1)}(s) &= \comm{\bH_1^{\text{d},(0)}(0)}{\bH_2^{\text{od},(1)}(s)}_{ab}^{ij}  \\
                    &= \sum_t [P(ab) f_a^t(0) v_{tb}^{ij,(1)}(s) - P(ij) f_t^i(0) v_{ab}^{tj,(1)}(s) ] \notag \\
                        &= \sum_t [ P(ab) \epsilon_a \delta_{at}v_{tb}^{ij,(1)}(s) - P(ij) \epsilon_i \delta_{it} v_{ab}^{tj,(1)}(s) ] \notag \\
                        &= P(ab) \epsilon_a v_{ab}^{ij,(1)}(s) - P(ij) \epsilon_i  v_{ab}^{ij,(1)}(s)  \notag\\
                        &= \left( \epsilon_a + \epsilon_b - \epsilon_i  - \epsilon_j \right) v_{ab}^{ij,(1)} \notag
\end{align}

We can now compute the first order contribution to Eq.~\eqref{eq:flowEquation}. We have seen that $\eta$ has no zeroth order contribution so
\begin{equation}
  \dv{\bH^{(1)}(s)}{s} = \comm{\boldsymbol{\eta}^{(1)}(s)}{\bH^{(0)}(s)}
\end{equation}
We start with the scalar contribution, \ie the PT1 energy
\begin{equation}
  \dv{E_0^{(1)}(s)}{s} = \mel{\phi}{\comm{\boldsymbol{\eta}_1^{(1)}(s)}{\bH_1^{(0)}(s)}}{\phi} + \mel{\phi}{\comm{\boldsymbol{\eta}_2^{(1)}(s)}{\bH_2^{(0)}(s)}}{\phi}
\end{equation}
where the second term is equal to zero.
Thus we have
\begin{align}
  \label{eq:diffEqScalPT1}
  \dv{E_0^{(1)}(s)}{s} &= \sum_{ip}\eta_i^{p,(1)}f_p^i(0) - \eta_p^{i,(1)}f_i^p(0) \\
                       &= \sum_i \epsilon_i(\eta_i^{i,(1)} - \eta_i^{i,(1)}) \notag \\
                       &= 0 \notag
\end{align}
Using the exact same reasoning as above we can show that there is only of the four terms of the one-body part of the commutator that is non-zero.
\begin{align}
  \dv{f_a^{i,(1)}(s)}{s} &= \comm{\boldsymbol{\eta}_1^{(1)}(s)}{\bH_1^{(0)}(s)}_a^i \\
                         &= \sum_r \eta_a^{r,(1)}(s)f_r^i(0) - \eta_r^{i,(1)}f_a^r(0)  \notag \\
                         &= (\epsilon_i - \epsilon_a) \eta_a^{i,(1)}(s) \notag \\
  &= -(\epsilon_i - \epsilon_a) ^2f_a^{i,(1)}(s)  \notag
\end{align}
The derivation for the two-body part to first order is once again similar
\begin{align}
  \dv{v_{ab}^{ij,(1)}(s)}{s} &= \comm{\boldsymbol{\eta}_2^{(1)}(s)}{\bH_1^{(0)}(s)}_{ab}^{ij} \\
                             &= -(\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b) ^2v_{ab}^{ij,(1)}(s)  \notag
\end{align}
These three differential equations can be integrated to obtain the analytical form of the Hamiltonian coefficients up to first order of perturbation theory.
\begin{align}
  E_0^{(1)}(s) &= E_0^{(1)}(0) = 0 \\
  f_a^{i,(1)}(s) &= f_a^{i,(1)}(0) e^{-s (\Delta_a^i )^2} = 0 \\
  v_{ab}^{ij,(1)}(s) &= v_{ab}^{ij,(1)}(0) e^{-s (\Delta_{ab}^{ij})^2} = \aeri{ij}{ab} e^{-s (\Delta_{ab}^{ij})^2}
\end{align}

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Second order Hamiltonian}
%%%%%%%%%%%%%%%%%%%%%%

Now turning to the differential equations, we start by computing the scalar part of Eq.~\eqref{eq:flowEquation}, \ie the differential equation for the second order energy.
\begin{align}
  \dv{E_0^{(2)}(s)}{s} & = \mel{\phi}{\comm{\boldsymbol{\eta}_1^{(2)}(s)}{\bH_1^{(0)}(s)}}{\phi} + \mel{\phi}{\comm{\boldsymbol{\eta}_2^{(2)}(s)}{\bH_2^{(0)}(s)}}{\phi} \\
  &+ \mel{\phi}{\comm{\boldsymbol{\eta}_1^{(1)}(s)}{\bH_1^{(1)}(s)}}{\phi} + \mel{\phi}{\comm{\boldsymbol{\eta}_2^{(1)}(s)}{\bH_2^{(1)}(s)}}{\phi}
\end{align}
The two first terms are equal to zero for the same reason as the PT1 scalar differential equation (see Eq.~\eqref{eq:diffEqScalPT1}).
In addition, the one-body hamiltonian has no first order contribution so 
\begin{align}
  \dv{E_0^{(2)}(s)}{s} &=  \mel{\phi}{\comm{\boldsymbol{\eta}_2^{(1)}(s)}{\bH_2^{(1)}(s)}}{\phi} \\
                       &= \frac{1}{4} \sum_{i j} \sum_{a b}\left(\eta_{i j}^{ab,(1)} H_{ab}^{i j,(1)} - H_{i j}^{ab,(1)} \eta_{ab}^{i j,(1)}\right) \\
                       &=\frac{1}{4} \sum_{i j} \sum_{a b}\left(\eta_{i j}^{ab,(1)}  - \eta_{ab}^{i j,(1)}\right) v_{ab}^{ij,(1)} \\
                       &= \frac{1}{4} \sum_{i j} \sum_{a b}\left(\Delta_{ab}^{ij}v_{ij}^{ab,(1)} - (-\Delta_{ab}^{ij} v_{ab}^{ij,(1)}) \right) \aeri{ij}{ab} e^{-s (\Delta_{ab}^{ij})^2} \\
                       &= \frac{1}{2} \sum_{i j} \sum_{a b} \Delta_{ab}^{ij} (v_{ab}^{ij,(1)})^2 \\
                       &= \frac{1}{2} \sum_{i j} \sum_{a b} \Delta_{ab}^{ij} \aeri{ij}{ab}^2 e^{-2s (\Delta_{ab}^{ij})^2}
\end{align}
After integration, using the initial condition $E_0^{(2)}(0)=0$, we obtain
\begin{equation}
  E_0^{(2)}(s) = \frac{1}{4} \sum_{i j} \sum_{a b} \frac{\aeri{ij}{ab}^2}{\Delta_{ab}^{ij}}\left(1-e^{-2s (\Delta_{ab}^{ij})^2}\right)
\end{equation}

We can continue the derivation further than Evangelista's paper.
To compute the second order contribution to the Hamiltonian coefficients, we first need to compute the second order contribution to $\boldsymbol{\eta}(s)$.
\begin{equation}
  \boldsymbol{\eta}^{(2)}(s) = \comm{\bH^{\text{d},(0)}(0)}{\bH^{\text{od},(2)}(s)} + \comm{\bH^{\text{d},(1)}(s)}{\bH^{\text{od},(1)}(s)}
\end{equation}
The expressions for the first commutator are computed analogously to the one of the previous subsection.
We focus on deriving expressions for the second term.
The one-body part of $\bH^{\text{od},(1)}(s)$ is equal to zero so two of the four terms contributing to the one-body part of  $\comm{\bH^{\text{d},(1)}(s)}{\bH^{\text{od},(1)}(s)}$ are zero.
In addition, the term $\comm{A_1}{B_2}$ is equal to zero as well because the coefficients $(A_1)_{i}^a$ are zero (see expression in Appendix).
So we have
\begin{align}
  \eta_a^{i,(2)}(s) &=  \comm{\bH_1^{\text{d},(0)}(0)}{\bH_1^{\text{od},(2)}(s)}_a^i + \comm{\bH_2^{\text{d},(1)}(s)}{\bH_2^{\text{od},(1)}(s)}_a^i \\
  &= (\epsilon_a - \epsilon_i)f_a^{i,(2)}(s) + \dots
\end{align}
Need to continue this derivation but this not needed for EPT2.

%=================================================================%
\section{The various flavors of MBPT}
\label{sec:mbpt}
% =================================================================%

The central equation of MBPT in practice is the following
\begin{equation}
  \label{eq:quasipart_eq}
  \bF{}{} + \bSig(\omega) = \omega \mathbb{1}.
\end{equation}
However, in order to use it we need to rely on approximations of the self-energy $\bSig(\omega)$.
                        
%%%%%%%%%%%%%%%%%%%%%%
\subsection{Self-energies and quasiparticle equations}
\label{sec:folded}
%%%%%%%%%%%%%%%%%%%%%%

In the following, we will focus on the GF(2), GW and GT approximations.
The GF($n$) formalism is defined such that the self-energy includes every diagram up to $n$-th order of M\"oller-Plesset perturbation theory.
\begin{align}
  \label{eq:GF2_selfenergy}
  \Sig{pq}{GF(2)}(\omega) &= \frac{1}{2} \sum_{ija} \frac{\aeri{pa}{ij}\aeri{qa}{ij}}{\omega + \eps _a -\eps_i -\eps_j - \ii \eta} \notag \\
  &+ \frac{1}{2} \sum_{iab} \frac{\aeri{pi}{ab}\aeri{qi}{ab}}{\omega + \eps _i -\eps_a -\eps_b + \ii \eta} \notag  
\end{align}
On the other hand, the GW self-energy is obtained by taking the RPA polarizability and removing the vertex correction in the exact definition of the self-energy.
\begin{equation}
  \label{eq:GW_selfenergy}
  \Sig{pq}{\GW}(\omega) = \sum_{im} \frac{\sERI{pi}{m} \sERI{qi}{m}}{\omega - \e{i}{} + \Om{m}{\dRPA} - \ii \eta} + \sum_{am} \frac{\sERI{pa}{m} \sERI{qa}{m}}{\omega - \e{a}{} - \Om{m}{\dRPA} + \ii \eta} \notag
\end{equation}
with
\begin{equation}
  \label{eq:GW_sERI}
    \sERI{pq}{m} = \sum_{ia} \eri{pq}{ia} \qty( \bX{m}{\dRPA} + \bY{m}{\dRPA} )_{ia} \notag
\end{equation}
Finally, the GT approximation corresponds to another approximation to the polarizability than in GW, namely the one coming from pp-hh-RPA
The corresponding self-energies read as
\begin{equation}
  \label{eq:GT_selfenergy}
  \Sig{pq}{\GT}(\omega) = \sum_{im} \frac{\eri{pi}{\chi^{N+2}_m}\eri{qi}{\chi^{N+2}_m}}{\omega + \e{i}{} - \Om{m}{N+2} + \ii \eta} + \sum_{am} \frac{\eri{pa}{\chi^{N-2}_m}\eri{qa}{\chi^{N-2}_m}}{\omega + \e{a}{} - \Om{m}{N-2} - \ii \eta} \notag
\end{equation}
with
\begin{align}
	\label{eq:GT_sERI}
  \eri{pq}{\chi^{N+2}_m} &= \sum_{c<d} \aeri{pq}{cd} \bX{cd,m}{N+2} + \sum_{k<l} \aeri{pq}{kl} \bY{kl,m}{N+2} \notag \\
  \eri{pq}{\chi^{N-2}_m} &= \sum_{c<d} \aeri{pq}{cd} \bX{cd,m}{N-2}  + \sum_{k<l} \aeri{pq}{kl} \bY{kl,m}{N-2} \notag
\end{align}
The two RPA problems giving the eigenvectors needed to build the GW and GT self-energies are given in Appendix~\ref{sec:rpa}.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{The unfolded equations}
\label{sec:unfolded}
%%%%%%%%%%%%%%%%%%%%%%

Following Schirmer for the GF(2) case or Bintrim \etal, the non-linear quasi-particle equations for each approximations can be unfolded in larger linear problems
\begin{equation}
  \label{eq:unfolded_equation}
  \bH  \bc_{(s)} = \epsilon_s \bc_{(s)}
\end{equation}
where $\bH$ depends on the approximation chosen for the self-energy.
For the three approximations considered here, the three matrices $\bH$ share the general form
\begin{equation}
  \label{eq:unfolded_matrice}
  \bH =
  \begin{pmatrix}
    \bF{}{}				                &	\bV{}{\text{2h1p}}		&	\bV{}{\text{2p1h}}	\\
    \T{(\bV{}{\text{2h1p}})}	&	\bC{}{\text{2h1p}}		&	\bO				\\
    \T{(\bV{}{\text{2p1h}})}	&	\bO					&	\bC{}{\text{2p1h}}	\\
  \end{pmatrix}
\end{equation}
The expression of the coupling blocks $\bV{}{}$ and the diagonal blocks $\bC{}{}$ in the different cases is given below.
\begin{itemize}
\item \textbf{GF(2)}
  \begin{align}
    V^\text{2h1p}_{p,klc} & = \frac{1}{\sqrt{2}}\aeri{pc}{kl}
    &
      V^\text{2p1h}_{p,kcd} & = \frac{1}{\sqrt{2}}\aeri{pk}{dc} \\
    C^\text{2h1p}_{ija,klc} & = \qty( \e{i}{} + \e{j}{} - \e{a}{}) \delta_{jl} \delta_{ac} \delta_{ik} 
    &
      C^\text{2p1h}_{iab,kcd} & = \qty( \e{a}{} + \e{b}{} - \e{i}{}) \delta_{ik} \delta_{ac} \delta_{bd} \notag
  \end{align}
\item \textbf{GW}
  \begin{align}
    V^\text{2h1p}_{p,klc} & = \eri{pc}{kl}
    &
      V^\text{2p1h}_{p,kcd} & = \eri{pk}{dc} \notag \\
    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} \notag & &
  \end{align}
\item \textbf{GT}
  \begin{align}
    V^\text{2h1p}_{p,klc} & = \aeri{pc}{kl}
    &
      V^\text{2p1h}_{p,kcd}&= \aeri{pk}{cd} \notag \\
    C^\text{2h1p}_{ija,klc} &= \qty[ \qty( \e{i}{} + \e{j}{} - \e{a}{}) \delta_{jl} \delta_{ac} - \aeri{ij}{kl} ] \delta_{ac} & & \\
    C^\text{2p1h}_{iab,kcd} &= \qty[ \qty( \e{a}{} + \e{b}{} - \e{i}{}) \delta_{ik} \delta_{ac} + \aeri{ab}{cd} ] \delta_{ik} & & \notag
  \end{align}
\end{itemize}

\textbf{\textcolor{red}{That would be nice to add electron-hole T matrix to see if it also correspond to one term that can be found in the CI below.}}

\textbf{\textcolor{red}{That would be nice to add the self-energies of the various post-RPA correction.}}

% =================================================================%
\section{A first quantization approach to SRG for MBPT}
\label{sec:matrix_srg}
%=================================================================%

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Initial conditions}
%%%%%%%%%%%%%%%%%%%%%%

Finding a second quantized effective Hamiltonian for MBPT is far from being trivial so we start the project with matrix perturbation theory.
A general upfolded MBPT matrix can be written as
\begin{equation}
  \label{eq:H_MBPT}
  H =
         \begin{pmatrix}
         \bF{}{} & \bV{}{} \\
         \bV{}{\dagger} & \bC{}{} 
         \end{pmatrix}
\end{equation}
Using SRG language, we define the diagonal and off-diagonal parts as 
\begin{equation}
  \label{eq:H_MBPT_partitioning}
  H(0) =
  \begin{pmatrix}
    \bF{}{} & \bO \\
    \bO & \bC{\text{d}}{} 
  \end{pmatrix}
  + \lambda
  \begin{pmatrix}
    \bO & \bV{}{} \\
    \bV{}{\dagger} & \bC{\text{od}}{}
  \end{pmatrix}
\end{equation}
which gives the following conditions
\begin{align}
  \bHd{0}(0) &=   \begin{pmatrix}
    \bF{}{} & \bO \\
    \bO & \bC{\text{d}}{} 
  \end{pmatrix}  & \bHod{0}(0) &= \bO \notag \\
  \bHd{1}(0) &= \bO & \bHod{1}(0) &=   \begin{pmatrix}
    \bO & \bV{}{} \\
    \bV{}{\dagger} & \bC{\text{od}}{} \notag
  \end{pmatrix}  
\end{align}
At this point, we aren't sure if the off-diagonal part of $\bC{}{}$ should be included or not in the off-diagonal part of the Hamiltonian $\bH_\text{od}$.
In the following derivation, we choose to do so because the other case can be obtained simply by taking $\bC{\text{od}}{} = \boldsymbol{0}$ and $\bC{\text{d}}{} = \bC{}{}$.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Zeroth order Hamiltonian}
%%%%%%%%%%%%%%%%%%%%%%

The zeroth-order commutator of the Wegner generator therefore gives
\begin{equation}
  \bEta{0} = \comm{\bHd{0}}{\bHod{0}} = \bO
\end{equation}
and similarly
\begin{equation}
  \dv{\bH^{(0)}}{s} = \comm{\bEta{0}}{\bH^{(0)}} = \bO
\end{equation}
Finally, we have
\begin{equation}
  \color{red}{\boxed{
    \color{black}{\bH^{(0)}(s) = \bH^{(0)}(0)}
    }}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%
\subsection{First order Hamiltonian}
%%%%%%%%%%%%%%%%%%%%%%

Now turning to the first-order contribution to the MBPT matrix, we start by computing the first order part of the Wegner generator.

\begin{align}
  &\bEta{1} = \comm{\bHd{0}}{\bHod{1}} \\
  &= \begin{pmatrix}
    \bO & \bF{}{(0)}\bV{}{(1)} - \bV{}{(1)}\bC{\text{d}}{(0)}\\
    \bC{\text{d}}{(0)}\bV{}{(1),\dagger} - \bV{}{(1),\dagger} \bF{}{(0)} & \bC{\text{d}}{(0)} \bC{\text{od}}{(1)} - \bC{\text{od}}{(1)} \bC{\text{d}}{(0)} \notag
  \end{pmatrix}  
\end{align}

\begin{align}
  \dv{\bH^{(1)}}{s} &= \comm{\bEta{1}}{\bHd{0}} =  \begin{pmatrix}
        \dv{\bF{}{(1)}}{s} & \dv{\bV{}{(1)}}{s} \\
        \dv{\bV{}{(1),\dagger}}{s} & \dv{\bC{}{(1)}}{s}
      \end{pmatrix}  \\
  \dv{\bF{}{(1)}}{s}  &= \bO  \\
  \dv{\bV{}{(1)}}{s} &= 2 \bF{}{(0)}\bV{}{(1)}\bC{\text{d}}{(0)} - (\bF{}{(0)})^2\bV{}{(1)} - \bV{}{(1)}(\bC{\text{d}}{(0)})^2 \\
  \dv{\bV{}{(1),\dagger}}{s} &= 2 \bC{\text{d}}{(0)}\bV{}{(1),\dagger}\bF{}{(0)} - \bV{}{(1),\dagger}(\bF{}{(0)})^2 - (\bC{\text{d}}{(0)})^2\bV{}{(1),\dagger} \\
  \dv{\bC{}{(1)}}{s}  &= 2 \bC{\text{d}}{(0)}\bC{\text{od}}{(1)}\bC{\text{d}}{(0)}- (\bC{\text{d}}{(0)})^2\bC{\text{od}}{(1)} - \bC{\text{od}}{(1)}(\bC{\text{d}}{(0)})^2
\end{align}
The last two equations can be solved differently depending on the form of $\bF{}{}$ and $\bC{}{}$.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Second order Hamiltonian}
%%%%%%%%%%%%%%%%%%%%%%

Recalling that $\bHod{0} = \bO$ and $\bHd{1} = \bO$, we derive 
\begin{align}
  &\bEta{2} = \comm{\bHd{0}}{\bHod{2}} + \comm{\bHd{1}}{\bHod{1}} \\
  &= \comm{\bHd{0}}{\bHod{2}} \notag \\
  &= \begin{pmatrix}
    \bO & \bF{}{(0)}\bV{}{(2)} - \bV{}{(2)}\bC{\text{d}}{(0)}\\
    \bC{\text{d}}{(0)}\bV{}{(2),\dagger} - \bV{}{(2),\dagger}\bF{}{(0)} & \bC{\text{d}}{(0)} \bC{\text{od}}{(2)} - \bC{\text{od}}{(2)} \bC{\text{d}}{(0)} \notag
  \end{pmatrix}  
\end{align}

\begin{align}
  \dv{\bH^{(2)}}{s} &= \comm{\bEta{2}}{\bHd{0}} + \comm{\bEta{1}}{\bHd{1}}  \\
  &=  \begin{pmatrix}
        \dv{\bF{}{(2)}}{s} & \dv{\bV{}{(2)}}{s} \\
        \dv{\bV{}{(2),\dagger}}{s} & \dv{\bC{}{(2)}}{s} 
      \end{pmatrix} \notag \\
  \dv{\bF{}{(2)}}{s}  &= \bF{}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bF{}{(0)} - 2 \bV{}{(1)}\bC{\text{d}}{(0)}\bV{}{(1),\dagger}\\
  \dv{\bC{}{(2)}}{s}  &= 2 \bC{\text{d}}{(0)}\bC{\text{od}}{(2)}\bC{\text{d}}{(0)}- (\bC{\text{d}}{(0)})^2\bC{\text{od}}{(2)} - \bC{\text{od}}{(2)}(\bC{\text{d}}{(0)})^2  \\
  &-2 \bC{\text{d}}{(1)}\bC{\text{od}}{(0)}\bC{\text{d}}{(1)} + (\bC{\text{d}}{(1)})^2\bC{\text{od}}{(0)} + \bC{\text{od}}{(0)}(\bC{\text{d}}{(1)})^2 \notag \\
  &+ \bC{\text{d}}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bC{\text{d}}{(0)} - 2 \bV{}{(1)}\bF{}{(0)}\bV{}{(1),\dagger} \notag \\
  \dv{\bV{}{(2)}}{s} &= 2 \bF{}{(0)}\bV{}{(2)}\bC{\text{d}}{(0)} - (\bF{}{(0)})^2\bV{}{(2)} - \bV{}{(2)}(\bC{\text{d}}{(0)})^2 \\
  &- 2 \bV{}{(1)} \bC{\text{d}}{(0)} \bC{\text{od}}{(1)}  + \bF{}{(0)} \bV{}{(1)} \bC{\text{od}}{(1)} + \bV{}{(1)} \bC{\text{od}}{(1)} \bC{\text{d}}{(0)} \notag \\
  \dv{\bV{}{(2),\dagger}}{s} &= 2 \bC{\text{d}}{(0)}\bV{}{(2),\dagger}\bF{}{(0)} - \bV{}{(2),\dagger}(\bF{}{(0)})^2 - (\bC{\text{d}}{(0)})^2\bV{}{(2),\dagger} \\
  &- 2 \bC{\text{od}}{(1)} \bC{\text{d}}{(0)} \bV{}{(1),\dagger} +  \bC{\text{od}}{(1)} \bV{}{(1),\dagger}  \bF{}{(0)} + \bC{\text{d}}{(0)} \bC{\text{od}}{(1)} \bV{}{(1),\dagger} \notag
\end{align}

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Downfolding the SRG-transformed matrix}
%%%%%%%%%%%%%%%%%%%%%%

In order to choose what to do with $\bC{\text{od}}{}$ we look at the downfolded SRG quasiparticle equation.

\begin{equation}
  \label{eq:H_SRGMBPT}
  H(s) =
  \begin{pmatrix}
    \bF{}{(0)}(0) + \bF{}{(2)}(s) & \bV{}{(1)}(s)  + \bV{}{(2)}(s) \\
    \bV{}{(1),\dagger}(s)  + \bV{}{(2),\dagger}(s) & \bC{}{(0)}(0) +\bC{}{(2)}(s) 
  \end{pmatrix}
\end{equation}

\begin{widetext}
  \begin{equation}
    \left\{
      \begin{aligned}
        (\bF{}{(0)}(0) + \bF{}{(2)}(s)) \bR^{1h/1p} + (\bV{}{(1)}(s) + \bV{}{(2)}(s)) \bR^{2h1p/2p1h} &= \omega \bR^{1h/1p} \\
        (\bV{}{(1),\dagger}(s) + \bV{}{(2),\dagger}(s))  \bR^{1h/1p} + (\bC{}{(0)}(0) +\bC{}{(2)}(s) ) \bR^{2h1p/2p1h}&=   \omega \bR^{2h1p/2p1h}             
      \end{aligned}
    \right.
  \end{equation}
  
  \begin{equation}
    (\bF{}{(0)}(0) + \bF{}{(2)}(s)) + (\bV{}{(1)}(s) + \bV{}{(2)}(s)) (\omega \mathbb{1} - \bC{\text{d}}{(0)}(0) - \bC{\text{od}}{(1)}(s)  -\bC{}{(2)}(s)  )^{-1}  (\bV{}{(1),\dagger}(s) + \bV{}{(2),\dagger}(s)) \bR^{1h/1p} = \omega \bR^{1h/1p}   
  \end{equation}

  If we want to truncate the quasiparticle equation to the second order we obtain 

    \begin{equation}
    (\bF{}{(0)}(0) + \bF{}{(2)}(s)) + \bV{}{(1)}(s)(\omega \mathbb{1} - \bC{\text{d}}{(0)}(0))^{-1} \bV{}{(1),\dagger}(s)  \bR^{1h/1p} = \omega \bR^{1h/1p}   
  \end{equation}
\end{widetext}
So if we choose to put the off-diagonal part of $\bC{}{}$ in the off-diagonal $\bH{}{}$ we see that the off diagonal part of $\bC{}{}$ is not present in the second order quasi-particle equation.
We believe that this is not desirable.
In the following, we will integrate order by order the differential equations obtained above in the case $\bC{\text{od}}{} = \boldsymbol{0}$ and $\bC{\text{d}}{} = \bC{}{}$.
The expression in the other case are given in Appendix~\ref{sec:diagC}.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{Integrating order by order}
%%%%%%%%%%%%%%%%%%%%%%

In the following, upper case indices correspond to the 2h1p and 2p1h sectors while lower case indices correspond to the 1h and 1p sectors.  Also the $\Delta\eps_R$ corresponds to the diagonal elements of the 2h1p and 2p1h sectors.

\subsubsection{First order}

\begin{align}
  \dv{\bF{}{(1)}}{s}  &= \bO  \Longleftrightarrow \bF{}{(1)}(s) = \bF{}{(1)}(0) \Longleftrightarrow \color{red}{\boxed{\color{black}{\bF{}{(1)}(s)= \bO}}} \\
  \dv{\bV{}{(1)}}{s} &= 2 \bF{}{(0)}\bV{}{(1)}\bC{}{(0)} - (\bF{}{(0)})^2\bV{}{(1)} - \bV{}{(1)}(\bC{}{(0)})^2 \\
  \dv{\bV{}{(1),\dagger}}{s} &= 2 \bC{}{(0)}\bV{}{(1),\dagger}\bF{}{(0)} - \bV{}{(1),\dagger}(\bF{}{(0)})^2 - (\bC{}{(0)})^2\bV{}{(1),\dagger} \\
  \dv{\bC{}{(1)}}{s}  &= \bO  \Longleftrightarrow \bC{}{(1)}(s) = \bC{}{(1)}(0) \Longleftrightarrow \color{red}{\boxed{\color{black}{\bC{}{(1)}(s)= \bO}}} 
\end{align}
The differential equation for the coupling blocks can be solved in the GF(2) case because in this case $\bC{}{}$ is diagonal (see Appendix~\ref{sec:diagC}).
However, in the general case this matrix differential equation is not trivial to solve.
In practice, we often resort to the \GOWO~or \evGW~schemes which implies that we only need to solve the above equations for $\bF{}{} = \epsilon_p$.
\begin{align}
  \dv{\bV{}{(1)}}{s} &= \bV{}{(1)} (2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2) \\
  \dv{\bV{}{(1),\dagger}}{s} &=  (2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2) \bV{}{(1),\dagger} \\
\end{align}
These matrix differential equations can be solved if we know how to diagonalize $2 \epsilon_p \bC{}{(0)} - \epsilon_p^2- (\bC{}{(0)})^2$.
We know how to diagonalize $\bC{}{(0)}$ so we know how to diagonalize polynomial of $\bC{}{(0)}$.

\textcolor{red}{\textbf{TODO Give analytical expression for the different cases.}}

\subsubsection{Second order}

\begin{align}
  \dv{\bF{}{(2)}}{s}  &= \bF{}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bF{}{(0)} - 2 \bV{}{(1)}\bC{}{(0)}\bV{}{(1),\dagger} \\
  \dv{\bC{}{(2)}}{s}  &= \bC{}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bC{}{(0)} - 2 \bV{}{(1)}\bF{}{(0)}\bV{}{(1),\dagger} \\
  \dv{\bV{}{(2)}}{s} &= 2 \bF{}{(0)}\bV{}{(2)}\bC{}{(0)} - (\bF{}{(0)})^2\bV{}{(2)} - \bV{}{(2)}(\bC{}{(0)})^2 \\
  \dv{\bV{}{(2),\dagger}}{s} &= 2 \bC{}{(0)}\bV{}{(2),\dagger}\bF{}{(0)} - \bV{}{(2),\dagger}(\bF{}{(0)})^2 - (\bC{}{(0)})^2\bV{}{(2),\dagger}
\end{align}

The two first equations can be solved by simple integrations.
The two last equations admit the same solutions as the first order coupling blocks differential equations with different initial conditions.

% =================================================================%
\section{Towards second quantized effective Hamiltonians for MBPT?}
\label{sec:second_quant_mbpt}
%=================================================================%

The many-body perturbation theory formalism and its various approximations are naturally derived using time-dependent Feynman diagrams.
These derivation are quite different from wave function methods like configuration interaction (CI) and coupled-cluster (CC) which are naturally expressed in second quantization.
One can study the link between these formalisms by expanding the MBPT Feynman diagrams into time-independent Goldstone diagrams and then compare them to the ones that appear in WFT.
However, that would be valuable to extend this connection by expressing the MBPT approximations in the second quantization formalism.
This is the aim of this section.

%%%%%%%%%%%%%%%%%%%%%%
\subsection{The IP/EA CI}
\label{sec:ip_ea_ci}
%%%%%%%%%%%%%%%%%%%%%%

We start by expressing the electronic Hamiltonian in the IP/EA basis to compare its expressions to the matrices of Sec.~\ref{sec:unfolded}.

\appendix

%=================================================================%
\section{Matrix elements of $C=[A, B]_{1,2}$}
%=================================================================%

An operator $A$ containing at most two-body terms may be written in normal ordered form with respect to the reference $\Phi$ as

$$
A=A_{0}+A_{1}+A_{2},
$$

where $A_{0}$ is a scalar, and

$$
\begin{gathered}
A_{1}=\sum_{p q} A_{p}^{q}\left\{\hat{a}_{q}^{p}\right\}, \\
A_{2}=\frac{1}{4} \sum_{p q r s} A_{p q}^{r s}\left\{\hat{a}_{r s}^{p q}\right\},
\end{gathered}
$$

with the second quantization operator written compactly as $\hat{a}_{q}^{p}=\hat{a}_{p}^{\dagger} \hat{a}_{q}$ and $\hat{a}_{r s}^{p q}=\hat{a}_{p}^{\dagger} \hat{a}_{q}^{\dagger} \hat{a}_{s} \hat{a}_{r}$. The commutator $C=[A, B]_{1,2}$ contains contributions from the following terms:

$$
\begin{gathered}
C_{0}=\left\langle\Phi\left|\left[A_{1}, B_{1}\right]\right| \Phi\right\rangle+\left\langle\Phi\left|\left[A_{2}, B_{2}\right]\right| \Phi\right\rangle, \\
C_{p}^{q}=\left[A_{1}, B_{1}\right]_{p}^{q}+\left[A_{1}, B_{2}\right]_{p}^{q}-\left[B_{1}, A_{2}\right]_{p}^{q}+\left[A_{2}, B_{2}\right]_{p}^{q}, \\
C_{p q}^{r s}=\left[A_{1}, B_{2}\right]_{p q}^{r s}-\left[B_{1}, A_{2}\right]_{p q}^{r s}+\left[A_{2}, B_{2}\right]_{p q}^{r s},
\end{gathered}
$$

where the unique contributions to the matrix elements are

$$
\begin{gathered}
\left\langle\Phi\left|\left[A_{1}, B_{1}\right]\right| \Phi\right\rangle=\sum_{p} \sum_{i}\left(A_{i}^{p} B_{p}^{i}-B_{i}^{p} A_{p}^{i}\right), \\
\left\langle\Phi\left|\left[A_{2}, B_{2}\right]\right| \Phi\right\rangle=\frac{1}{4} \sum_{i j} \sum_{a b}\left(A_{i j}^{a b} B_{a b}^{i j}-B_{i j}^{a b} A_{a b}^{i j}\right), \\
{\left[A_{1}, B_{1}\right]_{p}^{q}=\sum_{r}\left(A_{p}^{r} B_{r}^{q}-B_{p}^{r} A_{r}^{q}\right),} \\
{\left[A_{1}, B_{2}\right]_{p}^{q}=\sum_{i} \sum_{a} A_{i}^{a} B_{p a}^{q i}-A_{a}^{i} B_{p i}^{q a},} \\
{\left[A_{2}, B_{2}\right]_{p}^{q}=\frac{1}{2} \sum_{i j} \sum_{a}\left(A_{a p}^{i j} B_{i j}^{a q}-A_{i j}^{a q} B_{a p}^{i j}\right)} \\
+\frac{1}{2} \sum_{i} \sum_{a b}\left(A_{i p}^{a b} B_{a b}^{i q}-A_{a b}^{i q} B_{i p}^{a b}\right), \\
{\left[A_{1}, B_{2}\right]_{p q}^{r s}=\sum_{t}\left[P(p q) A_{p}^{t} B_{t q}^{r s}-P(r s) A_{t}^{r} B_{p q}^{t s}\right],} \\
{\left[A_{2}, B_{2}\right]_{p q}^{r s}=\frac{1}{2} \sum_{a b}\left(A_{p q}^{a b} B_{a b}^{r s}-A_{a b}^{r s} B_{p q}^{a b}\right)} \\
-\frac{1}{2} \sum_{i j}\left(A_{p q}^{i j} B_{i j}^{r s}-A_{i j}^{r s} B_{p q}^{i j}\right) \\
+\sum_{i} \sum_{a} P(p q) P(r s)\left[A_{p i}^{r a} B_{q a}^{s i}-A_{p a}^{r i} B_{q i}^{s a}\right] .
\end{gathered}
$$

In these equations $P(r s)$ is the antisymmetric permutation operator.

%=================================================================%
\section{The ph- and pp-RPA problems}
\label{sec:rpa}
%=================================================================%

%=================================================================%
\section{Perturbative matrix coefficients for $C^{(0)}$ diagonal}
\label{sec:diagC}
%=================================================================%

\begin{align}
  (\dv{\bV{}{(1)}}{s})_{pQ} &= (2 \bF{}{(0)}\bV{}{(1)}\bC{\text{d}}{(0)} - (\bF{}{(0)})^2\bV{}{(1)} - \bV{}{(1)}(\bC{\text{d}}{(0)})^2 )_{pQ}\\
                            &= \sum_{rS} 2 f^{(0)}_{pr} v^{(1)}_{rS}c^{(0)}_{SQ} - \sum_{rs} f^{(0)}_{pr} f^{(0)}_{rs} v^{(1)}_{sQ} - \sum_{RS} v^{(1)}_{pR} c^{(0)}_{RS}c^{(0)}_{SQ} \\
                            &=  \sum_{rS} 2 \epsilon^{(0)}_p\delta_{pr} v^{(1)}_{rS}\Delta\epsilon^{(0)}_Q\delta_{SQ} \\
                            &- \sum_{rs} \epsilon^{(0)}_p\delta_{pr} \epsilon^{(0)}_r\delta_{rs} v^{(1)}_{sQ} \\
                            &- \sum_{RS} v^{(1)}_{pR} \Delta\epsilon^{(0)}_R\delta_{RS}  \Delta\epsilon^{(0)}_Q\delta_{SQ} \\
                            &= (2 \epsilon^{(0)}_p\Delta\epsilon^{(0)}_Q - (\epsilon^{(0)}_p)^2 - (\Delta\epsilon^{(0)}_Q )^2) v^{(1)}_{pQ} \\
  \dv{v^{(1)}_{pQ}}{s}  &= - (\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_Q )^2 v^{(1)}_{pQ} \\
  &\color{red}{\boxed{\color{black}{v^{(1)}_{pQ}(s) = v^{(1)}_{pQ}(0) e^{-s(\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_Q )^2} }}} 
\end{align}
Note the close similarity with Evangelista's expressions for the off-diagonal part at first order!
\begin{align}
  (\dv{\bC{}{(1)}}{s})_{PQ} &= (2 \bC{\text{d}}{(0)}\bC{\text{od}}{(1)}\bC{\text{d}}{(0)}- (\bC{\text{d}}{(0)})^2\bC{\text{od}}{(1)} - \bC{\text{od}}{(1)}(\bC{\text{d}}{(0)})^2)_{PQ} \\
                            &= \sum_{RS} 2 c^{(0)}_{PR} c^{(1)}_{RS} c^{(0)}_{SQ} - c^{(0)}_{PR} c^{(0)}_{RS} c^{(1)}_{SQ} - c^{(1)}_{PR} c^{(0)}_{RS} c^{(0)}_{SQ} \\
                            &= 2 \Delta\epsilon^{(0)}_Pc^{(1)}_{PQ}\Delta\epsilon^{(0)}_Q - (\Delta\epsilon^{(0)}_P)^2  c^{(1)}_{PQ} -  c^{(1)}_{PQ} (\Delta\epsilon^{(0)}_Q)^2 \\
                            &= - (\Delta\epsilon^{(0)}_P - \Delta\epsilon^{(0)}_Q )^2 c^{(1)}_{PQ} \\
  &\color{red}{\boxed{\color{black}{c^{(1)}_{PQ}(s) = c^{(1)}_{PQ}(0) e^{-s(\Delta\epsilon^{(0)}_P - \Delta\epsilon^{(0)}_Q )^2} }}} 
\end{align}

\begin{align}
 &(\dv{\bF{}{(2)}}{s})_{pq}  = (\bF{}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bF{}{(0)} - 2 \bV{}{(1)}\bC{\text{d}}{(0)}\bV{}{(1),\dagger})_{pq}  \notag \\
                            &= \sum_{rS} f^{(0)}_{pr} v^{(1)}_{rS} v^{(1),\dagger}_{Sq} + \sum_{Rs}   v^{(1)}_{pR} v^{(1),\dagger}_{Rs} f^{(0)}_{sq} - 2\sum_{RS} v^{(1)}_{pR} c^{(0)}_{RS} v^{(1),\dagger}_{Sq} \notag \\
                            &= \sum_{S} \eps^{(0)}_{p} v^{(1)}_{pS} v^{(1)}_{qS} + \sum_{R} \eps^{(0)}_{q} v^{(1)}_{pR} v^{(1)}_{qR} - 2\sum_{R} \Delta\eps^{(0)}_R v^{(1)}_{pR} v^{(1)}_{qR} \notag \\
                            &= \sum_R (\eps^{(0)}_{p}  + \eps^{(0)}_{q} - 2 \Delta\eps^{(0)}_R) v^{(1)}_{pR} v^{(1)}_{qR} \notag \\
 &= \sum_R (\eps^{(0)}_{p}  + \eps^{(0)}_{q} - 2 \Delta\eps^{(0)}_R) v^{(1)}_{pR}(0) v^{(1)}_{qR}(0) e^{-s [ (\eps^{(0)}_p - \Delta\eps^{(0)}_R)^2+ (\eps^{(0)}_q - \Delta\eps^{(0)}_R)^2]} \notag \\
 &f^{(2)}_{pq}(s) = \notag \\
 &\color{red}{\boxed{\color{black}{- \sum_R\frac{ v^{(1)}_{pR}(0) v^{(1)}_{qR}(0)  (\eps^{(0)}_{p}  + \eps^{(0)}_{q} - 2 \Delta\eps^{(0)}_R)}{(\eps^{(0)}_p - \Delta\eps^{(0)}_R)^2+ (\eps^{(0)}_q - \Delta\eps^{(0)}_R)^2}(1 - e^{-s [ (\eps^{(0)}_p - \Delta\eps^{(0)}_R)^2+ (\eps^{(0)}_q - \Delta\eps^{(0)}_R)^2]})}}} \notag
\end{align}

\begin{align}
  (\dv{\bV{}{(2)}}{s})_{pQ} &=  (2 \bF{}{(0)}\bV{}{(2)}\bC{\text{d}}{(0)} - (\bF{}{(0)})^2\bV{}{(2)} - \bV{}{(2)}(\bC{\text{d}}{(0)})^2 \\
  & - 2 \bV{}{(1)} \bC{\text{d}}{(0)} \bC{\text{od}}{(1)}  + \bF{}{(0)} \bV{}{(1)} \bC{\text{od}}{(1)} + \bV{}{(1)} \bC{\text{od}}{(1)} \bC{\text{d}}{(0)})_{pQ} \notag \\
  v^{(2)}_{pQ}(s) &= v^{(2)}_{pQ}(0) e^{-s(\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_Q )^2} + \text{Non-homogeneous solution} \notag \\
  v^{(2)}_{pQ}(s) &= \text{Non-homogeneous solution}
\end{align}

\end{document}