\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{\alert}[1]{\textcolor{red}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#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{\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{\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}{\boldsymbol{F}}
\newcommand{\bU}{\boldsymbol{U}}
\newcommand{\bV}[2]{\boldsymbol{V}_{#1}^{#2}}
\newcommand{\bW}{\boldsymbol{W}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\bY}{\boldsymbol{Y}}
\newcommand{\bZ}{\boldsymbol{Z}}
\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}
%=================================================================%

The aim of this document is two-fold.
First, we want to re-derive (in details) the perturbative analysis of the similarity renormalisation group (SRG) formalism applied to the non-relativistic electronic Hamiltonian.
In a second time, we want to apply the same formalism to the unfolded GW Hamiltonian.
To do so, we first need to find a second quantization effective Hamiltonian for Green function theory.
Before jumping into these analysis, we do a brief presentation of the SRG formalism.

%=================================================================%
\section{The similarity renormalisation group}
%=================================================================%

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 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 electronic Hamiltonian}
%=================================================================%

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(s)_{ij}^{ab}\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) + F(0) + \la V(0)
\end{equation}
In addition, we know the following initial conditions.
We use the HF basis set of the reference such that $F^{\text{od}}(0) = 0$ and $F^{\mathrm{d}}(0)=\delta_{pq}\epsilon_p$
Therefore, we have
\begin{align}
  \bH^\text{d}(0)&=E_0(0) + F^{\mathrm{d}}(0) + \la V^{\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}
%%%%%%%%%%%%%%%%%%%%%%

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.


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{\Delta_{ab}^{ij}}{\aeri{ij}{ab}}\left(1-e^{-2s (\Delta_{ab}^{ij})^2}\right)
\end{equation}

%=================================================================%
\section{The unfolded Green's function}
% =================================================================%

%%%%%%%%%%%%%%%%%%%%%%
\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{}{} 
  \end{pmatrix}
  + \lambda
  \begin{pmatrix}
    \bO & \bV{}{} \\
    \bV{}{\dagger} & \bO 
  \end{pmatrix}
\end{equation}
which gives the following conditions
\begin{align}
  \bHd{0}(0) &=   \begin{pmatrix}
    \bF & \bO \\
    \bO & \bC{}{} 
  \end{pmatrix}  & \bHod{0}(0) &= \bO \\
  \bHd{1}(0) &= \bO & \bHod{1}(0) &=   \begin{pmatrix}
    \bO & \bV{}{} \\
    \bV{}{\dagger} & \bO 
  \end{pmatrix}  
\end{align}

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

The zero-th 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)}\bF^{(0)}\\
    \bC{}{(0)}\bV{}{(1),\dagger} - \bV{}{(1),\dagger}\bC{}{(0)} & \bO 
  \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  \Longleftrightarrow \color{red}{\boxed{\color{black}{\bF^{(1)}= \bO}}} \\
  \dv{\bC{}{(1)}}{s}  &= \bO \Longleftrightarrow \color{red}{\boxed{\color{black}{\bC{}{(1)}= \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}
\end{align}
The two last equations can be solved differently depending on the form of $\bF$ and $\bC{}{}$.
\subsubsection*{Diagonal $\bC{}{(0)}$}
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.

\begin{align}
  (\dv{\bV{}{(1)}}{s})_{pQ} &= (2 \bF^{(0)}\bV{}{(1)}\bC{}{(0)} - (\bF^{(0)})^2\bV{}{(1)} - \bV{}{(1)}(\bC{}{(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!

\subsubsection*{Non-diagonal $\bC{}{(0)}$}
We follow the same development as before
\begin{align}
  (\dv{\bV{}{(1)}}{s})_{pQ} &= (2 \bF^{(0)}\bV{}{(1)}\bC{}{(0)} - (\bF^{(0)})^2\bV{}{(1)} - \bV{}{(1)}(\bC{}{(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} c^{(0)}_{SQ}  \\
                            &- \sum_{rs} \epsilon^{(0)}_p\delta_{pr} \epsilon^{(0)}_r\delta_{rs} v^{(1)}_{sQ} \\
                            &- \sum_{RS} v^{(1)}_{pR} c^{(0)}_{RS}   c^{(0)}_{SQ} \\
                            &= - (\epsilon^{(0)}_p)^2v^{(1)}_{pQ}+  \sum_{S} 2 \epsilon^{(0)}_p v^{(1)}_{pS} c^{(0)}_{SQ}  - \sum_{RS} v^{(1)}_{pR} c^{(0)}_{RS}   c^{(0)}_{SQ} 
\end{align}
We obtain a set of coupled differential equations which seems far from being trivial to solve.
In order to simplify the problem we consider the case when $\bF = \eps_p$.
\begin{align}
  \dv{\bV{}{(1)}}{s} &= 2 \bF^{(0)}\bV{}{(1)}\bC{}{(0)} - (\bF^{(0)})^2\bV{}{(1)} - \bV{}{(1)}(\bC{}{(0)})^2 \\
                     &= 2 \eps_p\bV{}{(1)}\bC{}{(0)} - (\eps_p)^2\bV{}{(1)} - \bV{}{(1)}(\bC{}{(0)})^2  \\
                     &= \bV{}{(1)} (\eps_p\mathbb{1} - \bC{}{(0)})^2
\end{align}
Now to solve this matrix differential equation, we just need to diagonalize $(\eps_p \mathbb{1} - \bC{}{(0)})^2$.
Fortunately, this can be easily done because the eigenvalues of $\bC{}{(0)}$ are known to be the shifted RPA eigenvalues and the eigenvectors are given in Bintrim 2021.

\textbf{\color{red}{IDEA: Can we put the non-diagonal part of C in the off-diag H?}}

%%%%%%%%%%%%%%%%%%%%%%
\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}} \\
  &= \begin{pmatrix}
    \bO & \bF^{(0)}\bV{}{(2)} - \bV{}{(2)}\bF^{(0)}\\
    \bC{}{(0)}\bV{}{(2),\dagger} - \bV{}{(2),\dagger}\bC{}{(0)} & \bO 
  \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} \\
  \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}

Once again the integration of these equations is much simpler if $\bC{}{(0)}$ is diagonal.

\subsubsection*{Diagonal $\bC{}{(0)}$}

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

A similar derivation should give (\textbf{\textcolor{red}{TO CHECK}})

\begin{align}
 &c^{(2)}_{PQ}(s) = \notag \\
 &\color{red}{\boxed{\color{black}{- \sum_r \frac{\Delta\eps^{(0)}_{P}  + \Delta\eps^{(0)}_{Q} - 2 \eps^{(0)}_r}{(\Delta\eps^{(0)}_P - \eps^{(0)}_r)^2+ (\Delta\eps^{(0)}_Q - \eps^{(0)}_r)^2}(1 - e^{-s [ (\Delta\eps^{(0)}_P - \eps^{(0)}_r)^2+ (\Delta\eps^{(0)}_P - \eps^{(0)}_r)^2]})}}} \notag
\end{align}

\begin{align}
  &\dv{v^{(2)}_{pQ}}{s}  = - (\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_Q )^2 v^{(2)}_{pQ} \\
  &\color{red}{\boxed{\color{black}{v^{(2)}_{pQ}(s) = v^{(2)}_{pQ}(0) e^{-s(\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_Q )^2} = 0 }}} 
\end{align}


\subsubsection*{Non-diagonal $\bC{}{(0)}$}


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

\end{document}