# Enhanced One–Step Neville Algorithm with Access to the Convergence Rate

Ulrich D. Jentschura  Department of Physics and LAMOR, Missouri University of
Science and Technology, Rolla, Missouri 65409, USA  Ludovico T. Giorgini
Nordita, Royal Institute of Technology and Stockholm University, Stockholm 106
91, Sweden

###### Abstract

The recursive Neville algorithm allows one to calculate interpolating
functions recursively. Upon a judicious choice of the abscissas used for the
interpolation (and extrapolation), this algorithm leads to a method for
convergence acceleration. For example, one can use the Neville algorithm in
order to successively eliminate inverse powers of the upper limit of the
summation from the partial sums of a given, slowly convergent input series.
Here, we show that, for a particular choice of the abscissas used for the
extrapolation, one can replace the recursive Neville scheme by a simple one-
step transformation, while also obtaining access to subleading terms for the
transformed series after convergence acceleration. The matrix-based, unified
formulas allow one to estimate the rate of convergence of the partial sums of
the input series to their limit. In particular, Bethe logarithms for hydrogen
are calculated to 100 decimal digits.

###### Contents

  1. I Introduction
  2. II Enhanced Neville Algorithm
    1. II.1 Formulation of the Algorithm
    2. II.2 Neville Algorithm
    3. II.3 Universal Formula
  3. III Comparison to the Other Methods
  4. IV Numerical Examples
    1. IV.1 Slowly Convergent Model Series
    2. IV.2 Bethe Logarithm
  5. V Conclusions
  6. VI Acknowledgements
  7. A Coefficients for Higher–Order Terms
  8. B Bethe Logarithms for Excited States

##  I Introduction

Often in physics, one faces problems connected with slowly convergent,
nonalternating series. The convergence acceleration of nonalternating, slowly
convergent sums is notoriously problematic. Examples include angular momentum
expansions in both nonrelativistic as well as relativistic atomic physics
calculations JeMoSo1999 ; JeMoSoWe1999 , series representations of Bethe
logarithms (Chap. 4 of Ref. JeAd2022book ), or, logarithmic sums over
eigenvalues required for the calculation of the functional determinant in
$O(N)$ theories GiEtAl2020 ; GiEtAl2022 ; GiEtAl2024i . Nonalternating, slowly
convergent series belong to a class series which are notoriously difficult to
accelerate We1989 ; BrRZ1991 ; CaEtAl2007 . Traditional methods like the
Aitken $\Delta^{2}$ process Ai1926 , the Shanks transformation Sh1955 , and
Wynn’s epsilon algorithm Wy1956eps  (which calculates Padé approximants), and
other methods, all encounter numerical instabilities when applied to
nonalternating slowly convergent series.

In Refs. NP1961 ; JeMoSoWe1999 ; AkSaJeBeSoMo2003 , an algorithm was described
which first converts the nonalternating input series to an alternating series,
whose convergence is accelerated via nonlinear sequence transformations We1989
; JeMoSoWe1999 ; AkSaJeBeSoMo2003 . In fact, the first step of the two-step
condensation, as described in Ref. JeMoSoWe1999 , consists of the van
Wijngaarden transformation NP1961 , which rearranges the nonalternating input
series into an alternating series, but it requires the calculation of
individual terms of the input series of very high order in the summation
indices. It is followed by a nonlinear sequence transformations Si1979 ;
Si1980 ; Si1986 ; We1989 ; We1996c ; Si2003 , leading to what is known as the
combined nonlinear-condensation transformation. Sometimes, however, the option
of the calculation of terms of very high order in the summation indices is not
available. This problem could occur in various contexts, e.g., because, on
lattices, we cannot calculate eigenvalues with a high principal quantum number
GiEtAl2020 ; GiEtAl2022 ; GiEtAl2024i .

For the notoriously problematic case of a nonalternating, slowly convergent
input series, which cannot be practically converted into an alternating series
by the van Wijngaarden rearrangement transformation NP1961 , established
methods like Aitken $\Delta^{2}$ process Ai1926  and Wynn’s Epsilon Algorithm
Wy1956eps  are often used. One notes, according to Sec. 2.2.7 of Ref.
CaEtAl2007 , that the Aitken $\Delta^{2}$ process and Wynn’s Epsilon Algorithm
are special cases of the general E-algorithm Br1980E ; Ha1979E . Of course,
the Aitken $\Delta^{2}$ process and Wynn’s Epsilon Algorithm constitute viable
and widely used convergence acceleration methods, but numerical instability
can be encountered in higher orders of the transformation.

In view of the numerical instabilities encountered in the convergence
acceleration of non-alternating series, it is desirable to use asymptotic
information about the input series to the maximum extent possible, in order to
achieve the maximum convergence acceleration already in low orders of the
transformation. One possibility to achieve this rapid convergence is to use
the Neville algorithm, which, in its most general form, constitutes an
interpolation algorithm Ne1934 . It can be used for convergence acceleration,
which is tantamount to extrapolation, if one uses the algorithm in order to
calculate the interpolating polynomial outside of the domain of interpolation.
One first interprets the partial sums of the input series as values of a
function at abscissas which tend to a limiting value (say, zero) as the number
of terms of the input series is increased. One possibility is to use
$x=1/(i+1)$ for the partial sums $s_{i}$ of the input series
($i=0,1,2,\dots$). One then calculates the value of the interpolating function
at the limiting value of the abscissa (in our example case, at $x=0$) which is
outside of the domain $x\in(1/(n+1),1)$ used for the interpolation (here, $n$
is the maximum summation index used for the input series in a given order $n$
of the transformation). The Neville algorithm uses information from all
partial sums $(s_{0},\dots,s_{n})$ and hence, about the asymptotic structure
of the input series, in order to construct the interpolating polynomial in $x$
and can therefore lead to highly efficient convergence acceleration. Its usual
formulation is based on a recursive lozenge scheme which uses a three-term
recurrence relation. Here, we aim to enhance the Neville algorithm by
replacing the three-term recurrence relation with unifying analytic formulas
which also give access to subleading asymptotics of the behavior of the
partial sums $s_{n}$ of the input series for large $n$.

This paper is organized as follows. In Sec. II, we describe out formulation of
the enhance Neville algorithm. The idea behind our formulation is discussed in
Sec. II.1, while a comparison to the Neville algorithm is presented in Sec.
II.2. We find universal formulas which allow access to the subleading terms
for high $n$, in Sec. II.3. The performance of the algorithm is compared to
other established methods in Sec. III. Numerical examples are discussed in
Sec. IV, and conclusions are reserved for Sec. V. Appendix A lists some
higher-order coefficients for the subleading asymptotic behavior of the input
series, and Appendix B gives numerical values for hydrogen Bethe logarithms
obtained using our method.

##  II Enhanced Neville Algorithm

###  II.1 Formulation of the Algorithm

One starts with the following relations for the partial sums $s_{n}$ of the
input series, and the remainder terms $r_{n}$,

$s_{n}=\sum_{k=0}^{n}a_{k}=s_{\infty}+r_{n}\,,\qquad
r_{n}=-\sum_{k=n+1}^{\infty}a_{k}\,.$ (1)

(Here, we start the terms from index $k=0$, by convention We1989 .)

Let us assume that the terms of the infinite series possess an asymptotic
expansion in inverse integer powers of the summation index,

$a_{k}=\frac{A}{k^{2}}+\frac{B}{k^{3}}+\mathcal{O}(k^{-4})\,.$ (2)

