% REVTeX 3.0
% For APS galley-mode format:
%\documentstyle[aps]{revtex}
%\documentstyle[eqsecnum,aps]{revtex}
%\documentstyle[eqsecnum,aps,twocolumn]{revtex}
 \documentstyle[aps,twocolumn]{revtex}
%\documentstyle[eqsecnum,aps,amsfonts]{revtex}
%
% For preprint format:
%\documentstyle[preprint,aps]{revtex}
%\documentstyle[preprint,eqsecnum,aps]{revtex}
%\tighten
%
\begin{document}
% \draft command makes pacs numbers print
\draft
\title{Improving the convergence and estimating the accuracy of\\
summation approximants of $1/D$ expansions for Coulombic systems}
% repeat the \author\address pair as needed
\author{Melchior O. Elout, David Z. Goodson,
Carl D. Elliston,\cite{columbia} Shi-Wei Huang,
and Alexei V. Sergeev\cite{vavilov}}
\address{Department of Chemistry, Southern Methodist University,
Dallas, Texas  75275}
\author{Deborah K. Watson}
\address{Department of Physics and Astronomy, University of
Oklahoma, Norman, Oklahoma  73019}
\date{\today}
\maketitle
\begin{abstract}
The convergence of large-order expansions in $\delta=1/D$, where
$D$ is the dimensionality of coordinate space, for energies
$E(\delta)$ of Coulomb systems is strongly affected by singularities
at $\delta=1$ and $\delta=0$.  Pad\'e-Borel approximants with
modifications that completely remove the singularities at
$\delta=1$ and remove the dominant singularty at $\delta=0$
are demonstrated.  A renormalization of the interelectron
repulsion is found to move the dominant singularity of the
Borel function $F(\delta)=\sum_jE'_j/j!$, where $E'_j$ are
the the expansion coefficients of the energy with singularity
structure removed at $\delta=1$, farther from the origin and
thereby accelerate summation convergence.  The ground-state
energies of He and H$_2^{\,+}$ are used as test cases.  The
new methods give significant improvement over previous
summation methods.  {\it Shifted} Borel summation using
$F_m(\delta)=\sum_jE'_j/\Gamma(j+1-m)$ is considered.  The
standard deviation of results calculated with different values
of the shift parameter $m$ is proposed as a measure of
summation accuracy.
\end{abstract}
% insert suggested PACS numbers in braces on next line
\pacs{02.30.Lt, 31.15.Md}

% body of paper here
\narrowtext

\section{Introduction}

Large-order perturbation theory in $1/D$, where $D$ is the
dimensionality of coordinate space,
was proposed by Mlodinow \cite{mlod1,mlod2} as an
alternative to variational methods and coupling-constant
perturbation theories for solving the Schr\"odinger equation.
This ``dimensional perturbation theory'' has
some appealing features.  The first-order theory gives
a semiquantitative description of electron correlation
\cite{merwe,goscinski,drh,jgl,gscf1,zphys,kais1,book1,book2}
and other subtle many-body effects \cite{kais2,kais3}
and it is relatively easy to calculate the perturbation
expansion to very high order, using recursion relations
\cite{vainberg,dunn}.  $1/D$ expansions have now
been calculated through 30th order for two-electron atoms \cite{hegs}
and 10th order for 3-electron atoms \cite{li}.  Large-order
expansions have also been calculated for the H atom in an
external field \cite{bender,stark,popovrev,twofields,germannlett}
and for the H$_2^{\, +}$ molecule \cite{popovrev,h2p,huang}.

The expansions for energies of Coulombic systems are, in
general, divergent, due to singularities in the energy function
$E(\delta)$, with $\delta=1/D$ treated as a continuous
complex variable.  One type of singularity is a
pole in $E(\delta)$ that results from the fact that for certain
integer values of $D$ the expectation value of the Coulomb
potential can diverge at particle coalescences
\cite{doren}.  This singularity can be accurately
characterized and removed either by subtraction \cite{dorencpl}
or by rescaling \cite{drh,rescaling}, making partial sums
of the expansion appear rapidly convergent at low orders.
A second characteristic type of singularity is a complicated
branch point or essential singularity at $\delta=0$, which
leads to divergence at higher orders \cite{hegs}.
This singularity can in principle be explained by an instanton
analysis \cite{lipatov,zinn,divergence}.  The zeroth-order limit
of the perturbation theory corresponds to localization at an
extremum of an effective potential that is the sum of the
3-dimensional Coulomb potential and a centrifugal potential
that comes from the dimensional continuation of the
kinetic energy operator.  For applications to bound
systems, one generally chooses an extremum that corresponds
to a stable minimum.  However, there typically exist sub-barrier
trajectories in imaginary time that allow for tunneling out of
the minimum.  This causes the zeroth-order limit to be
singular, with the energy expansion coefficients diverging as a
factorial.

Summation methods in which the Coulombic poles are removed
from $E(\delta)$ and then the remainder expansion is summed
with Pad\'e approximants or Pad\'e-Borel approximants
yield apparently convergent results at all orders
\cite{hegs,summation,hexs}.  Our purpose here is to refine
those methods to improve the rate of convergence.
Our strategy is to renormalize the Hamiltonian, by redefining
the dimensional continuation in such a way that the effect of the
singularity at the origin is mitigated, and to improve the
the techniques for modeling the singularities.
In addition, we develop a method for estimating the accuracy
of a given summation approximant.
He and H$_2^{\,+}$ are used as test cases.


\section{Summation approximants}
\label{sec:summation}

