By upfolding the frequency-dependent $GW$ quasiparticle equation, we explain the appearance of multiple solutions and unphysical discontinuities in various physical quantities computed within the $GW$ approximation.
By considering the $GW$ self-energy as an effective Hamiltonian, the appearance of these multiple solutions and discontinuities can be directly related to the intruder states problem.
A simple and efficient regularization procedure is proposed to avoid such issues.
The $GW$ approximation of many-body perturbation theory \cite{Hedin_1965,Martin_2016} allows to compute accurate charged excitation (\ie, ionization potentials and electron affinities) in solids and molecules. \cite{Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019}
Its popularity in the molecular structure community is rapidly increasing \cite{Ke_2011,Bruneval_2012,Bruneval_2013,Bruneval_2015,Blase_2016,Bruneval_2016, Bruneval_2016a,Koval_2014,Hung_2016,Blase_2018,Boulanger_2014,Li_2017,Hung_2016,Hung_2017,vanSetten_2015,vanSetten_2018,vanSetten_2015, Maggio_2017,vanSetten_2018,Richard_2016,Gallandi_2016,Knight_2016,Dolgounitcheva_2016,Bruneval_2015,Krause_2015,Govoni_2018,Caruso_2016} thanks to its relatively low cost and somehow surprising accuracy for weakly-correlated systems. \cite{Bruneval_2021}
In recent studies, \cite{Loos_2018b,Veril_2018,Berger_2021,DiSabatino_2021} we discovered that one can observe (unphysical) irregularities and/or discontinuities in the energy surfaces of several key quantities (ionization potential, electron affinity, fundamental gap, total and correlation energies, as well as vertical excitation energies) even in the weakly correlated regime.
These issues were discovered in Ref.~\onlinecite{Loos_2018b} while studying a model two-electron system \cite{Seidl_2007,Loos_2009a,Loos_2009c} and they were further investigated in Ref.~\onlinecite{Veril_2018}, where we provided additional evidences and explanations of these undesirable features in real molecular systems.
In particular, we show that each branch of the self-energy $\Sigma$ is associated with a distinct quasiparticle solution, and that each switch between solutions implies a significant discontinuity in the quasiparticle energy due to the transfer of weight between two solutions of the quasiparticle equation. \cite{Veril_2018}
It was shown that these problems could be alleviated by using a static COHSEX self-energy \cite{Berger_2021} and by considering a fully self-consistent scheme. \cite{DiSabatino_2021}
However, none of these solutions is completely satisfying as a static approximation of the self-energy can induce significant loss of accuracy and fully self-consistent calculations can be quite challenging terms of implementation and cost.
In the present article, we provide further physical insights into the origin of these discontinuities by highlighting, in particular, the role of intruder states.
These new insights allow us to provide a cheap and efficient regularization scheme in order to avoid these issues.
We consider {\GOWO}\cite{Strinati_1980,Hybertsen_1985a,Hybertsen_1986,Godby_1988,Linden_1988,Northrup_1991,Blase_1994,Rohlfing_1995,Shishkin_2007} for the sake of simplicity but the same analysis can be performed in the case of (partially) self-consistent schemes.\cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016,Gui_2018,Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016}
Here, for the sake of simplicity, we consider a restricted Hartree-Fock (HF) starting point but the present analysis can be straightforwardly extended to a Kohn-Sham (KS) starting point.
Within the {\GOWO} approximation, in order to obtain the quasiparticle energies and the corresponding satellites, one solve, for each spatial orbital $p$, the following (non-linear) quasiparticle equation
Within the Tamm-Dancoff approximation, the screened two-electron integrals are given by
\begin{equation}
\ERI{pq}{m} = \sum_{ia}\ERI{pq}{ia} X_{ia,m}^\RPA
\end{equation}
where $\Om{m}{\RPA}$ and $\bX{m}{\RPA}$ are respectively the $m$th eigenvalue and eigenvector of the RPA problem in the Tamm-Dancoff approximation, \ie,
are two-electron integrals over the HF (spatial) orbitals $\MO{p}(\br)$.
As a non-linear equation, Eq.~\eqref{eq:qp_eq} has many solutions $\eps{p,s}{\GW}$ and their corresponding weight is given by the value of the so-called renormalization factor
In a well-behaved case, one of the solution (the so-called quasiparticle) $\eps{p}{\GW}\equiv\eps{p,s=0}{\GW}$ has a large weight $Z_{p}\equiv Z_{p,s=0}$.
The non-linear quasiparticle equation \eqref{eq:qp_eq} can be transformed into a larger linear problem via an upfolding process where the 2h1p and 2p1h sectors
The size of this eigenvalue problem is $N =1+ N^\text{2h1p}+ N^\text{2p1h}=1+ O^2 V + O V^2$, and this eigenvalue problem has to be solved for each orbital that one wants to correct.
Because the renormalization factor corresponds to the projection of the vector $\bc{}{(p,s)}$ onto the reference space, the weight of a solution $(p,s)$ is given by the the first coefficient of their corresponding eigenvector $\bc{}{(p,s)}$, \ie,
\begin{equation}
Z_{p,s} = \qty[ c_{1}^{(p,s)} ]^{2}
\end{equation}
It is important to understand that diagonalizing $\bH^{(p)}$ in Eq.~\eqref{eq:Hp} is completely equivalent to solving the quasiparticle equation \eqref{eq:qp_eq}.
The main difference between the two approaches is that, by diagonalizing Eq.~\eqref{eq:Hp}, one has directly access to the eigenvectors associated with each quasiparticle and satellites.
One can see this downfolding process as the construction of an frequency-dependent effective Hamiltonian where the reference (zeroth-order) space is composed by a single determinant of the 1h or 1p sector and the external (first-order) space by all the singly-excited states built from to this reference electronic configuration.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{An illustrative example}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Multiple solution issues in $GW$ appears all the time, especially for orbitals that are far in energy from the Fermi level.
Therefore, such issues are ubiquitous when one wants to compute core ionized states for example.
In order to illustrate the appearance and the origin of these multiple solutions, we consider the hydrogen molecule in the 6-31G basis set which corresponds to a system with 2 electrons and 4 spatial orbitals (one occupied and three virtual).
This example was already considered in our previous work \cite{Veril_2018} but here we provide further insights on the origin of the appearances of these multiple solutions.
Quasiparticle energies (left), correlation part of the self-energy (center) and renormalization factor (right) as functions of the internuclear distance $\RHH$ (in \si{\angstrom}) for various orbitals of \ce{H2} at the {\GOWO}@HF/6-31G level.
}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Figure \ref{fig:H2} shows the evolution of the quasiparticle energies at the {\GOWO}@HF/6-31G level as a function on the internuclear distance $\RHH$.
As one can see there are two problematic regions showing obvious discontinuities around $\RHH=\SI{1}{\AA}$ and $\RHH=\SI{1}{\AA}$
Note that these issues do not only appear at the {\GOWO} level but also at the partially self-consistent levels such as ev$GW$\cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016,Gui_2018} and qs$GW$. \cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016}
However, fully self-consistent $GW$ methods \cite{Stan_2006,Stan_2009,Rostgaard_2010,Caruso_2012,Caruso_2013,Caruso_2013a,Caruso_2013b,Koval_2014,Wilhelm_2018} where one considers not only the quasiparticle solutions but also the satellites at each iteration are free of these irregularities. \cite{DiSabatino_2021}
The $T$-matrix-based formalism as well as second-order Green's function (or second Born) scheme \cite{SzaboBook,Casida_1989,Casida_1991,Stefanucci_2013,Ortiz_2013, Phillips_2014,Phillips_2015,Rusakov_2014,Rusakov_2016,Hirata_2015,Hirata_2017} exhibit the same problems.
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).}