The idea is that, if $n$ is sufficiently large, we can use the asymptotic
approximation for the remainder term,

$r_{n}=-\sum_{k=n+1}^{\infty}a_{k}=-\sum_{k=n+1}^{\infty}\left[\frac{A}{k^{2}}+\frac{B}{k^{3}}+\mathcal{O}(k^{-4})\right]\,.$
(3)

One can use the relation

$\sum_{k=n+1}^{\infty}\frac{1}{k^{a}}=\frac{(-1)^{a}}{(a-1)!}\,\left.\frac{\partial^{a-1}\psi(z)}{\partial
z^{a-1}}\right|_{z=n+1}\,,$ (4)

where $\psi(z)$ is the logarithmic derivative of the Gamma function, to expand
the remainder term $r_{n}$ in inverse power of $n$. The first terms in this
expansion read as follows,

$r_{n}=\frac{A}{n}+\frac{B-A}{2n^{2}}+\frac{A-3B+2C}{6n^{3}}+\mathcal{O}(n^{-4})\,.$
(5)

One then tries to eliminate the asymptotic terms of the remainder function,
which are powers of $1/n$. To this end, one interprets the first $(n+1)$
partial sums $s_{i}$ with $i\in(0,\ldots,n)$ as approximations to the values
of a function $f(x)$ at the abscissa values $x=x_{i}=1/i$,

$s_{i}\approx f(x=x_{i}=1/(i+1))\,,\qquad i\in(0,\ldots,n)\,,$ (6)

so that $\lim_{\epsilon\to 0^{+}}f(\epsilon)=s_{\infty}$. If we assume that
$f(x)$ can be expanded as a power series, then we can write

$f(x_{i})=\sum_{j=0}^{\infty}c_{j}\,x_{i}^{j}\approx\sum_{j=0}^{n}c_{j}\,x_{i}^{j}=c_{0}+\frac{c_{1}}{(i+1)}+\frac{c_{2}}{(i+1)^{2}}+\ldots+\frac{c_{n}}{(i+1)^{n}}\,,$
(7)

with $n+1$ unknown coefficients $c_{j}$ ($j=0,\dots,n$). Given $n+1$ partial
sums $s_{i}$ of the input series $\\{s_{i}\\}_{n=0}^{\infty}$, the idea of the
algorithm is to solve the system of equations

$s_{i}=\sum_{j=0}^{n}c_{j}(n)\,\left(\frac{1}{i+1}\right)^{j}\,,\qquad
i,j\in(0,\ldots,n)\,,$ (8)

for the coefficients $c_{j}(n)$. One then assumes that the value of the
infinite series is recovered as $n$ is increased,

$s=s_{\infty}=\lim_{n\to\infty}c_{0}(n)\,,$ (9)

where $c_{0}(n)$ is the zeroth-order coefficient of the interpolating
polynomial $\sum_{j=0}^{n}c_{j}(n)\,x^{j}$, which, at the abscissas
$x=x_{i}=1/(i+1)$, coincides with the partial sums $s_{i}$ of the input
series. In view of the fact that $c_{0}(n)$ is identical to the value of the
interpolating polynomial at argument $x=0$ while all abscissas $x_{i}$ are
positive, the latter step constitutes an extrapolation which can lead to
convergence acceleration CaEtAl2007 .

###  II.2 Neville Algorithm

For the calculation of the coefficient $c_{0}(n)$ of the interpolating
polynomial, one can use the Neville algorithm Ne1934 . Given $n+1$ partial
sums $s_{i}$ of the input series $\\{s_{i}\\}_{n=0}^{\infty}$, the Neville
algorithm refines the approximation of $s_{\infty}$ by iteratively
interpolating between points defined by the sequence $s_{i}$. Recursive
schemes have been described in Refs. GaGu1974 ; vT1994 . Adapted to our choice
of abscissas $x_{i}=1/(i+1)$, one may formulate a lozenge scheme as follows.
One starts the recursion with the values

$s^{0}_{i}=s_{i},\qquad i\in(0,\ldots,n)\,,$ (10)

and uses a three-term recursion,

$s^{m}_{i}=\frac{(i+1)s^{m-1}_{i}-(i+1-m)\,s^{m-1}_{i-1}}{m},\qquad 1\leq
m\leq n\,,\qquad m\leq i\leq n\,.$ (11)

One increases $m$, thereby decreasing the allowable values of $i$ in the
process, until, for $m=n$, only one possible value is left for $i$, namely,
$i=m=n$. Finally, the desired coefficient $c_{0}(n)$ is recovered as follows,

$c_{0}(n)=s^{n}_{n}\,.$ (12)

While the recursive scheme is very helpful, it does not establish a direct
connection between each term of the partial sums used for extrapolation and
the final estimated limit. Hence, no additional information can be derived for
the asymptotic behavior of the partial sums $s_{n}$ for large $n$. In
particular, while the recursive scheme enables the estimation of the limit
$\lim_{x_{i}\to 0^{+}}f(x_{i})=s_{\infty}$, it does not provide insights into
the behavior of the interpolating function $f(x_{i})$ near $x_{i}=0$, e.g.,
the rate at which the limit is approached.

Let us briefly clarify a notation used in the following. Our goal is to
estimate $s_{\infty}$ utilizing the initial $n$ terms of the series. We denote
the estimate of the limit after the convergence acceleration transformation by
the symbol $\mathcal{T}_{n}$. Let now $\mathcal{T}$ be the transformation
which takes as input the partial sums $(s_{0},s_{1},\ldots,s_{n})$, and by
which one obtains $\mathcal{T}_{n}$,

$\mathcal{T}_{n}=\mathcal{T}(s_{0},s_{1},\ldots,s_{n}).$ (13)

For example, if $\mathcal{T}$ is the Neville algorithm as given in the
adaptation (11), then $\mathcal{T}_{n}=c_{0}(n)=s_{n}^{n}$. For the other
methods discussed here, the definition of $\mathcal{T}_{n}$ is adapted in the
obvious way.

###  II.3 Universal Formula

In order to obtain universal formulas which allow access to subleading terms
(that serve to measure the rate of convergence of the input series), it is
advantageous to map the problem onto a system of linear equations. One
interprets the coefficients $c_{j}(n)$ of the polynomial as the elements of a
vector $\vec{c}(n)=(c_{0}(n),...,c_{n}(n))$ with $n+1$ elements,

$s_{i}=\sum_{j=0}^{n}c_{j}(n)\,\left(\frac{1}{i+1}\right)^{j}=\sum_{j=0}^{n}M_{ij}(n)\,c_{j}(n)\,,\qquad
M_{ij}(n)=\left(\frac{1}{i+1}\right)^{j}\,,\qquad i,j\in(0,\ldots,n)\,.$ (14)

The elements $M_{ij}(n)$ can be interpreted as components of a matrix of
dimension $(n+1)\times(n+1)$, in which case one has
$\vec{s}=\mathbbm{M}(n)\cdot\vec{c}(n)$, where $\vec{s}=(s_{0},...,s_{n})$ is
the vector of the first $n+1$ partial sums of the input series. Therefore, the
solution can be obtained by matrix inversion,

$\vec{c}(n)=[\mathbbm{M}(n)]^{-1}\cdot\vec{s}\,,\qquad\vec{c}=(c_{0}(n),...,c_{n}(n))\,.$
(15)

