Effective Summation and Interpolation of Series by Self-Similar Root Approximants

Abstract

We describe a simple analytical method for effective summation of series, including divergent series. The method is based on self-similar approximation theory resulting in self-similar root approximants. The method is shown to be general and applicable to different problems, as is illustrated by a number of examples. The accuracy of the method is not worse, and in many cases better, than that of Padé approximants, when the latter can be defined.

1. Introduction

In numerous cases of applied mathematics and mathematical physics the solutions to problems can only be represented as series derived by means of some kind of perturbation theory or iterative procedure. A great majority of such series is even divergent, having meaning only as asymptotic series for an infinitesimally small expansion variable. While the considered problems often require to consider finite values of this variable, sometimes even very large values. The standard way of treating such asymptotic series, for the purpose of their extrapolation to the finite values of the variable, is by invoking the Padé approximants [1]. The latter, however, exhibit several deficiencies limiting their applicability, as is discussed in References [1,2], for instance, such a notorious deficiency as the appearance of spurious poles. Another weak point is the ambiguity of choosing one of the Padé approximants P_M/N from the table of many admissible, for each series of order k, variants satisfying the condition M + N + 1 = k. Also, in the limit of a large variable x the approximant P_M/N(x) behaves as x^M−N. Hence, only integer powers of x are allowed. It is possible to improve the results by employing the modified Padé approximants [3], corresponding to the power $P_{M/N}^{\mathrm{\gamma }}$, with choosing the appropriate value of γ satisfying the large-variable limit.

In the present paper, we show that it is possible to formulate a general method for effectively extrapolating and interpolating asymptotic series. The method enjoys the following advantages: (i) It is unambiguously defined for each given series of order k; (ii) It allows for the treatment of large-variable behavior of any type, whether with integer, rational, or irrational powers; (iii) Being more general, it is not less accurate than the method of the Padé approximants, when the latter exist, in many cases, being more accurate.

In the great majority of realistic situations, only a few terms of asymptotic expansions are available. Therefore, in the examples below, we do not consider very large series, showing that even several terms allow us to derive quite accurate approximations.

2. Self-Similar Root Approximants

Suppose we are interested in finding a real function f(x) of a real variable x. This function is defined by a complicated equation that cannot be solved exactly. However, applying a kind of perturbation theory, we can derive the small-variable behavior of this function

$$f(x)\simeq f_k(x)(x\to 0),$$

(1)

represented by asymptotic series, with the k-th order expansion

$$f_k(x)=f_0(x)\left(1+{\displaystyle \sum _{n=1}^k{a}_n{x}^n}\right),$$

(2)

where

$$f_0(x)=A{x}^{\alpha }.$$

(3)

Sometimes, the large-variable behavior of the function

$$f(x)\simeq f^{(p)}(x)(x\to \infty )$$

(4)

is also known and can be represented by an expansion over 1/x as

$$f^{(p)}(x)=f_{\infty }(x)\left(1+{\displaystyle \sum _{n=1}^p\frac{b_n}{x^n}}\right),$$

(5)

with

$$f_{\infty }(x)=B{x}^{\beta }.$$

(6)

For what follows, it is convenient to deal with the ratio f(x)/f₀(x), which at small variable x → 0 behaves as

$$\frac{f(x)}{f_0(x)}\simeq \frac{f_k(x)}{f_0(x)}=1+{\displaystyle \sum _{n=1}^k{a}_n{x}^n},$$

(7)

and at large values of the variable x → ∞ it tends to

$$\frac{f(x)}{f_0(x)}\simeq \frac{f_{\infty }(x)}{f_0(x)}=\frac{B}{A}{x}^{\beta -\alpha }.$$

(8)

The extrapolation of the small-variable expansions to the large-variable region can be done by means of self-similar approximation theory [4–8]. In this approach, the transfer from a k-th order approximation, say, a small-variable expansion, to the higher orders of approximations are treated as the motion with respect to the approximation order k playing the role of discrete time. Constructing a dynamical system, whose trajectory is bijective to the sequence of approximations, makes it feasible to find a fixed point representing the sought function. The convergence to the fixed point is governed by control functions. The self-similar approximation theory combines the methods of optimal control theory, dynamical theory, and renormalization-group approach. We shall not go into the details and mathematical justification of the self-similar approximation theory that has been thoroughly expounded in References [4–8], but we shall use some of its consequences.