The asymptotic expansion in $1/D$ for the electronic energy of
an atom or molecule has the form
\begin{equation}
E(\delta)\sim \delta^2\sum_{j=0}^k E_j\delta^j,
\label{expansion0}
\end{equation}
where $\delta=1/D$.  There are various algorithms for
computing the expansion coefficients $E_j$, but the most
efficient, especially for problems with more than one
degree of freedom, seems to be the matrix method of
Dunn {\it et al.} \cite{dunn}.  The direct evaluation of
Eq.~(\ref{expansion0}) at $\delta=1/3$ using partial sums,
from truncation at given $k$, typically gives slow
convergence at low orders and rapid divergence at high order.

For ground states of Coulombic systems, $E$ has an
expansion about $D=1$ in the form \cite{doren}
\begin{equation}
E(\delta)\sim a_{-2}\, (D-1)^{-2}
+a_{-1}\, (D-1)^{-1}+\cdots .
\label{expansion1}
\end{equation}
For a one-electron atom with nuclear charge $Z$, the first
term in Eq.~(\ref{expansion1}) constitutes the exact solution,
with $a_{-2}=2Z^2$.  For two-electron atoms $a_{-2}$ is equal to
the ground-state energy of a two-electron Schr\"odinger equation
in which the Coulomb potentials have been replaced with delta
functions \cite{doren,herrstil}.  The exact solution for this
equation was obtained by Rosenthal \cite{rosenthal} for
arbitrary $Z$.