The inversion of the matrix $\mathbbm{M}(n)$ can lead to significant numerical
instability and loss of accuracy. This is because large numbers can cause
overflow or underflow issues in floating-point arithmetic, while the entries
of $\mathbbm{M}(n)$ span many orders of magnitude, which leads to unfavorable
condition numbers for higher $n$. This can exacerbate rounding errors. For
this reason, it is advantageous to derive analytic formulas for the first five
rows of the inverted matrix $[\mathbbm{M}(n)]^{-1}$, which leads to analytic
formulas for the coefficients $c_{0}(n)$, $c_{1}(n)$, $c_{2}(n)$, $c_{3}(n)$,
and $c_{4}(n)$. We find that

$\displaystyle c_{j}(n)=$
$\displaystyle\;\sum_{i=0}^{n}\frac{(-1)^{n-i}\,(i+1)^{n}}{\Gamma(i+1)\,\Gamma(n-i+1)}\,\mathcal{P}_{j}(i,n)\,s_{i}\,,\qquad
j\in(0,1,2,3,4)\,,$ (16a) $\displaystyle\mathcal{P}_{0}(i,n)=$
$\displaystyle\;1\,,$ (16b) $\displaystyle 2\,\mathcal{P}_{1}(i,n)=$
$\displaystyle\;2i-3n-n^{2}\,,$ (16c) $\displaystyle 24\mathcal{P}_{2}(i,n)=$
$\displaystyle\;24i+24i^{2}-26n-36in+9n^{2}-12in^{2}+14n^{3}+3n^{4}\,,$ (16d)
$\displaystyle 48\mathcal{P}_{3}(i,n)=$
$\displaystyle\;48i+96i^{2}+48i^{3}-48n-124in-72i^{2}n+26n^{2}$
$\displaystyle\;-6in^{2}-24i^{2}n^{2}+29n^{3}+28in^{3}-n^{4}+6in^{4}-5n^{5}-n^{6}\,.$
(16e) $\displaystyle 5760\mathcal{P}_{4}(i,n)=$
$\displaystyle\;5760i+17280i^{2}+17280i^{3}+5760i^{4}-5712n-20640in-23520i^{2}n-8640i^{3}n$
$\displaystyle\;+3380n^{2}+2400in^{2}-3600i^{2}n^{2}-2880i^{3}n^{2}+3660n^{3}+6840in^{3}+3360i^{2}n^{3}$
$\displaystyle\;-385n^{4}+600in^{4}+720i^{2}n^{4}-888n^{5}-600in^{5}-130n^{6}-120in^{6}+60n^{7}+15n^{8}\,.$
(16f)

Results for the coefficients $c_{j}(n)$ with $j>4$ are given in Appendix A.
The results given in Eq. (16) facilitate the estimation of the functional
relationship of the $f(x_{i})$ for small abscissas $x_{i}$, thereby providing
insight the convergence rate of the partial sums of the input series toward
its limit. Specifically, one notes the relationship

$\displaystyle s_{n}=$
$\displaystyle\;s_{\infty}-r_{n}\approx\frac{c_{1}}{n+1}+\frac{c_{2}}{(n+1)^{2}}+\dots\,,\qquad
n\to\infty\,,\qquad c_{j}\equiv\lim_{n\to\infty}c_{j}(n)\,,$ (17)
$\displaystyle r_{n}=$
$\displaystyle\;s_{\infty}-s_{n}\approx-\frac{c_{1}}{n+1}-\frac{c_{2}}{(n+1)^{2}}-\dots\,,\qquad
n\to\infty\,,$ (18)

where we are pedantic with respect to the sign of the remainder term $r_{n}$.

Figure 1:  We investigate the convergence of the enhanced Neville
transformation (green curve) for the input series (26), as measured by the
quantity $\chi(n)$ defined in Eq. (27). The quantity $\chi(n)$ roughly
decreases by unity with every iteration of the algorithm, indicating the gain
of roughly one converged decimal with every order of the harvest
transformation. The convergence of the enhanced Neville transformation is
compared with the convergence of the Wynn epsilon algorithm (blue curves), of
the Aitken $\Delta^{2}$ process (red curves), and of the series itself (black
curves). Figure 2:  The figures illustrates the convergence of sequences
$d_{j}(n)$ toward $c_{j}$ for $j=1,2,3,4$, according to Eq. (31). The cases
$j=1,2,3,4$ are treated in Fig. (a), (b), (d), and (d), respectively. The
first 150 $d_{j}(n)$ with $n=0,\dots,149$ are compared to the value
$c_{j}=d_{j}(\infty)$, the latter being denoted by a black diamond. One notes
that $c_{1}=-1$ is integer-valued, while the first 50 decimals of the
coefficients $c_{2}$, $c_{3}$ and $c_{4}$ are rational numbers given in Eq.
(29).

##  III Comparison to the Other Methods

Some remarks may be in order with respect to the comparison of the method
outlined above to other established methods like the Aitken’s $\Delta^{2}$
process and Wynn’s epsilon algorithm We1989 ; BrRZ1991 ; CaEtAl2007 . These
are two of the most well-known and frequently used algorithms for series
convergence acceleration. Both of them are based on the Shanks transformation
Sh1955 , which produces Padé approximants if the input data are the partial
sums of a power series We1989 . The former produces an efficient recursive
scheme to compute conveniently such transformation. The formula for the
epsilon algorithm is given by:

$\epsilon_{k+1}^{(n)}=\epsilon_{k-1}^{(n+1)}+\frac{1}{\epsilon_{k}^{(n+1)}-\epsilon_{k}^{(n)}},\qquad\mbox{$-1<k<n-1$}\,,$
(19)

with $k$ indicating the iteration level and $n$ the position within that
level. The initial conditions are $\epsilon_{-1}^{(n)}=0$ for all $n$, and
$\epsilon_{0}^{(n)}=s_{n}$, where $s_{n}$ is the $n$th term of the sequence
being accelerated. The result of the transformation of the partial sums
$(s_{0},\dots,s_{n})$ is $\mathcal{T}_{n}=\epsilon^{(0)}_{n}$, which describes
a sequence of upper-diagonal (odd $n$) and diagonal (even $n$) Padé
approximants to the input series. Aitken’s $\Delta^{2}$ process consists
instead of a recursive iteration of the first-order Shanks transformation:

$\tilde{s}_{n}=s_{n}-\frac{(\Delta s_{n})^{2}}{\Delta^{2}s_{n}},$ (20)

where $\tilde{s}_{n}$ is the improved estimate after applying the process to
the original sequence $\tilde{s}_{n}$, $\Delta s_{n}=s_{n+1}-s_{n}$ is the
first difference of the sequence, and $\Delta^{2}s_{n}=\Delta s_{n+1}-\Delta
s_{n}$ is the second difference of the sequence. One iteration of the process
transforms the sequence of partial sums $(s_{0},\dots,s_{n})$ into the
sequence $(\tilde{s}_{0},\dots,\tilde{s}_{n-1})$. After $n$ iterations, one
obtains a sequence with a single element $\bar{s}_{0}$, where the bar over the
symbol is equivalent to an $n$-fold tilde. The result of the iterated process
is $\mathcal{T}_{n}=\bar{s}_{0}$.

Our assumptions, formulated in Eq. (5), imply that the input series exhibits
logarithmic convergence. Specifically, one has the asymptotic relation

$\lim_{n\to\infty}\frac{s_{n+1}-s}{s_{n}-s}=\lim_{n\to\infty}\frac{r_{n+1}}{r_{n}}=1.$
(21)