Employing this theory for the purpose of interpolation between the small-variable and large-variable regions, it is possible to come [9–11] to the self-similar root approximant

$$\frac{f_k^{*}(x)}{f_0(x)}={\left({\left(\dots {\left(1+A_{1}x\right)}^{n_1}+A_2{x}^2\right)}^{n_2}+\dots +A_k{x}^k\right)}^{n_k}.$$

(9)

A theorem has been proved [12] stating that all parameters A_i and powers n_i of approximant (9) are uniquely defined through the large-variable form (5).

However, the root approximant (9) cannot be uniquely defined through the small-variable expansion (2). This hinders the applicability of the approximant (9), since in the majority of cases, the small-variable expansion is better known, providing a number of terms, while the knowledge of the large-variable behavior is limited by just a single term (6), often even without precise data for the amplitude B. In order to extend the applicability of approximant (9) to be uniquely defined through the small-variable expansion, it is necessary to impose some constraints on the powers n_j. Such a straightforward constraint is the requirement that all parameters A_j of approximant (9) be involved in the definition of the large-variable limit, which implies the relation

$$n_j=\frac{j+1}{j}(j=1,2,\dots ,k-1),$$

(10)

with n_k = β − α. By expanding Eq. (9) in powers of x, it is easy to prove that all parameters A_j are uniquely defined through the coefficients a_j of small-variable expansion (2). In addition, we can require the validity of the limiting form (6), which improves accuracy.

The self-similar root approximant (9), with conditions (6) and (10), whose parameters A_j are uniquely defined by the accuracy-through-order procedure and are expressed through the coefficients a_j of the small-variable expansion (2), can be called, for short, the root approximant. In the following sections, we demonstrate that this root approximant provides quite accurate approximations for different problems, uniformly extrapolating the small-variable expansion (2), valid for x → 0, to the whole region of x ∈ [0, ∞].

3. Illustration by Simple Examples

Before going to more complicated problems, we show the efficiency of the method by simple cases.

3.1. Hard-Core Scattering Problem

Let us start the illustration of the method from the problem considered by Baker and Gammel [3]. When calculating the scattering length of a repulsive square-well potential, one meets the integral

$$S(x)={\displaystyle {\int }_0^x{\left(\frac{\sin t}{t^3}-\frac{\cos t}{t^2}\right)}^{2}dt,}$$

whose limit, as x → ∞, equals π/15. Baker and Gammel state that this integral cannot be correctly evaluated by the standard Padé method. To solve the problem, they suggest a modified method employing a power of the Padé approximant. We show below that such integrals can easily be treated by means of the root approximants.

The small-variable expansion of this integral reads as

$$S(x)\simeq \frac{x}{9}-\frac{x^3}{135}+\frac{x^5}{2625}-\frac{4x^7}{297675}+\frac{2x^9}{5893965}-\frac{x^{11}}{166080925}+\frac{x^{13}}{10672286625}.$$

Comparing this with form (2), we have S₀(x) = x/9. Since expansion (2) is in powers of x², we construct the root approximants (9) using x² as a variable. Thus, the root approximant of third order is

$$S_3^{*}(x)=\frac{x}{9}{\left({\left({\left(1+A_1{x}^2\right)}^2+A_2{x}^4\right)}^{3/2}+A_3{x}^6\right)}^{-1/6},$$

where the parameters are

$$A_1=0.133333,A_2=0.012952,A_3=\mathrm{0.016907.}$$

To fourth order,

$$S_4^{*}(x)=\frac{x}{9}{\left({\left({\left({\left(1+A_1{x}^2\right)}^2+A_2{x}^4\right)}^{3/2}+A_3{x}^6\right)}^{4/3}+A_4{x}^8\right)}^{-1/8},$$

where

$$A_1=0.133333,A_2=0.012952,A_3=0.002757,A_4=\mathrm{0.004636.}$$

To fifth order,

