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

%============================================================%
%%% NEWCOMMANDS %%%
% ============================================================%

\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{\ant}[1]{\textcolor{green}{#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 $\lambda$ is such that
\begin{equation}
  \bH(0) = E_0(0) + \bF{}{}(0) + \lambda \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) + \lambda \bV{}{\mathrm{d}}(0) & \bH^\text{od}(0)&= \lambda \bV{}{\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}
  \label{eq:SRG_MP2}
  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}
\PFL{Not quite. You're missing the eigenvectors to make it a non-linear eigenvalue problem.}
However, in order to use it we need to rely on approximations of the dynamical 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{equation}
  \label{eq:GF2_selfenergy}
  \Sigma_{pq}^{\text{GF(2)}}(\omega) 
  = \sum_{ija} \frac{W_{pa,ij}^{\text{GF(2)}}W_{qa,ij}^{\text{GF(2)}}}{\omega + \epsilon _a -\epsilon_i -\epsilon_j - \ii \eta} 
  + \sum_{iab} \frac{W_{pi,ab}^{\text{GF(2)}}W_{qi,ab}^{\text{GF(2)}}}{\omega + \epsilon _i -\epsilon_a -\epsilon_b + \ii \eta} 
\end{equation}
with
\begin{equation}
  \label{eq:GF2_sERI}
  W^{\GF}_{pq,rs}= \frac{1}{\sqrt{2}}\aeri{pq}{rs} 
\end{equation}
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}
  \Sigma_{pq}^{\GW}(\omega) 
  = \sum_{iv} \frac{W_{pi,v}^{\GW} W_{qi,v}^{\GW}}{\omega - \epsilon_i + \Omega_{v}^{\dRPA} - \ii \eta} 
  + \sum_{av} \frac{W_{pa,v}^{\GW}W_{qa,v}^{\GW}}{\omega - \epsilon_a - \Omega_{v}^{\dRPA} + \ii \eta} \notag
\end{equation}
with
\begin{equation}
  \label{eq:GW_sERI}
  W_{pq,v}^{\GW} = \sum_{ia}\eri{pi}{qa}\qty( \bX_{v}^{\dRPA} + \bY_{v}^{\dRPA} )_{ia}
\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}
  \Sigma_{pq}^{\GT}(\omega) 
  = \sum_{iv} \frac{W_{pi,v}^{N+2} W_{qi,v}^{N+2}}{\omega + \epsilon_i - \Omega_{v}^{N+2} + \ii \eta} 
  + \sum_{av} \frac{W_{pa,v}^{N-2} W_{qa,v}^{N-2}}{\omega + \epsilon_a - \Omega_{v}^{N-2} - \ii \eta} \notag
\end{equation}
with
\begin{align}
	\label{eq:GT_sERI}
  W_{pq,v}^{N+2} & = \sum_{c<d} \aeri{pq}{cd} \bX_{cd,m}^{N+2} + \sum_{k<l} \aeri{pq}{kl} \bY_{kl,m}^{N+2} \notag \\
  W_{pq,v}^{N-2} & = \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 approximation 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}}	\\
    (\bV{}{\text{2h1p}})^{\mathsf{T}}	&	\bC{}{\text{2h1p}}		&	\bO				\\
    (\bV{}{\text{2p1h}})^{\mathsf{T}}	&	\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}
    \label{eq:GF2_unfolded}
    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( \epsilon_i + \epsilon_j - \epsilon_a) \delta_{jl} \delta_{ac} \delta_{ik} 
    \\
    C^\text{2p1h}_{iab,kcd} & = \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} \delta_{bd} 
  \end{align}
\item \textbf{GW}
  \begin{align}
    \label{eq:GW_unfolded}
    V^\text{2h1p}_{p,klc} & = \eri{pc}{kl}
    \\
    V^\text{2p1h}_{p,kcd} & = \eri{pk}{dc}
    \\
    C^\text{2h1p}_{ija,klc} & = \qty[ \qty( \epsilon_i + \epsilon_j - \epsilon_a) \delta_{jl} \delta_{ac} - \eri{jc}{al} ] \delta_{ik}
    \\
    C^\text{2p1h}_{iab,kcd} & = \qty[ \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} + \eri{ak}{ic} ] \delta_{bd} 
  \end{align}
\item \textbf{GT}
  \begin{align}
    \label{eq:GT_unfolded}
    V^\text{2h1p}_{p,klc} & = \aeri{pc}{kl}
    \\
    V^\text{2p1h}_{p,kcd} & = \aeri{pk}{cd}
    \\
    C^\text{2h1p}_{ija,klc} & = \qty[ \qty( \epsilon_i + \epsilon_j - \epsilon_a) \delta_{jl} \delta_{ac} - \aeri{ij}{kl} ] \delta_{ac} 
    \\
    C^\text{2p1h}_{iab,kcd} & = \qty[ \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} + \aeri{ab}{cd} ] \delta_{ik} 
  \end{align}
\end{itemize}
The downfolding procedure to obtain the $GW$ self-energy is derived in details in Appendix~\ref{sec:downfolding}.

\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{Order by order differential equations}
%%%%%%%%%%%%%%%%%%%%%%

%///////////////////////////%
\subsubsection{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}
  \bH =
         \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}
  \bH(0) =
  \begin{pmatrix}
    \bF{}{} & \bO \\
    \bO & \bC{}{} 
  \end{pmatrix}
  + \lambda
  \begin{pmatrix}
    \bO & \bV{}{} \\
    \bV{}{\dagger} & \bO
  \end{pmatrix}
\end{equation}
we could also have defined it like this
\begin{equation}
  \label{eq:H_MBPT_partitioning}
  \bH(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}
However, in the end this alternative choice was not judicious and the reason is explained why in Appendix~\ref{sec:partitioning}.

The initial conditions corresponding to the first partitioning are
\begin{align}
  \bHd{0}(0) &=   \begin{pmatrix}
                    \bF{}{} & \bO \\
                    \bO & \bC{}{} 
                  \end{pmatrix}  &
                                   \bHod{0}(0) &=\begin{pmatrix}
                                                   \bO & \bO \\
                                                   \bO & \bO 
                                                 \end{pmatrix} \notag \\
  \bHd{1}(0) &=\begin{pmatrix}
                 \bO & \bO \\
                 \bO & \bO 
               \end{pmatrix} &
                               \bHod{1}(0) &=   \begin{pmatrix}
                                                  \bO & \bV{}{} \\
                                                  \bV{}{\dagger} & \bO \notag
                                                \end{pmatrix}   \notag 
\end{align}


%///////////////////////////%
\subsubsection{Zeroth order differential equations}
%///////////////////////////%

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}

%///////////////////////////%
\subsubsection{First order differential equations}
%///////////////////////////%

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{}{(0)}\\
    \bC{}{(0)}\bV{}{(1),\dagger} - \bV{}{(1),\dagger} \bF{}{(0)} & \bO \notag
  \end{pmatrix}  
\end{align}

\begin{align}
  \dv{\bH^{(1)}}{s} &= \comm{\bEta{1}}{\bHd{0}}  \\
  \dv{\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
\end{align}

The last two equations can be solved differently depending on the form of $\bF{}{}$ and $\bC{}{}$.

%///////////////////////////%
\subsubsection{Second order differential equations}
%///////////////////////////%

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{2}}{\bHod{0}} \\
  &= \comm{\bHd{0}}{\bHod{2}} \notag \\
  &= \begin{pmatrix}
    \bO & \bF{}{(0)}\bV{}{(2)} - \bV{}{(2)}\bC{}{(0)}\\
    \bC{}{(0)}\bV{}{(2),\dagger} - \bV{}{(2),\dagger}\bF{}{(0)} & \bO \notag
  \end{pmatrix}  