This observation indicates that the Shanks transform, which is predicated on
the assumption that, in the sequence’s tail (for large $n$), the ratio

$\frac{s_{n+1}-s}{s_{n}-s}\approx\frac{s_{n+2}-s}{s_{n+1}-s}\approx\lambda,$
(22)

remains constant, may not yield accurate results. In fact, advantage can be
taken in Eq. (22) for a linearly convergent series (i.e., for the case
$\lambda<1$). In this case, the sequence terms are expected to follow

$s_{n}=\sigma_{0}+\sigma_{1}\alpha_{1}^{n}+O(\alpha_{2}^{n}),\quad|\alpha_{2}|<|\alpha_{1}|,$
(23)

indicating that the error in approximating Eq. (22) diminishes geometrically
with $n$. Specifically,

$\frac{s_{n+1}-s}{s_{n}-s}-\frac{s_{n+2}-s}{s_{n+1}-s}=\mathcal{O}\left(\frac{\alpha_{2}}{\alpha_{1}}\right)^{n}.$
(24)

However, for a series whose terms are represented by Eq. (7), the error
reduction is merely polynomial in $n$:

$\frac{s_{n+1}-s}{s_{n}-s}-\frac{s_{n+2}-s}{s_{n+1}-s}=\mathcal{O}\left(\frac{1}{n}\right),$
(25)

which leads to a much slower rate of convergence for large $n$ as compared to
the factor $\lambda=\alpha_{2}/\alpha_{1}<1$. Consequently, methodologies
dependent on such as Aitken’s $\Delta^{2}$ process Ai1926 , the Shanks
transformation Sh1955 , or Wynn’s $\epsilon$ algorithm Wy1956eps , demonstrate
suboptimal performance for these series. For our particular model example,
given in Eq. (26), these statements are illustrated in Fig. 1.

##  IV Numerical Examples

###  IV.1 Slowly Convergent Model Series

We discuss the application of the enhanced Neville algorithm to the infinite
series with partial sums

$s_{n}=\sum_{k=0}^{n}a_{k}\,,\qquad
a_{k}=\frac{4}{\pi\,(k+1)^{2}}\,\arctan\left(\frac{k+2}{k+3}\right)\,,$ (26)

where the prefactor is chosen so that $a_{k}\approx 1/k^{2}$ for $k\to\infty$.
We calculated the enhanced Neville transformations, up to order $n=100$, and
the determination of the convergence according to the formula

$\chi(n)=\frac{\ln(|\mathcal{T}_{n}-\mathcal{T}_{n+1}|)}{\ln(10)}\,,$ (27)

with $\mathcal{T}_{n}$ defined in Eq. (13). The quantity $\chi(n)$ measures
the apparent convergence of the transforms as the order of the transformation
is increased, by calculating the number of apparently converged decimal
digits.

In Fig. 1, we compare the convergence of the enhanced Neville transformation
to Aitken’s $\Delta^{2}$ process Ai1926 , the Shanks transformation Sh1955 ,
Wynn’s epsilon algorithm Wy1956eps  and alongside the series itself. The
enhanced Neville scheme is shown to outperform the other convergence
acceleration algorithms, securing approximately one additional decimal of
convergence with each transformation order. In contrast, the other methods
progress much more slowly. This observation is tied to the asymptotic
structure of the input series, according to Eq. (2).

Our analysis, based on the enhanced Neville algorithm, yields a 50-decimal
approximation of the series limit:

$s_{\infty}=\lim_{n\to\infty}c_{0}(n)=1.31279\,49538\,25865\,79196\,34865\,83906\,40442\,73891\,27574\,77554.$
(28)

Furthermore, we obtain the results

$\displaystyle c_{1}=\lim_{n\to\infty}c_{1}(n)=$ $\displaystyle\;-1\,,\qquad
c_{2}=\lim_{n\to\infty}c_{2}(n)=\frac{1}{2}+\frac{1}{\pi}=0.81830\,98861\,83790\dots\,,$
(29a) $\displaystyle c_{3}=\lim_{n\to\infty}c_{3}(n)=$
$\displaystyle\;-\frac{1}{6}-\frac{2}{\pi}=-0.80328\,64390\,34248\dots\,,\qquad
c_{4}=\lim_{n\to\infty}c_{4}(n)=\frac{37}{12\pi}=0.98145\,54824\,00021\dots\,.$
(29b)

The next higher coefficients are $c_{5}=\frac{1}{30}-\frac{131}{30\pi}$,
$c_{6}=\frac{268}{45\pi}$, $c_{7}=-\frac{1}{42}-\frac{549}{70\pi}$,
$c_{8}=\frac{98347}{10000\pi}$, $c_{9}=-\frac{1}{30}-\frac{9253}{840\pi}$, and
$c_{10}=\frac{66011}{6300\pi}$. These results have been verified to about 100
digits against numerical data obtained with the help of the enhanced Neville
algorithm described here. They can be obtained analytically based on the
formula

$s_{n}=\sum_{k=0}^{n}a_{k}=s_{\infty}-\sum_{k=n+1}^{\infty}a_{k}\,,$ (30)

while realizing that the terms $a_{k}$ in the latter sum from $k=n+1$ to
$k=\infty$ can be expanded for large summation index $k$. We define the
approximations $d_{j}(n)\approx c_{j}$ as follows,

$d_{j}(n)=(n+1)^{j}\left(s_{n}-\sum_{r=0}^{j-1}\frac{c_{r}}{(n+1)^{r}}\right),\qquad\lim_{n\to\infty}d_{j}(n)=c_{j}\,,\qquad
j=1,2,3,4\,.$ (31)

It is instructive to plot the sequence of the $d_{j}(n)$ against $1/n$, for
$j=1,2,3,4$, as shown in Fig. 2, for the model series given in Eq. (26). The
value at infinite $n$ (where $1/n\to 0$) is given by
$\lim_{n\to\infty}d_{j}(n)=c_{j}$. The trend of the data presented in Fig. 2
is consistent with the convergence of the approximants $d_{j}(n)\to c_{j}$ for
$j=1,2,3,4$, for large $n$.

###  IV.2 Bethe Logarithm

We investigate the series of partial sums $s_{n}$, where

$b_{n}=\sum_{k=0}^{n}a_{k}\,,\qquad
a_{k}=\frac{16(k+2)}{(k+1)^{2}\,(k+3)^{2}}\,\Phi\left(-\frac{3+k}{1+k},1,2k+4\right)\,,$
(32)

where $\Phi$ is the Lerch phi transcendent, with
$\Phi(z,s,a)=\sum_{n=0}^{\infty}z^{n}/(n+a)^{s}$. According to Chap. 4 of Ref.
JeAd2022book , the limit $s_{\infty}$ of the input series is related to the
Bethe logarithm Be1947  $\ln k_{0}(1S)$ of the hydrogen ground-state by the
following relation,

$b_{\infty}=\ln k_{0}(1S)-10\,\ln(2)+2\,\zeta(2)+1\,.$ (33)

The sum $b_{\infty}$ and hence, the Bethe logarithm $\ln k_{0}(1S)$ of the
ground state of hydrogen, is calculated in the same manner as the example
treated in Sec. IV.1. The enhance Neville algorithm obtains the 50-figure
result

$\ln
k_{0}(1S)=2.98412\,85557\,65497\,61075\,97770\,90013\,79796\,99751\,80566\,17002$
(34)

in the 58th order of transformation. It is somewhat surprising that the
coefficients $c_{j}$ with $j=1,2,3,4$ are rational numbers,