$$S_5^{*}(x)=\frac{x}{9}{\left({\left({\left({\left({\left(1+A_1{x}^2\right)}^2+A_2{x}^4\right)}^{3/2}+A_3{x}^6\right)}^{4/3}+A_4{x}^8\right)}^{5/4}+A_5{x}^{10}\right)}^{-1/10},$$

where

$$A_1=0.133333,A_2=0.012952,A_3=0.002757,A_4=0.000578,A_5=\mathrm{0.001285.}$$

And to sixth order,

$$S_6^{*}(x)=\frac{x}{9}{\left({\left({\left({\left({\left({\left(1+A_1{x}^2\right)}^2+A_2{x}^4\right)}^{3/2}+A_3{x}^6\right)}^{4/3}+A_4{x}^8\right)}^{5/4}+A_5{x}^{10}\right)}^{6/5}+A_6{x}^{12}\right)}^{-1/12},$$

where

$$\begin{array}{l}{A}_1=0.133333,A_2=0.012952,A_3=0.002757,\\ A_4=0.000578,A_5=0.000137,A_6=\mathrm{0.000356.}\end{array}$$

All these approximants converge to π/15, as x → ∞. The higher the approximant order, the faster the convergence.

3.2. Debye Function

The n-th order Debye function is defined [13] through the integral representation

$$D(n,x)\equiv \frac{n}{x^n}{\displaystyle {\int }_0^x\frac{t^n}{e^t-1}dt.}$$

For |x| < 2π and n ≥ 1, it possesses the expansion

$$D(n,x)\simeq 1-\frac{n}{2(n+1)}x+n{\displaystyle \sum _{k=1}^{\infty }\frac{B_{2k}}{(2k+n)(2k)!}{x}^{2k},}$$

in which B_2k are Bernoulli numbers. At large x and Re n > 0, one has

$$D(n,x)\simeq \frac{C_n}{x^n}(x\to \infty ,\mathrm{Re}n>0),$$

where

$$C_n\equiv n\mathrm{\Gamma }(n+1)\zeta (n+1).$$

Below, we consider the case of n = 3, corresponding to the Debye function

$$D(x)\equiv D(3,x)=\frac{3}{x^3}{\displaystyle {\int }_0^x\frac{t^3}{e^t-1}dt.}$$

The small-variable expansion for the latter takes the form

$$D(x)\simeq 1-\frac{3}{8}x+{\displaystyle \sum _{k=1}^{\infty }{a}_{2k}{x}^{2k}}(x\to 0),$$

in which

$$a_{2k}=\frac{B_{2k}}{(2k+3)(2k)!}.$$

While the large-variable behavior is given by the expression

$$D(x)\simeq \frac{C_3}{x^3}(x\to \infty ),$$

with

$$C_3=\frac{{\pi }^4}{5}=19.481818.$$

Constructing the root approximant

$$D_5^{*}(x)=\frac{x}{9}{\left({\left({\left({\left({\left(1+A_{1}x\right)}^2+A_2{x}^2\right)}^{3/2}+A_3{x}^3\right)}^{4/3}+A_4{x}^4\right)}^{5/4}+A_5{x}^5\right)}^{-3/5},$$

we compare it with the exact numerical values of the function D(x) and find that $D_5^{*}(x)$ approximates well this function in the whole region of x ∈ [0, ∞], with the maximal error of 15% at x = 5. The best two-point Padé approximant of the same order, P_1/4(x), is less accurate, yielding the maximal error of 33% at x = 15.

3.3. Fermi-Dirac Integral

The general form of the j-th order Fermi-Dirac integral is

$$F(j,x)=\frac{1}{\mathrm{\Gamma }(j+1)}{\displaystyle {\int }_0^{\infty }\frac{t^j}{e^{t-x}+1}dt.}$$

Its asymptotic expansions are known [14].

For concreteness, let us consider the zero-order case that reduces to the function

$$F(x)\equiv F(0,x)=\ln (1+e^x).$$

At small x, this function tends to ln 2, and at large x, we have

$$F(x)\simeq x(x\to \infty ).$$

The root approximant

$$F_5^{*}(x)=\ln 2{\left({\left({\left({\left({\left(1+A_{1}x\right)}^2+A_2{x}^2\right)}^{3/2}+A_3{x}^3\right)}^{4/3}+A_4{x}^4\right)}^{5/4}+A_5{x}^5\right)}^{1/5},$$

