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.
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.
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 which are lines or curves where the function ``jumps'' from one value to another.
Critical points are singularities which lie on the real axis and where the nature of the function undergoes a sudden transition.
However, these do not clearly belong to a given class of singularities and they cannot be rigorously classified as they have more complicated functional forms.
Exact energies for the Hubbard dimer ($U=4t$) as functions of $\lambda$ on the real axis (\subref{subfig:FCI_real}) and in the complex plane (\subref{subfig:FCI_cplx}).
Only the interacting closed-shell singlets are plotted in the complex plane, becoming degenerate at the EP (black dot).
To illustrate the concepts discussed throughout this article, we consider the symmetric Hubbard dimer at half filling, \ie\ with two opposite-spin fermions.
Analytically solvable models are essential in theoretical chemistry and physics as their mathematical simplicity compared to realistic systems (e.g., atoms and molecules) allows new concepts and methods to be
easily tested while retaining the key physical phenomena.
To illustrate the formation of an EP, we scale the off-diagonal coupling strength by introducing the complex parameter $\lambda$ through the transformation $t \to\lambda t$ to give the parameterised Hamiltonian $\hH(\lambda)$.
While the open-shell triplet ($E_{\text{T}}$) and singlet ($E_{\text{S}}$) are independent of $\lambda$, the closed-shell singlet ground state ($E_{-}$) and doubly-excited state ($E_{+}$) couple strongly to form an avoided crossing at $\lambda=0$ (see Fig.~\ref{subfig:FCI_real}).
At non-zero values of $U$ and $t$, these closed-shell singlets can only become degenerate at a pair of complex conjugate points in the complex $\lambda$ plane
Crucially, the energy surface becomes non-analytic at $\lambda_{\text{EP}}$ and a square-root singularity forms with two branch cuts running along the imaginary axis from $\lambda_{\text{EP}}$ to $\pm\i\infty$ (see Fig.~\ref{subfig:FCI_cplx}).
On the real $\lambda$ axis, these EPs lead to the singlet avoided crossing at $\lambda=\Re(\lambda_{\text{EP}})$.
The ``shape'' of this avoided crossing is related to the magnitude of $\Im(\lambda_{\text{EP}})$, with smaller values giving a ``sharper'' interaction.
Remarkably, the existence of these square-root singularities means that following a complex contour around an EP in the complex $\lambda$ plane will interconvert the closed-shell ground and excited states (see Fig.~\ref{subfig:FCI_cplx}).
As a result, completely encircling an EP leads to the interconversion of the two interacting states, while a second complete rotation returns the two states to their original energies.
Additionally, the wave functions pick up a geometric phase in the process, and four complete loops are required to recover their starting forms.\cite{MoiseyevBook}
In the Hartree--Fock (HF) approximation, the many-electron wave function is approximated as a single Slater determinant $\Psi^{\text{HF}}(\vb{x}_1,\ldots,\vb{x}_N)$, where $\vb{x}=(\sigma,\vb{r})$ is a composite vector gathering spin and spatial coordinates.
This Slater determinant is defined as an antisymmetric combination of $\Ne$ (real-valued) occupied one-electron spin-orbitals $\phi_p(\vb{x})$, which are, by definition, eigenfunctions of the one-electron Fock operator
From hereon, $i$ and $j$ denote occupied orbitals, $a$ and $b$ denote unoccupied (or virtual) orbitals, while $p$, $q$, $r$, and $s$ denote arbitrary orbitals.
i) the Epstein-Nesbet (EN) partitioning which consists in taking the diagonal elements of $\hH$ as the zeroth-order Hamiltonian. \cite{Nesbet_1955,Epstein_1926}
ii) the weak correlation partitioning in which the one-electron part is consider as the unperturbed Hamiltonian $\hH^{(0)}$ and the two-electron part is the perturbation operator $\hV$, and
iii) the strong coupling partitioning where the two operators are inverted compared to the weak correlation partitioning. \cite{Seidl_2018}
Convergence of the RMP series as a function of the perturbation order $n$ for the Hubbard dimer at $U/t =3.5$ (where $r_c > 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$.
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$.
%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.
\textit{``In the limit of large order, the series coefficients become equivalent to the Taylor series coefficients of the singularity closest to the origin. ''}
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.
Furthermore, the two-state model of Olsen \etal{}\cite{Olsen_2000} was simply too minimal to understand the complexity of
divergences caused by the MP critical point.
% BASIS SET DEPENDENCE (INCLUDE?)
%Finally, it was shown that $\beta$ singularities are very sensitive to changes in the basis set but not to the bond stretching.
%On the contrary, $\alpha$ singularities are relatively insensitive to the basis sets but very sensitive to bond stretching.
%According to Goodson, \cite{Goodson_2004} the singularity structure of stretched molecules is difficult because there is more than one significant singularity.
%This is consistent with the observation of Olsen and coworkers \cite{Olsen_2000} on the \ce{HF} molecule at equilibrium and stretched geometries.
%To the best of our knowledge, the effect of bond stretching on singularities, its link with spin contamination and symmetry breaking of the wave function has not been as well understood as the ionization phenomenon and its link with diffuse functions.
%In the 2000's, Sergeev and Goodson \cite{Sergeev_2005, Sergeev_2006} analysed this problem from a more mathematical point of view by looking at the whole singularity structure where Olsen and collaborators were trying to find the dominant singularity causing the divergence. \cite{Olsen_1996,Olsen_2000,Olsen_2019}
%Singularities of $\alpha$-type are related to large avoided crossing between the ground and low-lying excited states, whereas $\beta$ singularities come from a sharp avoided crossing between the ground state and a highly diffuse state.
%They succeeded to explain the divergence of the series caused by $\beta$ singularities following previous work by Stillinger. \cite{Stillinger_2000}
%This reasoning is done on the exact Hamiltonian and energy, \ie, the Hamiltonian in the complete Hilbert space, this is the exact energy which exhibits this singularity on the negative real axis.
%However, in a finite basis set which does not span the complete Hilbert space, one can prove that, for a Hermitian Hamiltonian, the singularities of $E(\lambda)$ occurs in complex conjugate pairs with non-zero imaginary parts.
%Sergeev and Goodson proved, \cite{Sergeev_2005} as predicted by Stillinger, \cite{Stillinger_2000} that in a finite basis set the critical point on the real axis is modelled by a cluster of sharp avoided crossings with diffuse functions, equivalently by a cluster of $\beta$ singularities in the negative half plane.
%This explains that Olsen \textit{et al.}, because they used a simple two-state model, only observed the first singularity of this cluster of singularities causing the divergence. \cite{Olsen_2000}
% RELATIONSHIP TO QUANTUM PHASE TRANSITION
When a Hamiltonian is parametrised by a variable such as $\lambda$, the existence of abrupt changes in the
eigenstates as a function of $\lambda$ indicate the presence of a zero-temperature quantum phase transition (QPT).%
Meanwhile, as an avoided crossing becomes increasingly sharp, the corresponding EPs move increasingly close to the real axis, eventually forming a critical point.
%In the previous section, we saw that a careful analysis of the structure of the Hamiltonian allows us to predict the existence of a critical point.
%In a finite basis set, this critical point is model by a cluster of $\beta$ singularities.
%It is now well known that this phenomenon is a special case of a more general phenomenon.
%Indeed, theoretical physicists proved that EPs close to the real axis are connected to \textit{quantum phase transitions} (QPTs). \cite{Heiss_1988,Heiss_2002,Borisov_2015,Sindelka_2017,CarrBook,Vojta_2003,SachdevBook,GilmoreBook}
%In quantum mechanics, the Hamiltonian is almost always dependent of, at least, one parameter.
%In some cases the variation of a parameter can lead to abrupt changes at a critical point.
%These QPTs exist both for ground and excited states as shown by Cejnar and coworkers. \cite{Cejnar_2005,Cejnar_2007,Caprio_2008,Cejnar_2009,Sachdev_2011,Cejnar_2015,Cejnar_2016, Macek_2019,Cejnar_2020}
%A ground-state QPT is characterised by the derivatives of the ground-state energy with respect to a non-thermal control parameter. \cite{Cejnar_2009, Sachdev_2011}
%The transition is called discontinuous and of first order if the first derivative is discontinuous at the critical parameter value.
%Otherwise, it is called continuous and of $m$th order (with $m \ge 2$) if the $m$th derivative is discontinuous.
%A QPT can also be identify by the discontinuity of an appropriate order parameter (or one of its derivatives).
%
%The presence of an EP close to the real axis is characteristic of a sharp avoided crossing.
%Yet, at such an avoided crossing, eigenstates change abruptly.
%Although it is now well understood that EPs are closely related to QPTs, the link between the type of QPT (ground state or excited state, first or higher order) and EPs still need to be clarified.
%One of the major obstacles that one faces in order to achieve this resides in the ability to compute the distribution of EPs.
%The numerical assignment of an EP to two energies on the real axis is very difficult in large dimensions.
%Hence, the design of specific methods are required to get information on the location of EPs.
%Following this idea, Cejnar \textit{et al.}~developed a method based on a Coulomb analogy giving access to the density of EP close to the real axis. \cite{Cejnar_2005, Cejnar_2007}
%More recently Stransky and coworkers proved that the distribution of EPs is characteristic of the QPT order. \cite{Stransky_2018}
%In the thermodynamic limit, some of the EPs converge towards a critical point $\lambda_\text{c}$ on the real axis.
%They showed that, within the interacting boson model, \cite{Lipkin_1965} EPs associated to first- and second-order QPT behave differently when the number of particles increases.
%The position of these singularities converge towards the critical point on the real axis at different rates (exponentially and algebraically for the first and second orders, respectively) with respect to the number of particles.
%
%Moreover, Cejnar \textit{et al.}~studied the so-called shape-phase transitions of the interaction boson model from a QPT point of view. \cite{Cejnar_2000, Cejnar_2003, Cejnar_2007a, Cejnar_2009}
%The phase of the ensemble of $s$ and $d$ bosons is characterised by a dynamical symmetry.
%When a parameter is continuously modified the dynamical symmetry of the system can change at a critical value of this parameter, leading to a deformed phase.
%They showed that at this critical value of the parameter, the system undergoes a QPT.
%For example, without interaction the ground state is the spherical phase (a condensate of $s$ bosons) and when the interaction increases it leads to a deformed phase constituted of a mixture of $s$ and $d$ bosons states.
%In particular, we see that the transition from the spherical phase to the axially symmetric one is analog to the symmetry breaking of the wave function of the hydrogen molecule when the bond is stretched. \cite{SzaboBook}
%It seems like our understanding of the physics of spatial and/or spin symmetry breaking in HF theory can be enlightened by QPT theory.
%Indeed, the second derivative of the HF ground-state energy is discontinuous at the point of spin symmetry-breaking which means that the system undergo a second-order QPT.
%
%Moreover, the $\beta$ singularities introduced by Sergeev and coworkers to describe the EPs modelling the formation of a bound cluster of electrons are actually a more general class of singularities.
%The EPs close to the real axis (the so-called $\beta$ singularities) are connected to QPT because they result from a sharp avoided crossings at which the eigenstates change quickly.
%However, the $\alpha$ singularities arise from large avoided crossings.
%Thus, they cannot be connected to QPT.
%The avoided crossings generating $\alpha$ singularities generally involve the ground state and low-lying doubly-excited states.
%Those excited states have a non-negligible contribution to the exact FCI solution because they have (usually) the same spatial and spin symmetry as the ground state.
%We believe that $\alpha$ singularities are connected to states with non-negligible contribution in the CI expansion thus to the dynamical part of the correlation energy, while $\beta$ singularities are linked to symmetry breaking and phase transitions of the wave function, \ie, to the multi-reference nature of the wave function thus to the static part of the correlation energy.
The inability of Taylor series to model properly the energy function $E(\lambda$) can be simply understood by the fact that one aims at modelling a complicated function with potentially poles and singularities by a simple polynomial of finite order.
Nonetheless, the description of complex energy functions can be significantly improved thanks to Pad\'e approximant, \cite{Pade_1892} and related techniques. \cite{BakerBook,BenderBook}
(with $b_0=1$), where the coefficients of the polynomials $A(\lambda)$ and $B(\lambda)$ are determined by collecting terms according to power of $\lambda$.
Pad\'e approximants are extremely useful in many areas of physics and chemistry \cite{Loos_2013,Pavlyukh_2017,Tarantino_2019,Gluzman_2020} as they can model poles, which appears at the locations of the roots of $B(\lambda)$.
However, they are unable to model functions with square-root branch points (which are ubiquitous in the singularity structure of a typical perturbative treatment) and more complicated functional forms appearing at critical points (where the nature of the solution undergoes a sudden transition) for example.
Figure \ref{fig:PadeRMP} illustrates the improvement brought by diagonal (\ie, $d_A = d_B$) Pad\'e approximants as compared to the usual Taylor expansion in cases where the RMP series of the Hubbard dimer converges ($U/t =3.5$) and diverges ($U/t =4.5$).
More quantitatively, Table \ref{tab:PadeRMP} gathers estimates of the RMP ground-state energy at $\lambda=1$ provided by various truncated Taylor series and Pad\'e approximants for these two values of the ratio $U/t$.
While the truncated Taylor series converges laboriously to the exact energy at $U/t =3.5$ when one increases the truncation degree, the Pad\'e approximants yield much more accurate results with, additionally, a rather good estimate of the radius of convergence of the RMP series.
For $U/t =4.5$, the struggles of the truncated Taylor expansions are magnified and the Pad\'e approximants still provide quite accurate energies even outside the radius of convergence of the RMP series.
In a nutshell, the idea behind quadratic approximant is to model the singularity structure of the energy function $E(\lambda)$ via a generalised version of the square-root singularity expression \cite{Mayer_1985,Goodson_2011,Goodson_2019}
and substituting $E(\lambda$) by its $n$th-order expansion and the polynomials by their respective expressions \eqref{eq:PQR} yields $n+1$ linear equations for the coefficients $p_k$, $q_k$, and $r_k$ (where we are free to assume that $q_0=1$).
A quadratic approximant, characterised by the label $[d_P/d_Q,d_R]$, generates, by construction, $n_\text{bp}=\max(2d_p,d_q+d_r)$ branch points at the roots of the polynomial $P^2(\lambda)-4 Q(\lambda) R(\lambda)$.
Note that, by construction, a quadratic approximant has only two branches which hampers the faithful description of more complicated singularity structures.
As shown in Ref.~\onlinecite{Goodson_2000a}, quadratic approximants provide convergent results in the most divergent cases considered by Olsen and collaborators \cite{Christiansen_1996,Olsen_1996} and Leininger \etal\cite{Leininger_2000}
For the RMP series of the Hubbard dimer, the $[0/0,0]$ and $[1/0,0]$ quadratic approximant are quite poor approximation, but its $[1/0,1]$ version already model perfectly the RMP energy function by predicting a single pair of EPs at $\lambda_\text{EP}=\pm i 4t/U$.
We can anticipate that the singularity structure of the UMP energy function is going to be much more challenging to model properly, and this is indeed the case as the UMP energy function contains three branches (see Figs.~\ref{subfig:UMP_3} and \ref{subfig:UMP_7}).
However, by ramping up high enough the degree of the polynomials, one is able to get both, as shown in Fig.~\ref{fig:QuadUMP} and Table \ref{tab:QuadUMP}, accurate estimates of the radius of convergence of the UMP series and of the ground-state energy at $\lambda=1$, even in cases where the convergence of the UMP series is painfully slow (see Fig.~\ref{subfig:UMP_cvg}).
Figure \ref{fig:QuadUMP} evidences that the Pad\'e approximants are trying to model the square root singularity by placing a pole on the real axis (for [3/3]) or just off the real axis (for [4/4]).
Thanks to greater flexibility, the quadratic approximants are able to model nicely the avoided crossing and the location of the singularities.
Besides, they provide accurate estimates of the ground-state energy at $\lambda=1$ (see Table \ref{tab:QuadUMP}).
\caption{Estimate of the radius of convergence $r_c$ of the UMP energy function provided by various resummation techniques at $U/t =3$ and $7$.
The truncation degree of the Taylor expansion $n$ of $E(\lambda)$ and the number of branch points $n_\text{bp}=\max(2d_p,d_q+d_r)$ generated by the quadratic approximants are also reported.
An interesting point raised in Ref.~\onlinecite{Goodson_2019} suggests that low-order quadratic approximants might struggle to model the correct singularity structure when the energy function has poles in both the positive and negative half-planes.
In such a scenario, the quadratic approximant will have the tendency to place its branch points in-between, potentially introducing singularities quite close to the origin.
A simple potential cure for this consists in applying a judicious transformation (like a bilinear conformal mapping) which does not affect the points at $\lambda=0$ and $\lambda=1$. \cite{Feenberg_1956}
Recently, Mih\'alka \textit{et al.} studied the partitioning effect on the convergence properties of Rayleigh-Schr\"odinger perturbation theory by considering the MP and the EN partitioning as well as an alternative partitioning \cite{Mihalka_2017a} (see also Ref.~\onlinecite{Surjan_2000}).
Taking as an example (in particular) the water molecule at equilibrium and at stretched geometries, they could estimate the radius of convergence via a quadratic Pad\'e approximant and convert divergent perturbation expansions to convergent ones in some cases thanks to a judicious choice of the level shift parameter.
In a subsequent study by the same group, \cite{Mihalka_2017b} they use analytic continuation techniques to resum divergent MP series \cite{Goodson_2011} taking again as an example the water molecule in a stretched geometry.
In a nutshell, their idea consists in calculating the energy of the system for several values of $\lambda$ for which the MP series is rapidly convergent (\ie, for $\abs{\lambda} < r_c$), and to extrapolate the final energy to the physical system at $\lambda=1$ via a polynomial- or Pad\'e-based fit.
However, the choice of the functional form of the fit remains a subtle task.
This technique was first generalised by using complex scaling parameters and applying analytic continuation by solving the Laplace equation, \cite{Surjan_2018} and then further improved thanks to Cauchy's integral formula \cite{Mihalka_2019}
which states that the value of the energy can be computed at $\lambda=a$ inside the complex contour $\gamma$ only by the knowledge of its values on the same contour.
Their method consists in refining self-consistently the values of $E(\lambda)$ computed on a contour going through the physical point at $\lambda=1$ and encloses points of the ``trusted'' region (where the MP series is convergent). The shape of this contour is arbitrary but no singularities are allowed inside the contour to ensure $E(\lambda)$ is analytic.
When the values of $E(\lambda)$ on the so-called contour are converged, Cauchy's integrals formula \eqref{eq:Cauchy} is invoked to compute the values at $E(\lambda=1)$ which corresponds to the final estimate of the FCI energy.
The authors illustrate this protocol on the dissociation curve of \ce{LiH} and the stretched water molecule showing encouraging results. \cite{Mihalka_2019}
In order to model accurately chemical systems, one must choose, in a ever growing zoo of methods, which computational protocol is adapted to the system of interest.
This choice can be, moreover, motivated by the type of properties that one is interested in.
That means that one must understand the strengths and weaknesses of each method, \ie, why one method might fail in some cases and work beautifully in others.
We have seen that for methods relying on perturbation theory, their successes and failures are directly connected to the position of EPs in the complex plane.
Exhaustive studies have been performed on the causes of failure of MP perturbation theory.
First, it was understood that, for chemical systems for which the HF Slater determinant is a poor approximation to the exact wave function, MP perturbation theory fails too.
Such systems can be, for example, molecules where the exact ground-state wave function is dominated by more than one configuration, \ie, multi-reference systems.
For instance, it has been shown that systems considered as well understood (\eg, \ce{Ne}) can exhibit divergent behaviour when the basis set is augmented with diffuse functions.
Later, these erratic behaviours of the perturbation series were investigated and rationalised in terms of avoided crossings and singularities in the complex plane.
It was shown that the singularities can be classified in two families.
The first family includes $\alpha$ singularities resulting from a large avoided crossing between the ground state and a low-lying doubly-excited states.
The $\beta$ singularities, which constitutes the second family, are artefacts generated by the incompleteness of the Hilbert space, and they are directly connected to an ionisation phenomenon occurring in the complete Hilbert space.
We have found that the $\beta$ singularities modelling the ionisation phenomenon described by Sergeev and Goodson are actually part of a more general class of singularities.
Indeed, those singularities close to the real axis are connected to quantum phase transitions and symmetry breaking, and theoretical physics have demonstrated that the behaviour of the EPs depends of the type of transitions from which the EPs result (first or higher orders, ground state or excited state transitions).
To conclude, this work shows that our understanding of the singularity structure of the energy is still incomplete but we hope that it opens new perspectives for the understanding of the physics of EPs in electronic structure theory.
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).