Perturbation theory in the complex plane: exceptional points and where to find them

Within the Born–Oppenheimer approximation, the exact molecular Hamiltonian with N electrons and M (clamped) nuclei is defined for a given nuclear framework as

$$\begin{equation}\hat{H}\left(\mathbf{R}\right)=-\frac{1}{2}\sum\limits _{i}^{N}{\mathbf{\nabla }}_{i}^{2}-\sum\limits _{i}^{N}\sum\limits _{A}^{M}\frac{{Z}_{A}}{\left\vert {\mathbf{r}}_{i}-{\mathbf{R}}_{A}\right\vert }+\sum\limits _{i{< }j}^{N}\frac{1}{\left\vert {\mathbf{r}}_{i}-{\mathbf{r}}_{j}\right\vert },\end{equation} \tag{ 1 }$$

where r_i defines the position of the ith electron, R_A and Z_A are the position and charge of the Ath nucleus respectively, and R = (R₁, ..., R_M ) is a collective vector for the nuclear positions. The first term represents the kinetic energy of the electrons, while the two following terms account for the electron–nucleus attraction and the electron–electron repulsion.

The exact many-electron wave function at a given nuclear geometry Ψ(R) corresponds to the solution of the (time-independent) Schrödinger equation

$$\begin{equation}\hat{H}\left(\mathbf{R}\right)\enspace {\Psi}\left(\mathbf{R}\right)=E\left(\mathbf{R}\right)\enspace {\Psi}\left(\mathbf{R}\right),\end{equation} \tag{ 2 }$$

with the eigenvalues E(R) providing the exact energies. The energy E(R) can be considered as a 'one-to-many' function since each input nuclear geometry yields several eigenvalues corresponding to the ground and excited states of the exact spectrum. However, exact solutions to equation (2) are only possible in the simplest of systems, such as the one-electron hydrogen atom and some specific two-electron systems with well-defined mathematical properties [63–66]. In practice, approximations to the exact Schrödinger equation must be introduced, including perturbation theories and the Hartree–Fock approximation considered in this review. In what follows, we will drop the parametric dependence on the nuclear geometry and, unless otherwise stated, atomic units will be used throughout.

To illustrate the concepts discussed throughout this article, we consider the symmetric Hubbard dimer at half filling, i.e., 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.

Using the (localised) site basis, the Hilbert space of the Hubbard dimer comprises the four configurations

$$\begin{equation*}\left\vert {\mathcal{L}}^{{\uparrow}}{\mathcal{L}}^{{\downarrow}}\right\rangle ,\left\vert {\mathcal{L}}^{{\uparrow}}{\mathcal{R}}^{{\downarrow}}\right\rangle ,\left\vert {\mathcal{R}}^{{\uparrow}}{\mathcal{L}}^{{\downarrow}}\right\rangle ,\left\vert {\mathcal{R}}^{{\uparrow}}{\mathcal{R}}^{{\downarrow}}\right\rangle ,\end{equation*}$$

where ${\mathcal{L}}^{\sigma }$ (${\mathcal{R}}^{\sigma }$) denotes an electron with spin σ on the left (right) site. The exact, or full configuration interaction (FCI), Hamiltonian is then

$$\begin{equation}\mathbf{H}=\left(\begin{matrix}\hfill U\hfill & \hfill -t\hfill & \hfill -t\hfill & \hfill 0\hfill \\ \hfill -t\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -t\hfill \\ \hfill -t\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -t\hfill \\ \hfill 0\hfill & \hfill -t\hfill & \hfill -t\hfill & \hfill U\hfill \end{matrix}\right),\end{equation} \tag{ 3 }$$

where t is the hopping parameter and U is the on-site Coulomb repulsion. We refer the interested reader to references [67, 68] for more details about this system. The parameter U controls the strength of the electron correlation. In the weak correlation regime (small U), the kinetic energy dominates and the electrons are delocalised over both sites. In the large-U (or strong correlation) regime, the electron repulsion term becomes dominant and the electrons localise on opposite sites to minimise their Coulomb repulsion. This phenomenon is often referred to as Wigner crystallisation [11].

To illustrate the formation of an EP, we scale the off-diagonal coupling strength by introducing the complex parameter λ through the transformation t → λt to give the parameterised Hamiltonian $\hat{H}\left(\lambda \right)$. When λ is real, the Hamiltonian is Hermitian with the distinct (real-valued) (eigen)energies

$$\begin{equation}{E}_{\mp }=\frac{1}{2}\left(U\mp \sqrt{{\left(4\lambda t\right)}^{2}+{U}^{2}}\right),\end{equation} \tag{ 4a }$$

$$\begin{equation}{E}_{\text{T}}=0,\end{equation} \tag{ 4b }$$

$$\begin{equation}{E}_{\text{S}}=U.\end{equation} \tag{ 4c }$$

While the open-shell triplet (E_T) and singlet (E_S) are independent of λ, the closed-shell singlet ground state (E₋) and doubly-excited state (E₊) couple strongly to form an avoided crossing at λ = 0 (see figure 1(a)). In contrast, when λ is complex, the energies may become complex-valued, with the real components shown in figure 1(b). Although the imaginary component of the energy is linked to resonance lifetimes elsewhere in non-Hermitian quantum mechanics [24], its physical interpretation in the current context is unclear. Throughout this work, we will generally consider and plot only the real component of any complex-valued energies.

Zoom In Zoom Out Reset image size

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 λ plane

$$\begin{equation}{\lambda }_{\text{EP}}={\pm}\mathrm{i}\frac{U}{4t},\end{equation} \tag{ 5 }$$

with energy

$$\begin{equation}{E}_{\text{EP}}=\frac{U}{2}.\end{equation} \tag{ 6 }$$

These λ values correspond to so-called EPs and connect the ground and excited states in the complex plane. Crucially, the ground- and excited-state wave functions at an EP become identical rather than just degenerate. Furthermore, the energy surface becomes non-analytic at λ_EP and a square-root singularity forms with two branch cuts running along the imaginary axis from λ_EP to ±i∞ (see figure 1(b)). Along these branch cuts, the real components of the energies are equivalent and appear to give a seam of intersection, but a strict degeneracy is avoided because the imaginary components are different.

On the real λ axis, these EPs lead to the singlet avoided crossing at $\lambda =\mathrm{R}\mathrm{e}\left({\lambda }_{\text{EP}}\right)$. The 'shape' of this avoided crossing is related to the magnitude of $\mathrm{I}\mathrm{m}\left({\lambda }_{\text{EP}}\right)$, with smaller values giving a 'sharper' interaction. In the limit U/t → 0, the two EPs converge at λ_EP = 0 to create a conical intersection with a gradient discontinuity on the real axis. This gradient discontinuity defines a critical point in the ground-state energy, where a sudden change occurs in the electronic wave function, and can be considered as a zero-temperature quantum phase transition (QPT) [29, 69–75].

Remarkably, the existence of these square-root singularities means that following a complex contour around an EP in the complex λ plane will interconvert the closed-shell ground and excited states (see figure 1(b)). This behaviour can be seen by expanding the radicand in equation (4a) as a Taylor series around λ_EP to give

$$\begin{equation}{E}_{{\pm}}\approx {E}_{\text{EP}}{\pm}\sqrt{32{t}^{2}{\lambda }_{\text{EP}}}\sqrt{\lambda -{\lambda }_{\text{EP}}}.\end{equation} \tag{ 7 }$$

Parametrising the complex contour as λ(θ) = λ_EP + R exp(iθ) gives the continuous energy pathways

$$\begin{equation}{E}_{{\pm}}\left(\theta \right)\approx {E}_{\text{EP}}{\pm}\sqrt{32{t}^{2}{\lambda }_{\text{EP}}R}\enspace \mathrm{exp}\left(\mathrm{i}\theta /2\right)\end{equation} \tag{ 8 }$$

such that E_±(2π) = E_∓(0) and E_±(4π) = E_±(0). 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 can pick up a geometric phase in the process, and four complete loops are required to recover their starting forms [24].

To locate EPs in practice, one must simultaneously solve

$$\begin{equation}\mathrm{det}\left[\hat{H}\left(\lambda \right)-E\hat{I}\right]=0,\end{equation} \tag{ 9a }$$

$$\begin{equation}\frac{\partial }{\partial E}\enspace \mathrm{det}\left[\hat{H}\left(\lambda \right)-E\hat{I}\right]=0,\end{equation} \tag{ 9b }$$

where $\hat{I}$ is the identity operator [76]. Equation (9a) is the well-known secular equation providing the (eigen)energies of the system. If the energy is also a solution of equation (9b), then this energy value is at least two-fold degenerate. These degeneracies can be conical intersections between two states with different symmetries for real values of λ [35], or EPs between two states with the same symmetry for complex values of λ.

One of the most common routes to approximately solving the Schrödinger equation is to introduce a perturbative expansion of the exact energy. Within Rayleigh–Schrödinger perturbation theory, the time-independent Schrödinger equation is recast as

$$\begin{equation}\hat{H}\left(\lambda \right){\Psi}\left(\lambda \right)=\left({\hat{H}}^{\left(0\right)}+\lambda \hat{V}\right){\Psi}\left(\lambda \right)=E\left(\lambda \right){\Psi}\left(\lambda \right),\end{equation} \tag{ 10 }$$

