The family of Green's function methods based on the $GW$ approximation has gained popularity in the electronic structure theory thanks to its accuracy in weakly correlated systems combined with its cost-effectiveness.
Despite this, self-consistent versions still pose challenges in terms of convergence.
A recent study \href{https://doi.org/10.1063/5.0089317}{[J. Chem. Phys. 156, 231101 (2022)]} has linked these convergence issues to the intruder-state problem.
In this work, a perturbative analysis of the similarity renormalization group (SRG) approach is performed on Green's function methods.
The SRG formalism enables us to derive, from first principles, the expression of a new, naturally Hermitian form of the static self-energy that can be employed in quasiparticle self-consistent $GW$ (qs$GW$) calculations.
The resulting SRG-based regularized self-energy significantly accelerates the convergence of qs$GW$ calculations.
One-body Green's functions provide a natural and elegant way to access the charged excitation energies of a physical system. \cite{CsanakBook,FetterBook,Martin_2016,Golze_2019}
The non-linear Hedin equations consist of a closed set of equations leading to the exact interacting one-body Green's function and, therefore, to a wealth of properties such as the total energy, density, ionization potentials, electron affinities, as well as spectral functions, without the explicit knowledge of the wave functions associated with the neutral and charged states of the system. \cite{Hedin_1965}
Unfortunately, solving exactly Hedin's equations is usually out of reach and one must resort to approximations.
In particular, the $GW$ approximation, \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019,Bruneval_2021} which has been first introduced in the context of solids \cite{Strinati_1980,Strinati_1982a,Strinati_1982b,Hybertsen_1985,Hybertsen_1986,Godby_1986,Godby_1987,Godby_1987a,Godby_1988,Blase_1995} and is now widely applied to molecular systems, \cite{Rohlfing_1999a,Horst_1999,Puschnig_2002,Tiago_2003,Rocca_2010,Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018,Blase_2018,Liu_2020,Li_2017,Li_2019,Li_2020,Li_2021,Blase_2020,Holzer_2018a,Holzer_2018b,Loos_2020e,Loos_2021,McKeon_2022} yields accurate charged excitation energies for weakly correlated systems \cite{Hung_2017,vanSetten_2015,vanSetten_2018,Caruso_2016,Korbel_2014,Bruneval_2021} at a relatively low computational cost. \cite{Foerster_2011,Liu_2016,Wilhelm_2018,Forster_2021,Duchemin_2019,Duchemin_2020,Duchemin_2021}
The $GW$ method approximates the self-energy $\Sigma$ which relates the exact interacting Green's function $G$ to a non-interacting reference version $G_0$ through a Dyson equation of the form
Approximating $\Sigma$ as the first-order term of its perturbative expansion with respect to the screened Coulomb potential $W$ yields the so-called $GW$ approximation \cite{Hedin_1965,Martin_2016}
Diagrammatically, $GW$ corresponds to a resummation of the (time-dependent) direct ring diagrams via the computation of the random-phase approximation (RPA) polarizability \cite{Ren_2012,Chen_2017} and is thus particularly well suited for weak correlation.
Despite a wide range of successes, many-body perturbation theory has well-documented limitations. \cite{Kozik_2014,Stan_2015,Rossi_2015,Tarantino_2017,Schaefer_2013,Schaefer_2016,Gunnarsson_2017,vanSetten_2015,Maggio_2017a,Duchemin_2020}
For example, modeling core-electron spectroscopy requires core ionization energies which have been proven to be challenging for routine $GW$ calculations. \cite{vanSetten_2018,Golze_2018,Golze_2020,Li_2022}
Many-body perturbation theory can also be used to access optical excitation energies through the Bethe-Salpeter equation. \cite{Salpeter_1951,Strinati_1988,Blase_2018,Blase_2020} However, the accuracy is not yet satisfying for triplet excited states, where instabilities often occur. \cite{Bruneval_2015,Jacquemin_2017a,Jacquemin_2017b,Holzer_2018a}
Therefore, even if $GW$ offers a good trade-off between accuracy and computational cost, some situations might require higher precision.
Unfortunately, defining a systematic way to go beyond $GW$ via the inclusion of vertex corrections has been demonstrated to be a tricky task. \cite{Baym_1961,Baym_1962,DeDominicis_1964a,DeDominicis_1964b,Bickers_1989a,Bickers_1989b,Bickers_1991,Hedin_1999,Bickers_2004,Shirley_1996,DelSol_1994,Schindlmayr_1998,Morris_2007,Shishkin_2007b,Romaniello_2009a,Romaniello_2012,Gruneis_2014,Hung_2017,Maggio_2017b,Mejuto-Zaera_2022}
For example, Lewis and Berkelbach have shown that naive vertex corrections can even worsen the quasiparticle energies with respect to $GW$. \cite{Lewis_2019}
We refer the reader to the recent review by Golze and co-workers (see Ref.~\onlinecite{Golze_2019}) for an extensive list of current challenges in many-body perturbation theory.
Many-body perturbation theory also suffers from the infamous intruder-state problem,\cite{Andersson_1994,Andersson_1995a,Roos_1995,Forsberg_1997,Olsen_2000,Choe_2001} where they manifest themselves as solutions of the quasiparticle equation with non-negligible spectral weights.
These multiple solutions hinder the convergence of partially self-consistent schemes, \cite{Veril_2018,Forster_2021,Monino_2022} such as quasiparticle self-consistent $GW$\cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016} (qs$GW$) and eigenvalue-only self-consistent $GW$\cite{Shishkin_2007a,Blase_2011b,Marom_2012,Kaplan_2016,Wilhelm_2016} (ev$GW$).
The simpler one-shot $G_0W_0$ scheme \cite{Strinati_1980,Hybertsen_1985a,Hybertsen_1986,Godby_1988,Linden_1988,Northrup_1991,Blase_1994,Rohlfing_1995,Shishkin_2007a} is also impacted by these intruder states, leading to discontinuities in a variety of physical quantities including charged and neutral excitation energies, correlation and total energies.\cite{Loos_2018b,Veril_2018,Loos_2020e,Berger_2021,DiSabatino_2021,Monino_2022,Scott_2023}
In a recent study, Monino and Loos showed that the discontinuities could be removed by the introduction of a regularizer inspired by the similarity renormalization group (SRG) in the quasiparticle equation. \cite{Monino_2022}
Encouraged by the recent successes of regularization schemes in many-body quantum chemistry methods, such as in single- and multi-reference perturbation theory, \cite{Lee_2018a,Shee_2021,Evangelista_2014b,ChenyangLi_2019a,Battaglia_2022} the present work investigates the application of the SRG formalism in $GW$-based methods.
In particular, we focus here on the possibility of curing the qs$GW$ convergence problems using the SRG.
The SRG formalism 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 introduced in quantum chemistry by White \cite{White_2002} before being explored in more detail by Evangelista and coworkers in the context of multi-reference electron correlation theories. \cite{Evangelista_2014b,ChenyangLi_2015, ChenyangLi_2016,ChenyangLi_2017,ChenyangLi_2018,ChenyangLi_2019a,Zhang_2019,ChenyangLi_2021,Wang_2021,Wang_2023}
The SRG has also been successful in the context of nuclear structure theory, where it was first developed as a mature computational tool thanks to the work of several research groups.
The SRG transformation aims at decoupling an internal (or reference) space from an external space while incorporating information about their coupling in the reference space.
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.
By definition, intruder states have energies that are close to the reference energy, and therefore are the last to be decoupled.
By stopping the SRG transformation once all external configurations except the intruder states have been decoupled,
correlation effects between the internal and external spaces can be incorporated (or folded) without the presence of intruder states.
The goal of this manuscript is to determine if the SRG formalism can effectively address the issue of intruder states in many-body perturbation theory, as it has in other areas of electronic and nuclear structure theory.
This open question will lead us to an intruder-state-free static approximation of the self-energy derived from first-principles that can be employed in partially self-consistent $GW$ calculations.
Note that throughout the manuscript we focus on the $GW$ approximation but the subsequent derivations can be straightforwardly applied to other self-energies such as the one derived from the GF(2) \cite{Casida_1989,Casida_1991,SzaboBook,Stefanucci_2013,Ortiz_2013,Phillips_2014,Phillips_2015,Rusakov_2014,Rusakov_2016,Hirata_2015,Hirata_2017,Backhouse_2021,Backhouse_2020b,Backhouse_2020a,Pokhilko_2021a,Pokhilko_2021b,Pokhilko_2022} (or second Born) or $T$-matrix \cite{Liebsch_1981,Bickers_1989a,Bickers_1991,Katsnelson_1999,Katsnelson_2002,Zhukov_2005,vonFriesen_2010,Romaniello_2012,Gukelberger_2015,Muller_2019,Friedrich_2019,Biswas_2021,Zhang_2017,Li_2021b,Loos_2022} (or Bethe-Goldstone) approximations.
The central equation of many-body perturbation theory based on Hedin's equations is the so-called dynamical and non-hermitian quasiparticle equation which, within the $GW$ approximation, reads
Throughout the manuscript, the indices $p,q,r,s$ are general orbitals while $i,j,k,l$ and $a,b,c,d$ refers to occupied and virtual orbitals, respectively.
The self-energy can be physically understood as a correction to the Hartree-Fock (HF) problem (represented by $\bF$) accounting for dynamical screening effects.
Similarly to the HF case, Eq.~\eqref{eq:quasipart_eq} has to be solved self-consistently but the dynamical and non-hermitian nature of $\bSig^{\GW}(\omega)$, as well as its functional form, makes it much more challenging to solve from a practical point of view.
The matrix elements of $\bSig^{\GW}(\omega)$ have the following closed-form expression \cite{Hedin_1999,Tiago_2006,Bruneval_2012,vanSetten_2013,Bruneval_2016}
The diagonal matrix $\boldsymbol{\Omega}$ contains the positive eigenvalues of the RPA problen defined in Eq.~\eqref{eq:full_dRPA} and its elements $\Omega_\nu$ appear in Eq.~\eqref{eq:GW_selfenergy}.
In the Tamm-Dancoff approximation (TDA), one sets $\bB=\bO$ in Eq.~\eqref{eq:full_dRPA} which reduces to a Hermitian eigenvalue problem of the form $\bA\bX=\bX\bOm$ (hence $\bY=0$).
As mentioned above, because of the frequency dependence of the self-energy, solving exactly the quasiparticle equation \eqref{eq:quasipart_eq} is a rather complicated task.
The most popular strategy is the one-shot (perturbative) $G_0W_0$ scheme, where the self-consistency is completely abandoned, and the off-diagonal elements of Eq.~\eqref{eq:quasipart_eq} are neglected.
The previous equations are non-linear with respect to $\omega$ and therefore have multiple solutions $\epsilon_{p,s}$ for a given $p$ (where the index $s$ is numbering solutions).
One obvious drawback of the one-shot scheme mentioned above is its starting point dependence.
Indeed, in Eq.~\eqref{eq:G0W0} we choose to rely on HF orbital energies but this is arbitrary and one could have chosen Kohn-Sham energies (and orbitals) instead.
As commonly done, one can even ``tune'' the starting point to obtain the best possible one-shot $GW$ quasiparticle energies. \cite{Korzdorfer_2012,Marom_2012,Bruneval_2013,Gallandi_2015,Caruso_2016,Gallandi_2016}
Alternatively, one may solve iteratively the set of quasiparticle equations \eqref{eq:G0W0} to reach convergence of the quasiparticle energies, leading to the partially self-consistent scheme named ev$GW$.
However, if one of the quasiparticle equations does not have a well-defined quasiparticle solution, reaching self-consistency can be challenging, if not impossible.
Even at convergence, the starting point dependence is not totally removed as the quasiparticle energies still depend on the initial set of orbitals. \cite{Marom_2012}
In order to update both the orbitals and their corresponding energies, one must consider the off-diagonal elements in $\bSig^{\GW}(\omega)$.
To avoid solving the non-Hermitian and dynamic quasiparticle equation defined in Eq.~\eqref{eq:quasipart_eq}, one can resort to the qs$GW$ scheme in which $\bSig^{\GW}(\omega)$ is replaced by a static approximation $\bSig^{\qsGW}$.
Then, the qs$GW$ equations are solved via a standard self-consistent field procedure similar to the HF algorithm where $\bF$ is replaced by $\bF+\bSig^{\qsGW}$.
Various choices for $\bSig^{\qsGW}$ are possible but the most popular is the following Hermitian approximation
which was first 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 functions by Ismail-Beigi. \cite{Ismail-Beigi_2017}
One of the main results of the present manuscript is the derivation, from first principles, of an alternative static Hermitian form for the $GW$ self-energy.
For example, the $\ii\eta$ term in the denominators of Eq.~\eqref{eq:GW_selfenergy}, sometimes referred to as a broadening parameter linked to the width of the quasiparticle peak, is similar to the usual imaginary-shift regularizer employed in various other theories plagued by the intruder-state problem. \cite{Surjan_1996,Forsberg_1997,Monino_2022,Battaglia_2022}.
However, this $\eta$ parameter stems from a regularization of the convolution to obtain $\Sigma$ and should theoretically be set to zero. \cite{Martin_2016}
Several other regularizers are possible \cite{Stuck_2013,Rostam_2017,Lee_2018a,Evangelista_2014b,Shee_2021} and in particular, it was shown in Ref.~\onlinecite{Monino_2022} that a regularizer inspired by the SRG had some advantages over the imaginary shift.
Nonetheless, 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.
where the flow parameter $s$ controls the extent of the decoupling and is related to an energy cutoff $\Lambda=s^{-1/2}$ that avoids states with energy denominators smaller than $\Lambda$ to be decoupled from the reference space, hence avoiding potential intruders.
To solve the flow equation at a lower cost than the one associated with the diagonalization of the initial Hamiltonian, one must introduce an approximate form for $\boldsymbol{\eta}(s)$.
This implies that the matrix elements of the off-diagonal part decrease in a monotonic way throughout the transformation.
Moreover, the coupling coefficients associated with the highest-energy determinants are removed first as we shall evidence in the perturbative analysis below.
Then, as performed in Sec.~\ref{sec:srggw}, one can collect order by order the terms in Eq.~\eqref{eq:flowEquation} and solve analytically the low-order differential equations.
By applying the SRG to $GW$, our aim is to gradually remove the coupling between the quasiparticle and the satellites resulting in a renormalized quasiparticle equation.
However, to do so, one must identify the coupling terms in Eq.~\eqref{eq:quasipart_eq}, which is not straightforward.
and the corresponding coupling blocks read [see Eq.~(\ref{eq:GW_sERI})]
\begin{align}
W^\text{2h1p}_{p,i\nu}& = W_{p,i\nu},
&
W^\text{2p1h}_{p,a\nu}& = W_{p,a\nu}.
\end{align}
The usual $GW$ non-linear equation can be obtained by applying L\"owdin partitioning technique \cite{Lowdin_1963} to Eq.~\eqref{eq:GWlin} yielding \cite{Bintrim_2021}
We refer to Ref.~\onlinecite{Bintrim_2021} for a detailed discussion of the up/downfolding processes of the $GW$ equations (see also Refs.~\onlinecite{Tolle_2022,Scott_2023}).
As can be readily seen in Eq.~\eqref{eq:GWlin}, the blocks $\bW^\text{2h1p}$ and $\bW^\text{2p1h}$ are coupling the 1h and 1p configuration to the 2h1p and 2p1h configurations.
Therefore, it is natural to define, within the SRG formalism, the diagonal and off-diagonal parts of the $GW$ effective Hamiltonian as
Once the closed-form expressions of the low-order perturbative expansions are known, they can be inserted in Eq.~\eqref{eq:downfolded_sigma} to define a renormalized version of the quasiparticle equation.
The choice of Wegner's generator in the flow equation [see Eq.~\eqref{eq:flowEquation}] implies that the off-diagonal correction is of order $\order*{\lambda}$ while the correction to the diagonal block is at least $\order*{\lambda^2}$.
Therefore, the zeroth-order Hamiltonian is independent of $s$ and we have
and, thanks to the diagonal structure of $\bF^{(0)}$ (which is a consequence of the HF starting point) and $\bC^{(0)}$, the differential equation for the coupling block in Eq.~\eqref{eq:W1} is easily solved and yields
At $s=0$ the elements $W_{p,q\nu}^{(1)}(0)$ are equal to the screened two-electron integrals defined in Eq.~\eqref{eq:GW_sERI}, while for $s\to\infty$, they tend to zero.
Therefore, $W_{p,q\nu}^{(1)}(s)$ are genuine renormalized two-electron screened integrals.
It is worth noting the close similarity of the first-order elements with the ones derived by Evangelista in Ref.~\onlinecite{Evangelista_2014b} in the context of single- and multi-reference perturbation theory (see also Ref.~\onlinecite{Hergert_2016}).
As can be readily seen above, $\bF^{(2)}$ is the only second-order block of the effective Hamiltonian contributing to the second-order SRG quasiparticle equation.
Schematic evolution of the quasiparticle equation as a function of the flow parameter $s$ in the case of the dynamic SRG-$GW$ flow (magenta) and the static SRG-qs$GW$ flow (cyan). \ANT{Maybe we should replace dynamic by full?}
Therefore, the SRG flow continuously transforms the dynamical self-energy $\widetilde{\bSig}^{\GW}(\omega; s)$ into a static correction $\widetilde{\bF}^{(2)}(s)$.
As illustrated in Fig.~\ref{fig:flow}, this transformation is done gradually starting from the states that have the largest denominators in Eq.~\eqref{eq:static_F2}.
Because the large-$s$ limit of Eq.~\eqref{eq:GW_renorm} is purely static and hermitian, the new alternative form of the self-energy reported in Eq.~\eqref{eq:static_F2} can be naturally used in qs$GW$ calculations to replace Eq.~\eqref{eq:sym_qsgw}.
Unfortunately, as we shall discuss further in Sec.~\ref{sec:results}, as $s\to\infty$, self-consistency is once again quite difficult to achieve, if not impossible.
However, one can define a more flexible new static self-energy, which will be referred to as SRG-qs$GW$ in the following, by discarding the dynamic part in Eq.~\eqref{eq:GW_renorm}.
Note that the SRG static approximation is naturally Hermitian as opposed to the usual case [see Eq.~(\ref{eq:sym_qsGW})] where it is enforced by brute-force symmetrization.
Another important difference is that the SRG regularizer is energy dependent while the imaginary shift is the same for every denominator of the self-energy.
It is well-known that in traditional qs$GW$ calculations, increasing the imaginary shift $\ii\eta$ to ensure convergence in difficult cases is most often unavoidable.
Similarly, in SRG-qs$GW$, one might need to decrease the value of $s$ to ensure convergence.
The fact that the $s\to\infty$ static limit does not always 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}$.
It is instructive to plot the functional form of both regularizing functions, this is done in Fig.~\ref{fig:plot}.
The surfaces are plotted for a regularizing parameter value of $\eta=1$, where we have set $s=1/2\eta^2$ such that the first order Taylor expansion around $(0,0)$ of both functional form are equal.
One can observe that the SRG surface is much smoother than its qs counterpart.
This is due to the fact that the SRG functional has less irregularities for $\eta=0$, in fact there is a single singularity at $x=y=0$.
To conclude this section, the case of discontinuities will be briefly discussed.
Indeed, it has been previously mentioned that intruder states are responsible for both the poor convergence of qs$GW$ and discontinuities in physical quantities at the $\GOWO$ level.
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.
will reverse the situation and now a finite value of $t$ will be well-designed to avoid discontinuities in the renormalized dynamic part.
The dynamic part after the change of variable is actually closely related to the SRG-inspired regularizer introduced by Monino and Loos in Ref.~\onlinecite{Monino_2022}.
\titou{Our set is composed by XX closed-shell organic molecules, displayed in Fig.~??, with singlet ground states.}
Following the same philosophy as the \textsc{quest} database for neutral excited states, \cite{Loos_2020d,Veril_2021} their geometries have been optimized at the CC3 level \cite{Christiansen_1995b,Koch_1997} in the aug-cc-pVTZ basis set using the \textsc{cfour} program. \cite{CFOUR}
The reference CCSD(T) principal ionization potentials (IPs) and electron affinities (EAs) have been obtained using Gaussian 16 \cite{g16} with default parameters, that is, within the restricted and unrestriced HF formalism for the neutral and charged species, respectively.
The two qs$GW$ variants considered in this work have been implemented in an in-house program, named \textsc{quack}. \cite{QuAcK}
The $GW$ implementation closely follows the one of \textsc{molgw}. \cite{Bruneval_2016}
In all $GW$ calculations, we use the aug-cc-pVTZ cartesian basis set and self-consistency is performed on all (occupied and virtual) orbitals, including core orbitals.
The maximum size of the DIIS space \cite{Pulay_1980,Pulay_1982} and the maximum number of iterations were set to 5 and 64, respectively.
In practice, one may achieve convergence, in some cases, by adjusting these parameters or by using an alternative mixing scheme.
However, in order to perform a black-box comparison, these parameters have been fixed to these default values.
Principal IP of the water molecule in the aug-cc-pVTZ basis set as a function of the flow parameter $s$ for the SRG-qs$GW$ method with and without TDA.
Principal IP of the \ce{Li2}, \ce{LiH} and \ce{BeO} in the aug-cc-pVTZ basis set as a function of the flow parameter $s$ for the SRG-qs$GW$ method with and without TDA.
Reference values (HF, qs$GW$ with and without TDA) are also reported as dashed lines.
Figure \ref{fig:fig2} shows the error of various methods for the principal IP with respect to the CCSD(T) reference value.
The HF IP (dashed black line) is overestimated; this is a consequence of the missing correlation and the lack of orbital relaxation for the cation, a result that is well understood. \cite{SzaboBook,Lewis_2019}
At $s=0$, the IP is equal to its HF counterpart as expected from the discussion of Sec.~\ref{sec:srggw}.
For $s\to\infty$, the IP reaches a plateau at an error that is significantly smaller than their $s=0$ starting point.
Even more, the value associated with this plateau is slightly more accurate than its qs$GW$ counterpart.
However, the SRG-qs$GW$ error do not decrease smoothly between the initial HF value and the $s\to\infty$ limit as for small $s$ values it is actually worst than the HF starting point.
where $k$ is the index of the highest molecular MO (HOMO).
The first term is the zeroth order IP and the two following terms come from the 2h1p and 2p1h coupling, respectively.
The denominator of the last term is always positive while the 2h1p term is negative.
When $s$ is increased, the first states that will be decoupled from the HOMO will be the 2p1h ones because their energy difference with the HOMO is larger than the ones of the 2h1p block.
Therefore, for small $s$ only the last term of Eq.~\eqref{eq:2nd_order_IP} will be partially included resulting in a positive correction to the IP.
As soon as $s$ is large enough to decouple the 2h1p block as well the IP will start to decrease and eventually go below the $s=0$ initial value as observed in Fig.~\ref{fig:fig2}.
The left panel of Fig.~\ref{fig:fig3} shows the results for the Lithium dimer, \ce{Li2} is an interesting case because the HP IP is actually underestimated.
Finally, the beryllium oxide is considered a prototypical example of a molecular system difficult to converge because of intruder states. \cite{vanSetten_2015,Veril_2018,Forster_2021}
However, as we will see in the next subsection these are just particular molecular systems and on average the RPA polarizability performs better than the TDA one.
\caption{First ionization potential (left) and first electron attachment (right) in eV calculated using $\Delta$CCSD(T) (reference), HF, $G_0W_0$@HF, qs$GW$ and SRG-qs$GW$. The statistical descriptors are computed for the errors with respect to the reference. \ANT{Maybe change the values of SRG with the one for s=1000}}
The test set considered in this study is based on the GW20 set introduced by Lewis and Berkelbach, \cite{Lewis_2019} which is made of the 20 smallest atoms and molecules of the GW100 benchmark set. \cite{vanSetten_2015}
This set has been augmented with the MgO and O3 molecules (which are part of GW100 as well) because they are known to possess intruder states. \cite{vanSetten_2015,Forster_2021}
As mentioned previously the HF IPs are overestimated with a mean signed error (MSE) of \SI{0.64}{\electronvolt} and a mean absolute error (MAE) of \SI{0.74}{\electronvolt}.
Performing a one-shot $G_0W_0$ calculation on top these mean-field results allows to divided by more than two the MSE and MAE, \SI{0.26}{\electronvolt} and \SI{0.32}{\electronvolt}, respectively.
However, there are still outliers with quite large errors, for example the IP of the dinitrogen is overestimated by \SI{1.56}{\electronvolt}.
Self-consistency can mitigate the error of the outliers as the maximum error at the qs$GW$ level is now \SI{0.56}{\electronvolt} and the standard deviation error (SDE) is decreased from \SI{0.39}{\electronvolt} for $G_0W_0$ to \SI{0.18}{\electronvolt} for qs$GW$.
In addition, the MSE and MAE (\SI{0.24}{\electronvolt}/\SI{0.25}{\electronvolt}) are also slightly improved with respect to $G_0W_0$@HF.
Now turning to the new results of this manuscript, \ie the alternative self-consistent scheme SRG-qs$GW$.
Table~\ref{tab:tab1} shows the SRG-qs$GW$ values for $s=100$.
The statistical descriptors corresponding to the alternative static self-energy are all improved with respect to qs$GW$.
Of course these are slight improvements but this is done with no additional computational cost and can be implemented really quickly just by changing the form of the static approximation.
The evolution of the statistical descriptors with respect to the various methods considered in Table~\ref{tab:tab1} is graphically illustrated by Fig.~\ref{fig:fig4}.
The decrease of the MSE and SDE correspond to a shift of the maximum toward zero and a shrink of the distribution width, respectively.
Indeed, up to $s=10^3$ self-consistency can be attained without any problems (mean and max number of iterations = n for s=100).
For $s=10^4$, convergence could not be attained for 12 molecules out of 22, meaning that some intruder states were included in the static correction for this value of $s$.
This is illustrated in the case of \ce{H2O} in the upper panel Fig.
However, this is not a problem as the MAE is already well converged before the intruder states are added to the SRG-qs$GW$ static self-energy (lower panel).
On the other hand, the qs$GW$ convergence with respect to $\eta$ is more difficult to evaluate.
The whole set considered in this work could be converged for $\eta=0.1$.
However, as soon as we decrease $\eta$ self-consistency could not be attained for the whole set of molecules using the black-box convergence parameters (see Sec.~\ref{sec:comp_det}).
Unfortunately, the convergence of the IP is not as tight as in the SRG case.
The values of the IP that could be converged for $\eta=0.01$ can vary between $10^{-3}$ and $10^{-1}$ with respect to $\eta=0.1$.
% \begin{table}
% \caption{First ionization potential in eV calculated using $G_0W_0^{\text{TDA}}$@HF, qs$GW^{\text{TDA}}$ and SRG-qs$GW^{\text{TDA}}$. The statistical descriptors are computed for the errors with respect to the reference.}
\caption{First electron attachment in eV calculated using $\Delta$CCSD(T) (reference), HF, $G_0W_0$@HF, qs$GW$ and SRG-qs$GW$. The statistical descriptors are computed for the errors with respect to the reference.}
Finally, we compare the performance of HF, $G_0W_0$@HF, qs$GW$ and SRG-qs$GW$ again but for the principal electron attachement (EA) energies.
The raw results are given in Tab.~\ref{tab:tab2} while the corresponding histograms of the error distribution are plotted in Fig.~\ref{fig:fig6}.
The HF EA are understimated in averaged with some large outliers while $G_0W_0$@HF mitigates the average error there are still large outliers.
The performance of the two qs$GW$ schemes are quite similar for EA, \ie a MAE of \SI{\sim 0.1}{\electronvolt} and the error of the outliers is reduced with respect to $G_0W_0$@HF.
\ANT{Maybe we should mention that some EA are not chemically meaningful.}
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{Non-TDA $GW$ and $GW_{\text{TDHF}}$ equations}
% \label{sec:nonTDA}
% %%%%%%%%%%%%%%%%%%%%%%
% The matrix elements of the $GW$ self-energy within the TDA are the same as in Eq.~\eqref{eq:GW_selfenergy} but the screened integrals are now defined as
% where $\bX$ is the eigenvector matrix of the TDA particle-hole dRPA problem obtained by setting $\bB= \bO$ in Eq.~\eqref{eq:full_dRPA}, \ie
% \begin{equation}
% \label{eq:TDA_dRPA}
% \bA \bX = \bX \boldsymbol{\Omega}.
% \end{equation}
% Defining an unfold version of this equation that 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}
% which already depends on the screened integrals and therefore require the knowledge of the eigenvectors of the dRPA problem defined in Eq.~\eqref{eq:full_dRPA}.
% Within the TDA the renormalized matrix elements have the same $s$ dependence as in the RPA case but the $s=0$ screened integrals $W_{p,q\nu}$ and eigenvalues $\Omega_\nu$ are replaced by the ones of Eq.~\eqref{eq:GWnonTDA_sERI} and \eqref{eq:TDA_dRPA}, respectively.