saving work, GW and GF(2) up to first order, working on second order

This commit is contained in:
Antoine Marie 2022-11-08 14:18:17 +01:00
parent 4fd96bc8d4
commit 9676c5ac60

View File

@ -13,126 +13,15 @@
]{hyperref} ]{hyperref}
\urlstyle{same} \urlstyle{same}
\newcommand{\ie}{\textit{i.e.}} %============================================================%
\newcommand{\eg}{\textit{e.g.}} %%% NEWCOMMANDS %%%
\newcommand{\etal}{\textit{et al.}} % ============================================================%
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\input{Commands}
\usepackage[normalem]{ulem} \usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}} \newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\ant}[1]{\textcolor{green}{#1}} \newcommand{\ant}[1]{\textcolor{green}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{r\ed}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\QP}{\textsc{quantum package}}
\newcommand{\T}[1]{#1^{\intercal}}
\newcommand{\Sig}[2]{\Sigma_{#1}^{#2}}
\newcommand{\dRPA}{\text{dRPA}}
% coordinates
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bx}{\boldsymbol{x}}
\newcommand{\dbr}{d\br}
\newcommand{\dbx}{d\bx}
% methods
\newcommand{\GW}{\text{$GW$}}
\newcommand{\GT}{\text{$GT$}}
\newcommand{\evGW}{ev$GW$}
\newcommand{\qsGW}{qs$GW$}
\newcommand{\GOWO}{$G_0W_0$}
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\xc}{\text{xc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\co}{\text{c}}
\newcommand{\x}{\text{x}}
\newcommand{\KS}{\text{KS}}
\newcommand{\HF}{\text{HF}}
\newcommand{\RPA}{\text{RPA}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
\newcommand{\sERI}[2]{(#1|#2)}
\newcommand{\e}[2]{\epsilon_{#1}^{#2}}
%
\newcommand{\Ne}{N}
\newcommand{\Norb}{K}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}
% operators
\newcommand{\hH}{\Hat{H}}
\newcommand{\hS}{\Hat{S}}
\newcommand{\ani}[1]{\hat{a}_{#1}}
\newcommand{\cre}[1]{\hat{a}_{#1}^\dagger}
\newcommand{\no}[2]{\mleft\{ \hat{a}_{#1}^{#2}\mright\} }
% energies
\newcommand{\Enuc}{E^\text{nuc}}
\newcommand{\Ec}[1]{E_\text{c}^{#1}}
\newcommand{\EHF}{E^\text{HF}}
% orbital energies
\newcommand{\eps}{\epsilon}
\newcommand{\reps}{\Tilde{\epsilon}}
% Matrix elements
\newcommand{\SigC}{\Sigma^\text{c}}
\newcommand{\rSigC}{\Tilde{\Sigma}^\text{c}}
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\SO}[1]{\psi_{#1}}
\newcommand{\eri}[2]{\braket{#1}{#2}}
\newcommand{\aeri}[2]{\mel{#1}{}{#2}}
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\rbra}[1]{(#1|}
\newcommand{\rket}[1]{|#1)}
% Matrices
\newcommand{\bO}{\boldsymbol{0}}
\newcommand{\bI}{\boldsymbol{1}}
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\bSig}{\boldsymbol{\Sigma}}
\newcommand{\bSigC}{\boldsymbol{\Sigma}^{\text{c}}}
\newcommand{\be}{\boldsymbol{\epsilon}}
\newcommand{\bOm}{\boldsymbol{\Omega}}
\newcommand{\bA}{\boldsymbol{A}}
\newcommand{\bB}{\boldsymbol{B}}
\newcommand{\bC}[2]{\boldsymbol{C}_{#1}^{#2}}
\newcommand{\bD}{\boldsymbol{D}}
\newcommand{\bF}[2]{\boldsymbol{F}_{#1}^{#2}}
\newcommand{\bR}{\boldsymbol{R}}
\newcommand{\bU}{\boldsymbol{U}}
\newcommand{\bV}[2]{\boldsymbol{V}_{#1}^{#2}}
\newcommand{\bW}{\boldsymbol{W}}
\newcommand{\bX}[2]{\boldsymbol{X}_{#1}^{#2}}
\newcommand{\bY}[2]{\boldsymbol{Y}_{#1}^{#2}}
\newcommand{\bZ}[2]{\boldsymbol{Z}_{#1}^{#2}}
\newcommand{\bc}{\boldsymbol{c}}
% orbitals, gaps, etc
\newcommand{\IP}{I}
\newcommand{\EA}{A}
\newcommand{\HOMO}{\text{HOMO}}
\newcommand{\LUMO}{\text{LUMO}}
\newcommand{\Eg}{E_\text{g}}
\newcommand{\EgFun}{\Eg^\text{fund}}
\newcommand{\EgOpt}{\Eg^\text{opt}}
\newcommand{\EB}{E_B}
% shortcuts for greek letters
\newcommand{\si}{\sigma}
\newcommand{\la}{\lambda}
\newcommand{\RHH}{R_{\ce{H-H}}}
\newcommand{\ii}{\mathrm{i}}
\newcommand{\bEta}[1]{\boldsymbol{\eta}^{(#1)}(s)}
\newcommand{\bHd}[1]{\bH_\text{d}^{(#1)}}
\newcommand{\bHod}[1]{\bH_\text{od}^{(#1)}}
% addresses % addresses
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France} \newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
@ -185,7 +74,7 @@ Therefore, the transformed Hamiltonian
depends on a flow parameter $s$. depends on a flow parameter $s$.
The resulting Hamiltonian possess up to $N$-body operators with $N$ the number of particle. The resulting Hamiltonian possess up to $N$-body operators with $N$ the number of particle.
\begin{equation} \begin{equation}
\bH(s) = E_0(s) + \bF(s) + \bV{}{}(s) + \bW(s) + \dots \bH(s) = E_0(s) + \bF{}{}(s) + \bV{}{}(s) + \bW(s) + \dots
\end{equation} \end{equation}
In the following, we will truncate every contribution superior to two-body operators. 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 We can easily derive an evolution equation for this Hamiltonian by taking the derivative of $\bH(s)$. This gives
@ -240,15 +129,15 @@ Hence, we define the off-diagonal Hamiltonian as
\end{equation} \end{equation}
Note that each coefficients depend on $s$. Note that each coefficients depend on $s$.
The perturbative parameter $\la$ is such that The perturbative parameter $\lambda$ is such that
\begin{equation} \begin{equation}
\bH(0) = E_0(0) + \bF{}{}(0) + \la \bV{}{}(0) \bH(0) = E_0(0) + \bF{}{}(0) + \lambda \bV{}{}(0)
\end{equation} \end{equation}
In addition, we know the following initial conditions. 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$. 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 Therefore, we have
\begin{align} \begin{align}
\bH^\text{d}(0)&=E_0(0) + \bF{}{\mathrm{d}}(0) + \la \bV{}{\mathrm{d}}(0) & \bH^\text{od}(0)&= \la V^{\mathrm{od}}(0) \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} \end{align}
Now, we want to compute the terms at each order of the following development Now, we want to compute the terms at each order of the following development
\begin{equation} \begin{equation}
@ -398,30 +287,30 @@ 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. The GF($n$) formalism is defined such that the self-energy includes every diagram up to $n$-th order of M\"oller-Plesset perturbation theory.
\begin{align} \begin{align}
\label{eq:GF2_selfenergy} \label{eq:GF2_selfenergy}
\Sig{pq}{GF(2)}(\omega) &= \frac{1}{2} \sum_{ija} \frac{\aeri{pa}{ij}\aeri{qa}{ij}}{\omega + \eps _a -\eps_i -\eps_j - \ii \eta} \notag \\ \Sigma_{pq}^{GF(2)}(\omega) &= \frac{1}{2} \sum_{ija} \frac{\aeri{pa}{ij}\aeri{qa}{ij}}{\omega + \epsilon _a -\epsilon_i -\epsilon_j - \ii \eta} \notag \\
&+ \frac{1}{2} \sum_{iab} \frac{\aeri{pi}{ab}\aeri{qi}{ab}}{\omega + \eps _i -\eps_a -\eps_b + \ii \eta} \notag &+ \frac{1}{2} \sum_{iab} \frac{\aeri{pi}{ab}\aeri{qi}{ab}}{\omega + \epsilon _i -\epsilon_a -\epsilon_b + \ii \eta} \notag
\end{align} \end{align}
On the other hand, the GW self-energy is obtained by taking the RPA polarizability and removing the vertex correction in the exact definition of the self-energy. 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} \begin{equation}
\label{eq:GW_selfenergy} \label{eq:GW_selfenergy}
\Sig{pq}{\GW}(\omega) = \sum_{im} \frac{\sERI{pi}{m} \sERI{qi}{m}}{\omega - \e{i}{} + \Om{m}{\dRPA} - \ii \eta} + \sum_{am} \frac{\sERI{pa}{m} \sERI{qa}{m}}{\omega - \e{a}{} - \Om{m}{\dRPA} + \ii \eta} \notag \Sigma_{pq}^{\GW}(\omega) = \sum_{iv} \frac{\ceri{pi}{v} \ceri{qi}{v}}{\omega - \epsilon_i + \Omega_{v}^{\dRPA} - \ii \eta} + \sum_{am} \frac{\ceri{pa}{v} \ceri{qa}{v}}{\omega - \epsilon_a - \Omega_{v}^{\dRPA} + \ii \eta} \notag
\end{equation} \end{equation}
with with
\begin{equation} \begin{equation}
\label{eq:GW_sERI} \label{eq:GW_sERI}
\sERI{pq}{m} = \sum_{ia} \eri{pq}{ia} \qty( \bX{m}{\dRPA} + \bY{m}{\dRPA} )_{ia} \notag \ceri{pq}{m} = \sum_{ia} \eri{pi}{qa} \qty( \bX_{m}^{\dRPA} + \bY_{m}^{\dRPA} )_{ia} \notag
\end{equation} \end{equation}
Finally, the GT approximation corresponds to another approximation to the polarizability than in GW, namely the one coming from pp-hh-RPA 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 The corresponding self-energies read as
\begin{equation} \begin{equation}
\label{eq:GT_selfenergy} \label{eq:GT_selfenergy}
\Sig{pq}{\GT}(\omega) = \sum_{im} \frac{\eri{pi}{\chi^{N+2}_m}\eri{qi}{\chi^{N+2}_m}}{\omega + \e{i}{} - \Om{m}{N+2} + \ii \eta} + \sum_{am} \frac{\eri{pa}{\chi^{N-2}_m}\eri{qa}{\chi^{N-2}_m}}{\omega + \e{a}{} - \Om{m}{N-2} - \ii \eta} \notag \Sigma_{pq}^{\GT}(\omega) = \sum_{im} \frac{\eri{pi}{\chi^{N+2}_m}\eri{qi}{\chi^{N+2}_m}}{\omega + \epsilon_i - \Omega_{m}^{N+2} + \ii \eta} + \sum_{am} \frac{\eri{pa}{\chi^{N-2}_m}\eri{qa}{\chi^{N-2}_m}}{\omega + \epsilon_a - \Omega_{m}^{N-2} - \ii \eta} \notag
\end{equation} \end{equation}
with with
\begin{align} \begin{align}
\label{eq:GT_sERI} \label{eq:GT_sERI}
\eri{pq}{\chi^{N+2}_m} &= \sum_{c<d} \aeri{pq}{cd} \bX{cd,m}{N+2} + \sum_{k<l} \aeri{pq}{kl} \bY{kl,m}{N+2} \notag \\ \eri{pq}{\chi^{N+2}_m} &= \sum_{c<d} \aeri{pq}{cd} \bX_{cd,m}^{N+2} + \sum_{k<l} \aeri{pq}{kl} \bY_{kl,m}^{N+2} \notag \\
\eri{pq}{\chi^{N-2}_m} &= \sum_{c<d} \aeri{pq}{cd} \bX{cd,m}{N-2} + \sum_{k<l} \aeri{pq}{kl} \bY{kl,m}{N-2} \notag \eri{pq}{\chi^{N-2}_m} &= \sum_{c<d} \aeri{pq}{cd} \bX_{cd,m}^{N-2} + \sum_{k<l} \aeri{pq}{kl} \bY_{kl,m}^{N-2} \notag
\end{align} \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}. The two RPA problems giving the eigenvectors needed to build the GW and GT self-energies are given in Appendix~\ref{sec:rpa}.
@ -442,39 +331,43 @@ For the three approximations considered here, the three matrices $\bH$ share the
\bH = \bH =
\begin{pmatrix} \begin{pmatrix}
\bF{}{} & \bV{}{\text{2h1p}} & \bV{}{\text{2p1h}} \\ \bF{}{} & \bV{}{\text{2h1p}} & \bV{}{\text{2p1h}} \\
\T{(\bV{}{\text{2h1p}})} & \bC{}{\text{2h1p}} & \bO \\ (\bV{}{\text{2h1p}})^{\mathsf{T}} & \bC{}{\text{2h1p}} & \bO \\
\T{(\bV{}{\text{2p1h}})} & \bO & \bC{}{\text{2p1h}} \\ (\bV{}{\text{2p1h}})^{\mathsf{T}} & \bO & \bC{}{\text{2p1h}} \\
\end{pmatrix} \end{pmatrix}
\end{equation} \end{equation}
The expression of the coupling blocks $\bV{}{}$ and the diagonal blocks $\bC{}{}$ in the different cases is given below. The expression of the coupling blocks $\bV{}{}$ and the diagonal blocks $\bC{}{}$ in the different cases is given below.
\begin{itemize} \begin{itemize}
\item \textbf{GF(2)} \item \textbf{GF(2)}
\begin{align} \begin{align}
\label{eq:GF2_unfolded}
V^\text{2h1p}_{p,klc} & = \frac{1}{\sqrt{2}}\aeri{pc}{kl} V^\text{2h1p}_{p,klc} & = \frac{1}{\sqrt{2}}\aeri{pc}{kl}
& &
V^\text{2p1h}_{p,kcd} & = \frac{1}{\sqrt{2}}\aeri{pk}{dc} \\ V^\text{2p1h}_{p,kcd} & = \frac{1}{\sqrt{2}}\aeri{pk}{dc} \\
C^\text{2h1p}_{ija,klc} & = \qty( \e{i}{} + \e{j}{} - \e{a}{}) \delta_{jl} \delta_{ac} \delta_{ik} C^\text{2h1p}_{ija,klc} & = \qty( \epsilon_i + \epsilon_j - \epsilon_a) \delta_{jl} \delta_{ac} \delta_{ik}
& &
C^\text{2p1h}_{iab,kcd} & = \qty( \e{a}{} + \e{b}{} - \e{i}{}) \delta_{ik} \delta_{ac} \delta_{bd} \notag C^\text{2p1h}_{iab,kcd} & = \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} \delta_{bd} \notag
\end{align} \end{align}
\item \textbf{GW} \item \textbf{GW}
\begin{align} \begin{align}
\label{eq:GW_unfolded}
V^\text{2h1p}_{p,klc} & = \eri{pc}{kl} V^\text{2h1p}_{p,klc} & = \eri{pc}{kl}
& &
V^\text{2p1h}_{p,kcd} & = \eri{pk}{dc} \notag \\ V^\text{2p1h}_{p,kcd} & = \eri{pk}{dc} \notag \\
C^\text{2h1p}_{ija,klc} &= \qty[ \qty( \e{i}{} + \e{j}{} - \e{a}{}) \delta_{jl} \delta_{ac} - \eri{jc}{al} ] \delta_{ik} & & C^\text{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( \e{a}{} + \e{b}{} - \e{i}{}) \delta_{ik} \delta_{ac} + \eri{ak}{ic} ] \delta_{bd} \notag & & C^\text{2p1h}_{iab,kcd} &= \qty[ \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} + \eri{ak}{ic} ] \delta_{bd} \notag & &
\end{align} \end{align}
\item \textbf{GT} \item \textbf{GT}
\begin{align} \begin{align}
\label{eq:GT_unfolded}
V^\text{2h1p}_{p,klc} & = \aeri{pc}{kl} V^\text{2h1p}_{p,klc} & = \aeri{pc}{kl}
& &
V^\text{2p1h}_{p,kcd}&= \aeri{pk}{cd} \notag \\ V^\text{2p1h}_{p,kcd}&= \aeri{pk}{cd} \notag \\
C^\text{2h1p}_{ija,klc} &= \qty[ \qty( \e{i}{} + \e{j}{} - \e{a}{}) \delta_{jl} \delta_{ac} - \aeri{ij}{kl} ] \delta_{ac} & & \\ C^\text{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( \e{a}{} + \e{b}{} - \e{i}{}) \delta_{ik} \delta_{ac} + \aeri{ab}{cd} ] \delta_{ik} & & \notag C^\text{2p1h}_{iab,kcd} &= \qty[ \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} + \aeri{ab}{cd} ] \delta_{ik} & & \notag
\end{align} \end{align}
\end{itemize} \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 electron-hole T matrix to see if it also correspond to one term that can be found in the CI below.}}
@ -646,7 +539,7 @@ The expression in the other case are given in Appendix~\ref{sec:diagC}.
\subsection{Integrating order by order} \subsection{Integrating order by order}
%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%
In the following, upper case indices correspond to the 2h1p and 2p1h sectors while lower case indices correspond to the 1h and 1p sectors. Also the $\Delta\eps_R$ corresponds to the diagonal elements of the 2h1p and 2p1h sectors. 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\epsilon_R$ corresponds to the diagonal elements of the 2h1p and 2p1h sectors.
\subsubsection{First order} \subsubsection{First order}
@ -660,13 +553,83 @@ The differential equation for the coupling blocks can be solved in the GF(2) cas
However, in the general case this matrix differential equation is not trivial to solve. However, in the general case this matrix differential equation is not trivial to solve.
In practice, we often resort to the \GOWO~or \evGW~schemes which implies that we only need to solve the above equations for $\bF{}{} = \epsilon_p$. In practice, we often resort to the \GOWO~or \evGW~schemes which implies that we only need to solve the above equations for $\bF{}{} = \epsilon_p$.
\begin{align} \begin{align}
\label{eq:matrixdiffeq}
\dv{\bV{}{(1)}}{s} &= \bV{}{(1)} (2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2) \\ \dv{\bV{}{(1)}}{s} &= \bV{}{(1)} (2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2) \\
\dv{\bV{}{(1),\dagger}}{s} &= (2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2) \bV{}{(1),\dagger} \\ \dv{\bV{}{(1),\dagger}}{s} &= (2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2) \bV{}{(1),\dagger}
\end{align} \end{align}
These matrix differential equations can be solved if we know how to diagonalize $2 \epsilon_p \bC{}{(0)} - \epsilon_p^2- (\bC{}{(0)})^2$.
We know how to diagonalize $\bC{}{(0)}$ so we know how to diagonalize polynomial of $\bC{}{(0)}$.
\textcolor{red}{\textbf{TODO Give analytical expression for the different cases.}} \begin{itemize}
\item \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}
\item \textbf{GW}
In order to solve the matrix differential equation Eq.~(\ref{eq:matrixdiffeq}), we need to diagonalize $2 \epsilon_p \bC{}{(0)} - \epsilon_p^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}
Therefore the eigenvalues of $2 \epsilon_p \bC{}{(0)} - \epsilon_p^2\mathbb{1} - (\bC{}{(0)})^2$ are $2 \epsilon_p \Omega_{(q,v)} - \epsilon_p^2 - \Omega_{q,v}^2 = -(\epsilon_p - \Omega_{(q,v)})^2$.
And finally the analytical expressions for the GW coupling blocks at first order are
\begin{align}
\bV{}{\hhp,(1)}(s) &= \bU \text{diag}(e^{-(\epsilon_p - \Omega_{(i,v)})^2s})\bU^{-1}\bV{}{\hhp,(1)}(0) \\
\bV{}{\pph,(1)}(s) &= \bU \text{diag}(e^{-(\epsilon_p - \Omega_{(a,v)})^2s})\bU^{-1}\bV{}{\pph,(1)}(0)
\end{align}
Therefore the downfolded SRG self-energy is
\begin{align}
&\bSig(\omega)^{\hhp} = \bV{}{\hhp,(1)}(0) \bU^{\hhp} \text{diag}(e^{-(\epsilon_p - \Omega_{(i,v)})^2s}) \\
&\text{diag}(\frac{1}{\omega - \epsilon_i + \Omega^{(v)}}) \text{diag}(e^{-(\epsilon_p - \Omega_{(i,v)})^2s}) (\bU^{\hhp})^{-1} (\bV{}{\hhp,(1)}(0))^{\mathsf{T}} \notag \\
&= \bV{}{\hhp,(1)}(0) \bU^{\hhp} \text{diag}(\frac{e^{-2(\epsilon_p - \Omega_{(i,v)})^2s}}{\omega - \epsilon_i + \Omega^{(v)}}) (\bU^{\hhp})^{-1} (\bV{}{\hhp,(1)}(0))^{\mathsf{T}} \notag
\end{align}
Renaming the index $p$ as $r$ and using a similar matrix product as in Appendix~\ref{sec:downfolding}, we finally obtain
\begin{align}
\label{eq:SRGGW_selfenergy}
\Sigma(\omega)_{pq} &= \sum_{(i,v)} \frac{M_{ip}^{(v)}M_{iq}^{(v)}}{\omega - \epsilon_i + \Omega^{(v)}}e^{-2(\epsilon_r - \Omega_{(i,v)})^2s} \notag \\
&+ \sum_{(a,v)} \frac{M_{ap}^{(v)}M_{aq}^{(v)}}{\omega - \epsilon_a - \Omega^{(v)}}e^{-2(\epsilon_r - \Omega_{(a,v)})^2s}
\end{align}
where the $\pph$ part has been obtained by an analog derivation.
\item \textbf{GF(2)}
The GF(2) case is much easier because $\bC{}{(0)}$ is already diagonal so the matrices $\bU$ and $\bU^{-1}$ are equal to the identity.
A compeltely analog derivation gives
\begin{align}
\label{eq:SRGGF2_selfenergy}
\Sigma_{pq}^{GF(2)}(\omega) &= \frac{1}{2} \sum_{ija} \frac{\aeri{pa}{ij}\aeri{qa}{ij}}{\omega + \epsilon _a -\epsilon_i -\epsilon_j - \ii \eta} e^{-2(\epsilon_r - \Delta_{ij}^a)^2s}\notag \\
&+ \frac{1}{2} \sum_{iab} \frac{\aeri{pi}{ab}\aeri{qi}{ab}}{\omega + \epsilon _i -\epsilon_a -\epsilon_b + \ii \eta} e^{-2(\epsilon_r - \Delta_{i}^{ab})^2s} \notag
\end{align}
\item \textbf{GT}
\textcolor{red}{\textbf{TODO Give analytical expression for the GT case.}}
\end{itemize}
\subsubsection{Second order} \subsubsection{Second order}
@ -753,6 +716,111 @@ In these equations $P(r s)$ is the antisymmetric permutation operator.
\label{sec:rpa} \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{Perturbative matrix coefficients for $C^{(0)}$ diagonal} \section{Perturbative matrix coefficients for $C^{(0)}$ diagonal}
\label{sec:diagC} \label{sec:diagC}
@ -780,11 +848,11 @@ Note the close similarity with Evangelista's expressions for the off-diagonal pa
\begin{align} \begin{align}
&(\dv{\bF{}{(2)}}{s})_{pq} = (\bF{}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bF{}{(0)} - 2 \bV{}{(1)}\bC{\text{d}}{(0)}\bV{}{(1),\dagger})_{pq} \notag \\ &(\dv{\bF{}{(2)}}{s})_{pq} = (\bF{}{(0)}\bV{}{(1)}\bV{}{(1),\dagger} + \bV{}{(1)}\bV{}{(1),\dagger}\bF{}{(0)} - 2 \bV{}{(1)}\bC{\text{d}}{(0)}\bV{}{(1),\dagger})_{pq} \notag \\
&= \sum_{rS} f^{(0)}_{pr} v^{(1)}_{rS} v^{(1),\dagger}_{Sq} + \sum_{Rs} v^{(1)}_{pR} v^{(1),\dagger}_{Rs} f^{(0)}_{sq} - 2\sum_{RS} v^{(1)}_{pR} c^{(0)}_{RS} v^{(1),\dagger}_{Sq} \notag \\ &= \sum_{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_{S} \epsilon^{(0)}_{p} v^{(1)}_{pS} v^{(1)}_{qS} + \sum_{R} \epsilon^{(0)}_{q} v^{(1)}_{pR} v^{(1)}_{qR} - 2\sum_{R} \Delta\epsilon^{(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 (\epsilon^{(0)}_{p} + \epsilon^{(0)}_{q} - 2 \Delta\epsilon^{(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 \\ &= \sum_R (\epsilon^{(0)}_{p} + \epsilon^{(0)}_{q} - 2 \Delta\epsilon^{(0)}_R) v^{(1)}_{pR}(0) v^{(1)}_{qR}(0) e^{-s [ (\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_R)^2+ (\epsilon^{(0)}_q - \Delta\epsilon^{(0)}_R)^2]} \notag \\
&f^{(2)}_{pq}(s) = \notag \\ &f^{(2)}_{pq}(s) = \notag \\
&\color{red}{\boxed{\color{black}{- \sum_R\frac{ v^{(1)}_{pR}(0) v^{(1)}_{qR}(0) (\eps^{(0)}_{p} + \eps^{(0)}_{q} - 2 \Delta\eps^{(0)}_R)}{(\eps^{(0)}_p - \Delta\eps^{(0)}_R)^2+ (\eps^{(0)}_q - \Delta\eps^{(0)}_R)^2}(1 - e^{-s [ (\eps^{(0)}_p - \Delta\eps^{(0)}_R)^2+ (\eps^{(0)}_q - \Delta\eps^{(0)}_R)^2]})}}} \notag &\color{red}{\boxed{\color{black}{- \sum_R\frac{ v^{(1)}_{pR}(0) v^{(1)}_{qR}(0) (\epsilon^{(0)}_{p} + \epsilon^{(0)}_{q} - 2 \Delta\epsilon^{(0)}_R)}{(\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_R)^2+ (\epsilon^{(0)}_q - \Delta\epsilon^{(0)}_R)^2}(1 - e^{-s [ (\epsilon^{(0)}_p - \Delta\epsilon^{(0)}_R)^2+ (\epsilon^{(0)}_q - \Delta\epsilon^{(0)}_R)^2]})}}} \notag
\end{align} \end{align}
\begin{align} \begin{align}