\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{\Ne}{N}
\newcommand{\Norb}{K}
\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}}}
\newcommand{\ii}{\mathrm{i}}

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

\begin{document}	

\title{Unphysical Discontinuities, Intruder States and Regularization in $GW$ Methods}

\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 recasting the non-linear frequency-dependent $GW$ quasiparticle equation into a linear eigenvalue problem, we explain the appearance of multiple solutions and unphysical discontinuities in various physical quantities computed within the $GW$ approximation.
Considering the $GW$ self-energy as an effective Hamiltonian, these issues 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, electron affinities and fundamental gaps) in solids and molecules. \cite{Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019}
Its popularity in the molecular electronic 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 computational cost \cite{Foerster_2011,Liu_2016,Wilhelm_2018,Forster_2021,Duchemin_2021} and somehow surprising accuracy for weakly correlated systems. \cite{Korbel_2014,vanSetten_2015,Caruso_2016,Hung_2017,vanSetten_2018,Bruneval_2021}

The idea behind the $GW$ approximation 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,Loos_2020e,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}
Multiple solution issues in $GW$ appears frequently, \cite{vanSetten_2015,Maggio_2017,Duchemin_2020} especially for orbitals that are energetically far from the Fermi level, such as in core ionized states. \cite{Golze_2018,Golze_2020}

In addition to obvious irregularities on potential energy surfaces that hampers the accurate determination of properties such as equilibrium bond lengths and harmonic vibrational frequencies, \cite{Loos_2020e,Berger_2021} one direct consequence of these discontinuities is the difficulty to converge (partially) self-consistent $GW$ calculations as the self-consistent procedure jumps erratically from one solution to the other even if convergence accelerator techniques such as DIIS are employed. \cite{Pulay_1980,Pulay_1982,Veril_2018}

It was shown that these problems can be tamed 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 $GW$ scheme, \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 solution but also the satellites at each iteration. \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 propose 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 starting point.
Throughout this article, $p$ and $q$ are general (spatial) orbitals, $i$, $j$, $k$, and $l$ denotes occupied orbitals, $a$, $b$, $c$, and $d$ are vacant orbitals, while $m$ labels single excitations $i \to a$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Downfolding: 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 is constituted by a hole (h) and a particle (p) term, as follows
\begin{equation}
\label{eq:SigC}
	\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 (that we enforce here for the sake of simplicity), 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 random-phase approximation (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}{ia} = \iint \MO{p}(\br_1) \MO{q}(\br_1) \frac{1}{\abs{\br_1 - \br_2}} \MO{i}(\br_2) \MO{a}(\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 weights are given by the value of the following renormalization factor
\begin{equation}
\label{eq:Z}
	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}$ has a large weight $Z_{p}$.
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}
which physically shows that the mean-field solution of unit weight is ``scattered'' by the effect of correlation in many solutions of smaller weights.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Upfolding: the linear $GW$ problem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The non-linear quasiparticle equation \eqref{eq:qp_eq} can be \textit{exactly} 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,Riva_2022}
For each orbital $p$, this yields a linear eigenvalue problem of the form
\begin{equation}
	\bH^{(p)} \cdot \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 $1 + O^2 V + O V^2$ (where $O$ and $V$ are the number of occupied and virtual orbitals, respectively), and this eigenvalue problem has to be solved for each orbital that one wishes to correct.
Note, however, that the blocks $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$ do not need to be recomputed for each orbital.
Because the renormalization factor \eqref{eq:Z} corresponds to the projection of the vector $\bc{}{(p,s)}$ onto the reference (or internal) space, the weight of a solution $(p,s)$ is given by the first coefficient of their corresponding eigenvector $\bc{}{(p,s)}$, \ie, 
\begin{equation}
\label{eq:Z_proj}
	Z_{p,s} = \qty[ c_{1}^{(p,s)} ]^{2}
\end{equation}

It is paramount to understand that diagonalizing $\bH^{(p)}$ [see Eq.~\eqref{eq:Hp}] is completely equivalent to solving the quasiparticle equation \eqref{eq:qp_eq}.
This can be further illustrated 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.

One can see this downfolding process as the construction of a frequency-dependent effective Hamiltonian where the internal space is composed by a single determinant of the 1h or 1p sector and the external (or outer) space by all the 2h1p and 2p1h configurations. \cite{Dvorak_2019a,Dvorak_2019b,Bintrim_2021a} 
The main mathematical difference between the two approaches is that, by diagonalizing Eq.~\eqref{eq:Hp}, one has directly access to the internal and external components of the eigenvectors associated with each quasiparticle and satellite, and not only their projection in the reference space as shown by Eq.~\eqref{eq:Z_proj}.