To remove the effects of the poles one can either rescale
the expansion \cite{drh,rescaling}, by multiplying it by
$(1-\delta)^2$, or explicitly subtract out the $\delta$ expansions
of the poles \cite{dorencpl,summation}.  If the residues are
known exactly then subtraction works best.  We consider
two approaches for the case where $a_{-2}$ is known but
$a_{-1}$ is not known: subtracting the second-order pole
and rescaling to remove the first-order pole,
\begin{mathletters}
\begin{eqnarray}
&&E(\delta) \sim \delta^2\left[
 a_{-2}(1-\delta)^{-2} + (1-\delta)^{-1}E'(\delta)\right],
\label{rescaledhybrid}
\\
&&E'(\delta)=\sum_{j=0}^k \delta^j\, E'_j,\quad
  E'_j=E_j-E_{j-1}-a_{-2},
\end{eqnarray}
\label{reshybeqs}
\end{mathletters}

\noindent
or subtraction of both poles,
\begin{mathletters}
\begin{eqnarray}
&&E(\delta) \sim \delta^2\!\left[
 a_{-2}(1-\delta)^{-2}\! + a_{-1}(1-\delta)^{-1}\!
 +E'(\delta)\right],\!
\label{hybrid}
\\
&&E'(\delta)=\sum_{j=0}^k \delta^j E'_j,\quad
  E'_j=E_j\! -(j+1)a_{-2}\! -a_{-1},
\end{eqnarray}
\label{hybrideqs}
\end{mathletters}

\noindent
with
\begin{equation}
a_{-1}=\lim_{\delta\to 1}\sum_{j=0}^k\delta^j\,
 \left[ E_j-(j+1)a_{-2}\right].
\label{hybridam1}
\end{equation}
Eq.~(\ref{hybridam1}) can be evaluated using Pad\'e summation
\cite{summation}.  The use of approximants based on
Eqs.~(\ref{hybrideqs}) has been called {\it hybrid}
summation \cite{dorencpl}, since it is the
combination of an expansion about $\delta=1$ and an expansion
about $\delta=0$.  We will refer to the approach in
Eqs.~(\ref{reshybeqs}) as {\it rescaled hybrid} summation.

It is possible to obtain convergent results by summing the
expansions for $E'$ with Pad\'e approximants or with Pad\'e-Borel
approximants.  The latter give a somewhat better convergence
rate.  Here we employ a more general technique \cite{zinn} that
will be useful not just for summing the energy expansion but also
for analyzing the singularity structure at $\delta=0$ and for
estimating summation accuracy.
Consider a function $F_m(\delta)$ with asymptotic
expansion
\begin{equation}
F_m(\delta)\sim\sum_{j=j_0}^k \delta^j F_{m,j},\quad
F_{m,j}=E'_j/\Gamma(j+1-m),
\label{Fexpansion}
\end{equation}
where $j_0$ is the smallest integer greater than or
equal to $m$.  This function is related to $E'(\delta)$ by the
transform
\begin{eqnarray}
&&E'(\delta)=\sum_{j=0}^{j_0-1} E'_j\delta^j
\nonumber\\
&&
 +\; e^{-(m-1)i\xi}\!\!\lim_{x\to\infty}
\int_0^{xe^{i\xi}}\!\!\!\!\!
x^{-m}F_m\big(\delta xe^{i\xi}\big)\exp\!\big( -xe^{i\xi}\big) dx,
\label{boreltransform}
\end{eqnarray}
with the angle $\xi$ chosen so that $F(\delta x e^{i\xi})$
is nonsingular along the path of integration and
$-\pi/2 <\xi <\pi/2$.  If $m=0$ then
Eq.~(\ref{boreltransform}) is the usual Borel transform.
For the function $F_m$ in Eq.~(\ref{boreltransform}) we will
substitute Pad\'e approximants of the expansion in
Eq.~(\ref{Fexpansion}).  For $m\ne 0$ we will refer to the
summation method as {\it shifted} Pad\'e-Borel summation.


\section{Renormalization} 
\label{sec:renormalize}

\subsection{Charge renormalization}
\label{sec:chargerenorm}

The standard approach to analytic continuation of the
Hamiltonian to arbitrary $D$ is to continue the kinetic energy
operator to $D$ dimensions while keeping the
physical, 3-dimensional, definition for the potential energy
operator \cite{herrstil,highJ}.  If the energy
is expressed in atomic units multiplied by $Z^2/D^2$, then
the effective potential for an $S$ state of the two-electron
atom at large $D$ is
\begin{equation}
V_{\rm eff}(r_1,r_2,\theta)=V_{\rm centr}+V_{\rm Coul},
\end{equation}
\begin{equation}
V_{\rm centr}=\left({1-6\delta+8\delta \over 8\sin^2\theta}
                       -{1\over 8}\right)
\left({1\over r_1^2}+{1\over r_2^2}\right),
\label{Vcentr}
\end{equation}
\begin{equation}
V_{\rm Coul}=-r_1^{-1}\! -r_2^{-1}\!
+\lambda\, (r_1^2+r_2^2-2r_1r_2\cos\theta)^{-1/2}\! .\!
\label{VCoul}
\end{equation}

\noindent
$V_{\rm centr}$ is the centrifugal potential from the $D$-dimensional
kinetic energy operator while $V_{\rm Coul}$ is the
3-dimensional Coulomb potential, with $\lambda = 1/Z$.
This effective potential has a stable global minimum, with
$r_1=r_2$, as long as $\lambda$ is less than a critical value
$\lambda_c\approx 0.8144$.  For H$^-$, with $\lambda=1$,
this symmetric extremum is a saddle point \cite{merwepr,dorenhyd},
and if this point is used for the zeroth order of the perturbation
theory, then the coefficients of the energy expansion will be
complex numbers.  Similarly, the symmetric large-$D$ extrema for
the H$_2^{\,+}$ and H$_2$ molecules also correspond to
saddle points \cite{merwepr,frantz,tan}.

A strategy for calculating $1/D$ expansions for systems that
are unstable in the large-$D$ limit was proposed in
Ref.~\cite{traynor}, using the fact that the definition of
the dimensional continuation of the Hamiltonian
is arbitrary as long as it gives the correct operator at
$D=3$.  If we renormalize the repulsive part of the
Coulomb potential, by replacing it with a continuous function of
$\delta$ such that the repulsion becomes weaker as $\delta$
approaches zero, it is possible to ensure that the symmetric
extremum will be a minimum.  This approach is closely
related to the renormalized perturbation theory developed by
Killingbeck for anharmonic oscillators \cite{killingbeck}.

Charge renormalization was applied to large-order perturbation
theory for H$^-$ in Ref.~\cite{watson}.  The parameter $\lambda$
in Eq.~(\ref{VCoul}) was replaced with the linear polynomial
\begin{equation}
\lambda(\delta)=\lambda_0+\lambda_1\delta,
\quad \lambda_1=3\, (Z^{-1}-\lambda_0)
\label{linearlambda}
\end{equation}
with $\lambda_0$ treated as a free parameter.  $\lambda_1$ is
defined so that $\lambda(1/3)=1/Z$.  (This modification of
$\lambda$ does not affect the scale factor in the energy
units.  We continue to use atomic units multiplied by
$Z^2/D^2$, where $Z$ is independent of $D$.)  The hybrid
expansion, Eq.~(\ref{hybrid}), with Pad\'e-Borel summation of
$E'$ gave convergent results for H$^-$ for any value of $\lambda_0$
less than $\lambda_c$.  The fastest convergence was for
$\lambda_0\approx (2Z)^{-1}$, which gave a more accurate result
at given order for the ground-state energy of H$^-$ than had 
been obtained \cite{hegs} for the ground-state energy of He
with unrenormalized $\lambda$.  This rapid convergence for
H$^-$ suggests that the use of a $\delta$-dependent
potential energy operator could be used to improve the
convergence even for systems that are stable at large $D$.

Figure~\ref{fig1} shows the effect on the convergence for He
of setting $\lambda_0=(2Z)^{-1}$.  Since there is a fair amount
of scatter in the convergence from order to order, we have
plotted cubic polynomial fits to more clearly show the
convergence patterns.  The dashed curves correspond to
hybrid approximants, based on Eq.~(\ref{hybrid}), while the solid
curves correspond to rescaled hybrid approximants, based on
Eq.~(\ref{rescaledhybrid}).  The expansion of
$E'$ was summed with unshifted Pad\'e-Borel approximants.
We define the ``accuracy'' of a
summation approximant $S_k$ at given order $k$ as
\begin{equation}
\alpha(k)=-\log_{10}|(S_k-E_{\rm exact})/E_{\rm exact}|,
\label{alphadefn}
\end{equation}
which is a continuous measure of the number of accurate digits.

Several characteristic trends are evident.  First, it is
clear that the rescaled hybrid summation is the better of
the two summation methods.  It is possible to
somewhat improve the accuracy of the unscaled hybrid approximants
by applying Shanks extrapolation, with some subjective
judgement, to the Pad\'e-approximant sequences for $a_{-1}$,
but the results remain less dependable than those from the
rescaled hybrid method.  The general shapes of the rescaled
hybrid curves are representative of those over a wide range
of $\lambda_0$.  They begin with an almost linear increase
in the number of accurate digits.  Then, upon reaching an
accuracy plateau of between approximately 6 to 8 digits the
rate of convergence slows significantly.  For $\lambda_0<0.8/Z$
there is a return to a more rapid convergence rate, beginning
at around 25th order.  Figure~\ref{fig2} shows the initial
slope of $\alpha(k)$ as a function of $\lambda_0$.
The increasing rate of convergence is balanced at very low
order by the fact that the zeroth-order accuracy decreases
for $\lambda_0$ far from the physical value.  (The zeroth-order
result is equal to the exact energy for $\lambda_0\approx 0.46$,
which incidently corresponds to a dip in Fig.~\ref{fig2}.)
Nevertheless, there is for the He atom a clear advantage to
choosing $\lambda_0$ in the general neighborhood of $0.5/Z$.

The functional form chosen for $\lambda(\delta)$ is
essentially arbitrary.  In addition to the linear polynomial,
Eq.~(\ref{linearlambda}), we have considered a quadratic polynomial,
\begin{equation}
\lambda(\delta)=\lambda_0+\lambda_1\delta+\lambda_2\delta^2,
\end{equation}
\begin{equation}
\lambda_1={9\over 2Z}-4\lambda_0-{1\over 2}\lambda(1),
\quad
\lambda_2=3\lambda_0+{3\over 2}\lambda(1)
-{9\over 2Z},
\end{equation}

\noindent
with $\lambda(1)$ and $\lambda_0=\lambda(0)$ treated as
free parameters.  This makes it possible to control the behavior
at both of the solvable limits, $\delta=0$ and $\delta=1$.
The position of the dominant singularities in the Borel
function depend only on the value of $\lambda_0$.  Therefore,
the rate of convergence of the summation approximants depends
less strongly on $\lambda(1)$ than on $\lambda(0)$.  We find
in practice that certain choices of $\lambda(1)$ do improve
the convergence at low and intermediate orders, especially
for H$^-$, for which the linear $\lambda(\delta)$ can imply a
very large value for $\lambda(1)$.


\subsection{Effect on singularity structure}
\label{sec:singularities}

Analysis of the expansion coefficients for the ground-state
energy of two-electron atoms has shown \cite{hegs} that the
dominant singularities of $F_0(\delta)$ (i.e., the singularities
nearest to the origin of the complex $\delta$ plane) are a
complex-conjugate pair of square-root branch points in the
negative half-plane.  Let $-\beta e^{\pm ib}$ be their
locations.  Then the coefficients of the asymptotic expansion
in $\delta$ of the expression
$G \left(\delta+\beta e^{ib}\right)^{1/2}
+\, G^*\left(\delta+\beta e^{-ib}\right)^{1/2}$,
where $G$ is a constant, will in the limit of large order be
equal to the coefficients $F_{0,j}$ of the asymptotic expansion
of $F_0$, Eq.~(\ref{Fexpansion}).  Thus, in the limit of large
$j$ we have
\begin{eqnarray}
&&F_{0,j} \sim 2\, |G|\, \beta^{1/2}(-1)^j
\nonumber\\
&&\quad\times\;
{\Gamma\left(j-{1\over 2}\right)
 \over\Gamma\left(-{1\over 2}\right)\Gamma(j+1)}\,\beta^{-j}
\cos\big[\, (j-{\scriptstyle {1\over 2}})b
       - \arg\, G\,\big].
\label{F0j}
\end{eqnarray}
Now compare this to the function
\begin{equation}
f(\delta)=
\gamma \left(\delta+\beta e^{ib}\right)^{-1}
+\,\gamma^* \left(\delta+\beta e^{-ib}\right)^{-1},
\label{Borelpoles}
\end{equation}
which has the expansion coefficients
\begin{equation}
f_j = 2\, |\gamma |\,\beta^{-1}(-1)^j\beta^{-j}
\cos\big[\, (j+1)b -\arg\gamma\,\big].
\end{equation}
If we set
$\gamma=2^{-1}\pi^{-1/2}\beta^{3/2}|G|e^{i(\arg G+3b /2)}$,
then $f_j\sim F_{3/2,j}$, which implies that the dominant
singularities of $F_{3/2}(\delta)$ will be first-order poles
at the same locations as the dominant branch points in
$F_0(\delta)$.

The effect on convergence of the renormalization
can be understood in terms of the effect on the locations of
the dominant singularities of $F_m(\delta)$.
A convenient method for determining these locations is to
examine the singularities of the Pad\'e approximants of
$F_m$.  Previously \cite{hegs} we used the quadratic Pad\'e
approximants of $F_0$ for this purpose, noting the positions
of the branch points nearest the origin.  Once the branch
points are approximately identified in this manner,
the result can be refined by noting the positions of the
corresponding poles in the linear Pad\'e approximants of
$F_{3/2}$, which stabilize somewhat faster.

We find that the positions of these singularities depend
only on the value of $\lambda_0$.  In general, there appear
to be 3 distinct complex-conjugate pairs of singularities
in the Borel functions.  For $\lambda_0$
close to $\lambda_c$ the dominant singularities are in the
positive half plane, but as $\lambda_0$ decreases this pair
rapidly moves away from the origin.  For $\lambda_0<0.8$
the dominant singularities are the pair in the negative
half plane.  Figure~\ref{fig3} shows the 
real and imaginary parts of this pair.  The steady movement
away from the origin as $\lambda_0$ decreases is presumably
related to the improvement in convergence.
In the limit $\lambda_0\to 0$ the singularities in the Borel
function must disappear, since $E'(\delta)$ is then nonsingular.
It is interesting to note that the mechanism for this
disappearance is for the prefactor multiplying the dominant
singularities to go to zero, rather than for the positions of
the singularities to move out to infinity.
The third pair of singularities has an imaginary
part of $\pm\pi$ with a small real part that depends weakly
on $\lambda_0$.  This is far enough from the origin that
these can be expected to have no significant effect on the
rate of convergence.


\subsection{Effect on roundoff error}

The recursive computation of the energy expansion coefficients
$E_k$ is somewhat unstable to roundoff error
\cite{dunn,hegs}, on account of their factorial divergence.
However, in the
hydrogenic limit $Z\to\infty$ (i.e. $\lambda\to 0$) the energy
expansion is convergent, with only a linear increase in
$E_k$ with $k$.  At small but nonzero $\lambda$
the perturbation theory is still divergent but the growth of
the expansion coefficients is nearly linear until very high
order.  This apparently comes about through massive cancelations
of intermediate quantities that grow factorially.  The result
is extreme loss of precision near the hydrogenic limit.

The number of significant figures in the $E_k$
decreases approximately linearly with $k$.  It was stated
previously \cite{hegs,hexs} that roundoff error
in the $E_k$ becomes manifest in the $S_k$ at about
the value of $k$ for which the precision of the $E_k$
intersects the accuracy of the $S_k$.  More careful analysis
reveals that this rule is generally valid but somewhat
conservative.  The actual propagation of error in the summation
approximants is complicated, with some approximants
much less sensitive than others and thereby retaining their
accuracy beyond the intersection.  Such approximants can be
identified empirically by adding to the
$E_k$ random error of a given magnitude.

For the renormalized calculations we find that the rate of
precision loss depends mainly on the value of the parameter
$\lambda_0$, increasing as $\lambda_0$ decreases.  We have
used quadruple-precision arithmetic (32 decimal digits) for
computing the $E_k$ in the calculations reported here,
so that the roundoff error is always less than the convergence
accuracy.


\subsection{Renormalization for diatomic molecules}

Within the Born-Oppenheimer approximation the Hamiltonian is
not homogeneous, since the internuclear distance $R$ is treated
as a parameter rather than as a dynamical variable.  Therefore,
transformation to atomic units multiplied by $D^{-2}$ replaces
$R$ with the $D$-dependent parameter $\tilde{R}=R/D^2$.  In a
recent study of H$_2^{\,+}$ it was shown that renormalization of
this $D$ dependence is a simple and effective way to improve the
convergence of the $1/D$ expansion \cite{huang}.  One
replaces $\tilde{R}$ with
\begin{equation}
\tilde{R}={2(3+\beta)\over (D-1)(D+\beta)}\, R,
\end{equation}
with $\beta$ treated as a free parameter.  (The factor of
$D-1$ is included to give a reasonable $D\to 1$ limit.)  An
important effect of this renormalization is to ensure that the
$D\to\infty$ limit corresponds to an effective potential with a
single global minimum.  Otherwise the molecule
undergoes a symmetry breaking transition to a double-well problem
as $R$ increases, which significantly degrades the convergence of
the large-order expansion.  For many-electron diatomic molecules a
renormalization of the interelectron repulsion, completely
analogous to the procedure for
atoms in Section \ref{sec:chargerenorm}, could be used to attenuate
the effect on convergence of the interelectron repulsion in addition
to the $R$ renormalization that prevents the symmetry breaking.


\section{Estimating the accuracy of the summation approximants}
\label{sec:estimate}

We have calculated the summation approximants $S_k$ using a range
of integer and half-integer values for the shift parameter
$m$ in the Borel function $F_m$.  The best results tend to
correspond to $m=0$, 1/2, 1, and 3/2, but which of
these four is best at given order appears to be random.
Therefore, to obtain an estimate of the accuracy of the
summation at given order $k$ we compute the standard deviation
$\sigma_k$ from the mean $\bar{S}_k$ of the four different $S_k$.

Figure~\ref{fig4} shows results for the He ground state energy
with $\lambda_0=0.5/Z$.  The dotted curve shows the
actual accuracy of the mean.  The solid curve shows
$-\log_{10}|\sigma_k/\bar{S}_k|$.  The largest
disagreements between the actual accuracy and this estimate
occurs at very high orders, where the mean tends to be
significantly more accurate than the individual approximants
and the estimate is too conservative.


\section{Darboux-Borel approximants}

The singularity structure of the Borel function $F(\delta)$ for
H$_2^{\,+}$ is similar to that for the two-electron
atom except that the dominant singularity is typically a
square-root branch point that lies on the negative
real axis \cite{h2p,divergence}.  Because this problem is
separable, in ellipsoidal coordinates, it is
straightforward to derive an exact solution for the
location $\delta_0$ of the branch point \cite{divergence}.
Since this singularity is the nearest to the origin in the
complex plane, it follows from Darboux's theorem \cite{darboux}
that the expression
$g(\delta)(1-\delta/\delta_0)^{1/2}+h(\delta)$,
where $g$ and $h$ are nonsingular for $|\delta|\le |\delta_0|$,
has the same expansion coefficients in the limit of
large order as the actual expansion of $F$.  This suggests the
use of Darboux approximants \cite{copenhagen},
\begin{equation}
F_{[L,N/M]}(\delta)=
{R_N(\delta)\over Q_M(\delta)}(1-\delta/\delta_0)^{1/2}
+{P_L(\delta)\over Q_M(\delta)},
\label{darbouxapprox}
\end{equation}
where $P_L$, $Q_M$, and $R_N$ are polynomials of degrees
$L$, $M$, and $N$, respectively.  These polynomials can be
determined by setting $F_{[L,M/N]}$ equal to the expansion
of $F$, expanding the square root, multiplying through by
$Q_M$, and then solving the resulting set of linear equations.

Figure~\ref{fig5} compares accuracies of
Darboux-Borel and Pad\'e-Borel approximants for the
ground-state energy of H$_2^{\,+}$.  These results
correspond to internuclear distance $R=2$, which is very
close to the equilibrium bond length.
(The renormalized $D$-dimensional Hamiltonian for this
problem is described in Ref.~\cite{huang}.)
There is noticeable improvement in accuracy from the
Darboux approximants.  This is seen as well for
other values of $R$.


\section{Discussion}
\label{sec:discussion}

Our central results are the following: 1)~Summation accuracy
depends significantly on removing the Coulombic poles at
the $D\to 1$ limit, and rescaling removes the first-order
pole more effectively than subtraction using an approximate
residue.  2)~Renormalization
of the $D$-dimensional Hamiltonian, by multiplying the
interelectron potential with a polynomial in $\delta$
can significantly improve the rate of convergence.
3)~Comparison of results from shifted Borel transforms,
given by Eqs.~(\ref{Fexpansion}) and (\ref{boreltransform}),
yields a standard deviation $\sigma_k$ that provides a
reliable estimate of the summation accuracy at given order $k$.
A reasonable summation procedure is to carry out a series
of inexpensive moderate-order calculations using a quadratic
$\lambda(\delta)$ over the full range of $\lambda_0$ values,
with some arbitrary value of $\lambda(1)$, in order to determine
an approximate range of $\lambda_0$ that can be expected to yield
the fastest convergence.  In the case of He, for example, a plot
such as in Fig.~\ref{fig2} would suggest $0.1\le\lambda_0\le 0.3$.
Then additional moderate-order calculations could be used to
determine an optimal pair of $\lambda_0$ and $\lambda(1)$
to be used in the expensive large-order calculation.
One could smooth the convergence of the
approximant sequence by using the mean $\bar{S}_k$ of various
shifted approximants, eliminating $\bar{S}_k$ that have large
$\sigma_k$.  For example, if $\sigma_{27}$ is greater than
$\sigma_{26}$, as in Fig.~\ref{fig4}, then $\bar{S}_{26}$
would be retained as the 27th-order result.