where

$$\begin{array}{l}{A}_1=0.721348,A_2=0.360674,A_3=0.390257,\\ A_4=0.410334,A_5=4.294519,\end{array}$$

provides an accurate approximation for the function F(x) in the whole region of x ∈ [0, ∞], the maximal error being 5%. The two-point Padé approximant P_3/2(x) is slightly less accurate, with the maximal error of 6%.

3.4. Fekete-Szegö Problem

The problem of maximizing the absolute value of a functional in subclasses of normalized functions is called the Fekete-Szegö problem [15,16]. The Fekete-Szegö functional is bounded by the function

$$f(x)=1+2\exp \left(-\frac{2x}{1-x}\right),$$

where 0 < x < 1.

In order to consider the interval [0, ∞], as in other examples, we can use the change of the variable

$$x=\frac{z}{1+z},z=\frac{x}{1-x}.$$

Then z → ∞ as x → 1. Expanding F(z) at small z gives

$$F(z)\equiv f(x(z))\simeq 3-4z+4z^2-\frac{8}{3}{z}^3+\frac{4}{3}{z}^4-\frac{8}{15}{z}^5.$$

The root approximant $F_3^{*}(z)$ uniformly approximates the function F(z) on the interval z ∈ [0, ∞], with the maximal error about 10%. The two-point Padé approximant P_2/2(z) is worse, having the maximal error twice larger than the root approximant $F_3^{*}(z)$.

4. Some Useful Tricks

It is important to mention some tricks allowing for the convenient use of the method. Below, we discuss the interchange of small-variable and large -variable limits and the problem of dealing with logarithms.

4.1. Inversion of Expansions

In the above examples, we have considered functions, whose expansions are better known for the small-variable limit, while a few, or just a single term, are available in the large-variable limit. But generally, the small-variable and large-variable limits are interchangeable. In those cases, when the large-variable expansion in powers of 1/x provides a number of terms and this expansion enjoys better convergence properties, it is possible to inverse the small-variable limit to the large-variable limit by using the variable change x = 1/t. Then, instead of the function f(x), we consider the function

$$F(t)\equiv f\left(\frac{1}{t}\right),t=\frac{1}{x}.$$

(11)

The small-variable limit (1) becomes the large-variable limit

$$F(t)\simeq F^{(k)}(t)\equiv f_k\left(\frac{1}{t}\right)(t\to \infty ),$$

(12)

in which

$$F^{(k)}(t)=F_{\infty }(t)\left(1+{\displaystyle \sum _{n=1}^k\frac{a_n}{t^n}}\right),$$

(13)

with

$$F_{\infty }(t)\equiv f_0\left(\frac{1}{t}\right)=A{t}^{-\alpha }.$$

Conversely, the large-variable behavior (4) transforms to the small-variable behavior

$$F(t)\simeq F_p(t)\equiv f^{(p)}\left(\frac{1}{t}\right)(t\to 0),$$

(14)

in which

$$F_p(t)=F_0(t)\left(1+{\displaystyle \sum _{n=1}^p{b}_n{t}^n}\right),$$

(15)

where

$$F_0(t)\equiv f_{\infty }\left(\frac{1}{t}\right)=B{t}^{-\beta }.$$

After this change of the variable it is straightforward to employ the same procedure of constructing the root approximants, as is explained in Section 2.

More generally, it is possible to use the change of the variable t = 1/x^s, with a positive power s > 0, so that again t → ∞, when x → 0.

4.2. Example of Inversion

As an illustration of the inversion procedure, we give below a typical example, discussing it rather briefly, since the whole method of constructing the root approximants is the same as before.

Let us consider the partition function of the so-called zero-dimensional oscillator, or the generating functional of zero-dimensional φ⁴ field theory, which is defined through the integral

$$I(x)=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\exp \left(-{\phi }^2-x{\phi }^4\right)d\phi ,}$$

where x plays the role of a coupling parameter. In the weak-coupling limit, one has [17] the asymptotic expansion

$$I(x)\simeq 1+{\displaystyle \sum _{n=1}^{\infty }{a}_n{x}^n}(x\to 0),$$