where ${\hat{H}}^{\left(0\right)}$ is a zeroth-order Hamiltonian and $\hat{V}=\hat{H}-{\hat{H}}^{\left(0\right)}$ represents the perturbation operator. Expanding the wave function and energy as power series in λ as

$$\begin{equation}{\Psi}\left(\lambda \right)=\sum\limits _{k=0}^{\infty }{\lambda }^{k}\enspace {{\Psi}}^{\left(k\right)},\end{equation} \tag{ 11a }$$

$$\begin{equation}E\left(\lambda \right)=\sum\limits _{k=0}^{\infty }{\lambda }^{k}\enspace {E}^{\left(k\right)},\end{equation} \tag{ 11b }$$

solving the corresponding perturbation equations up to a given order n, and setting λ = 1 then yields approximate solutions to equation (2).

Mathematically, equation (11b) corresponds to a Taylor series expansion of the exact energy around the reference system λ = 0. The energy of the target 'physical' system is recovered at the point λ = 1. However, like all series expansions, equation (11b) has a radius of convergence r_c. When r_c < 1, the Rayleigh–Schrödinger expansion will diverge for the physical system. The value of r_c can vary significantly between different systems and strongly depends on the particular decomposition of the reference and perturbation Hamiltonians in equation (10) [61]. From complex analysis [56], the radius of convergence for the energy can be obtained by looking for the non-analytic singularities of E(λ) in the complex λ plane. This property arises from the following theorem [77]:

'The Taylor series about a point z₀ of a function over the complex z plane will converge at a value z₁ if the function is non-singular at all values of z in the circular region centred at z₀ with radius $\left\vert {z}_{1}-{z}_{0}\right\vert$. If the function has a singular point z_s such that $\left\vert {z}_{\mathrm{s}}-{z}_{0}\right\vert {< }\left\vert {z}_{1}-{z}_{0}\right\vert$, then the series will diverge when evaluated at z₁.'

As a result, the radius of convergence for a function is equal to the distance from the origin of the closest singularity in the complex plane, referred to as the 'dominant' singularity. This singularity may represent a pole of the function, or a branch point (e.g., square-root or logarithmic) in a multi-valued function.

For example, the simple function

$$\begin{equation}f\hspace{-1pt}\left(x\right)=\frac{1}{1+{x}^{4}}.\end{equation} \tag{ 12 }$$

is smooth and infinitely differentiable for $x\in \mathbb{R}$, and one might expect that its Taylor series expansion would converge in this domain. However, this series diverges for x ⩾ 1. This divergence occurs because f(x) has four poles in the complex plane (e^iπ/4, e^−iπ/4, e^i3π/4, and e^−i3π/4) with a modulus equal to 1, demonstrating that complex singularities are essential to fully understand the series convergence on the real axis [56].

The radius of convergence for the perturbation series equation (11b) is therefore dictated by the magnitude ${r}_{\mathrm{c}}=\left\vert {\lambda }_{\mathrm{c}}\right\vert$ of the singularity in E(λ) that is closest to the origin. Note that when $\left\vert \lambda \right\vert ={r}_{\mathrm{c}}$, one cannot a priori predict if the series is convergent or not. For example, the series ${\sum }_{k=1}^{\infty }{\lambda }^{k}/k$ diverges at λ = 1 but converges at λ = −1.

Like the exact system in section 2.2, the perturbation energy E(λ) represents a 'one-to-many' function with the output elements representing an approximation to both the ground and excited states. The most common singularities on E(λ) therefore correspond to non-analytic EPs in the complex λ plane where two states become degenerate. Additional singularities can also arise at critical points of the energy. A critical point corresponds to the intersection of two energy surfaces where the eigenstates remain distinct but a gradient discontinuity occurs in the ground-state energy. In contrast, at a square-root branch point, both the energies and the associated wave functions of the intersecting surfaces become identical. Later we will demonstrate how the choice of reference Hamiltonian controls the position of these EPs, and ultimately determines the convergence properties of the perturbation series.

In the HF approximation, the many-electron wave function is approximated as a single Slater determinant Ψ_HF(x₁, ..., x_N ), where x = (σ, r) is a composite vector gathering spin and spatial coordinates. This Slater determinant is defined as an antisymmetric combination of N (real-valued) occupied one-electron spin–orbitals ϕ_p (x), which are, by definition, eigenfunctions of the one-electron Fock operator

$$\begin{equation}\hat{f}\left(\mathbf{x}\right){\phi }_{p}\left(\mathbf{x}\right)=\left[\hat{h}\left(\mathbf{x}\right)+{\hat{v}}_{\text{HF}}\left(\mathbf{x}\right)\right]{\phi }_{p}\left(\mathbf{x}\right)={{\epsilon}}_{p}{\phi }_{p}\left(\mathbf{x}\right).\end{equation} \tag{ 13 }$$

Here the (one-electron) core Hamiltonian is

$$\begin{equation}\hat{h}\left(\mathbf{x}\right)=-\frac{{\mathbf{\nabla }}^{2}}{2}+\sum\limits _{A}^{M}\frac{{Z}_{A}}{\left\vert \mathbf{r}-{\mathbf{R}}_{A}\right\vert }\end{equation} \tag{ 14 }$$

and

$$\begin{equation}{\hat{v}}_{\text{HF}}\left(\mathbf{x}\right)=\sum\limits _{i}^{N}\left[{\hat{J}}_{i}\left(\mathbf{x}\right)-{\hat{K}}_{i}\left(\mathbf{x}\right)\right]\end{equation} \tag{ 15 }$$

is the HF mean-field electron–electron potential with

$$\begin{equation} {\hat{J}}_{i}\left(\mathbf{x}\right){\phi }_{j}\left(\mathbf{x}\right)=\left(\int {\phi }_{i}\left({\mathbf{x}}^{\prime }\right)\frac{1}{\left\vert \mathbf{r}-{\mathbf{r}}^{\prime }\right\vert }{\phi }_{i}\left({\mathbf{x}}^{\prime }\right)d{\mathbf{x}}^{\prime }\right){\phi }_{j}\left(\mathbf{x}\right),\end{equation} \tag{ 16a }$$

$$\begin{equation} {\hat{K}}_{i}\left(\mathbf{x}\right){\phi }_{j}\left(\mathbf{x}\right)=\left(\int {\phi }_{i}\left({\mathbf{x}}^{\prime }\right)\frac{1}{\left\vert \mathbf{r}-{\mathbf{r}}^{\prime }\right\vert }{\phi }_{j}\left({\mathbf{x}}^{\prime }\right)d{\mathbf{x}}^{\prime }\right){\phi }_{i}\left(\mathbf{x}\right),\end{equation} \tag{ 16b }$$

defining the Coulomb and exchange operators (respectively) in the spin–orbital basis [3]. The HF energy is then defined as

$$\begin{equation}{E}_{\text{HF}}=\frac{1}{2}\sum\limits _{i}^{N}\left({h}_{i}+{f}_{i}\right),\end{equation} \tag{ 17 }$$

with the corresponding matrix elements

$$\begin{equation}{h}_{i}=\left\langle {\phi }_{i}\right\vert \hat{h}\left\vert {\phi }_{i}\right\rangle ,\qquad {f}_{i}=\left\langle {\phi }_{i}\right\vert \hat{f}\left\vert {\phi }_{i}\right\rangle .\end{equation} \tag{ 18 }$$

The optimal HF wave function is identified by using the variational principle to minimise the HF energy. For any system with more than one electron, the resulting Slater determinant is not an eigenfunction of the exact Hamiltonian $\hat{H}$. However, it is by definition an eigenfunction of the approximate many-electron HF Hamiltonian constructed from the one-electron Fock operators as

$$\begin{equation}{\hat{H}}_{\text{HF}}=\sum\limits _{i}^{N}f\hspace{-1pt}\left({\mathbf{x}}_{i}\right).\end{equation} \tag{ 19 }$$

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.

In the most flexible variant of real HF theory (generalised HF) the one-electron orbitals can be complex-valued and contain a mixture of spin-up and spin-down components [78, 79]. However, the application of HF theory with some level of constraint on the orbital structure is far more common. Forcing the spatial part of the orbitals to be the same for spin-up and spin-down electrons leads to restricted HF (RHF) method, while allowing different orbitals for different spins leads to the so-called unrestricted HF (UHF) approach. [80]. The advantage of the UHF approximation is its ability to correctly describe strongly correlated systems, such as antiferromagnetic phases [81] or the dissociation of the hydrogen dimer [82]. However, by allowing different orbitals for different spins, the UHF wave function is no longer required to be an eigenfunction of the total spin operator ${\hat{\mathcal{S}}}^{2}$, leading to 'spin-contamination'.

In the Hubbard dimer, the HF energy can be parametrised using two rotation angles θ^α and θ^β as

$$\begin{align}{E}_{\text{HF}}\left({\theta }^{\enspace \alpha },{\theta }^{\enspace \beta }\right)& =-t\enspace \left(\mathrm{sin}\enspace {\theta }^{\enspace \alpha }+\mathrm{sin}\enspace {\theta }^{\enspace \beta }\right)\\ & \quad +\frac{U}{2}\left(1+\mathrm{cos}\enspace {\theta }^{\enspace \alpha }\enspace \mathrm{cos}\enspace {\theta }^{\enspace \beta }\right),\end{align} \tag{ 20 }$$