We find that the best two-parameter fit of summation accuracy
vs. order, for two-electron atoms and for H$_2^{\,+}$, is
$\alpha(k)=\alpha_0+\alpha_1k^{1/2}$.  However, there appears
to be a systematic deviation from this behavior, with a leveling
off of accuracy at intermediate orders followed by an increase
at higher orders, which led us to use cubic polynomials
for the fits in Fig.~\ref{fig1}.
For He the plateau is at approximately 8 significant figures.
This effect was noted in Ref.~\cite{hegs}, for
unrenormalized expansions, and was attributed to the need to
model more than just the dominant singularity in the Borel
function.  For He one obtains approximately 3 significant
figures by explicitly including the poles at $\delta=1$,
and apparently, 7 figures from modeling also the dominant
branch point in the Borel function.  This explanation is
supported by the behavior of the Darboux-Borel summation
for H$_2^{\,+}$ in Fig.~\ref{fig5}, which explicitly includes the
square-root branch point.  In that case there is rapid convergence
to almost 7 figures at 14th order but then it is not until
29th order that the convergence resumes.  The plateau for the
Pad\'e-Borel results is less distinct, presumably because
at intermediate orders the Pad\'e approximants are
still refining the description of the branch point while also
attempting to decipher other singularity structure.
By order 27 the Pad\'e-Borel accuracy is about equal to the
Darboux-Borel accuracy.