\end{align}

\begin{align}
  \dv{\bH^{(2)}}{s} &= \comm{\bEta{2}}{\bH^{(0)}} + \comm{\bEta{1}}{\bH^{(1)}} + \comm{\bEta{0}}{\bH^{(2)}}  \\
  &= \comm{\bEta{2}}{\bHd{0}} + \comm{\bEta{1}}{\bHod{1}} \notag
\end{align}

\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),\dagger}\bV{}{(1)} + \bV{}{(1),\dagger}\bV{}{(1),}\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}

%///////////////////////////%
\subsubsection{Third order differential equations}
% ///////////////////////////%

\begin{equation}
  \bEta{3} = \comm{\bHd{0}}{\bHod{3}} + \comm{\bHd{2}}{\bHod{1}} 
\end{equation}

We will show after that $\bV{}{(2)} = \bO$, and hence $\bEta{2} = 0$, but we already use this result to simplify derivations

\begin{equation}
  \dv{\bH^{(3)}}{s} = \comm{\bEta{3}}{\bHd{0}}  + \comm{\bEta{1}}{\bHd{2}} 
\end{equation}

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

%///////////////////////////%
\subsubsection{Forth order differential equations}
% ///////////////////////////%

\begin{equation}
  \bEta{4} = \comm{\bHd{0}}{\bHod{4}} + \comm{\bHd{2}}{\bHod{2}} + \comm{\bHd{3}}{\bHod{1}} 
\end{equation}

\begin{equation}
  \dv{\bH^{(4)}}{s} = \comm{\bEta{4}}{\bHd{0}} + \comm{\bEta{2}}{\bHd{2}} + \comm{\bEta{1}}{\bHd{3}} 
\end{equation}

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

%///////////////////////////%
\subsubsection{First order Hamiltonian elements}
%///////////////////////////%

\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 $\bC{}{\GF}(0)$ is diagonal.
Inspired by this fact we transform $\bC{}{(0)}$ to its diagonal representation using
$\bC{}{(0)} = \bU \bD^{(0)} \bU^{-1}$ and insert it in the differential equation for $\bV{}{(1)}$
\begin{align}
  \dv{\bV{}{(1)}}{s} &= 2 \bF{}{(0)}\bV{}{(1)}\bU \bD^{(0)} \bU^{-1}- (\bF{}{(0)})^2\bV{}{(1)} - \bV{}{(1)}\bU (\bD^{(0)})^2 \bU^{-1}\\
  \dv{\bV{}{(1)}}{s}\bU &= 2 \bF{}{(0)}\bV{}{(1)}\bU \bD^{(0)} - (\bF{}{(0)})^2\bV{}{(1)} \bU - \bV{}{(1)}\bU (\bD^{(0)})^2 \\
  \dv{\bW^{(1)}}{s} &= 2 \bF{}{(0)}\bW^{(1)} \bD^{(0)} - (\bF{}{(0)})^2\bW^{(1)} - \bW^{(1)} (\bD^{(0)})^2
\end{align}
where in the last line we have defined the matrix of screened integral $\bW^{(1)}=\bV{}{(1)} \bU$.
Note that in the GF(2) case $\bU = \mathbb{1}$ and thus $\bW^{(1)}=\bV{}{(1)}$.
\begin{align}
  &(\dv{\bW^{(1)}}{s})_{p,(q,v)} = (2 \bF{}{(0)}\bW^{(1)} \bD^{(0)} - (\bF{}{(0)})^2\bW^{(1)} - \bW^{(1)} (\bD^{(0)})^2)_{p,(q,v)} \notag \\
  &= \sum_{r,(s,x)} 2 F_{pr}^{(0)}W_{r,(s,x)}^{(1)}D_{(s,x),(q,v)}^{(0)} - \sum_{r,s} F_{pr}^{(0)} F_{rs}^{(0)}W_{s,(q,v)}^{(1)} \notag \\
  &- \sum_{(r,x),(s,y)} W_{p,(r,x)}^{(1)}D_{(r,x),(s,y)}^{(0) } D_{(s,y),(q,v)}^{(0)} \notag \\
  &= 2 F_{pp}^{(0)}W_{p,(q,v)}^{(1)}D_{(q,v),(q,v)}^{(0)} - (F_{pp}^{(0)})^2 W_{p,(q,v)}^{(1)} \notag - W_{p,(q,v)}^{(1)} (D_{(q,v),(q,v)}^{(0)})^2 \\
  &(\dv{\bW^{(1)}}{s})_{p,(q,v)} = - (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 W_{p,(q,v)}^{(1)} 
\end{align}
This equation can be integrated to give
\begin{equation}
  \color{red}{\boxed{\color{black}{W_{p,(q,v)}^{(1)}(s) = W_{p,(q,v)}^{(1)}(0) e^{- (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 s}}}} 
\end{equation}

%///////////////////////////%
\subsubsection{Second order Hamiltonian elements}
%///////////////////////////%

\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. 

\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} \\
                     &= \bF{}{(0)}\bW^{(1)}\bW^{(1),\dagger} + \bW^{(1)}\bW^{(1),\dagger}\bF{}{(0)} - 2 \bW^{(1)}\bD^{(0)}\bW^{(1),\dagger} \notag
\end{align}

\begin{align}
  \dv{F_{pq}^{(2)}}{s} &= F_{pp}^{(0)}(\bW^{(1)}\bW^{(1),\dagger} )_{pq} + (\bW^{(1)}\bW^{(1),\dagger})_{pq}F_{qq}^{(0)}\\
                       &- 2 (\bW^{(1)}\bC{}{(0)}\bW^{(1),\dagger})_{pq}
\end{align}

\begin{itemize}
\item \textbf{GW}
  
  In the GW case this evaluates to
\begin{align}
  &\dv{F_{pq}^{(2)}}{s} = \sum_{r,v}(\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)})W^{(1)}_{p,(r,v)}W^{(1),\dagger}_{(r,v),q} - 2 W^{(1)}_{p,(r,v)}C^{(0)}_{(r,v),(r,v)}W^{(1),\dagger}_{(r,v),q} \notag \\
  &= \sum_{r,v}(\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2 C^{(0)}_{(r,v),(r,v)})W^{(1)}_{p,(r,v)}W^{(1),\dagger}_{(r,v),q} \notag \\
  &= \sum_{r,v}(\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2 (\epsilon_r^{(0)} \pm \Omega_v))W^{(0)}_{p,(r,v)} \notag \\
  &\times W^{(0),\dagger}_{(r,v),q} e^{-(\epsilon_p - \epsilon_r \pm \Omega_v)^2s}   e^{-(\epsilon_q - \epsilon_r \pm \Omega_v)^2s}  \notag
\end{align}

