% Manuscript for Advances in Quantum Chemistry, Goscinski volume.
\documentclass[12pt]{article}
\usepackage{graphicx}
\input epsf
\title{Singularity Structure of M{\o}ller-Plesset Perturbation
Theory\protect\footnote{This article is dedicated to Osvaldo
Goscinski in recognition of his fundamental contributions to
quantum chemistry.}}
\author{David Z. Goodson and Alexey V. Sergeev\\
Department of Chemistry and Biochemistry\\
University of Massachusetts Dartmouth\\
North Dartmouth, MA\ \ 02747, USA\\
}
\begin{document}
\maketitle
\begin{abstract}
M{\o}ller-Plesset perturbation theory expresses the energy as a function
$E(z)$ of a perturbation parameter, $z$. This function contains
singular points in the complex $z$-plane that affect the convergence
of the perturbation series. A review is given of what is known in
advance about the singularity structure of $E(z)$ from functional
analysis of the Schr\"odinger equation, and of techniques for empirically
analyzing the singularity structure using large-order perturbation series.
The physical significance of the singularities is discussed. They fall
into two classes, which behave differently in response to changes in
basis set or molecular geometry. One class consists of complex-conjugate
square-root branch points that connect the ground state to a low-lying
excited state. The other class consists of a critical point on the
negative real $z$-axis, corresponding to a dissociation phenomenon.
These two kinds of singularities are characterized and contrasted using
quadratic summation approximants. A new classification scheme for
M{\o}ller-Plesset perturbation series is proposed, based on the
relative positions in the $z$-plane of the two classes of singularities.
Possible applications of this singularity analysis to practical problems
in quantum chemistry are described.
\end{abstract}
\newpage
\tableofcontents
\section{Introduction}
The starting point for most {\it ab initio} calculations of molecular
energies is the Hartree-Fock (HF) approximation. Therefore, a central
problem of quantum chemistry is the construction of an extrapolation
from the HF energy to the true energy eigenvalue of the Schr\"odinger
equation. A particularly straightforward approach, at least in
principle, is a perturbation
theory proposed by M{\o}ller and Plesset \cite{mp} in which the
HF wavefunction is taken as the zeroth-order approximation for the
eigenfunction. The theory can be formulated by partitioning the
Hamiltonian according to \cite{cremerreview}
\begin{equation}
H=H_0+\left( H_{\rm phys}-H_0\right) z,
\label{MPpartition}
\end{equation}
where $H_0$ is the sum of one-electron Fock operators, $H_{\rm phys}$
is the Schr\"odinger Hamiltonian, and $z$ is a perturbation parameter.
The energy is then obtained as a power series
$E(z)=E_0+E_1z+E_2z^2+\cdots$. Thus, in M{\o}ller-Plesset (MP)
theory the energy is a function of $z$, in the complex
$z$-plane, such that $E(0)$ is equal to the
sum of HF orbital energies and $E(1)$ is the extrapolation to the
physical energy.
Traditionally, $E(z)$ is calculated by partial summation, that is,
the power series is truncated at some given order and then evaluated
at $z=1$. Truncation at order $n$ yields the ``MP$n$'' approximation
to the energy. Thus, $E(z)$ is evaluated as a polynomial.
The power series is an asymptotic series, which is a rigorously
correct solution only in the $z\to 0$ limit \cite{benderorszag}.
The true functional form of $E$ is much more complicated than
a polynomial. In
particular, it has a rich structure of singular points, that is, points
in an infinitesimal neighborhood of which the first derivative does
not exist. Since a polynomial is a functional form that is
nonsingular at any finite $z$, it cannot describe the true function
in the neighborhood of a singularity or outside the convergence
radius determined by the nearest singularity to the origin. This
is the cause of poor convergence seen in recent high-order MP$n$
studies \cite{dunningmp,olsen1,leininger}.
Knowledge of the singularity structure of $E(z)$
can be used to predict the convergence behavior of an MP
series \cite{dzghopt,mp4let,dzgijqc}.
The singularity positions in the MP energy function also
affect the convergence of the corresponding coupled-cluster (CC)
sequence \cite{dzgcc,zheng}. Since CC theory is another
method that extrapolates from the HF approximation to the exact
solution, this is not surprising. Although the function $E(z)$
corresponds to the ground state energy, its singularity structure
contains information about excited states. If the singularities are
not all reasonably far from the points $z=0$ and $z=1$, then this
indicates that a single-determinant HF wavefunction will not be a
suitable starting point for describing the true wavefunction.
Our purpose here is to review what is known in advance about the
functional form of $E(z)$ and then to describe methods for
obtaining detailed information about the singularities from analysis
of MP series. In Section~\ref{section_analysis} we present existence
theorems that tell us what singularity structure to expect. Then
in Section~\ref{section_FCI} we describe methods for ``empirical''
studies in which information about singularity structure is extracted
from high-order series obtained from full configuration-interaction (FCI)
calculations. In Section~\ref{section_examples} we analyze some
representative examples and clarify the connection between singularity
structure and convergence ``class.'' MP series
convergence is often described in terms of a classification scheme
developed by Schmidt {\it et al.} \cite{handyfeenberg}.
Series with monotonic convergence are put in class~A while series
in which the $E_i$ alternate in sign are put in class~B.
We suggest a modified classification scheme that
explains the behavior of molecules with triple bonds or with distorted
geometries, which have previously been considered
anomolous \cite{olsen1,leininger,cremerjpc}.
Finally, we discuss possible applications of
MP singularity analysis to practical problems in quantum chemistry.
\section{Functional Analysis of the Energy}
\label{section_analysis}
\subsection{Singularities of the Energy Function}
\label{section_energysingularities}
MP perturbation theory is generally used only for the ground state
energy. The zeroth-order approximation for the wavefunction is taken
to be the lowest-energy HF determinant. Nevertheless, the function
$E(z)$ obtained from this perturbation series, when considered as a
function over the complex $z$ plane, contains within it the full
energy spectrum of eigenstates with the same symmetry as the ground
state. This is the consequence of a theorem presented by
Katz \cite{katz}:
\medskip
\it
\noindent
Let $E^{(0)}$, $E^{(1)}$, $E^{(2)}$,\dots be the spectrum of energy
functions of Eq.~(\ref{MPpartition}) corresponding to a given set of all
conserved quantities (e.g., angular momentum). For any $i$ and $j$
there exist complex-conjugate pairs of branch points $z_{ij}$, $z_{ij}^*$
in the complex $z$ plane at which $E^{(i)}$ is equal to $E^{(j)}$.
In the neighborhood of~$z_{ij}$
\begin{equation}
E^{(i)}(z) = b_{ij}-c_{ij}(1-z/z_{ij})^{1/2},\quad
E^{(j)}(z) = b_{ij}+c_{ij}(1-z/z_{ij})^{1/2},
\label{katzeq}
\end{equation}
where $b_{ij}$ and $c_{ij}$ are constants, and similarly for $z_{ij}^*$.
\medskip
\rm
\noindent
Now consider the ground state function $E^{(0)}(z)$. Starting at
the origin, trace a path in the $z$ plane that circles about the
point $z_{01}$. This point is a branch point singularity, which
means that a 360$^\circ$ circuit will lead to a new Riemann sheet
corresponding to the function $E^{(1)}(z)$. Similarly, following
$E^{(0)}(z)$ on a path that circles $z_{02}$ leads to $E^{(2)}(z)$,
and so on.
This phenomenon was studied by Baker \cite{bakerrmp}
using a model two-body problem for which an exact solution for
$E(z)$ can be obtained. Depending on the choice of parameters, the
complex-conjugate branch points connecting $E^{(0)}$ and $E^{(1)}$ could
be in the negative half-plane or in the positive half-plane. If branch
cuts were specified from these branch points out to infinity, then there
were no other branch points. However, if the branch cut was drawn
perpendicular to the real axis, joining the two points, then some of
the other branch points predicted by Katz, connecting the ground state
with higher excited states, were found. The situation is shown
schematically in Fig.~1. If one follows $E^{(0)}(z)$ on a path from
the origin along the real axis, there will be an avoided
crossing as $z$ passes between the 0-1 branch points and no
additional avoided crossings with higher states. However, a path
that temporarily leaves the real axis to avoid the 0-1 branch points
will show no
interaction between $E^{(0)}$ and $E^{(1)}$ but will have an avoided
crossing with $E^{(2)}$. Since Katz's analysis does not depend on the
specific form of the potential energy function,
we expect this behavior to be generic.
\begin{figure}
\noindent
\quad\epsffile{singsf1b.ps}
\epsffile{singsf1d.ps}
\noindent
\quad\epsffile{singsf1a.ps}
\epsffile{singsf1c.ps}
\caption{\label{fig:fig1} Schematic representation of energy functions
$E^{(0)}$, $E^{(1)}$, and $E^{(2)}$ for two paths in the complex $z$-plane.
Branch cuts (dashed lines) are shown for the branch points $z_{01}$ and
$z_{01}^*$.}
\end{figure}
These square-root branch points are not the only kind of singularities
in the energy function. Baker \cite{bakerrmp} noted that a
many-fermion system in which there are both attractive and repulsive
terms in the potential will undergo a spontaneous autoionization at some
point on the real $z$ axis. He suggested this is analogous to a
critical point in an $(E,z)$ phase diagram. For $z$ beyond the
critical value $z_c$, there can be no eigenstate in which the system
is bound. $z_c$ is a branch-point singularity, although presumably of
a more complicated form than Eq.~(\ref{katzeq}).
This was studied by Stillinger \cite{stillinger} for the
specific case of an atom. He noted that a negative real $z$, in
Eq.~(\ref{MPpartition}), implies that the $r_{12}^{-1}$ interelectron
potential energy terms in $H_{\rm phys}$ are negative, and hence
electrons are mutually attractive. At the same time the repulsive
mean field central potential in $H_0$ is, on the negative real axis,
increased by a factor of $(1+|z|)$, which works against the nucleus-electron
attraction. For sufficiently large $|z|$ it will be the case that a
bound electron cluster infinitely separated from the nucleus will have
a lower energy than any state in which the nucleus is bound.
Thus, at $z_c$ the nucleus separates from the atom.
The effect of $z_c$ on the behavior of $E^{(0)}(z)$ on the
real axis will be rather different from that of the Katz singularities.
On a path from the origin along the negative real axis $E^{(0)}(z)$
acquires an imaginary part as it passes through the singularity. This
is because the eigenstate then corresponds to a state in the
scattering continuum. We expect that the derivative of the function
along this path will be continuous.\footnote{The discontinuity of the
derivative occurs on paths through $z_c$ that leave the real axis.}
This is due to the fact that when $z$ is only slightly
beyond the critical point there exists a long-range tunneling barrier
in the potential, which yields an exponentially small resonance width.
Similar behavior is seen in the $1/Z$ expansion, where
$Z$ is the nuclear charge, for the two-electron atom \cite{morganpra}.
\subsection{Convergence Classes}
\label{section_convergenceclass}
It is common practice to classify MP series in terms of their
convergence patterns. With class~A systems the $E_i$ series coefficients
are all of the same sign. The series show monotonic but often
rather slow convergence. Class~B series, except at the lowest orders,
have $E_i$ that alternate in sign. They often converge rapidly
but sometimes diverge. A physical interpretation of this classification
system was developed by Cremer and He \cite{cremerjpc}. They examined
the terms in the wavefunction that lead to the sign patterns in the $E_i$
and concluded that the class~B sign alternation is caused by orbital
crowding, as in systems with highly electronegative atoms and closed
shells, while class~A behavior results from uncrowded orbitals, as in
alkanes and radicals.\footnote{Cremer and coworkers \protect\cite{forsberg}
have recently proposed a new nomenclature in which the former classes A
and B are now referred to as type~I and II, respectively, based on the
orbital structure, and with class designations A through E describing
a variety of actual convergence patterns seen in practice.}
In Ref.~\cite{dzghopt} it was suggested that these two convergence
patterns could be explained in terms of singularity structure.
The sign alternation seen for class~B would result if the dominant
singularity lies on the negative real axis, while class~A would
result if the dominant singularity is on the positive real axis.
The former is consistent with the Stillinger critical point.
The Katz analysis predicts complex-conjugate
pairs of branch points displaced from the real axis. This
would seem to be inconsistent with the class~A
sign pattern, since a complex conjugate pair,
$z_{01}$, $z_{01}^*$ in the positive half plane
will give a periodic sign pattern with regions in which successive
$E_i$ are negative alternating with regions in which they are positive.
However, if the imaginary part of these singularities is small then the
period, given by \cite{morganpra,dzgdpt}
\begin{equation}
n_0=\pi/\arctan\left( |{\rm Im}\, z_{01}/{\rm Re}\, z_{01}|\right) ,
\end{equation}
will be sufficiently large that the series will appear to converge
monotonically up to very high orders.
A connection between singularity analysis and the Cremer-He analysis can
perhaps be made as follows. If the valence electron orbitals are
crowded into a relatively small region of space, then the attractive
interaction between the electrons (for negative $z$) is especially
strong, as is the mean-field repulsive perturbation. Both of these
effects favor autoionization at small $|z|$. It is not
surprising that of the class~B systems for which high-order MP
series have been calculated, the one with the dominant singularity
closest to the origin is F$^-$, with the aug-cc-pVDZ basis set.
It is interesting that the position of $z_c$
is strongly dependent on the basis set. In the case of F$^-$
if the cc-pVDZ basis, which lacks the diffuse functions of the augmented
basis sets, is used instead, then the singularity on the negative real axis
seems to disappear \cite{olsen1,dzghopt,olsenmodel}. Apparently, diffuse
functions are needed to describe the autoionization.
In contrast, the position of the dominant singularity for class~A systems
is relatively insensitive to the basis set. This is a Katz branch point
pair connecting the ground state with the first excited state. The
interaction between these two states can be described reasonably well
with a small basis set. The effect of additional polarization functions,
representing highly excited configurations, is small. However, changes
in molecular geometries can have a large effect on the class~A
singularity \cite{zheng}, with bond stretching causing the branch points
to move in toward the physical point $z=1$. The distance of these
branch points from $z=1$ is related to the energy gap between the
states at $z=1$. As a bond is stretched the energy gap decreases and
the point of avoided crossing along the positive real axis moves closer
to $z=1$.
It is important to note that a path along the {\it negative} real axis
could also show an avoided crossing between the ground state and
the first excited state. Below, in Section~\ref{section_complicated},
we present empirical evidence indicating that such a Katz
branch point pair, in the negative half plane, is responsible for
the worst examples of MP series divergence. Thus, there are two
distinct kinds of class~B systems. We suggest that
classic class~B behavior, with alternating signs in the $E_i$, is due
to a critical point while the complicated sign patterns
in systems such as the C$_2$ molecule
or in molecules with bonds stretched well beyond
the equilibrium distance are due to the simultaneous
effects of Katz degeneracies in
both the positive and negative half planes.
\section{Extracting Singularity Structure from
Perturbation Series}
\label{section_FCI}
There are two general approaches for obtaining information about singularity
structure from a perturbation series \cite{pearce}: asymptotic methods
and approximant methods. The former category includes such familiar
procedures as the ratio test and the $n$th-root test as well as some more
modern techniques \cite{ninham,hunter}, all of which are based on
a theorem of Darboux that states that in the limit of large order the
series coefficients of the function in question become equivalent to
the corresponding coefficients of the Taylor series of the {\it dominant}
singularity, that is, the singularity closest to the origin. The
advantage of these methods is that we have for them rigorous convergence
proofs. A disadvantage is that they only provide information about the
one singularity closest to the origin. More distant singularities can
sometimes be studied using conformal mappings \cite{pearce}, but, even
for the dominant singularity, convergence can be significantly slowed
by interference from other singularities.
For analysis of MP series we have found approximant methods to be more
useful than asymptotic methods.
The idea is to construct a model function (a
{\it summation approximant}) containing parameters that are chosen so
that its Taylor series agrees with the coefficients of the perturbation
series up to some given order. Although there are no generally
applicable convergence theorems for approximants,
if the functional form of the approximant is a reasonable model for
the true functional form then the convergence can be rapid.
Approximants can simultaneously fit more than one singularity.
Furthermore, they provide an analytical continuation to regions
of the complex plane where partial sums are poorly convergent or
divergent.
A systematic approach to constructing approximants was proposed by
Pad\'e \cite{pade,gravesmorris}. A set of polynomials $A^{(m)}_k (z)$,
where $k$ indicates the degree of the polynomial, is defined by the
asymptotic relation
\begin{equation}
\sum_{m=0}^M A^{m}_{k_m} E^m = {\cal O}(z^K),\quad
K=M+\sum_{m=0}^M k_m.
\label{padeapprox}
\end{equation}
The notation ${\cal O}(z^K)$ means that the coefficients of the
Taylor series of the left-hand side are equal to zero up through
order $K-1$. $E$ in Eq.~(\ref{padeapprox}) represents the perturbation
series of the energy function. Collecting terms according to the
power of $z$ yields $K$ linear equations for the $K+1$ coefficients of the
polynomials. The final coefficient is set according to an arbitrary
normalization condition, typically, $A^{(M)}(0)=1$. Once the polynomial
coefficients are obtained, Eq.~(\ref{padeapprox}) is solved for $E(z)$.
With $M=1$ we have rational approximants,
$E(z)=A^{(1)}(z)/A^{(0)}(z)$. The application of these to perturbation
theories of the Schr\"odinger equation was pioneered by
Brand\"as and Goscinski \cite{brandaesgoscinski}. They proposed using
approximants to determine singularity positions. However, the
singularity structure of rational approximants is not of the form
we expect for the MP energy. They are single-valued functions,
containing poles, at the zeros of $A^{(0)}$, but no branch points.
Given enough series
coefficients, they can model a branch cut with combinations of
poles and zeros \cite{dzghopt,bakerpade} but they cannot model more
than one branch.
More appropriate are quadratic
approximants \cite{dzghopt,bakerpade,shafer,jordan,popov,sergeev,quadvib},
with $M=2$, which have the double-valued solution
\begin{equation}
E={1\over 2}\left[ { A^{(1)}\over A^{(2)} } \pm
{ 1\over A^{(2)} }\sqrt{ (A^{(1)})^2-4A^{(0)}A^{(2)} } \right] ,
\label{quadapprox}
\end{equation}
with a square-root branch point as specified by the Katz theorem.
It has two branches ($\pm$) connected by branch points at the values
of $z$ at which the discriminant polynomial, $(A^{(1)})^2-4A^{(0)}A^{(2)}$,
is zero. In principle, one branch describes the ground state and the
other, the first excited state. More complicated branch points
are in practice fit by these approximants with clusters of square-root
branch points. Approximants with $M>2$ can describe additional branches
of the function, although the accuracy for higher
branches tends to be poor unless the series in known to quite high
order \cite{drazin,algebra,tipping,fernandez}. Differential approximants,
with the powers of $E$ in Eq.~(\ref{padeapprox}) replaced by derivatives,
can describe a wider variety of singularity types \cite{gravesmorris,khan}.
Other kinds of approximants can also be used. For example,
Olsen {\it et al.} \cite{olsenmodel} analyzed MP series convergence
using a $2\times 2$ matrix eigenvalue equation \cite{wilson,finley},
which implicitly incorporates a square-root branch point.
It is of course possible simply to explicitly construct an approximant
as an arbitrary function with the singularity structure that $E(z)$
is expected to have. We suggest, for example, approximants of the
form\footnote{This is an example of an
Hermite-Pad\'e approximant \protect\cite{bakerpade}.}
\begin{equation}
E(z)=
\pm P_A(z)
\left[ \left( 1-z/z_{01}\right)\left( 1-z/z_{01}^*\right)\right]^{1/2}
+P_B(z)B(z-z_c)+P_C(z),
\label{singfit}
\end{equation}
where $P_A$, $P_B$, and $P_C$ are polynomials and
\begin{equation}
B(z)=\int_0^\infty (1+zt)^{3/2} e^{-t}dt,
\label{Bfunction}
\end{equation}
The function $B(z-z_c)$, which can be expressed in terms of an
exponential integral, has the behavior we expect for autoionization
in that it has a branch point at $z_c$ at which
the function becomes complex but with the derivative continuous along
the real axis. The coefficients of the polynomials can be chosen
so that the Taylor series agrees with the perturbation series.
Note that the gap between the energy of the ground state and that of the
first excited state is given by
\begin{equation}
E^{(+)}-E^{(-)}
= 2 P_A(z)
\left[ \left( 1-z/z_{01}\right)\left( 1-z/z_{01}^*\right)\right]^{1/2}.
\label{gap}
\end{equation}
The value of the gap at $z=0$ is simply $2P_A(0)$. Since
$z=0$ corresponds to the Hartree-Fock solution, this value is easily
calculated. By thus constraining $P_A$, we can include in the analysis
additional information about the singularity structure beyond the
information we obtain from the perturbation series. The quadratic
approximant, Eq.~(\ref{quadapprox}), can be similarly constrained.
At present, the highest order for which direct calculations of MP series
have been carried out is only six \cite{cremerijqc}. However,
MP series to arbitrary order can be computed from elements of the
Hamiltonian matrix that are obtained in the course of an FCI
calculation \cite{bartletthopt,handyhopt}.
As a result, high-order MP series are now available in the literature
for a variety of systems with as many as 10 correlated valence
electrons \cite{olsen1,leininger}. In principle, the
singularity structure of $E(z)$ can be obtained from an FCI calculation
directly, without calculating the perturbation series, by
analyzing the characteristic polynomial as a function of $z$. This
is a polynomial of very high degree, which makes a full numerical
analysis of the $z$ dependence difficult. However, given an
initial estimate of a singularity position one could refine the
result with a local analysis of the characteristic polynomial in the
neighborhood of that point.
\section{Examples}
\label{section_examples}
\subsection{Class A}
Consider the molecule BH. Since the boron atom is four electrons short
of a full octet, this should be a clear example of class~A
according to the Cremer-He criteria \cite{cremerjpc}. The MP series has been
calculated \cite{olsen1,leininger}
through about 15th order for several of the correlation-consistent
basis sets \cite{ccbases}
and all the $E_i$ have the same sign.
\begin{table}
\caption{Singularity analysis of the energy function of the BH molecule
from quadratic approximants. The ``weight'' corresponds to $|c_{ij}|$
of Eq.~(\protect\ref{katzeq}).
\bigskip}
\begin{tabular}{|c|c|c|c|l|}
\hline
&Branch point in&&Branch point in&\\
&positive&&negative&\\
Basis set&half-plane&Weight&half-plane&Weight\\ \hline
cc-pVDZ& $1.65\pm 0.43 i$&0.4&$-3.6$&\ 0.001\\
cc-pVTZ& $1.60\pm 0.40 i$&0.3&$-3.5$&\ 0.01\\
cc-pVQZ& $1.60\pm 0.31 i$&0.4&$-3.5$&\ 0.04\\
aug-cc-pVDZ& $1.55\pm 0.41 i$&0.1&$-2.5$&\ 0.001\\
aug-cc-pVTZ& $1.59\pm 0.36 i$&0.4&$-3$&\ 0.01\\
aug-cc-pVQZ& $1.61\pm 0.35 i$&0.4&$-2.1$&\ 0.002\\
\hline
\end{tabular}
\label{tab:tab1}
\end{table}
In principle, the zeros of
the discriminant polynomial in Eq.~(\ref{quadapprox}) correspond
to branch points of $E(z)$. Since the number of zeros increases
rapidly with the order of the series, we expect most of these branch
points to be spurious, and in practice most are unstable from
order to order and are often quite distant from the origin or appear as
nearly coincident pairs. For BH, the branch points that are stable from
order to order are shown in Table~\ref{tab:tab1}. They are consistent with
the predictions of the functional analysis, with a complex-conjugate
pair $z_{01}$, $z_{01}^*$ with positive real part and a negative
real $z_c$. The former are the dominant singularities, as expected
for a class~A system. The $z_{01}$ values imply that the period of
the sign pattern of the $E_i$ is somewhere between 13 and 17. Since
the highest-order coefficients probably contain significant
roundoff error, these values are not inconsistent with the lack of
apparent sign alternation in the series.
The basis sets represent two families, cc-pVXZ and aug-cc-pVXZ. These
are equivalent except that the latter are augmented with diffuse
functions. The presence of diffuse functions seems to shift $z_c$
toward the origin but has little effect on $z_{01}$.
\begin{figure}
\epsffile{singsf2.ps}
\caption{Energy\protect\cite{energynote} in $E_{\rm h}$
for BH with the aug-cc-pVQZ basis as a function of
real $z$, from MP series of Leininger {\it et al.} \protect\cite{leininger}.
The solid curves are the two branches of the quadratic approximant
with polynomial degrees 5, 5, and 5, with the values corresponding to
the real part of the branch points marked by crosses. The dashed curve
is from the rational approximant with degrees 8 and 8. The filled
circle shows the physical point.}
\label{fig:z01BH}
\end{figure}
In Fig.~\ref{fig:z01BH} we plot the two branches of a quadratic
approximant as a function of $z$ on the real axis, for BH with
the largest basis set. A broad avoided crossing can be
seen at ${\rm Re}\, z_{01}$. Also shown is the result from a rational
approximant, which agrees very well with the quadratic approximant
for the lower branch up to the avoided crossing. It then jumps to
the higher branch after passing through the branch cut between
$z_{01}$ and $z_{01}^*$. Very similar behavior is seen with the
other five basis sets.
\subsection{Class B}
In Table~\ref{tab:tab1} we also list a singularity for BH on the negative
real axis, which presumably corresponds to the critical point. However,
the convergence is rather poor since it is much farther from the origin
than are the dominant singularities. For a better example we consider
the class~B system F$^-$. Christiansen {\it et al.} \cite{olsen1} found
that the MP series for F$^-$ with the aug-cc-pVDZ basis is divergent.
The only stable branch point from quadratic approximants is
at $-0.60$, which presumably corresponds to the critical point, $z_c$.
\begin{figure}
\epsffile{singsf3.ps}
\caption{Energy\protect\cite{energynote} in $E_{\rm h}$
for F$^-$ with the aug-cc-pVDZ basis as a function of
real $z$, from MP series of Ref.~\protect\cite{olsen1}.
The solid curve is the quadratic approximant with degrees
6, 6, and 6. When this approximant is complex, the real part
is shown as the dash-dot curve. The dashed curve
is the rational approximant with degree 9 for the numerator and
10 for the denominator.}
\label{fig:zcFm}
\end{figure}
Figure~\ref{fig:zcFm} shows $E(z)$ for negative real $z$. Up until
$z_c$ the rational and quadratic approximants are in agreement.
After passing through $z_c$ the quadratic approximant gains an
imaginary part. The rational approximant, unable to produce a branch
point, instead maps out a branch cut with alternating poles and
zeros \cite{bakerpade} along the real axis starting slightly beyond~$z_c$.
In Section~\ref{section_energysingularities} we argued that for a path
along the real axis the derivative at $z_c$ should be continuous. By
construction, the singularity from a quadratic approximant is a square-root
branch point, of the form $P(z_c)(1-z/z_c)^{1/2}$ where $P$ is nonsingular
at $z_c$. This expression
has a discontinuous derivative along any path.
However, the quadratic approximant in Fig.~\ref{fig:zcFm} restricts
the effects of this discontinuity to the immediate neighborhood of $z_c$
by making the weight of the singularity, $|P(z_c)|$, rather small, in
this case, 0.03. For this reason, the discontinuity of the derivative
is apparent only on close examination of the curve. Similar behavior
is seen for $z_c$ of BH, in which case, as shown in Table~\ref{tab:tab1},
the weight of $z_c$ is significantly smaller than the weight of $z_{01}$.
Additional evidence for a qualitative difference between the singularity
on the negative real axis and the complex-conjugate pair in the positive
half-plane is the fact that the quadratic approximants show more
clustering of spurious branch points around the former, suggesting
that isolated square-root branch points are insufficient to accurately
model the function.
\subsection{Systems with Complicated Singularity Structure}
\label{section_complicated}
Molecules such as N$_2$, C$_2$, and CN$^+$, and in general
molecules with bonds stretched from the equilibrium geometry,
do not fit well into the A/B classification
scheme \cite{olsen1,leininger,cremerjpc,olsenmodel,cremerijqc}.
They show irregularities in the sign patterns and poor convergence.
This has been attributed to the presence of
significant singularity structure in
both half-planes \cite{dzghopt,dzgijqc,dzgcc}.
For C$_2$, quadratic approximants show stable branch points at
$1.23\pm 0.35i$ and $-0.95\pm 0.33i$.
The behavior of $E(z)$ in the vicinity of these
points in shown in Fig.~\ref{fig:C2}. It is qualitatively similar to
the behavior at $z_{01}$ for BH. In the case of C$_2$ we see no
singularity on the negative real axis. This may be due to the fact
that the basis set (cc-pVDZ) lacks diffuse functions.
\begin{figure}
\epsffile{singsf4.ps}
\caption{Energy\protect\cite{energynote} in $E_{\rm h}$
of C$_2$ for real $z$ with the cc-pVDZ basis.
The MP series is from Ref.~\protect\cite{leininger}.
Solid curves are the branches of the quadratic approximant
with degrees 6, 6, and 6. Values corresponding to
the real part of the branch points are marked by crosses. The dashed curve
is the rational approximant with degrees 9 and 10. The filled
circle shows the physical point.}
\label{fig:C2}
\end{figure}
The anomolous systems have in common a relatively small weight for the
Hartree-Fock reference determinant in the $z=1$ wavefunction, which
suggests that Katz branch points will be especially significant.
Stretching a bond decreases the energy gap, which
shifts avoided crossings to lower values of $|z|$, causing branch point
positions to vary significantly along the potential energy surface.
Olsen {\it et al.} \cite{olsenmodel}
found that the negative real singularity for the hydrogen
fluoride molecule, determined from the large-order series, also
shifts toward the origin as the bond is stretched, but
the effect is small. An analysis \cite{zheng}
of MP4 series found a significant shift of singularities
toward $z=1$ in the positive half-plane
and toward $z=0$ in the negative half-plane
as the distances increase toward dissociation.
The MP4 analysis could not determine the imaginary part of
the singular points, and therefore could not distinguish the nature of
the singularity in the negative~half-plane.
\section{Discussion}
The MP energy function $E(z)$ is expected to have a
complex-conjugate pair of square-root branch points in the positive
half plane, a branch point on the negative real axis corresponding to
a critical point, and, possibly, a complex-conjugate pair of square-root
branch points in the negative half-plane. The complex-conjugate
branch points in either half-plane result from the same physical
phenomenon, avoided crossing between the ground state and a low-lying
excited state. They are insensitive to changes in the basis set but
shift significantly as bonds are stretched. In contrast, the critical
point represents an autoionization process and its position can depend
strongly on the basis.
The A/B classification scheme, based on series
coefficient sign patterns, is inadequate for describing this situation.
It seems more sensible to classify MP series according to singularity
structure. We suggest that the complex-conjugate branch points,
regardless of the real part, be referred to as ``class~$\alpha$''
singularities while the critical point be referred to as ``class~$\beta$.''
An atom or molecule in which the class $\beta$ singularity dominates
the perturbation series can be called a class~$\beta$ system.
If the class~$\alpha$ singularities in the positive half-plane
dominate, then we have a class~$\alpha$ system. Systems in which
the positive class~$\alpha$ singularities are approximately the same
distance from the origin as the other class~$\alpha$ or the class~$\beta$
singularities can be called, respectively, class~$\alpha\alpha$ or
$\alpha\beta$.
Since for most applications to systems of chemical interest it is not
feasible at present to compute MP series beyond 4th or 5th order,
methods will have to be developed for estimating singularity
positions from low-order series. For clear class~$\alpha$ or
class~$\beta$ systems, quadratic approximants at MP4 can
give a reasonably accurate result \cite{dzghopt}. Difficulties
arise with $\alpha\alpha$ and $\alpha\beta$ systems because the approximants
tend to model the multiple singularities with a single branch point
half way between them \cite{dzgijqc}. One approach for dealing with
this is to use a conformal mapping to shift the singular points so
that one will always be somewhat closer to the origin than the
other \cite{pearce}. This is the basic idea behind the \hbox{MP4-q$\lambda$}
method \cite{mp4let}, which makes quadratic approximants in
most cases a dependable
tool for summing the series. In contrast, conformal mappings are less
effective with partial summation, since then one must shift
all singularities {\it away} from the origin.
An important potential application to practical quantum chemistry
is to use singularity positions as a diagnostic criterion for
predicting the accuracy of a calculation. This has already
been implemented to some extent, using the \hbox{MP4-q$\lambda$}
singularity positions to identify cases in which $\alpha\alpha$ and
$\alpha\beta$ effects adversely affect energy accuracy. The
positions of singularities in the negative half-plane seems to
be an excellent diagnostic of the accuracy of the MP energy
summation \cite{mp4let,dzgijqc}. Coupled-cluster calculations with
the CCSD(T) method are relatively insensitive to the class~$\beta$
singularity but are affected by $\alpha\alpha$ singularity structure.
The \hbox{MP4-q$\lambda$} singularity analysis can apparently
distinguish to some extent between $\alpha\alpha$ and $\alpha\beta$
systems, and thereby serve as an accuracy diagnostic for
coupled-cluster methods \cite{dzgcc}.
Another application is to multireference perturbation
theory (MRPT). A key question in choosing the reference set
of eigenstates is the extent of interaction
between different states \cite{hose,zarrabian,freed,fuelscher}.
Positions of the complex-conjugate branch points
provide a measure of this \cite{warken}. One could
monitor singularity structure along a path on the potential
energy hypersurface, changing the reference set as needed to
eliminate harmful branch points. At present, software is readily
available only for 3rd-order MRPT \cite{mrpt3}, which is likely not
sufficient for reliable singularity analysis. A 4th-order
MRPT with a small reference set might benefit significantly from
a quadratic approximant analysis, both for summation and for
determining branch point~positions.
\begin{thebibliography}{99}
\bibitem{mp}
C. M{\o}ller and M. S. Plesset, Phys. Rev. {\bf 46}, 618 (1934).
\bibitem{cremerreview}
D. Cremer, in {\it Encyclopedia of Computational Chemistry},
P. v. R. Schleyer, {\it et al.} eds. (Wiley, New York, 1998), pp.~1706-1735.
\bibitem{benderorszag}
C. M. Bender and S. A. Orszag,
{\it Advanced Mathematical Methods for Scientists and Engineers}
(McGraw-Hill, New York, 1978).
\bibitem{dunningmp}
K. A. Petereson and T. H. Dunning, Jr., J. Phys. Chem. {\bf 99}, 3898 (1995);
T. H. Dunning, Jr. and K. A. Peterson, J. Chem. Phys. {\bf 108}, 4761 (1998).
\bibitem{olsen1}
O. Christiansen, J. Olsen, P. J{\o}rgensen, H. Koch, and P.-A. Malmqvist,
Chem. Phys. Lett. {\bf 261}, 369 (1996);
J. Olsen, O. Christiansen, H. Koch, and P. J{\o}rgensen, J. Chem. Phys.
{\bf 105}, 5082 (1996).
\bibitem{leininger}
M. L. Leininger, W. D. Allen, H. F. Schaefer III, and C. D. Sherrill,
J. Chem. Phys. {\bf 112}, 9213 (2000).
\bibitem{dzghopt}
D. Z. Goodson, J. Chem. Phys. {\bf 112}, 4901 (2000).
\bibitem{mp4let}
D. Z. Goodson, J. Chem. Phys. {\bf 113}, 6461 (2000).
\bibitem{dzgijqc}
D. Z. Goodson, Int. J. Quantum Chem. {\bf 92}, 35 (2003).
\bibitem{dzgcc}
D. Z. Goodson, J. Chem. Phys. {\bf 116}, 6948 (2002).
\bibitem{zheng}
D. Z. Goodson and M. Zheng, Chem. Phys. Lett. {\bf 365}, 396 (2002).
\bibitem{handyfeenberg}
C. Schmidt, M. Warken, and N. C. Handy, Chem. Phys. Lett. {\bf 211},
272 (1993).
\bibitem{cremerjpc}
D. Cremer and Z. He, J. Phys. Chem. {\bf 100}, 6173 (1996).
\bibitem{katz}
A. Katz, Nucl. Phys. {\bf 29}, 353 (1962).
\bibitem{bakerrmp}
G. A. Baker, Jr., Rev. Mod. Phys. {\bf 43}, 479 (1971).
\bibitem{stillinger}
F. H. Stillinger, J. Chem. Phys. {\bf 112}, 9711 (2000).
\bibitem{morganpra}
J. D. Baker, D. E. Freund, R. N. Hill, and J. D. Morgan III,
Phys. Rev. A {\bf 41}, 1247 (1990).
\bibitem{forsberg}
B. Forsberg, Z. He, Y. He, and D. Cremer, Int. J. Quantum Chem. {\bf 76},
306 (2000).
\bibitem{dzgdpt}
D. Z. Goodson, M. L\'opez-Cabrera, D. R. Herschbach, and J. D. Morgan III,
J. Chem. Phys. {\bf 97}, 8481 (1992).
\bibitem{olsenmodel}
J. Olsen, P. J{\o}rgensen, T. Helgaker, and O. Christiansen,
J. Chem. Phys. {\bf 112}, 9736 (2000).
\bibitem{pearce}
C. J. Pearce, Adv. Phys. {\bf 27}, 89 (1978); and references therein.
\bibitem{ninham}
D. W. Ninham, J. Math. Phys. {\bf 4}, 679 (1963).
\bibitem{hunter}
C. Hunter and B. Guerrieri, SIAM J. Appl. Math. {\bf 39}, 248 (1980).
\bibitem{pade}
H. Pad\'e, Ann. de l'Ecole Normale Sup. 3i\'eme S\'erie {\bf 9},
Suppl. 1 (1892).
\bibitem{gravesmorris}
G. A. Baker, Jr., and P. Graves-Morris, {\it Pad\'e Approximants}
(Cambridge Univ. Press, Cambridge, 1996); and references therein.
\bibitem{brandaesgoscinski}
E. Brand\"as and O. Goscinski, Phys. Rev. A {\bf 1}, 552 (1970);
Int. J. Quantum Chem. {\bf 4}, 5718 (1970).
\bibitem{bakerpade}
G. A. Baker, Jr., {\it The Essentials of Pad\'e Approximants}
(Academic, New York, 1975).
\bibitem{shafer}
R. E. Shafer, SIAM J. Math. Anal. {\bf 11}, 447 (1975).
\bibitem{jordan}
K. D. Jordan, Int. J. Quantum Chem. Symp. {\bf 9}, 325 (1975);
Chem. Phys. {\bf 9}, 199 (1975).
\bibitem{popov}
V. D. Va\u{\i}nberg, V. D. Mur, V. S. Popov, and A. V. Sergeev,
Pis'ma Zh. Eksp. Teor. Fiz. {\bf 44}, 9 (1986) [JETP Lett. {\bf 44},
9 (1986)].
\bibitem{sergeev}
A. V. Sergeev, J. Phys. A: Math. Gen. {\bf 28}, 4157 (1995).
\bibitem{quadvib}
D. Z. Goodson and A. V. Sergeev, J. Chem. Phys. {\bf 110}, 8205 (1999).
\bibitem{drazin}
P. G. Drazin and Y. Tourigny, SIAM J. Appl. Math. {\bf 56}, 1-18 (1996).
\bibitem{algebra}
A. V. Sergeev and D. Z. Goodson, J. Phys. A: Math. Gen. {\bf 31}, 4301 (1998).
\bibitem{tipping}
F. M. Fern\'andez and R. H. Tipping, J. Mol. Struct. Theochem {\bf 488},
157 (1999).
\bibitem{fernandez}
F. M. Fern\'andez, and C. G. Diaz, Eur. Phys. J. D {\bf 15}, 41 (2001);
C. G. Diaz and F. M. Fern\'andez, J. Mol. Struct. Theochem {\bf 541}, 31
(2001).
\bibitem{khan}
M. A. H. Khan, J. Comp. App. Math. {\bf 149}, 457 (2002).
\bibitem{wilson}
S. Wilson, K. Jankowski, and J. Paldus, Int. J. Quantum Chem. {\bf 28},
525 (1985).
\bibitem{finley}
J. P. Finley, R. K. Chaudhuri, and K. F. Freed, J. Chem. Phys. {\bf 103},
4990 (1995).
\bibitem{cremerijqc}
Z. He and D. Cremer, Int. J. Quantum Chem. {\bf 59}, 15 (1996);
31 (1996); 57 (1996); 71 (1996).
\bibitem{bartletthopt}
W. D. Laidig, G. Fitzgerald, and R. J. Bartlett, Chem. Phys. Lett.
{\bf 113}, 151 (1985).
\bibitem{handyhopt}
P. J. Knowles, K. Somasundram, N. C. Handy, and K. Hirao,
Chem. Phys. Lett. {\bf 113}, 8 (1985);
N. C. Handy, P. J. Knowles, and K. Somasundram, Theor. Chim. Acta
{\bf 68}, 87 (1985).
\bibitem{ccbases}
T. H. Dunning, Jr., J. Chem. Phys. {\bf 90}, 1007 (1989);
R. A. Kendall, T. H. Dunning, Jr., and R. J. Harrison,
{\it ibid.} {\bf 96}, 6769 (1992).
\bibitem{energynote}
The plotted energy here is
$E_0+[E(z)-E_0]/z\sim (E_0+E_1)+E_2z+E_3z^2+\cdots$.
\bibitem{hose}
G. Hose, J. Chem. Phys. {\bf 84}, 4505 (1986).
\bibitem{zarrabian}
S. Zarrabian and R. J. Bartlett, Chem. Phys. Lett. {\bf 153}, 133 (1988);
S. Zarrabian and J. Paldus, Int. J. Quantum Chem. {\bf 38}, 761 (1990).
\bibitem{freed}
J. P. Finley and K. F. Freed, J. Chem. Phys. {\bf 102}, 1306 (1995); and
references therein.
\bibitem{fuelscher}
J. Olsen and M. P. F\"ulscher, Chem. Phys. Lett. {\bf 326}, 225 (2000).
\bibitem{warken}
M. Warken, J. Chem. Phys. {\bf 103}, 5554 (1995).
\bibitem{mrpt3}
H.-J. Werner, Molec. Phys. {\bf 89}, 645 (1996).
\end{thebibliography}
\end{document}