SRGGW/Manuscript/SRGGW.tex

652 lines
39 KiB
TeX

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts,bbold}
\usepackage[version=4]{mhchem}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}
\usepackage[
colorlinks=true,
citecolor=blue,
breaklinks=true
]{hyperref}
\urlstyle{same}
\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{A Similarity Renormalization Group Approach To Many-Body Perturbation 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}
%=================================================================%
One-body Green's functions provide a natural and elegant way to access the charged excitations energies of a physical system. \cite{Martin_2016,Golze_2019}
The one-body non-linear Hedin's equations give a recipe to obtain the exact interacting one-body Green's function and therefore the exact ionization potentials and electron affinities. \cite{Hedin_1965}
Unfortunately, fully solving Hedin's equations is out of reach and one must resort to approximations.
In particular, the $GW$ approximation, \cite{Hedin_1965} which has first been mainly used in the context of solids \cite{Strinati_1980,Strinati_1982,Hybertsen_1985,Hybertsen_1986,Godby_1986,Godby_1987,Godby_1987a,Godby_1988,Blase_1995} and is now widely used for molecules as well \ant{ref?}, provides fairly accurate results for weakly correlated systems\cite{Hung_2017,vanSetten_2015,vanSetten_2018,Caruso_2016,Korbel_2014,Bruneval_2021} at a low computational cost. \cite{Foerster_2011,Liu_2016,Wilhelm_2018,Forster_2021,Duchemin_2021}
The $GW$ approximation is an approximation for the self-energy $\Sigma$ which role is to relate the exact interacting Green's function $G$ to a non-interacting reference one $G_0$ through the Dyson equation
\begin{equation}
\label{eq:dyson}
G = G_0 + G_0\Sigma G.
\end{equation}
The self-energy encapsulates all the Hartree-exchange-correlation effects which are not taken in account in the reference system.
%Throughout this manuscript the references are chosen to be the Hartree-Fock (HF) ones so that the self-energy only account for the missing correlation.
Approximating $\Sigma$ as the first order term of its perturbation expansion with respect to the screened interaction $W$ gives the so-called $GW$ approximation. \cite{Hedin_1965, Martin_2016}
Alternatively one could choose to define $\Sigma$ as the $n$-th order expansion in terms of the bare Coulomb interaction leading to the GF($n$) class of approximations. \cite{Hirata_2015,Hirata_2017}
The GF(2) approximation is also known as the second Born approximation. \ant{ref ?}
Despite a wide range of successes, many-body perturbation theory is not flawless.
It has been shown that a variety of physical quantities such as charged and neutral excitations energies or correlation and total energies computed within many-body perturbation theory exhibits some discontinuities. \cite{Veril_2018,Loos_2018b}
Even more worrying these discontinuities can happen in the weakly correlated regime where $GW$ is thought to be valid.
These discontinuities are due to a transfer of spectral weight between two solutions of the quasi-particle equation. \cite{Monino_2022}
This is another occurrence of the infamous intruder-state problem. \cite{Roos_1995,Olsen_2000,Choe_2001} \ant{more ref}
In addition, systems for which two quasi-particle solutions have a similar spectral weight are known to be particularly difficult to converge for partially self-consistent $GW$. \cite{Forster_2021}
In a recent study, Monino and Loos showed that these discontinuities could be removed by introduction of a regularizer inspired by the similarity renormalisation group (SRG) in the quasi-particle equation. \cite{Monino_2022}
Encouraged by this result, this work will investigate the application of the SRG formalism to many-body perturbation theory in its $GW$ and GF(2) variants.
The SRG has been developed independently by Wegner \cite{Wegner_1994} and Glazek and Wilson \cite{Glazek_1993,Glazek_1994} in the context of condensed matter systems and light-front quantum field theories, respectively.
This formalism has been been introduced in quantum chemistry by White \cite{White_2002} before being explored in more details by Evangelista and his co-workers in the context of multi-reference electron correlation theories. \cite{Evangelista_2014b,Li_2015, Li_2016, Li_2017, Li_2018, Li_2019a}
The SRG has also been successful in the context of nuclear theory, \cite{Bogner_2007,Tsukiyama_2011,Tsukiyama_2012,Hergert_2013,Hergert_2016} see Ref.\onlinecite{Hergert_2016a} for a recent review in this field. \ant{Maybe search for recent papers of T. Duguet as well.}
The SRG transformation aims at decoupling a reference space from an external space while folding information about the coupling in the reference space.
This is often during such decoupling that intruder states appear. \ant{ref}
However, SRG is particularly well-suited to avoid them because the decoupling of each external configuration is inversely proportional to its energy difference with the reference space.
Because intruder states have energies really close to the reference energies they will be the last ones decoupled.
Therefore the SRG continuous transformation can be stopped once every external configurations except the intruder ones have been decoupled.
Doing so, it gives a way to fold in information about the coupling in the reference space while avoiding intruder states.
The aim of this manuscript is to investigate whether SRG can treat the intruder-state problem in many-body perturbation theory as successfully as it has been in other fields.
We begin by reviewing the $GW$ formalism in Sec.~\ref{sec:gw} and then briefly review the SRG formalism in Sec.~\ref{sec:srg}.
Section~\ref{sec:theoretical} is concluded by a perturbative analysis of the SRG formalism applied to $GW$ (see Sec.~\ref{sec:srggw}).
The computational details of our implementation are provided in Sec.~\ref{sec:comp_det} before turning to the results section.
This section starts by
%=================================================================%
\section{Theoretical background}
\label{sec:theoretical}
%=================================================================%
%%%%%%%%%%%%%%%%%%%%%%
\subsection{Many-body perturbation theory in the GW approximation}
\label{sec:gw}
%%%%%%%%%%%%%%%%%%%%%%
The central equation of many-body perturbation theory based on Hedin's equations is the so-called quasi-particle equation
\begin{equation}
\label{eq:quasipart_eq}
\left[ \bF + \bSig(\omega = \epsilon_p) \right] \psi_p = \epsilon_p \psi_p,
\end{equation}
where $\bF$ is the Fock matrix, \cite{SzaboBook} and $\bSig(\omega)$ is the self-energy, both are $K \times K$ matrices with $K$ the number of one-body basis functions.
The self-energy can be physically understood as a dynamical screening correction to the Hartree-Fock (HF) problem represented by $\bF$.
Similarly to the HF case, this equation needs to be solved self-consistently.
Note that $\bSig$ is dynamical, \ie it depends on both the eigenvalues $\epsilon_p$ and eigenvectors $\psi_p$ while $\bF$ depends only on the eigenvectors.
Because of this $\omega$ dependence, fully solving this equation is a rather complicated task, hence several approximate solving schemes has been developed.
The most popular one is probably the one-shot scheme, known as $G_0W_0$ if the self-energy is the $GW$ one, in which the off-diagonal elements of Eq.~(\ref{eq:quasipart_eq}) are neglected and the self-consistency is abandoned.
In this case, there are $K$ quasi-particle equations which read as
tr\begin{equation}
\label{eq:G0W0}
\epsilon_p^{\HF} + \Sigma_{p}(\omega) - \omega = 0,
\end{equation}
where $\Sigma_{p}(\omega)$ are the diagonal elements of $\bSig$ and $\epsilon_p^{\HF}$ are the HF orbital energies.
The previous equation is non-linear with respect to $\omega$ and therefore can have multiple solutions $\epsilon_{p,s}$ for a given $p$.
These solutions can be characterised by their spectral weight defined as the renormalisation factor $Z_{p,s}$
\begin{equation}
\label{eq:renorm_factor}
0 \leq Z_{p,s} = \left[ 1 - \pdv{\Sigma_{p}(\omega)}{\omega}\bigg|_{\omega=\epsilon_{p,s}} \right]^{-1} \leq 1.
\end{equation}
The solution with the largest weight is referred to as the quasi-particle solution while the others are known as satellites or shake-up solutions.
However, in some cases Eq.~(\ref{eq:G0W0}) can have two (or more) solutions with similar weights and the quasi-particle solution is not well-defined.
In fact, these cases are related to the discontinuities and convergence problems discussed earlier because the additional solutions with large weights are the previously mentioned intruder states.
One obvious flaw of the one-shot scheme mentioned above is its starting point dependence.
Indeed, in Eq.~(\ref{eq:G0W0}) we chose to use the HF orbital energies but this is arbitrary and one could have chosen Kohn-Sham orbitals for example.
Therefore, one could try to optimise the starting point to obtain the best one-shot energies possible. \cite{Korzdorfer_2012,Marom_2012,Bruneval_2013,Gallandi_2015,Caruso_2016, Gallandi_2016}
Alternatively, one could solve this equation self-consistently leading to the eigenvalue only self-consistent scheme. \cite{Shishkin_2006,Blase_2011,Marom_2012,Kaplan_2016,Wilhelm_2016}
To do so the energy of the quasi-particle solution of the previous iteration is used to build Eq.~(\ref{eq:G0W0}) and then this equation is solved for $\omega$ again until convergence is reached.
However, if the quasi-particle solution is not well-defined, self-consistency can be quite difficult, if not impossible, to reach.
Even if self-consistency has been reached, the starting point dependence has not been totally removed because the results still depend on the starting molecular orbitals. \cite{Marom_2012}
To update both the energies and the molecular orbitals, one needs to take into account the off-diagonal elements in Eq.~(\ref{eq:quasipart_eq}).
To take into account the effect of off-diagonal elements without fully solving the quasi-particle equation, one can resort to the quasi-particle self-consistent (qs) scheme in which $\bSig(\omega)$ is replaced by a static approximation $\bSig^{\qs}$.
The algorithm to solve the qs problem is totally analog to the HF case with $\bF$ replaced by $\bF + \bSig^{\qs}$.
Various choice for $\bSig^\qs$ are possible but the most popular one is the following Hermitian approximation
\begin{equation}
\label{eq:sym_qsgw}
\Sigma_{pq}^\qs = \frac{1}{2}\Re\left(\Sigma_{pq}(\epsilon_p) + \Sigma_{pq}(\epsilon_q) \right).
\end{equation}
This form has first been introduced by Faleev and co-workers \cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007} before being derived as the effective Hamiltonian that minimizes the length of the gradient of the Klein functional for non-interacting Green's function. \cite{Ismail-Beigi_2017}
One of the main results of this manuscript is the derivation from first principles of an alternative static Hermitian form, this will be done in the next section.
In this case as well self-consistency can be difficult to reach in cases where multiple solutions have large spectral weights.
Multiple solutions arise due to the $\omega$ dependence of the self-energy.
Therefore, by suppressing this dependence the static qs approximation relies on the fact that there is one well-defined quasi-particle solution.
If it is not the case, the qs scheme will oscillates between the solutions with large weights. \cite{Forster_2021}
Convergence problems arise when a shake-up configuration has an energy similar to the associated quasi-particle solution, this satellite is the so-called intruder state.
The intruder state problem can be dealt with by introducing \textit{ad hoc} regularizers.
The $\ii eta$ term that is usually added in the denominator of the self-energy (see below) is the usual imaginary-shift regularizer used in various other theories flawed by intruder states. \cite{Battaglia_2022} \ant{more ref...}
Various other regularizers are possible and in particular one of us has shown that a regularizer inspired by the SRG had some advantages, in the $GW$ case, over the imaginary shift one. \cite{Monino_2022}
But it would be more rigorous, and more instructive, to obtain this regularizer from first principles by applying the SRG formalism to many-body perturbation theory.
This is the aim of this work.
Therefore if we apply it, the SRG would gradually remove the coupling between the quasi-particle and the satellites resulting in a renormalized quasi-particle.
However, to do so one needs to identify the coupling terms in Eq.~(\ref{eq:quasipart_eq}), which is not straightforward.
The way around this problem is to transform Eq.~(\ref{eq:quasipart_eq}) to its upfolded version and the coupling terms will elegantly appear in the process.
From now on, we will restrict ourselves to the $GW$ in the Tamm-Dancoff approximation (TDA) case but the same derivation could be done for the non-TDA $GW$ and GF(2) self-energies.
The corresponding formula are given in Appendix~\ref{sec:nonTDA} and \ref{sec:GF2}, respectively.
The upfolded $GW$ quasi-particle equation is the following
\begin{equation}
\label{eq:GWlin}
\begin{pmatrix}
\bF & \bV^{\text{2h1p}} & \bV^{\text{2p1h}} \\
(\bV^{\text{2h1p}})^{\mathrm{T}} & \bC^{\text{2h1p}} & \bO \\
(\bV^{\text{2p1h}})^{\mathrm{T}} & \bO & \bC^{\text{2p1h}} \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX \\
\bY^{\text{2h1p}} \\
\bY^{\text{2p1h}} \\
\end{pmatrix}
=
\begin{pmatrix}
\bX \\
\bY^{\text{2h1p}} \\
\bY^{\text{2p1h}} \\
\end{pmatrix}
\cdot
\boldsymbol{\epsilon},
\end{equation}
where $\boldsymbol{\epsilon}$ is a diagonal matrix collecting the quasi-particle and satellite energies, the 2h1p and 2p1h matrix elements are
\begin{subequations}
\begin{align}
C^\text{2h1p}_{ija,klc} & = \qty[ \qty( \epsilon_{i}^{\GW} + \epsilon_{j}^{\GW} - \epsilon_{a}^{\GW}) \delta_{jl} \delta_{ac} - \eri{jc}{al} ] \delta_{ik} ,
\\
C^\text{2p1h}_{iab,kcd} & = \qty[ \qty( \epsilon_{a}^{\GW} + \epsilon_{b}^{\GW} - \epsilon_{i}^{\GW}) \delta_{ik} \delta_{ac} + \eri{ak}{ic} ] \delta_{bd},
\end{align}
\end{subequations}
and the corresponding coupling blocks read
\begin{align}
V^\text{2h1p}_{p,klc} & = \eri{pc}{kl},
&
V^\text{2p1h}_{p,kcd} & = \eri{pk}{dc}.
\end{align}
Throughout the manuscript $p,q,r,s$ indices are used for general orbitals while $i,j,k,l$ and $a,b,c,d$ refers to occupied and virtual orbitals, respectively.
The indices $v$ and $w$ will be used for neutral excitations, \ie composite indices $v=(ia)$.
The usual $GW$ non-linear equation can be obtained by applying L\"odwin partitioning technique to Eq.~(\ref{eq:GWlin}) \cite{Lowdin_1963,Bintrim_2021}
\begin{equation}
\label{eq:GWnonlin}
\left( \bF + \bSig(\omega) \right) \bX = \omega \bX,
\end{equation}
with
\begin{align}
\bSig(\omega) &= \bV^{\hhp} \left(\omega \mathbb{1} - \bC^{\hhp}\right)^{-1} (\bV^{\hhp})^{\mathsf{T}} \\
&+ \bV^{\pph} \left(\omega \mathbb{1} - \bC^{\pph})^{-1} (\bV^{\pph}\right)^{\mathsf{T}}, \notag
\end{align}
which can be further developed to give
\begin{equation}
\label{eq:GW_selfenergy}
\Sigma_{pq}(\omega)
= \sum_{iv} \frac{W_{p,(i,v)} W_{q,(i,v)}}{\omega - \epsilon_i + \Omega_{v} - \ii \eta}
+ \sum_{av} \frac{W_{p,(a,v)}W_{q,(a,v)}}{\omega - \epsilon_a - \Omega_{v} + \ii \eta},
\end{equation}
with the screened integrals defined as
\begin{equation}
\label{eq:GW_sERI}
W_{p,(q,v)} = \sum_{ia}\eri{pi}{qa}\qty( \bX_{v})_{ia},
\end{equation}
where $\bX$ is the matrix of eigenvectors of the particle-hole direct RPA (dRPA) problem in the TDA defined as
\begin{equation}
\bA \bX = \boldsymbol{\Omega} \bX,
\end{equation}
with
\begin{equation}
A^\dRPA_{ij,ab} = (\epsilon_i - \epsilon_a) \delta_{ij}\delta_{ab} + \eri{ib}{aj}.
\end{equation}
$\boldsymbol{\Omega}$ is the diagonal matrix of eigenvalues and its elements $\Omega_v$ appear in Eq.~(\ref{eq:GW_selfenergy}).
Equations~(\ref{eq:GWlin}) and~(\ref{eq:GWnonlin}) have exactly the same solutions but one is linear and the other not.
The price to pay for this linearity is that the size of the matrix in the former equation is $\mathcal{O}(K^3)$ while it is $\mathcal{O}(K)$ in the latter one.
We refer to Ref.~\onlinecite{Bintrim_2021} for a detailed discussion of the up/downfolding processes of the $GW$ equations (see also Chapter 8 of Ref.~\onlinecite{Schirmer_2018} for the GF(2) case).
As can be readily seen in Eq.~\eqref{eq:GWlin}, the blocks $V^\text{2h1p}$ and $ V^\text{2p1h}$ are coupling the 1h and 1p configuration to the dressed 2h1p and 2p1h configurations.
Therefore, these blocks will be the target of our SRG transformation but before going in more details we will review the SRG formalism.
%%%%%%%%%%%%%%%%%%%%%%
\subsection{The similarity renormalization group}
\label{sec:srg}
%%%%%%%%%%%%%%%%%%%%%%
The similarity renormalization group method aims at continuously transforming an Hamiltonian to a diagonal form, or more often to a block-diagonal form.
Therefore, the transformed Hamiltonian
\begin{equation}
\label{eq:SRG_Ham}
\bH(s) = \bU(s) \, \bH \, \bU^\dag(s),
\end{equation}
depends on a flow parameter $s$, such that $\bH(s=0)$ is the initial untransformed Hamiltonian and $\bH(s=\infty)$ is the (block)-diagonal Hamiltonian.
An evolution equation for $\bH(s)$ can be easily obtained by deriving Eq~(\ref{eq:SRG_Ham}) with respect to $s$.
This gives the flow equation
\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 defining such an approximation, we need to define what are the blocks to suppress to obtain a block-diagonal Hamiltonian.
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}
and by definition we have the following condition on $\bH^\text{od}$
\begin{equation}
\bH^\text{od}(s=\infty) = \boldsymbol{0}.
\end{equation}
In this work, we will use Wegner's canonical generator which is defined as \cite{Wegner_1994}
\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}
If this generator is used, the following condition is verified \cite{Kehrein_2006}
\begin{equation}
\label{eq:derivative_trace}
\dv{s}\text{Tr}\left[ \bH^\text{od}(s)^2 \right] \leq 0
\end{equation}
which implies that the matrix elements of the off-diagonal part will decrease in a monotonic way.
Even more, the coupling coefficients associated with the highest energy determinants are removed first as will be evidenced by the perturbative analysis after.
The main flaw of this generator is that it generates a stiff set of ODE which is therefore difficult to solve numerically. \ant{ref}
However, here we will not tackle the full SRG problem but only consider analytical low-order perturbative expressions so we will not be affected by this problem. \cite{Evangelista_2014, Hergert_2016}
Let us now perform the perturbative analysis of the SRG equations.
For $s=0$, the initial problem is
\begin{equation}
\bH(0) = \bH^\text{d}(0) + \lambda ~ \bH^\text{od}(0)
\end{equation}
where $\lambda$ is the usual perturbation parameter and the off-diagonal part of the Hamiltonian has been defined as the perturber.
For finite values of $s$, we have the following perturbation expansion of the Hamiltonian
\begin{equation}
\label{eq:perturbation_expansionH}
\bH(s) = \bH^{(0)}(s) + \lambda ~ \bH^{(1)}(s) + \lambda^2 \bH^{(2)}(s) + \dots
\end{equation}
Hence, the generator $\boldsymbol{\eta}(s)$ admits a perturbation expansion as well.
Then, one can collect order by order the terms in Eq.~(\ref{eq:flowEquation}) and solve analytically the low-order differential equations.
%%%%%%%%%%%%%%%%%%%%%%
\subsection{Renormalised GW}
\label{sec:srggw}
%%%%%%%%%%%%%%%%%%%%%%
Finally, the SRG formalism exposed above will be applied to $GW$.
First, one needs to defined the diagonal and off-diagonal parts of the $GW$ effective Hamiltonian.
As hinted at the end of section~\ref{sec:gw}, the diagonal and off-diagonal parts will be defined as
\begin{align}
\label{eq:diag_and_offdiag}
\bH^\text{d}(s) &=
\begin{pmatrix}
\bF & \bO & \bO \\
\bO & \bC^{\text{2h1p}} & \bO \\
\bO & \bO & \bC^{\text{2p1h}} \\
\end{pmatrix}
&
\bH^\text{od}(s) &=
\begin{pmatrix}
\bO & \bV^{\text{2h1p}} & \bV^{\text{2p1h}} \\
(\bV^{\text{2h1p}})^{\mathrm{T}} & \bO & \bO \\
(\bV^{\text{2p1h}})^{\mathrm{T}} & \bO & \bO \\
\end{pmatrix}
\end{align}
where we have omitted the $s$ dependence of the matrix elements for the sake of brevity.
Then, the aim of this section is to solve order by order the flow equation [see Eq.~(\ref{eq:flowEquation})] knowing that the initial conditions 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}
where we have defined the matrix $\bC$ and $\bV$ that collects the 2h1p and 2p1h channels for the sake of conciseness.
Then, the perturbative expansions can be inserted in Eq.~(\ref{eq:GWlin}) before downfolding to obtain a renormalised quasi-particle equation.
In particular, in this manuscript the focus will be on the second-order renormalized quasi-particle equation.
%///////////////////////////%
\subsubsection{Zero-th order matrix elements}
%///////////////////////////%
There is only one zero-th order term in the right hand side of the flow equation
\begin{equation}
\dv{\bH^{(0)}}{s} = \comm{\comm{\bHd{0}}{\bHod{0}}}{\bH^{(0)}},
\end{equation}
and performing the block matrix products gives the following system of equations
\begin{subequations}
\begin{align}
\dv{\bF^{(0)}}{s} &= \bO \\
\dv{\bC^{(0)}}{s} &= \bO \\
\dv{\bV^{(0),\dagger}}{s} &= 2 \bC^{(0)}\bV^{(0),\dagger}\bF^{(0)} - \bV^{(0),\dagger}(\bF^{(0)})^2 - (\bC^{(0)})^2\bV^{(0),\dagger} \\
\dv{\bV^{(0)}}{s} &= 2 \bF^{(0)}\bV^{(0)}\bC^{(0)} - (\bF^{(0)})^2\bV^{(0)} - \bV^{(0)}(\bC^{(0)})^2
\end{align}
\end{subequations}
where the $s$ dependence of $\bV^{(0)}$ and $\bV^{(0),\dagger}$ has been dropped in the last two equations.
$\bF^{(0)}$ and $\bC^{(0)}$ do not depend on $s$ as a consequence of the first two equations.
The last equation can be solved by introducing $\bU$ the matrix that diagonalizes $\bC^{(0)} = \bU \bD^{(0)} \bU^{-1}$ such that the differential equation for $\bV^{(0)}$ becomes
\begin{equation}
\dv{\bW^{(0)}}{s} = 2 \bF^{(0)}\bW^{(0)} \bD^{(0)} - (\bF^{(0)})^2\bW^{(0)} - \bW^{(0)} (\bD^{(0)})^2
\end{equation}
where $\bW^{(0)}= \bV^{(0)} \bU$.
Note that the matrix $\bU$ is also used in the downfolding process of Eq.~(\ref{eq:GWlin}). \cite{Bintrim_2021}
Due to the diagonal structure of $\bF^{(0)}$ and $\bD^{(0)}$, these equations can be easily solved and give
\begin{equation}
W_{p,(q,v)}^{(0)}(s) = W_{p,(q,v)}^{(0)}(0) e^{- (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 s}
\end{equation}
Due to the initial conditions $\bV^{(0)}(0) = \bO$, we have $\bW^{(0)}(s)=\bO$ and therefore $\bV^{(0)}(s)=\bO=\bV^{(0)}(0) $.
The two first equations of the system are trivial and finally we have
\begin{equation}
\bH^{(0)}(s) = \bH^{(0)}(0)
\end{equation}
which shows that the zero-th order matrix elements are independent of $s$.
The matrix elements of $\bU$ and $\bD$ are
\begin{align}
U_{(p,v),(q,w)}^{(0)} &= \delta_{pq} \bX_{v,w} \\
D_{(p,v),(q,w)}^{(0)} &= \left(\epsilon_p + \text{sign}(\epsilon_p-\epsilon_F)\Omega_v\right)\delta_{pq}\delta_{vw}
\end{align}
where $\epsilon_F$ is the Fermi level.
%///////////////////////////%
\subsubsection{First order matrix elements}
%///////////////////////////%
Knowing that $\bHod{0}(s)=\bO$, the first order flow equation is
\begin{equation}
\dv{\bH^{(1)}}{s} = \comm{\comm{\bHd{0}}{\bHod{1}}}{\bHd{0}}
\end{equation}
which gives the same system of equations as in the previous subsection except that $\bV^{(0)}$ and $\bV^{(0),\dagger}$ should be replaced by $\bV^{(1)}$ and $\bV^{(1),\dagger}$.
Once again the two first equations are easily solved
\begin{align}
\bF^{(1)}(s) &= \bF^{(1)}(0) = \bO & \bC^{(1)}(s) &= \bC^{(1)}(0) = \bO.
\end{align}
and the first order coupling elements are given by (up to a multiplication by $\bU^{-1}$)
\begin{align}
W_{p,(q,v)}^{(1)}(s) &= W_{p,(q,v)}^{(1)}(0) e^{- (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 s} \\
W_{p,(q,v)}^{(1)}(s) &= \left( \sum_{ia}\eri{pi}{qa}\qty( \bX_{v})_{ia} \right) e^{- (\epsilon_p - \epsilon_q - \text{sign}(\epsilon_q-\epsilon_F)\Omega_v)^2 s}
\end{align}
Note that at $s=0$ the elements $W_{p,(q,v)}^{(1)}(0)$ are equal to the two-electron screened integrals defined in Eq.~(\ref{eq:GW_sERI}) and that for $s\to\infty$ they go to zero.
Therefore, $W_{p,(q,v)}^{(1)}(s)$ are renormalized two-electrons screened integrals.
Note the close similarity of the first-order element expressions with the ones of Evangelista in Ref.~\onlinecite{Evangelista_2014b} obtained in a second quantization formalism (see also Ref.~\onlinecite{Hergert_2016}).
%///////////////////////////%
\subsubsection{Second-order matrix elements}
% ///////////////////////////%
The second-order renormalised quasi-particle equation is given by
\begin{equation}
\label{eq:GW_renorm}
\left( \widetilde{\bF}(s) + \widetilde{\bSig}(\omega; s) \right) \bX = \omega \bX
\end{equation}
with
\begin{align}
\widetilde{\bF}(s) &= \bF^{(0)}+\bF^{(2)}(s)\\
\label{eq:srg_sigma}
\widetilde{\bSig}(\omega; s) &= \bV^{(1)}(s) \left(\omega \mathbb{1} - \bC^{(0)}\right)^{-1} (\bV^{(1)}(s))^{\dagger}
\end{align}
As can be readily seen above, $\bF^{(2)}$ is the only second-order block of the effective Hamiltonian contributing to the second-order SRG quasi-particle equation.
Collecting every second-order terms and performing the block matrix products results in the following differential equation for $\bF^{(2)}$
\begin{equation}
\label{eq:diffeqF2}
\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} .
\end{equation}
This can be solved by simple integration along with the initial condition $\bF^{(2)}=\bO$ to give
\begin{equation}
F_{pq}^{(2)}(s) = \sum_{r,v} \frac{\Delta_{prv}+ \Delta_{qrv}}{\Delta_{prv}^2 + \Delta_{qrv}^2} W_{p,(r,v)} W^{\dagger}_{(r,v),q}\left(1 - e^{-(\Delta_{prv}^2 + \Delta_{qrv}^2) s}\right).
\end{equation}
with $\Delta_{pqv} = \epsilon_p - \epsilon_q - \text{sign}(\epsilon_q-\epsilon_F)\Omega_v$.
At $s=0$, this second-order correction is null while for $s\to\infty$ it tends towards the following static limit
\begin{equation}
\label{eq:static_F2}
F_{pq}^{(2)}(\infty) = \frac{\Delta_{prv}+ \Delta_{qrv}}{\Delta_{prv}^2 + \Delta_{qrv}^2} W_{p,(r,v)} W^{\dagger}_{(r,v),q}.
\end{equation}
Note that in the $s\to\infty$ limit the dynamic part of the self-energy [see Eq.~\eqref{eq:srg_sigma}] tends to zero.
Therefore, the SRG flow gradually transforms the dynamic degrees of freedom of $\bSig(\omega)$ in static ones, starting from the ones that have the largest denominators in Eq.~(\ref{eq:static_F2}).
Interestingly, the static limit, \ie $s\to\infty$ limit, of Eq.~(\ref{eq:GW_renorm}) defines an alternative qs$GW$ approximation to the one defined by Eq.~(\ref{eq:sym_qsgw}).
Yet, both are closely related as they share the same diagonal terms.
Also, note that the hermiticity is naturally enforced in the SRG static approximation as opposed to the symmetrized case.
However, as will be discussed in more details in the results section, the convergence of the qs$GW$ scheme using $\widetilde{\bF}(\infty)$ is very poor.
This is similar to the symmetric case when the imaginary shift $\ii \eta$ is set to zero.
Indeed, in qs$GW$ calculation using the symmetrized static form, increasing $\eta$ to ensure convergence in difficult cases is most often unavoidable. \ant{ref fabien ?}
Therefore, we will define the SRG-qs$GW$ static effective Hamiltonian as
\begin{align}
\label{eq:SRG_qsGW}
\Sigma_{pq}^{\text{SRG}}(s) &= \epsilon_p \delta_{pq} + \sum_{r,v} \frac{\left(\epsilon_{p} + \epsilon_{q} - 2 (\epsilon_r \pm \Omega_v)\right) W_{p,(r,v)} W_{q,(r,v)}}{(\epsilon_p - \epsilon_r \pm \Omega_v)^2 + (\epsilon_q - \epsilon_r \pm \Omega_v)^2} \notag \\
&\times \left(1 - e^{-(\epsilon_p - \epsilon_r \pm \Omega_v)^2s} e^{-(\epsilon_q - \epsilon_r \pm \Omega_v)^2s}\right)
\end{align}
which depends on one regularising parameter $s$ analogously to $eta$ in the usual case.
The fact that the $s\to\infty$ static limit does not well converge when used in a qs$GW$ calculation could have been predicted because in this limit even the intruder states have been included in $\tilde{\bF}$.
Therefore, we should use a value of $s$ large enough to include almost every states but small enough to avoid intruder states.
To conclude this section, we will discuss the case of discontinuities.
Indeed, we have previously said that intruder states are responsible for both the poor convergence of qs$GW$ and discontinuities in physical quantities at the $\GOWO$ level.
So is it possible to remove discontinuities by using the SRG machinery developed above?
In fact, not directly because discontinuities are due to intruder states in the dynamic part while we have seen just above that a finite value of $s$ is well-designed to avoid the intruder states in the static part.
However, doing a change of variable such that
\begin{align}
e^{-s} &= 1-e^{-t} & 1 - e^{-s} &= e^{-t}
\end{align}
hence choosing a finite value of $t$ is well-designed to avoid discontinuities in the dynamic.
In fact, the dynamic part after the change of variable is closely related to the SRG-inspired regularizer introduced by Monino and Loos in Ref.~\onlinecite{Monino_2022}.
%=================================================================%
\section{Computational details}
\label{sec:comp_det}
% =================================================================%
The two qs$GW$ variants considered in this work has been implemented in a in-house program.
The $GW$ implementation closely follows the one of mol$GW$. \cite{Bruneval_2016}
The geometry have been optimized at the CC3 level in the aug-cc-pvtz basis set without frozen core using the CFOUR program.
The reference CCSD(T) IP energies have obtained using default parameters of Gaussian 16.
This means that the cations used an unrestricted HF reference while the neutral ground-state energies have been obtained in a restricted formalism.
%=================================================================%
\section{Results}
\label{sec:results}
%=================================================================%
%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
%%%%%%%%%%%%%%%%%%%%%%
Here comes the conclusion.
%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Data availability statement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The data that supports the findings of this study are available within the article.% and its supplementary material.
%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{SRGGW}
%%%%%%%%%%%%%%%%%%%%%%%%
\appendix
%%%%%%%%%%%%%%%%%%%%%%
\section{Non-TDA $GW$ equations}
\label{sec:nonTDA}
%%%%%%%%%%%%%%%%%%%%%%
The $GW$ self-energy without TDA is the same as in Eq.~(\ref{eq:GW_selfenergy}) but the screened integrals are now defined as
\begin{equation}
\label{eq:GWnonTDA_sERI}
W_{p,(q,v)} = \sum_{ia}\eri{pi}{qa}\qty( \bX_{v} + \bY_{v})_{ia},
\end{equation}
where $\bX$ and $\bY$ are the matrix of eigenvectors of the full particle-hole dRPA problem defined as
\begin{equation}
\label{eq:full_dRPA}
\begin{pmatrix}
\bA & \bB \\
- \bB & \bA \\
\end{pmatrix}
\begin{pmatrix}
\bX \\
\bY \\
\end{pmatrix} = \boldsymbol{\Omega}
\begin{pmatrix}
\bX \\
\bY \\
\end{pmatrix},
\end{equation}
with
\begin{align}
A_{ij,ab} &= (\epsilon_i - \epsilon_a) \delta_{ij}\delta_{ab} + \eri{ib}{aj}, \\
B_{ij,ab} &= \eri{ij}{ab}.
\end{align}
$\boldsymbol{\Omega}$ is the diagonal matrix of eigenvalues. Note that $\boldsymbol{\Omega}$ in this case has the same size as in the TDA because we consider only the positive excitations of the full dRPA problem.
Defining an unfold version of this equation which does not require a diagonalization of the RPA problem before unfolding is a tricky task (see supplementary material of Ref.~\onlinecite{Bintrim_2021}).
However, because we will eventually downfold again the upfolded matrix, we can use the following matrix \cite{Tolle_2022}
\begin{equation}
\label{eq:nonTDA_upfold}
\begin{pmatrix}
\bF & \bW^{\text{2h1p}} & \bW^{\text{2p1h}} \\
(\bW^{\text{2h1p}})^{\mathrm{T}} & \bD^{\text{2h1p}} & \bO \\
(\bW^{\text{2p1h}})^{\mathrm{T}} & \bO & \bD^{\text{2p1h}} \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX \\
\bY^{\text{2h1p}} \\
\bY^{\text{2p1h}} \\
\end{pmatrix}
=
\begin{pmatrix}
\bX \\
\bY^{\text{2h1p}} \\
\bY^{\text{2p1h}} \\
\end{pmatrix}
\cdot
\boldsymbol{\epsilon},
\end{equation}
which already depends on the screened integrals and therefore require the knowledge of the eigenvectors of the dRPA problem defined in Eq.~(\ref{eq:full_dRPA}).
where $\boldsymbol{\epsilon}$ is a diagonal matrix collecting the quasi-particle and satellite energies, the 2h1p and 2p1h matrix elements are
\begin{subequations}
\begin{align}
D^\text{2h1p}_{ija,klc} & = \delta_{ik}\delta_{jl} \delta_{ac} \qty[\epsilon_i - \Omega_{ja}] ,
\\
D^\text{2p1h}_{iab,kcd} & = \delta_{ik}\delta_{ac} \delta_{bd} \qty[\epsilon_a + \Omega_{ib}] ,
\end{align}
\end{subequations}
and the corresponding coupling blocks read
\begin{align}
W^\text{2h1p}_{p,klc} & = \sum_{ia}\eri{pi}{ka} \qty( \bX_{lc} + \bY_{lc})_{ia} \\
W^\text{2p1h}_{p,kcd} & = \sum_{ia}\eri{pi}{ca} \qty( \bX_{kd} + \bY_{kd})_{ia}
\end{align}
Using the SRG on this matrix instead of Eq.~(\ref{eq:GWlin}) gives the same expression for $\bW^{(1)}$, $\bF^{(2)}$ and $\bSig^{\text{SRG}}$ but now the screened integrals are the one of Eq.~\eqref{eq:GWnonTDA_sERI} and the eigenvalues $\Omega$ and eigenvectors $\bX$ and $\bY$ are the ones of the full RPA problem defined in Eq.~\eqref{eq:full_dRPA}.
%%%%%%%%%%%%%%%%%%%%%%
\section{GF(2) equations}
\label{sec:GF2}
%%%%%%%%%%%%%%%%%%%%%%
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 matrix elements of its second-order version read as
\begin{align}
\label{eq:GF2_selfenergy}
\Sigma_{pq}^{\text{GF(2)}}(\omega)
&= \frac{1}{\sqrt{2}} \sum_{ija} \frac{\aeri{pa}{ij}\aeri{qa}{ij}}{\omega + \epsilon _a -\epsilon_i -\epsilon_j - \ii \eta} \\
&+ \frac{1}{\sqrt{2}} \sum_{iab} \frac{\aeri{pi}{ab}\aeri{qi}{ab}}{\omega + \epsilon _i -\epsilon_a -\epsilon_b + \ii \eta}
\end{align}
This self-energy can be upfolded similarly to the $GW$ case and one obtain the following ``super-matrix''
\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{}{}$ is given below.
\begin{align}
\label{eq:GF2_unfolded}
V^\text{2h1p}_{p,ija} & = \frac{1}{\sqrt{2}}\aeri{pa}{ij}
\\
V^\text{2p1h}_{p,iab} & = \frac{1}{\sqrt{2}}\aeri{pi}{ab}
\\
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}
Note that this matrix is exactly the ADC(2) matrix for charged excitations.
The fact that the integrals are not screened in GF(2) manifests itself in the fact that the $\bC$ matrices are already diagonal.
Applying the SRG formalism to this matrix is completely analog to the derivation exposed in the main text.
We only give the analytical expressions of the matrix elements needed for the second-order SRG-GF(2) quasiparticle equations.
\begin{equation}
(V^\text{2h1p}_{p,ija})^{(1)}(s) = \frac{1}{\sqrt{2}}\aeri{pa}{ij} e^{- (\epsilon_p + \epsilon_a - \epsilon_i - \epsilon_j)^2 s}
\end{equation}
\begin{equation}
(V^\text{2h1p}_{p,iab})^{(1)}(s) = \frac{1}{\sqrt{2}}\aeri{pi}{ab} e^{- (\epsilon_p + \epsilon_i - \epsilon_a - \epsilon_b)^2 s}
\end{equation}
We define $ \Delta_{pq,rs} = \epsilon_p + \epsilon_q - \epsilon_r - \epsilon_s $
\begin{align}
F_{pq}^{(2)}(s) &= \sum_{ria} \frac{\epsilon_{p} + \epsilon_{q} - 2 (\epsilon_r \pm \Omega_v)}{(\epsilon_p - \epsilon_r \pm \Omega_v)^2 + (\epsilon_q - \epsilon_r \pm \Omega_v)^2} W_{p,(r,v)} \notag \\
&\times W^{\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).
\end{align}
\end{document}