which can be integrated as
\begin{align}
  F_{pq}^{(2)}(s) &= -\sum_{r,v} \frac{\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2 (\epsilon_r^{(0)} \pm \Omega_v)}{(\epsilon_p^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2 + (\epsilon_q^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2} W^{(0)}_{p,(r,v)} \notag \\
  &\times W^{(0),\dagger}_{(r,v),q} e^{-(\epsilon_p^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2s}   e^{-(\epsilon_q^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2s}  + \text{Cte} \notag
\end{align}

The constant is determined as
\begin{align}
  &F_{pq}^{(2)}(0) = 0 \notag \\
  &= - \sum_{r,v} \frac{\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2 (\epsilon_r^{(0)} \pm \Omega_v)}{(\epsilon_p^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2 + (\epsilon_q^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2} W^{(0)}_{p,(r,v)} W^{(0),\dagger}_{(r,v),q} + \text{Cte} \notag
\end{align}
Which finally gives
\begin{align}
  F_{pq}^{(2)}(s) &= \sum_{r,v} \frac{\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2 (\epsilon_r^{(0)} \pm \Omega_v)}{(\epsilon_p^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2 + (\epsilon_q^{(0)} - \epsilon_r^{(0)} \pm \Omega_v)^2} W^{(0)}_{p,(r,v)} \notag \\
                  &\times W^{(0),\dagger}_{(r,v),q}\left(1 - e^{-(\epsilon_p - \epsilon_r \pm \Omega_v)^2s}   e^{-(\epsilon_q - \epsilon_r \pm \Omega_v)^2s}\right) \notag
\end{align}

\item \textbf{GF(2)}

  The expression for the GF(2) case is
\begin{align}
  F_{pq}^{(2)}(s) &= \frac{1}{2} \sum_{ija} \frac{\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2\Delta_{ij}^a}{(\epsilon_p^{(0)}  - \Delta_{ij}^a)^2 + (\epsilon_q^{(0)}  - \Delta_{ij}^a)^2} \aeri{pa}{ij}\notag \\
                  &\times\aeri{qa}{ij}  \left(1 - e^{-(\epsilon_p - \Delta_{ij}^a)^2s}   e^{-(\epsilon_q - \Delta_{ij}^a)^2s}\right) \notag \\
  & + \frac{1}{2} \sum_{iab} \frac{\epsilon_{p}^{(0)} + \epsilon_{q}^{(0)} - 2\Delta_{i}^{ab}}{(\epsilon_p^{(0)}  - \Delta_{ij}^a)^2 + (\epsilon_q^{(0)}  - \Delta_{i}^{ab})^2} \aeri{pi}{ab}\notag \\
                  &\times\aeri{qi}{ab}  \left(1 - e^{-(\epsilon_p - \Delta_{i}^{ab})^2s}   e^{-(\epsilon_q - \Delta_{i}^{ab})^2s}\right) \notag
\end{align}

\end{itemize}

The equation for $\bV{}{(2)}$ and $\bV{}{(2),\dagger}$ are the same that for their first order counterpart.
So we define $\bW^{(2)} = \bV{}{(2)} \bU$ and we obtain
\begin{equation}
  W_{p,(q,v)}^{(2)}(s) = W_{p,(q,v)}^{(2)}(0) e^{- (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 s} = 0
\end{equation}
Therefore, we have
\begin{align}
  &\color{red}{\boxed{\color{black}{\bV{}{(2)}(s)= \bO}}} & &\color{red}{\boxed{\color{black}{\bV{}{(2),\dagger}(s)= \bO}}} 
\end{align}

Finally, the elements of $\bC{}{(2)}$ can be obtained exactly like their $\bF{}{(2)}$ counterpart just by swapping the role of $\bF{}{(0)}$ and $\bC{}{(0)}$ in the formula.

%///////////////////////////%
\subsubsection{Third order Hamiltonian elements}
%///////////////////////////%

\begin{align}
  \dv{\bF{}{(3)}}{s}  &= \bO  \Longleftrightarrow \bF{}{(3)}(s) = \bF{}{(3)}(0) \Longleftrightarrow \color{red}{\boxed{\color{black}{\bF{}{(3)}(s)= \bO}}} \\
  \dv{\bC{}{(3)}}{s}  &= \bO  \Longleftrightarrow \bC{}{(1)}(s) = \bC{}{(1)}(0) \Longleftrightarrow \color{red}{\boxed{\color{black}{\bC{}{(3)}(s)= \bO}}} 
\end{align}

The homogeneous equations for $\bV{}{(3)}$ and $\bV{}{(3), \dagger}$ give homogeneous solutions equal to zero for the same reasons as $\bV{}{(2)}$ and $\bV{}{(2), \dagger}$.
Therefore the differential equations can be simplified to
\begin{align}
   \dv{\bV{}{(3)}}{s} &=  2 \bF{}{(0)}\bV{}{(1)}\bC{}{(2)} + 2 \bF{}{(2)}\bV{}{(1)}\bC{}{(0)} \notag \\
                    &- \left( \bF{}{(0)}\bF{}{(2)} + \bF{}{(2)}\bF{}{(0)} \right)\bV{}{(1)} \notag \\
                    &- \bV{}{(1)}\left( \bC{}{(0)}\bC{}{(2)} + \bC{}{(2)}\bC{}{(0)} \right)  \notag \\
  \dv{\bV{}{(3),\dagger}}{s} &= 2 \bC{}{(0)}\bV{}{(1),\dagger}\bF{}{(2)} + 2 \bC{}{(2)}\bV{}{(1),\dagger}\bF{}{(0)} \notag \\
                    &- \bV{}{(1),\dagger}\left( \bF{}{(0)}\bF{}{(2)} + \bF{}{(2)}\bF{}{(0)} \right) \notag \\
                    &-\left( \bC{}{(0)}\bC{}{(2)} + \bC{}{(2)}\bC{}{(0)} \right)\bV{}{(1),\dagger}  \notag 
\end{align}
and can be solved by integration of the right side using the previous analytic formula for Hamiltonian elements.

\begin{align}
  \text{Part B}: &(\dv{\bW^{(3)}}{s})_{pR} = - \left[\left( \bF{}{(0)}\bF{}{(2)} + \bF{}{(2)}\bF{}{(0)} \right)\bW^{(1)}\right]_{pR} \notag \\
                 &= - \sum_q\left( \bF{}{(0)}\bF{}{(2)} + \bF{}{(2)}\bF{}{(0)} \right)_{pq}W^{(1)}_{qR} \notag \\
                 &= - \sum_q(\epsilon_p+\epsilon_q)F_{pq}^{(2)}W^{(1)}_{qR} \notag \\
                 &= - \sum_{qS} (\epsilon_p+\epsilon_q) \frac{\Delta_{pS}+\Delta_{qS}}{\Delta_{pS}^2+\Delta_{qS}^2}W^{(1)}_{pS}(0)W^{(1)}_{qS}(0) W^{(1)}_{qR}(0) \notag \\
                 &\times (1 - e^{-(\Delta_{pS}^2+\Delta_{qS}^2)s})e^{-\Delta_{qR}^2s}   \notag \\
  \text{Part B}: & (\bW^{(3)}(s))_{pR} = - \sum_{qS} (\epsilon_p+\epsilon_q) \frac{\Delta_{pS}+\Delta_{qS}}{\Delta_{pS}^2+\Delta_{qS}^2}W^{(1)}_{pS}(0)W^{(1)}_{qS}(0) W^{(1)}_{qR}(0) \notag \\
  &\times \left(-\frac{e^{-\Delta_{qR}^2s} }{\Delta_{qR}^2} + \frac{e^{-(\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2)s}}{\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2}\right)
\end{align}

\begin{align}
  \text{Part C}: &(\dv{\bW^{(3)}}{s})_{pR} = - \left[\bW^{(1)}\left( \bC{}{(0)}\bC{}{(2)} + \bC{}{(2)}\bC{}{(0)} \right)\right]_{pR} \notag \\
                 &= - \sum_S W^{(1)}_{pS} \left( \bD^{(0)}\bD^{(2)} + \bD^{(2)}\bD^{(0)} \right)_{SR}\notag \\
                 &= - \sum_SW^{(1)}_{pS}(D_R+D_S)D_{SR}^{(2)}\notag \\
                 &= - \sum_{Sq}W^{(1)}_{pS}(0)e^{-\Delta_{pS}^2s} (D_R+D_S)\frac{-\Delta_{qS}-\Delta_{qR}}{\Delta_{qS}^2+\Delta_{qR}^2} \notag \\
                 &\times W^{(1)}_{qS}(0)W^{(1)}_{qR}(0)(1-e^{-(\Delta_{qS}^2+\Delta_{qR}^2)s})) \notag \\
  \text{Part C}: & (\bW^{(3)}(s))_{pR} = + \sum_{Sq} (D_R+D_S) \frac{\Delta_{qS}+\Delta_{qR}}{\Delta_{qS}^2+\Delta_{qR}^2} W^{(1)}_{pS}(0)W^{(1)}_{qS}(0) W^{(1)}_{qR}(0) \notag \\
  &\times \left(-\frac{e^{-\Delta_{pS}^2s} }{\Delta_{pS}^2} + \frac{e^{-(\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2)s}}{\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2}\right) \notag 
\end{align}

\begin{align}
  \text{Part A}: &(\dv{\bW^{(3)}}{s})_{pR}  = (2 \bF{}{(0)}\bW^{(1)}\bD^{(2)} + 2 \bF{}{(2)}\bW^{(1)}\bD^{(0)})_{pR} \notag \\
                 &= 2\sum_{qS} \epsilon_p\delta_{pq} W_{qS}^{(1)}D_{SR}^{(2)} + F_{pq}^{(2)}W_{qS}^{(1)}D_S\delta_{SR} \notag \\
                 &= 2\sum_S \epsilon_pW_{pS}^{(1)}(0) e^{-\Delta_{pS}^2 s} D_{SR}^{(2)} \notag \\
                 &+ 2\sum_q D_S W_{qR}^{(1)}(0) e^{-\Delta_{qR}^2 s} F_{pq}^{(2)} \notag \\
                 &=2\sum_{qS} \epsilon_pW_{pS}^{(1)}(0) e^{-\Delta_{pS}^2 s} \frac{-\Delta_{qS}-\Delta_{qR}}{\Delta_{qS}^2+\Delta_{qR}^2} \notag \\
                 &\times W^{(1)}_{qS}(0)W^{(1)}_{qR}(0)(1-e^{-(\Delta_{qS}^2+\Delta_{qR}^2)s})) \notag \\
                 &+ 2\sum_{qS} D_S W_{qR}^{(1)}(0) e^{-\Delta_{qR}^2 s} \frac{\Delta_{pS}+\Delta_{qS}}{\Delta_{pS}^2+\Delta_{qS}^2}\notag \\
                 &\times W^{(1)}_{pS}(0)W^{(1)}_{qS}(0)  (1 - e^{-(\Delta_{pS}^2+\Delta_{qS}^2)s}) \notag
\end{align}

\begin{align}
  & (\bW^{(3)}(s))_{pR} = \sum_{qS} (2D_S - \epsilon_p - \epsilon_q) \frac{\Delta_{pS}+\Delta_{qS}}{\Delta_{pS}^2+\Delta_{qS}^2}W^{(1)}_{pS}(0)W^{(1)}_{qS}(0) W^{(1)}_{qR}(0) \notag \\
  &\times \left(-\frac{e^{-\Delta_{qR}^2s} }{\Delta_{qR}^2} + \frac{e^{-(\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2)s}}{\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2}\right) \notag \\
  &+ \sum_{Sq} (D_R+D_S - 2\epsilon_p) \frac{\Delta_{qS}+\Delta_{qR}}{\Delta_{qS}^2+\Delta_{qR}^2} W^{(1)}_{pS}(0)W^{(1)}_{qS}(0) W^{(1)}_{qR}(0) \notag \\
  &\times \left(-\frac{e^{-\Delta_{pS}^2s} }{\Delta_{pS}^2} + \frac{e^{-(\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2)s}}{\Delta_{pS}^2+\Delta_{qS}^2+\Delta_{qR}^2}\right) \notag 
\end{align}

%///////////////////////////%
\subsubsection{Forth order Hamiltonian elements}
%///////////////////////////%

%%%%%%%%%%%%%%%%%%%%%%
\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)} + \bF{}{(2)}(s) + \bF{}{(4)}(s) & \bV{}{(1)}(s)  + \bV{}{(3)}(s) \\
    \bV{}{(1),\dagger}(s)  + \bV{}{(3),\dagger}(s) & \bC{}{(0)} +\bC{}{(2)}(s) + \bC{}{(4)}(s)
  \end{pmatrix}
\end{equation}

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

  Then the quasiparticle equation truncated to zeroth, second and forth order is
  \begin{align}
    \left[ \bF{}{(0)} \right] \bR^{1h/1p} &= \omega \bR^{1h/1p}    \\
    \left[ \bF{}{(0)} + \bF{}{(2)}(s) + \bV{}{(1)}(s) (\omega \mathbb{1} - \bC{}{(0)} )^{-1}  \bV{}{(1),\dagger}(s)  \right] \bR^{1h/1p} &= \omega \bR^{1h/1p}    \\
    \left[ (\bF{}{(0)} + \bF{}{(2)}(s) + \bF{}{(4)}(s)) + (\bV{}{(1)}(s) + \bV{}{(3)}(s)) (\omega \mathbb{1} - \bC{}{(0)} - \bC{}{(2)}(s))^{-1}  (\bV{}{(1),\dagger}(s) + \bV{}{(3),\dagger}(s)) \right] \bR^{1h/1p} &= \omega \bR^{1h/1p}   
  \end{align}
  
\end{widetext}

To zeroth order in the coupling we recover the HF quasiparticle energies which makes sense.
Now turning to the first non-zero correction to the quasiparticle equation which is the second order equation

\begin{align}
  \bSig_c^{(2)} (\omega) &= \bV{}{\hhp,(1)} \bU^{\hhp} \text{diag}(\frac{1}{\omega - D_{(i,v),(i,v)}}) (\bU^{\hhp})^{-1} (\bV{}{\hhp,(1)})^{\mathsf{T}} \notag \\
  &+ \bV{}{\pph,(1)} \bU^{\pph} \text{diag}(\frac{1}{\omega - D_{(a,v),(a,v)}}) (\bU^{\pph})^{-1} (\bV{}{\pph,(1)})^{\mathsf{T}} \notag \\
                        &= \bW^{\hhp,(1)}  \text{diag}(\frac{1}{\omega - D_{(i,v),(i,v)}}) (\bW^{\hhp,(1)})^{\mathsf{T}} \notag \\
  &+ \bW^{\pph,(1)}  \text{diag}(\frac{1}{\omega - D_{(a,v),(a,v)}}) (\bW^{\pph,(1)})^{\mathsf{T}} \notag 
 \end{align}

 \begin{itemize}

 \item \textbf{GF(2)}

   In the GF(2) case we have $D_{(i,v),(i,v)} = D_{ija,ija} = \epsilon_i + \epsilon_j  - \epsilon_a$ and $D_{(a,v),(a,v)} = D_{iab,iab} = \epsilon_a + \epsilon_b - \epsilon_i $ and the $W$ matrix elements have been defined in Eq.~(\ref{eq:GF2_sERI}).
   \begin{align}
     (\bSig_c^{(2)} (\omega,s))_{pq} &= \sum_{ija} W_{p,i[ja]}^{(1)}   \frac{1}{\omega - D_{ija,ija}}(W^{\mathsf{T}})_{i[ja],q}^{(1)} \notag \\
                                  &+ \sum_{iab} W_{p,[ia]b}^{(1)}  \frac{1}{\omega - D_{iab,iab}}(W^{\mathsf{T}})_{[ia]b,q}^{(1)} \notag \\
                                  &= \sum_{ija}  \frac{W_{pa,ij}^{(1)}  W_{qa,ij}^{(1)} }{\omega - (\epsilon_i + \epsilon_j  - \epsilon_a)} \notag \\
                                   &+ \sum_{iab} \frac{W_{pi,ab}^{(1)}  W_{qi,ab}^{(1)}}{\omega - (\epsilon_a + \epsilon_b  - \epsilon_i)} \notag \\
                                   &= \sum_{ija}  \frac{W_{pa,ij}^{(1)}(0)  W_{qa,ij}^{(1)}(0) }{\omega - (\epsilon_i + \epsilon_j  - \epsilon_a)}  e^{-(\epsilon_p - \Delta_{ij}^a)^2s}e^{-(\epsilon_q - \Delta_{ij}^a)^2s} \notag \\
     &+ \sum_{iab}  \frac{W_{pi,ab}^{(1)}(0)  W_{qi,ab}^{(1)}(0) }{\omega - (\epsilon_a + \epsilon_b  - \epsilon_i)}  e^{-(\epsilon_p - \Delta_{i}^{ab})^2s}e^{-(\epsilon_q - \Delta_{i}^{ab})^2s} \notag
 \end{align}

 \item \textbf{GW}

   A similar derivation gives
   
   \begin{align}
     \label{eq:SRGGW_selfenergy}
     (\Sigma_c^{\GW}(\omega,s))_{pq} &= \sum_{iv} \frac{W_{pi,v}^{(1)}(0)  W_{qi,v}^{(1)}(0)}{\omega - \epsilon_i + \Omega_{v}^{\dRPA} - \ii \eta}e^{-(\epsilon_p - \epsilon_i + \Omega_v)^2s}   e^{-(\epsilon_q - \epsilon_i + \Omega_v)^2s}   \notag \\
     &+ \sum_{av} \frac{W_{pa,v}^{(1)}(0) W_{qa,v}^{(1)}(0)}{\omega - \epsilon_a - \Omega_{v}^{\dRPA} + \ii \eta} e^{-(\epsilon_p - \epsilon_a - \Omega_v)^2s}e^{-(\epsilon_q - \epsilon_a - \Omega_v)^2s}  \notag
 \end{align}

 \item \textbf{GT}
   
 \textcolor{red}{\textbf{TODO Give analytical expression for the GT case.}}

\end{itemize}

Now that we have analytical expression for the second order correlation part of the self-energy as well as for $\bF{}{(2)}$, we have every ingredients for the second order quasi-particle equation.

In the previous formula we can see that the diagonal elements at $s \to \infty$ correspond to the same values as in the usual diagonal static approximation.
Therefore, the SRG as a true renormalization group method removes the coupling $\bV{}{}$ while incorporating some of its physics in the non-coupled problem $\bF{}{}$.
This formalism gives us a rationalization of the diagonal static approximation from a RG perspective.

In addition, this gives us an unambiguous way to define a non-diagonal static approximation.
On the other hand the usual non-diagonal static self-energy used in qsGW, 
\begin{equation}
  (\Sigma_c^{qsGW})_{pq} = \frac{\Sigma_c(\epsilon_p)_{pq} + \Sigma_c(\epsilon_q)_{qp}}{2}
\end{equation}
is not unambiguously defined because this is not the only possible choice starting from the $GW$ self-energy.
If we define $x=\epsilon_p^{(0)} - \epsilon_r^{(0)} \pm \Omega_v$ and $y=\epsilon_q^{(0)} - \epsilon_r^{(0)} \pm \Omega_v$ then the SRG qsGW is of the form $(x+y)/(x^2 + y^2)$ while the usual qsGW is of the form $(x+y)/2xy$.
Note that both diagonal are of the form $1/x$ which is consistent with the fact that the SRG diagonal correspond to the usual static diagonal.

Even more, the SRG formalism defines a hierarchy of static approximation by considering higher and higher perturbation order for $\bSig$.
This hierarchy could be compared to another hierarchy of static approximation obtained by perturbing the static self-energy by its difference to its dynamic counterpart.

One of the con of the static approximation is that we loose information about the satellites and this is also true for the SRG also when the coupling has been totally removed.
However, SRG allows us to stop at a finite value $s$ corresponding to a renormalized coupling but the coupling is still present.
Therefore the satellites can still be observed due to the non-linearity induced by the renormalized coupling.

%=================================================================%
\section{Discussion on the choice of regularizers}
\label{sec:regularizers}
%=================================================================%

Before discussing the regularizers in MBPT, we start by analyzing behavior of the second order SRG correction to the energy derived in Eq.~(\ref{eq:SRG_MP2}). The equation is written here again for readability and also we used $s = 1/\Lambda^2$ where $\Lambda$ is the characteristic cutoff energy.
\begin{equation}
  E_0^{(2)}(s= 1/\Lambda^2) = \frac{1}{4} \sum_{i j} \sum_{a b} \frac{\aeri{ij}{ab}^2}{\Delta_{ab}^{ij}}\left(1-e^{-2\left(\frac{\Delta_{ab}^{ij}}{\Lambda}\right)^2}\right)
\end{equation}
For $\Lambda \to 0$ the SRG-MP2 energy is equal to the MP2 one while  for $\Lambda  \to \infty$ the SRG-MP2 energy goes to zero.
For a finite value of $s$, hence a finite energy cutoff $\Lambda$, then the term of the sum with $\Delta_{ab}^{ij} < \Lambda$ are almost zero. Therefore a small cutoff removes only the divergent $1/\Delta_{ab}^{ij}$.

A similar analysis can be done about the regularized correlation self-energy introduced by Monino and Loos. Here we discuss only the GW self-energy but without loss of generality.
\begin{align}
  \label{eq:GW_selfenergy_regularized}
  (\Sigma_c^{\GW}(\omega,\Lambda))_{pq} &= \sum_{iv} \frac{W_{pi,v}^{\GW} W_{qi,v}^{\GW}}{\omega - \Omega_{i,v}^{\dRPA} - \ii \eta}\left(1-e^{-2\left( \frac{\omega - \Omega_{i,v}^{\dRPA} - \ii \eta}{\Lambda} \right)^2}\right) \notag \\
  &+ \sum_{av} \frac{W_{pa,v}^{\GW}W_{qa,v}^{\GW}}{\omega - \Omega_{a,v}^{\dRPA} + \ii \eta}\left(1-e^{-2\left( \frac{\omega - \Omega_{a,v}^{\dRPA} + \ii \eta}{\Lambda} \right)^2}\right) \notag
\end{align}
Similarly to the SRG-MP2 case, here we see that for small $\Lambda$ only the divergent denominators are removed.
These denominators are precisely the one responsible of the discontinuities which explain the good results of regularized GW to remove discontinuities seen in Enzo's paper.
Enzo also showed that for larger $\Lambda$ the quasiparticle energies are still smooth yet less accurate.
This is because in addition to the divergent denominators we are removing more and more terms from the self-energy which leads to worse energies.

Finally, we discuss the renormalized correlation self-energy introduced in this work.
\begin{align}
     (\Sigma_c^{\GW}(\omega,\Lambda))_{pq} &= \sum_{iv} \frac{W_{pi,v}^{\GW}W_{qi,v}^{\GW}}{\omega -\Omega_{i,v}^{\dRPA} - \ii \eta}e^{-\left( \frac{\epsilon_p - \Omega_{i,v}^{\dRPA} }{\Lambda} \right)^2}   e^{-\left( \frac{\epsilon_p - \Omega_{i,v}^{\dRPA} }{\Lambda} \right)^2}    \notag \\
     &+ \sum_{av} \frac{W_{pa,v}^{\GW} W_{qa,v}^{\GW}}{\omega - \Omega_{a,v}^{\dRPA} + \ii \eta} e^{-\left( \frac{\epsilon_p - \Omega_{a,v}^{\dRPA}}{\Lambda} \right)^2}e^{-\left( \frac{\epsilon_q - \Omega_{a,v}^{\dRPA} }{\Lambda} \right)^2}  \notag
 \end{align}
 In this case the situation is reversed, \ie the divergent denominators will be the last removed when $\Lambda$ is increased.
 Therefore the renormalized self-energy seems not to be the good strategy to remove discontinuities.
 However, it defines SRG-PT2 approximations to the quasiparticle energies which have the same pros as the SRG-MP2 discussed above.
 Schematically, the determinants included in $\Tilde{F}$ and $\Tilde{\Sigma}$ are like in the following figure where blue means that the determinant is included.

 \begin{figure}
   \centering
   \includegraphics[width=0.9\linewidth]{Figures/renormalizedF.pdf}
   \caption{Determinants included at a finite value of $s$ according to the diagonal denominators of $\Tilde{F}$}
   \label{fig:fig_1}
 \end{figure}


%=================================================================%
\section{Preliminary results}
\label{sec:results}
%=================================================================%

  
%%%%%%%%%%%%%%%%%%%%%% 
\subsection{$\eta = 0$}
%%%%%%%%%%%%%%%%%%%%%% 

\begin{widetext}
  
  \begin{table}
    \caption{Symmetrized static form DIIS = 1}
    \begin{tabular}{rrrrr}
      & He & Ne & H2 & Ar \\
      NbIter & 2.0000 & 13.0000 & 115.0000 & 6.0000 \\
      HOMO & -24.3588 & -21.0664 & -16.1493 & -15.3200 \\
      EHF & -2.8552 & -128.4885 & -1.1280 & -526.8000 \\
      RPA & -0.0460 & -0.2271 & -0.0453 & -0.1775 \\
    \end{tabular}
  \end{table}
  
  \begin{table}
    \caption{SRG static for  DIIS = 1}
    \begin{tabular}{rrrrrrr}
    & He & Ne & H2 & Ar & LiF & HCl \\
      NbIter & 2.0000 & 13.0000 & 51.0000 & 4.0000 & 207.0000 & 163.0000 \\
      HOMO & -24.3516 & -20.9991 & -16.2258 & -15.3139 & -10.9319 & -12.3324 \\
      EHF & -2.8552 & -128.4887 & -1.1286 & -526.8002 & -106.9429 & -460.0891 \\
      RPA & -0.0460 & -0.2271 & -0.0459 & -0.1774 & -0.2481 & -0.1935 \\
    \end{tabular}
  \end{table}

  \begin{table}
    \caption{Symmetrized static form DIIS = 3}
    \begin{tabular}{rrrrrrrrrr}
      & He & Ne & H2 & LiH & HF & Ar & LiF & F2 \\
      NbIter & 2.0000 & 4.0000 & 11.0000 & 58.0000 & 227.0000 & 4.0000 & 234.0000 & 132.0000 \\
      HOMO & -24.3588 & -21.0662 & -16.1493 & -7.9603 & -15.5864 & -15.3198 & -11.0257 & -15.6825 \\
      EHF & -2.8552 & -128.4885 & -1.1280 & -7.9834 & -100.0101 & -526.7999 & -106.9414 & -198.6606 \\
      RPA & -0.0460 & -0.2271 & -0.0453 & -0.0424 & -0.2447 & -0.1775 & -0.2493 & -0.4709 \\
    \end{tabular}
  \end{table}

  
  \begin{table}
    \caption{SRG static form DIIS = 3}
    \begin{tabular}{rrrrrrrrrrr}
      & He & Ne & H2 & LiH & HF & Ar & LiF & HCl & BF & F2 \\
      NbIter & 2.0000 & 4.0000 & 38.0000 & 248.0000 & 91.0000 & 4.0000 & 28.0000 & 141.0000 & 140.0000 & 52.0000 \\
      HOMO & -24.3516 & -20.9990 & -16.2258 & -7.8883 & -15.6498 & -15.3138 & -10.9314 & -12.3324 & -10.8823 & -15.7551 \\
      EHF & -2.8552 & -128.4887 & -1.1286 & -7.9833 & -100.0190 & -526.8000 & -106.9429 & -460.0891 & -124.1050 & -198.6830 \\
      RPA & -0.0460 & -0.2271 & -0.0459 & -0.0425 & -0.2436 & -0.1774 & -0.2484 & -0.1936 & -0.3130 & -0.4684 \\
    \end{tabular}
  \end{table}
\end{widetext}

%%%%%%%%%%%%%%%%%%%%%%
\subsection{$\eta = 0.01$}
%%%%%%%%%%%%%%%%%%%%%%

\begin{widetext}
  
  \begin{table}
    \caption{Symmetrized static form DIIS = 1}
    \begin{tabular}{rrrrrrr}
      & He & Ne & H2 & LiH & Ar & LiF \\
      NbIter & 2.0000 & 13.0000 & 60.0000 & 22.0000 & 6.0000 & 29.0000 \\
      HOMO & -24.3588 & -21.0664 & -16.1499 & -7.9504 & -15.3200 & -11.0513 \\
      EHF & -2.8552 & -128.4885 & -1.1280 & -7.9834 & -526.8000 & -106.9418 \\
      RPA & -0.0460 & -0.2271 & -0.0453 & -0.0424 & -0.1775 & -0.2493 \\
    \end{tabular}
  \end{table}
  
  \begin{table}
    \caption{SRG static for  DIIS = 1}
    \begin{tabular}{rrrrrrrrr}
      & He & Ne & H2 & LiH & Ar & LiF & HCl & CH4 \\
      NbIter & 2.0000 & 13.0000 & 47.0000 & 14.0000 & 4.0000 & 240.0000 & 19.0000 & 172.0000 \\
      HOMO & -24.3516 & -20.9991 & -16.2258 & -7.8880 & -15.3139 & -10.9312 & -12.3324 & -14.3818 \\
      EHF & -2.8552 & -128.4887 & -1.1286 & -7.9833 & -526.8002 & -106.9428 & -460.0891 & -40.1982 \\
      RPA & -0.0460 & -0.2271 & -0.0459 & -0.0423 & -0.1774 & -0.2485 & -0.1935 & -0.2289 \\
    \end{tabular}
  \end{table}

  \begin{table}[h!]
    \caption{Symmetrized static form DIIS = 3}
    \begin{tabular}{rrrrrrrrrrr}
                    & He            & Ne                & H2             & Li2             & LiH          & HF              & Ar                & H2O         & LiF                & HCl  \\ 
      NbIter & 2.0000     & 4.0000        & 9.0000     & 12.0000   & 20.0000 & 18.0000    & 4.0000       & 33.0000  & 12.0000      & 34.0000  \\
      HOMO & -24.3588 & -21.0662    & -16.1499 & -5.3026    & -7.8579  & -15.5895   & -15.3198   & -12.1571 & -11.0512    & -12.3535 \\
      EHF      & -2.8552    & -128.4885 & -1.1280    & -14.8685 & -7.9826  & -100.0089 & -526.7999 & -76.0202 & -106.9418 & -460.0862 \\
      RPA      & -0.0460    & -0.2271      & -0.0453    & -0.0382   & -0.0429  & -0.2454      & -0.1775      & -0.2507   & -0.2493      & -0.1946 \\
      \hline
                    & BeO          & CO               & N2             & & BH3          & NH3        & BF               & BN              & SH2 \\
      NbIter & 27.0000 & 18.0000 & 18.0000 & & 15.0000 & 33.0000 & 17.0000 & 26.0000 & 13.0000 \\
      HOMO & -10.1480 & -14.1411 & -15.4355 & & -13.3479 & -10.5387 & -10.6645 & -11.6532 & -10.0743 \\
      EHF      & -89.3969 & -112.6821 & -108.9450 & & -26.3889 & -56.1865 & -124.0993 & -78.8715 & -398.6901 \\
      RPA      & -0.2867 & -0.3626 & -0.3538 & & -0.1504 & -0.2446 & -0.3107 & -0.3085 & -0.2056 \\
    \end{tabular}
  \end{table}

  
  \begin{table}[h!]
    \caption{SRG static form DIIS = 3}
    \begin{tabular}{rrrrrrrrrrrrrrrrrrrrr}
      Mol        & He             & Ne              & H2            & Li2             & LiH          & HF               & Ar                 & H2O        & LiF                & HCl             \\ 
      NbIter   & 2.0000     & 4.0000       & 8.0000    & 9.0000     & 12.0000 & 23.0000     & 4.0000       & 10.0000  & 7.0000         & 17.0000    \\
      HOMO & -24.3516 & -20.9990   & -16.2258 & -5.2716    & -7.8879  & -15.6498    & -15.3138   & -12.1953 & -10.9314    & -12.3324   \\
      EHF      & -2.8552    & -128.4887 & -1.1286   & -14.8692 & -7.9833   & -100.0190 & -526.8000 & -76.0262 & -106.9429 & -460.0891 \\
      RPA      & -0.0460    & -0.2271     & -0.0459    & -0.0383   & -0.0423   & -0.2436      & -0.1774     & -0.2492    & -0.2485      & -0.1935      \\
      \hline
      Mol       & BeO & CO & N2 & CH4 & BH3 & NH3 & BF & BN & SH2 & F2 \\
      NbIter  & 13.0000 & 12.0000 & 53.0000 & 39.0000 & 113.0000 & 35.0000 & 8.0000 & 80.0000 & 8.0000 & 19.0000 \\
      HOMO & -10.0058 & -14.1362 & -15.4820 & -14.3811 & -13.3194 & -10.5563 & -10.8822 & -11.5763 & -10.0370 & -15.7552 \\
      EHF       & -89.4032 & -112.6875 & -108.9521 & -40.1982 & -26.3904 & -56.1949 & -124.1050 & -78.8792 & -398.6937 & -198.6830 \\
      RPA       & -0.2848 & -0.3614 & -0.3508 & -0.2288 & -0.1503 & -0.2429 & -0.3121 & -0.3125 & -0.2043 & -0.4684 \\      
    \end{tabular}
  \end{table}
\end{widetext}

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

The particle-hole RPA is defined as
\begin{equation}
  \begin{pmatrix}
    \bA^\phRPA & \bB^\phRPA \\
    - \bB^\phRPA & \bA^\phRPA \\
  \end{pmatrix}
  \begin{pmatrix}
    \bX \\
    \bY \\
  \end{pmatrix} = \boldsymbol{\Omega}
  \begin{pmatrix}
    \bX \\
    \bY \\
  \end{pmatrix}
\end{equation}

with

\begin{align}
  A^\phRPA_{ij,ab} &= (\epsilon_i - \epsilon_a) \delta_{ij}\delta_{ab} + \aeri{ib}{aj} \\
  B^\phRPA_{ij,ab} &= \aeri{ij}{ab}
\end{align}

The equation above correspond to the RPAx version while the direct RPA version used in GW has the same element but with non-antisymmetrized Coulomb integrals.

The particle-particle RPA reads as
\begin{equation}
  \begin{pmatrix}
    \bA^\ppRPA & \bB^\ppRPA \\
    - (\bB^\ppRPA)^\mathsf{T} & - \bC{}{\ppRPA} \\
  \end{pmatrix}
  \begin{pmatrix}
    \bX \\
    \bY \\
  \end{pmatrix} = \boldsymbol{\Omega}
  \begin{pmatrix}
    \bX \\
    \bY \\
  \end{pmatrix}
\end{equation}

with

\begin{align}
  A_{ab,cd}^\ppRPA &= \delta_{ab}\delta_{cd}(\epsilon_a + \epsilon_b) + \aeri{ab}{cd} \\
  B_{ab,ij}^\ppRPA &= \aeri{ab}{ij} \\
  C_{ij,kl}^\ppRPA &= -\delta_{ik}\delta_{jl}(\epsilon_i + \epsilon_j) + \aeri{ij}{kl}
\end{align}

%=================================================================%
\section{Downfolding procedure of the linear GW matrix}
\label{sec:downfolding}
%=================================================================%

\begin{equation}
\left\{ \begin{aligned} 
    \bF{}{}\bR^{\hp} + \bV{}{\hhp} \bR^{\hhp} + \bV{}{\pph} \bR^{\pph} &= E \bR^{\hp}	\\
    (\bV{}{\hhp})^{\mathsf{T}}\bR^{\hp} + \bC{}{\hhp} \bR^{\hhp} &= E \bR^{\hhp}	\\
    (\bV{}{\pph})^{\mathsf{T}}\bR^{\hp} + \bC{}{\pph} \bR^{\pph} &= E \bR^{\pph}
\end{aligned} \right.
\end{equation}

\begin{equation}
  \left( \bF{}{} + \bSig(E) \right) \bR^{\hp} = E \bR^{\hp}
\end{equation}

\begin{align}
  \bSig(E) &= \bV{}{\hhp} (E \mathbb{1} - \bC{}{\hhp})^{-1}  (\bV{}{\hhp})^{\mathsf{T}} \\
           &+ \bV{}{\pph} (E \mathbb{1} - \bC{}{\pph})^{-1}  (\bV{}{\pph})^{\mathsf{T}} \notag
\end{align}

We introduce the matrices of eigenvectors $\bU^{\hhp}$ and $\bU^{\pph}$ such that
\begin{align}
  \bC{}{\hhp} &= \bU^{\hhp} \text{diag}(\epsilon_i - \Omega^{(v)}) (\bU^{\hhp})^{-1} \\
  \bC{}{\pph} &= \bU^{\pph} \text{diag}(\epsilon_a + \Omega^{(v)}) (\bU^{\pph})^{-1} 
\end{align}
with
\begin{align}
  U_{i[ja],(k,v)}^{\hhp} &= \bX_{ja}^{(v)}\delta_{pi} & U_{[ia]b,(c,v)}^{\pph} &= \bX_{ia}^{(v)}\delta_{bp}
\end{align}

Therefore we have
\begin{align}
  (E \mathbb{1} - \bC{}{\hhp})^{-1} &= [\bU^{\hhp} (E \mathbb{1} - \text{diag}(\epsilon_i - \Omega^{(v)}))^{-1} (\bU^{\hhp})^{-1}]^{-1} \notag \\
                                    &= \bU^{\hhp} (\text{diag}(E - \epsilon_i + \Omega^{(v)}))^{-1}(\bU^{\hhp})^{-1} \notag \\
                                    &= \bU^{\hhp} \text{diag}(\frac{1}{E - \epsilon_i + \Omega^{(v)}}) (\bU^{\hhp})^{-1}  
\end{align}

\begin{equation}
  \bSig(E)^{\hhp} = \bV{}{\hhp} \bU^{\hhp} \text{diag}(\frac{1}{E - \epsilon_i + \Omega^{(v)}}) (\bU^{\hhp})^{-1} (\bV{}{\hhp})^{\mathsf{T}}
\end{equation}

We have
\begin{align}
  &(\bV{}{\hhp} \bU^{\hhp})_{p,(k,v)} = \sum_{i[ja]} v_{p,i[ja]} U_{i[ja],(k,v)}  \\
  &= \sum_{i[ja]} \eri{pa}{ij} \bX_{ja}^{(v)}\delta_{ik} \notag \\
  &= \sum_{[ja]} \eri{pa}{kj} \bX_{ja}^{(v)} \notag \\
  &(\bV{}{\hhp} \bU^{\hhp} \text{diag}(\frac{1}{E - \epsilon_i + \Omega^{(v)}}) )_{p,(k,v)} \\
  &=\sum_{l,(w)} \left(\sum_{[ja]} \eri{pa}{lj} \bX_{ja}^{(w)}\right) \text{diag}(\frac{1}{E - \epsilon_i + \Omega^{(v)}})_{(l,w),(k,v)} \delta_{lk}\delta_{(w),(v)} \notag \\
  &= \left(\sum_{[ja]} \eri{pa}{kj} \bX_{ja}^{(v)}\right) \frac{1}{E - \epsilon_k + \Omega^{(v)}}  \notag \\
  &(\bSig(E)^{\hhp})_{p,q} \\
  &= \sum_{(k,v)} \left(\sum_{[ja]} \eri{pa}{kj} \bX_{ja}^{(v)}\right) \frac{1}{E - \epsilon_k + \Omega^{(v)}} \left(\sum_{[ja]} \eri{qa}{kj} \bX_{ja}^{(v)}\right) \notag \\
  &= \sum_{(k,v)} \frac{M_{kp}^{(v)}M_{kq}^{(v)}}{E - \epsilon_k + \Omega^{(v)}} \notag
\end{align}

%=================================================================%
\section{Alternative partitioning}
\label{sec:partitioning}
%=================================================================%

In this section we discuss the alternative partitioning of the Hamiltonian mentioned previously.

%///////////////////////////%
\subsubsection{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}

%///////////////////////////%
\subsubsection{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{widetext}
  
  \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}

% \section{Old stuff}

% 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_r$.
% \begin{align}
%   \label{eq:matrixdiffeq}
%   \dv{\bV{r}{(1)}}{s} &= \bV{}{(1)} (2 \epsilon_r \bC{}{(0)} - \epsilon_r^2\mathbb{1} - (\bC{}{(0)})^2) \\
%   \dv{\bV{r}{(1),\dagger}}{s} &=  (2 \epsilon_r \bC{}{(0)} - \epsilon_r^2\mathbb{1} - (\bC{}{(0)})^2) \bV{}{(1),\dagger} 
% \end{align}
% These equations can be solved if 

% \textbf{Matrix differential equation}
%   A general differential matrix equation of the form
%   \begin{equation}
%     \dv{\bX}{s} = \bA \bX
%   \end{equation}
%   admits a solution
%   \begin{align}
%     \bX(s) &= c_1e^{\lambda_1 s}\bU_1 + c_2e^{\lambda_2 s}\bU_2 + \dots + c_ne^{\lambda_n s}\bU_n  \\
%            &= \bU \text{diag}(e^{\lambda_i s}) \bc \notag
%   \end{align}
%   where $c_i$ are coefficients to determine, $\lambda_i$ are the eigenvalues of $\bA$ and $\bU_i$ the corresponding eigenvectors.
%   The $c_i$ are determined by the initial condition which gives
%   \begin{equation}
%     \label{eq:solution_eqdiff}
%     \bX(s) = \bU \text{diag}(e^{\lambda_i s}) \bU^{-1}\bX(0)
%   \end{equation}

%    In order to solve the matrix differential equation Eq.~(\ref{eq:matrixdiffeq}), we need to diagonalize $2 \epsilon_r \bC{}{(0)} - \epsilon_r^2\mathbb{1} - (\bC{}{(0)})^2$ which is a polynomial in $\bC{}{(0)}$ so the polynomial admits the same eigenvectors $\bC{}{(0)}$.
%   The elements of the matrix $\bC{}{(0)}$ in the GW case are given in Eq.~(\ref{eq:GW_unfolded}) and can be written equivalently as
%    \begin{align}
%     \label{eq:GW_unfolded}
%     C^\text{2h1p}_{i[ja],k[lc]} &= \epsilon_i  \delta_{jl} \delta_{ac} \delta_{ik} - A_{ja,lc}^{\phRPA}\delta_{ik} \notag \\
%     C^\text{2p1h}_{[ia]b,[kc]d} &= \epsilon_i  \delta_{ik} \delta_{ac} \delta_{bd} + A_{ia,kc}^{\phRPA}\delta_{bd} \notag
%    \end{align}
%    So the matrix $\bC{}{(0)}$ is a diagonal block matrix with each block corresponding to a shifted RPA problem.
%    Because we know the eigenvectors of the RPA problem we can buil the eigenvectors of $\bC{}{(0)}$ as
%    \begin{align}
%      U_{i[ja],(p,v)} &= \bX_{ja}^{(v)}\delta_{pi} & U_{[ia]b,(p,v)} &= \bX_{ia}^{(v)}\delta_{bp}
%    \end{align}
%    where the eigenvectors are indexed by a collective index $(n,v)$ where $n$ refers to the block number and $v$ refers to the eigenvector inside the block.
%    The corresponding eigenvalues are
%    \begin{align}
%      \Omega_{(i,v)} &= \epsilon_i - \Omega_v & \Omega_{(a,v)} &= \epsilon_a + \Omega_v 
%    \end{align}


\end{document}