We have shown results only for ground states, but in
general we expect similar singularity structure, and hence
similar convergence behavior, for the lowest eigenstate of
any given symmetry.  With other excited states there are
branch points in the {\it energy} function $E(\delta)$
\cite{hexs}.  Although not the dominant
singularities, in practice they dominate the convergence
behavior until very high order.  For example, it has been
estimated \cite{hexs} that the asymptotic behavior due to the
singularity at the origin for the $1s2s$ $^1S$ state will not
assert itself until 88th order.  In that case Borel summation
is not very effective, and it is more appropriate to use a
summation method designed to model the apparent dominating
singularities.

Recent 10th-order calculations for the 3-electron atom \cite{li}
show convergence behavior almost identical to that for
2-electron atoms.  It remains to be seen if this will
hold as well for 4-electron atoms.  Loeser \cite{jgl} has
formulated the $1/D$ expansion for $S$ states of all atoms in the
periodic table.  However, with 4 electrons there are 10
independent internal coordinates at large $D$ but only 9 at
$D=3$.  Therefore, it is not clear that Loeser's theory will
converge to the physical solution when taken to large order.  Such
coordinate redundancies proliferate quadratically as the number
of electrons increases.  An alternative formulation that
avoids this situation is ``partial'' dimensional continuation,
applied recently in a study of triatomic molecular
rotation spectra \cite{suvernev}.  In that case one
body-fixed axis was allowed to rotate in $D$ dimensions
while the other axis was treated as rotating about the first
in only 3 dimensions.  Some sort of related procedure could,
in principle, give a $D$-dependent Hamiltonian
with the physical number of coordinates.