where we have introduced occupied ${\phi }_{1}^{\sigma }$ and unoccupied ${\phi }_{2}^{\sigma }$ molecular orbitals for the spin-σ electrons as

$$\begin{equation}{\phi }_{1}^{\sigma }=\mathrm{cos}\left(\frac{{\theta }^{\sigma }}{2}\right){\mathcal{L}}^{\sigma }+\mathrm{sin}\left(\frac{{\theta }^{\sigma }}{2}\right){\mathcal{R}}^{\sigma },\end{equation} \tag{ 21a }$$

$$\begin{equation}{\phi }_{2}^{\sigma }=-\mathrm{sin}\left(\frac{{\theta }^{\sigma }}{2}\right){\mathcal{L}}^{\sigma }+\mathrm{cos}\left(\frac{{\theta }^{\sigma }}{2}\right){\mathcal{R}}^{\sigma }\end{equation} \tag{ 21b }$$

Equations (20), (21a) and (21b) are valid for both RHF and UHF. In the weak correlation regime 0 ⩽ U ⩽ 2t, the angles which minimise the HF energy, i.e., ∂E_HF/∂θ^σ = 0, are

$$\begin{equation}{\theta }_{\text{RHF}}^{\enspace \alpha }={\theta }_{\text{RHF}}^{\enspace \beta }=\pi /2,\end{equation} \tag{ 22 }$$

giving the molecular orbitals

$$\begin{equation}{\phi }_{1,\text{RHF}}^{\sigma }=\frac{{\mathcal{L}}^{\sigma }+{\mathcal{R}}^{\sigma }}{\sqrt{2}},\qquad {\phi }_{2,\text{RHF}}^{\sigma }=\frac{{\mathcal{L}}^{\sigma }-{\mathcal{R}}^{\sigma }}{\sqrt{2}},\end{equation} \tag{ 23 }$$

and the ground-state RHF energy (figure 2)

$$\begin{equation}{E}_{\text{RHF}}\equiv {E}_{\text{HF}}\left({\theta }_{\text{RHF}}^{\enspace \alpha },{\theta }_{\text{RHF}}^{\enspace \beta }\right)=-2t+\frac{U}{2}.\end{equation} \tag{ 24 }$$

Here, the molecular orbitals respectively transform according to the ${{\Sigma}}_{\text{g}}^{+}$ and ${{\Sigma}}_{\text{u}}^{+}$ irreducible representations of the D_∞h point group that represents the symmetric Hubbard dimer. We can therefore consider these as symmetry-pure molecular orbitals. However, in the strongly correlated regime U > 2t, the closed-shell orbital restriction prevents RHF from modelling the correct physics with the two electrons on opposite sites.

Zoom In Zoom Out Reset image size

As the on-site repulsion is increased from 0, the HF approximation reaches a critical value at U = 2t where an alternative 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 ⩾ 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 [82]. For U ⩾ 2t, the optimal orbital rotation angles for the UHF orbitals become

$$\begin{equation}{\theta }_{\text{UHF}}^{\enspace \alpha }=\mathrm{arctan}\left(-\frac{2t}{\sqrt{{U}^{2}-4{t}^{2}}}\right),\end{equation} \tag{ 25a }$$

$$\begin{equation}{\theta }_{\text{UHF}}^{\enspace \beta }=\mathrm{arctan}\left(+\frac{2t}{\sqrt{{U}^{2}-4{t}^{2}}}\right),\end{equation} \tag{ 25b }$$

with the corresponding UHF ground-state energy (figure 2)

$$\begin{equation}{E}_{\text{UHF}}\equiv {E}_{\text{HF}}\left({\theta }_{\text{UHF}}^{\enspace \alpha },{\theta }_{\text{UHF}}^{\enspace \beta }\right)=-\frac{2{t}^{2}}{U}.\end{equation} \tag{ 26 }$$

The molecular orbitals of the lower-energy UHF solution do not transform as an irreducible representation of the D_∞h point group and therefore break spatial symmetry. Allowing different orbitals for the different spins also means that the overall wave function is no longer an eigenfunction of the ${\mathcal{S}}^{2}$ operator and can be considered to break spin symmetry. This combined spatial and spin symmetry-breaking occurs for all U ⩾ 2t. Furthermore, time-reversal symmetry dictates that this UHF wave function must be degenerate with its spin-flipped counterpart, obtained by swapping ${\theta }_{\text{UHF}}^{\enspace \alpha }$ and ${\theta }_{\text{UHF}}^{\enspace \beta }$ in equations (25a) and (25b). 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 [83]. Spatial 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 [80, 84].

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 [85, 86]. Alternatively, the non-linear terms arising from the Coulomb and exchange operators can be considered as a perturbation from the core Hamiltonian (14) by introducing the transformation U → λ U, giving the parametrised Fock operator

$$\begin{equation}\hat{f}\left(\mathbf{x};\lambda \right)=\hat{h}\left(\mathbf{x}\right)+\lambda \enspace {\hat{v}}_{\text{HF}}\left(\mathbf{x}\right).\end{equation} \tag{ 27 }$$

The orbitals in the reference problem λ = 0 correspond to the symmetry-pure eigenfunctions of the one-electron core Hamiltonian, while self-consistent solutions at λ = 1 represent the orbitals of the true HF solution.

For real λ, the self-consistent HF energies at given (real) U and t values in the Hubbard dimer directly mirror the energies shown in figure 2, with coalesence points at

$$\begin{equation}{\lambda }_{\text{c}}={\pm}\frac{2t}{U}.\end{equation} \tag{ 28 }$$

In contrast, when λ becomes complex, the HF equations become non-Hermitian and each HF solution can be analytically continued for all λ values using the holomorphic HF approach [87–89]. Remarkably, the coalescence point in this analytic continuation emerges as a quasi-EP on the real λ axis (figure 3), where the different HF solutions become equivalent but not self-orthogonal [36]. By analogy with perturbation theory, the regime where this quasi-EP occurs within λ_c ⩽ 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 λ = 1. For example, in the Hubbard dimer with U > 2t, one finds λ_c < 1 and the symmetry-pure orbitals do not provide a good representation of the HF ground state. In contrast, U < 2t yields λ_c > 1 and corresponds to the regime where the HF ground state is correctly represented by symmetry-pure orbitals.

Zoom In Zoom Out Reset image size

We have recently shown that the complex scaled Fock operator (28) also allows states of different symmetries to be interconverted by following a well-defined contour in the complex λ-plane [36]. In particular, by slowly varying λ in a similar (yet different) manner to an adiabatic connection in density-functional theory (DFT) [90–92], 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 quasi-EPs in determining the behaviour of the HF approximation. Furthermore, the complex-scaled Fock operator can be used routinely to construct analytic continuations of HF solutions beyond the points where real HF solutions coalesce and vanish [93].

In electronic structure, the HF Hamiltonian (19) is often used as the zeroth-order Hamiltonian to define Møller–Plesset (MP) perturbation theory [10]. This approach can recover a large proportion of the electron correlation energy [94–96], and provides the foundation for numerous post-HF approximations. With the MP partitioning, the parametrised perturbation Hamiltonian becomes

$$\begin{align}\hat{H}\left(\lambda \right)& =\sum\limits _{i}^{N}\left[-\frac{{\mathbf{\nabla }}_{i}^{2}}{2}-\sum\limits _{A}^{M}\frac{{Z}_{A}}{\left\vert {\mathbf{r}}_{i}-{\mathbf{R}}_{A}\right\vert }\right]\\ & \quad +\left(1-\lambda \right)\sum\limits _{i}^{N}{\hat{v}}_{\text{HF}}\left({\mathbf{x}}_{i}\right)+\lambda \sum\limits _{i{< }j}^{N}\frac{1}{\left\vert {\mathbf{r}}_{i}-{\mathbf{r}}_{j}\right\vert }.\end{align} \tag{ 29 }$$

Any set of orbitals can be used to define the HF Hamiltonian, although either the RHF or UHF orbitals are usually chosen to define the RMP or UMP series respectively. The MP energy at a given order n (i.e., MPn) is then defined as

$$\begin{equation}{E}_{\text{MP}n}=\sum\limits _{k=0}^{n}{E}_{\text{MP}}^{\left(k\right)},\end{equation} \tag{ 30 }$$

where ${E}_{\text{MP}}^{\left(k\right)}$ is the kth-order MP correction and

$$\begin{equation}{E}_{\text{MP}1}={E}_{\text{MP}}^{\left(0\right)}+{E}_{\text{MP}}^{\left(1\right)}={E}_{\text{HF}}.\end{equation} \tag{ 31 }$$

The second-order MP2 energy correction is given by

$$\begin{equation}{E}_{\text{MP}}^{\left(2\right)}=\frac{1}{4}\sum\limits _{ij}\sum\limits _{ab}\frac{{\left\vert \left\langle ij\right\vert \left\vert ab\right\rangle \right\vert }^{2}}{{{\epsilon}}_{i}+{{\epsilon}}_{j}-{{\epsilon}}_{a}-{{\epsilon}}_{b}},\end{equation} \tag{ 32 }$$

where $\left\langle pq\right\vert \left\vert rs\right\rangle =\left\langle pq\vert rs\right\rangle -\left\langle pq\vert sr\right\rangle$ are the anti-symmetrised two-electron integrals in the molecular spin–orbital basis [97]