in which the coefficients are

$$a_n=\frac{{\left(-1\right)}^n}{\sqrt{\pi }n!}\mathrm{\Gamma }\left(2n+\frac{1}{2}\right).$$

For instance

$$a_1=-\frac{3}{4},a_2=\frac{105}{32},a_3=-\frac{3465}{128},$$

and so on.

The strong-coupling expansion reads as

$$I(x)\simeq 1.022765x^{-1/4}-0.345684x^{-3/4}+0.127846x^{-5/4}(x\to \infty ).$$

here the strong-coupling expansion provides a number of terms. Moreover, the absolute values of the coefficients in this expansion diminish with increasing order, contrary to the coefficients a_n in the weak-coupling expansion, which grow as nⁿ with increasing order n. This makes the strong-coupling expansion more suitable for constructing root approximants.

Resorting to the change of the variable x = 1/t⁴, we consider the function J(t) ≡ I(1/t⁴) and follow the scheme of the previous section. We define the root approximants $J_k^{*}(t)$ that give us the approximants $I_k^{*}(x)=J_k^{*}\left(1/x^{1/4}\right)$ for the sought function. Found in that way approximant $I_3^{*}(x)$ has the maximal error of 5% for the whole range of x ∈ [0, ∞]. For comparison, the Padé approximant P_1/2(x) has the maximal error of about 20%, which is much less accurate.

4.3. Dealing with Logarithms

It is worth paying attention to the problem of series involving logarithms, which often appear in physics applications. Such series do not yield any complication for the method of root approximants described here. There are two equivalent ways of treating such series. Thus, if a series contains the terms with xⁿ, xⁿ⁺¹, and with xⁿ ln x, then it is admissible to consider as the terms of one order either those containing xⁿ and xⁿ ln x or the terms xⁿ⁺¹ and xⁿ ln x.

As an illustration, let us consider, e.g., the typical form of such a series involving logarithms as that one arising in the Nambu-Iona Lasinio model [18] and leading to the function

$$f(x)=x\left[\sqrt{1+x^2}-x^2\ln \left(\frac{1+\sqrt{1+x^2}}{x}\right)\right],$$

where x plays the role of mass. At asymptotically small x, it follows

$$f(x)\simeq x+\left(\frac{1}{2}-\ln 2+\ln x\right)x^3(x\to 0).$$

While at large x, one has

$$f(x)\simeq \frac{2}{3}-\frac{1}{5x^2}+\frac{3}{28x^4}(x\to \infty ).$$

Keeping in mind the dependence of the last expansion on 1/x², it is convenient to use the variablez = 1/x². The root approximant, satisfying the required limits, has the form

$$f_4^{*}(x)=\frac{2}{3}{\left({\left({\left(1+A_{1}z\right)}^2+A_2{z}^2\right)}^{3/2}+A_3{z}^2\ln (1+z)+A_4{z}^3\right)}^{-1/6},$$

with all parameters uniquely defined by the given expansions. This expression approximates well the initial function f(x), with the maximal error of 2% at x ≈ 2. Contrary to this, the best Padé approximant of the same order has the error of 11% at x ≈ 1.5.

5. Ground-State Energy of Electron Gas

Important and not trivial problems arise when studying the properties of charged systems [19–21]. Here we show how our method works for the case of homogeneous electron systems.

5.1. One-Dimensional Electron Gas

The Hartree-Fock part of the uniform electron energy is well known. The problem arises in calculating the correlation energy. The latter is usually presented in a reduced dimensionless form ε(r_s) as a function of the Seitz radius r_s. High-density expansion for one-dimensional uniform electron gas [22] corresponds to small r_s, when for the correlation energy one has

$$\epsilon (r_s)\simeq C+0.00845r_s(r_s\to 0),$$

where

$$C=-\frac{{\pi }^2}{360}=-\mathrm{0.027416.}$$

The low-density expansion [22] implies large r_s, when

$$\epsilon (r_s)\simeq \frac{b_1}{r_s}+\frac{b_2}{r_s^{3/2}}(r_s\to \infty ),$$

where

$$b_1=-\left(\ln \sqrt{2\pi }-\frac{3}{4}\right)=-0.168939,b_2=\mathrm{0.359933.}$$