$\begin{split}c_{1}=&\,0,\qquad c_{2}=\,0,\qquad c_{3}=-\frac{4}{3},\qquad
c_{4}=\frac{27}{4},\qquad c_{5}=-\frac{703}{30},\qquad
c_{6}=\frac{3329}{48},\\\ c_{7}=&-\frac{63163}{336},\qquad
c_{8}=\frac{184961}{384},\qquad c_{9}=-\frac{569323}{480},\qquad
c_{10}=\frac{7256477}{2560}\,.\end{split}$ (35)

The results given in Eq. (35) currently constitute conjectures which we have
verified to at least 100 decimals. The fact that the coefficients
$c_{1},c_{2}=0$ is expected from the observation that
$\lim_{k\to\infty}a_{k}\sim\frac{1}{k^{4}}$, which then leads to
$\lim_{k\to\infty}b_{k}=b_{\infty}+\mathcal{O}\left(\frac{1}{k^{3}}\right)$.
Values of the Bethe logarithms $\ln k_{0}(1S)$, $\ln k_{0}(2S)$ and $\ln
k_{0}(2P)$ accurate to 100 digits are presented in Appendix B. These could be
used for the search of possible analytic formulas using the PSLQ algorithm
FeBa1992 ; BaPl1997 ; FeBaAr1999 ; BaBr2001 . The PSLQ algorithm can be used
to search for analytic expressions of accurately known numerical quantities in
terms of known constants such as the Euler constant $\gamma_{E}=0.57721\dots$,
various Riemann zeta functions, powers of $\pi$, and multiplicative
combinations of these constants. We can report that we have carried out a
limited set of searches with the same constants that were used in Eq. (A11) of
Ref. GiEtAl2020  without success.

##  V Conclusions

For slowly convergent series which possess an asymptotic expansion compatible
with Eq. (2), the enhanced Neville algorithm described here harvests the
asymptotic structure to eliminate the maximum number of asymptotic terms in
the expansion in inverse powers of $n$, in any given order to the
transformation. We take advantage of the asymptotic structure of the terms to
be summed, according to Eq. (2). The general structure of the transformation
has been described in Sec. II.1 and and compared to the recursive Neville
algorithm in Sec. II.2. Our matrix-based formulation of the problem in Sec.
II.1 enables us to not only derive formulas for the limit of the series, but
also for subleading terms, according to Eq. (16).

In Sec. III, we have compared the performance of the enhanced Neville
algorithm to that of Wynn’s $\epsilon$ algorithm Wy1956eps , which calculates
Padé approximants, and Aitken’s $\Delta^{2}$ process Ai1926 , including the
iterated Aitken process. We find that our method provides numerically superior
results, and harvests the asymptotic structure of the input series, according
to Eq. (2). Furthermore, in comparison to the combined nonlinear-condensation
transformation, described in Refs JeMoSoWe1999 ; AkSaJeBeSoMo2003 , the
enhanced Neville transformation eliminates the need to “sample” the input
series at very high orders, which is otherwise necessitated by the van
Wijngaarden transformation NP1961 , which is the first step of the combined
nonlinear-condensation transformation.

For the Bethe logarithm, discussed in Sec. IV.2, we find that the subleading
terms of the partial sums $b_{n}$ of the series given in Eq. (32) can be
expressed in terms of rational coefficients,

$b_{n}=\ln
k_{0}(1S)-10\,\ln(2)+2\,\zeta(2)+1-\frac{4}{3(n+1)^{3}}-\frac{27}{4(n+1)^{4}}+\mathcal{O}\left(\frac{1}{n^{5}}\right)\,,$
(36)

for large $n$. For the definition of $b_{n}$, see Eq. (32). The coefficients
of the subleading terms ($-4/3$ and $-27/4$) have been determined based on our
enhanced Neville algorithm.

Coefficients for the derivation of higher-order subleading asymptotics are
given in Appendix A. Furthermore, Bethe logarithms for the $2S$ and $2P$
states are discussed in Appendix B, and numerical results are given to 100
decimal digits. With the enhanced Neville algorithm presented here, the
accuracy of the calculation of the Bethe logarithms can easily be enhanced
beyond the 100 decimals indicated in Appendix B, facilitating a possible
search for analytic representations using the PSLQ algorithm FeBa1992 ;
BaPl1997 ; FeBaAr1999 ; BaBr2001 .

Finally, as a possible direction for future research, we mention possible
adaptations of the methods discussed here, for input series of a different
asymptotic structure as compared to Eq. (2), for example, input series whose
terms exhibit inverse half-integer powers of the summation index in asymptotic
limits.

##  VI Acknowledgements

The authors acknowledge helpful conversations with Ernst Joachim Weniger and
Jean Zinn–Justin. This research was supported by the National Science
Foundation (grant NSF PHY–2110294). LTG gratefully acknowledges support from
the Swedish Research Council (Vetenskapsradet) Grant No. 638-2013-9243.

##  Appendix A Coefficients for Higher–Order Terms

We report in this Appendix the expressions for the polynomials
$\mathcal{P}_{k}(i,n)$ with $5\leq k\leq 10$. For $k=5$, the result is

$\displaystyle 11520\mathcal{P}_{5}(i,n)=$
$\displaystyle\;11520i+46080i^{2}+69120i^{3}+46080i^{4}+11520i^{5}-11520n-52704in-88320i^{2}n$
$\displaystyle-64320i^{3}n-17280i^{4}n+6768n^{2}+11560in^{2}-2400i^{2}n^{2}-12960i^{3}n^{2}-5760i^{4}n^{2}$
$\displaystyle+7652n^{3}+21000in^{3}+20400i^{2}n^{3}+6720i^{3}n^{3}-680n^{4}+430in^{4}+2640i^{2}n^{4}+1440i^{3}n^{4}$
$\displaystyle-2085n^{5}-2976in^{5}-1200i^{2}n^{5}-395n^{6}-500in^{6}-240i^{2}n^{6}+198n^{7}+120in^{7}$
$\displaystyle+70n^{8}+30in^{8}-5n^{9}-3n^{10}\,.$ (37)

For $k=6$, we have the result

$\displaystyle 2903040\mathcal{P}_{6}(i,n)=$
$\displaystyle\;2903040i+14515200i^{2}+29030400i^{3}+29030400i^{4}+14515200i^{5}+2903040i^{6}$
$\displaystyle-2914560n-16184448in-35538048i^{2}n-38465280i^{3}n-20563200i^{4}n-4354560i^{5}n$
$\displaystyle+1667232n^{2}+4618656in^{2}+2308320i^{2}n^{2}-3870720i^{3}n^{2}-4717440i^{4}n^{2}-1451520i^{5}n^{2}$
$\displaystyle+1942136n^{3}+7220304in^{3}+10432800i^{2}n^{3}+6834240i^{3}n^{3}+1693440i^{4}n^{3}-97020n^{4}$
$\displaystyle-63000in^{4}+773640i^{2}n^{4}+1028160i^{3}n^{4}+362880i^{4}n^{4}-523446n^{5}-1275372in^{5}$
$\displaystyle-1052352i^{2}n^{5}-302400i^{3}n^{5}-146727n^{6}-225540in^{6}-186480i^{2}n^{6}-60480i^{3}n^{6}$
$\displaystyle+44070n^{7}+80136in^{7}+30240i^{2}n^{7}+30177n^{8}+25200in^{8}+7560i^{2}n^{8}+406n^{9}$
$\displaystyle-1260in^{9}-2205n^{10}-756in^{10}-126n^{11}+63n^{12}\,.$ (38)