$$\begin{equation}\left\langle pq\vert rs\right\rangle =\iint d{\mathbf{x}}_{1}\enspace d{\mathbf{x}}_{2}\frac{{\phi }_{p}^{{\ast}}\left({\mathbf{x}}_{1}\right){\phi }_{q}^{{\ast}}\left({\mathbf{x}}_{2}\right){\phi }_{r}\left({\mathbf{x}}_{1}\right){\phi }_{s}\left({\mathbf{x}}_{2}\right)}{\left\vert {\mathbf{r}}_{1}-{\mathbf{r}}_{2}\right\vert }.\end{equation} \tag{ 33 }$$

While most practical calculations generally consider only the MP2 or MP3 approximations, higher order terms can be computed to understand the convergence of the MPn series [15]. A priori, there is no guarantee that this series will provide the smooth convergence that is desirable for a systematically improvable theory. In fact, when the reference HF wave function is a poor approximation to the exact wave function, for example in multi-configurational systems, MP theory can yield highly oscillatory, slowly convergent, or catastrophically divergent results [15, 16, 19, 21, 22, 98]. Furthermore, the convergence properties of the MP series can depend strongly on the choice of restricted or unrestricted reference orbitals.

Although practically convenient for electronic structure calculations, the MP partitioning is not the only possibility and alternative partitionings have been considered [99] including: i) the Epstein–Nesbet (EN) partitioning which consists in taking the diagonal elements of $\hat{H}$ as the zeroth-order Hamiltonian [100, 101], ii) the weak correlation partitioning in which the one-electron part is consider as the unperturbed Hamiltonian ${\hat{H}}^{\left(0\right)}$ and the two-electron part is the perturbation operator $\hat{V}$, and iii) the strong coupling partitioning where the two operators are inverted compared to the weak correlation partitioning [102, 103]. While an in-depth comparison of these different approaches can offer insight into their relative strengths and weaknesses for various situations, we will restrict our current discussion to the convergence properties of the MP expansion.

Among the most desirable properties of any electronic structure technique is the existence of a systematic route to increasingly accurate energies. In the context of MP theory, one would like a monotonic convergence of the perturbation series towards the exact energy such that the accuracy increases as each term in the series is added. If such well-behaved convergence can be established, then our ability to compute individual terms in the series becomes the only barrier to computing the exact correlation in a finite basis set. Unfortunately, the computational scaling of each term in the MP series increases with the perturbation order, and practical calculations must rely on fast convergence to obtain high-accuracy results using only the lowest order terms.

MP theory was first introduced to quantum chemistry through the pioneering works of Bartlett et al in the context of many-body perturbation theory [104], and Pople and co-workers in the context of determinantal expansions [105, 106]. Early implementations were restricted to the fourth-order MP4 approach that was considered to offer state-of-the-art quantitative accuracy [106, 107]. However, it was quickly realised that the MP series often demonstrated very slow, oscillatory, or erratic convergence, with the UMP series showing particularly slow convergence [13–15]. For example, RMP5 is worse than RMP4 for predicting the homolytic barrier fission of He₂ ²⁺ using a minimal basis set, while the UMP series monotonically converges but becomes increasingly slow beyond UMP5 [16]. The first examples of divergent MP series were observed in the N₂ and F₂ diatomics, where low-order RMP and UMP expansions give qualitatively wrong binding curves [17].

The divergence of RMP expansions for stretched bonds can be easily understood from two perspectives [20]. Firstly, the exact wave function becomes increasingly multi-configurational as the bond is stretched, and the RHF wave function no longer provides a qualitatively correct reference for the perturbation expansion. Secondly, the energy gap between the occupied and unoccupied orbitals associated with the stretch becomes increasingly small at larger bond lengths, leading to a divergence, for example, in the MP2 correction (32). In contrast, the origin of slow UMP convergence is less obvious as the reference UHF energy remains qualitatively correct at large bond lengths and the orbital degeneracy is avoided. Furthermore, this slow convergence can also be observed in molecules with a UHF ground state at the equilibrium geometry (e.g., CN⁻), suggesting a more fundamental link with spin-contamination in the reference wave function [18].

Using the UHF framework allows the singlet ground state wave function to mix with triplet wave functions, leading to spin contamination where the wave function is no longer an eigenfunction of the ${\hat{\mathcal{S}}}^{2}$ operator. The link between slow UMP convergence and this spin-contamination was first systematically investigated by Gill et al using the minimal basis H₂ model [19]. In this work, the authors identified that the slow UMP convergence arises from its failure to correctly predict the amplitude of the low-lying double excitation. This erroneous description of the double excitation amplitude has the same origin as the spin-contamination in the reference UHF wave function, creating the first direct link between spin-contamination and slow UMP convergence [19]. Lepetit et al later analysed the difference between perturbation convergence using the UMP and EN partitionings [21]. They argued that the slow UMP convergence for stretched molecules arises from (i) the fact that the MP denominator (see equation (32)) tends to a constant value instead of vanishing, and (ii) the slow convergence of contributions from the singly-excited configurations that strongly couple to the doubly-excited configurations and first appear at fourth-order [21]. Drawing these ideas together, we believe that slow UMP convergence occurs because the single excitations must focus on removing spin-contamination from the reference wave function, limiting their ability to fine-tune the amplitudes of the higher excitations that capture the correlation energy.

A number of spin-projected extensions have been derived to reduce spin-contamination in the wave function and overcome the slow UMP convergence. Early versions of these theories, introduced by Schlegel [108, 109] or Knowles and Handy [110, 111], exploited the 'projection-after-variation' philosophy, where the spin-projection is applied directly to the UMP expansion. These methods succeeded in accelerating the convergence of the projected MP series and were considered as highly effective methods for capturing the electron correlation at low computational cost [111]. However, the use of projection-after-variation leads to gradient discontinuities in the vicinity of the UHF symmetry-breaking point, and can result in spurious minima along a molecular binding curve [108, 110]. More recent formulations of spin-projected perturbations theories have considered the 'variation-after-projection' framework using alternative definitions of the reference Hamiltonian [112, 113]. These methods yield more accurate spin-pure energies without gradient discontinuities or spurious minima.

The behaviour of the RMP and UMP series observed in H₂ can also be illustrated by considering the analytic Hubbard dimer with a complex-valued perturbation strength. In this system, the stretching of the H–H bond is directly mirrored by an increase in the ratio U/t. Using the ground-state RHF reference orbitals leads to the parametrised RMP Hamiltonian

$$\begin{equation}{\mathbf{H}}_{\text{RMP}}\left(\lambda \right)=\left(\begin{matrix}\hfill -2t+U-\lambda U/2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \lambda U/2\hfill \\ \hfill 0\hfill & \hfill U-\lambda U/2\hfill & \hfill \lambda U/2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \lambda U/2\hfill & \hfill U-\lambda U/2\hfill & \hfill 0\hfill \\ \hfill \lambda U/2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 2t+U-\lambda U/2\hfill \end{matrix}\right),\end{equation} \tag{ 34 }$$

which yields the ground-state energy

$$\begin{equation}{E}_{-}\left(\lambda \right)=U-\frac{\lambda U}{2}-\frac{1}{2}\sqrt{{\left(4t\right)}^{2}+{\lambda }^{2}{U}^{2}}.\end{equation} \tag{ 35 }$$

From this expression, the EPs can be identified as λ_EP = ±i4t/U, giving the radius of convergence

$$\begin{equation}{r}_{\text{c}}=\left\vert \frac{4t}{U}\right\vert .\end{equation} \tag{ 36 }$$

Remarkably, these EPs are identical to the exact EPs discussed in section 2.2. The Taylor expansion of the RMP energy can then be evaluated to obtain the kth-order MP correction

$$\begin{equation}{E}_{\text{RMP}}^{\left(k\right)}=U{\delta }_{0,k}-\frac{1}{2}\frac{{U}^{k}}{{\left(4t\right)}^{k-1}}\left(\begin{matrix}\hfill 1/2\hfill \\ \hfill k/2\hfill \end{matrix}\right).\end{equation} \tag{ 37 }$$

The RMP series is convergent at λ = 1 for U = 3.5t with r_c > 1, as illustrated for the individual terms at each perturbation order in figure 4(b). In contrast, for U = 4.5t one finds r_c < 1, and the RMP series becomes divergent at λ = 1. The corresponding Riemann surfaces for U = 3.5t and 4.5t are shown in figures 4(a) and (c), respectively, with the single EP at λ_EP (black dot). We illustrate the surface $\left\vert \lambda \right\vert =1$ using a vertical cylinder of unit radius to provide a visual aid for determining if the series will converge at the physical case λ = 1. For the divergent case, λ_EP lies inside this unit cylinder, while in the convergent case λ_EP 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. Note that, when λ_EP lies on the unit cylinder, we cannot a priori determine whether the perturbation series will converge or not.

