\documentclass[aps,prb,reprint,noshowkeys,linenumbers,superscriptaddress]{revtex4-1} \usepackage{subcaption} \usepackage{bm,graphicx,tabularx,array,booktabs,dcolumn,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,siunitx,mhchem} \usepackage[utf8]{inputenc} \usepackage[T1]{fontenc} \usepackage{txfonts} \usepackage[normalem]{ulem} \definecolor{hughgreen}{RGB}{0, 128, 0} \newcommand{\titou}[1]{\textcolor{red}{#1}} \newcommand{\hugh}[1]{\textcolor{hughgreen}{#1}} \newcommand{\trash}[1]{\textcolor{red}{\sout{#1}}} \newcommand{\trashHB}[1]{\textcolor{orange}{\sout{#1}}} \usepackage[ colorlinks=true, citecolor=blue, linkcolor=blue, filecolor=blue, urlcolor=blue, breaklinks=true ]{hyperref} \urlstyle{same} \newcommand{\ctab}{\multicolumn{1}{c}{---}} \newcommand{\latin}[1]{#1} %\newcommand{\latin}[1]{\textit{#1}} \newcommand{\ie}{\latin{i.e.}} \newcommand{\eg}{\latin{e.g.}} \newcommand{\etal}{\textit{et al.}} \newcommand{\mc}{\multicolumn} \newcommand{\fnm}{\footnotemark} \newcommand{\fnt}{\footnotetext} \newcommand{\mcc}[1]{\multicolumn{1}{c}{#1}} \newcommand{\mr}{\multirow} % operators \newcommand{\bH}{\mathbf{H}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bh}{\mathbf{h}} \newcommand{\bQ}{\mathbf{Q}} \newcommand{\bSig}{\mathbf{\Sigma}} \newcommand{\br}{\mathbf{r}} \newcommand{\bp}{\mathbf{p}} \newcommand{\cP}{\mathcal{P}} \newcommand{\cS}{\mathcal{S}} \newcommand{\cT}{\mathcal{T}} \newcommand{\cC}{\mathcal{C}} \newcommand{\PT}{\mathcal{PT}} \newcommand{\EPT}{E_{\PT}} \newcommand{\laPT}{\lambda_{\PT}} \newcommand{\EEP}{E_\text{EP}} \newcommand{\laEP}{\lambda_\text{EP}} \newcommand{\Ne}{N} % Number of electrons \newcommand{\Nn}{M} % Number of nuclei \newcommand{\hI}{\Hat{I}} \newcommand{\hH}{\Hat{H}} \newcommand{\hS}{\Hat{S}} \newcommand{\hT}{\Hat{T}} \newcommand{\hW}{\Hat{W}} \newcommand{\hV}{\Hat{V}} \newcommand{\hc}[2]{\Hat{c}_{#1}^{#2}} \newcommand{\hn}[1]{\Hat{n}_{#1}} \newcommand{\n}[1]{n_{#1}} \newcommand{\Dv}{\Delta v} \newcommand{\ra}{\rightarrow} \newcommand{\up}{\uparrow} \newcommand{\dw}{\downarrow} \newcommand{\updot}{% \mathrel{\ooalign{\hfil$\vcenter{ \hbox{$\scriptscriptstyle\bullet$}}$\hfil\cr$\uparrow$\cr} }% } \newcommand{\dwdot}{% \mathrel{\ooalign{\hfil$\vcenter{ \hbox{$\scriptscriptstyle\bullet$}}$\hfil\cr$\downarrow$\cr} }% } \newcommand{\vac}{% \mathrel{\ooalign{\hfil$\vcenter{ \hbox{$\scriptscriptstyle\bullet$}}$\hfil\cr$ $\cr} }% } \newcommand{\uddot}{% \mathrel{\ooalign{\hfil$\vcenter{ \hbox{$\scriptscriptstyle\bullet$}}$\hfil\cr$\uparrow\downarrow$\cr} }% } % Center tabularx columns \newcolumntype{Y}{>{\centering\arraybackslash}X} % HF rotation angles \newcommand{\ta}{\theta_{\alpha}} \newcommand{\tb}{\theta_{\beta}} % Some constants \renewcommand{\i}{\mathrm{i}} % Imaginary unit \newcommand{\e}{\mathrm{e}} % Euler number \newcommand{\rc}{r_{\text{c}}} \newcommand{\lc}{\lambda_{\text{c}}} \newcommand{\lep}{\lambda_{\text{EP}}} % Some energies \newcommand{\Emp}{E_{\text{MP}}} % Blackboard bold \newcommand{\bbR}{\mathbb{R}} \newcommand{\bbC}{\mathbb{C}} \newcommand{\Lup}{\mathcal{L}^{\uparrow}} \newcommand{\Ldown}{\mathcal{L}^{\downarrow}} \newcommand{\Lsi}{\mathcal{L}^{\sigma}} \newcommand{\Rup}{\mathcal{R}^{\uparrow}} \newcommand{\Rdown}{\mathcal{R}^{\downarrow}} \newcommand{\Rsi}{\mathcal{R}^{\sigma}} \newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France.} \newcommand{\UCAM}{Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, U.K.} \newcommand{\UOX}{Physical and Theoretical Chemical Laboratory, Department of Chemistry, University of Oxford, Oxford, OX1 3QZ, U.K.} \begin{document} \title{Perturbation Theory in the Complex Plane: Exceptional Points and Where to Find Them} \author{Antoine \surname{Marie}} \affiliation{\LCPQ} \author{Hugh G.~A.~\surname{Burton}} \email{hugh.burton@chem.ox.ac.uk} \affiliation{\UOX} \author{Pierre-Fran\c{c}ois \surname{Loos}} \email{loos@irsamc.ups-tlse.fr} \affiliation{\LCPQ} \begin{abstract} In this review, we explore the extension of quantum chemistry in the complex plane and its link with perturbation theory. We observe that the physics of a quantum system is intimately connected to the position of energy singularities in the complex plane, known as exceptionnal points. After a presentation of the fundamental notions of quantum chemistry in the complex plane, such as the mean-field Hartree--Fock approximation and Rayleigh-Schr\"odinger perturbation theory, and their illustration with the ubiquitous (symmetric) Hubbard dimer at half filling, we provide a historical overview of the various research activities that have been performed on the physics of singularities. In particular, we highlight the seminal work of several research groups on the convergence behaviour of perturbative series obtained within M{\o}ller--Plesset perturbation theory and its apparent link with quantum phase transitions. \end{abstract} \maketitle \raggedbottom \tableofcontents %%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} \label{sec:intro} %%%%%%%%%%%%%%%%%%%%%%% Due to the ubiquitous influence of processes involving electronic states in physics, chemistry, and biology, their faithful description from first principles has been one of the grand challenges faced by theoretical chemists since the dawn of computational chemistry. Accurately predicting ground- and excited-state energies (hence excitation energies) is particularly valuable in this context, and it has concentrated most of the efforts within the community. An armada of theoretical and computational methods have been developed to this end, each of them being plagued by its own flaws. The fact that none of these methods is successful in every chemical scenario has encouraged chemists to carry on the development of new methodologies, their main goal being to get the most accurate energies (and properties) at the lowest possible computational cost in the most general context. One common feature of all these methods is that they rely on the notion of quantised energy levels of Hermitian quantum mechanics, in which the different electronic states of a molecule or an atom are energetically ordered, the lowest being the ground state while the higher ones are excited states. Within this quantised paradigm, electronic states look completely disconnected from one another. Many current methods study excited states using only ground-state information, creating a ground-state bias that leads to incorrect excitation energies. However, one can gain a different perspective on quantisation extending quantum chemistry into the complex domain. In a non-Hermitian complex picture, the energy levels are \textit{sheets} of a more complicated topological manifold called \textit{Riemann surface}, and they are smooth and continuous \textit{analytic continuation} of one another. In other words, our view of the quantised nature of conventional Hermitian quantum mechanics arises only from our limited perception of the more complex and profound structure of its non-Hermitian variant. \cite{MoiseyevBook,BenderPTBook} The realisation that ground and excited states both emerge from one single mathematical structure with equal importance suggests that excited-state energies can be computed from first principles in their own right. One could then exploit the structure of these Riemann surfaces to develop methods that directly target excited-state energies without needing ground-state information. \cite{Burton_2019,Burton_2019a} By analytically continuing the electronic energy $E(\lambda)$ in the complex domain (where $\lambda$ is a coupling parameter), the ground and excited states of a molecule can be smoothly connected. This connection is possible because by extending real numbers to the complex domain, the ordering property of real numbers is lost. Hence, electronic states can be interchanged away from the real axis since the concept of ground and excited states has been lost. Amazingly, this smooth and continuous transition from one state to another has recently been experimentally realised in physical settings such as electronics, microwaves, mechanics, acoustics, atomic systems and optics. \cite{Bittner_2012,Chong_2011,Chtchelkatchev_2012,Doppler_2016,Guo_2009,Hang_2013,Liertzer_2012,Longhi_2010,Peng_2014, Peng_2014a,Regensburger_2012,Ruter_2010,Schindler_2011,Szameit_2011,Zhao_2010,Zheng_2013,Choi_2018,El-Ganainy_2018} Exceptional points (EPs) are branch point singularities where two (or more) states become exactly degenerate. \cite{MoiseyevBook,Heiss_1988,Heiss_1990,Heiss_1999,Berry_2011,Heiss_2012,Heiss_2016,Benda_2018} They are the non-Hermitian analogs of conical intersections, \cite{Yarkony_1996} which are ubiquitous in non-adiabatic processes and play a key role in photochemical mechanisms. In the case of auto-ionising resonances, EPs have a role in deactivation processes similar to conical intersections in the decay of bound excited states. \cite{Benda_2018} Although Hermitian and non-Hermitian Hamiltonians are closely related, the behaviour of their eigenvalues near degeneracies is starkly different. For example, encircling non-Hermitian degeneracies at EPs leads to an interconversion of states, and two loops around the EP are necessary to recover the initial energy. \cite{MoiseyevBook,Heiss_2016,Benda_2018} Additionally, the wave function picks up a geometric phase (also known as Berry phase \cite{Berry_1984}) and four loops are required to recover the initial wave function. In contrast, encircling Hermitian degeneracies at conical intersections only introduces a geometric phase while leaving the states unchanged. More dramatically, whilst eigenvectors remain orthogonal at conical intersections, at non-Hermitian EPs the eigenvectors themselves become equivalent, resulting in a \textit{self-orthogonal} state. \cite{MoiseyevBook} More importantly here, although EPs usually lie off the real axis, these singular points are intimately related to the convergence properties of perturbative methods and avoided crossing on the real axis are indicative of singularities in the complex plane. \cite{BenderBook,Olsen_1996,Olsen_2000,Olsen_2019,Mihalka_2017a,Mihalka_2017b,Mihalka_2019} \titou{The use of non-Hermitian Hamiltonians in quantum chemistry has a long history; these Hamiltonians have been used extensively as a method for describing metastable resonance phenomena. \cite{MoiseyevBook} Through a complex-scaling of the electronic or atomic coordinates,\cite{Moiseyev_1998} or by introducing a complex absorbing potential,\cite{Riss_1993,Ernzerhof_2006,Benda_2018} outgoing resonance states are transformed into square-integrable wave functions that allow the energy and lifetime of the resonance to be computed. We refer the interested reader to the excellent book of Moiseyev for a general overview. \cite{MoiseyevBook}} \titou{Discussion around the different types of singularities in complex analysis. At a singular point, a function and/or its derivatives becomes infinite or undefined (hence non analytic). One very common type of singularities (belonging to the family of isolated singularities) are poles where the function behaves $1/(\lambda - \lambda_c)^n$ where $n \in \mathbb{N}^*$ is the order of the pole. Another class of singularities are branch points resulting from a multi-valued function such as a square root or a logarithm function and usually implying the presence of so-called branch cuts where the function ``jumps'' from one value to another. Yet another family of singularities are formed by critical points which lie on the real axis, have more complicated functional forms and where the nature of the function undergoes a sudden transition. } %%%%%%%%%%%%%%%%%%%%%%% \section{Exceptional Points in Electronic Structure} \label{sec:EPs} %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% \subsection{Time-Independent Schr\"odinger Equation} \label{sec:TDSE} %%%%%%%%%%%%%%%%%%%%%%% Within the Born-Oppenheimer approximation, the exact molecular Hamiltonian with $\Ne$ electrons and $\Nn$ (clamped) nuclei is defined for a given nuclear framework as \begin{equation}\label{eq:ExactHamiltonian} \hH(\vb{R}) = - \frac{1}{2} \sum_{i}^{\Ne} \grad_i^2 - \sum_{i}^{\Ne} \sum_{A}^{\Nn} \frac{Z_A}{\abs{\vb{r}_i-\vb{R}_A}} + \sum_{i2t$, the closed-shell orbital restriction prevents RHF from modelling the correct physics with the two electrons on opposing sites. %%% FIG 3 (?) %%% % Analytic Continuation of HF %%%%%%%%%%%%%%%%% \begin{figure*}[t] \begin{subfigure}{0.49\textwidth} \includegraphics[height=0.65\textwidth,trim={0pt 0pt 0pt -35pt},clip]{HF_cplx_angle} \subcaption{\label{subfig:UHF_cplx_angle}} \end{subfigure} \begin{subfigure}{0.49\textwidth} \includegraphics[height=0.65\textwidth]{HF_cplx_energy} \subcaption{\label{subfig:UHF_cplx_energy}} \end{subfigure} \caption{% (\subref{subfig:UHF_cplx_angle}) Real component of the UHF angle $\ta^{\text{UHF}}$ for $\lambda \in \bbC$. Symmetry-broken solutions correspond to individual sheets and become equivalent at the \textit{quasi}-EP $\lambda_{\text{c}}$ (black dot). The RHF solution is independent of $\lambda$, giving the constant plane at $\pi/2$. (\subref{subfig:UHF_cplx_energy}) The corresponding HF energy surfaces show a non-analytic point at the \textit{quasi}-EP. \label{fig:HF_cplx}} \end{figure*} %%%%%%%%%%%%%%%%% As the on-site repulsion is increased from 0, the HF approximation reaches a critical value at $U=2t$ where a symmetry-broken UHF solution appears with a lower energy than the RHF one. Note that the RHF wave function remains a genuine solution of the HF equations for $U \ge 2t$, but corresponds to a saddle point of the HF energy rather than a minimum. This critical point is analogous to the infamous Coulson--Fischer point identified in the hydrogen dimer.\cite{Coulson_1949} For $U \ge 2t$, the optimal orbital rotation angles for the UHF orbitals become \begin{subequations} \begin{align} \ta^\text{UHF} & = \arctan (-\frac{2t}{\sqrt{U^2 - 4t^2}}), \label{eq:ta_uhf} \\ \tb^\text{UHF} & = \arctan (+\frac{2t}{\sqrt{U^2 - 4t^2}}), \label{eq:tb_uhf} \end{align} \end{subequations} with the corresponding UHF ground-state energy (Fig.~\ref{fig:HF_real}) \begin{equation} E_\text{UHF} \equiv E_\text{HF}(\ta^\text{UHF}, \tb^\text{UHF}) = - \frac{2t^2}{U}. \end{equation} Time-reversal symmetry dictates that this UHF wave function must be degenerate with its spin-flipped pair, obtained by swapping $\ta^{\text{UHF}}$ and $\tb^{\text{UHF}}$ in Eqs.~\eqref{eq:ta_uhf} and \eqref{eq:tb_uhf}. This type of symmetry breaking is also called a spin-density wave in the physics community as the system ``oscillates'' between the two symmetry-broken configurations. \cite{GiulianiBook} Symmetry breaking can also occur in RHF theory when a charge-density wave is formed from an oscillation between the two closed-shell configurations with both electrons localised on one site or the other.\cite{StuberPaldus,Fukutome_1981} %============================================================% \subsection{Self-Consistency as a Perturbation} %OR {Complex adiabatic connection} %============================================================% % INTRODUCE PARAMETRISED FOCK HAMILTONIAN The inherent non-linearity in the Fock eigenvalue problem arises from self-consistency in the HF approximation, and is usually solved through an iterative approach.\cite{Roothaan_1951,Hall_1951} Alternatively, the non-linear terms arising from the Coulomb and exchange operators can be considered as a perturbation from the core Hamiltonian \eqref{eq:Hcore} by introducing the transformation $U \rightarrow \lambda\, U$, giving the parametrised Fock operator \begin{equation} \Hat{f}(\vb{x} ; \lambda) = \Hat{h}(\vb{x}) + \lambda\, \Hat{v}_\text{HF}(\vb{x}). \end{equation} The orbitals in the reference problem $\lambda=0$ correspond to the symmetry-pure eigenfunctions of the one-electron core Hamiltonian, while self-consistent solutions at $\lambda = 1$ represent the orbitals of the true HF solution. % INTRODUCE COMPLEX ANALYTIC-CONTINUATION For real $\lambda$, the self-consistent HF energies at given (real) $U$ and $t$ values in the Hubbard dimer directly mirror the energies shown in Fig.~\ref{fig:HF_real}, with coalesence points at \begin{equation} \lambda_{\text{c}} = \pm \frac{2t}{U}. \label{eq:scaled_fock} \end{equation} In contrast, when $\lambda$ becomes complex, the HF equations become non-Hermitian and each HF solutions can be analytically continued for all $\lambda$ values using the holomorphic HF approach.\cite{Hiscock_2014,Burton_2016,Burton_2018} Remarkably, the coalescence point in this analytic continuation emerges as a \textit{quasi}-EP on the real $\lambda$ axis (Fig.~\ref{fig:HF_cplx}), where the different HF solutions become equivalent but not self-orthogonal.\cite{Burton_2019} By analogy with perturbation theory, the regime where this \textit{quasi}-EP occurs within $\lambda_{\text{c}} \le 1$ can be interpreted as an indication that the symmetry-pure reference orbitals no longer provide a qualitatively accurate representation for the true HF ground state at $\lambda = 1$. For example, in the Hubbard dimer with $U > 2t$, one finds $\lambda_{\text{c}} < 1$ and the symmetry-pure orbitals do not provide a good representation of the HF ground state. In contrast, $U < 2t$ yields $\lambda_{\text{c}} > 1$ and corresponds to the regime where the HF ground state is correctly represented by symmetry-pure orbitals. % COMPLEX ADIABATIC CONNECTION We have recently shown that the complex scaled Fock operator Eq.~\eqref{eq:scaled_fock} also allows states of different symmetries to be interconverted by following a well-defined contour in the complex $\lambda$-plane.\cite{Burton_2019} In particular, by slowly varying $\lambda$ in a similar (yet different) manner to an adiabatic connection in density-functional theory,\cite{Langreth_1975,Gunnarsson_1976,Zhang_2004} a ground-state wave function can be ``morphed'' into an excited-state wave function via a stationary path of HF solutions. This novel approach to identifying excited-state wave functions demonstrates the fundamental role of \textit{quasi}-EPs in determining the behaviour of the HF approximation. %\titou{In a recent paper, \cite{Burton_2019} using holomorphic Hartree--Fock (h-HF) \cite{Hiscock_2014,Burton_2018,Burton_2016} as an analytic continuation of conventional HF theory, we have demonstrated, on a simple model, that one can interconvert states of different symmetries and natures by following well-defined contours in the complex $\lambda$-plane, where $\lambda$ is the strength of the electron-electron interaction (see Fig.~\ref{fig:iAC}). %In particular, by slowly varying $\lambda$ in a similar (yet different) manner to an adiabatic connection in density-functional theory, \cite{Langreth_1975,Gunnarsson_1976,Zhang_2004} one can ``morph'' a ground-state wave function into an excited-state wave function via a stationary path of HF solutions. \cite{Seidl_2018} %In such a way, we could obtain a doubly-excited state wave function starting from the ground state wave function, a process which is not as easy as one might think. \cite{Gilbert_2008,Thom_2008,Shea_2018} %One of the fundamental discovery we made was that Coulson-Fischer points (where multiple symmetry-broken solutions coalesce) play a central role and can be classified as \textit{quasi}-exceptional points, as the wave functions do not become self-orthogonal. %The findings reported in Ref.~\onlinecite{Burton_2019} represent the very first study of non-Hermitian quantum mechanics for the exploration of multiple solutions at the HF level. %It perfectly illustrates the deeper topology of electronic states revealed using a complex-scaled electron-electron interaction. %Through the introduction of non-Hermiticity, we have provided a more general framework in which the complex and diverse characteristics of multiple solutions can be explored and understood.} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{M{\o}ller--Plesset Theory in the Complex Plane} \label{sec:MP} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %=====================================================% \subsection{Background Theory} %=====================================================% In electronic structure, the HF Hamiltonian \eqref{eq:HFHamiltonian} is often used as the zeroth-order Hamiltonian to define M\o{}ller--Plesset (MP) perturbation theory.\cite{Moller_1934} This approach can recover a large proportion of the electron correlation energy,\cite{Lowdin_1955a,Lowdin_1955b,Lowdin_1955c} and provides the foundation for numerous post-HF approximations. With the MP partitioning, the parametrised perturbation Hamiltonian becomes \begin{multline}\label{eq:MPHamiltonian} \hH(\lambda) = \sum_{i}^{N} \qty[ - \frac{\grad_i^2}{2} - \sum_{A}^{M} \frac{Z_A}{\abs{\vb{r}_i-\vb{R}_A}} ] \\ + (1-\lambda) \sum_{i}^{N} v^{\text{HF}}(\vb{x}_i) + \lambda\sum_{i 1$) and $4.5$ (where $r_c < 1$). The Riemann surfaces associated with the exact energies of the RMP Hamiltonian \eqref{eq:H_RMP} are also represented for these two values of $U/t$ as functions of $\lambda$. \label{fig:RMP}} \end{figure*} The behaviour of the RMP and UMP series observed in \ce{H2} can also be illustrated by considering the analytic Hubbard dimer with a complex-valued perturbation strength. In this system, the stretching of the \ce{H\bond{-}H} bond is directly mirrored by an increase in the electron correlation $U/t$. Using the ground-state RHF reference orbitals leads to the parametrised RMP Hamiltonian \begin{widetext} \begin{equation} \label{eq:H_RMP} \bH_\text{RMP}\qty(\lambda) = \begin{pmatrix} -2t + U - \lambda U/2 & 0 & 0 & \lambda U/2 \\ 0 & U - \lambda U/2 & \lambda U/2 & 0 \\ 0 & \lambda U/2 & U - \lambda U/2 & 0 \\ \lambda U/2 & 0 & 0 & 2t + U - \lambda U/2 \\ \end{pmatrix}, \end{equation} \end{widetext} which yields the ground-state energy \begin{equation} \label{eq:E0MP} E_{-}(\lambda) = U - \frac{\lambda U}{2} - \frac{1}{2} \sqrt{(4t)^2 + \lambda ^2 U^2}. \end{equation} From this expression, the EPs can be identified as $\lep = \pm \i 4t / U$, giving the radius of convergence \begin{equation} \rc = \abs{\frac{4t}{U}}. \end{equation} Remarkably, these EPs are identical to the exact EPs discussed in Sec.~\ref{sec:example}. The Taylor expansion of the RMP energy can then be evaluated to obtain the $k$th-order MP correction \begin{equation} E_\text{RMP}^{(k)} = U \delta_{0,k} - \frac{1}{2} \frac{U^k}{(4t)^{k-1}} \mqty( 1/2 \\ k/2). \end{equation} %with %\begin{equation} % E_{\text{MP}n}(\lambda) = \sum_{k=0}^n E_\text{MP}^{(k)} \lambda^k. %\end{equation} %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % RADIUS OF CONVERGENCE PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[htb] \includegraphics[width=\linewidth]{RadConv} \caption{ Radius of convergence $r_c$ for the RMP ground state (red), the UMP ground state (blue), and the UMP excited state (orange) series as functions of the ratio $U/t$. \label{fig:RadConv}} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%% The RMP series is convergent for $U = 3.5\,t$ with $\rc > 1$, as illustrated for the individual terms at each order of perturbation in Fig.~\ref{subfig:RMP_cvg}. In contrast, for $U = 4.5t$ one finds $\rc < 1$, and the RMP series becomes divergent. The corresponding Riemann surfaces for $U = 3.5\,t$ and $4.5\,t$ are shown in Figs.~\ref{subfig:RMP_3.5} and \ref{subfig:RMP_4.5}, respectively, with the single EP at $\lep$ (black dot) and the radius of convergence indicated by the vertical cylinder of unit radius. For the divergent case, the $\lep$ lies inside this cylinder of convergence, while in the convergent case $\lep$ lies outside this cylinder. In both cases, the EP connects the ground state with the doubly-excited state, and thus the convergence behaviour for the two states using the ground-state RHF orbitals is identical. % HGAB: This cannot be relevant here as the single-excitations don't couple to either ground or excited state. %The convergent and divergent series start to differ at fourth order, corresponding to the lowest-order contribution %%from the single excitations.\cite{Lepetit_1988} %%% FIG 3 %%% \begin{figure*} \begin{subfigure}{0.32\textwidth} \includegraphics[height=0.75\textwidth]{fig3a} \subcaption{\label{subfig:UMP_3} $U/t = 3$} \end{subfigure} % \begin{subfigure}{0.32\textwidth} \includegraphics[height=0.75\textwidth]{fig3b} \subcaption{\label{subfig:UMP_cvg}} \end{subfigure} % \begin{subfigure}{0.32\textwidth} \includegraphics[height=0.75\textwidth]{fig3c} \subcaption{\label{subfig:UMP_7} $U/t = 7$} \end{subfigure} \caption{ Convergence of the UMP series as a function of the perturbation order $n$ for the Hubbard dimer at $U/t = 3$ and $7$. The Riemann surfaces associated with the exact energies of the UMP Hamiltonian \eqref{eq:H_UMP} are also represented for these two values of $U/t$ as functions of $\lambda$. \label{fig:UMP}} \end{figure*} The behaviour of the UMP series is more subtle than the RMP series as the spin-contamination in the wave function introduces additional coupling between the singly- and doubly-excited configurations. Using the ground-state UHF reference orbitals in the Hubbard dimer yields the parametrised UMP Hamiltonian \begin{widetext} \begin{equation} \label{eq:H_UMP} \bH_\text{UMP}\qty(\lambda) = \begin{pmatrix} -2t^2 \lambda/U & 0 & 0 & 2t^2 \lambda/U \\ 0 & U - 2t^2 \lambda/U & 2t^2\lambda/U & 2t \sqrt{U^2 - (2t)^2} \lambda/U \\ 0 & 2t^2\lambda/U & U - 2t^2 \lambda/U & -2t \sqrt{U^2 - (2t)^2} \lambda/U \\ 2t^2 \lambda/U & 2t \sqrt{U^2 - (2t)^2} \lambda/U & -2t \sqrt{U^2 - (2t)^2} \lambda/U & 2U(1-\lambda) + 6t^2\lambda/U \\ \end{pmatrix}. \end{equation} \end{widetext} While a closed-form expression for the ground-state energy exists, it is cumbersome and we eschew reporting it. Instead, the radius of convergence of the UMP series can be obtained numerically as a function of $U/t$, as shown in Fig.~\ref{fig:RadConv}. These numerical values reveal that the UMP ground-state series has $\rc > 1$ for all $U/t$ and always converges. However, in the strong correlation limit (large $U$), this radius of convergence tends to unity, indicating that the corresponding UMP series becomes increasingly slow. Furthermore, the doubly-excited state using the ground-state UHF orbitals has $\rc < 1$ for almost any value of $U/t$, reaching the limiting value of $1/2$ for $U/t \rightarrow \infty$, and thus the excited-state UMP series will always diverge. % DISCUSSION OF UMP RIEMANN SURFACES The convergence behaviour can be further elucidated by considering the full structure of the UMP energies in the complex $\lambda$-plane. These Riemann surfaces are illustrated for $U = 3t$ and $7t$ alongside the perturbation terms at each order in Fig.~\ref{subfig:UMP_cvg}. At $U = 3t$, the RMP series is convergent, while RMP becomes divergent for $U=7t$. The ground-state UMP expansion is convergent in both cases, although the rate of convergence is significantly slower for larger $U/t$ as the radius of convergence becomes increasingly close to one (Fig.~\ref{fig:RadConv}). % EFFECT OF SYMMETRY BREAKING As the UHF orbitals break the molecular symmetry, new coupling terms emerge between the electronic states that cause fundamental changes to the structure of EPs in the complex $\lambda$-plane. For example, while the RMP energy shows only one EP between the ground state and the doubly-excited state (Fig.~\ref{fig:RMP}), the UMP energy has two EPs: one connecting the ground state with the singly-excited open-shell singlet, and the other connecting this single excitation to the doubly-excited second excitation (Fig.~\ref{fig:UMP}). This new ground-state EP always appears outside the unit cylinder and guarantees convergence of the ground-state energy. However, the excited-state EP is moved within the unit cylinder and causes the convergence of the excited-state UMP series to deteriorate. Our interpretation of this effect is that the symmetry-broken orbital optimisation has redistributed the strong coupling between the ground- and doubly-excited states into weaker couplings between all states, and has thus sacrificed convergence of the excited-state series so that the ground-state convergence can be maximised. Since the UHF ground state already provides a good approximation to the exact energy, the ground-state sheet of the UMP energy is relatively flat and the corresponding EP in the Hubbard dimer always lies outside the unit cylinder. The slow convergence observed in stretched \ce{H2}\cite{Gill_1988} can then be seen as this EP moves increasingly close to the unit cylinder at large $U/t$ and $\rc$ approaches one (from above). Furthermore, the majority of the UMP expansion in this regime is concerned with removing spin-contamination from the wave function rather than improving the energy. It is well-known that the spin-projection needed to remove spin-contamination can require non-linear combinations of highly-excited determinants,\cite{Lowdin_1955c} and thus it is not surprising that this process proceeds very slowly as the perturbation order is increased. %The convergence of the UMP as a function of the ratio $U/t$ is shown in Fig.~\ref{subfig:UMP_cvg} for two specific values: the first ($U = 3t$) is well within the RMP convergence region, while the second ($U = 7t$) falls outside. %Note that in the case of UMP, there are now two pairs of EPs as the open-shell singlet now couples strongly with both the ground and doubly-excited states. %This has the clear tendency to move away from the origin the EP dictating the convergence of the ground-state energy, while deteriorating the convergence properties of the excited-state energy. %For $U = 3t$ (see Fig.~\ref{subfig:UMP_3}), the ground-state energy is remarkably flat since the UHF energy is already a pretty good estimate of the exact energy thanks to the symmetry-breaking process. %Most of the UMP expansion is actually correcting the spin-contamination in the wave function. %For $U = 7t$ (see Fig.~\ref{subfig:UMP_7}), we are well towards the strong correlation regime, where we see that the UMP series is slowly convergent while RMP diverges. %We see a single EP on the ground-state surface which falls just outside (maybe on?) the radius of convergence. %An EP this close to the radius of convergence gives an increasingly slow convergence of the UMP series, but it will converge eventually as observed in Fig.~\ref{subfig:UMP_cvg}. %On the other hand, there is an EP on the excited energy surface that is well within the radius of convergence. %We can therefore say that the use of a symmetry-broken UHF wave function can retain a convergent ground-state perturbation series %at the expense of a divergent excited-state perturbation series. (Note: the orbitals are not optimised for excited-state here). %In contrast, the RMP expansion was always convergent for the open-shell excited state (which was a single CSF) while %the radius of convergence for the doubly-excited state was identical to the ground-state as this was the only exceptional point. %==========================================% \subsection{Classifying Types of Convergence} % Behaviour} % Further insights from a two-state model} %==========================================% % CREMER AND HE As computational implementations of higher-order MP terms improved, the systematic investigation of convergence behaviour in a broader class of molecules became possible. Cremer and He introduced an efficient MP6 approach and used it to analyse the RMP convergence of 29 atomic and molecular systems with respect to the FCI energy.\cite{Cremer_1996} They established two general classes: ``class A'' systems that exhibit monotonic convergence; and ``class B'' systems for which convergence is erratic after initial oscillations. %Their system set contains stretched molecules as well as molecules at their equilibrium geometry for various basis sets. By analysing the different cluster contributions to the MP energy terms, they proposed that class A systems generally include well-separated and weakly correlated electron pairs, while class B systems are characterised by dense electron clustering in one or more spatial regions.\cite{Cremer_1996} %\textit{``Class A systems are characterised by electronic structures with well-separated electron pairs while class B systems are characterized by electronic structures with electron clustering in one or more regions.''} %Moreover, they analysed the contribution of the triple (T) excitations to the MP4, MP5 and MP6 energies next to the single, double and quadruple (SDQ) excitations contribution. In class A systems, they showed that the majority of the correlation energy arises from pair correlation, with little contribution from triple excitations. On the other hand, triple excitations have an important contribution in class B systems, including providing orbital relaxation, and these contributions lead to oscillations of the total correlation energy. %This observation on the contribution to the MP$n$ energy corroborates the electronic structure discussed above. Using these classifications, Cremer and He then introduced simple extrapolation formulas for estimating the exact correlation energy $\Delta E$ using terms up to MP6\cite{Cremer_1996} \begin{subequations} \begin{align} \Delta E_{\text{A}} &= \Emp^{(2)} + \Emp^{(3)} + \Emp^{(4)} + \frac{\Emp^{(5)}}{1 - (\Emp^{(6)} / \Emp^{(5)})}, \\[5pt] \Delta E_{\text{B}} &= \Emp^{(2)} + \Emp^{(3)} + \qty(\Emp^{(4)} + \Emp^{(5)}) \exp(\Emp^{(6)} / \Emp^{(5)}). \end{align} \end{subequations} %As one can only compute the first terms of the MP series, a smart way of getting more accurate results is to use extrapolation formula, \ie, estimating the limit of the series with only few terms. %Cremer and He proved that using specific extrapolation formulas of the MP series for class A and class B systems improves the precision of the results compared to the formula used without resorting to classes. \cite{Cremer_1996} These class-specific formulas reduced the mean absolute error from the FCI correlation energy by a factor of four compared to previous class-independent extrapolations, highlighting how one can leverage a deeper understanding of MP convergence to improve estimates of the correlation energy at lower computational costs. In Section~\ref{sec:Resummation}, we consider more advanced extrapolation routines that take account of EPs in the complex $\lambda$-plane. %The mean absolute deviation taking the FCI correlation energies as reference is $0.3$ millihartree with the class-specific formula whereas the deviation increases to 12 millihartree using the general formula. %Even if there were still shaded areas in their analysis and that their classification was incomplete, the work of Ref.~\onlinecite{Cremer_1996} clearly evidenced that understanding the origin of the different modes of convergence could potentially lead to a more rationalised use of MP perturbation theory and, hence, to more accurate correlation energy estimates. In the late 90's, Olsen \etal\ discovered an even more concerning behaviour of the MP series. \cite{Olsen_1996} They showed that the series could be divergent even in systems that were considered to be well understood, such as \ce{Ne} or the \ce{HF} molecule. \cite{Olsen_1996, Christiansen_1996} Cremer and He had already studied these two systems and classified them as \textit{class B} systems.\cite{Cremer_1996} However, Olsen and co-workers performed their analysis in larger basis sets containing diffuse functions, finding that the corresponding MP series becomes divergent at (very) high order. The discovery of this divergent behaviour is particularly worrying as large basis sets are required to get meaningful and accurate energies.\cite{Loos_2019d,Giner_2019} Furthermore, diffuse functions are particularly important for anions and/or Rydberg excited states, where the wave function is inherently more diffuse than the ground state.\cite{Loos_2018a,Loos_2020a} Olsen \etal\ investigated the causes of these divergences and the different types of convergence by analysing the relation between the dominant singularity (\ie, the closest singularity to the origin) and the convergence behaviour of the series.\cite{Olsen_2000} Their analysis is based on Darboux's theorem: \cite{Goodson_2011} \begin{quote} \textit{``In the limit of large order, the series coefficients become equivalent to the Taylor series coefficients of the singularity closest to the origin. ''} \end{quote} Following this theory, a singularity in the unit circle is designated as an intruder state, with a front-door (or back-door) intruder state if the real part of the singularity is positive (or negative). Using their observations in Ref.~\onlinecite{Olsen_1996}, Olsen and collaborators proposed a simple method that performs a scan of the real axis to detect the avoided crossing responsible for the dominant singularities in the complex plane. \cite{Olsen_2000} By modelling this avoided crossing using a two-state Hamiltonian, one can obtain an approximation for the dominant singularities as the EPs of the $2\times2$ matrix \begin{equation} \label{eq:Olsen_2x2} \underbrace{\mqty(\alpha & \delta \\ \delta & \beta )}_{\bH} = \underbrace{\mqty(\alpha + \alpha_{\text{s}} & 0 \\ 0 & \beta + \beta_{\text{s}} )}_{\bH^{(0)}} + \underbrace{\mqty( -\alpha_{\text{s}} & \delta \\ \delta & - \beta_{\text{s}})}_{\bV}, \end{equation} where the diagonal matrix is the unperturbed Hamiltonian matrix $\bH^{(0)}$ with level shifts $\alpha_{\text{s}}$ and $\beta_{\text{s}}$, and $\bV$ represents the perturbation. The authors first considered molecules with low-lying doubly-excited states with the same spatial and spin symmetry as the ground state. \cite{Olsen_2000} In these systems, the exact wave function has a non-negligible contribution from the doubly-excited states, and thus the low-lying excited states are likely to become intruder states. For \ce{CH_2} in a diffuse, yet rather small basis set, the series is convergent at least up to the 50th order, and the dominant singularity lies close (but outside) the unit circle, causing slow convergence of the series. These intruder-state effects are analogous to the EP that dictates the convergence behaviour of the RMP series for the Hubbard dimer (Fig.~\ref{fig:RMP}). Furthermore, the authors demonstrated that the divergence for \ce{Ne} is due to a back-door intruder state that arise when the ground state undergoes sharp avoided crossings with highly diffuse excited states. %They used their two-state model on this avoided crossings and the model was actually predicting the divergence of the series. %They concluded that the divergence of the series was due to the interaction with a highly diffuse excited state. This divergence is related to a more fundamental critical point in the MP energy surface that we will discuss in Section~\ref{sec:MP_critical_point}. Finally, Ref.~\onlinecite{Olsen_1996} proved that the extrapolation formulas of Cremer and He \cite{Cremer_1996} are not mathematically motivated when considering the complex singularities causing the divergence, and therefore cannot be applied for all systems. For example, \ce{HF} contains both back-door intruder states and low-lying doubly-excited states that result in alternating terms up to 10th order. The series becomes monotonically convergent at higher orders since the two pairs of singularities are approximately the same distance from the origin. More recently, this two-state model has been extended to non-symmetric Hamiltonians as\cite{Olsen_2019} \begin{equation} \underbrace{\mqty(\alpha & \delta_1 \\ \delta_2 & \beta)}_{\bH} = \underbrace{\mqty(\alpha & 0 \\ 0 & \beta + \gamma )}_{\bH^{(0)}} + \underbrace{\mqty( 0 & \delta_2 \\ \delta_1 & - \gamma)}_{\bV}. \end{equation} This extension allows various choices of perturbation to be analysed, including coupled cluster perturbation expansions \cite{Pawlowski_2019a,Pawlowski_2019b,Pawlowski_2019c,Pawlowski_2019d,Pawlowski_2019e} and other non-Hermitian perturbation methods. Note that new forms of perturbation expansions only occur when the sign of $\delta_1$ and $\delta_2$ differ. Using these non-Hermitian two-state model, the convergence of a perturbation series can be characterised according to a so-called ``archetype'' that defines the overall ``shape'' of the energy convergence.\cite{Olsen_2019} For Hermitian Hamiltonians, these archetypes can be subdivided into five classes (zigzag, interspersed zigzag, triadic, ripples, and geometric), while two additional archetypes (zigzag-geometric and convex-geometric) are observed in non-Hermitian Hamiltonians. %Other features characterising the convergence behaviour of a perturbation method are its rate of convergence, its length of recurring period, and its sign pattern. The geometric archetype appears to be the most common for MP expansions,\cite{Olsen_2019} but the ripples archetype corresponds to some of the early examples of MP convergence. \cite{Handy_1985,Lepetit_1988,Leininger_2000} The three remaining Hermitian archetypes seem to be rarely observed in MP perturbation theory. In contrast, the non-Hermitian coupled cluster perturbation theory,% \cite{Pawlowski_2019a,Pawlowski_2019b,Pawlowski_2019c,Pawlowski_2019d,Pawlowski_2019e} exhibits a range of archetypes including the interspersed zigzag, triadic, ripple, geometric, and zigzag-geometric forms. This analysis highlights the importance of the primary critical point in controlling the high-order convergence, regardless of whether this point is inside or outside the complex unit circle. \cite{Handy_1985,Olsen_2000} %======================================= \subsection{M{\o}ller--Plesset Critical Point} \label{sec:MP_critical_point} %======================================= % STILLINGER INTRODUCES THE CRITICAL POINT In the early 2000's, Stillinger reconsidered the mathematical origin behind the divergent series with odd-even sign alternation.\cite{Stillinger_2000} This type of convergence behaviour corresponds to Cremer and He's class B systems with closely spaced electron pairs and includes \ce{Ne}, \ce{HF}, \ce{F-}, and \ce{H2O}.\cite{Cremer_1996} Stillinger proposed that these series diverge due to a dominant singularity on the negative real $\lambda$ axis, corresponding to a multielectron autoionisation threshold.\cite{Stillinger_2000} To understand Stillinger's argument, consider the parametrised MP Hamiltonian in the form \begin{multline} \label{eq:HamiltonianStillinger} \hH(\lambda) = \sum_{i}^{\Ne} \Bigg[ \overbrace{-\frac{1}{2}\grad_i^2 - \sum_{A}^{\Nn} \frac{Z_A}{\abs{\vb{r}_i-\vb{R}_A}}}^{\text{independent of $\lambda$}} \\ + \underbrace{(1-\lambda)v^{\text{HF}}(\vb{x}_i)}_{\text{repulsive for $\lambda < 1$}} + \underbrace{\lambda\sum_{i