The root approximant, enjoying the same expansions, but valid for arbitrary r_s reads as

$${\epsilon }_3^{*}(r_s)=-\frac{{\pi }^2}{360}{\left({\left({\left(1+A_1{r}_s\right)}^{3/2}+A_2{r}_s^2\right)}^{5/4}+A_3{r}_s^3\right)}^{-1/3},$$

with the parameters

$$A_1=0.493150,A_2=0.056122,A_3=\mathrm{0.004274.}$$

Comparing the prediction of the root approximant with the data from diffusion Monte Carlo calculations [22] in the interval 0 < r_s < 20, we find that the maximal error of ${\epsilon }_3^{*}$ is 8%. Padé approximants give the errors between 2% and 10%. Thus, $P_{1/2}\left(\sqrt{r_s}\right)$ has the error of 2%, while $P_{0/3}\left(\sqrt{r_s}\right)$ has the maximal error of 10%. The Cioslowski interpolation method [23] results [22] in a better accuracy of 1%. However, this method includes an additional parameter that is fitted from numerical Monte Carlo calculations. While our aim has been in constructing good approximations without fitting parameters, being based only on asymptotic expansions. The principal importance of avoiding fitting parameters is crucial for those problems where no exact numerical data are available.

5.2. Two-Dimensional Electron Gas

Correlation energy of a homogenous two-dimensional electron gas was studied in several articles, e.g., in References [24–31]. In high-density limit (small r_s), the ground-state energy reads [31] as

$$E_0(r_s)\simeq \frac{c_{-2}}{r_s^2}+\frac{c_{-1}}{r_s}+\epsilon (r_s)(r_s\to 0),$$

where the first two terms constitute the Hartree-Fock energy, with

$$c_{-2}=\frac{1}{2},c_{-1}=-\frac{4\sqrt{2}}{3\pi }.$$

And the last term is the correlation energy

$$\epsilon (r_s)\simeq c_0+{c^{\prime }}_1{r}_s\ln r_s(r_s\to 0),$$

with the coefficients

$$c_0=-0.192495,{c^{\prime }}_1=-\sqrt{2}\left(\frac{10}{3\pi }-1\right)=-\mathrm{0.0863136.}$$

In the low-density limit (large r_s) the asymptotic expansion for the correlation energy can be written [26] as

$$\epsilon (r_s)\simeq \frac{b_1}{r_s}+\frac{b_2}{r_s^{3/2}}+\frac{b_3}{r_s^2}(r_s\to \infty ),$$

where

$$b_1=-0.472189,b_2=0.4964,b_3=-\mathrm{0.5297.}$$

For intermediate r_s, there have been suggested [27,28,30] several phenomenological expressions with parameters fitted from Monte Carlo calculations. Thus, Gori-Giorgi et al. [28] suggested the form

$$\epsilon (r_s)=A_0+\left(B_0{r}_s+C_0{r}_s^2+D_0{r}_s^3\right)\ln \left(1+\frac{1}{E_0{r}_s+F_0{r}_s^{3/2}+G_0{r}_s^2+H_0{r}_s^3}\right),$$

with the parameters

$$\begin{array}{c}{A}_0=-0.1925,B_1=0.0863136,C_0=0.057234,D_0=\mathrm{0.003362896.}\\ E_0=1.0022,F_0=-0.02069,G_0=0.34,H_0=\mathrm{0.01747.}\end{array}$$

This expression can be used as a numerical result for estimating the accuracy of approximate analytic formulas.

The root approximant, satisfying all asymptotic expansions reads as

$${\epsilon }_5^{*}(r_s)=\frac{b_1}{r_s}{\left({\left({\left(1+\frac{A_1}{\sqrt{r_s}}\right)}^2+\frac{A_1}{r_s}\right)}^{3/2}+\frac{A_3}{r_s}\ln \left(1+\frac{1}{\sqrt{r_s}}\right)+\frac{A_4}{r_s^{3/2}}+\frac{A_5}{r_s^2}\right)}^{-1/2},$$

where the parameters are

$$\begin{array}{l}{b}_1=-0.472189,A_1=0.700849,A_2=2.723702,\\ A_3=10.792193,A_4=-5.764339,A_5=\mathrm{6.017150.}\end{array}$$