Zoom In Zoom Out Reset image size

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{equation}{\mathbf{H}}_{\text{UMP}}\left(\lambda \right)=\left(\begin{matrix}\hfill -2{t}^{2}\lambda /U\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 2{t}^{2}\lambda /U\hfill \\ \hfill 0\hfill & \hfill U-2{t}^{2}\lambda /U\hfill & \hfill 2{t}^{2}\lambda /U\hfill & \hfill 2t\sqrt{{U}^{2}-{\left(2t\right)}^{2}}\lambda /U\hfill \\ \hfill 0\hfill & \hfill 2{t}^{2}\lambda /U\hfill & \hfill U-2{t}^{2}\lambda /U\hfill & \hfill -2t\sqrt{{U}^{2}-{\left(2t\right)}^{2}}\lambda /U\hfill \\ \hfill 2{t}^{2}\lambda /U\hfill & \hfill 2t\sqrt{{U}^{2}-{\left(2t\right)}^{2}}\lambda /U\hfill & \hfill -2t\sqrt{{U}^{2}-{\left(2t\right)}^{2}}\lambda /U\hfill & \hfill 2U\left(1-\lambda \right)+6{t}^{2}\lambda /U\hfill \end{matrix}\right).\end{equation} \tag{ 38 }$$

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 figure 5. These numerical values reveal that the UMP ground-state series has r_c > 1 for all U/t and always converges. However, in the strong correlation limit (large U/t), this radius of convergence tends to unity, indicating that the convergence of the corresponding UMP series becomes increasingly slow. Furthermore, the doubly-excited state using the ground-state UHF orbitals has r_c < 1 for almost any value of U/t, reaching the limiting value of 1/2 for U/t → ∞. Hence, the excited-state UMP series will always diverge.

Zoom In Zoom Out Reset image size

The convergence behaviour can be further elucidated by considering the full structure of the UMP energies in the complex λ-plane (see figures 6(a) and (c)). These Riemann surfaces are illustrated for U = 3t and 7t alongside the perturbation terms at each order in figure 6(b). At U = 3t, the RMP series is convergent at λ = 1, while RMP becomes divergent at λ = 1 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 (figure 5).

Zoom In Zoom Out Reset image size

As the UHF orbitals break the spatial and spin symmetry, new coupling terms emerge between the electronic states that cause fundamental changes to the structure of EPs in the complex λ-plane. For example, while the RMP energy shows only one EP between the ground and doubly-excited states (figure 4), the UMP energy has two pairs of complex-conjugate 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 (figure 6). 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 H₂ [19] can then be seen as this EP moves increasingly close to the unit cylinder at large U/t and r_c 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 [96], and thus it is not surprising that this process proceeds very slowly as the perturbation order is increased.

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 [114]. They established two general classes: 'class A' systems that exhibit monotonic convergence; and 'class B' systems for which convergence is erratic after initial oscillations. 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 [114]. 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 orbital relaxation to doubly-excited configurations, and these contributions lead to oscillations of the total correlation energy.

Using these classifications, Cremer and He then introduced simple extrapolation formulas for estimating the exact correlation energy ΔE using terms up to MP6 [114]

$$\begin{equation} {\Delta}{E}_{\text{A}}={E}_{\text{MP}}^{\left(2\right)}+{E}_{\text{MP}}^{\left(3\right)}+{E}_{\text{MP}}^{\left(4\right)}+\frac{{E}_{\text{MP}}^{\left(5\right)}}{1-\left({E}_{\text{MP}}^{\left(6\right)}/{E}_{\text{MP}}^{\left(5\right)}\right)},\end{equation} \tag{ 39a }$$

$$\begin{equation} {\Delta}{E}_{\text{B}}={E}_{\text{MP}}^{\left(2\right)}+{E}_{\text{MP}}^{\left(3\right)}+\left({E}_{\text{MP}}^{\left(4\right)}+{E}_{\text{MP}}^{\left(5\right)}\right)\mathrm{exp}\left({E}_{\text{MP}}^{\left(6\right)}/{E}_{\text{MP}}^{\left(5\right)}\right).\end{equation} \tag{ 39b }$$

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 4, we consider more advanced extrapolation routines that take account of EPs in the complex λ-plane.

In the late 90s, Olsen et al discovered an even more concerning behaviour of the MP series [57]. They showed that the series could be divergent even in systems that were considered to be well understood, such as Ne or the HF molecule [57, 115]. Cremer and He had already studied these two systems and classified them as class B systems [114]. 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 [116, 117]. 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 [118, 119].

Olsen et al investigated the causes of these divergences and the different types of convergence by analysing the relation between the dominant singularity (i.e., the closest singularity to the origin) and the convergence behaviour of the series [58]. Their analysis is based on Darboux's theorem: [77]

'In the limit of large order, the series coefficients become equivalent to the Taylor series coefficients of the singularity closest to the origin'.

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 reference [57], 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 [58]. By modelling this avoided crossing using a two-state Hamiltonian, one can obtain an approximation for the dominant singularities as the EPs of the two-state matrix

$$\begin{equation}\underset{\mathbf{H}}{\underbrace{\left(\begin{matrix}\hfill \alpha \hfill & \hfill \delta \hfill \\ \hfill \delta \hfill & \hfill \beta \hfill \end{matrix}\right)}}=\underset{{\mathbf{H}}^{\left(0\right)}}{\underbrace{\left(\begin{matrix}\hfill \alpha +{\alpha }_{\text{s}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \beta +{\beta }_{\text{s}}\hfill \end{matrix}\right)}}+\underset{\mathbf{V}}{\underbrace{\left(\begin{matrix}\hfill -{\alpha }_{\text{s}}\hfill & \hfill \delta \hfill \\ \hfill \delta \hfill & \hfill -{\beta }_{\text{s}}\hfill \end{matrix}\right)}},\end{equation} \tag{ 40 }$$

where the diagonal matrix is the unperturbed Hamiltonian matrix H⁽⁰⁾ with level shifts α_s and β_s, and V 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 [58]. 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 CH₂ in a diffuse, yet rather small basis set, the series is convergent at λ = 1 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 (figure 4). Furthermore, the authors demonstrated that the divergence for Ne is due to a back-door intruder state that arise when the ground state undergoes sharp avoided crossings with highly diffuse excited states. This divergence is related to a more fundamental critical point in the MP energy surface that we will discuss in section 3.5.

Finally, reference [57] proved that the extrapolation formulas of Cremer and He [114] (see equations (39a) and (39b)) are not mathematically motivated when considering the complex singularities causing the divergence, and therefore cannot be applied for all systems. For example, the HF molecule 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 [59]

$$\begin{equation}\underset{\mathbf{H}}{\underbrace{\left(\begin{matrix}\hfill \alpha \hfill & \hfill {\delta }_{1}\hfill \\ \hfill {\delta }_{2}\hfill & \hfill \beta \hfill \end{matrix}\right)}}=\underset{{\mathbf{H}}^{\left(0\right)}}{\underbrace{\left(\begin{matrix}\hfill \alpha \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \beta +\gamma \hfill \end{matrix}\right)}}+\underset{\mathbf{V}}{\underbrace{\left(\begin{matrix}\hfill 0\hfill & \hfill {\delta }_{2}\hfill \\ \hfill {\delta }_{1}\hfill & \hfill -\gamma \hfill \end{matrix}\right)}}.\end{equation} \tag{ 41 }$$

This extension allows various choices of perturbation to be analysed, including coupled cluster perturbation expansions [120–124] and other non-Hermitian perturbation methods. Note that new forms of perturbation expansions only occur when the sign of δ₁ and δ₂ differ. Using this 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 [59]. 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. The geometric archetype appears to be the most common for MP expansions [59], but the ripples archetype corresponds to some of the early examples of MP convergence [15, 21, 22, 98]. The three remaining Hermitian archetypes seem to be rarely observed in MP perturbation theory. In contrast, the non-Hermitian coupled cluster perturbation theory [120–124], exhibits a range of archetypes including the interspersed zigzag, triadic, ripple, geometric, and zigzag-geometric forms. This analysis highlights the importance of the primary singularity in controlling the high-order convergence, regardless of whether this point is inside or outside the complex unit circle [15, 58].

In the early 2000s, Stillinger reconsidered the mathematical origin behind the divergent series with odd–even sign alternation [125]. This type of convergence behaviour corresponds to Cremer and He's class B systems with closely spaced electron pairs and includes Ne, HF, F⁻, and H₂O [114]. Stillinger proposed that these series diverge due to a dominant singularity on the negative real λ axis, corresponding to a multielectron autoionisation threshold [125]. To understand Stillinger's argument, consider the parametrised MP Hamiltonian in the form

$$\begin{align}\hat{H}\left(\lambda \right)& =\sum\limits _{i}^{N}\left[\text{independent}\;\text{of}\;\lambda {\overbrace{-\frac{1}{2}{\mathbf{\nabla }}_{i}^{2}-\sum\limits _{A}^{M}\frac{{Z}_{A}}{\left\vert {\mathbf{r}}_{i}-{\mathbf{R}}_{A}\right\vert }}}\right.\\ & \quad +\left.\underset{\text{repulsive}\;\text{for}\;\lambda {< }1}{\underbrace{\left(1-\lambda \right){\hat{v}}_{\text{HF}}\left({\mathbf{x}}_{i}\right)}}+\underset{\text{attractive}\;\text{for}\;\lambda {< }0}{\underbrace{\lambda \sum\limits _{i{< }j}^{N}\frac{1}{\vert {\mathbf{r}}_{i}-{\mathbf{r}}_{j}\vert }}}\right].\end{align} \tag{ 42 }$$

The mean-field potential ${\hat{v}}_{\text{HF}}$ essentially represents a negatively charged field with the spatial extent controlled by the extent of the HF orbitals, usually located close to the nuclei. When λ is negative, the mean-field potential becomes increasingly repulsive, while the explicit two-electron Coulomb interaction becomes attractive. There is therefore a negative critical point λ_c where it becomes energetically favourable for the electrons to dissociate and form a bound cluster at an infinite separation from the nuclei [125]. This autoionisation effect is closely related to the critical point for electron binding in two-electron atoms (see reference [126]). Furthermore, a similar set of critical points exists along the positive real axis, corresponding to single-electron ionisation processes [127]. While these critical points are singularities on the real axis, their exact mathematical form is difficult to identify and remains an open question.

To further develop the link between the critical point and types of MP convergence, Sergeev and Goodson investigated the relationship with the location of the dominant singularity that controls the radius of convergence [128]. They demonstrated that the dominant singularity in class A systems corresponds to an EP with a positive real component, where the magnitude of the imaginary component controls the oscillations in the signs of successive MP terms [129, 130]. In contrast, class B systems correspond to a dominant singularity on the negative real λ axis representing the MP critical point. The divergence of class B systems, which contain closely spaced electrons (e.g., F⁻), can then be understood as the HF potential ${\hat{v}}_{\text{HF}}$ is relatively localised and the autoionization is favoured at negative λ values closer to the origin. With these insights, they regrouped the systems into new classes: i) α singularities which have 'large' imaginary parts, and ii) β singularities which have very small imaginary parts [128, 131].