The element $\eps{p}{\HF}$ of $\bH^{(p)}$ [see Eq.~\eqref{eq:Hp}] corresponds to the (approximate) relative energy of the $(\Ne\pm1)$-electron reference determinant (compared to the $\Ne$-electron HF determinant) while the eigenvalues of the blocks $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$, which are $\eps{i}{\HF} - \Om{m}{\RPA}$ and $\eps{a}{\HF} + \Om{m}{\RPA}$ respectively, provide an estimate of the relative energy of the 2h1p and 2p1h determinants.
In some situations, one (or several) of these determinants from the external space may become of similar energy than the reference determinant, resulting in a vanishing denominator in the self-energy \eqref{eq:SigC}.
Hence, these two diabatic electronic configurations may cross and form an avoided crossing, and this outer-space determinant may be labeled as an intruder state.
As we shall see below, discontinuities, which are ubiquitous in molecular systems, arise in such scenarios.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{An illustrative 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 two-electron system with four spatial orbitals (one occupied and three virtuals).
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 discontinuities.
The downfolded and upfolded {\GOWO} schemes have been implemented in the electronic structure package QuAcK \cite{QuAcK} which is freely available at \url{https://github.com/pfloos/QuAcK}.
We denote as $\ket*{1\Bar{1}}$ the $\Ne$-electron ground-state determinant where the orbital 1 is occupied by one spin-up and one spin-down electron.
Similarly notations will be employed for the $(\Ne\pm1)$-electron configurations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	FIGURE 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
	\includegraphics[width=0.47\linewidth]{fig1a}
	\hspace{0.05\linewidth}
	\includegraphics[width=0.47\linewidth]{fig1b}
	\caption{
	\label{fig:H2}
	Selection of quasiparticle and satellite energies $\eps{p,s}{\GW}$ (solid lines) and their renormalization factor $Z_{p,s}$ (dashed lines) as functions of the internuclear distance $\RHH$ (in \si{\angstrom}) for the HOMO$+2$ ($p=3$) and HOMO$+3$ ($p=4$) orbitals of \ce{H2} at the {\GOWO}@HF/6-31G level.
	The quasiparticle solution (which corresponds to the solution with the largest weight) is represented as a thicker line.
}
\end{figure*}	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figure \ref{fig:H2} shows the evolution of the quasiparticle energy, energetically close-by satellites and their corresponding weights at the {\GOWO}@HF level as a function on the internuclear distance $\RHH$.
One can easily diagnose two problematic regions showing obvious discontinuities around $\RHH = \SI{1.2}{\angstrom}$ for the HOMO$+2$ ($p = 3$) and $\RHH = \SI{0.5}{\angstrom}$ for the HOMO$+3$ ($p = 4$).

Let us first look more closely at the region around $\RHH = \SI{1.2}{\angstrom}$.
As one can see, an avoided crossing is formed between two solutions of the quasiparticle equation.
Inspection of their corresponding eigenvectors reveals that the $(\Ne+1)$-electron determinants involved are the reference 1p determinant $\ket*{1\Bar{1}3}$ and an excited $(\Ne+1)$-electron of configuration $\ket*{12\Bar{2}}$ which becomes lower in energy than the reference determinant for $\RHH > \SI{1.2}{\angstrom}$.
The algorithm diabatically follows the reference determinant $\ket*{1\Bar{1}3}$, while it should be adiabatically switching to the $\ket*{12\Bar{2}}$ determinant which becomes lower in energy for $\RHH > \SI{1.1}{\angstrom}$.
This is not the $\Ne$-electron situation where one has to check if it is multireference, but for the $\Ne\pm1$-electron situations.

As similar scenario is at the play in the region around $\RHH = \SI{0.5}{\angstrom}$ but it now involves 3 electronic configurations: the $\ket*{1\Bar{1}4}$ reference determinant as well as two external determinants of configuration $\ket*{1\Bar{?}?}$ and $\ket*{1\Bar{?}?}$.
These forms two avoided crossings in rapid successions, which involves two discontinuities in the energy surface.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introducing regularized $GW$ methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
One way to alleviate these issues and to massively improve the convergence properties of self-consistent $GW$ calculations is to resort to a regularization of the self-energy without altering too much the quasiparticle energies.

From a general perspective, a regularized $GW$ self-energy reads 
\begin{equation}
\begin{split}
	\rSigC{p}(\omega;\eta) 
	& = \sum_{im} 2\ERI{pi}{m}^2 f_\eta(\omega - \eps{i}{\HF} + \Om{m}{\RPA})
	\\
	& + \sum_{am} 2\ERI{pa}{m}^2 f_\eta(\omega - \eps{a}{\HF} - \Om{m}{\RPA})
\end{split}
\end{equation}
where various choices for the ``regularizer'' $f_\eta$ are possible.
The main purpose of $f_\eta$ is to ensure that $\rSigC{p}(\omega;\eta)$ remains finite even if one of the denominators goes to zero.

The most common and well-established way of regularizing $\Sigma$ is via the simple regularizer $f_\eta(\Delta) = (\Delta \pm \eta)^{-1}$, a strategy somehow related to the imaginary shift used in multiconfigurational perturbation theory. \cite{Forsberg_1997}
Other choices are legitimate like the regularizers considered by Head-Gordon and coworkers within second-order M{\o}ller-Plesset theory. \cite{Lee_2018a,Shee_2021}

Our investigations have shown that the following regularizer
\begin{equation}
	f_\eta(\Delta) = \frac{1-e^{-2\Delta^2/\eta^2}}{\Delta}
\end{equation}
derived from the (second-order) perturbative analysis of the similarity renormalization group equations \cite{Wegner_1994,Glazek_1994,White_2002,Evangelista_2014} is particularly convenient and effective in the present context.
Increasing $\eta$ gradually integrates out states with denominators $\Delta$ larger than $\eta$ while the states with $\Delta \ll \eta$ are not decoupled from the reference space, hence avoiding intruder state problems. \cite{Li_2019a}
Of course, by construction, we have
\begin{equation}
	\lim_{\eta\to0} \rSigC{p}(\omega;\eta) = \SigC{p}(\omega) 
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	FIGURE 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
%	\includegraphics[width=\linewidth]{fig1a}
%	\includegraphics[width=\linewidth]{fig1b}
	\caption{
	\label{fig:H2_reg}
	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_reg} evidences how the regularization of the $GW$ self-energy diabatically linked the two solutions to get rid of the discontinuities.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Concluding remarks}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%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{
The authors thank Pina Romaniello for insightful discussions.
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{ufGW}
%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