For $k=7$, the result is

$\displaystyle 5806080\mathcal{P}_{7}(i,n)=$
$\displaystyle\;5806080i+34836480i^{2}+87091200i^{3}+116121600i^{4}+87091200i^{5}+34836480i^{6}$
$\displaystyle+5806080i^{7}-5806080n-38198016in-103444992i^{2}n-148006656i^{3}n-118056960i^{4}n$
$\displaystyle-49835520i^{5}n-8709120i^{6}n+3303936n^{2}+12571776in^{2}+13853952i^{2}n^{2}-3124800i^{3}n^{2}$
$\displaystyle-17176320i^{4}n^{2}-12337920i^{5}n^{2}-2903040i^{6}n^{2}+3783456n^{3}+18324880in^{3}$
$\displaystyle+35306208i^{2}n^{3}+34534080i^{3}n^{3}+17055360i^{4}n^{3}+3386880i^{5}n^{3}-147000n^{4}$
$\displaystyle-320040in^{4}+1421280i^{2}n^{4}+3603600i^{3}n^{4}+2782080i^{4}n^{4}+725760i^{5}n^{4}-920780n^{5}$
$\displaystyle-3597636in^{5}-4655448i^{2}n^{5}-2709504i^{3}n^{5}-604800i^{4}n^{5}-313950n^{6}-744534in^{6}$
$\displaystyle-824040i^{2}n^{6}-493920i^{3}n^{6}-120960i^{4}n^{6}+27573n^{7}+248412in^{7}+220752i^{2}n^{7}$
$\displaystyle+60480i^{3}n^{7}+65187n^{8}+110754in^{8}+65520i^{2}n^{8}+15120i^{3}n^{8}+14457n^{9}-1708in^{9}$
$\displaystyle-2520i^{2}n^{9}-5397n^{10}-5922in^{10}-1512i^{2}n^{10}-1729n^{11}-252in^{11}+273n^{12}$
$\displaystyle+126in^{12}+63n^{13}-9n^{14}.$ (39)

With $C_{8}=1393459200$, we have the following formula for $k=8$,

$\displaystyle C_{8}\,\mathcal{P}_{8}(i,n)=$
$\displaystyle\;1393459200i+9754214400i^{2}+29262643200i^{3}+48771072000i^{4}+48771072000i^{5}+29262643200i^{6}$
$\displaystyle+9754214400i^{7}+1393459200i^{8}-1387653120n-10560983040in-33994321920i^{2}n$
$\displaystyle-60348395520i^{3}n-63855267840i^{4}n-40294195200i^{5}n-14050713600i^{6}n-2090188800i^{7}n$
$\displaystyle+807277824n^{2}+3810170880in^{2}+6342174720i^{2}n^{2}+2574996480i^{3}n^{2}-4872268800i^{4}n^{2}$
$\displaystyle-7083417600i^{5}n^{2}-3657830400i^{6}n^{2}-696729600i^{7}n^{2}+891826560n^{3}+5306000640in^{3}$
$\displaystyle+12871461120i^{2}n^{3}+16761669120i^{3}n^{3}+12381465600i^{4}n^{3}+4906137600i^{5}n^{3}+812851200i^{6}n^{3}$
$\displaystyle-67571600n^{4}-112089600in^{4}+264297600i^{2}n^{4}+1205971200i^{3}n^{4}+1532563200i^{4}n^{4}$
$\displaystyle+841881600i^{5}n^{4}+174182400i^{6}n^{4}-204569280n^{5}-1084419840in^{5}-1980740160i^{2}n^{5}$
$\displaystyle-1767588480i^{3}n^{5}-795432960i^{4}n^{5}-145152000i^{5}n^{5}49594888n^{6}-254036160in^{6}-376457760i^{2}n^{6}$
$\displaystyle-316310400i^{3}n^{6}-147571200i^{4}n^{6}-29030400i^{5}n^{6}-1304520n^{7}+66236400in^{7}+112599360i^{2}n^{7}$
$\displaystyle+67495680i^{3}n^{7}+14515200i^{4}n^{7}+6310455n^{8}+42225840in^{8}+42305760i^{2}n^{8}+19353600i^{3}n^{8}$
$\displaystyle+3628800i^{4}n^{8}+5741280n^{9}+3059760in^{9}-1014720i^{2}n^{9}-604800i^{3}n^{9}+383204n^{10}$
$\displaystyle-2716560in^{10}-1784160i^{2}n^{10}-362880i^{3}n^{10}-825840n^{11}-475440in^{11}-60480i^{2}n^{11}-76790n^{12}$
$\displaystyle+95760in^{12}+30240i^{2}n^{12}+57120n^{13}+15120in^{13}+1260n^{14}-2160in^{14}-1800n^{15}+135n^{16}.$
(40)

With $C_{9}=2786918400$ and $\mathcal{P}_{9}\equiv\mathcal{P}_{9}(i,n)$, the
result reads, for $k=9$,

$\displaystyle C_{9}\,\mathcal{P}_{9}=$
$\displaystyle\;2786918400i+22295347200i^{2}+78033715200i^{3}+156067430400i^{4}+195084288000i^{5}$
$\displaystyle+156067430400i^{6}+78033715200i^{7}+22295347200i^{8}+2786918400i^{9}-2786918400n$
$\displaystyle-23897272320in-89110609920i^{2}n-188685434880i^{3}n-248407326720i^{4}n-208298926080i^{5}n$
$\displaystyle-108689817600i^{6}n-32281804800i^{7}n-4180377600i^{8}n+1642567680n^{2}+9234897408in^{2}$
$\displaystyle+20304691200i^{2}n^{2}+17834342400i^{3}n^{2}-4594544640i^{4}n^{2}-23911372800i^{5}n^{2}-21482496000i^{6}n^{2}$
$\displaystyle-8709120000i^{7}n^{2}-1393459200i^{8}n^{2}+1839755520n^{3}+12395654400in^{3}+36354923520i^{2}n^{3}$
$\displaystyle+59266260480i^{3}n^{3}+58286269440i^{4}n^{3}+34575206400i^{5}n^{3}+11437977600i^{6}n^{3}+1625702400i^{7}n^{3}$
$\displaystyle-193851264n^{4}-359322400in^{4}+304416000i^{2}n^{4}+2940537600i^{3}n^{4}+5477068800i^{4}n^{4}$
$\displaystyle+4748889600i^{5}n^{4}+2032128000i^{6}n^{4}+348364800i^{7}n^{4}-486633040n^{5}-2577978240in^{5}$
$\displaystyle-6130320000i^{2}n^{5}-7496657280i^{3}n^{5}-5126042880i^{4}n^{5}-1881169920i^{5}n^{5}-290304000i^{6}n^{5}$
$\displaystyle-53606320n^{6}-607262096in^{6}-1260987840i^{2}n^{6}-1385536320i^{3}n^{6}-927763200i^{4}n^{6}$
$\displaystyle-353203200i^{5}n^{6}-58060800i^{6}n^{6}+40638488n^{7}+129863760in^{7}+357671520i^{2}n^{7}$
$\displaystyle+360190080i^{3}n^{7}+164021760i^{4}n^{7}+29030400i^{5}n^{7}-6368192n^{8}+97072590in^{8}$
$\displaystyle+169063200i^{2}n^{8}+123318720i^{3}n^{8}+45964800i^{4}n^{8}+7257600i^{5}n^{8}-252565n^{9}+17602080in^{9}$
$\displaystyle+4090080i^{2}n^{9}-3239040i^{3}n^{9}-1209600i^{4}n^{9}+5530785n^{10}-4666712in^{10}-9001440i^{2}n^{10}$
$\displaystyle-4294080i^{3}n^{10}-725760i^{4}n^{10}-77996n^{11}-2602560in^{11}-1071840i^{2}n^{11}-120960i^{3}n^{11}$
$\displaystyle-873524n^{12}+37940in^{12}+252000i^{2}n^{12}+60480i^{3}n^{12}+34690n^{13}+144480in^{13}+30240i^{2}n^{13}$
$\displaystyle+61670n^{14}-1800in^{14}-4320i^{2}n^{14}-6212n^{15}-3600in^{15}-1620n^{16}+270in^{16}+315n^{17}-15n^{18}.$