The existence of the MP critical point can also explain why the divergence observed by Olsen et al in the Ne atom and the HF molecule occurred when diffuse basis functions were included [57]. Clearly diffuse basis functions are required for the electrons to dissociate from the nuclei, and indeed using only compact basis functions causes the critical point to disappear. While a finite basis can only predict complex-conjugate branch point singularities, the critical point is modelled by a cluster of sharp avoided crossings between the ground state and high-lying excited states [127]. Alternatively, Sergeev et al demonstrated that the inclusion of a 'ghost' atom also allows the formation of the critical point as the electrons form a bound cluster occupying the ghost atom orbitals [127]. This effect explains the origin of the divergence in the HF molecule as the fluorine valence electrons jump to the hydrogen at a sufficiently negative λ value [127]. Furthermore, the two-state model of Olsen and collaborators [58] was simply too minimal to understand the complexity of divergences caused by the MP critical point.

When a Hamiltonian is parametrised by a variable such as λ, the existence of abrupt changes in the eigenstates as a function of λ indicate the presence of a zero-temperature QPT [29, 69–75]. Meanwhile, as an avoided crossing becomes increasingly sharp, the corresponding EPs move increasingly close to the real axis. When these points converge on the real axis, they effectively 'annihilate' each other and no longer behave as EPs. Instead, they form a 'critical point' singularity that resembles a conical intersection, and the convergence of a pair of complex-conjugate EPs on the real axis is therefore diagnostic of a QPT [132, 133].

Since the MP critical point corresponds to a singularity on the real λ axis, it can immediately be recognised as a QPT with respect to varying the perturbation parameter λ. However, a conventional QPT can only occur in the thermodynamic limit, which here is analogous to the complete basis set limit [134]. The MP critical point and corresponding β singularities in a finite basis must therefore be modelled by pairs of complex-conjugate EPs that tend towards the real axis, exactly as described by Sergeev et al [127] In contrast, α singularities correspond to large avoided crossings that are indicative of low-lying excited states which share the symmetry of the ground state [128], and are thus not manifestations of a QPT. Notably, since the exact MP critical point corresponds to the interaction between a bound state and the continuum, its functional form is more complicated than a conical intersection and remains an open question.

The simplified site basis of the Hubbard dimer makes explicitly modelling the ionisation continuum impossible. Instead, we can use an asymmetric version of the Hubbard dimer [67, 68] where we consider one of the sites as a 'ghost atom' that acts as a destination for ionised electrons being originally localised on the other site. To mathematically model this scenario in this asymmetric Hubbard dimer, we introduce a one-electron potential − on the left site to represent the attraction between the electrons and the model 'atomic' nucleus, where we define > 0. The reference Slater determinant for a doubly-occupied atom can be represented using RHF orbitals (see equation (21)) with ${\theta }_{\text{RHF}}^{\enspace \alpha }={\theta }_{\text{RHF}}^{\enspace \beta }=0$, which corresponds to strictly localising the two electrons on the left site. With this representation, the parametrised asymmetric RMP Hamiltonian becomes

$$\begin{equation}{\mathbf{H}}_{\text{asym}}\left(\lambda \right)=\left(\begin{matrix}\hfill 2\left(U-{\epsilon}\right)-\lambda U\hfill & \hfill -\lambda t\hfill & \hfill -\lambda t\hfill & \hfill 0\hfill \\ \hfill -\lambda t\hfill & \hfill \left(U-{\epsilon}\right)-\lambda U\hfill & \hfill 0\hfill & \hfill -\lambda t\hfill \\ \hfill -\lambda t\hfill & \hfill 0\hfill & \hfill \left(U-{\epsilon}\right)-\lambda U\hfill & \hfill -\lambda t\hfill \\ \hfill 0\hfill & \hfill -\lambda t\hfill & \hfill -\lambda t\hfill & \hfill \lambda U\hfill \end{matrix}\right).\end{equation} \tag{ 43 }$$

For the ghost site to perfectly represent ionised electrons, the hopping term between the two sites must vanish (i.e., t = 0). This limit corresponds to the dissociative regime in the asymmetric Hubbard dimer as discussed in reference [68], and the RMP energies become

$$\begin{equation}{E}_{-}=2\left(U-{\epsilon}\right)-\lambda U,\end{equation} \tag{ 44a }$$

$$\begin{equation}{E}_{\text{S}}=\left(U-{\epsilon}\right)-\lambda U,\end{equation} \tag{ 44b }$$

$$\begin{equation}{E}_{+}=U\lambda ,\end{equation} \tag{ 44c }$$

as shown in figure 7(a) (dashed lines). The RMP critical point then corresponds to the intersection E₋ = E₊, giving the critical λ value

$$\begin{equation}{\lambda }_{\text{c}}=1-\frac{{\epsilon}}{U}.\end{equation} \tag{ 45 }$$

Clearly the radius of convergence ${r}_{\text{c}}=\left\vert {\lambda }_{\text{c}}\right\vert$ is controlled directly by the ratio /U, with a convergent RMP series at λ = 1 occurring for > 2U. The on-site repulsion U controls the strength of the HF potential localised around the 'atomic site', with a stronger repulsion encouraging the electrons to be ionised at a less negative value of λ. Large U can be physically interpreted as strong electron repulsion effects in electron dense molecules. In contrast, smaller gives a weaker attraction to the atomic site, representing strong screening of the nuclear attraction by core and valence electrons, and again a less negative λ is required for ionisation to occur. Both of these factors are common in atoms on the right-hand side of the periodic table, e.g., F, O, Ne. Molecules containing these atoms are therefore often class β systems with a divergent RMP series due to the MP critical point [128, 131].

Zoom In Zoom Out Reset image size

The critical point in the exact case t = 0 is represented by the gradient discontinuity in the ground-state energy on the negative real λ axis (figure 7(a): solid lines), mirroring the behaviour of a QPT [134]. The autoionisation process is manifested by a sudden drop in the 'atomic site' electron density ρ_atom (figure 8). However, in practical calculations performed with a finite basis set, the critical point is modelled as a cluster of branch points close to the real axis. The use of a finite basis can be modelled in the asymmetric dimer by making the second site a less idealised destination for the ionised electrons with a non-zero (yet small) hopping term t. Taking the value t = 0.1 (figure 7(a): dashed lines), the critical point becomes an avoided crossing with a complex-conjugate pair of EPs close to the real axis (figure 7(b)). In contrast to the exact critical point with t = 0, the ground-state energy remains smooth through this avoided crossing, with a more gradual drop in the atomic site density. In the limit t → 0, these EPs approach the real axis (figure 7(c)) and the avoided crossing becomes a gradient discontinuity, mirroring Sergeev's discussion on finite basis set representations of the MP critical point [131].

Zoom In Zoom Out Reset image size

Returning to the symmetric Hubbard dimer, we showed in section 3.3 that the slow convergence of the strongly correlated UMP series was due to a complex-conjugate pair of EPs just outside the radius of convergence. These EPs have positive real components and small imaginary components (see figure 6), suggesting a potential connection to MP critical points and QPTs (see section 3.5). For λ > 1, the HF potential becomes an attractive component in Stillinger's Hamiltonian displayed in equation (42), while the explicit electron–electron interaction becomes increasingly repulsive. Closed-shell critical points along the positive real λ axis may then represent points where the two-electron repulsion overcomes the attractive HF potential and a single electron dissociates from the molecule (see reference [131]).

In contrast, spin symmetry-breaking in the UMP reference creates different HF potentials for the spin-up and spin-down electrons. Consider one of the two reference UHF solutions where the spin-up and spin-down electrons are localised on the left and right sites respectively. The spin-up HF potential will then be a repulsive interaction from the spin-down electron density that is centred around the right site (and vice-versa). As λ becomes greater than 1 and the HF potentials become attractive, there will be a sudden driving force for the electrons to swap sites. This swapping process can also be represented as a double excitation, and thus an avoided crossing will occur for λ ⩾ 1 (figure 9(a)). While this appears to be an avoided crossing between the ground and first-excited state, the presence of an earlier excited-state avoided crossing means that the first-excited state qualitatively represents the reference double excitation for λ > 1/2. We can visualise this swapping process by considering the difference in the electron density on the left and right sites, defined for each spin as