The error of this approximant is about 5%.

6. Systems with Spherical Symmetry

Finite quantum systems often enjoy spherical symmetry. Below, we consider two examples of such systems that are important for applications.

6.1. Energy of Harmonium Atoms

An N-electron harmonium atom is described by the Hamiltonian

$$\widehat{H}=\frac{1}{2}{\displaystyle \sum _{i=1}^N\left(-{\nabla }_i^2+{\omega }^2{r}_i^2\right)}+\frac{1}{2}{\displaystyle \sum _{i\ne j}^N\frac{1}{r_{ij}}},$$

where dimensionless units are employed and

$$r_i\equiv \left|r_i\right|,r_{ij}\equiv \left|r_i-r_j\right|.$$

This Hamiltonian provides a rather realistic modeling of trapped ions, quantum dots, and some other finite systems, such as atomic nuclei and metallic grains [32]. This is why the energy of harmonium atoms has been intensively studied [33–37]. Here we show that root approximants give a good approximation for the energy of such systems. We consider the ground-state energy of a two-electron harmonium.

At a shallow harmonic potential, the energy can be expanded [23] in powers of ω, so that

$$E(\omega )\simeq E_k(\omega )(\omega \to 0),$$

with the truncated series

$$E_k(\omega )={\displaystyle \sum _{n=0}^k{c}_n}{\omega }^{\left(2+n\right)/3}.$$

For instance, to third order, we get

$$E_3(\omega )=c_0{\omega }^{2/3}+c_1\omega +c_2{\omega }^{4/3},$$

with the coefficients

$$c_0=\frac{3}{2^{4/3}}=1.19055,c_1=\frac{1}{2}\left(3+\sqrt{3}\right)=2.36603,c_2=\frac{7}{36}{2}^{-2/3}=\mathrm{0.122492.}$$

And for a rigid potential, the energy is approximated [23] as

$$E(\omega )\simeq E^{(p)}(\omega )(\omega \to \infty ),$$

where

$$E^{(p)}(\omega )={\displaystyle \sum _{n=0}^P{b}_n{\omega }^{(2-n)/2}}.$$

To fourth order, one has

$$E^{(4)}(\omega )=b_0\omega +b_1{\omega }^{1/2}+b_2+b_3{\omega }^{-1/2},$$

where

$$\begin{array}{c}{b}_0=3,b_1=\sqrt{\frac{2}{\pi }}=0.797885,b_2=-\frac{2}{\pi }\left(1-\frac{\pi }{2}+\ln 2\right)=-0.077891,\\ b_3={\left(\frac{2}{\pi }\right)}^{3/2}\left[2-2G-\frac{3}{2}\pi +\left(\pi +3\right)\ln 2+\frac{3}{2}{\left(\ln 2\right)}^2-\frac{{\pi }^2}{24}\right]=0.0112528,\end{array}$$

with the Catalan constant

$$G\equiv {\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(2n+1\right)}^2}}=\mathrm{0.91596559.}$$

The root approximant, respecting all given small-ω, as well as large-ω expansions, is

$$E_6^{*}(\omega )=c_0{\omega }^{2/3}{\left({\left({\left({\left({\left({\left(1+A_1{\omega }^{1/3}\right)}^{1/2}+A_2{\omega }^{2/3}\right)}^{3/4}+A_3\omega \right)}^{5/6}+A_4{\omega }^{4/3}\right)}^{7/8}+A_5{\omega }^{5/3}\right)}^{9/10}+A_6{\omega }^2\right)}^{1/6},$$

with the parameters

$$\begin{array}{c}{c}_0=1.19055,A_1=48.4532,A_2=564.108,\\ A_3=1088.39,A_4=1221.08,A_5=796.791,A_6=256.\end{array}$$

We estimate the accuracy of the root approximant comparing it with the numerical data from Reference [38] and find that its maximal error is only 0.9%. Note that Padé approximants cannot be used in the case of harmonium, since the small-variable and large-variable asymptotic expansions are incompatible.

6.2. Energy of Two-Electron Spherium

The two-electron spherium is a system consisting of two electrons that are confined to the surface of a sphere of radius R. The ground-state energy of the system [23,39] possesses the small-radius expansion

