\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts,siunitx}
\usepackage[version=4]{mhchem}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}

\usepackage[
	colorlinks=true,
    citecolor=blue,
    breaklinks=true
	]{hyperref}
\urlstyle{same}

\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{r\ed}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}

\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\QP}{\textsc{quantum package}}
\newcommand{\T}[1]{#1^{\intercal}}

% coordinates
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bx}{\boldsymbol{x}}
\newcommand{\dbr}{d\br}
\newcommand{\dbx}{d\bx}

% methods
\newcommand{\GW}{\text{$GW$}}	
\newcommand{\evGW}{ev$GW$}	
\newcommand{\qsGW}{qs$GW$}	
\newcommand{\GOWO}{$G_0W_0$}	
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\xc}{\text{xc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\co}{\text{c}}
\newcommand{\x}{\text{x}}
\newcommand{\KS}{\text{KS}}
\newcommand{\HF}{\text{HF}}
\newcommand{\RPA}{\text{RPA}}

% 
\newcommand{\Norb}{N}
\newcommand{\Nocc}{O}
\newcommand{\Nvir}{V}

% operators
\newcommand{\hH}{\Hat{H}}
\newcommand{\hS}{\Hat{S}}

% energies
\newcommand{\Enuc}{E^\text{nuc}}
\newcommand{\Ec}[1]{E_\text{c}^{#1}}
\newcommand{\EHF}{E^\text{HF}}

% orbital energies
\newcommand{\eps}[2]{\epsilon_{#1}^{#2}}
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}

% Matrix elements
\newcommand{\Sig}[2]{\Sigma_{#1}^{#2}}
\newcommand{\SigC}[1]{\Sigma^\text{c}_{#1}}
\newcommand{\rSigC}[1]{\widetilde{\Sigma}^\text{c}_{#1}}
\newcommand{\SigX}[1]{\Sigma^\text{x}_{#1}}
\newcommand{\SigXC}[1]{\Sigma^\text{xc}_{#1}}
\newcommand{\MO}[1]{\phi_{#1}}
\newcommand{\SO}[1]{\psi_{#1}}
\newcommand{\ERI}[2]{(#1|#2)}
\newcommand{\rbra}[1]{(#1|}
\newcommand{\rket}[1]{|#1)}


% Matrices
\newcommand{\bO}{\boldsymbol{0}}
\newcommand{\bI}{\boldsymbol{1}}
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\bvc}{\boldsymbol{v}}
\newcommand{\bSig}[1]{\boldsymbol{\Sigma}^{#1}}
\newcommand{\be}{\boldsymbol{\epsilon}}
\newcommand{\bOm}[1]{\boldsymbol{\Omega}^{#1}}
\newcommand{\bA}[2]{\boldsymbol{A}_{#1}^{#2}}
\newcommand{\bB}[2]{\boldsymbol{B}_{#1}^{#2}}
\newcommand{\bC}[2]{\boldsymbol{C}_{#1}^{#2}}
\newcommand{\bV}[2]{\boldsymbol{V}_{#1}^{#2}}
\newcommand{\bX}[2]{\boldsymbol{X}_{#1}^{#2}}
\newcommand{\bY}[2]{\boldsymbol{Y}_{#1}^{#2}}
\newcommand{\bZ}[2]{\boldsymbol{Z}_{#1}^{#2}}
\newcommand{\bc}[2]{\boldsymbol{c}_{#1}^{#2}}

% orbitals, gaps, etc
\newcommand{\IP}{I}
\newcommand{\EA}{A}
\newcommand{\HOMO}{\text{HOMO}}
\newcommand{\LUMO}{\text{LUMO}}
\newcommand{\Eg}{E_\text{g}}
\newcommand{\EgFun}{\Eg^\text{fund}}
\newcommand{\EgOpt}{\Eg^\text{opt}}
\newcommand{\EB}{E_B}

\newcommand{\RHH}{R_{\ce{H-H}}}


% addresses
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}

\begin{document}	

\title{Unphysical Discontinuities in $GW$ Methods and the Role of Intruder States}

\author{Enzo \surname{Monino}}
	\affiliation{\LCPQ}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
	\email{loos@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}

\begin{abstract}
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 state problem.
A simple and efficient regularization procedure is proposed to avoid such issues.
%\bigskip
%\begin{center}
%	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
%\end{center}
%\bigskip
\end{abstract}

\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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 growing \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}

The idea behind the $GW$ approximation (or Green's function-based methods in general) is to recast the many-body problem into a set of non-linear one-body equations. The introduction of the self-energy $\Sigma$ links the non-interacting Green's function $G_0$ to its fully-interacting version $G$ via the following Dyson equation:
\begin{equation}
	G = G_0 + G_0 \Sigma G
\end{equation}
Electron correlation is then explicitly incorporated into one-body quantities via a sequence of self-consistent steps known as Hedin's equations. \cite{Hedin_1965}
%which connect $G$, the irreducible vertex function $\Gamma$, the irreducible polarizability $P$, the dynamically-screened Coulomb interaction $W$, and $\Sigma$ through a set of five equations.
%\begin{subequations}
%\begin{align}
%	\label{eq:G}
%	& G(12) = G_0(12) + \int G_\text{H}(13) \Sigma(34) G(42) d(34),
%	\\
%	\label{eq:Gamma}
%	& \Gamma(123) = \delta(12) \delta(13) 
%	\notag 
%	\\
%	& \qquad \qquad		+ \int \fdv{\Sigma(12)}{G(45)} G(46) G(75) \Gamma(673) d(4567),
%	\\
%	\label{eq:P}
%	& P(12) = - i \int G(13) \Gamma(324) G(41) d(34),
%	\\
%	\label{eq:W}
%	& W(12) = v(12) + \int v(13) P(34) W(42) d(34),
%	\\
%	\label{eq:Sig}
%	& \Sigma(12) = i \int G(13) W(14) \Gamma(324) d(34),
%\end{align}
%\end{subequations}
%where $v$ is the bare Coulomb interaction, $\delta(12)$ is Dirac's delta function and $(1)$ is a composite coordinate gathering spin, space and time variables $(\sigma_1,\boldsymbol{r}_1,t_1)$.

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 showed 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 Coulomb-hole plus screened-exchange (COHSEX) \cite{Hedin_1965,Hybertsen_1986,Hedin_1999,Bruneval_2006} self-energy \cite{Berger_2021} or 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 in accuracy and fully self-consistent calculations can be quite challenging in terms of implementation and cost.

In the present article, via an upfolding process of the non-linear $GW$ equation, \cite{Bintrim_2021a} we provide further physical insights into the origin of these discontinuities by highlighting, in particular, the role of intruder states.
Inspired by regularized electronic structure theories, \cite{Lee_2018a,Evangelista_2014b} these new insights allow us to provide a cheap and efficient regularization scheme in order to avoid these issues.

Here, we consider the one-shot {\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}
Moreover, we consider a restricted Hartree-Fock (HF) starting point but it can be straightforwardly extended to a Kohn-Sham (KS) starting point.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Downfold: The non-linear $GW$ problem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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
\begin{equation}
\label{eq:qp_eq}
	\eps{p}{\HF} + \SigC{p}(\omega) - \omega = 0
\end{equation}
where $\eps{p}{\HF}$ is the $p$th HF orbital energy and the correlation part of the {\GOWO} self-energy reads
\begin{equation}
	\SigC{p}(\omega) 
	= \sum_{im} \frac{2\ERI{pi}{m}^2}{\omega - \eps{i}{\HF} + \Om{m}{\RPA}}
	+ \sum_{am} \frac{2\ERI{pa}{m}^2}{\omega - \eps{a}{\HF} - \Om{m}{\RPA}}
\end{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, 
\begin{equation}
	\bA{}{\RPA} \cdot \bX{m}{\RPA} = \Om{m}{\RPA} \bX{m}{\RPA}
\end{equation}
with
\begin{equation}
	A_{ia,jb}^{\RPA} = (\eps{a}{\HF} - \eps{i}{\HF}) \delta_{ij} \delta_{ab} + \ERI{ia}{bj}
\end{equation}
and 
\begin{equation}
	\ERI{pq}{rs} = \iint \MO{p}(\br_1) \MO{q}(\br_1) \frac{1}{\abs{\br_1 - \br_2}} \MO{r}(\br_2) \MO{s}(\br_2) d\br_1 \dbr_2
\end{equation}
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
\begin{equation}
	0 \le Z_{p,s} = \qty[ 1 - \eval{\pdv{\SigC{p}(\omega)}{\omega}}_{\omega = \eps{p,s}{\GW}} ]^{-1} \le 1
\end{equation} 
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}$.
Note that we have the following important conservation rules \cite{Martin_1959,Baym_1961,Baym_1962}
\begin{align}
	\sum_{s} Z_{p,s} & = 1
	&
	\sum_{s} Z_{p,s} \eps{p,s}{\GW} & = \eps{p}{\HF}
\end{align}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Upfolding: the linear $GW$ problem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
are upfolded from the 1h and 1p sectors. \cite{Bintrim_2021a}
For each orbital $p$, this yields a linear eigenvalue problem of the form
\begin{equation}
	\bH^{(p)} \bc{}{(p,s)} = \eps{p,s}{\GW} \bc{}{(p,s)}
\end{equation}
with
\begin{equation}
\label{eq:Hp}
	\bH^{(p)} = 
	\begin{pmatrix}
		\eps{p}{\HF}		&	\bV{p}{\text{2h1p}}	&	\bV{p}{\text{2p1h}}
		\\
		\T{(\bV{p}{\text{2h1p}})}	&	\bC{}{\text{2h1p}}			&	\bO
		\\
		\T{(\bV{p}{\text{2p1h}})}	&	\bO				&	\bC{}{\text{2p1h}}	
	\end{pmatrix}
\end{equation}
where
\begin{align}
	C^\text{2h1p}_{ija,kcl} & = \qty[ \qty( \eps{i}{\HF} + \eps{j}{\HF} - \eps{a}{\HF}) \delta_{jl} \delta_{ac} - 2 \ERI{ja}{cl} ] \delta_{ik} 
	\\
	C^\text{2p1h}_{iab,kcd} & = \qty[ \qty( \eps{a}{\HF} + \eps{b}{\HF} - \eps{i}{\HF}) \delta_{ik} \delta_{ac} + 2 \ERI{ai}{kc} ] \delta_{bd} 
\end{align}
and the corresponding coupling blocks read
\begin{align}
	V^\text{2h1p}_{p,klc} & = \sqrt{2} \ERI{pk}{cl}
	&
	V^\text{2p1h}_{p,kcd} & = \sqrt{2} \ERI{pd}{kc}
\end{align}
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.
Note, however, that the block $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$ do not need to be recomputed for each orbital.
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}.
This can be further illustrate by expanding the secular equation associated with Eq.~\eqref{eq:Hp}
\begin{equation}
	\det[ \bH^{(p)} - \omega \bI ] = 0 
\end{equation}
and comparing it with Eq.~\eqref{eq:qp_eq} by setting
\begin{multline}
	\SigC{p}(\omega) 
	= \bV{p}{\text{2h1p}} \cdot \qty(\omega \bI - \bC{}{\text{2h1p}} )^{-1} \cdot \T{\qty(\bV{p}{\text{2h1p}})} 
	\\
	+ \bV{p}{\text{2p1h}} \cdot \qty(\omega \bI - \bC{}{\text{2p1h}} )^{-1} \cdot \T{\qty(\bV{p}{\text{2p1h}})} 
\end{multline}
where $\bI$ is the identity matrix.
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 a 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 2h1p and 2p1h configurations.  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	FIGURE 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
%	\includegraphics[width=\linewidth]{fig1a}
%	\includegraphics[width=\linewidth]{fig1b}
	\caption{
	\label{fig:H2}
	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}$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introducing regularized $GW$ methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

One way to hamper such issues is to resort to regularization of the $GW$ self-energy.
Of course, the way of regularizing the self-energy is not unique but here we consider 3 different ways directly imported from MP2 theory.

This helps greatly convergence for (partially) self-consistent $GW$ methods.

The regularized self-energy is 
\begin{equation}
\begin{split}
	\rSigC{p}(\omega;s) 
	& = \sum_{im} 2\ERI{pi}{m}^2 f_s(\omega - \eps{i}{\HF} + \Om{m}{\RPA})
	\\
	& + \sum_{am} 2\ERI{pa}{m}^2 f_s(\omega - \eps{a}{\HF} - \Om{m}{\RPA})
\end{split}
\end{equation}
where various choices for the ``regularizer'' $f_s$ are possible.
Our investigation have shown that
\begin{equation}
	f_s(\Delta) = \frac{1-e^{-2s\Delta^2}}{\Delta}
\end{equation}
is a very convenient form and has been derived from the flow equation within driven-similarity renormalization group \cite{Evangelista_2014}
Of course, by construction, we have
\begin{equation}
	\lim_{s\to\infty} \rSigC{p}(\omega;s) = \SigC{p}(\omega) 
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Concluding remarks}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.

%%%%%%%%%%%%%%%%%%%%%%%%
\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).}
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{ufGW}
%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}