$$\begin{equation}{\Delta}{\rho }^{\sigma }={\rho }_{\mathcal{R}}^{\sigma }-{\rho }_{\mathcal{L}}^{\sigma },\end{equation} \tag{ 46 }$$

where ${\rho }_{\mathcal{L}}^{\sigma }$ (${\rho }_{\mathcal{R}}^{\sigma }$) is the spin-σ electron density on the left (right) site. This density difference is shown for the UMP ground-state at U = 5t in figure 10 (solid lines). Here, the transfer of the spin-up electron from the right site to the left site can be seen as λ passes through 1 (and similarly for the spin-down electron).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 'sharpness' of the avoided crossing is controlled by the correlation strength U/t. For small U/t, the HF potentials will be weak and the electrons will delocalise over the two sites, both in the UHF reference and the exact wave function. This delocalisation dampens the electron swapping process and leads to a 'shallow' avoided crossing (solid lines in figure 9(a)) that corresponds to EPs with non-zero imaginary components (figure 9(b)). As U/t becomes larger, the HF potentials become stronger and the on-site repulsion dominates the hopping term to make electron delocalisation less favourable. In other words, the electrons localise on individual sites to form a Wigner crystal. These effects create a stronger driving force for the electrons to swap sites until, eventually, this swapping occurs suddenly at λ = 1, as shown for U = 50t in figure 10 (dashed lines). In this limit, the ground-state EPs approach the real axis (figure 9(c)) and the avoided crossing creates a gradient discontinuity in the ground-state energy (dashed lines in figure 9(a)). We therefore find that, in the strong correlation limit, the symmetry-broken ground-state EP becomes a new type of MP critical point and represents a QPT as the perturbation parameter λ is varied. Remarkably, this argument explains why the dominant UMP singularity lies so close, but always outside, the radius of convergence (see figure 5).

The failure of a Taylor series for correctly modelling the MP energy function E(λ) arises because one is trying to model a complicated function containing multiple branches, branch points, and singularities using a simple polynomial of finite order. A truncated Taylor series can only predict a single sheet and does not have enough flexibility to adequately describe functions such as the MP energy. Alternatively, the description of complex energy functions can be significantly improved by introducing Padé approximants [136], and related techniques [56, 137].

A Padé approximant can be considered as the best approximation of a function by a rational function of given order. More specifically, a [d_A /d_B ] Padé approximant is defined as

$$\begin{equation}{E}_{\left[{d}_{A}/{d}_{B}\right]}\left(\lambda \right)=\frac{A\left(\lambda \right)}{B\left(\lambda \right)}=\frac{{\sum }_{k=0}^{{d}_{A}}{a}_{k}\enspace {\lambda }^{k}}{1+{\sum }_{k=1}^{{d}_{B}}{b}_{k}\enspace {\lambda }^{k}},\end{equation} \tag{ 47 }$$

where the coefficients of the polynomials A(λ) and B(λ) are determined by collecting and comparing terms for each power of λ with the low-order terms in the Taylor series expansion. Padé approximants are extremely useful in many areas of physics and chemistry [138–141] as they can model poles, which appear at the roots of B(λ). However, they are unable to model functions with square-root branch points (which are ubiquitous in the singularity structure of perturbative methods) and more complicated functional forms appearing at critical points (where the nature of the solution undergoes a sudden transition). Despite this limitation, the successive diagonal Padé approximants (i.e., d_A = d_B ) often define a convergent perturbation series in cases where the Taylor series expansion diverges.

Figure 11 illustrates the improvement provided by diagonal Padé approximants 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 1 gathers estimates of the RMP ground-state energy at λ = 1 provided by various truncated Taylor series and Padé approximants for these two values of the ratio U/t. While the truncated Taylor series converges laboriously to the exact energy as the truncation degree increases at U/t = 3.5, the Padé approximants yield much more accurate results. Furthermore, the distance of the closest pole to the origin $\left\vert {\lambda }_{\text{c}}\right\vert$ in the Padé approximants indicate that they provide a relatively good approximation to the position of the true branch point singularity in the RMP energy. For U/t = 4.5, the Taylor series expansion performs worse and eventually diverges, while the Padé approximants still offer relatively accurate energies and recovers a convergent series.

Zoom In Zoom Out Reset image size

Table 1. RMP ground-state energy estimate at λ = 1 of the Hubbard dimer provided by various truncated Taylor series and Padé approximants at U/t = 3.5 and 4.5. We also report the distance of the closest pole to the origin $\left\vert {\lambda }_{\text{c}}\right\vert$ provided by the diagonal Padé approximants.

		$\left\vert {\lambda }_{\text{c}}\right\vert$		E−(λ = 1)
Method	Degree	U/t = 3.5	U/t = 4.5	U/t = 3.5	U/t = 4.5
Taylor	2			−1.015 63	−1.015 63
	3			−1.015 63	−1.015 63
	4			−0.869 08	−0.615 17
	5			−0.869 08	−0.615 17
	6			−0.925 18	−0.868 58
Padé	[1/1]	2.29	1.78	−1.611 11	−2.642 86
	[2/2]	2.29	1.78	−0.821 24	−0.484 46
	[3/3]	1.73	1.34	−0.919 95	−0.819 29
	[4/4]	1.47	1.14	−0.905 79	−0.748 66
	[5/5]	1.35	1.05	−0.907 78	−0.762 77
Exact		1.14	0.89	−0.907 54	−0.760 40

We can expect the UMP energy function to be much more challenging to model properly as it contains three connected branches (see figures 6(a) and (c)). Figure 12 and table 2 indicate that this is indeed the case. In particular, figure 12 illustrates that the Padé approximants are trying to model the square root branch point that lies close to λ = 1 by placing a pole on the real axis (e.g., [3/3]) or with a very small imaginary component (e.g., [4/4]). The proximity of these poles to the physical point λ = 1 means that any error in the Padé functional form becomes magnified in the estimate of the exact energy, as seen for the low-order approximants in table 2. However, with sufficiently high degree polynomials, one obtains accurate estimates for the position of the closest singularity and the ground-state energy at λ = 1, even in cases where the convergence of the UMP series is incredibly slow (see figure 6(b)).

Zoom In Zoom Out Reset image size

Table 2. Estimate for the distance of the closest singularity (pole or branch point) to the origin $\left\vert {\lambda }_{\text{c}}\right\vert$ in the UMP energy function of the Hubbard dimer provided by various truncated Taylor series and approximants at U/t = 3 and 7. The truncation degree of the Taylor expansion n of E(λ) and the number of branch points n_bp = max(2d_p , d_q + d_r ) generated by the quadratic approximants are also reported.

				$\left\vert {\lambda }_{\text{c}}\right\vert$		E−(λ = 1)
Method		n	nbp	U/t = 3	U/t = 7	U/t = 3	U/t = 7
Taylor			2				−0.740 74	−0.291 55
		3				−0.781 89	−0.296 90
		4				−0.822 13	−0.302 25
		5				−0.857 69	−0.307 58
		6				−0.888 82	−0.312 89
Padé	[1/1]	2		9.000	49.00	−0.750 00	−0.291 67
	[2/2]	4		0.974	1.003	0.750 00	−17.93 75
	[3/3]	6		1.141	1.004	−1.108 96	−1.498 56
	[4/4]	8		1.068	1.003	−0.853 96	−0.335 96
	[5/5]	10		1.122	1.004	−0.972 54	−0.355 13
Quadratic	[2/1, 2]	6	4	1.086	1.003	−1.010 09	−0.534 72
	[2/2, 2]	7	4	1.082	1.003	−1.005 53	−0.534 63
	[3/2, 2]	8	6	1.082	1.001	−1.005 68	−0.524 73
	[3/2, 3]	9	6	1.071	1.002	−0.999 73	−0.531 02
	[3/3, 3]	10	6	1.071	1.002	−0.999 66	−0.531 03
(Pole-free)	[3/0, 2]	6	6	1.059	1.003	−1.137 12	−0.571 99
	[3/0, 3]	7	6	1.073	1.002	−1.003 35	−0.531 13
	[3/0, 4]	8	6	1.071	1.002	−1.000 74	−0.531 16
	[3/0, 5]	9	6	1.070	1.002	−1.000 42	−0.531 14
	[3/0, 6]	10	6	1.070	1.002	−1.000 39	−0.531 13
Exact				1.069	1.002	−1.000 00	−0.531 13

Quadratic approximants are designed to model the singularity structure of the energy function E(λ) via a generalised version of the square-root singularity expression [77, 135, 142]

$$\begin{equation}{E}_{\left[{d}_{P}/{d}_{Q},{d}_{R}\right]}\left(\lambda \right)=\frac{1}{2Q\left(\lambda \right)}\left[P\left(\lambda \right){\pm}\sqrt{{P}^{2}\left(\lambda \right)-4Q\left(\lambda \right)R\left(\lambda \right)}\right],\end{equation} \tag{ 48 }$$

with the polynomials