Probably a more serious obstacle to applying the $1/D$
expansion to many-electron systems is the rapid increase
in computational cost with increasing numbers of internal
coordinates \cite{dunn}.
With larger systems it will be necessary to introduce
approximations, such as separability assumptions \cite{gscf2}
or basis-set truncations.  The tradeoff will be between
a small systematic error, from a large-order inexact theory, and
a random, but possibly larger, error from truncating the exact
theory at lower order.

\acknowledgements
This work was supported by grants from the National Science
Foundation and the Robert. A. Welch Foundation.

\begin{references}

\bibitem[*]{columbia} Current address: Dept. of Physics,
Columbia University, 538 W. 120th St., New York, NY  10027.

\bibitem[\dagger]{vavilov} Permanent address: S. I. Vavilov State
Optical Institute, 199034 St. Petersburg, Russian Federation.

\bibitem{mlod1} L. D. Mlodinow, {\it The Large $N$ Expansion
in Quantum Mechanics}, Ph.D. Dissertation, University of
California, Berkeley (1981).

\bibitem{mlod2} L. D. Mlodinow, in {\it Progress in
Particle and Nuclear Physics}, vol. 8, edited by D. Wilkinson
(Pergamon Press, Oxford, 1982), pp. 387-399.

\bibitem{merwe} P. du T. van der Merwe, J. Chem. Phys.
{\bf 81}, 5976 (1984); {\bf 82}, 5293 (1985).

