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 naturally static and Hermitian form of the 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, slightly improves the overall accuracy, and is straightforward to implement in existing code.
The one-body Green's function provides 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 electronic 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$ involves 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 \cite{Golze_2019} for an extensive list of current challenges in Green's function methods.
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 eigenvalue-only self-consistent $GW$\cite{Shishkin_2007a,Blase_2011b,Marom_2012,Kaplan_2016,Wilhelm_2016} (ev$GW$) and quasiparticle self-consistent $GW$\cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016} (qs$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 and/or irregularities 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, in the quasiparticle equation, of a regularizer inspired by the similarity renormalization group (SRG). \cite{Monino_2022}
Encouraged by this study and 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,Coveney_2023} the present work investigates the application of the SRG formalism in $GW$-based methods.
The SRG formalism has been developed independently by Wegner \cite{Wegner_1994} in the context of condensed matter systems and Glazek \& Wilson \cite{Glazek_1993,Glazek_1994} in light-front quantum field theory.
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 these because the decoupling of each external configuration is inversely proportional to its energy difference with the reference space.
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 second-order Green's function \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 the $T$-matrix approximation.\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}
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$ refer 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(\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(\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 problem 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) $GW$ scheme, $G_0W_0$, 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,z}$ for a given $p$ (where the index $z$ is numbering solutions).
The solution with the largest weight $Z_p \equiv Z_{p,z=0}$ is referred to as the quasiparticle while the others are known as satellites (or shake-up transitions).
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}
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(\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}$.
which was first introduced by Faleev and co-workers \cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Lei_2022} before being derived by Ismail-Beigi as the effective Hamiltonian that minimizes the length of the gradient of the Klein functional for non-interacting Green's functions. \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 qs$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 is required to define the Fourier transformation between time and energy representation 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,Coveney_2023} 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}$.
For a given value of $s$, only states with energy difference (with respect to the reference space) greater than $\Lambda$ are decoupled from the reference space, hence avoiding potential intruders.
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.
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}
Equations \eqref{eq:GWlin} and \eqref{eq:quasipart_eq} yield exactly the same quasiparticle and satellite energies but one is linear and the other is not.
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_2023,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}$. \cite{Hergert_2016}
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
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).
Therefore, the SRG flow continuously transforms the dynamical self-energy $\widetilde{\bSig}(\omega; s)$ into a static correction $\widetilde{\bF}^{(2)}(s)$.
As illustrated in Fig.~\ref{fig:flow} (magenta curve), this transformation is done gradually starting from the states that have the largest denominators in Eq.~\eqref{eq:static_F2}.
For a fixed value of the energy cutoff $\Lambda$, if $\abs*{\Delta_{pr}^{\nu}}\gg\Lambda$, then $W_{pr}^{\nu} e^{-(\Delta_{pr}^{\nu})^2 s}\approx0$, meaning that the state is decoupled from the 1h and 1p configurations, while, for $\abs*{\Delta_{pr}^{\nu}}\ll\Lambda$, we have $W_{pr}^{\nu}(s)\approx W_{pr}^{\nu}$, that is, the state remains coupled.
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} (see cyan curve in Fig.~\ref{fig:flow}).
Note that the static SRG-qs$GW$ approximation defined in Eq.~\eqref{eq:SRG_qsGW} is straightforward to implement in existing code and is naturally Hermitian as opposed to the usual case [see Eq.~\eqref{eq:sym_qsGW}] where it is enforced by brute-force symmetrization.
Indeed, the fact that SRG-qs$GW$ calculations do not always converge in the large-$s$ limit is expected as, in this limit, potential intruder states have been included.
Therefore, one should use a value of $s$ large enough to include as many states as possible but small enough to avoid intruder states.
These have been 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 $(x,y)=(0,0)$ of both functional forms is equal.
The convergence properties and the accuracy of both static approximations are quantitatively gauged in Sec.~\ref{sec:results}.
To conclude this section, we briefly discussed the case of discontinuities mentioned in Sec.~\ref{sec:intro}.
Indeed, it has been previously mentioned that intruder states are responsible for both the poor convergence of qs$GW$ and discontinuities in physical quantities. \cite{Loos_2018b,Veril_2018,Loos_2020e,Berger_2021,DiSabatino_2021,Monino_2022,Scott_2023}
Is it then possible to rely on the SRG machinery to remove discontinuities?
However, as we have seen just above the functional form of the renormalized equation makes it possible to choose $s$ such that there is no intruder states in its static part.
Performing a bijective transformation of the form,
on the renormalized quasiparticle equation \eqref{eq:GW_renorm} reverses the situation and makes it possible to choose $t$ such that there is no intruder states in the dynamic part, hence removing discontinuities.
\textcolor{red}{The intruder-state free dynamic part makes it possible to define a SRG-$G_0W_0$ and SRG-ev$GW$.
The main manuscript focus on SRG-qs$GW$ but the performance of SRG-$G_0W_0$ and SRG-ev$GW$ are discussed in the {\SupInf} for the sake of completeness.}
Note that, after this transformation, the form of the regularizer is actually closely related to the SRG-inspired regularizer introduced by Monino and Loos in Ref.~\onlinecite{Monino_2022}.
Our set of systems is composed by closed-shell compounds that correspond to the 50 smallest atoms and molecules (in terms of the number of electrons) of the $GW$100 benchmark set. \cite{vanSetten_2015}
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/aug-cc-pVTZ basis level \cite{Christiansen_1995b,Koch_1997} using the \textsc{cfour} program. \cite{CFOUR}
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 $\eta$ value has been set to \num{e-3} for the conventional $G_0W_0$ calculations (where we eschew linearizing the quasiparticle equation) while, for the qs$GW$ calculations, $\eta$ has been chosen as the largest value where one successfully converges the 50 systems composing the test set.
The $\Delta$CCSD(T) principal ionization potentials (IPs) and electron affinities (EAs) have been obtained using \textsc{gaussian 16}\cite{g16} (with default parameters) within the restricted and unrestricted formalism for the neutral and charged species, respectively.
%The first step will be to analyze in depth the behavior of the two static self-energy approximations in a few test cases.
%Then, the accuracy of the principal IPs and EAs produced by the qs$GW$ and SRG-qs$GW$ schemes are statistically gauged over the test set of molecules described in Sec.~\ref{sec:comp_det}.
Error [with respect to $\Delta$CCSD(T)] in the principal IP of water in the aug-cc-pVTZ basis set as a function of the flow parameter $s$ for SRG-qs$GW$ (green curve).
Error [with respect to $\Delta$CCSD(T)] in the principal IP of \ce{Li2}, \ce{LiH}, and the principal EA of \ce{F2} in the aug-cc-pVTZ basis set as a function of the flow parameter $s$ for the SRG-qs$GW$ method (green curves).
Figure \ref{fig:fig3} shows the error in the principal IP [with respect to the $\Delta$CCSD(T) reference value] as a function of the flow parameter in SRG-qs$GW$ (green curve).
The IP at the HF level (cyan line) is too large; this is a consequence of the missing correlation and the lack of orbital relaxation in the cation, a result that is well understood. \cite{SzaboBook,Lewis_2019}
The usual qs$GW$ scheme (blue line) brings a quantitative improvement as the IP is now within \SI{0.3}{\eV} of the reference value.
The first term of the right-hand side of Eq.~\eqref{eq:2nd_order_IP} is the zeroth-order IP and the following two terms originate from the 2h1p and 2p1h coupling, respectively.
As $s$ increases, the first states that decouple from the HOMO are the 2p1h configurations because their energy difference with respect to the HOMO is larger than the ones associated with the 2h1p block.
Therefore, for small $s$, only the last term of Eq.~\eqref{eq:2nd_order_IP} is partially included, resulting in a positive correction to the IP.
As soon as $s$ is large enough to decouple the 2h1p block, the IP starts decreasing and eventually goes below the initial value at $s=0$, as observed in Fig.~\ref{fig:fig3}.
%\ANT{I don't know if we should remove this paragraph and the TDA curves in Fig 3 and 4 or not...}
%In addition, the qs$GW$ and SRG-qs$GW$ methods based on a TDA screening (dubbed qs$GW^\TDA$ and SRG-qs$GW^\TDA$) are also considered in Fig.~\ref{fig:fig3}.
%The TDA values are now underestimating the IP, unlike their RPA counterparts.
%For both static self-energies, the TDA leads to a slight increase in the absolute error.
Next, the flow parameter dependence of SRG-qs$GW$ is investigated for the principal IP of two additional molecular systems as well as the principal EA of \ce{F2}.
The left panel of Fig.~\ref{fig:fig4} shows the results for the lithium dimer, \ce{Li2}, which is an interesting case because, unlike in water, HF underestimates the reference IP.
Now turning to lithium hydride, \ce{LiH} (see middle panel of Fig.~\ref{fig:fig4}), we see that the qs$GW$ IP is actually worse than the fairly accurate HF value.
Finally, we also consider the evolution with respect to $s$ of the principal EA of \ce{F2} that is displayed in the right panel of Fig.~\ref{fig:fig4}.
% Finally, beryllium oxide, \ce{BeO}, is considered as it is a well-known example where it is particularly difficult to converge self-consistent $GW$ calculations due to the presence of intruder states. \cite{vanSetten_2015,Veril_2018,Forster_2021}
% The SRG-qs$GW$ calculations could be converged without any issue even for large $s$ values.
% The convergence properties of both schemes will be systematically investigated in the next subsection.
% Once again, a plateau is attained and the corresponding value is slightly more accurate than its qs$GW$ counterpart.
%Note that, for \ce{LiH} and \ce{BeO}, the TDA actually improves the accuracy compared to the RPA-based qs$GW$ schemes.
%However, as we shall see in Sec.~\ref{sec:SRG_vs_Sym}, these are special cases, and, on average, the RPA polarizability performs better than its TDA version.
%Also, SRG-qs$GW^\TDA$ is better than qs$GW^\TDA$ in the three cases of Fig.~\ref{fig:fig4}.
%However, it is not a general rule.
%Therefore, it seems that the effect of the TDA cannot be systematically predicted.
Histogram of the errors [with respect to $\Delta$CCSD(T)] for the principal IP of the $GW$50 test set calculated using HF, $G_0W_0$@HF, qs$GW$, and SRG-qs$GW$.
All calculations are performed with the aug-cc-pVTZ basis.
As previously mentioned, the HF approximation overestimates the IPs with a mean signed error (MSE) of \SI{0.56}{\eV} and a mean absolute error (MAE) of \SI{0.69}{\eV}.
Performing a $G_0W_0$ calculation on top of this mean-field starting point, $G_0W_0$@HF, reduces by more than a factor two the MSE and MAE, \SI{0.29}{\eV} and \SI{0.33}{\eV}, respectively.
Self-consistency mitigates the error of the outliers as the MAE at the qs$GW$ level is now \SI{0.57}{\eV} and the standard deviation of the error (SDE) is decreased from \SI{0.31}{\eV} for $G_0W_0$@HF to \SI{0.18}{\eV} for qs$GW$.
Of course, these are small improvements but this is done with no additional computational cost and it can be easily implemented in existing code by changing the form of the static self-energy.
The evolution of the statistical descriptors with respect to the various methods considered in Table \ref{tab:tab1} is graphically illustrated in Fig.~\ref{fig:fig4}.
The decrease of the MSE and SDE correspond to a shift of the maximum of the distribution toward zero and a contraction of the distribution width, respectively.
On the other hand, the qs$GW$ convergence behavior is more erratic as shown by the blue curve of Fig.~\ref{fig:fig6} where we report the variation of the qs$GW$ MAE as a function of $\eta=\sqrt{1/(2s)}$.
At $\eta=\num{e-2}$ ($s=\num{5e3}$), convergence could not be reached for 13 molecules while 2 systems were already problematic at $\eta=\num{5e-2}$ ($s=200$).
For example, out of the 37 molecules that could be converged for $\eta=\num{e-2}$, the variation of the IP with respect to $\eta=\num{5e-2}$ can go up to \SI{0.1}{\eV}.
Histogram of the errors [with respect to $\Delta$CCSD(T)] for the principal EA of the $GW$50 test set calculated using HF, $G_0W_0$@HF, qs$GW$, and SRG-qs$GW$.
Finally, we compare the performance of HF, $G_0W_0$@HF, qs$GW$, and SRG-qs$GW$ again but for the principal EAs of $GW$50.
The raw data are reported in Table \ref{tab:tab1} while the corresponding histograms of the error distribution are plotted in Fig.~\ref{fig:fig7}.
The HF EAs are, on average, underestimated with a MAE of \SI{0.31}{\eV} and some clear outliers: \SI{-2.03}{\eV} for \ce{F2} and \SI{1.04}{\eV} for \ce{CH2O}, for example.
$G_0W_0$@HF mitigates the average error (MAE equals to \SI{0.16}{\eV}) but the minimum and maximum error values are not satisfactory.
The performance of the two qs$GW$ schemes are quite similar for EAs with MAEs of the order of \SI{0.1}{\eV}.
These two partially self-consistent methods reduce also the minimum errors but, interestingly, they do not decrease the maximum error compared to HF.
Note that a positive EA indicates a bounded anion state, which can be accurately described by the methods considered in this study.
However, a negative EA suggests a resonance state, which is beyond the scope of the methods used in this study, including the $\Delta$CCSD(T) reference.
As such, it is not advisable to assign a physical interpretation to these values.
Nonetheless, it is possible to compare $GW$-based and $\Delta$CCSD(T) values in such cases, provided that the comparison is limited to a given basis set.
The present manuscript applies the similarity renormalization group (SRG) to the $GW$ approximation of many-body perturbation theory, which is known to be plagued by intruder states.
The problems caused by intruder states in many-body perturbation theory are numerous but here we focus on the convergence issues caused by them.
SRG's central equation is the flow equation, which is usually solved numerically but can be solved analytically for low perturbation order.
Applying this approach in the $GW$ context yields closed-form renormalized expressions for the Fock matrix elements and the screened two-electron integrals.
These renormalized quantities lead to a regularized $GW$ quasiparticle equation, referred to as SRG-$GW$, which is the main result of this work.
By isolating the static component of SRG-$GW$, we obtain an alternative Hermitian and intruder-state-free self-energy that can be used in the context of qs$GW$ calculations.
This new variant is called SRG-qs$GW$.
Additionally, we demonstrate how SRG-$GW$ can effectively resolve the discontinuity problems that arise in $GW$ due to intruder states.
This provides a first-principles justification for the SRG-inspired regularizer proposed in Ref.~\onlinecite{Monino_2022}.
We first study the flow parameter dependence of the SRG-qs$GW$ IPs for a few test cases.
The results show that the IPs gradually evolve from the HF starting point at $s=0$ to a plateau value for $s\to\infty$ that is much closer to the $\Delta$CCSD(T) reference than the HF initial value.
For small values of the flow parameter, the SRG-qs$GW$ IPs are actually worse than their starting point.
Therefore, it is advisable to use the largest possible value of $s$, similar to qs$GW$ calculations where one needs to use the smallest possible $\eta$ value.
Next, we gauge the accuracy of the SRG-qs$GW$ principal IP for a test set of 50 atoms and molecules (referred to as $GW$50).
The results show that, on average, SRG-qs$GW$ is slightly better than its qs$GW$ parent.
Despite the fact that the increase in accuracy is relatively modest, it comes with no additional computational cost and is straightforward to implement, as only the expression of the static self-energy needs to be modified.
Moreover, SRG-qs$GW$ calculations are much easier to converge than their traditional qs$GW$ counterparts thanks to the intruder-state-free nature of SRG-qs$GW$.
Finally, the principal EAs of the $GW$50 set are also investigated.
It is found that the performances of qs$GW$ and SRG-qs$GW$ are quite similar in this case.
However, it should be noted that most of the anions of the $GW$50 set are resonance states, and the associated physics cannot be accurately described by the methods considered in this study.
Therefore, test sets of molecules with bound anions, such as this one of organic electron-acceptor molecules, \cite{Richard_2016,Gallandi_2016,Knight_2016,Dolgounitcheva_2016} and their accompanying accurate reference values are greatly valuable to the many-body perturbation theory community.
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.