$$\begin{equation}P\left(\lambda \right)=\sum\limits _{k=0}^{{d}_{P}}{p}_{k}{\lambda }^{k},Q\left(\lambda \right)=\sum\limits _{k=0}^{{d}_{Q}}{q}_{k}{\lambda }^{k},R\left(\lambda \right)=\sum\limits _{k=0}^{{d}_{R}}{r}_{k}{\lambda }^{k},\end{equation} \tag{ 49 }$$

defined such that d_P + d_Q + d_R = n − 1, and n is the truncation order of the Taylor series of E(λ). Recasting equation (48) as a second-order expression in E(λ), i.e.,

$$\begin{equation}Q\left(\lambda \right){E}^{2}\left(\lambda \right)-P\left(\lambda \right)E\left(\lambda \right)+R\left(\lambda \right)\sim \mathcal{O}{\lambda }^{n+1},\end{equation} \tag{ 50 }$$

and substituting E(λ) by its nth-order expansion and the polynomials by their respective expressions (49) yields n + 1 linear equations for the coefficients p_k , q_k , and r_k (where we are free to assume that q₀ = 1). A quadratic approximant, characterised by the label [d_P /d_Q , d_R ], generates, by construction, n_bp = max(2d_p , d_q + d_r ) branch points at the roots of the polynomial P²(λ) − 4Q(λ)R(λ) and d_q poles at the roots of Q(λ).

Generally, the diagonal sequence of quadratic approximant, i.e., [0/0, 0], [1/0, 0], [1/0, 1], [1/1, 1], [2/1, 1], is of particular interest as the order of the corresponding Taylor series increases on each step. However, while a quadratic approximant can reproduce multiple branch points, it can only describe a total of two branches. This constraint can hamper the faithful description of more complicated singularity structures such as the MP energy surface. Despite this limitation, reference [129] demonstrates that quadratic approximants provide convergent results in the most divergent cases considered by Olsen and collaborators [57, 115] and Leininger et al [98].

As a note of caution, reference [135] suggests that low-order quadratic approximants can struggle to correctly model the singularity structure when the energy function has poles in both the positive and negative half-planes. In such a scenario, the quadratic approximant will tend to place its branch points in-between, potentially introducing singularities quite close to the origin. The remedy for this problem involves applying a suitable transformation of the complex plane (such as a bilinear conformal mapping) which leaves the points at λ = 0 and λ = 1 unchanged [143].

For the RMP series of the Hubbard dimer, the [0/0, 0] and [1/0, 0] quadratic approximants are quite poor approximations, but the [1/0, 1] version perfectly models the RMP energy function by predicting a single pair of EPs at λ_EP = ±i4t/U. This is expected from the form of the RMP energy (see equation (35)), which matches the ideal target for quadratic approximants. Furthermore, the greater flexibility of the diagonal quadratic approximants provides a significantly improved model of the UMP energy in comparison to the Padé approximants or Taylor series. In particular, these quadratic approximants provide an effective model for the avoided crossings (figure 12) and an improved estimate for the distance of the closest branch point to the origin. Table 2 shows that they provide remarkably accurate estimates of the ground-state energy at λ = 1.

While the diagonal quadratic approximants provide significantly improved estimates of the ground-state energy, we can use our knowledge of the UMP singularity structure to develop even more accurate results. We have seen in previous sections that the UMP energy surface contains only square-root branch cuts that approach the real axis in the limit U/t → ∞. Since there are no true poles on this surface, we can obtain more accurate quadratic approximants by taking d_q = 0 and increasing d_r to retain equivalent accuracy in the square-root term (see equation (48)). Figure 13 illustrates this improvement for the pole-free [3/0, 4] quadratic approximant compared to the [3/2, 2] approximant with the same truncation degree in the Taylor expansion. Clearly, modelling the square-root branch point using d_q = 2 has the negative effect of introducing spurious poles in the energy, while focussing purely on the branch point with d_q = 0 leads to a significantly improved model. Table 2 shows that these pole-free quadratic approximants provide a rapidly convergent series with essentially exact energies at low order.

Zoom In Zoom Out Reset image size

Finally, to emphasise the improvement that can be gained by using either Padé, diagonal quadratic, or pole-free quadratic approximants, we collect the energy and error obtained using only the first 10 terms of the UMP Taylor series in table 3. The accuracy of these approximants reinforces how our understanding of the MP energy surface in the complex plane can be leveraged to significantly improve estimates of the exact energy using low-order perturbation expansions.

While the Padé and quadratic approximants can yield a convergent series representation in cases where the standard MP series diverges, there is no guarantee that the rate of convergence will be fast enough for low-order approximations to be useful. However, these low-order partial sums or approximants often contain a remarkable amount of information that can be used to extract further information about the exact result. The Shanks transformation presents one approach for extracting this information and accelerating the rate of convergence of a sequence [56, 144].

Consider the partial sums ${S}_{n}={\sum }_{k=0}^{n}{s}_{k}$ defined from the truncated summation of an infinite series $S={\sum }_{k=0}^{\infty }{s}_{k}$. If the series converges, then the partial sums will tend to the exact result

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {S}_{n}=S.\end{equation} \tag{ 51 }$$

The Shanks transformation attempts to generate increasingly accurate estimates of this limit by defining a new series as

$$\begin{equation}T\left({S}_{n}\right)=\frac{{S}_{n+1}{S}_{n-1}-{S}_{n}^{2}}{{S}_{n+1}-2{S}_{n}+{S}_{n-1}}.\end{equation} \tag{ 52 }$$

This series can converge faster than the original partial sums and can thus provide greater accuracy using only the first few terms in the series. However, it is only designed to accelerate converging partial sums with the approximate form S_n ≈ S + αβⁿ . Furthermore, while this transformation can accelerate the convergence of a series, there is no guarantee that this acceleration will be fast enough to significantly improve the accuracy of low-order approximations.

To the best of our knowledge, the Shanks transformation has never previously been applied to accelerate the convergence of the MP series. We have therefore applied it to the convergent Taylor series, Padé approximants, and quadratic approximants for RMP and UMP in the symmetric Hubbard dimer. The UMP approximants converge too slowly for the Shanks transformation to provide any improvement, even in the case where the quadratic approximants are already very accurate. In contrast, acceleration of the diagonal Padé approximants for the RMP cases can significantly improve the estimate of the energy using low-order perturbation terms, as shown in table 4. Even though the RMP series diverges at U/t = 4.5, the combination of diagonal Padé approximants with the Shanks transformation reduces the absolute error in the best energy estimate to 0.002% using only the first 10 terms in the Taylor series. This remarkable result indicates just how much information is contained in the first few terms of a perturbation series, even if it diverges.

Recently, Mihálka et al have studied the effect of different partitionings, such as MP or EN theory, on the position of branch points and the convergence properties of Rayleigh–Schrödinger perturbation theory [61] (see also references [145–147]). Taking the equilibrium and stretched water structures as an example, they estimated the radius of convergence using quadratic Padé approximants. The EN partitioning provided worse convergence properties than the MP partitioning, which is believed to be because the EN denominators are generally smaller than the MP denominators. To remedy the situation, they showed that introducing a suitably chosen level shift parameter can turn a divergent series into a convergent one by increasing the magnitude of these denominators [61]. However, like the UMP series in stretched H₂ [21], the cost of larger denominators is an overall slower rate of convergence.

In a later study by the same group, they used analytic continuation techniques to resum a divergent MP series such as a stretched water molecule [60]. Any MP series truncated at a given order n can be used to define the scaled function

$$\begin{equation}{E}_{\text{MP}n}\left(\lambda \right)=\sum\limits _{k=0}^{n}{\lambda }^{k}{E}_{\text{MP}}^{\left(k\right)}.\end{equation} \tag{ 53 }$$

Reliable estimates of the energy can be obtained for values of λ where the MP series is rapidly convergent (i.e., for $\left\vert \lambda \right\vert {< }{r}_{\text{c}}$), as shown in figure 14 for the RMP10 series of the symmetric Hubbard dimer with U/t = 4.5. These values can then be analytically continued using a polynomial- or Padé-based fit to obtain an estimate of the exact energy at λ = 1. However, choosing the functional form for the best fit remains a difficult and subtle challenge.

Zoom In Zoom Out Reset image size

This technique was first generalised using complex scaling parameters to construct an analytic continuation by solving the Laplace equations [148]. It was then further improved by introducing Cauchy's integral formula [62]

$$\begin{equation}E\left(\lambda \right)={\frac{1}{2\pi \mathrm{i}}\oint }_{\mathcal{C}}\frac{E\left({\lambda }^{\prime }\right)}{{\lambda }^{\prime }-\lambda },\end{equation} \tag{ 54 }$$

which states that the value of the energy can be computed at λ inside the complex contour $\mathcal{C}$ using only the values along the same contour. Starting from a set of points in a 'trusted' region where the MP series is convergent, their approach self-consistently refines estimates of the E(λ') values on a contour that includes the physical point λ = 1. The shape of this contour is arbitrary, but there must be no branch points or other singularities inside the contour. Once the contour values of E(λ') are converged, Cauchy's integral formula equation (54) can be invoked to compute the value at E(λ = 1) and obtain a final estimate of the exact energy. The authors illustrate this protocol for the dissociation curve of LiH and the stretched water molecule and obtained encouragingly accurate results [62].