Finally, with $C_{10}=367873228800$ and
$\mathcal{P}_{10}\equiv\mathcal{P}_{10}(i,n)$, we have the following result
for $k=10$,

$\displaystyle C_{10}\,\mathcal{P}_{10}=$
$\displaystyle\;367873228800i+3310859059200i^{2}+13243436236800i^{3}+30901351219200i^{4}+46352026828800i^{5}$
$\displaystyle+46352026828800i^{6}+30901351219200i^{7}+13243436236800i^{8}+3310859059200i^{9}+367873228800i^{10}$
$\displaystyle-370660147200n-3522313175040in-14917040455680i^{2}n-36669077913600i^{3}n$
$\displaystyle-57696244531200i^{4}n-60285225369600i^{5}n-41842514165760i^{6}n-18608254156800i^{7}n$
$\displaystyle-4813008076800i^{8}n-551809843200i^{9}n+211314216960n^{2}+1435825391616in^{2}+3899225696256i^{2}n^{2}$
$\displaystyle+5034352435200i^{3}n^{2}+1747653304320i^{4}n^{2}-3762781102080i^{5}n^{2}-5991990681600i^{6}n^{2}$
$\displaystyle-3985293312000i^{7}n^{2}-1333540454400i^{8}n^{2}-183936614400i^{9}n^{2}+252909144576n^{3}$
$\displaystyle+1879074109440in^{3}+6435076285440i^{2}n^{3}+12621996288000i^{3}n^{3}+15516933949440i^{4}n^{3}$
$\displaystyle+12257714810880i^{5}n^{3}+6073740288000i^{6}n^{3}+1724405760000i^{7}n^{3}+214592716800i^{8}n^{3}$
$\displaystyle-12322523136n^{4}-73018923648in^{4}-7247644800i^{2}n^{4}+428333875200i^{3}n^{4}+1111124044800i^{4}n^{4}$
$\displaystyle+1349826508800i^{5}n^{4}+895094323200i^{6}n^{4}+314225049600i^{7}n^{4}+45984153600i^{8}n^{4}-77131325216n^{5}$
$\displaystyle-404528688960in^{5}-1149495367680i^{2}n^{5}-1798761000960i^{3}n^{5}-1666196421120i^{4}n^{5}-924952089600i^{5}n^{5}$
$\displaystyle-286634557440i^{6}n^{5}-38320128000i^{7}n^{5}-18618732528n^{6}-87234630912in^{6}-246608991552i^{2}n^{6}$
$\displaystyle-349341189120i^{3}n^{6}-305355536640i^{4}n^{6}-169087564800i^{5}n^{6}-54286848000i^{6}n^{6}-7664025600i^{7}n^{6}$
$\displaystyle+13184713168n^{7}+22506296736in^{7}+64354656960i^{2}n^{7}+94757731200i^{3}n^{7}+69195962880i^{4}n^{7}$
$\displaystyle+25482885120i^{5}n^{7}+3832012800i^{6}n^{7}+3778269000n^{8}+11972980536in^{8}+35129924280i^{2}n^{8}$
$\displaystyle+38594413440i^{3}n^{8}+22345424640i^{4}n^{8}+7025356800i^{5}n^{8}+958003200i^{6}n^{8}-2732570786n^{9}$
$\displaystyle+2290135980in^{9}+2863365120i^{2}n^{9}+112337280i^{3}n^{9}-587220480i^{4}n^{9}-159667200i^{5}n^{9}$
$\displaystyle-185879199n^{10}+114057636in^{10}-1804196064i^{2}n^{10}-1755008640i^{3}n^{10}-662618880i^{4}n^{10}$
$\displaystyle-95800320i^{5}n^{10}+555823886n^{11}-353833392in^{11}-485020800i^{2}n^{11}-157449600i^{3}n^{11}$
$\displaystyle-15966720i^{4}n^{11}-37142193n^{12}-110297088in^{12}+38272080i^{2}n^{12}+41247360i^{3}n^{12}+7983360i^{4}n^{12}$
$\displaystyle-65980420n^{13}+23650440in^{13}+23063040i^{2}n^{13}+3991680i^{3}n^{13}+9237162n^{14}+7902840in^{14}$
$\displaystyle-807840i^{2}n^{14}-570240i^{3}n^{14}+3778940n^{15}-1295184in^{15}-475200i^{2}n^{15}-860970n^{16}-178200in^{16}$
$\displaystyle+35640i^{2}n^{16}-48378n^{17}+41580in^{17}+29205n^{18}-1980in^{18}-2970n^{19}+99n^{20}.$

##  Appendix B Bethe Logarithms for Excited States

We supplement the results presented in Sec. III with a somewhat more detailed
discussion. The Lerch $\Phi$ transcendent AbSt1972  is defined as

$\Phi(z,s,a)=\sum_{n=0}^{\infty}\frac{z^{n}}{(n+a)^{s}}\,.$ (41)

A formula valid in a larger domain of the complex plane is

$\Phi(z,s,a)=\frac{1}{1-z}\;\sum_{n=0}^{\infty}\left(\frac{-z}{1-z}\right)^{n}\sum_{k=0}^{n}\frac{(-1)^{k}}{(a+k)^{s}}\,{n\choose
k}\,.$ (42)

The Lerch transcendent can be thought of as a generalization of the $\zeta$
function, because $\Phi(1,s,1)=\zeta(s)$. The Lerch transcendent, for the case
$s=1$, can be written as a hypergeometric function,
$\Phi(z,1,a)={}_{2}F_{1}(1,a,1+a,z)/a$. According to Chap. 4 of Ref.
JeAd2022book , hydrogenic Bethe logarithms of $1S$, $2S$ and $2P$ states can
be written as rather compact sums over Lerch $\Phi$ functions, and evaluated
to essentially arbitrary accuracy using the convergence acceleration
algorithms outlined here. We have for the ground state,

$\displaystyle\ln k_{0}(1S)=$
$\displaystyle\;10\,\ln(2)-2\zeta(2)-1+\sum_{k=2}^{\infty}\frac{16\,k}{(k-1)^{2}\,(k+1)^{2}}\,\Phi\left(\frac{1+k}{1-k},1,2\,k\right)$
$\displaystyle=$
$\displaystyle\;2.98412\,85557\,65497\,61075\,97770\,90013\,79796\,99751\,80566\,17002\,$
$\displaystyle\,\;\;\;\,00048\,15926\,13924\,06576\,62306\,75532\,86860\,62013\,30404\,72249\,.$
(43)

This result adds 50 figures to the result communicated in Eq. (34). The
logarithmic sum, for the $2S$ state, is given by