\bibitem{goscinski} O. Goscinski and V. Mujica, Int. J. Quantum
Chem. {\bf 29}, 897 (1986).

\bibitem{drh} D. R. Herschbach, J. Chem. Phys. {\bf 84}, 838
(1986).

\bibitem{jgl} J. G. Loeser, J. Chem. Phys. {\bf 86}, 5635 (1987).

\bibitem{gscf1} D. Z. Goodson and D. R. Herschbach, J. Chem.
Phys. {\bf 86}, 4997 (1987).

\bibitem{zphys} D. R. Herschbach, J. G. Loeser, and D. K. Watson,
Z. Phys. D {\bf 10}, 195 (1988).

\bibitem{kais1} S. Kais and D. R. Herschbach, J. Chem. Phys.
{\bf 100}, 4367 (1994).

\bibitem{book1} D. R. Herschbach, in {\it Dimensional
Scaling in Chemical Physics}, edited by D. R. Herschbach, J.
Avery, and O. Goscinski (Kluwer, Dordrecht, 1993), pp. 7-80;
and references therein.

\bibitem{book2} J. G. Loeser and D. R. Herschbach, in
{\it New Methods in Quantum Theory}, edited by C. A. Tsipis,
V. S. Popov, D. R. Herschbach, and J. S. Avery (Kluwer,
Dordrecht, 1996), pp. 1-33; and references therein.

\bibitem{kais2} P. Serra and S. Kais, Chem. Phys. Lett.
{\bf 260}, 302 (1996); Phys. Rev. A {\bf 55}, 238 (1997).

\bibitem{kais3} J. P. Neirotti, P. Serra, and S. Kais,
Phys. Rev. Lett., in press.

\bibitem{vainberg} V. M. Va\u\i nberg, V. D. Mur, V. S. Popov,
A. V. Sergeev, and A. V. Shcheblykin, Teor. Mat. Fiz. {\bf 74},
399 (1988) [Theo. Math. Phys. {\bf 74}, 269 (1988)].

\bibitem{dunn} M. Dunn, T. C. Germann, D. Z. Goodson, C. A. Traynor,
J. D. Morgan III, D. K. Watson, and D. R. Herschbach, J. Chem. Phys.
{\bf 101}, 5987 (1994).

\bibitem{hegs} D. Z. Goodson, M. L\'opez-Cabrera, D. R. Herschbach,
and J. D. Morgan III, J. Chem. Phys. {\bf 97}, 8481 (1992).

\bibitem{li} M. O. Elout and D. Z. Goodson, unpublished.

\bibitem{bender} C. M. Bender, L. D. Mlodinow, and N.
Papanicolaou, Phys. Rev. A {\bf 25}, 1305 (1982).

\bibitem{stark} V. M. 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{popovrev} V. D. Mur, V. S. Popov and A. V. Sergeev, Zh.
Eksp. Teor. Fiz. {\bf 97}, 32 (1990) [Sov. Phys. JETP {\bf 70},
16 (1990)]; and references therein.

\bibitem{twofields} V. M. Va\u\i nberg, V. S. Popov, and
A. V. Sergeev, Zh. Eksp. Teor. Fiz. {\bf 98}, 847 (1990)
[Sov. Phys. JETP {\bf 71}, 470 (1990)].

\bibitem{germannlett} T. C. Germann, D. R. Herschbach, M.
Dunn, and D. K. Watson, Phys. Rev. Lett. {\bf 74}, 658 (1995).

\bibitem{h2p} M. L\'opez-Cabrera, D. Z. Goodson, D. R.
Herschbach, and J. D. Morgan III, Phys. Rev. Lett. {\bf 68}, 1992
(1992).

\bibitem{huang} S.-W. Huang, D. Z. Goodson, M. L\'opez-Cabrera,
and T. C. Germann, submitted to Phys. Rev. A.

\bibitem{doren} D. J. Doren and D. R. Herschbach, Phys. Rev.
A {\bf 34}, 2654 (1986);  J. Chem. Phys. {\bf 87}, 433 (1987).

\bibitem{dorencpl} D. J. Doren and D. R. Herschbach, Chem.
Phys. Lett. {\bf 118}, 115 (1985).

\bibitem{rescaling} M. L\'opez-Cabrera, A. L. Tan, and J. G.
Loeser, J. Phys. Chem. {\bf 97}, 2467 (1993).

\bibitem{lipatov} L. N. Lipatov, Pis'ma Zh. Eksp. Teor. Fiz.
{\bf 25}, 116 (1977) [JETP Lett. {\bf 25}, 79 (1977)].

\bibitem{zinn} J. Zinn-Justin, Phys. Rep. {\bf 109}, 109 (1981).

\bibitem{divergence} V. S. Popov and A. V. Sergeev, Phys. Lett. A
{\bf 193}, 165 (1994); Zh. Eksp. Teor. Fiz. {\bf 105}, 568 (1994)
[JETP {\bf 78}, 303 (1994)].