$$E(R)\simeq \frac{1}{R}+c_0+c_{1}R+c_2{R}^2+c_3{R}^3(R\to 0),$$

in which

$$\begin{array}{c}{c}_0=4\ln 2-3=-0.22741128,c_1=8{\left(\ln 2\right)}^2-40\ln 2+24=0.11773689,\\ c_2=-0.05027560,c_3=\mathrm{0.01395783.}\end{array}$$

The coefficients c₂ and c₃ can also be expressed in closed forms that, however, are too much cumbersome [39], because of which we give here only their numerical values.

In the large-radius limit, the energy has the expansion

$$E(R)\simeq \frac{1}{2R}+\frac{1}{2R^{3/2}}-\frac{1}{8R^2}-\frac{1}{128R^{5/2}}(R\to \infty ).$$

The root approximant can be writen in the form

$$E_5^{*}(R)=\frac{1}{R}+c_0{\left({\left({\left({\left({\left(1+A_{1}R\right)}^{3/2}+A_2{R}^2\right)}^{5/4}+A_3{R}^3\right)}^{7/6}+A_4{R}^4\right)}^{9/8}+A_5{R}^5\right)}^{-1/5},$$

where

$$\begin{array}{c}{A}_1=1.05188915,A_2=0.56453530,A_3=0.36000617,\\ A_4=0.12606787,A_5=\mathrm{0.01946301.}\end{array}$$

Comparing this expression with numerical data [39], we find that the maximal error occurs at R = 20, being only 0:1%. The best Padé approximant $P_{5/5}\left(\sqrt{R}\right)$ is much less accurate, having the maximal error, also at R = 20, but an order larger, 1:5%.

7. Discussion

We have described a simple and general method for interpolating functions between their small-variable and large-variable asymptotic expansions. The method is based on the construction of self-similar root approximants enjoying the general form

$$f_k^{*}(x)=f_0(x){\left({\left({\left(\dots {\left(1+A_{1}x\right)}^{n_1}+A_2{x}^2\right)}^{n_2}+A_3{x}^3\right)}^{n_3}+\dots +A_k{x}^k\right)}^{n_k}.$$

All parameters A_i can be uniquely defined through the corresponding asymptotic expansions. By changing the variable, it is easy to inverse the expansions between the small-variable and large-variable limits.

Our aim has been to suggest a method that would involve no fitting parameters. This is especially important in those complicated cases, where numerical data in the whole region of the variable are not available. The absence of fitting parameters makes our aproach different from other intrepolation methods, such as the Cioslowski method [23].

We have demonstrated the method of root approximants by several examples, whose structure is typical for many applications, including the hard-core scattering problem, Debye function, Fermi-Dirac integral, Fekete-Szegö problem, zero-dimensional oscillator, homogeneous electron gas, harmonium atom, and spherium.

We have analyzed several more problems, e.g., the interpolation of the polaron mass between weak-coupling and strong-coupling limits studied earlier by the Feynman variational procedure [40] and by other methods [41–44]. Our approach provides approximations, whose accuracy is comparable or better than that of other methods, being at the same time more simple.

Generally, the suggested method provides the accuracy not worse than the method of Padé approximants and in the majority of cases is more accurate than the latter.

Except the root approximants of the general form (9), we also have considered additive approximants represented by the sums

$$f_{M/N}^{*}(x)={\displaystyle \sum _{i=1}^{(M+N)/2}{A}_i{\left(1+B_{i}x\right)}^{n_i}}.$$

This type of expressions can be considered either as additive root approximants or an additive variant resulting from self-similar factor approximants [45].

For example, in the case of one-dimensional electron gas, the correlation energy is approximated as

$${\epsilon }_{2/2}^{*}(r_s)=A_1{\left(1+B_1{r}_s\right)}^{-1}+A_2{\left(1+B_2{r}_s\right)}^{-3/2},$$

with the parameters

$$A_1=-0.044941,A_2=0.017526,B_1=0.266023,B_2=\mathrm{0.133344.}$$

This expression has the maximal error of 11%. However, a more detailed analysis of such additive approximants requires a separate investigation, which is out of the scope of the present paper.