$\displaystyle\ln k_{0}(2S)=$
$\displaystyle\;-\frac{545}{36}+\frac{16}{3}\,\ln(2)-14\,\zeta(2)+24\,\zeta(3)+\sum_{k=3}^{\infty}\frac{1024\,k\,(k-1)\,(k+1)}{(k-2)^{3}\,(k+1)^{3}}\,\Phi\left(\frac{2+k}{2-k},1,2\,k\right)$
$\displaystyle=$
$\displaystyle\;2.81176\,98931\,20563\,51521\,97427\,85941\,63611\,28935\,51470\,29732$
$\displaystyle\,\;\;\;\,41909\,18696\,96453\,24020\,20118\,89106\,87017\,48612\,02831\,24031\,,$
(44)

whereas for the $2P$ state, it reads

$\displaystyle\ln k_{0}(2P)=$
$\displaystyle\;-\frac{3437}{2916}+\frac{3280}{2187}\,\ln(2)-\frac{14}{3}\,\zeta(2)+\frac{136}{9}\,\zeta(3)-\frac{64}{3}\,\zeta(4)$
$\displaystyle\;+\sum_{k=3}^{\infty}\frac{256\,k^{3}\,(11\,k^{2}-12)}{3\,(k-2)^{4}\,(k+2)^{4}}\,\Phi\left(\frac{2+k}{2-k},1,2\,k\right)$
$\displaystyle=$
$\displaystyle\;-0.03001\,67086\,30212\,90244\,36757\,10951\,14406\,39409\,33044\,23103$
$\displaystyle\,\;\;\;\;\;\;\;\,04668\,98525\,32719\,44796\,89622\,57183\,26244\,10312\,70799\,73828\,.$
(45)

While these formulas do not provide direct, analytic results for Bethe
logarithms, they still illustrate that one can think of the Bethe logarithm as
a sum of logarithms, zeta functions, and generalized zeta functions.

## References

  * (1) U. D. Jentschura, P. J. Mohr, and G. Soff, Calculation of the Electron Self-Energy for Low Nuclear Charge, Phys. Rev. Lett. 82, 53–56 (1999).
  * (2) U. D. Jentschura, P. J. Mohr, G. Soff, and E. J. Weniger, Convergence acceleration via combined nonlinear-condensation transformations, Comput. Phys. Commun. 116, 28–54 (1999).
  * (3) U. D. Jentschura and G. S. Adkins, Quantum Electrodynamics: Atoms, Lasers and Gravity (World Scientific, Singapore, 2022).
  * (4) L. T. Giorgini, U. D. Jentschura, E. M. Malatesta, G. Parisi, T. Rizzo, and J. Zinn-Justin, Two–Loop Corrections to the Large–Order Behavior of Correlation Functions in the One–Dimensional $N$–Vector Model, Phys. Rev. D 101, 125001 (2020).
  * (5) L. T. Giorgini, U. D. Jentschura, E. M. Malatesta, G. Parisi, T. Rizzo, and J. Zinn-Justin, Correlation Functions of the Anharmonic Oscillator: Numerical Verification of Two–Loop Corrections to the Large–Order Behavior, Phys. Rev. D 105, 105012 (2022).
  * (6) L. T. Giorgini et al., Instantons in $\phi^{4}$ Theories: Transseries, Virial Theorems and Numerical Aspects, in preparation (2024).
  * (7) E. J. Weniger, Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series, Comput. Phys. Rep. 10, 189–371 (1989).
  * (8) C. Brezinski and M. Redivo-Zaglia, Extrapolation Methods (North-Holland, Amsterdam, 1991).
  * (9) E. Caliceti, M. Meyer-Hermann, P. Ribeca, A. Surzhykov, and U. D. Jentschura, From Useful Algorithms for Slowly Convergent Series to Physical Predictions Based on Divergent Perturbative Expansions, Phys. Rep. 446, 1–96 (2007).
  * (10) A. C. Aitken, On Bernoulli’s numerical solution of algebraic equations, Proc. Roy. Soc. Edinburgh 46, 289–305 (1926).
  * (11) D. Shanks, Non-linear transformations of divergent and slowly convergent sequences, J. Math. and Phys. (Cambridge, Mass.) 34, 1–42 (1955).
  * (12) P. Wynn, On a device for computing the $e_{m}(S_{n})$ transformation, Math. Tables Aids Comput. 10, 91–96 (1956).
  * (13) C. W. Clenshaw, E. T. Goodwin, D. W. Martin, G. F. Miller, F. W. J. Olver, and J. H. Wilkinson, (National Physical Laboratory), Modern Computing Methods, Notes on Applied Science, 2 ed. (H. M. Stationary Office, London, 1961), Vol. 16.
  * (14) S. V. Aksenov, M. A. Savageau, U. D. Jentschura, J. Becher, G. Soff, and P. J. Mohr, Application of the combined nonlinear-condensation transformation to problems in statistical analysis and theoretical physics, Comput. Phys. Commun. 150, 1–20 (2003).
  * (15) A. Sidi, Convergence properties of some nonlinear sequence transformations, Math. Comput. 33, 315–326 (1979).
  * (16) A. Sidi, Analysis of the convergence of the $T$-transformation for power series, Math. Comput. 35, 833–850 (1980).
  * (17) A. Sidi, Borel summability and converging factors for some everywhere divergent series, SIAM J. Math. Anal. 17, 1222–1232 (1986).
  * (18) E. J. Weniger, Computation of the Whittaker function of the second kind by summing its divergent asymptotic series with the help of nonlinear sequence transformations, Comput. Phys. 10, 496–503 (1996).
  * (19) A. Sidi, Practical Extrapolation Methods (Cambridge University Press, Cambridge, 2003).
  * (20) C. Brezinski, A general extrapolation algorithm, Numer. Math. 35, 175–180 (1980).
  * (21) T. Havie, Generalized Neville type extrapolation schemes, B.I.T. 19, 204–213 (1979).
  * (22) E. H. Neville, Iterative Interpolation, J. Indian Math. Soc. 20, 87–120 (1934).
  * (23) D. S. Gaunt and A. J. Guttmann, in Phase Transitions and Critical Phenomena 3, edited by C. Domb and M. S. Green (Academic Press, London, 1974), p. 181.
  * (24) A. H. van Tuyl, Acceleration of Convergence of a Family of Logarithmically Convergent Series, Math. Comput. 63, 229–246 (1994).
  * (25) H. A. Bethe, The Electromagnetic Shift of Energy Levels, Phys. Rev. 72, 339–341 (1947).
  * (26) H. R. P. Ferguson and D. H. Bailey, A Polynomial Time, Numerically Stable Integer Relation Algorithm, RNR Techn. Rept. RNR–91–032 (1992).
  * (27) D. H. Bailey and S. Plouffe, in Organic Mathematics: Proceedings of the Workshop Held in Burnaby, BC, edited by J. Borwein, P. Borwein, L. Jörgenson, and R. Corless (American Mathematical Society, Philadelphia, PA, 1997), pp. 73–88.
  * (28) H. R. P. Ferguson, D. H. Bailey, and S. Arno, Analysis of PSLQ, an integer relation finding algorithm, Math. Comput. 68, 351–369 (1999).
  * (29) D. H. Bailey and D. Broadhurst, Integer Relation Detection, Math. Comput. 70, 1719–1736 (2001).
  * (30) M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, 10 ed. (National Bureau of Standards, Washington, D. C., 1972).