\bibitem{summation} D. Z. Goodson and D. R. Herschbach, Phys.
Rev. A {\bf 46}, 5428 (1992).

\bibitem{hexs} D. Z. Goodson and D. K. Watson, Phys. Rev. A {\bf 48},
2668 (1993).

\bibitem{herrstil} D. R. Herrick and F. H. Stillinger,
Phys. Rev. A {\bf 11}, 42 (1975).

\bibitem{rosenthal} C. M. Rosenthal, J. Chem. Phys. {\bf 55},
2474 (1971).

\bibitem{highJ} M. Dunn and D. K. Watson, Ann. Phys. (NY)
{\bf 251}, 266, 319 (1996); Few-Body Systems {\bf 21}, 187 (1996).

\bibitem{merwepr} P. du T. van der Merwe, Phys. Rev. A {\bf 36},
3446 (1987).

\bibitem{dorenhyd} D. J. Doren and D. R. Herschbach, J. Phys.
Chem. {\bf 92}, 1816 (1988).

\bibitem{frantz} D. D. Frantz and D. R. Herschbach, Chem. Phys.
{\bf 126}, 59 (1988).

\bibitem{tan} A. L. Tan and J. G. Loeser, in
{\it Dimensional Scaling in Chemical Physics}, edited by D. R.
Hersch\-bach, J. Avery, and O. Goscinski (Kluwer Academic, Dordrecht,
1993), pp. 230-255.

\bibitem{traynor} C. A. Traynor and D. Z. Goodson, J. Phys. Chem.
{\bf 97}, 2464 (1993).

\bibitem{killingbeck} J. Killingbeck, J. Phys. A {\bf 14},
1005 (1981).

\bibitem{watson} D. K. Watson and D. Z. Goodson, Phys. Rev. A
{\bf 51}, R5 (1995).

\bibitem{baker} J. D. Baker, D. E. Freund, R. N. Hill, and J.
D. Morgan III, Phys. Rev. A {\bf 41}, 1247 (1990).

\bibitem{copenhagen} D. Z. Goodson and M. L\'opez-Cabrera, in
{\it Dimensional Scaling in Chemical Physics}, edited by D. R.
Hersch\-bach, J. Avery, and O. Goscinski (Kluwer Academic,
Dordrecht, 1993), pp. 275-314.

\bibitem{darboux} B. W. Ninham, J. Math. Phys. {\bf 4}, 679
(1963).

\bibitem{suvernev} A. A. Suvernev and D. Z. Goodson, J. Chem.
Phys. {\bf 107}, 4099 (1997).

\bibitem{gscf2} D. Z. Goodson, Phys. Rev. A {\bf 55}, 4155 (1997).

\end{references}

\begin{figure}
\caption{Cubic polynomial fits of the number of accurate digits
from summation approximants for the ground state energy of He,
as a function of order, with ``accuracy'' defined as
$-\log|\Delta E/E_{\rm exact}|$, where
$\Delta E=E_{\rm approx}-E_{\rm exact}$, with $E_{\rm exact}$
from Ref.~\protect\cite{baker}.  The dashed curves are
from unshifted hybrid Pad\'e-Borel approximants while the
solid curves are from unshifted rescaled hybrid Pad\'e-Borel
approximants.  Results are shown for $\lambda_0=1/2$
(unrenormalized) and for $\lambda_0=1/4$, as labeled,
$\lambda(\delta)=\lambda_0+\lambda_1\delta$.}
\label{fig1}
\end{figure}

\begin{figure}
\caption{The rate of convergence at low order for unshifted
rescaled hybrid Pad\'e-Borel approximants for the ground state
energy of He, as a function of the renormalization parameter
$\lambda_0$, using $\lambda(\delta)=\lambda_0+\lambda_1\delta$.
The rate of convergence is defined here as the parameter
$\alpha_1$ from the 3-parameter least-squares fit of the accuracy
$\alpha(k)=\alpha_0+\alpha_1k+\alpha_2k^2+\alpha_3k^3$
with $\alpha_0$ constrained to equal the accuracy of $S_0$}
\label{fig2}
\end{figure}

\begin{figure}
\caption{Real part (solid curve) and imaginary part (dashed
curve) of one of the dominant singularities in the Borel
functions of the two-electron isoelectronic series, as a
function of $\lambda_0$.  The other dominant
singularity is at the complex conjugate of this location.}
\label{fig3}
\end{figure}

\begin{figure}
\caption{Accuracy of rescaled hybrid Pad\'e-Borel approximants
for the ground state energy of He, as a function of order in
the perturbation expansion.  The points are determined using
Borel functions $F_m$ with the values of the shift parameter
as follows: $+$, $m=0$; $\Box$, $m=-1/2$; $\triangle$, $m=-1$;
$\Diamond$, $m=-3/2$.  The dotted curve shows
$-\log_{10}|(\bar{S}_k-E_{\rm exact})/E_{\rm exact}|$,
where $\bar{S}_k$ is the mean of the four summation
approximants at order $k$.  The solid curve shows
$-\log_{10}|\sigma_k/\bar{S}_k|$, where $\sigma_k$ is the
standard deviation from the mean of the four approximants
at order $k$.}
\label{fig4}
\end{figure}

\begin{figure}
\caption{Accuracy of rescaled hybrid Darboux-Borel
approximants ($\nabla$) and Pad\'e-Borel approximants ($*$)
for the ground state energy of H$_2^{\,+}$ with internuclear
distance $R=2$, as a function of order in the perturbation
expansion.  The Borel summation is unshifted.}
\label{fig5}
\end{figure}

\end{document}