Extended Comment on "One-Range Addition Theorems for Coulomb Interaction Potential and Its Derivatives" by I. I. Guseinov (Chem. Phys. Vol. 309 (2005), pp. 209 - 213)

TL;DR

This paper critically examines Guseinov's derivation of one-range addition theorems for the Coulomb potential, highlighting mathematical flaws and emphasizing the need for proper interpretation as generalized functions.

Contribution

The paper provides a rigorous critique of Guseinov's methods, demonstrating the invalidity of his addition theorems for the Coulomb potential and clarifying their proper mathematical context.

Findings

Guseinov's addition theorems are mathematically flawed.

The Coulomb potential does not belong to the Hilbert spaces used.

Proper interpretation requires viewing these as generalized function expansions.

Abstract

Addition theorems are principal tools that express a function $f (\bm{r} \pm \bm{r}')$ in terms of products of other functions that only depend on either $\bm{r}$ or $\bm{r}'$. The best known example of such an addition theorem is the Laplace expansion of the Coulomb potential which possesses a characteristic two-range form. Guseinov [Chem. Phys. {\bf 309}, 209 - 213 (2005)] derived one-range addition theorems for the Coulomb potential via the limit $\beta \to 0$ in previously derived one-range addition theorems for the Yukawa potential $\exp \bigl(-\beta | \bm{r}-\bm{r}'| \bigr) /| \bm{r}-\bm{r}'|$. At first sight, this looks like a remarkable achievement, but from a mathematical point of view, Guseinov's work is at best questionable and in some cases fundamentally flawed. One-range addition theorems are expansions in terms of functions that are complete and orthonormal in a given Hilbert space, but Guseinov replaced the complete and orthonormal functions by nonorthogonal Slater-type functions and rearranged the resulting expansions. This is a dangerous operation whose validity must be checked. It is shown that the one-center limit $\bm{r}' = \bm{0}$ of Guseinov's rearranged Yukawa addition theorems as well as of several other addition theorems does not exist. Moreover, the Coulomb potential does not belong to any of the Hilbert spaces implicitly used by Guseinov. Accordingly, one-range addition theorems for the Coulomb potential diverge in the mean. Instead, these one-range addition theorems have to interpreted as expansions of generalized functions in the sense of Schwartz that converge weakly in suitable functionals.

Equations

$$F_{\ell }^m(\mathbf{r})=f_{\ell }(r)Y_{\ell }^m(\mathbf{r}/r),$$

$$\frac{1}{|\mathbf{r}\pm {\mathbf{r}}^{\prime }|}$$

$$|{\mathbf{r}}_{<}|=\min (r,r^{\prime }),|{\mathbf{r}}_{>}|=\max (r,r^{\prime }).$$

$$f(\mathbf{r}\pm {\mathbf{r}}^{\prime })=\sum _{n=0}^{\infty }\frac{(\pm {\mathbf{r}}^{\prime }\cdot \nabla {)}^n}{n!}f(\mathbf{r})={\mathrm{e}}^{\pm {\mathbf{r}}^{\prime }\cdot \nabla }f(\mathbf{r}).$$

$${\mathrm{e}}^{{\mathbf{r}}^{\prime }\cdot \nabla }={\mathrm{e}}^{x^{\prime }\partial /\partial x}{\mathrm{e}}^{y^{\prime }\partial /\partial y}{\mathrm{e}}^{z^{\prime }\partial /\partial z}$$

$${\chi }_{N,L}^M(\beta ,\mathbf{r})=(\beta r{)}^{N-L-1}{\mathrm{e}}^{-\beta r}{\mathcal{Y}}_L^M(\beta \mathbf{r}).$$

$$\Vert u\Vert =\sqrt{(u|u)},$$

$$f_N=\sum _{n=0}^N{C}_n^{(N)}{\phi }_n$$

$$f=\sum _{n=0}^{\infty }({\phi }_n|f){\phi }_n$$

$$f=\sum _{n=0}^{\infty }{C}_n{\phi }_n$$

$$\Vert f{\Vert }^2=\sum _{n=0}^{\infty }|({\phi }_n|f){|}^2.$$

$$(f|g)_{2}\;=\;\int\bigl{[}f(\bm{r})\bigr{]}^{*}\,g(\bm{r})\,\mathrm{d}^{3}\bm{r}\,.$$

$$\Vert f{\Vert }_2=\sqrt{(f|f{)}_2}.$$

$$L^2({\mathbb{R}}^3)$$

$$\displaystyle\;=\;\bigl{\{}f\colon\mathbb{R}^{3}\to\mathbb{C}\big{|}\,\|f\|_{2}<\infty\bigr{\}}$$

$$(f|g)_{w,2}\;=\;\int\bigl{[}f(\bm{r})\bigr{]}^{*}\,w(\bm{r})\,g(\bm{r})\,\mathrm{d}^{3}\bm{r}\,.$$

$$\Vert f{\Vert }_{w,2}=\sqrt{(f|f{)}_{w,2}},$$

$$L_w^2({\mathbb{R}}^3)$$

$$\displaystyle\;=\;\bigl{\{}f\colon\mathbb{R}^{3}\to\mathbb{C}\big{|}\,\|f\|_{w,2}<\infty\bigr{\}}\,.$$

$$\langle f|g\rangle_{2,1}\;=\;\int\bigl{[}f(\bm{r})\bigr{]}^{*}\,\frac{\eta^{2}-\nabla^{2}}{2\eta^{2}}\,g(\bm{r})\,\mathrm{d}^{3}\bm{r}\,.$$

$$⟨⟨f⟩{⟩}_{2,1}=\sqrt{⟨f|f{⟩}_{2,1}}.$$

$$W_2^{(1)}({\mathbb{R}}^3)$$

$$\displaystyle\;=\;\Bigl{\{}f\colon\mathbb{R}^{3}\to\mathbb{C}\Bigm{|}\,\int\bigl{[}f(\bm{r})\bigr{]}^{*}\,\frac{\eta^{2}-\nabla^{2}}{2\eta^{2}}\,f(\bm{r})\,\mathrm{d}^{3}\bm{r}<\infty\Bigr{\}}$$

$$\displaystyle\;=\;\bigl{\{}f\colon\mathbb{R}^{3}\to\mathbb{C}\big{|}\,\langle\!\langle f\rangle\!\rangle_{2,1}<\infty\bigr{\}}\,.$$

$$f(\mathbf{r}\pm {\mathbf{r}}^{\prime })$$

$$C_{n,\ell }^m(f;\pm {\mathbf{r}}^{\prime })$$

$$\sum _{n\ell m}=\sum _{n=1}^{\infty }\sum _{\ell =0}^{n-1}\sum _{m=-\ell }^{\ell },$$

$$C_{n,\ell }^m(f;\pm {\mathbf{r}}^{\prime })$$

$$T_{n^{\prime }{\ell }^{\prime }{m}^{\prime }}^{n\ell m}(f;\pm )$$

$$f(\mathbf{r}\pm {\mathbf{r}}^{\prime })=\sum _{\begin{array}{c}n\ell m\\ n^{\prime }{\ell }^{\prime }{m}^{\prime }\end{array}}{T}_{n^{\prime }{\ell }^{\prime }{m}^{\prime }}^{n\ell m}(f;\pm ){\phi }_{n,\ell }^m(\mathbf{r}){\phi }_{n^{\prime },{\ell }^{\prime }}^{m^{\prime }}({\mathbf{r}}^{\prime }).$$

$$f(\mathbf{r}\pm {\mathbf{r}}^{\prime })$$

$$=\sum _{\begin{array}{c}n\ell m\\ n^{\prime }{\ell }^{\prime }{m}^{\prime }\end{array}}{\mathbf{T}}_{n^{\prime }{\ell }^{\prime }{m}^{\prime }}^{n\ell m}(f,w;\pm ){\psi }_{n,\ell }^m(\mathbf{r}){\psi }_{n^{\prime },{\ell }^{\prime }}^{m^{\prime }}({\mathbf{r}}^{\prime }).$$

$$f(\mathbf{r}\pm {\mathbf{r}}^{\prime })=\sum _{n\ell m}{\mathbf{C}}_{n,\ell }^m(f,w;\pm {\mathbf{r}}^{\prime }){\psi }_{n,\ell }^m(\mathbf{r}),$$

$${\mathbf{C}}_{n,\ell }^m(f,w;\pm {\mathbf{r}}^{\prime })$$

$$\displaystyle\qquad\;=\;\int\,\bigl{[}\psi_{n,\ell}^{m}(\bm{r})\bigr{]}^{*}\,w(\bm{r})\,f(\bm{r}\pm\bm{r}^{\prime})\,\mathrm{d}^{3}\bm{r}\,,$$

Full text

Keywords: Coulomb potential, addition theorems, exponentially decaying functions, orthogonal and nonorthogonal expansions, Hilbert spaces, generalized functions, weak convergence.

1. 1 Introduction
2. 2 Basic Hilbert Space Theory
3. 3 One-Range Addition Theorems
4. 4 Laguerre-Type Function Sets and Their Hilbert Spaces
5. 5 Addition Theorems for Exponentially Decaying Functions
6. 6 Guseinov’s One-Range Addition Theorems for Slater-Type Functions
7. 7 Weakly Convergent One-Range Addition Theorems for the Coulomb Potential
8. 8 Differentiation Techniques and the Spherical Tensor Gradient Operator
9. 9 Summary and Conclusions
10. A Spherical Harmonics and Gaunt Coefficients

1 Introduction

In many subfields of physics and physical chemistry, an essential step towards a solution of a problem consists in a separation of variables. Principal mathematical tools, that can accomplish such a separation of variables, are so-called addition theorems. These are expansions of a given function $f(\bm{r}\pm\bm{r}^{\prime})$ with $\bm{r},\bm{r}^{\prime}\in\mathbb{R}^{3}$ in products of other functions that only depend on either $\bm{r}$ or $\bm{r}^{\prime}$. A review of the relevant literature with an emphasis on electronic structure calculations can be found in [1, 2, 3]. Applications of addition theorems in classical physics are described in books by Martin [4] and Jones [5].

In atomic or molecular calculations, we are predominantly interested in irreducible spherical tensors

[TABLE]

that can be represented as products of a radial function $f_{\ell}(r)$ multiplied by a surface spherical harmonic $Y_{\ell}^{m}(\theta,\phi)$ defined by (A). The convenient orthonormality and completeness properties of spherical harmonics make it highly desirable that the functions of either $\bm{r}$ or $\bm{r}^{\prime}$, which occur in the expansion of $f(\bm{r}\pm\bm{r}^{\prime})$, are also irreducible spherical tensors of the type of (1.1). Thus, addition theorems are expansions in terms of products of spherical harmonics with arguments $\theta,\phi=\bm{r}/r$ and $\theta^{\prime},\phi^{\prime}=\bm{r}^{\prime}/r^{\prime}$, respectively.

The best known example of such an addition theorem is the Laplace expansion of the Coulomb potential:

[TABLE]

Here, ${\mathcal{Y}}_{\ell}^{m}$ and ${\mathcal{Z}}_{\ell}^{m}$ denote regular and irregular solid harmonics defined by (A.2) and (A.3), respectively.

The Laplace expansion (1) leads to a separation of the variables ${\bm{r}}$ and ${\bm{r}}^{\prime}$ and decouples integration variables in expectation values of the Coulomb potential. However, the right-hand side of (1) depends on ${\bm{r}}$ and ${\bm{r}}^{\prime}$ only indirectly via the two vectors $\bm{r}_{<}$ and $\bm{r}_{>}$ that satisfy $|\bm{r}_{<}|<|\bm{r}_{>}|$. Hence, the Laplace expansion has a two-range form, depending on the relative length of $\bm{r}$ and $\bm{r}^{\prime}$. This is a complication, which occurs frequently among addition theorems and which can lead to nontrivial technical problems.

As discussed in more details in [1, 2, 3], an addition theorem for $f(\bm{r}\pm\bm{r}^{\prime})$ can be derived by performing a three-dimensional Taylor expansion (see for example [6, p. 181]):

[TABLE]

Thus, the translation operator

[TABLE]

generates $f(\bm{r}+\bm{r}^{\prime})$ by performing a three-dimensional Taylor expansion of $f$ around $\bm{r}$ with shift vector $\bm{r}^{\prime}$. Since the variables $\bm{r}$ and $\bm{r}^{\prime}$ are separated, the series expansion (1.3) is indeed an addition theorem.

We could also expand $f$ around $\bm{r}^{\prime}$ and use $\bm{r}$ as the shift vector. This would produce an addition theorem for $f(\bm{r}\pm\bm{r}^{\prime})$ in which the roles of $\bm{r}$ and $\bm{r}^{\prime}$ are interchanged. Both approaches are mathematically legitimate and equivalent if $f$ is analytic in the sense of complex analysis at $\bm{r}$, $\bm{r}^{\prime}$, and $\bm{r}\pm\bm{r}^{\prime}$ for essentially arbitrary vectors $\bm{r},\bm{r}^{\prime}\in\mathbb{R}^{3}$. Unfortunately, this is normally not true. Most functions of interest in atomic and molecular electronic structure theory are either singular or not analytic at the origin $\bm{r}=(0,0,0)$. Obvious examples are the Coulomb potential, which is singular at the origin, or the $1s$ hydrogen eigenfunction, which possesses a cusp at the origin. In fact, all exponentially decaying function sets discussed in this article are not analytic at the origin.

Nonanalyticity is not a mathematical sophistication that can safely be ignored in practical applications. The radius of convergence of the Taylor expansion of a function is determined by the location of its nearest singularity or pole. If a function $f$ is not analytic at the origin, an addition theorem derived via (1.3) can only converge if $|\bm{r}|>|\bm{r}^{\prime}|$, which gives rise to a two-range form.

A two-range form can be quite inconvenient. The use of such an addition theorem in a multicenter integral ultimately leads to indefinite integrals of special functions, whose efficient and reliable evaluations can be very difficult (I struggled with these problems long ago in my diploma thesis [7], which was published in condensed form in [8]). In contrast, a $1s$ Gaussian $\exp(-\beta r^{2})$ with $\beta>0$ is analytic for all $\bm{r}\in\mathbb{R}^{3}$. This implies that its addition theorem, which can be obtained via (1.3), has a one-range form (see for example [9, Eq. (9)] or also [10, Eq. (17)]).

The fact that addition theorems derived via (1.3) are just rearranged three-dimensional Taylor series automatically implies that they converge pointwise.

It is one of the central results of modern mathematics that convergence is not an abstract or global property. It depends very much on how we measure it. Pointwise convergence is a very demanding requirement. Moreover, it is not really necessary if we only want to use addition theorems in multicenter integrals. Therefore, it makes sense to wonder, whether computational benefits can be achieved by relaxing the requirement of pointwise convergence and by demanding instead a weaker form of convergence.

This is indeed possible. Often, one-range addition theorems can be constructed that are essentially expansions of $f(\bm{r}\pm\bm{r}^{\prime})$ in terms of complete orthonormal function sets. Such an addition theorem converges in the mean with respect to the norm of the corresponding Hilbert space. This had already been done by several authors (see for example [11, 12, 13, 14, 15, 16, 17, 18, 19] and references therein). Guseinov also derived and applied one-range addition theorems in several articles [20, 21, 22].

In recent years, Guseinov and his coworkers made an industry out of one-range addition theorems of exponentially decaying functions. They published – in addition to several other articles on different topics – a very long list of articles on the derivation and application of one-range addition theorems [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 63, 69], that are often remarkably similar and that do not give due credit to the previous work of others.

At first sight, this long list of articles looks quite impressive. Unfortunately, from a mathematical point of view Guseinov’s treatment of one-range addition theorems is at best questionable and in some cases fundamentally flawed. Guseinov and his coworkers ignored basic facts of Hilbert space and approximation theory as well as all questions of convergence and existence. Such a cavalier attitude is not uncommon among scientists, but it is completely unacceptable in the context of addition theorems and multicenter integrals which are essentially a mathematical topic.

As discussed in this article, Guseinov’s approach leads to serious problems, which are particularly evident in the case of Guseinov’s derivation of a one-range addition theorem for the Coulomb potential [39]. However, as discussed in more details later, serious mathematical problems persist also in all other articles dealing with one-range addition theorem. Therefore, the title of this article is somewhat misleading. In this article, I will try to provide a reasonably balanced and detailed discussion of one-range addition theorems. Accordingly, my criticism is not limited to [39] and extends to all articles by Guseinov and coworkers dealing with one-range addition theorems [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 63, 69].

Guseinov’s starting point for his one-range addition theorems of the Coulomb potential was a class of one-range addition theorems [23, 24, 25, 26] for Slater-type functions which in unnormalized form can be expressed as follows:

[TABLE]

Here, $N$ is some kind of generalized principal quantum number, $L$ and $M$ are the usual angular momentum quantum numbers, and $\beta$ is a positive scaling parameter. In the vast majority of articles dealing with Slater-type functions, $N$ is assumed to be a positive integer $\geq 1$, i.e., $N\in\mathbb{N}$. However, several authors – among them Guseinov and coworkers [24, 26, 70, 27, 29, 30, 32, 33, 44, 53, 54, 55, 56, 57, 58, 59, 60, 61, 43, 66, 70, 71, 72, 73, 74, 75, 66] – considered also nonintegral principal quantum numbers $N\in\mathbb{R}\setminus\mathbb{N}$.

Guseinov’s one-range addition theorems for Slater-type functions with nonintegral principal quantum numbers [24, 25, 26] contain as the special case $N=L=M=0$ addition theorems for the Yukawa potential $\exp(-\beta r)/r$ [76], which may be viewed to be an exponentially screened Coulomb potential. Guseinov derived his one-range addition theorems for the Coulomb potential by exploiting the obvious relationship $1/r=\lim_{\beta\to 0}\exp(-\beta r)/r$ in the one-range addition theorems for the Yukawa potential.

At first sight, Guseinov’s derivation of a class of one-range addition theorems for the Coulomb potential looks like a remarkable achievement. So far, the Laplace expansion (1) has routinely been used in atomic and molecular electronic structure calculations, but its characteristic two-range form can easily lead to nontrivial technical problems. Obviously, a one-range addition theorem for the Coulomb potential should simplify subsequent integrations in multicenter integrals substantially.

As already remarked above, one-range addition theorems are essentially expansions of a given function $f(\bm{r}\pm\bm{r}^{\prime})$ belonging to a suitable Hilbert space in terms of functions that are complete and orthonormal in this Hilbert space. It is an essential feature of orthogonal expansions in general and of one-range addition theorems in special that they only converge in the mean with respect to the norm of the underlying Hilbert space, but not necessarily pointwise (see for example [77]).

The fact, that orthogonal expansions only converge in the mean if the function, which is to be expanded, is an element of the corresponding Hilbert space, has an obvious consequence: Since the Coulomb potential does not belong to any of the Hilbert spaces, which Guseinov implicitly used and which all involve an integration over the whole $\mathbb{R}^{3}$, Guseinov’s one-range addition theorems for the Coulomb potential all diverge in the mean.

But this is not the only problem. For the derivation of his addition theorems, Guseinov had first expanded Slater-type functions $\chi_{N,L}^{M}(\beta,\mathbf{r}\pm\bm{r}^{\prime})$ with arbitrary principal quantum numbers $N\in\mathbb{R}$ in terms of the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ defined by (4), yielding the addition theorems (6.1). The functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k=-1,0,1,2,\dots$, which Guseinov had introduced in [24, Eq. (1)], are – as discussed in Section 4 – complete and orthonormal in the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$ defined by (4). As discussed in more details in Section (4), Guseinov’s functions generalize some other, well established function sets also based on the generalized Laguerre polynomials.

For some reasons, which had never had been explained comprehensively and which I thus do not really understand, Guseinov considered it to be advantageous to replace in his addition theorems (6.1) his complete and orthonormal functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ by nonorthogonal Slater-type functions with integral principal quantum numbers and to rearrange the order of summations of his expansions. In this way, Guseinov obtained the expansions (6) of Slater-type functions with in general nonintegral principal quantum numbers in terms of Slater-type functions with integral principal quantum numbers located at a different center.

Due to the intrinsic complexity of addition theorems, it is by no means easy to decide whether such a rearrangement is legitimate and whether the resulting expansions (6) are mathematically meaningful. In general, this is an open question, but if the principal quantum number of $\chi_{N,L}^{M}(\beta,\mathbf{r}\pm\bm{r}^{\prime})$ is nonintegral, $N\in\mathbb{R}\setminus\mathbb{N}$, an affirmative answer is possible: The one-center limit $\bm{r}^{\prime}=\bm{0}$ in (6) does not exist, which means that Guseinov’s rearranged addition theorems (6) do not exist for the whole argument set $\mathbb{R}^{3}\times\mathbb{R}^{3}$.

This conclusion is highly consequential: The rearranged addition theorems (6) with $N=L=M=0$ were the starting point for Guseinov’s derivation of one-range addition theorems for the Coulomb potential. Thus, Guseinov obtained his one-range addition theorems for the Coulomb potential via addition theorems for the Yukawa potential that do not exist for the whole argument set.

Guseinov was not the first one who had derived a divergent expansion of the Coulomb potential in terms of a complete and orthonormal function set. Salmon, Birss, and Ruedenberg [78] derived a bipolar expansion of the Coulomb potential in terms of the Gaussian-type eigenfunctions of a three-dimensional isotropic harmonic oscillator defined by (4) which are complete and orthonormal in the Hilbert space $L^{2}(\mathbb{R}^{3})$ of square integrable functions defined by (2) (see for instance [18, Section V]). However, Silverstone and Kay [79] demonstrated that this bipolar expansion diverges. This observation was later confirmed by Ruedenberg and Salmon [80].

Because of their divergence, one might dismiss Guseinov’s one-range addition theorems of the Coulomb potential to be practically useless, just as the expansion of Salmon, Birss, and Ruedenberg [78] had been dismissed. However, this would be premature and the situation is actually more complicated but also less hopeless than it may appear at first sight. As discussed in Section 7, these addition theorems can be interpreted as weakly convergent expansions of generalized functions, that are mathematically meaningful and yield convergent results when used in suitable functionals.

It is quite obvious that Guseinov had failed to understand the mathematical theory behind one-range addition theorems. Thus is bad enough. However, the referees of Guseinov’s numerous recent articles on one-range addition theorems and related topics [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 63, 69] apparently also failed to understand this theory. Therefore, I will try to give a compact, but hopefully comprehensive treatment of one-range addition theorems and of the mathematical tools, which are needed for their derivation, interpretation, and application.

Section 2 gives a compact review of the for our purposes most fundamental properties of Hilbert spaces in the context of approximation theory and of orthogonal expansions. Section 3 discusses the derivation and the basic properties of one-range addition theorems. In Section 4, complete and orthonormal Laguerre-type function sets and their corresponding Hilbert spaces are discussed. Section 5 discusses the derivation of one-range addition theorems for Slater-type functions in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ using techniques that were applied by Filter and Steinborn [12] for the derivation of related addition theorems. Section 6 discusses Guseinov’s highly questionable derivation of one-range addition theorems for Slater-type functions. Section 7 discusses how divergent one-range addition theorems of the Coulomb potential can be interpreted as weakly convergent expansions that are mathematically meaningful in suitable functionals. Section 8 discusses alternatives to the differentiation techniques, which Guseinov had used for the generation of new one-range addition theorems and which are neither convenient from a technical point of view nor mathematically justified. This article is concluded by a short summary in Section 9. Finally, there is Appendix A listing some basic facts about spherical harmonics and Gaunt coefficients.

2 Basic Hilbert Space Theory

Let us assume that $\mathcal{V}$ is a vector space over the complex numbers $\mathbb{C}$ that possesses an inner product $(\cdot|\cdot)\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}$ (see for instance [81, p. 36]). If $\mathcal{V}$ is complete with respect to the norm $\|\cdot\|\colon\mathcal{V}\to\mathbb{R}_{+}$ defined by

[TABLE]

i.e., if every Cauchy sequence in $\mathcal{V}$ converges with respect to (2.1) to an element of $\mathcal{V}$, then $\mathcal{V}$ is called a Hilbert space. In the mathematical literature, convergence in the means with respect to the norm $\|\cdot\|$ is frequently called strong convergence, in contrast to weak convergence which will be discussed later in Section 7.

Hilbert spaces play a major role in various branches of mathematics and mathematical physics and in particular also in approximation theory. If $f$ is an element of some Hilbert space $\mathcal{H}$ and if $\{\varphi_{n}\}_{n=0}^{\infty}$ is linearly independent and complete in $\mathcal{H}$, then we can construct approximations

[TABLE]

to $f$, where $N$ is a finite integer. The coefficients $C_{n}^{(N)}$ are chosen in such a way that the mean square deviation $\|f-f_{N}\|^{2}=(f-f_{N}|f-f_{N})$ becomes minimal.

It is one of the central results of approximation theory that $\|f-f_{N}\|^{2}$ becomes minimal if $\{\varphi_{n}\}_{n=0}^{\infty}$ is an orthonormal function set satisfying $(\varphi_{n}|\varphi_{n^{\prime}})=\delta_{n{n^{\prime}}}$ for all indices $n$ and $n^{\prime}$ and if the coefficients are chosen according to $C_{n}^{(N)}=(\varphi_{n}|f)$ (see for example [82, Theorem 9 on p. 51]). Since the coefficients $(\varphi_{n}|f)$ do not depend on the truncation order $N$ in $f_{N}$, $f$ possesses an expansion

[TABLE]

in terms of the orthogonal function set $\{\varphi_{n}\}_{n=0}^{\infty}$ that converges in the mean with respect to the norm $\|\cdot\|$ of the Hilbert space $\mathcal{H}$.

This result emphasizes the central role of orthogonal expansions both in Hilbert spaces as well as in approximation theory. If the functions $\{\varphi_{n}\}_{n=0}^{\infty}$ are only normalized satisfying $(\varphi_{n}|\varphi_{n})=1$ but not orthogonal, then it is in general only possible to construct finite approximations $f_{N}$ of the type of (2.2) by minimizing the mean square deviation $\|f-f_{N}\|^{2}=(f-f_{N}|f-f_{N})$, but we cannot tacitly assume that the coefficients $C_{n}^{(N)}$ in $f_{N}$ converge to a well defined limit $C_{n}=C_{n}^{(\infty)}$ as $N\to\infty$. Accordingly, the existence of an expansion

[TABLE]

is not guaranteed if the functions $\{\varphi_{n}\}_{n=0}^{\infty}$ are not orthogonal. This fact is well documented both in the mathematical literature (see for example [82, Theorem 10 on p. 54] or [83, Section 1.4]) as well as in the literature on electronic structure calculations [84, 85, 86, 87, 88, 89]).

Thus, completeness of a function set $\{\varphi_{n}\}_{n=0}^{\infty}$ in a Hilbert space $\mathcal{H}$ does not suffice to guarantee the existence of formal expansions in terms of these functions: These expansions may or they may not exist. Horrifying examples of pathologies of nonorthogonal expansions can be found in [88, Section III.I].

There is another, practically very consequential aspect of orthogonal expansions of the type of (2.3); The coefficients $(\varphi_{n}|f)$ of orthogonal expansions satisfy Parseval’s equality (see for example [81, Eq. (II.2) on p. 45])

[TABLE]

Thus, the coefficients $(\varphi_{n}|f)$ are bounded in magnitude and they have to vanish as $n\to\infty$. This may well be the main reason why orthogonal expansions tend to be computationally well behaved. In the case of nonorthogonal expansions of the type of (2.4) it can instead happen that the expansion coefficients $C_{n}$ are unbounded, have alternating signs, and increase in magnitude with increasing index (see for example [11, Appendix E on pp. 162 - 164] or [88, Table I on p. 166]). Such a behavior can easily lead to a cancellation of significant digits and to a catastrophic accumulation of rounding errors.

If we remove a single function from a complete and orthonormal function set, it becomes incomplete. In the case of nonorthogonal expansions of the type of (2.4), the situation is much more complicated since nonorthogonal function sets are in general overcomplete as well as almost linearly dependent (see for example [87, 88] and references therein). Again, this can be the source of serious numerical problems.

Of course, there are situations in which nonorthogonal expansions offer computational advantages (see fore example the discussion in [90]). However, in the vast majority of all cases, orthogonal expansions are clearly superior. Consequently, one should should not voluntarily abandon the highly useful feature of orthogonality unless there are truly compelling reasons.

The natural inner product for effective one-particle wave functions $f,g\colon\mathbb{R}^{3}\to\mathbb{C}$ in atomic and molecular electronic structure calculations on the basis of the Hartree-Fock-Roothaan equations [91, 92, 93] is given by the integral

[TABLE]

As always in this article, integration extends over the whole $\mathbb{R}^{3}$. The corresponding norm satisfies

[TABLE]

Accordingly, the Hilbert space

[TABLE]

of square-integrable one-particle wave functions provides the natural setting for atomic and molecular electronic structure calculations in general as well as for the derivation of one-range addition theorems in special.

The whole formalism of inner products, norms, and function spaces can be generalized to include weight functions $w\colon\mathbb{R}^{3}\to\mathbb{R}_{+}$. If $w(\bm{r})\geq 0$ is such a positive weight function, we can define an inner product with respect to the weight function $w$ for functions $f,g\colon\mathbb{R}^{3}\to\mathbb{C}$ according to

[TABLE]

On the basis of the inner product (2.9), the norm of a function $f\colon\mathbb{R}^{3}\to\mathbb{C}$ with respect to the weight function $w$ is defined according to

[TABLE]

and the Hilbert space $L_{w}^{2}(\mathbb{R}^{3})$ of square integrable functions with respect to the weight function $w$ is defined via the norm (2.10) according to

[TABLE]

By augmenting the inner product (2.6) by a nontrivial weight function $w(\bm{r})\neq 1$ according to (2.9), it should at least in principle be possible to accomplish some fine-tuning in approximation processes. So, it looks like an obvious idea to use instead of $L^{2}(\mathbb{R}^{3})$ a weighted Hilbert $L_{w}^{2}(\mathbb{R}^{3})$ that is better adapted to the problem under consideration.

Unfortunately, there are some principal problems which must not be ignored. In general, we neither have $L_{w}^{2}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$ nor $L^{2}(\mathbb{R}^{3})\subset L_{w}^{2}(\mathbb{R}^{3})$. Thus, the two Hilbert spaces $L^{2}(\mathbb{R}^{3})$ and $L_{w}^{2}(\mathbb{R}^{3})$ are mathematically inequivalent. Accordingly, the expansion of a function $f\in L^{2}(\mathbb{R}^{3})$ in terms of a function set, that is complete and orthonormal in $L_{w}^{2}(\mathbb{R}^{3})$, does not necessarily lead to a convergent result. Moreover, expansions of a given function in terms of function sets, that are complete and orthonormal in either $L^{2}(\mathbb{R}^{3})$ or in $L_{w}^{2}(\mathbb{R}^{3})$, respectively, can also be computationally different, since they can have substantially different rates of convergence (see the simple example discussed in Section 4).

Even the Hilbert space $L^{2}(\mathbb{R}^{3})$ of square integrable functions is – loosely speaking – too large for electronic structure calculations on the basis of the Hartree-Fock-Roothaan equations. All effective one-particle wave functions must belong to $L^{2}(\mathbb{R}^{3})$, but the converse is not necessarily true, i.e., there are elements of $L^{2}(\mathbb{R}^{3})$ that cannot be used in electronic structure calculations. An example is Yukawa potential $\exp(-\beta r)/r$. It belongs to $L^{2}(\mathbb{R}^{3})$, but it cannot be used as a trial function in electronic structure calculations since the expectation values of the kinetic energy and of the Coulomb potential do not exist.

The necessary exclusion of unsuitable elements from $L^{2}(\mathbb{R}^{3})$ can be accomplished via some other generalizations of $L^{2}$ spaces, the so-called Sobolev spaces (see for example [94, 95, 96]), which are of considerable importance not only in electronic structure theory. The in the context of electronic structure calculations most important Sobolev space $W_{2}^{(1)}(\mathbb{R}^{3})$ can be defined via the following inner product for differentiable functions $f,g\colon\mathbb{R}^{3}\to\mathbb{C}$:

[TABLE]

Here, \nabla=\bigl{(}\partial/\partial x,\partial/\partial y,\partial/\partial z\bigr{)} is the three-dimensional gradient operator, and $\eta$ is a real constant. As usual, integration extends over the whole $\mathbb{R}^{3}$.

On the basis of the inner product (2.12), the following Sobolev-type norm of a differentiable function $f\colon\mathbb{R}^{3}\to\mathbb{C}$ can be defined:

[TABLE]

Both the inner product $\langle f|g\rangle_{2,1}$ as well as the norm $\langle\!\langle f\rangle\!\rangle_{2,1}$ depend on the real parameter $\eta$. Since it can be shown that the norms defined via (2.13) are equivalent for all $\eta\in\mathbb{R}\setminus{0}$ and thus give rise to the same topology, their dependence on $\eta$ is not explicitly indicated.

The Sobolev space $W_{2}^{(1)}(\mathbb{R}^{3})$ is defined via the norm (2.13) according to

[TABLE]

Obviously, the Sobolev space $W_{2}^{(1)}(\mathbb{R}^{3})$ is also a Hilbert space. In addition, $W_{2}^{(1)}(\mathbb{R}^{3})$ is also a proper subspace of the Hilbert space $L^{2}(\mathbb{R}^{3})$, i.e., $W_{2}^{(1)}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$. Nevertheless, approximation processes in $L^{2}(\mathbb{R}^{3})$ and $W_{2}^{(1)}(\mathbb{R}^{3})$, respectively, can differ substantially. Loosely speaking, we can say that an approximation scheme in $L^{2}(\mathbb{R}^{3})$ approximates $f$ in the mean, but in $W_{2}^{(1)}(\mathbb{R}^{3})$ both $f$ as well as $\nabla f$ are approximated in the mean. Thus, in Sobolev spaces certain possible pathologies can be avoided. This can have far-reaching consequences in quantum mechanical calculations which are essentially special approximation processes. In [85, Section 9], it was shown that completeness of a one-particle basis in $L^{2}(\mathbb{R}^{3})$ does not suffice to guarantee the convergence of computations based on the Rayleigh-Ritz variational principle. This is only guaranteed if the one-particle basis is complete in $W_{2}^{(1)}(\mathbb{R}^{3})$ (see also [97, 98]).

3 One-Range Addition Theorems

Let us assume that $f\in L^{2}(\mathbb{R}^{3})$ and that the functions $\{\varphi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ are complete and orthonormal in $L^{2}(\mathbb{R}^{3})$. It makes sense to assume that the functions $\{\varphi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ are irreducible spherical tensors of the type of (1.1). Thus, the index $n$ can be viewed to be some kind of generalized principal quantum number, and $\ell$ and $m$ are the usual angular momentum quantum numbers.

An addition theorem for $f(\bm{r}\pm\bm{r}^{\prime})$, which converges in the mean with respect to the norm of $L^{2}(\mathbb{R}^{3})$, can be constructed by expanding $f$ in terms of the orthonormal functions $\{\varphi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$:

[TABLE]

The summation limit in (3.1) depend on the exact definition of the function set $\{\varphi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$. This article always uses the convention

[TABLE]

which is in agreement with the usual convention for the bound state hydrogen eigenfunctions.

The expansion (3.1) is a one-range addition theorem since the variables ${\bm{r}}$ and ${\bm{r}}^{\prime}$ are completely separated: The dependence on $\bm{r}$ is entirely contained in the functions $\varphi_{n,\ell}^{m}(\bm{r})$, whereas $\bm{r}^{\prime}$ occurs only in the expansion coefficients $C_{n,\ell}^{m}(f;\pm\bm{r}^{\prime})$ which are overlap integrals.

If the overlap integrals $C_{n,\ell}^{m}(f;\pm\bm{r}^{\prime})$ can be expanded in terms of the functions $\varphi_{n,\ell}^{m}(\bm{r}^{\prime})$ according to

[TABLE]

then the addition theorem (3.1) assumes a completely symmetrical form:

[TABLE]

From a purely formal point of view, the derivation of one-range addition theorems is a triviality. The challenging part is the construction of computationally convenient mathematical expressions for the overlap integrals $C_{n,\ell}^{m}(f;\pm\bm{r}^{\prime})$ in (3.1) or the coefficients $T_{n^{\prime}\ell^{\prime}m^{\prime}}^{n\ell m}(f;\pm)$ in (3.4). In realistic applications, we cannot not tacitly assume that the use of one-range addition theorems in a multicenter integral necessarily leads to rapidly convergent expansions. Therefore, we must be able to compute $C_{n,\ell}^{m}(f;\pm\bm{r}^{\prime})$ and/or $T_{n^{\prime}\ell^{\prime}m^{\prime}}^{n\ell m}(f;\pm)$ efficiently and reliably even for possibly very large indices (see for instance the convergence rates reported in [99]).

Symmetrical one-range addition theorems of the kind of (3.4), which converge in the mean with respect to the norm of the Hilbert space $L^{2}(\mathbb{R}^{3})$, were constructed by Filter and Steinborn [12, Eqs. (5.11) and (5.12)] and later applied by Kranz and Steinborn [100] and by Trivedi and Steinborn [99].

As already remarked above, a nontrivial weight function $w(\bm{r})\neq 1$ in an inner product can give more weight to those regions of space in which $f$ is large, while deemphasizing the contribution from those regions in which $f$ is small. Accordingly, the inclusion of a suitable weight function $w$ can improve convergence. It is thus an in principle obvious idea to construct one-range addition theorems that converge with respect to the norm of a suitable weighted Hilbert space $L_{w}^{2}(\mathbb{R}^{3})$ defined in (2). Let us therefore assume that $f$ belongs to $L_{w}^{2}(\mathbb{R}^{3})$. Then we can construct a one-range addition theorem by expanding $f(\bm{r}\pm\bm{r}^{\prime})$ with respect to a function set $\{\psi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ that is complete in $L_{w}^{2}(\mathbb{R}^{3})$ and orthonormal with respect to the modified inner product (2.9):

[TABLE]

For the derivation of (3), we only have to replace (3.1) by

[TABLE]

and (3.3) by

[TABLE]

As discussed in Section 4, Guseinov and coworkers [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 50, 47, 48, 49, 51, 53, 54, 55, 56, 57, 59, 60, 61, 62] used for the derivation of one-range addition theorems and of related mathematical objects predominantly the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ defined in (4) that are complete and orthonormal in the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$ defined in (4). As long as the functions, which are to be expanded, belong to the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$, this is a completely legitimate approach that leads to addition theorems of the type of (3) which converge with respect to the norm (4.19) of $L_{r^{k}}^{2}(\mathbb{R}^{3})$.

Unfortunately, the discussion in Section 4 shows that the use of a weight function $r^{k}$ with $k\geq 1$ does not necessarily lead to better results. Moreover, it is very difficult or even practically impossible to decide on the basis of simple a priori considerations for which $k$ the weight function $r^{k}$ produces best results. Thus, it is by no means clear how we can actually profit from this additional degree of freedom.

There is also another annoying practical problem which can easily occur in multicenter integrals: Let us assume that $f(\bm{r}\pm\bm{r}^{\prime})$ is expanded in terms of a function set $\{\psi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ that is complete in a weighted Hilbert space $L_{w}^{2}(\mathbb{R}^{3})$ and orthonormal with respect to the corresponding inner product (2.9). If we use this one-range addition theorem in a multicenter integral containing $f(\bm{r}\pm\bm{r}^{\prime})$, the weight function $w(r)$, which makes $\{\psi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ orthonormal, normally does not occur there. Consequently, subsequent integrations involving $f(\bm{r}\pm\bm{r}^{\prime})$ are not necessarily simplified by orthogonality. Thus, it may be more difficult to apply a one-range addition theorems of the type of (3) based on a nontrivial weight function $w(\bm{r})\neq 1$ than an addition theorem of the type of (3.4). It is not the complexity of an addition theorem, that really matters, but the complexity of the series expansion, which we obtain by inserting an addition theorem into a multicenter integral.

If we compare one-range and two-range addition theorems, it is obvious that one-range addition theorems greatly facilitate subsequent integrations. Thus, they seem to be clearly superior to two-range addition theorems of the type of the Laplace expansion (1).

However, one-range addition theorems also have disadvantages, and a balanced assessment of the relative merits of one-range and two-range addition theorems is by no means easy. Firstly, one-range addition theorems usually have a more complicated structure than two-range addition theorems (typically, they contain one additional infinite summation). Secondly, one-range addition theorems normally do not converge pointwise, but only in the mean with respect to the norm of the underlying Hilbert space. This can be quite advantageous, since it makes it possible to expand functions with singularities and/or discontinuities. If, however, singularities and/or discontinuities are present, then a nonsmooth function has to be approximated by smooth functions. Very often, this leads to slow convergence.

The probably most severe disadvantage of one-range addition theorems compared to two-range addition theorems is that the approach sketched above cannot be applied to all functions of interest. As remarked in Section 2, inner products for effective one-particle wave functions in atomic and molecular electronic structure calculations all involve integrations over the whole $\mathbb{R}^{3}$. Thus, an effective one-particle wave function $f$ has to belong to the Hilbert space $L^{2}(\mathbb{R}^{3})$ or also to a suitable weighted Hilbert space $L_{w}^{2}(\mathbb{R}^{3})$. Many functions of considerable relevance such as the Coulomb potential or the irregular solid harmonic $\mathcal{Z}_{\ell}^{m}(\bm{r})$ are only locally integrable, not square integrable with respect to an integration over the whole $\mathbb{R}^{3}$ and thus do not belong to $L^{2}(\mathbb{R}^{3})$ or to any other weighted Hilbert spaces discussed in this article. If we nevertheless want to formulate one-range addition theorems for these functions by expanding them in terms of a complete and orthonormal function set, we are faced with the problem that these one-range addition theorems diverge in the mean.

4 Laguerre-Type Function Sets and Their Hilbert Spaces

The surface spherical harmonics $Y_{\ell}^{m}(\theta,\phi)$ are complete and orthonormal with respect to an integration over the surface of the unit sphere in $\mathbb{R}^{3}$ (an explicit proof can for instance be found in [101, Section III.7.6]). Since more complex Hilbert spaces can be constructed by forming tensor products of simpler Hilbert spaces (see for example [101, Section II.6.5]), we only have to find suitable radial functions that are complete and orthogonal with respect to an integration from [math] to $\infty$ (see also [86, Lemma 6 on p. 31]). Thus, we more or less automatically arrive at function sets based on the generalized Laguerre polynomials.

The generalized Laguerre polynomials $L_{n}^{(\alpha)}(x)$ with $\mathrm{Re}(\alpha)>-1$ are orthogonal polynomials associated with the integration interval $[0,\infty)$ and the weight function $w(x)=x^{\alpha}\mathrm{exp}^{-x}$. They possess the following explicit expressions [102, p. 240]:

[TABLE]

These expressions can be used to extend the definition of $L_{n}^{(\alpha)}(x)$ to complex values of $\alpha$. The generalized Laguerre polynomials can also be defined via the following Rodrigues relationship [102, p. 241]:

[TABLE]

The generalized Laguerre polynomials satisfy for $\mathrm{Re}(\alpha)>-1$ and $m,n\in\mathbb{N}_{0}$ the following orthogonality relationship [102, p. 241]:

[TABLE]

Henceforth, the condition $\mathrm{Re}(\alpha)>-1$, which is necessary for the existence of this and related integrals, will always be tacitly assumed.

The completeness of the generalized Laguerre polynomials in the weighted Hilbert space

[TABLE]

is a classic result of mathematical analysis (see for example [83, p. 33], [103, pp. 349 - 351], or [104, pp. 235 - 238]).

A different convention for Laguerre polynomials is frequently used in the quantum mechanical literature. For example, Bethe and Salpeter [105, Eq. (3.5)] define associated Laguerre functions \bigl{[}L_{n}^{m}(x)\bigr{]}_{\text{BS}} with $n,m\in\mathbb{N}_{0}$ via the Rodrigues-type relationships

[TABLE]

Comparison of (4.2) and (4.5b) implies:

[TABLE]

The convention of Bethe and Salpeter [105] is also used in the books by Condon and Shortley [106, Eqs. (6) and (9) on p. 115] and by Condon and Odabaşi [107, Eq. (2) on p. 189] as well as in numerous articles.

In my opinion, the use of associated Laguerre functions \bigl{[}L_{n+m}^{m}(x)\bigr{]}_{\text{BS}} defined by (4.5) is not recommendable. It follows from (4.6) that these functions cannot express generalized Laguerre polynomials $L_{n}^{(\alpha)}$ with nonintegral superscripts $\alpha$. This is both artificial and unnecessary. For example, the eigenfunctions $\Omega_{n,\ell}^{m}(\beta,\bm{r})$ of the Hamiltonian $\beta^{-2}\nabla^{2}-\beta^{2}r^{2}$ of the three-dimensional isotropic harmonic oscillator contain generalized Laguerre polynomials with half-integral superscripts (see for example [18, Eq. (5.4)] and references therein):

[TABLE]

These functions satisfy the orthonormality condition (see for example [18, Eq. (5.5)])

[TABLE]

and are complete and orthonormal in $L^{2}(\mathbb{R}^{3})$ (see also [86, pp. 36 - 37]).

In this article, only the mathematical notation for generalized Laguerre polynomials is used. Additional conventions used in physics were discussed by Kaijser and Smith [108, Footnote 1 on p. 48].

Guseinov and coworkers use the convention for Laguerre polynomials that is common in quantum mechanics (compare [107, Eq. (4) on p. 189] with [24, Eq. (3)]), although they give the book by Gradshteyn and Rhyzhik [109] as their primary reference on Laguerre polynomials. This is misleading: Gradshteyn and Rhyzhik use the mathematical definition of the generalized Laguerre polynomials (compare [109, Eq. (8.970.1) on p. 1061] with (4.2)).

In principle, every function set, that is complete and orthonormal in either $L^{2}(\mathbb{R}^{3})$ or also in a suitable weighted Hilbert space $L_{w}^{2}(\mathbb{R}^{3})$, can be used for the construction of one-range addition theorems. If we insist on an exponentially decaying function set that is complete and orthonormal in $L^{2}(\mathbb{R}^{3})$, then the following functions, which were introduced by Hylleraas [110, Footnote ∗ on p. 349], Shull and Löwdin [111], and Löwdin and Shull [112, Eq. (46)], are the most natural choice:

[TABLE]

Here, $\beta>0$ is a scaling parameter.

Lambda functions are orthonormal with respect to a integration over the whole $\mathbb{R}^{3}$:

[TABLE]

Accordingly, they are complete and orthonormal in $L^{2}(\mathbb{R}^{3})$.

Closely related to both bound-state hydrogen eigenfunctions and Lambda functions are the following functions which were already used in 1928 by Hylleraas [113, Eq. (25) on p. 478] and which are commonly called Coulomb Sturmians or simply Sturmians:

[TABLE]

If we replace in (4) $\beta$ by $Z/n$, we obtain the bound state eigenfunctions of a hydrogenlike ion with nuclear charge $Z$. These eigenfunctions are orthonormal with respect to an integration over the whole $\mathbb{R}^{3}$, but they are incomplete without the inclusion of the continuum eigenfunctions (see for example [17] and references therein).

Sturmians – or rather their three-dimensional Fourier transforms – occur also in Fock’s treatment of the hydrogen atom [114], albeit in a somewhat disguised form [18, Section VI].

It follows at once from (4.3) that the Sturmians satisfy [18, Eq. (4.7)]:

[TABLE]

Accordingly, the Sturmians are a complete orthogonal set in the weighted Hilbert space

[TABLE]

Unfortunately, the Hilbert space $L_{1/r}^{2}(\mathbb{R}^{3})$ is not necessarily suited for quantum mechanical calculations since we neither have $L^{2}(\mathbb{R}^{3})\subset L_{1/r}^{2}(\mathbb{R}^{3})$ nor $L_{1/r}^{2}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$. Thus, one might arrive at the premature conclusion that Sturmians are not particularly suited for our purposes.

If, however, we combine the differential equation satisfied by the Sturmians [18, Eq. (4.9)],

[TABLE]

with the orthogonality relationship (4.12), we see that the Sturmians satisfy the Sobolev-type orthonormality [18, Eq. (4.10)]

[TABLE]

Accordingly, Sturmians are complete and orthonormal in the Sobolev space $W_{2}^{(1)}(\mathbb{R}^{3}$.

In [24, Eq. (1)], Guseinov introduced a fairly large class of complete and orthonormal functions which, if the mathematical notation for generalized Laguerre polynomials is used, can be expressed as follows:

[TABLE]

The indices satisfy $n\in\mathbb{N}$, $k=-1,0,1,2,\dots$, $\ell\in\mathbb{N}_{0}\leq n-1$, $-\ell\leq m\leq\ell$, and the scaling parameter $\beta$ is positive. Guseinov and coworkers used these functions quite extensively in the context of one-range addition theorems [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 50, 47, 48, 49, 51, 53, 54, 55, 56, 57, 59, 60, 61, 62].

If Guseinov had used the mathematical definition (4.1) of the generalized Laguerre polynomials and not the unnecessarily restrictive convention (4.5), he probably would have noticed that the parameter $k=-1,0,1,2,\dots$ in his functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ need not be integral and can be generalized to real values $k\in[-1,\infty)$. One only has to replace the factorial $(n+\ell+k+1)!$ in (4) by the gamma function $\Gamma(n+\ell+k+2)$.

Guseinov’s functions satisfy the orthonormality relationship (compare also [25, Eq. (4)])

[TABLE]

Accordingly, Guseinov’s functions are a complete and orthonormal set in the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$ with $k=-1,0,1,2,\dots$, which is defined via the inner product

[TABLE]

and the norm

[TABLE]

according to

[TABLE]

For $k=-1$, the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ are apart from a slightly different normalization factor identical to the Sturmian functions defined by (4). Thus, they are complete and orthonormal in the Hilbert space $L_{1/r}^{2}(\mathbb{R}^{3})$ defined by (4) or – when used in combination with the inner product (2.12) with $\zeta=\beta$ – complete and orthogonal in the Sobolev space $W_{2}^{(1)}(\mathbb{R}^{3})$, which is a proper subspace of the Hilbert space $L^{2}(\mathbb{R}^{3})$. It is important to notice that $L_{1/r}^{2}(\mathbb{R}^{3})$ and $W_{2}^{(1)}(\mathbb{R}^{3})$ are not identical, although the Sturmians and thus Guseinov’s functions with $k=-1$ are complete and orthogonal in both spaces.

For $k=0$, the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ are identical to the Lambda functions defined by (4). Thus, they are complete and orthonormal in $L^{2}(\mathbb{R}^{3})$.

For $k=1,2,3,\dots$, the weighted Hilbert spaces $L_{r^{k}}^{2}(\mathbb{R}^{3})$ are genuinely different from the Hilbert space $L^{2}(\mathbb{R}^{3})$ and we neither have $L^{2}(\mathbb{R}^{3})\subset L_{r^{k}}^{2}(\mathbb{R}^{3})$ nor $L_{r^{k}}^{2}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$. The weight function $r^{k}$ in the inner product (4.18) admits functions to $L_{r^{k}}^{2}(\mathbb{R}^{3})$ that are too singular at the origin to belong to $L^{2}(\mathbb{R}^{3})$, and it excludes some functions belonging to $L^{2}(\mathbb{R}^{3})$ that decay slowly like a fixed power $r^{-\alpha}$ as $r\to\infty$.

If $f\in L_{r^{k}}^{2}(\mathbb{R}^{3})$, then it is guaranteed that the expansion

[TABLE]

in terms of Guseinov’s functions converges in the mean with respect to the norm (4.19) of $L_{r^{k}}^{2}(\mathbb{R}^{3})$. If, however, $f\in L^{2}(\mathbb{R}^{3})$, the convergence of this expansion with respect to the norm (4.19) of $L_{r^{k}}^{2}(\mathbb{R}^{3})$ is not guaranteed.

For fixed $k=-1,0,1,\dots$, the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ are by construction complete and orthogonal in the corresponding Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$, but for $k\neq 0$ the expansion (4.21) does not necessarily converge with respect to (4.19) for arbitrary functions $f$ belonging to the for our purposes most important Hilbert space $L(\mathbb{R}^{3})$ of square integrable functions. Thus, the convergence of (4.21) is only guaranteed if we have $f\in L^{2}(\mathbb{R}^{3})\cap L_{r^{k}}^{2}(\mathbb{R}^{3})$. This is a complications which we have to take into account whenever we do expansions in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k\neq 0$.

It is by no means clear whether and for which $k\neq 0$ a weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$ provides a proper setting for bound state electronic structure calculations. It certainly makes sense to be cautious and to avoid potential complications whenever possible. Therefore, I would refrain from using the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k\neq 0$ unless I see obvious computational benefits. Generalization for the sake of generalization is rarely a good idea.

In addition, the rate of convergence of the expansion of a given function in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ may depend quite strongly on the choice of $k$. For $k\geq 1$, the inner product (4.18) gives less weight to the region close to the origin and greater weight to the region away from the origin, which may or may not improve convergence.

The numerical consequences of a variation of $k$ can be studied by expanding the exponential

[TABLE]

in terms of Guseinov’s functions:

[TABLE]

The orthonormality of the spherical harmonics implies:

[TABLE]

This expression can be rewritten as follows:

[TABLE]

With the help of the integral [109, Eq. 7.414.8 on p. 850]

[TABLE]

we obtain

[TABLE]

and

[TABLE]

This result shows that the expansion (4.23) becomes trivial for $x=1$ because then only the term with $n=1$ is different from zero. However, it also shows that for $x\neq 1$ the rate of convergence of the expansion (4.23) decreases with increasing $k\in\mathbb{N}$.

This example does not imply that Guseinov’s functions necessarily produce inferior convergence rates for higher values of $k$. Most likely, other examples can be constructed for which higher values of $k$ improve convergence. However, this example shows that the inclusion of the nontrivial weight function $r^{k}$ in an inner product can have significant numerical consequences. Since it is not at all clear whether the weighted Hilbert spaces $L_{r^{k}}^{2}(\mathbb{R}^{3})$ with $k\neq 0$ provide a proper setting for bound state electronic structure calculations, I cannot recommend to routinely use Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k\neq 0$ for the construction of one-range addition theorems.

This assessment does not necessarily apply to Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k=-1$, which are up to a slightly different normalization factor identical to the Sturmians defined by (4). Sturmians are complete and orthogonal in the weighted Hilbert space $L_{1/r}^{2}(\mathbb{R}^{3})$, whose usefulness in electronic structure calculations is not clear since we neither have $L^{2}(\mathbb{R}^{3})\subset L_{1/r}^{2}(\mathbb{R}^{3})$ nor $L_{1/r}^{2}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$. But Sturmians are also complete and orthonormal in the Sobolev pace $W_{2}^{(1)}(\mathbb{R}^{3}$ defined by (2), which is a proper subspace of $L^{2}(\mathbb{R}^{3})$. As discussed in Section 2, a basis set has to be complete in $W_{2}^{(1)}(\mathbb{R}^{3}$ to guarantee the convergence of variational calculations. Completeness in $L^{2}(\mathbb{R}^{3})$ does not suffice.

Let me emphasize once more that Laguerre expansions converge in general only in the mean, but not necessarily pointwise (see for example [77]). Additional conditions, which a function has to satisfy in order to guarantee that its Laguerre expansion converges pointwise, were discussed by Szegö [115, Theorem 9.1.5 on p. 246] (see also [12, Appendix]).

5 Addition Theorems for Exponentially Decaying Functions

Pointwise convergent two-range addition theorems for all the commonly occurring exponentially decaying functions can be constructed via the two-range addition theorems [2, 116] of the so-called $B$ functions that were introduced by Filter and Steinborn [117, Eq. (2.14)]:

[TABLE]

Here, $\hat{k}_{n-1/2}$ is a reduced Bessel function of half-integral order $n-1/2$ defined by [118, Eqs. (3.1) and (3.2)]

[TABLE]

and $K_{\nu}(z)$ is a modified Bessel function of the second kind [102, p. 66]. If the order $\nu$ of a reduced Bessel function is half-integral, it can be expressed as an exponential multiplied by a terminating confluent hypergeometric series ${}_{1}F_{1}$ (see for example [119, Eq. (3.7)]).

$B$ functions are comparatively complicated mathematical objects, and it is not obvious why they should offer any advantages compared to other exponentially decaying function sets. However, $B$ functions possess a three-dimensional Fourier transform of remarkable simplicity:

[TABLE]

This is the most consequential and also the most often cited result of my PhD thesis [120, Eq. (7.1-6) on p. 160]. Later, (5) was published in [121, Eq. (3.7)]. Independently and almost simultaneously, (5) was also derived by Niukkanen [122, Eqs. (57) - (58)].

The exceptionally simple Fourier transform (5) gives $B$ functions a unique position among exponentially decaying functions. It also explains why other exponentially decaying functions like Slater-type functions with integral principal quantum numbers, bound state hydrogen eigenfunctions, and other functions based on generalized Laguerre polynomials such as Lambda functions, Sturmians, or Guseinov’s functions can all be expressed in terms of finite linear combinations of $B$ functions (details and further references can be found in [18, Section IV] or [2, Section 4]). Two-range addition theorems for the exponentially decaying functions mentioned above can be written down immediately by forming finite linear combinations of the corresponding addition theorems of $B$ functions [2, 116].

$B$ functions also turned out to be extremely useful for the derivation of one-range addition theorems. Filter and Steinborn [12, Eqs. (5.11) and (5.12)] derived symmetrical one-range addition theorems of the type of (3.4) for Lambda and $B$ functions by expanding them in terms of Lambda functions:

[TABLE]

These addition theorems contain Gaunt coefficients defined by (A). Filter and Steinborn [12, Eqs. (4.9), (4.10), and (4.33)] were able to derive explicit expressions for the coefficients $T_{n_{1}\ell_{1}}^{n_{2}\ell_{2}NL}$ and $a_{n_{1}\ell_{1}}^{n_{2}\ell_{2}N+LL}$ occurring in (5) and (5), respectively. In [18, Eq. (7.8)] it was shown that the explicit expression for the coefficients $T_{n_{1}\ell_{1}}^{n_{2}\ell_{2}NL}$ in the addition theorem (5) can also be derived via the weakly convergent expansion of the plane wave in terms of Lambda functions [18, Eq. (4.38)]. The same approach works also for the coefficients $a_{n_{1}\ell_{1}}^{n_{2}\ell_{2}N+LL}$ in the $B$ function addition theorem (5) (see [18, Eq. (7.7)]), and it should also work in the case of Slater-type functions with integral principal quantum numbers defined by (1.5), although this has apparently not been done yet.

The derivation of the addition theorems (5) and (5) by Filter and Steinborn was based on their remarkably compact convolution theorem of $B$ functions [123, Eq. (4.1)] which, however, can be derived more easily with the help of the Fourier transform (5) [121, Section V]:

[TABLE]

The abbreviation $\Delta\ell$ is defined by (A.9). For an application of this highly convenient expression for the construction of the addition theorem (5), Filter and Steinborn only had to express Lambda function in terms of $B$ functions. This can be accomplished with the help of the following expression [120, Eq. (3.3-35)] (see also [12, Eq. (3.17) and Ref. [23] on p. 2736]):

[TABLE]

This yields a finite sum representation of Lambda functions in terms of $B$ functions [12, Eq. (3.18)]:

[TABLE]

In this way, all overlap integrals of Lambda functions or between $B$ and Lambda functions, which occur according to (3.1) as intermediate steps in the derivation of the one-range addition theorems (5) and (5), can because of (5) be expressed in terms of $B$ functions.

Now, one only has to express all $B$ functions in terms of Lambda functions. This can be done with the help of the following relationship [124, Eq. (B.2) on p. 214]:

[TABLE]

As is well known, a generalized Laguerre polynomial with superscript $\beta$ can be expressed as a finite sum of generalized Laguerre polynomials with superscript $\alpha$ [102, p. 249]:

[TABLE]

This yields the following finite sum representation of a reduced Bessel function with half integral order in terms of generalized Laguerre polynomials with an essentially arbitrary superscripts $\alpha$:

[TABLE]

Filter and Steinborn [12, Eq. (3.13)] gave an expression which wrongly contains the additional factor $(-1)^{n}$.

Thus, we can express a $B$ function as a finite sum of Lambda functions (see also [12, Eq. (3.14)]):

[TABLE]

Starting from the Lambda function addition theorem (5), a symmetrical one-range addition theorem of the type of (3.4) for Slater-type functions with integral principal quantum numbers $n\in\mathbb{N}$ can be derived easily. As is well known, an integral power $x^{m}$ with $m\in\mathbb{N}_{0}$ can be expressed as a finite sum of generalized Laguerre polynomials [125, Eq. (2) on p. 207]:

[TABLE]

Thus, a Slater-type function with integral principal quantum number $n\in\mathbb{N}$ defined by (1.5) can be expressed as a finite sum of Lambda functions [12, Eq. (3.10)]:

[TABLE]

Accordingly, a symmetrical one-range addition theorem for Slater-type functions with integral principal quantum numbers $n\in\mathbb{N}$ in terms of Lambda functions can be written down immediately. The expansion coefficients of this addition theorem are just simple finite sums of the coefficients $T_{n_{1}\ell_{1}}^{n_{2}\ell_{2}NL}$ in (5).

Alternatively, one could just as well start from the $B$ function addition theorem (5) since a Slater-type function with an integral principal quantum number can also be expressed as a finite sum of $B$ functions [117, Eqs. (3.3) and (3.4)].

[TABLE]

If the principal quantum number $n$ is a positive integer, this expansion terminates because of the Pochhammer symbols $(-[n-\ell-1]/2)_{\sigma}$ and $(-[n-\ell]/2)_{\sigma}$ after a finite number of steps.

Thus, the expansion coefficients of the Slater addition theorem can also be expressed as simple finite sums of the coefficients $a_{n_{1}\ell_{1}}^{n_{2}\ell_{2}N+LL}$ in (5).

The finite sum formula (5.13) can be generalized to nonintegral powers $x^{\mu}$ with $\mu\in\mathbb{C}\setminus\mathbb{N}_{0}$. With the help of [109, Eq. (7.414.7) on p. 850]

[TABLE]

we obtain:

[TABLE]

If we set $\mu=m$ with $m\in\mathbb{N}_{0}$, the infinite series on the right-hand side terminates because of the Pochhammer symbol $(-m)_{n}$ and we obtain the finite sum (5.13). However, there is a fundamental difference between (5.13) and (5). The finite sum formula (5.13) is a relationship among polynomials. Therefore, it is valid pointwise for arbitrary $x\in\mathbb{C}$. In the case of the infinite series (5), we know that it converges in the mean in the radial Hilbert space (4), but we cannot assume that it converges pointwise. Moreover, the index dependence of the series coefficients on the right-hand side of (5) indicates that the convergence of this expansion is slow if it does not terminate. These convergence problems can be demonstrated by considering the special case $\mu=-1$ and $x\to 0$ in (5). Then, the left-hand side approaches $+\infty$. For an analysis of the behavior of the right-hand side of (5) as $x\to 0$, we use (4.1b) and obtain

[TABLE]

Inserting this into (4.1b) yields:

[TABLE]

Since $(1)_{n}=n!$, the series diverges for $x\to 0$ to $+\infty$, although each individual term of the series (5) with $\mu=-1$ is well behaved as $x\to 0$. Of course, pointwise convergence – or also divergence – of an expansion like (5) for nonintegral $\mu\in\mathbb{R}\setminus\mathbb{N}_{0}$ is not really important if we want to use it in integrals. However, the at best slow decay of the terms on the right-hand side of (5) indicates that expansions of a Slater-type function with nonintegral principal quantum number $n\notin\mathbb{N}$ in terms of Lambda functions or other exponentially decaying Laguerre-type function sets converge slowly.

By combining (5) with the Lambda function addition theorem (5), a symmetrical one-range addition theorem for a Slater-type functions with a nonintegral principal quantum numbers in terms of Lambda functions can be formulated easily. The only principal constraint on the nonintegral principal quantum number is that the corresponding Slater-type function must be square integrable, i.e., it must belong to the Hilbert space $L^{2}(\mathbb{R}^{3})$. However, the expansion coefficients of this one range addition theorem are now given by an infinite series involving the coefficients $T_{n_{1}\ell_{1}}^{n_{2}\ell_{2}NL}$ in (5). As remarked above, we have no a priori reason to assume that the convergence of this series would be rapid. So, unless somebody accomplishes a substantial simplification – for example by expressing an inner sum in closed form via a summation theorem of a generalized hypergeometric series with unit argument – we are confronted with a one-range addition theorem for Slater-type functions with nonintegral principal quantum numbers that is significantly more complicated and less suited for practical work than the corresponding addition theorem for Slater-type functions with integral principal quantum numbers.

Let us now assume that we have a symmetrical one-range addition theorem for Slater-type functions in terms of Lambda functions of the type of (5) or (5). Then we only have to replace the Lambda functions by Guseinov’s functions to obtain an expansion of a Slater-type function in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k=-1,0,1,2,\dots$ defined by (4). With the help of (5.10) we obtain the following expression for a Lambda function as a finite sum of Guseinov’s functions:

[TABLE]

In this way, expansions in terms of Lambda functions can be transformed easily to expansions in terms of Guseinov’s functions.

An inverse relationship – the expansion of Guseinov’s functions in terms of Lambda functions – can also be derived via (5.10):

[TABLE]

With the help of this finite sum, we could – starting from the addition theorem (5) for Lambda functions – immediately write down a symmetrical one-range addition theorem for Guseinov’s function. Moreover, all Lambda functions in this addition theorem could be replaced by Guseinov’s functions via (5).

The approach of Filter and Steinborn [12], which is based on the convolution theorem (5) of $B$ functions, can also be used to construct from the scratch symmetrical one-range addition theorems that are expansions in terms of Guseinov’s functions. Setting $w(r)=r^{k}$ and $\psi_{n,\ell}^{m}(\bm{r})=\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ in (3) yields the following one-range addition theorem:

[TABLE]

If $f\in L_{r^{k}}^{2}(\mathbb{R}^{3})$, this addition theorem converges in the mean according to the norm (4.18) of the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$.

For the derivation of an addition theorem of the type of (5.22) for Guseinov’s functions or $B$ functions via the convolution theorem (5) of $B$ functions, we have to express Guseinov’s functions in terms of $B$ functions. This can be done with the help of (5), yielding

[TABLE]

With the help of this relationship, two-range addition theorems of Guseinov’s function can be written down immediately by forming finite linear combinations of the corresponding addition theorems of $B$ functions [2, 116].

The weight function $r^{k}$ in (5.22) and (5.22), which is responsible for the orthonormality of Guseinov’s functions according to (4), can be absorbed with the help of [123, Eq. (6.1))]

[TABLE]

yielding

[TABLE]

If $s$ is an even integer, the Pochhammer symbol $(-s/2)_{\sigma}$ causes a truncation of the summation after a finite number of steps, and if $s$ is odd, this truncation is accomplished by the Pochhammer symbol $(-n-[s-1]/2)_{\sigma}$.

Finally, we have to express $B$ functions as finite linear combinations of Guseinov’s functions. With the help of (5), we obtain:

[TABLE]

With the help of (5), (5), and (5) it is possible to generalize the approach of Filter and Steinborn [12], which was based on the convolution theorem (5) of $B$ functions and which produced expansions in terms of Lambda functions, to one-range addition theorems for exponentially decaying functions that are expansions in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$.

We cannot expect that these addition theorems are as compact as the corresponding expansions in terms of Lambda functions, which can be obtained by setting $k=0$ in Guseinov’s functions. For example, (5) simplifies considerably for $k=-1$ and $k=0$ [18, Eqs. (4.19) and (4.20)]. Moreover, for $k=-1$ and $k=1,2,\dots$, the weight function $r^{k}$ has to be absorbed via (5), which produces an additional inner sum. So, unless some additional simplifications can be found, one-range addition theorems in terms of Guseinov’s functions will have a more complicated structure than the corresponding addition theorems in terms of Lambda functions.

This higher complexity is one reason why I am not interested in explicitly constructing symmetrical one-range addition theorems that are expansions in terms of Guseinov’s functions with $k\geq 1$. I doubt that these addition theorems would be useful enough to justify the effort. A second reason is that I have – as outlined in Section 4 – severe doubts whether the weighted Hilbert spaces $L_{r^{k}}^{2}(\mathbb{R}^{3})$ with $k\neq 0$ are really suited for bound state electronic structure calculations.

6 Guseinov’s One-Range Addition Theorems for Slater-Type

Functions

As discussed in Section 5, symmetrical one-range addition theorems of the type of (3.1) for exponentially decaying functions can be derived comparatively easily via the remarkably compact convolution theorem (5) of $B$ functions. However, Guseinov preferred to proceed differently.

In [23], Guseinov derived one-range addition theorems for Slater-type functions with integral principal quantum numbers by expanding them in terms of Sturmians and Lambda functions, and in [24] he did this for Slater-type functions with nonintegral principal quantum numbers. In [25], Guseinov introduced his functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ defined by (4) and used them for the construction of one-range addition theorems for Slater-type functions with integral principal quantum numbers. Since Sturmians and Lambda functions are special cases of Guseinov’s functions with $k=-1$ and $k=0$, respectively, these addition theorems generalize earlier expansions in terms of Lambda functions and Sturmians derived in [23]. In [26] Guseinov derived one-range addition theorems in terms of his functions for so-called “central and noncentral potentials” which are nothing but special Slater-type functions with integral and nonintegral principal quantum numbers. In [27], Guseinov provided a “unified treatment” of multicenter integrals of Slater-type functions with integral and nonintegral principal quantum numbers by expanding Slater-type functions in terms of his functions. In [28], Guseinov constructed addition theorems for his functions both in the coordinate as well as in the momentum representation by expanding his functions in terms of Slater-type functions with integral principal quantum numbers. In [29], Guseinov used his addition theorems for Slater-type functions to handle multicenter integrals of what he calls “central and noncentral interaction potentials” which are special Slater-type functions with integral and nonintegral principal quantum numbers. In [30], Guseinov provided again a “unified analytical treatment” of multicenter integrals of “central and noncentral interaction potentials” via one-range addition theorems of Slater-type functions derived with the help of his functions. In [31], Guseinov provided in this way a “unified analytical treatment” of multicenter nuclear attraction, electric field and electric field gradient integrals over Slater-type functions. In [32, 35], Guseinov provided another “unified treatment” of essentially the same integrals as in [31], but this type he emphasized the use of his functions. In [33], we find another “unified treatment” of multicenter integrals of Slater-type functions, but this time the potentials are called “integer and noninteger $u$ Yukawa-type screened Coulomb type potentials” (again, these potentials are nothing but special Slater-type functions with integral and nonintegral principal quantum numbers). In [34], Guseinov constructed one-range addition theorems for derivatives of Slater-type functions, and in [36, 37], he considered again multicenter integrals of “central and noncentral interaction potentials”. In [38], Guseinov provided another “unified treatment”, but this time of electronic attraction, electric field, and electric-field gradient multicenter integrals of “screened and nonscreened Coulomb potentials”. In [39], Guseinov derived one-range addition theorems for the Coulomb potential starting from his addition theorems for the Yukawa potential which is essentially a special Slater-type function $\chi_{N,L}^{M}$ with $N=L=M=0$. In [40, 50], Guseinov used the momentum space addition theorems derived in [28] for the construction of momentum space addition theorems for Slater-type functions. In [41, 42, 43, 45, 46], Guseinov considered one-range addition theorems for “Yukawa-like central and noncentral interaction potentials and their derivatives”, for derivatives of his functions, for “derivatives of integer and noninteger $u$ Coulomb-Yukawa type central and noncentral potentials”, and for “combined Coulomb and Yukawa like central and noncentral interaction potentials and their derivatives”, respectively. In [47], Guseinov considered one-range expansions for two-center charge densities of Slater-type functions with integral and nonintegral principal quantum numbers. Finally, in [48] Guseinov provided another “unified treatment” of expansion theorems and one-range addition theorems of complete orthonormal sets of functions in coordinate, momentum and four-dimensional spaces.

These examples show that Guseinov used quite a few different names for the functions he expanded, but with few exceptions he concentrated on one-range addition theorems for Slater-type functions with integral and nonintegral principal quantum numbers in the coordinate representation. Consequently. I will also focus on these addition theorems and will examine critically their derivation by Guseinov.

Guseinov exclusively derived unsymmetrical addition theorems of the type of (3.1) or (3.6), in which the two vectors $\bm{r}$ and $\bm{r}^{\prime}$ are treated differently. He never derived completely symmetrical addition theorems of the type of (3.4), (3), or (5.22).

Moreover, Guseinov expanded Slater-type functions in terms of functions with a different scaling parameter. Normally, this is a bad idea since overlap integrals with different scaling parameters are significantly more complicated than overlap integrals with equal scaling parameters (this is the reason why exclusively addition theorems with equal scaling parameters were discussed in Section 5). However, Guseinov needed this additional degree of freedom since he wanted to construct one-range addition theorems for the Coulomb potential by representing it as the limiting case of the Yukawa potential according to $1/r=\lim_{\beta\to 0}\exp(-\beta r)/r$.

Guseinov’s work on one-range addition theorems is based on expansions of Slater-type functions with integral or nonintegral principal quantum numbers in terms of his functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ with $k=-1,0,1,2,\dots$ and $\gamma>0$ (see for example [25, Eqs. (11) and (12)]):

[TABLE]

As remarked above, Guseinov did not derive completely symmetrical addition theorems of the type of (5.22) by expanding the overlap integrals (6.1) in terms of his functions. This is not completely satisfactory: In the Slater-type function $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$, the vectors $\bm{r}$ and $\bm{r}^{\prime}$ play an identical role, but not in the addition theorems (6.1).

Nevertheless, a definite assessment of the relative merits of symmetrical and unsymmetrical one-range addition theorems is not so easy. Since we cannot tacitly assume that one-range addition theorems necessarily lead to rapidly convergent expansions for multicenter integrals, the crucial question is whether it is easier to compute either overlap integrals like (6.1) or purely numerical coefficients as they occur in (3.3b), (3.7), and (5.22) both effectively and reliably for large indices.

It also makes a big difference if one has to integrate over only one of the two vectors $\bm{r}$ and $\bm{r}^{\prime}$ or over both of them. For example, in [39, Table I] Guseinov considered three-center one-electron nuclear attraction integrals

[TABLE]

where $f$ and $g$ are Slater-type functions, and the vectors $\bm{A}$, $\bm{B}$, and $\bm{C}$ are atomic centers. The integral (6) is the most complicated one-electron integral occurring in molecular electronic structure calculations based on the Hartree-Fock-Roothaan equations. By means of a shift of origin, one of the three vectors $\bm{A}$, $\bm{B}$, and $\bm{C}$ can be made to vanish. Thus, two addition theorems are needed to decouple the arguments in (6).

In the case of a one-electron integral like (6), the use of unsymmetrical one-range addition theorems of the type of (6.1) may well be feasible: It should not matter too much whether the remaining two of the three vectors $\bm{A}$, $\bm{B}$, and $\bm{C}$ occur in overlap integrals or in complete and orthonormal functions. Nevertheless, a comparison of Guseinov’s results [39, Table I] with those reported by Bouferguene and Jones [126], who had used two-range addition theorems, would have been of considerable interest.

The most difficult multicenter integrals occurring in the LCAO-MO approach are the notorious six-dimensional two-electron integrals

[TABLE]

which describe the Coulomb interaction of two in general nonclassical charge distributions $[f(\bm{r})]^{*}$ and $g(\bm{r}^{\prime})$ consisting of products of effective one-particle wave functions located at different atomic centers. We can only benefit in both integrations from the simplifying power of orthogonality if we use in (6.3) symmetrical one-range addition theorems of the type of (3.4), (3), or (5.22).

In principle, this applies also to multicenter integrals of those functions which Guseinov and coworkers had called “central and noncentral interaction potentials” [27, 29, 30, 33, 36, 37, 41, 42, 43, 45, 46, 56, 57, 61, 62] and which are nothing but special Slater-type functions with integral and nonintegral principal quantum numbers. In the case of one-electron integrals, unsymmetrical addition theorems of the type of (6.1) may well be sufficient, but if two-electron integrals of these interaction potentials should indeed be physically meaningful and have to be computed (see also [127]), I am skeptical about the feasibility of Guseinov’s approach based on the exclusive use of unsymmetrical addition theorems.

The central computational problem of Guseinov’s unsymmetrical one-range addition theorems (6.1) is the efficient and reliable evaluation of the overlap integrals (6.1) between Slater-type functions and the expansion functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ even for large indices.

Guseinov’s solution is simple but not necessarily good. It follows at once from the explicit expression (4.1) of the generalized Laguerre polynomials that Guseinov’s functions can be expressed as finite sums of Slater-type functions with integral principal quantum numbers:

[TABLE]

Accordingly, the overlap integrals (6.1) can be expressed as finite sums of overlap integrals of Slater-type functions:

[TABLE]

Using this in the addition theorems (6.1) yields:

[TABLE]

Superficially, this approach, which is mathematically completely legitimate, seems to have the advantage that existing programs for overlap integrals of Slater-type functions can be used for the evaluation of the overlap integrals (6.1) (this may be the reason for Guseinov’s approach). However, stability problems are likely. The explicit expression (4.1) of a generalized Laguerre polynomial $L_{n}^{(\alpha)}(x)$ in powers of $x$ becomes numerical unstable for larger values of the index $n$. Because of these stability problems, orthogonal polynomials are normally computed recursively.

So, as long as the opposite is not explicitly proved, it certainly makes sense to assume that an expression like (6.4) inherits the stability problems of (4.1). In addition, conventional programs for overlap integrals of Slater-type functions, as they are for instance used in semiempirical calculations, normally cannot be used in the case of (very) large principal and angular momentum quantum numbers. Thus, we need alternative expressions for the overlap integrals (6.1) not based on (6.4), which permit an efficient and reliable evaluation even for (very) large indices.

In his desire to reduce his whole formalism of one-range addition theorems to Slater-type functions, Guseinov even expressed the functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ on the right-hand side of (6) by Slater-type functions according to (6.4) (see for example [25, Eq. (14)]):

[TABLE]

This is still correct since the sum over $\nu^{\prime}$ is nothing but the function $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ in disguise. But Guseinov rearranged the order of summations in (6), obtaining expansions for Slater-type functions $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ with integral or nonintegral principal quantum numbers in terms of Slater-type functions $\chi_{n+\ell,\ell}^{m}(\beta,\bm{r})$ with integral principal quantum numbers located at a different center (see for example [25, Eq. (15)]):

[TABLE]

This step is potentially disastrous. A rearrangement of the order of summations of a double series is not always legitimate and can easily lead to a divergent result (see for example [128, Chapter V]). So, when going from (6) to (6), one cannot tacitly assume that the rearranged expansion (6) indeed converges.

As discussed in Section 2, a formal expansion of the type of (2.3) of a function $f$ in terms of a given function set $\{\varphi_{n}\}_{n=0}^{\infty}$ exists and converges in the mean if $f$ belongs to the underlying Hilbert space and if the expansion functions $\{\varphi_{n}\}_{n=0}^{\infty}$ are complete and orthonormal in this Hilbert space. If the expansion functions are only complete, but not orthogonal, it is possible to construct finite approximations of the type of (2.2) by minimizing the mean square deviation $\|f-f_{N}\|^{2}=(f-f_{N}|f-f_{N})$, but the existence of formal expansions of the type of (2.4) in terms of nonorthogonal functions is not guaranteed: These expansions may or may not exist.

Slater-type functions are complete in all the Hilbert space considered in this article (see for example [86, Section 4]), but not orthogonal. Thus, it is not clear whether Guseinov’s rearranged addition theorems (6) are mathematically meaningful. This has to be checked.

Convergence and existence problems due to a rearrangement can be illustrated via (5), which expresses a nonintegral power $x^{\mu}$ with $\mu\in\mathbb{C}\setminus\mathbb{N}_{0}$ as an infinite series of generalized Laguerre polynomials and which converges in the mean in the Hilbert space (4). If we insert the explicit expression (4.1) of the generalized Laguerre polynomials into (5) and interchange the order of summations, we obtain after some algebra:

[TABLE]

It looks as if we succeeded in constructing a power series expansion for the nonintegral power $x^{\mu}$. However, the generalized hypergeometric series ${}_{1}F_{0}$ with unit argument is the limiting case $y\to-1$ of the following binomial series satisfying [102, p. 38]

[TABLE]

We then obtain for the hypergeometric series in (6):

[TABLE]

Thus, the power series (6) does not exist and the Laguerre series (5) for $x^{\mu}$ with $\mu\in\mathbb{C}\setminus\mathbb{N}_{0}$ cannot be reformulated as a power series in $x$. Of course, this is not really surprising: The general power function $z^{\mu}$ with $z\in\mathbb{C}$ and $\mu\in\mathbb{C}\setminus\mathbb{N}_{0}$ is not analytic at $z=0$ in the sense of complex analysis. For $\mu=m$ with $m\in\mathbb{N}_{0}$, Taylor expansion of $z^{m}$ around $z=0$ yields the trivial identity $z^{m}=z^{m}$.

One-range additions theorems for exponentially decaying functions are fairly complicated mathematical objects. Consequently, explicit proofs of their convergence or divergence are very difficult. It is understandable that Guseinov was not interested in presenting such proofs, although he should have done so in order to justify his manipulations. Fortunately, valuable insight can in some cases be gained by considering not the addition theorems themselves, but their much simpler one-center limits.

Let us therefore assume that we succeeded in constructing either symmetrical addition theorems of the type of (5.22) or also unsymmetrical addition theorems involving overlap integrals for some function $f(\bm{r}\pm\bm{r}^{\prime})$ by expanding it in terms of Guseinov’s functions. We now consider the one-center limit by setting $\bm{r}^{\prime}=\bm{0}$. Then, our addition theorem must simplify to yield an identity for $f(\bm{r})$. Under fortunate circumstances, the mathematical nature of this identity can provide valuable insight.

First, we set $\beta=\gamma$ and $\bm{r}^{\prime}=\bm{0}$ in the addition theorems (6.1) for $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$. After the cancellation of all common factors we obtain expressions that are special cases of the expansion (5) expressing an in general nonintegral power $x^{\mu}$ as an infinite series of generalized Laguerre polynomials. As shown in (6) - (6), this series can only be reformulated as a power series in $x$ if $\mu$ is a nonnegative integer, $\mu=m$ with $m\in\mathbb{N}_{0}$, yielding the trivial identity $x^{m}=x^{m}$. If we have instead $\mu\in\mathbb{R}\setminus\mathbb{N}_{0}$, a rearranged power series in $x$ does not exist.

Thus, we can conclude that for $\beta=\gamma$, the one-center limit $\bm{r}^{\prime}=\bm{0}$ of the rearranged addition theorem (6) does not exist if the principal quantum number $N$ of the Slater-type function $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ is nonintegral, $N\in\mathbb{R}\setminus\mathbb{N}$.

Next, we assume $\beta\neq\gamma$ in the addition theorem (6.1) and set $\bm{r}^{\prime}=\bm{0}$. We then obtain after the cancellation of all common factors expressions that are special cases of the following expansion:

[TABLE]

This expansion, which can be derived with the help of (5), converges in the mean with respect to the norm of the weighted Hilbert space $L^{2}_{\mathrm{e}^{-x}x^{\alpha}}(\mathbb{R}_{+})$. For $u=0$, (6) simplifies to give (5). This can be shown easily with the help of Gauss’ summation theorem [102, p. 40].

If we insert the explicit expression (4.1) of the generalized Laguerre polynomials into (6) and interchange the order of summations, we also obtain a formal power series in $x$. Unfortunately, an analysis of the resulting power series becomes very difficult because of the terminating Gaussian hypergeometric series ${}_{2}F_{1}$ in (6) (probably, this would be a nontrivial research project in its own right). However, we can argue that the function $z^{\mu}\exp(uz)$ with $\mu,u,z\in\mathbb{C}$ is only analytic at $z=0$ in the sense of complex analysis if $\mu$ is a nonnegative integer, $\mu=m$ with $m\in\mathbb{N}_{0}$, yielding $z^{m}\exp(uz)=\sum_{n=0}^{\infty}u^{n}z^{m+n}/n!$. If $\mu$ is nonintegral, $\mu\in\mathbb{C}\setminus\mathbb{N}_{0}$, a power series around $z=0$ does not exist.

Thus, we can conclude that also for $\beta\neq\gamma$, the one-center limit $\bm{r}^{\prime}=\bm{0}$ of the rearranged addition theorems (6) for $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ does not exist if the principal quantum number $N$ is nonintegral, $N\in\mathbb{R}\setminus\mathbb{N}$.

From a mathematical point of view, a one-range addition theorem for a function $f(\bm{r}\pm\bm{r}^{\prime})$ is a mapping $\mathbb{R}^{3}\times\mathbb{R}^{3}\to\mathbb{C}$. We are interested in one-range addition theorems because we would like to have a unique (infinite series) representation of $f(\bm{r}\pm\bm{r}^{\prime})$ with separated variables $\bm{r}$ and $\bm{r}^{\prime}$ that is valid for the whole argument set $\mathbb{R}^{3}\times\mathbb{R}^{3}$. If we accept this premise, then we have to conclude that Guseinov’s manipulations, which produced the rearranged addition theorems (6) for $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$, are at least in the case of nonintegral principal quantum numbers $N\in\mathbb{R}\setminus\mathbb{N}$ a complete failure.

This observation does not rule out the possibility that an appropriate reinterpretation might make the rearranged addition theorems (6) with $N\in\mathbb{R}\setminus\mathbb{N}$ mathematically meaningful in a restricted sense as an approximation, although it does exist for the whole argument set $\mathbb{R}^{3}\times\mathbb{R}^{3}$. This has to be investigated. Moreover, it is still open whether the rearranged addition theorems (6) exists for the whole argument set $\mathbb{R}^{3}\times\mathbb{R}^{3}$ if the principal quantum number $N$ of the Slater-type function is a positive integer, $N\in\mathbb{N}$. Also this remains has to be investigated. In all cases, the burden of proof lies with Guseinov.

These observations should be particularly worrisome for those who advocate the use of Slater-type functions with nonintegral principal quantum numbers as basis functions in molecular electronic structure calculations. It is generally accepted that a good basis set should produce highly accurate approximations, but it should also lead to multicenter integrals that can be computed efficiently and accurately at tolerable computational costs. These two requirements have so far been mutually exclusive.

Slater-type functions with integral principal quantum numbers or other exponentially decaying function sets are well suited to produce highly accurate approximations in electronic structure calculations, but their multicenter integrals are notoriously difficult. In contrast, Gaussian functions are nonphysical. Consequently, large basis sets are needed to accomplish sufficiently accurate results. However, Gaussian functions have one, albeit decisive advantage: Their multicenter integrals can be computed comparatively easily.

If we nevertheless want to use physically better motivated exponentially decaying basis functions, we should concentrate on those functions that promise the simplest multicenter integrals, even if we should have to sacrifice some accuracy. Superficially, Slater-type functions with nonintegral principal quantum numbers look attractive since they promise somewhat better results than Slater-type functions with integral principal quantum numbers. But they achieve this at the cost of significantly more complicated multicenter integrals. Therefore, we would first need some fundamental mathematical breakthroughs to make Slater-type functions with nonintegral principal quantum numbers competitive with Slater-type functions with integral principal quantum numbers or with other exponentially decaying functions. As far as I can judge it, these breakthroughs are not in sight.

There is another, more general remark which I would like to make. Addition theorems of exponentially decaying functions as well as their multicenter integrals are fairly complicated mathematical objects. Accordingly, it is usually very difficult and often even practically impossible to justify nontrivial manipulations by rigorous mathematical proofs. This is highly embarrassing. The only way out, which I see, is to be very conservative and cautious in order to be on the safe side.

Let us assume that $f$ belongs to a suitable Hilbert space. If we expand $f(\bm{r}\pm\bm{r}^{\prime})$ in terms of a function set, which is complete and orthogonal in that Hilbert space, then the resulting one-range addition theorem converges in the mean with respect to the norm of the corresponding Hilbert space. Moreover, in many cases the Schwarz inequality (see for example [82, Eq. (6) on p. 31]) suffices to guarantee that the use of this addition theorem in a multicenter integral produces a convergent expansion.

However, we are no longer on the safe side if we replace the complete and orthogonal functions by other functions that are only complete but not orthogonal. Formally, this may yield another one-range addition theorem for $f(\bm{r}\pm\bm{r}^{\prime})$, but we do not know whether this expansion makes sense or not. This has to be proved. If we cannot prove this, then we should rather avoid such a possibly dangerous manipulation. It is not acceptable to ignore potential problems of that kind and hope for the best, as it was done by Guseinov.

7 Weakly Convergent One-Range Addition Theorems for the Coulomb

Potential

Guseinov [39] derived one-range addition theorems for the Coulomb potential via addition theorems for the Yukawa potential $\exp(-\beta r)/r$ [76], which can be viewed to be an exponentially screened Coulomb potential and which is also a special Slater-type function according to

[TABLE]

Guseinov’s idea was to set $N=L=M=0$ in his unsymmetrical one-range addition theorems for Slater-type functions discussed in Section 6 and to perform the limit $\beta\to 0$. Unfortunately, things are more complicated and Guseinov’s simplistic approach is fundamentally flawed.

For $k=0,1,2,\dots$, the Yukawa potential belongs to the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$ defined by (4), but not for $k=-1$ (see for example [19, p. 410]). This was apparently overlooked by Guseinov who used his functions with unspecified $k$. Accordingly, the addition theorems (6.1) with $N=L=M=0$ converge for $k=0,1,2,\dots$ in the mean with respect to the norm (4.19) of $L_{r^{k}}^{2}(\mathbb{R}^{3})$, but not for $k=-1$.

Guseinov replaced in his one-range addition theorems (6.1) the overlap integrals (6.1) involving his complete and orthonormal functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ by overlap integrals of Slater-type functions according to (6.5). This step, which yields the addition theorems (6), is mathematically legitimate, but numerically dubious. Also the next step – the replacement of the expansion functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ on the right-hand side of (6) by Slater-type functions according to (6.4), which yields the addition theorems (6) – is in the case of the Yukawa potential at least for $k=0,1,2\dots$ mathematically justified, but only as long as the order of the infinite summations in (6) is retained.

Guseinov’s final step – the rearrangement of the infinite summations in (6) which formally produces the addition theorems (6) with Slater-type functions as expansion functions – is not legitimate: The one-center limit of these addition theorems does not exist for $N=L=M=0$: If we set $\bm{r}^{\prime}=\bm{0}$ in the addition theorems (6.1) with $N=L=M=0$, we obtain expansions that are special cases of the Laguerre expansion (6) for $z^{\mu}\exp(uz)$ with $\mu=-1$. This Laguerre expansion cannot be reformulated as a power series in $z$ since $\exp(uz)/z$ is not analytic at $z=0$ in the sense of complex analysis. Accordingly, the rearranged addition theorem (6) with $N=L=M=0$ does not exist for the whole argument set $\mathbb{R}^{3}\times\mathbb{R}^{3}$.

Thus, Guseinov’s derivation [39] of one-range addition theorems for the Coulomb potential with Slater-type functions as expansion functions fails since his starting point – the unsymmetrical one-range addition theorems (6) with $N=L=M=0$ – do not exist for all arguments $\bm{r},\bm{r}^{\prime}\in\mathbb{R}^{3}$.

As a possible remedy, we could expand the Coulomb potential in terms of a suitable complete and orthonormal set of functions. If we expand $1/|\bm{r}-\bm{r}^{\prime}|$ in terms of Guseinov’s functions, we formally obtain the following symmetrical one-range addition theorems which are special cases of the general addition theorem (5.22):

[TABLE]

Analogous symmetrical addition theorems can also be derived for the Yukawa potential.

Instead of Guseinov’s functions, we could use in (7.2) or in analogous addition theorems for the Yukawa potential any other function set that is complete and orthonormal in an appropriate Hilbert space. An obvious, albeit Gaussian-type choice would be the eigenfunctions $\Omega_{n,\ell}^{m}(\beta,\bm{r})$ of the three-dimensional isotropic harmonic oscillator, which are defined by (4) and which are complete and orthonormal in $L^{2}(\mathbb{R}^{3})$.

In the case of (7.2) there is a principal problem that was apparently overlooked by Guseinov: The Coulomb potential does not belong to any of the Hilbert spaces considered in this article or implicitly used by Guseinov, since they all involve an integration over the whole $\mathbb{R}^{3}$. Thus, for $k=-1,0,1,2,\dots$ the addition theorems (7.2) diverge in the mean with respect to the norm (4.19) of $L_{r^{k}}^{2}(\mathbb{R}^{3})$.

It is tempting to conclude that all attempts of constructing one-range addition theorems for $1/|\bm{r}-\bm{r}^{\prime}|$ by expanding it in terms of function sets, that are complete and orthonormal with respect to an inner product involving an integration over the whole $\mathbb{R}^{3}$, are doomed. However, this conclusion is premature, and expansions of the type of (7.2) may well be our best chance of achieving our aim.

As discussed in Section 2, one-range addition theorems are based on Hilbert space theory and thus rely heavily on concepts from approximation theory. Loosely speaking, we may say that it is the purpose of approximation theory to provide convenient expressions – for example series expansions – that allow an efficient and reliable evaluation of a given mathematical object.

There is, however, a very important difference between the use of one-range addition theorems in multicenter integrals and more conventional applications of approximation theory: We are not interested in evaluating the functions that are represented by addition theorems. We only use these addition theorems to simplify the integrations in multicenter integral. Addition theorems are only intermediate results which ultimately produce series expansions for multicenter integrals.

The use of convergent expansions in integrals has the undeniable advantage that normally only comparatively mild assumptions are needed to guarantee that integration and summation can be interchanged and that the resulting expansions converge. Nevertheless, the use of convergent expansions in integrals is to some extent a luxury and not strictly necessary. We are free to use a divergent expansion in an integral and interchange integration and summation if we can guarantee that the resulting expansion converges to the correct result.

Obviously, such an approach gives us additional possibilities, but it would be naive to expect a free lunch: It is grossly negligent to use divergent series in integrals without knowing criteria of manageable complexity that guarantee the convergence of the resulting expansions.

As discussed in Section 1, two-range addition theorems of the type of the Laplace expansion (1) converge pointwise since they are essentially rearranged three-dimensional Taylor series in the Cartesian components of the shift vector. By relaxing the requirement of pointwise convergence, we arrive at one-range addition theorems which in general converge in the mean or strongly with respect to the norm of a suitable Hilbert space.

As a consequence, one-range addition theorems for functions, which do not belong to a given Hilbert space, diverge with respect to its corresponding norm. If we nevertheless insist on extending the formalism of one-range addition theorems to functions that do not belong to this Hilbert space, the requirement of convergence in the mean or strong convergence is too demanding. We have to replace it by a weaker type of converge that is related to and inspired by the theory of generalized functions in the sense of Schwartz [129].

Quantum theory is largely based on functionals, which are mappings from a function space to the complex numbers. Accordingly, the numerical values and the local properties of wave functions are – strictly speaking – irrelevant, and convergence of functionals, also called weak convergence, is all we really need. Stronger types of convergence such as pointwise convergence or convergence in the mean are in principle pure luxury.

Thus, the one-range addition theorems (7.2), which diverge for any $k=-1,0,1,2,\dots$ with respect to the norm (4.19) of the Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$, should be interpreted as weakly convergent expansions that produces correct results when used in suitable functionals.

In the context of multicenter integrals, weakly convergent expansions have largely been a terra incognita, but there is a detailed mathematical literature both on weak convergence in general (see for example [130, Chapters 10 and 11]) as well as on the expansion of generalized functions or distributions in terms of orthogonal polynomials (see for example [131, 132, 133] and references therein).

Weakly convergent expansions can be quite useful. In [18] weakly convergent expansions for $\exp(\mathrm{i}\bm{p}\cdot\bm{r})$ in terms of complete orthonormal and biorthogonal function sets were constructed. In some cases, these expansions simplify the evaluation of Fourier transforms, and they can also be used for the construction of one-range addition theorems (see [18, Section VII] or [19]).

Expansions for the plane wave, that closely resemble those derived in [18], were also constructed by Guseinov [28, Eqs. (45) - (46)]. Guseinov, who did not mention the article [18] in [28], either overlooked or ignored the obvious fact that that the plane wave does not belong to any of the Hilbert spaces which he implicitly used. Accordingly, Guseinov’s expansion diverge in the mean and can only converge weakly. Guseinov’s oversight is hard to understand because he had cited [18] in several other articles [25, 33, 34, 40, 41, 49, 52, 53, 57].

In my opinion, the essential features of weak convergence in contrast to strong convergence can be explained most easily via the Euclidean vector space $\mathbb{C}^{\infty}$ of infinite row or column vectors $\bm{u}=(u_{1},u_{2},\dots)$ with complex coefficients $u_{n}$. By equipping the vector space $\mathbb{C}^{\infty}$ with the inner product $(\bm{u}|\bm{v})=\sum_{n=1}^{\infty}[u_{n}]^{*}v_{n}$, we obtain the corresponding Hilbert space $\ell^{2}\subset\mathbb{C}^{\infty}$ with norm $\|\bm{u}\|=\sqrt{(\bm{u}|\bm{u})}$. The condition $\|\bm{u}\|<\infty$, which defines the Hilbert space $\ell^{2}$, can only be satisfied if the coefficients $u_{n}\in\mathbb{C}$ of $\bm{u}$ decay sufficiently fast as $n\to\infty$.

Let us now assume that some vector $\bm{w}\in\mathbb{C}^{\infty}$ cannot be normalized, i.e., $\|\bm{w}\|^{2}=(\bm{w}|\bm{w})=\sum_{n=1}^{\infty}|w_{n}|^{2}<\infty$ does not hold. Thus, $\bm{w}\notin\ell^{2}$, but this does not imply that all inner products $(\bm{w},\bm{u})$ with normalizable $\bm{u}\in\ell^{2}$ do not exist. If the coefficients $u_{n}$ of $\bm{u}$ decay sufficiently rapidly as $n\to\infty$, then the infinite series $(\bm{w}|\bm{u})=\sum_{n=1}^{\infty}[w_{n}]^{*}u_{n}$ may well converge to a finite result, although the series $(\bm{w}|\bm{w})=\sum_{n=1}^{\infty}|w_{n}|^{2}$ diverges.

It is important to note that $|(\bm{w}|\bm{u})|<\infty$ cannot hold for arbitrary $\bm{u}\in\ell^{2}$, but only for a suitably restricted subspace $\mathcal{W}\subset\ell^{2}\subset\mathbb{C}^{\infty}$ of normalizable vectors such that $\bm{\omega}\in\mathcal{W}$ implies $|(\bm{w}|\bm{\omega})|<\infty$. Accordingly, my interpretation of weak convergence resembles at least conceptually the theory of rigged Hilbert spaces or Gelfand triplets (see for example [134, 135] and references therein). A very readable account of rigged Hilbert spaces from the perspective of quantum mechanics and their relationship with Dirac’s bra and ket formalism can be found in the book by Ballentine [136, Chapter 1.4].

In principle, it is desirable to specify for a given $\bm{w}$ the whole set $\mathcal{W}$, but this may be difficult. In practice, it may well be sufficient to specify only a sufficiently large subset of $\mathcal{W}$. Let us for instance assume that the coefficients $w_{n}$ of a vector $\bm{w}\notin\ell^{2}$ are all finite and satisfy $|w_{n}|\sim n^{\beta}$ with $\beta\geq-1/2$ as $n\to\infty$. Thus, $\bm{w}$ cannot be normalized. If, however, the coefficients $\omega_{n}$ of a vector $\bm{\omega}\in\ell^{2}$ are all finite and satisfy $|\omega_{n}|\sim n^{\alpha}$ with $\alpha<-\beta-1$ as $n\to\infty$, then the infinite series $(\bm{w}|\bm{\omega})=\sum_{n=1}^{\infty}[w_{n}]^{*}\omega{}_{n}$ converges and the inner product $(\bm{w}|\bm{\omega})$ makes sense.

However, the condition $|\omega_{n}|\sim n^{\alpha}$ with $\alpha<-\beta-1$ as $n\to\infty$ does not suffice to specify the whole set $\mathcal{W}$ of vectors $\bm{\omega}$ with $|(\bm{w}|\bm{\omega})|<\infty$. It can happen that the series $\sum_{n=1}^{\infty}[w_{n}]^{*}\omega{}_{n}$ diverges but is summable. If summability techniques are included in our arsenal of mathematical techniques, we could even discard the otherwise essential requirement $\bm{\omega}\in\ell^{2}$ and try to make inner products $(\bm{w}|\bm{\omega})$ with both $\bm{w},\bm{\omega}\notin\ell^{2}$ mathematically meaningful. A (very) condensed review of the classic summability methods associated with the names of Cesàro, Abel, and Riesz can be found in Zayed’s book [137, Chapter 1.11.1]. More detailed treatments of linear summability methods can be found in specialized monographs such as the books by Boos [138], Hardy [139], Knopp [140], and Powell and Shah [141].

An interpretation of one-range addition theorems of the type of (7.2) as weakly convergent expansions requires that we first specify the functional, in which addition theorems are to be used. Then, regularity criteria have to be formulated whose validity guarantee the convergence of the resulting expansions for this functional.

The most difficult integrals, which occur in the LCAO-MO approach, are the notorious two-electron integrals $\mathcal{C}(f,g)$ defined by (6.3), whose evaluation can become extremely difficult if $f$ and $g$ are nonclassical two-center charge distributions of exponentially decaying functions.

A theoretical analysis of the integrals $\mathcal{C}(f,g)$ is by no means easy (see for example [142, Chapter 9]). In this context, it is instructive to replace in (6.3) the Coulomb potential by the Yukawa potential, yielding

[TABLE]

Obviously, we have $\mathcal{C}(f,g)=\lim_{\beta\to 0}\mathcal{Y}(f,g;\beta)$. Unfortunately, it is not guaranteed that this limiting process is continuous and produces a finite result. It can be shown with the help of Fourier transformation (see for example [143, Section II] and references therein) that $\mathcal{Y}(f,g;\beta)$ exists for arbitrary densities $f,g\in L^{2}(\mathbb{R}^{3})$, but this does not suffice to guarantee the existence of $\mathcal{C}(f,g)$. The densities $f$ and $g$ have to satisfy more sophisticated conditions than square integrability to guarantee this (see for example [142, Section 4.3 (Hardy-Littlewood-Sobolev inequality)] or [144, Example 3 (Sobolev’s inequality) on p. 31]).

If we want to apply weakly convergent addition theorems of the type of (7.2) for the evaluation of two-electron integrals, we need additional criteria that guarantee the convergence of the resulting expansions to the value of $\mathcal{C}(f,g)$ for physically and mathematically reasonable charge densities $f$ and $g$.

If we insert the addition theorems (7.2) into the integral (6.3), we obtain the following series expansions:

[TABLE]

If only a finite number of the integrals $\prescript{}{k}{\mathbf{F}}_{n,\ell}^{m}(\beta)$ and $\prescript{}{k}{\mathbf{G}}_{n^{\prime},\ell^{\prime}}^{m^{\prime}}(\beta)$ are nonzero, then the use of (7.2) in $\mathcal{C}(f,g)$ obviously produces the correct result.

In general, there will be infinitely many nonzero integrals $\prescript{}{k}{\mathbf{F}}_{n,\ell}^{m}(\beta)$ and $\prescript{}{k}{\mathbf{G}}_{n^{\prime},\ell^{\prime}}^{m^{\prime}}(\beta)$. In this case, asymptotic estimates for the coefficients $\prescript{}{k}{\mathbf{\Gamma}}_{n^{\prime}\ell^{\prime}m^{\prime}}^{n\ell m}(\beta)$ defined by (7.2) have to be constructed that hold in the case of large indices. Subsequently, asymptotic conditions on the decay of the integrals $\prescript{}{k}{\mathbf{F}}_{n,\ell}^{m}(\beta)$ and $\prescript{}{k}{\mathbf{G}}_{n^{\prime},\ell^{\prime}}^{m^{\prime}}(\beta)$, respectively, in the case of large indices can be formulated which guarantee the convergence of the series expansion (7.4).

It cannot be denied that such a proof of the convergence of the series expansions (7.4) – although most likely manageable – would be highly pedestrian. Nontrivial technical difficulties are also quite likely. Obviously, a more elegant convergence proof based on concepts of functional analysis and suitable functions spaces would be highly desirable. Unfortunately, this is currently not in sight.

There is also another problems. Kato [145] showed that the singularities of the potential of atomic and molecular Hamiltonians produce discontinuities of the wave functions commonly called Coulomb or correlation cusps (see also [146]). This implies that atomic or molecular wave functions do not possess continuous derivatives of arbitrary order at the locations of the nuclei. Consequently, exponentially decaying functions, which cannot be differentiated arbitrarily often at the origin, provide much better approximations than infinitely differentiable Gaussian-type functions. According to experience, a few exponentially decaying functions normally suffice to accurately model the discontinuities of atomic and molecular wave functions at the nuclei. But there is a price: Exponentially decaying functions are fairly complicated mathematical objects whose multicenter integrals are notoriously difficult.

It is quite likely that we also have to pay some price if we try to analyze the weak convergence of the addition theorems (7.2) in appropriate functionals such as the integral (6.3). As remarked above, weak convergence is related to Schwartz’s theory of generalized functions which applies to integrals containing a distribution multiplied by a so-called test function. To make these integrals well defined, the test functions must possess a sufficient amount of mathematical “niceness” in order to compensate the possibly highly irregular behavior of generalized functions such as the Dirac delta “function” or their derivatives.

In mathematics, it is common to use test functions with highly idealized properties because this greatly facilitates proofs. A frequently used test function space is the Schwartz space $\mathcal{S}(\mathbb{R}^{3})$ of rapidly decreasing and infinitely differentiable functions $\psi\colon\mathbb{R}^{3}\to\mathbb{C}$, which satisfy

[TABLE]

for all $k,m,n,u,v,w\in\mathbb{N}_{0}$ (see for example [81, p. 133]). The oscillator eigenfunctions $\Omega_{n,\ell}^{m}(\beta,\bm{r})$ defined by (4) or also other Gaussian-type functions belong to $\mathcal{S}(\mathbb{R}^{3})$, but none of the exponentially decaying function sets considered in this article.

Even more idealized are test function $\psi\colon\mathbb{R}^{3}\to\mathbb{C}$ belonging to the space $\mathcal{D}(\mathbb{R}^{3})$ consisting of infinitely differentiable functions with compact support (see for example [142, Chapter 6.2]).

Because of their infinite differentiability, test functions belonging to $\mathcal{S}(\mathbb{R}^{3})$ or $\mathcal{D}(\mathbb{R}^{3})$, respectively, are extremely convenient analytical tools. Their use greatly simply mathematical proofs. However, infinitely differentiable functions are not ideally suited for the representation of effective one-particle wave functions which are discontinuous and also decay exponentially (see for example [147, 148] and references therein).

The completeness of the oscillator eigenfunctions $\Omega_{n,\ell}^{m}(\beta,\bm{r})$ in $L^{2}(\mathbb{R}^{3})$ implies that they can approximate effective one-particle wave functions in the mean with an in principle unlimited precision, but one should not expect rapid convergence. If we instead use physically better motivated exponentially decaying functions, which are not infinitely differentiable, we can expect faster convergence, but all mathematical manipulations involving these functions will be (much) more difficult. Most likely, this also applies to the analysis of the properties of weakly convergent addition theorems of the type of (7.2).

Both the index $k$ as well as the scaling parameter $\beta$ of the expansion functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ influence the rate of convergence of the infinite series in (7.4). It is an obvious idea to try to maximize convergence by optimizing both $k$ and $\beta$. Unfortunately, our current level of understanding does not allow detailed predictions. Most likely, the “optimal” values of $k$ and $\beta$ depend quite strongly on the properties of the two charge densities $f$ and $g$. It is also possible that for given $f$ and $g$ only some values of $k$ will lead to a convergent series expansion (7.4).

The use of a nonzero $k$ can at least potentially lead to stability problems. For $k=0$, the integrals $\prescript{}{k}{\mathbf{F}}_{n,\ell}^{m}(\beta)$ and $\prescript{}{k}{\mathbf{G}}_{n^{\prime},\ell^{\prime}}^{m^{\prime}}(\beta)$ are the coefficients of expansions of $f$ and $g^{*}$ in terms of Lambda functions, which are according to (4.10) orthonormal with respect to an integration over the whole $\mathbb{R}^{3}$. Since the coefficients of orthogonal expansions satisfy Parseval’s equality (2.5), the integrals $\prescript{}{k}{\mathbf{F}}_{n,\ell}^{m}(\beta)$ and $\prescript{}{k}{\mathbf{G}}_{n^{\prime},\ell^{\prime}}^{m^{\prime}}(\beta)$ with $k=0$ are bounded in magnitude and vanish for large indices. For $k\neq 0$, Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ are not orthogonal with respect to an integration over the whole $\mathbb{R}^{3}$, but satisfy (4). Thus, $\prescript{}{k}{\mathbf{F}}_{n,\ell}^{m}(\beta)$ and $\prescript{}{k}{\mathbf{G}}_{n^{\prime},\ell^{\prime}}^{m^{\prime}}(\beta)$ with $k\neq 0$ are not the coefficients of orthogonal expansions and do not satisfy Parseval’s equality. Accordingly, it is not guaranteed that they are bounded in magnitude and that they vanish for large indices.

Let me emphasize that the convergence or divergence of a one-range addition theorem and the convergence or divergence of the resulting expansion for a multicenter integral are not directly related. In [39, pp. 212 - 213] Guseinov claimed that the convergence of his addition theorems for the Coulomb potential can be demonstrated via the convergence of the resulting series expansions for three-center integrals of the type of (6). This is of course wrong: Weak convergence does not imply convergence in the mean, let alone pointwise convergence (see for instance [149, §29]). Moreover, Guseinov should know that an observed agreement of different floating point computations up to a certain number of digits does not necessarily prove anything (see for example [150, Table 2]).

But even if the one-range addition theorems (7.2) or related addition theorems should lead to divergent series expansions for multicenter integrals, the resulting expansions can nevertheless be computationally useful. In the book by Bornemann, Laurie, Wagon, and Waldvogel [151, p. 225], there is the following instructive remark:

The question whether a series converges is largely irrelevant when the reason for using a series is to approximate its sum numerically.

Thus, even wildly divergent series can be used for computational purposes if suitable summation methods can be found.

I am aware of several predominantly theoretical articles dealing with the summation of divergent Laguerre expansions (see for example [152, 153] and references therein). However, from a practical point of view, it is probably more effective to use purely numerical summation methods. I am also skeptical about the classical summability methods mentioned above and discussed in books by Boos [138], Hardy [139], Knopp [140], and Powell and Shah [141]. Instead, I suggest to use nonlinear sequence transformations, which can also be applied in the case of slowly convergent sequences and series and which often achieve spectacular improvements of convergence. In that context, it may be of interest to note that the most recent (third) edition of the book Numerical Recipes [154] now also discusses nonlinear sequence transformations (for further details, see also [155]).

The best known and most often used sequence transformations are the so-called Padé approximants [156] which accomplish an acceleration of convergence or a summation by converting the partial sums of a power series to a doubly indexed sequence of rational functions. As documented by the long list of successful applications in the book by Baker and Graves-Morris [157], Padé approximants have become the standard tool in theoretical physics and in applied mathematics to overcome problems with slowly convergent or divergent power series.

It is, however, not so well known among non-specialists that alternative sequence transformations as for example Wynn’s epsilon and rho algorithm [158, 159], Brezinski’s theta algorithm [160], or Levin’s transformation [161] and later generalizations (see for example [162, 163, 164] and references therein) can at least for certain computational problems be much more effective than Padé approximants.

As emphasized already several times, we have no a priori reason to assume that one-range addition theorems necessarily lead to rapidly convergent expansions for multicenter integrals. Accordingly, our ability of evaluating infinite series effectively and reliably is crucial for the practical usefulness of addition theorems.

The conventional process of adding up the terms of a series successively is at least in principle able to produce approximations with unlimited accuracy for a convergent series, but it completely fails in the case of a divergent series. Moreover, adding up the terms successively is in far too many cases prohibitively inefficient. Therefore, it makes sense to try to use sequence transformations for the evaluation of series expansions for multicenter integrals whenever possible. This is not a new idea. The oldest article using sequence transformations for the evaluation of multicenter integrals, which I am aware of, was published in 1967 by Petersson and McKoy [165].

Personally, I became interested in sequence transformations during my PhD thesis [120], in which series expansions for multicenter integrals played a major role. Since it (too) often happened that my series expansions converged slowly (see for example [119, Table II]), it was a natural idea to speed up convergence with the help of sequence transformations. During my PhD thesis, I only knew linear series transformations as described in the classic, but now outdated book by Knopp [140], which turned out to be ineffective. I was completely ignorant of the more powerful nonlinear transformations, which often accomplish spectacular improvements of convergence. Only later, I used nonlinear sequence transformations with considerable success for the evaluation of multicenter integrals of exponentially decaying functions [166, 167, 168, 169, 170].

The usefulness of sequence transformations in the context of multicenter integrals is not limited to the evaluation of infinite series representations. Particularly noteworthy seems to be Safouhi’s approach. Starting from the Fourier transform (5) of $B$ functions, Safouhi converts complicated multicenter integrals of $B$ or Slater-type functions to multidimensional integral representations that have to be evaluated by numerical quadrature.

At first sight, this does not look like a good idea because the oscillatory nature of the multidimensional integral representations makes the straightforward application of conventional quadrature methods difficult. However, these problems can be overcome by combining quadrature schemes with suitable nonlinear sequence transformations. Based on previous work of Sidi [171] and of Levin and Sidi [172], Safouhi succeeded in developing some extrapolation techniques specially suited to his needs. This permits a remarkably efficient and reliable evaluation of complicated multicenter integrals via oscillatory integral representations (see for example [173, 174, 175, 176, 177, 178, 179] and references therein).

In my opinion, Safouhi’s work is a convincing demonstration of the practical usefulness of extrapolation and convergence acceleration techniques in electronic structure calculations. Safouhi’s work also shows that there is no reason for despair if existing sequence transformations turn out to be not powerful enough to solve certain problems. It may well be possible to construct new transformations that can do the job.

Sequence transformations can be useful also in completely different contexts. For example, I have applied sequence transformations successfully in such diverse fields as the evaluation of special functions and related objects [163, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189], the summation of strongly divergent quantum mechanical perturbation expansions [186, 190, 191, 192, 193, 194, 195, 182, 196, 184, 197, 198, 199, 200, 201, 202, 203], the prediction of unknown series coefficients [190, 194, 195, 201, 204], and the extrapolation of quantum chemical crystal orbital and cluster electronic structure calculations for oligomers to their infinite chain limits of stereoregular quasi-onedimensional organic polymers [205, 206, 207]. Many other applications of sequence transformations are listed in my Habilitation thesis [184] and in [164].

These examples should convince even a skeptical reader that sequence transformations are extremely useful computational tools, and that it is worth while to invest time and effort to understand their power as well as their shortcomings. Of course, sequence transformations are no panacea, and not all my attempts were successful. However, even my failures often turned out to be fruitful in the long run, since they frequently provided new insight which ultimately paved the way for the construction of new convergence acceleration and summation techniques [163, 164, 167, 183, 186, 187, 196, 208, 209].

In my opinion, it is hard to understand that there are still researchers like Guseinov who try to evaluate series expansions for multicenter integrals without the help of sequence transformations. Multicenter integrals of exponentially decaying functions are notoriously complicated, and a nice explicit expression for a multicenter integral does not necessarily permit its efficient and reliable evaluation. It is also necessary to employ powerful and sophisticated numerical techniques. The conventional process of adding up the terms of a series successively until convergence is achieved does not fall into this category.

The currently most complete reference on Padé approximants is the impressive book by Baker and Graves-Morris [157]. The most recent monograph on sequence transformations is the book by Sidi [210]. It contains a wealth of information and is undeniably useful for specialists, but presentation and choice of topics makes it basically unsuited for novices (compare also my book reviews [211, 212]). So, the best choice may well be the book by Brezinski and Redivo Zaglia [213]. Some people also like my long review [163]. As a first introduction for novices, I can recommend Appendix A of the book by Bornemann, Laurie, Wagon, and Waldvogel [151].

8 Differentiation Techniques and the Spherical Tensor Gradient

Operator

Often, it comparatively easy to obtain explicit analytical expressions for multicenter integrals over scalar functions, which are irreducible spherical tensors of rank zero with respect to their local (atomic) coordinate systems. If, however, the functions occurring in the multicenter integral are irreducible spherical tensors of higher ranks, the derivation of explicit expressions can easily become extremely difficult or even practically impossible.

It is also an empirical fact that it is usually much easier to differentiate than to integrate (integration is an art, but every fool and in particular also computer algebra systems like Maple or Mathematica can differentiate). Accordingly, it is an obvious idea to try to generate an explicit expression for a difficult multicenter integral over nonscalar functions by differentiating the simpler multicenter integral over scalar functions – preferably the simplest scalar functions – with respect to scaling parameters and/or nuclear coordinates (see for example [121, Section IV]). This approach is not restricted to the notoriously difficult multicenter integrals of exponentially decaying functions. It is also used in the case of multicenter integrals of the so-called spherical Gaussian functions (see for example [214, 215, 216, 217, 218] and references therein).

It is relatively easy to generate multicenter integrals of higher scalar functions by differentiating the simplest scalar functions with respect to their scaling parameters. In the case of a $1s$ Slater-type or Gaussian function, we can construct higher scalar functions easily by repeatedly using the relationships $\partial\exp(-\alpha r)/\partial\alpha=-r\exp(-\alpha r)$ or $\partial\exp(-\alpha r^{2})/\partial\alpha=-r^{2}\exp(-\alpha r^{2})$.

The generation of anisotropic functions, which are irreducible spherical tensors of rank $\ell$, from scalar functions is less straightforward, but also here an highly developed mathematical technology based on the differential operator $\mathcal{Y}_{\ell}^{m}(\nabla)$ is available (see for example [3] and [4, Section 3] and references therein).

In [39, Section 3], Guseinov used differentiation techniques for the derivation of more complicated addition theorems by differentiating his one-range addition theorems for $1/|\bm{r}-\bm{r}^{\prime}|$ with respect to the Cartesian components of either $\bm{r}=(x,y,z)$ or $\bm{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$. Guseinov pursued essentially analogous approaches in his articles [31, 32, 33, 34, 35, 38, 41, 42, 43, 45, 46, 56, 57, 60, 61].

It is not a new idea to generate more complicated addition theorems of anisotropic functions by differentiating simpler addition theorems of scalar functions. In [219, Sections IV and VI], it was shown that the two-range addition theorems of the irregular solid harmonic and the modified Helmholtz harmonic, respectively, can be derived by differentiating the two-range addition theorems for the Coulomb potential and the Yukawa potential, respectively. Similar ideas were pursued in [116].

Unfortunately, this does not imply that Guseinov’s approach is computationally efficient or at least mathematically sound. This Section first describes how irreducible spherical tensors can be differentiated comparatively easily with respect to the Cartesian components of their argument vectors. The second topic is less technical and refers to more fundamental mathematical problems: In [31, 32, 33, 34, 35, 38, 39, 41, 42, 43, 45, 46, 56, 57, 60, 61], Guseinov consistently ignored all questions of convergence and existence of his addition theorems obtained by differentiating simpler addition theorems.

This article considers exclusively addition theorems of irreducible spherical tensors of the type of (1.1). Moreover, all functions occurring in these addition theorems are also irreducible spherical tensors of a given rank $\ell\geq 0$. Irreducible spherical tensors $F_{\ell}^{m}(\bm{r}\pm\bm{r}^{\prime})$, that are of interest in the context of electronic structure calculations, can all be differentiated at least a finite number of times with respect to the Cartesian components of either $\bm{r}$ or $\bm{r}^{\prime}$. Because of the convenient orthonormality properties of the spherical harmonics, it is highly desirable to express the angular part of Cartesian differentiations in terms of spherical harmonics. The necessary algebra can be done, but in particular for large angular momentum quantum numbers of the irreducible spherical tensors and for large orders of the differential operators, we would be confronted with messy expressions and nontrivial technical problems. There is also the danger that delta function contributions, which can occur in multicenter integrals and addition theorems (see for example the articles by Pitzer Kern, and Liscomb [220] and by Kay, Todd, and Silverstone [221] or the book by Judd [222, Chapter 5.3]), are easily overlooked if we differentiate with respect to Cartesian components (see also [223, 224, 225, 226]).

A much more convenient approach is possible that completely avoids all differentiations of irreducible spherical tensors $F_{\ell}^{m}(\bm{r}\pm\bm{r}^{\prime})$ with respect to Cartesian components of the argument vectors and only requires differentiations with respect to the radial variables $r=|\bm{r}|$ and $r^{\prime}=|\bm{r}^{\prime}|$, respectively. This is accomplished by reformulating the original differential operators with unspecified transformation properties as finite sums of differential operators that are irreducible spherical tensors of a given integral rank.

Let us assume that $\mathcal{P}_{n}(\bm{r})$ is a polynomial of degree $n$ in the Cartesian components of $\bm{r}=(x,y,z)$:

[TABLE]

The regular solid harmonic ${\mathcal{Y}}_{\ell}^{m}(\bm{r})$ is according to (A) a homogeneous polynomial of degree $\ell$ in the Cartesian components of $\bm{r}=(x,y,z)$. Thus, the completeness and the orthonormality of the surface spherical harmonic $Y_{\ell}^{m}(\theta,\phi)$ with respect to an integration over the unit sphere in $\mathbb{R}^{3}$ implies that the polynomial $\mathcal{P}_{n}(\bm{r})$ can be expressed as a finite sum of solid harmonics in $\bm{r}$ multiplied by even powers of $r=|\bm{r}|$ (see also [227, §96 on pp. 147 - 148]):

[TABLE]

This is a relationship among polynomials in the Cartesian components of an essentially arbitrary three-dimensional vector $\bm{r}$. Thus, (8.2) also holds if we replace $\bm{r}$ by $\nabla=(\partial/\partial x,\partial/\partial y,\partial/\partial z)$, and any differential operator, which is polynomial of degree $n$ in the Cartesian components of $\nabla$, can be expressed as a finite sum of products of integral powers of the Laplacian $\nabla^{2}=\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}+\partial^{2}/\partial z^{2}$ multiplied by a solid harmonic in $\nabla$:

[TABLE]

We obtain an explicit expression for the spherical tensor gradient operator $\mathcal{Y}_{\ell}^{m}(\nabla)$ by replacing in (A) the Cartesian components of $\bm{r}$ by those of $\nabla$:

[TABLE]

In Martin’s book [4, p. 62], $\mathcal{Y}_{\ell}^{m}(\nabla)$ is called the Erdélyi operator. I prefer instead the name spherical tensor gradient operator because the expressions for products $\mathcal{Y}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})$, that will be mentioned later, generalize – as emphasized by Bayman [228, p. 2558] – the well known gradient formula in angular momentum theory (see for example [229, Chapters 5.7 and 5.9] or [230, Chapter II.11]).

In Martin’s book [4, Section 3], the mathematical theory of $\mathcal{Y}_{\ell}^{m}(\nabla)$ as well as numerous applications predominantly in classical physics are discussed. In [3], one also finds a detailed discussion of the mathematical properties of $\mathcal{Y}_{\ell}^{m}(\nabla)$ as well as numerous applications, but this time emphasis is on electronic structure theory.

A new differential operator is not necessarily a useful thing, let alone a major achievement. As remarked above, differentiating an irreducible spherical tensor $F_{\ell}^{m}(\bm{r})$ of the type of (1.1) with respect to the Cartesian components of $\bm{r}$ produces messy expressions. So, if we look at (8) and take into account (8), we may get the impression that we replaced something complicated – the left-hand side of (8) – by something even more complicated.

We thus arrive at the paradoxical conclusion that the differential operators $\mathcal{P}_{n}(\nabla)$ and ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ are practically useful only if it is not necessary to differentiate an irreducible spherical tensor $F_{\ell}^{m}(\bm{r})$ with respect to the Cartesian components of $\bm{r}$ via the defining explicit expression (8) and (8). Fortunately, this is the case. Differentiations with respect to Cartesian components can be avoided completely. The Laplacian $\nabla^{2}$ is an irreducible spherical tensor of rank zero, and ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ is just like the corresponding regular solid harmonic ${\mathcal{Y}}_{\ell}^{m}(\bm{r})$ an irreducible spherical tensor of rank $\ell$ (a formal proof can be found in [6, p. 312]). Consequently, we can hope for substantial computational and technical benefits if products involving ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ and other irreducible spherical tensors are handled via the powerful machinery of angular momentum coupling.

The first article on the differential operator ${\mathcal{Y}}_{\ell}^{m}(\nabla)$, which I am aware of, is due to Hobson who derived in 1892 a very consequential theorem on the differentiation of functions $f\colon\mathbb{R}^{n}\to\mathbb{C}$ [231, p. 67]. This theorem is also discussed in Hobson’s book [227, pp. 124 - 129], in Martin’s book [4, Section 3.5], and in [3, Section 3].

With the help of Hobson’s theorem, it is comparatively easy to obtain explicit expressions for the product of $\mathcal{Y}_{\ell}^{m}(\nabla)$ and a radially symmetric function $\varphi\colon\mathbb{R}^{3}\to\mathbb{}C$ that only depends on $r=|\bm{r}|$ (see for example [3, Eq. (3.4)]):

[TABLE]

Thus, ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ is a generating differential operator that transforms an irreducible spherical tensor of rank zero to an irreducible spherical tensor of rank $\ell$.

With the help of Hobson’s theorem it is also possible to construct explicit expressions for products ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})$. For that purpose, let us assume that the irreducible spherical tensor $F_{\ell_{2}}^{m_{2}}(\bm{r})$ satisfies a relation of the kind of (8.5), i.e., it can be generated by applying ${\mathcal{Y}}_{\ell_{2}}^{m_{2}}(\nabla)$ to a suitable scalar function $\Phi_{\ell_{2}}(r)$:

[TABLE]

In view of ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})={\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla){\mathcal{Y}}_{\ell_{2}}^{m_{2}}(\nabla)\,\Phi_{\ell_{2}}(r)$, we need an explicit expression of manageable complexity for the product ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla){\mathcal{Y}}_{\ell_{2}}^{m_{2}}(\nabla)$. This can be accomplished easily. Since Gaunt coefficients defined by (A) linearize the product of two spherical harmonics, multiplication of (A) by $r^{\ell_{1}+\ell_{2}}$ yields the linearization formula for the product of the regular solid harmonics $\mathcal{Y}_{\ell_{1}}^{m_{1}}(\bm{r})$ and ${Y}_{\ell_{2}}^{m_{2}}(\bm{r})$ (see for example [3, Eq. (3.6)]). In this linearization formula, we only have to replace the Cartesian components of $\bm{r}$ by those of $\nabla$ and obtain (see for example [3, Eq. (3.7)]):

[TABLE]

The abbreviation $\Delta\ell$, which is either a positive integer or zero, is defined by (A.9).

By combining (8.5), (8.6), and (8) we obtain (see for example [3, Eq. (3.8)]):

[TABLE]

In principle, (8) should suffice for our purposes since the scalar function $\Phi_{\ell_{2}}(r)$ in (8.6) can according to (8.5) be obtained from the scalar function $f_{\ell_{2}}(r)$ in (1.1) by repeated integration with respect to $r$. However, repeated integrations can at least potentially lead to nontrivial technical problems. Alternative expressions for the product ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})$ are thus desirable.

By systematically exploiting the tensorial nature of the spherical tensor gradient operator, the product ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})$ can be expressed as a finite linear combination of Gaunt coefficients, radial functions $\gamma_{\ell_{1}\ell_{2}}^{\ell}(r)$, and spherical harmonics [232, Eq. (4.7)]:

[TABLE]

The functions $\gamma_{\ell_{1}\ell_{2}}^{\ell}(r)$ in (8) can be obtained by differentiating the radial part $f_{\ell_{2}}(r)$ of the spherical tensor $F_{\ell_{2}}^{m_{2}}(\bm{r})$ with respect to $r=|\bm{r}|$ (see [232, Eqs. (3.29), (4.15)

• (4.18), and (4.24)] or [3, Eqs. (4.11) - (4.16)]):

[TABLE]

The abbreviations $\Delta l$, $\Delta l_{1}$, $\Delta l_{2}$, and $\sigma(\ell)$ are defined by (A.9) - (A.12).

Other expressions for the product ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})$ can be found in articles by Bayman [228], Santos [233], Stuart [234], Niukkanen [235], and Rashid [236].

There are some radially symmetric functions of considerable relevance in electronic structure theory that lead to remarkably simple expressions if ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ is applied to them via (8.5). The classic example is the Coulomb potential $1/r$. Hobson [231] showed that the irregular solid harmonic $\mathcal{Z}_{\ell}^{m}(\bm{r})=r^{-\ell-1}Y_{\ell}^{m}(\theta,\phi)$ is generated by applying ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ to $1/r$ (further details can be found in Hobson’s book [227, pp. 124 - 129]). In modern notation, Hobson’s result can be expressed as follows (see for example [121, Eq. (4.16)]):

[TABLE]

Thus, (8) and (8.16) imply that Cartesian derivatives of the Coulomb potential as considered by Guseinov in [39, Section 3] or by Guseinov and Mamedov [60] can be expressed by linear combinations of irregular solid harmonics $\mathcal{Z}_{\ell}^{m}$ multiplied by integral powers of the Laplacian $\nabla^{2}$. If a nonzero power of the Laplacian acts on an irregular solid harmonic, we obtain (see for example [237, Eq. (29)]):

[TABLE]

Here, $\delta_{\ell}^{m}$ is the so-called spherical delta function (see for example [237, Eq. (30)]):

[TABLE]

If we set in (8.17) $\ell=0$, we obtain the well known Poisson equation of a unit point charge:

[TABLE]

Thus, the spherical delta function $\delta_{\ell}^{m}$ can be viewed as a generalized solution of the Poisson equation of a unit multipole charge.

In [32, 33, 34, 35, 38, 41, 42, 43, 45, 46, 56, 57, 61], Guseinov considered Cartesian derivatives of the Yukawa potential and more generally of Slater-type functions. Also in these cases, is is advantageous to rewrite differentiations with respect to Cartesian components in tensorial form according to (8).

In this context, it is highly advantageous that Slater-type functions with integral principal quantum numbers can according to (5) be expressed as finite linear combinations of $B$ functions. The reason is that it is remarkably easy to apply the spherical tensor gradient operator to $B$ functions. For example, application of ${\mathcal{Y}}_{\ell}^{m}(\nabla)$ to a scalar $B$ function simply produces a nonscalar $B$ function [121, Eq. (4.12)]:

[TABLE]

It is also fairly easy to apply the spherical tensor gradient operator to a nonscalar $B$ function [232, Eq. (6.25)]:

[TABLE]

If we set $n=\ell=m=0$ in (5), we see that the Yukawa potential is also a special $B$ function:

[TABLE]

Thus, (8.20) implies that the application of the spherical tensor gradient operator to the Yukawa potential yields the so-called modified Helmholtz harmonic [219, Eq. (6.9)]:

[TABLE]

In the limit of vanishing screening ($\beta\to 0$), the modified Helmholtz harmonic approaches an irregular solid harmonic according to [143, Eq. (3.10)]

[TABLE]

Moreover, the modified Helmholtz harmonic satisfies the following exponentially screened variant of the anisotropic Poisson equation (8.17) [3, Eq. (6.18)]:

[TABLE]

In view of this relationship and also because of [232, Eq. (5.6)]

[TABLE]

which shows that the differential operator $1-\beta^{-2}\nabla^{2}$ of the modified Helmholtz equation acts as a ladder operator in the case of $B$ functions, it makes sense to define a distributional $B$ function as the following derivative of the three-dimensional Dirac delta function [232, Eq. (6.20)]:

[TABLE]

It is also easy to apply an integral power of the Laplacian to a $B$ function. The binomial expansion of $\beta^{-2\nu}\,\nabla^{2\nu}$ in powers of $1-\beta^{-2}\nabla^{2}$ in combination with (8.26) yields [232, Eq. (5.7)]:

[TABLE]

The $B$ function relationships given above show that the application of the spherical tensor gradient operator multiplied by integral powers of the Laplacian to the Yukawa potential leads to remarkably compact expressions. Since (8.28) also holds for distributional $B$ functions of the type of (8), it is almost trivially simple to keep track of delta function contributions.

If only first and second order derivatives with respect to the Cartesian components of the Coulomb or the Yukawa potential are needed, it is not really necessary to do the differentiations via (8), and it is also not too difficult to keep track of delta function contributions. If, however, differentiations of (very) high orders have to be done, the tensorial version (8) offers substantial advantages.

For example, the two-range addition theorems discussed in [1, 2] and in [3, Section 7] were derived via the following decomposition of the translation operator in (1.3) in terms of tensorial invariants:

[TABLE]

Apparently, this expansion was first published by Santos [233, Eq. (A.6)], who emphasized that it should be useful for the derivation of addition theorems, although he never used it for that purpose. However, the addition theorems derived in [1, 2] and in [3, Section 7] show that the tensorial decomposition (8) is indeed a practically useful mathematical tool, although it involves irreducible spherical tensors of arbitrary order.

After this essentially technical digression on Guseinov’s way of differentiating functions $f\colon\mathbb{R}^{3}\to\mathbb{C}$ with respect to the Cartesian components of the argument vector, more fundamental questions such as the existence and convergence of derivatives of one-range addition theorems will be discussed.

As emphasized already several times, two-range addition theorems are essentially rearranged three-dimensional Taylor expansions in the Cartesian components of the shift vector, which converge pointwise and in suitable open sets even uniformly. Accordingly, two-range addition theorems represent ordinary functions, which are locally smooth and can be differentiated. Thus, it is possible to generate more complicated two-range addition theorems by differentiating simpler ones (see for example [116, 219]).

In contrast, Guseinov’s one-range addition theorems discussed in Sections 6 and 7 – if they exist at all – either converge in the mean with respect to the norm of an appropriate Hilbert space or they converge weakly, i.e., they only make sense in suitable functionals restricted to subsets of certain Hilbert spaces. As is well known, smoothness is neither necessary for strong or weak convergence nor implied by it. Thus, it is not clear whether Guseinov’s one-range addition theorems represent sufficiently smooth functions that can be differentiated, and if yes, in which sense they can be differentiated.

These problems are not uncommon in mathematics. Often, it is comparatively easy to prove convergence in the mean, but very difficult or even impossible to prove pointwise convergence or even uniform convergence directly. Scenarios like this one obviously raise the question, under which conditions uniform convergence is implied by convergence in the mean. A condensed review of these issues and a discussion of some elementary results can be found in a recent article by Ford and Pennline [238]. Related topics, albeit with an emphasis on distribution theory, are discussed in the book by Strichartz [239, Chapter 8].

Guseinov’s derivatives of strongly convergent one-range addition theorems for Slater-type functions can be made mathematically meaningful if it can be shown that these addition theorems do not only converge in the mean, but also pointwise and uniformly.

Probably, there is some hope in the case of the one-range addition theorems (6.1), which are expansions in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ and which converge with respect to the norm (4.19) of the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$, if the principal quantum number $N$ of the Slater-type function $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ is a positive integer $N\in\mathbb{N}$. We would have to check whether the additional conditions discussed in Szegö’s book [115, Theorem 9.1.5 on p. 246] are satisfied in this case. Such an approach may also be successful in the case of the derivatives of the one-range addition theorems for Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ considered in [42, 43].

I am, however, very skeptical about all addition theorems for Slater-type functions $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ with nonintegral principal quantum number $N\in\mathbb{R}\setminus\mathbb{N}$.

The situation is more complicated if weakly convergent one-range addition theorems of the type of (7.2) are to be differentiated. Weakly convergent addition theorems are essentially expansions of generalized functions in the sense of Schwartz [129], whose derivatives also have to be interpreted in this sense. If $f$ is a generalized function and $\psi\colon\mathbb{R}^{3}\to\mathbb{C}$ is a suitable test function, derivatives of $f$ are defined by the following functional [129, Chapter II]:

[TABLE]

If $\psi\in\mathcal{S}(\mathbb{R}^{3})$ or $\psi\in\mathcal{D}(\mathbb{R}^{3})$, $\psi$ is infinitely differentiable and it is possible to define weak derivatives of arbitrary order of a generalized function $f$ via (8).

There is a detailed mathematical literature on derivatives of generalized functions and their interpretation and application. For example, there is a monograph by Ziemer on weakly differentiable functions [240]. Unfortunately, the mathematical literature is not particularly helpful if we want to apply derivatives of one-range addition theorems in multicenter integrals and if we insist on using exponentially decaying functions. Mathematicians almost exclusively use infinitely differentiable test functions belonging to either $\mathcal{S}(\mathbb{R}^{3})$ or $\mathcal{D}(\mathbb{R}^{3})$. However, Kato’s work on cusps [145] shows that atomic and molecular wave functions are not infinitely differentiable. This is also true if we approximate atomic and molecular wave functions in variational calculations by exponentially decaying functions.

If we want to use one-range addition theorems, which are derivatives of generalized functions, in multicenter integrals, we are confronted with the annoying problem that the remainder of the integrand, which assumes the role of the test function, not only has to decay sufficiently rapidly for large arguments, but also has to possess continuous derivatives of a sufficiently high order. Obviously, this complicates considerably our attempts to prove that the resulting expansions converge to the correct results. Mathematicians know why they prefer infinitely differentiable test functions belonging to either $\mathcal{S}(\mathbb{R}^{3})$ or $\mathcal{D}(\mathbb{R}^{3})$.

To illuminate these problems, let us apply the spherical tensor gradient operator $\mathcal{Y}_{L}^{M}(\nabla)=\mathcal{Y}_{L}^{M}(\nabla_{\bm{r}})$ to the one-range addition theorems (7.2) for the Coulomb potential. Then, we formally obtain with the help of (8.16) one-range addition theorems for the irregular solid harmonic:

[TABLE]

An essentially identical expression – the factor $(-1)^{L}$ would be missing – can be derived by differentiating $1/|\bm{r}-\bm{r}^{\prime}|$ instead with respect to $\bm{r}^{\prime}$.

If we use the addition theorems (8) in a multicenter integral, we essentially have two options: Either, we can try to construct an explicit expression for the product $\mathcal{Y}_{L}^{M}(\nabla)\,\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$, or we could use (8) and apply $\mathcal{Y}_{L}^{M}(\nabla)$ to the $\bm{r}$-dependent part of the remainder of the integrand, which assumes the role of the test function.

If we apply $\mathcal{Y}_{L}^{M}(\nabla)$ to $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ and use one of the numerous expressions for $\gamma_{\ell_{1}\ell_{2}}^{\ell}(r)$ defined by (8), it should be possible to construct a closed form expression for this product. However, the resulting functions of $\bm{r}$ will no longer be orthogonal and they may also have a (much) more complicated structure.

Alternatively, we can use (8) and shift $\mathcal{Y}_{L}^{M}(\nabla)$ to the remainder to the $\bm{r}$-dependent part of the integrand. If this happens to be a nonclassical two-center density, we have to apply the Leibniz theorem of the spherical tensor gradient operator (see [3, Eq. (4.22)] and references therein). Messy expressions are then likely, possibly involving singular or distributional contributions. In either case, nontrivial technical difficulties are to be expected. One must not forget that we still have to prove that the resulting series expansion for the multicenter integral converges to the correct result.

If we should indeed need weakly convergent one-range addition theorems of the type of (7.2) for the irregular solid harmonic $\mathcal{Z}_{L}^{M}(\bm{r}-\bm{r}^{\prime})$, then it is probably easier to construct it from the scratch by expanding $\mathcal{Z}_{L}^{M}$ in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$. I strongly suspect that in this case it would be much easier to prove convergence of the resulting expansion to the correct result.

One the basis of our current level of knowledge (or rather the lack of it), the use of derivatives of weakly convergent addition theorems in multicenter integrals would be purely experimental.

9 Summary and Conclusions

The efficient and reliable evaluation of multicenter integrals is among the oldest mathematical and computational problems of molecular electronic structure theory (a review of the older literature can be found in an article by Dalgarno [241]). In spite of the efforts of numerous researchers including Guseinov and coworkers, no completely satisfactory solution has been found yet. The situation is particularly unsatisfactory in the case of the notoriously difficult integrals of exponentially decaying functions, but even in the case of the much simpler integrals of Gaussians, there is still active research going on (see for example [10, 214, 215, 216, 217, 218, 242] and references therein).

Different centers in the integrand make the evaluation of an integral difficult. They prevent the straightforward separation of a multicenter into products of simpler integrals. As discussed in Section 1, addition theorems, which are expansions of a given function $f(\bm{r}\pm\bm{r}^{\prime})$ with $\bm{r},\bm{r}^{\prime}\in\mathbb{R}^{3}$ in products of other functions depending on either $\bm{r}$ or $\bm{r}^{\prime}$, are principal tools that can accomplish a separation of variables, albeit at the cost of infinite series expansions.

In electronic structure calculations, predominantly those addition theorems have played a major role that depend on $\bm{r}$ and $\bm{r}^{\prime}$ only indirectly via $\bm{r}_{<}$ and $\bm{r}_{>}$ and thus possess a two-range form. The prototype of such an addition theorem is the Laplace expansion (1) of the Coulomb potential. Unfortunately, the use of two-range addition theorems in multicenter integrals can easily lead to nontrivial technical problems. This explains why many authors have tried to construct one-range addition theorems that can be applied more easily.

Two-range addition theorems can be constructed by applying the translation operator $\mathrm{e}^{\bm{r}_{<}\cdot\nabla_{>}}$ in its tensorially invariant form (8) to $f(\bm{r}_{>})$, which is by assumption analytic at $\bm{r}_{>}$. Accordingly, two-range addition theorems are essentially rearranged three-dimensional Taylor expansions in the Cartesian components of the shift vector $\bm{r}_{<}$ that converge pointwise and in suitable open sets even uniformly. If $f$ is not analytic at the origin, its addition theorem converges only if it has a two-range form.

Classical and complex analysis is mainly concerned with power series that converge pointwise. However, pointwise convergence is a very demanding requirement, and its scope is too limited to cover all cases of interest. If we want to avoid the troublesome two-range form of addition theorems, we have to replace pointwise convergence by a weaker form of convergence.

One-range addition theorems can be constructed by exploiting Hilbert space theory, whose basic facts are reviewed in Section 2. In this approach, a function $f(\bm{r}\pm\bm{r}^{\prime})$ belonging to a suitable Hilbert space is according to (3.1) expanded in terms of a function set $\{\varphi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ that is complete and orthonormal in this Hilbert space. In general, such an expansion converges only in the mean, i.e., with respect to the norm of the underlying Hilbert space, but not pointwise.

The natural Hilbert space for electronic structure calculations based on effective one-particle wave functions is the Hilbert space $L^{2}(\mathbb{R}^{3})$ of square integrable functions defined by (2). As discussed in Section 3, it is also possible to construct one-range addition theorems that converge with respect to the norm of a suitable weighted Hilbert space $L_{w}^{2}(\mathbb{R}^{3})$ defined by (2). Here, $w(\bm{r})\geq 0$ is a suitable weight function. There is, however, a principle problem: In general, we neither have $L^{2}(\mathbb{R}^{3})\subset L_{w}^{2}(\mathbb{R}^{3})$ nor $L_{w}^{2}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$. Accordingly, it is by no means easy to find a nontrivial weight function $w(\bm{r})\neq 1$ that leads to substantial computational benefits and simultaneously avoids major drawbacks.

If we demand that the complete and orthonormal functions $\{\varphi_{n,\ell}^{m}(\bm{r})\}_{n,\ell,m}$ are also irreducible spherical tensors of the type of (1.1), we more or less automatically arrive at functions based on the generalized Laguerre polynomials. Section 4 discusses Lambda functions $\Lambda_{n,\ell}^{m}(\beta,\bm{r})$ defined by (4), Sturmians $\Psi_{n,\ell}^{m}(\beta,\bm{r})$ defined by (4), and Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$ with $k=-1,0,1,2,\dots$ defined by (4).

Guseinov’s functions satisfy the orthonormality relationship (4) involving the weight function $w(\bm{r})=r^{k}$. Accordingly, they are complete and orthonormal in the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$ with $k=-1,0,1,2,\dots$ defined by (4). Guseinov’s functions provide some unification since they contain for $k=0$ and $k=-1$, respectively, Lambda functions and Sturmians as special cases. However, this unification does not apply to approximation processes. For $k\neq 0$, we neither have $L_{r^{k}}^{2}(\mathbb{R}^{3})\subset L^{2}(\mathbb{R}^{3})$ nor $L^{2}(\mathbb{R}^{3})\subset L_{r^{k}}^{2}(\mathbb{R}^{3})$. Thus, $f\in L^{2}(\mathbb{R}^{3})$ does not imply that an expansion of $f$ in terms of Guseinov’s functions with $k\neq 0$ necessarily converges in $L_{r^{k}}^{2}(\mathbb{R}^{3})$.

One-range addition theorems for exponentially decaying functions had been considered by several researchers before Guseinov. The so far most compact results were obtained by Filter and Steinborn [12] who derived the symmetrical one-range addition theorems (5) and (5) for Lambda and $B$ functions, respectively, by expanding them in terms of Lambda functions, which are complete and orthonormal in $L^{2}(\mathbb{R}^{3})$. In Section 5, it is shown that the approach of Filter and Steinborn, which was based on the convolution theorem (5) of $B$ functions, works also in case of Slater-type functions. The approach of Filter and Steinborn can be generalized to expansions in terms of Guseinov’s functions. However, these addition theorems have a more complicated structure than the corresponding expansions in terms of Lambda functions.

Section 6 describes Guseinov’s mathematically dubious attempts of deriving one-range addition theorems for Slater-type functions. Guseinov first derived the expansions (6.1) of a Slater-type function $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ with integral or nonintegral principal quantum number $N$ in terms of his functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\gamma,\bm{r})$ with $\gamma>0$. Guseinov’s results are mathematically correct, although not necessarily optimal. In a later step, Guseinov replaced his orthonormal functions by nonorthogonal Slater-type functions with integral principal quantum numbers according to (6.4) and rearranged the order of summations in his addition theorems. In this way, Guseinov obtained the expansions (6) of $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$ in terms of Slater-type functions $\chi_{n,\ell}^{m}(\gamma,\bm{r})$ with integral principal quantum numbers $n\in\mathbb{N}$.

As reviewed in Section 2, a function $f$ belonging to a given Hilbert space can be expanded in terms of complete and orthogonal functions. Such an expansion converges with respect to the norm of the Hilbert space, and the coefficients of this expansion satisfy Parseval’s equality (2.5). Thus, the expansion coefficients are bounded in magnitude and they also vanish for large indices.

If the orthogonal functions are replaced by nonorthogonal functions, it is possible to construct finite approximations of the type of (2.2). Unfortunately, it is not guaranteed that an infinite expansion of the type of (2.4) exists, and if it does, it can happen that the coefficients of this expansion are unbounded and do not vanish for large indices.

Accordingly, Guseinov’s manipulations, which produced the rearranged addition theorem (6) from (6.1), are dangerous. Thus, he should have proved that his manipulations are legitimate and lead to a meaningful result.

Addition theorems of exponentially decaying functions are fairly complicated mathematical objects. Therefore, it is very difficult or practically even impossible to prove explicitly that a given addition theorem converges or diverges. This applies also to Guseinov’s rearranged addition theorem (6) for Slater-type functions $\chi_{N,L}^{M}(\beta,\bm{r}\pm\bm{r}^{\prime})$. However, as shown in Section 6, the one-center limit $\bm{r}^{\prime}=\bm{0}$ of Guseinov’s addition theorem (6) does not exist if the principal quantum $N$ is not a positive integer. Thus, for $N\in\mathbb{R}\setminus\mathbb{N}$ the rearranged addition theorem (6) does not exist for the whole argument set $\mathbb{R}^{3}\times\mathbb{R}^{3}$. The remaining question, whether the rearranged addition theorem (6) exists for $N\in\mathbb{N}$, cannot be decided in this way and is still open.

The nonexisting rearranged addition theorems (6) with $N=L=M=0$ were the starting point for Guseinov’s construction [39] of one-range addition theorems for the Coulomb potential, but this is not the only problem. The Coulomb potential does not belong to any of the Hilbert spaces implicitly used by Guseinov, since they all involve an integration over the whole $\mathbb{R}^{3}$. Accordingly, expansions of $1/|\bm{r}-\bm{r}^{\prime}|$ in terms of Guseinov’s functions $\prescript{}{k}{\Psi}_{n,\ell}^{m}(\beta,\bm{r})$, which formally yield the symmetrical one-range addition theorems (7.2), diverge for any $k=-1,0,1,2,\dots$ in the mean with respect to the norm of the weighted Hilbert space $L_{r^{k}}^{2}(\mathbb{R}^{3})$.

This observation seems to imply that all attempts of constructing a symmetrical one-range addition theorem for $1/|\bm{r}-\bm{r}^{\prime}|$ by expanding it in terms of functions, that are complete and orthonormal with respect to an inner product involving an integration over the whole $\mathbb{R}^{3}$, are futile. Although seemingly obvious, this conclusion is premature and the situation is better than it may look.

It is possible to use divergent one-range addition theorems like (7.2) in multicenter integrals in a mathematically rigorous way. However, convergence in the mean, which in the mathematical literature is frequently called strong convergence, is too demanding for that. We have to replace it by an even weaker form of convergence. Instead, we require that a possibly divergent one-range addition theorem produces meaningful results when used in a multicenter integral. Obviously, this approach is inspired by the theory of generalized functions in the sense of Schwartz [129].

Weak convergence of addition theorems is all we really need in variational electronic structure calculations. Pointwise convergence or convergence in the mean – although undeniably convenient – are not really necessary. Thus, the use of weakly convergent addition theorems offers new perspectives and simplifies or even solves certain technical problems. However, it would be overly optimistic to expect that we could get all these benefits without having to pay a price.

The key problem, which makes the application of weakly convergent addition theorems difficult, is the formulation of sufficiently simple criteria which guarantee the convergence of the resulting expansions for multicenter integrals. In that respect, a lot of work remains to be done.

Regularity conditions – or rather their absence – also play a major role in Section 8. In several articles, Guseinov had constructed more complicated one-range addition theorems by differentiating simpler one-range addition theorems with respect to the Cartesian components of their argument vectors. Guseinov’s one-range addition theorems – if they exist at all – either converge in the mean with respect to the norm of an appropriate weighted Hilbert space or they converge weakly, i.e., they only make sense as generalized functions in suitable functionals. It is thus not at all clear whether Guseinov’s addition represent functions in the ordinary sense that are locally smooth and can be differentiated. We need to know much more before we can safely apply derivatives of one-range addition theorems in multicenter integrals.

Section 8 also discusses a more convenient alternative to the troublesome differentiation of irreducible spherical tensors of the type of (1.1) with respect to the Cartesian components of their argument vectors. The alternative proposed here is based on the reformulation (8) of a differential operator as a finite sum of integral powers of the Laplacian multiplied by irreducible spherical tensor gradient operators $\mathcal{Y}_{\ell}^{m}(\nabla)$ defined by (8). This differential operator plays a decisive role in the derivation of two-range addition theorems as rearranged Taylor expansions [1, 2, 3].

At first sight, the use of (8) looks like a bad idea. However, $\mathcal{Y}_{\ell}^{m}(\nabla)$ is an irreducible spherical tensor of rank $\ell$. Consequently. products ${\mathcal{Y}}_{\ell_{1}}^{m_{1}}(\nabla)F_{\ell_{2}}^{m_{2}}(\bm{r})$ can be simplified considerably with the help of the powerful machinery of angular momentum coupling, yielding the finite linear combination (8) of Gaunt coefficients, radial functions, and spherical harmonics. Thus, no differentiations with respect to the Cartesian components of $\bm{r}$ are needed. We only have to differentiate with respect to the radial variable $r=|\bm{r}|$.

It is the intention of this article to demonstrate that the analytical tools of classical and complex analysis do not suffice for a successful treatment of one-range addition theorems. More modern and also more sophisticated concepts such as Hilbert spaces and generalized functions are indispensable.

In particular, we need different concepts of convergence. If a function $f$ belongs to a suitable Hilbert space, one-range addition theorems can be constructed that converge in the mean with respect to the norm of that Hilbert space, even if $f$ is not analytic everywhere in $\mathbb{R}^{3}$. If, however, $f$ does mot belong to that Hilbert space, one-range addition theorems can only converge weakly in the sense of generalized functions.

Analogous considerations apply also with respect to numerical techniques. In the context of multicenter integrals, addition theorems are essentially mathematical recipes that generate infinite series expansions. We have no a priori reason to assume that these series expansions converge rapidly. Thus, our assessment of the usefulness of addition theorems depends crucially on our ability of evaluating slowly convergent or even divergent series effectively and reliably.

The conventional process of adding up the terms of an infinite series successively until convergence is finally achieved is in far too many cases hopelessly inefficient. It is thus an obvious idea to use convergence acceleration and summation techniques. Of course, we do not know whether the currently known convergence acceleration and summation techniques, which were mentioned shortly in Section 7, are powerful enough to produce satisfactory results or whether it will be necessary to construct new and hopefully more powerful transformations. This has to be investigated.

Let me emphasize: If we want to evaluate even difficult multicenter integrals effectively and reliably via series expansion, then there is no alternative to the systematic and intelligent use of convergence acceleration and summation techniques.

It is dangerous to rely too much on completeness. The completeness of a function set in a Hilbert space only implies that any $f$ belonging to this Hilbert space can be expanded in terms of these functions and that this expansion converges in the mean. Unfortunately, convergence is a very weak statement. In particular, the convergence of an infinite series does not imply that this series is numerically useful. A simple and yet striking example is provided by the Dirichlet series $\zeta(s)=\sum_{n=0}^{\infty}(n+1)^{-s}$ for the Riemann zeta function. This series converges if $\mathrm{Re}(s)>1$, but is notorious for extremely slow convergence if $\mathrm{Re}(s)$ is only slightly larger than 1. For example, in [207, p. 194] it was estimated that in the order of $10^{600}$ terms of the Dirichlet series would be needed to compute $\zeta(1.01)$ with an accuracy of 6 decimal digits.

Appendix A Spherical Harmonics and Gaunt Coefficients

If we choose the phase convention of Condon and Shortley [106, Chapter III.4], the spherical harmonic can be expressed as follows [6, p. 69]:

[TABLE]

Here, $P_{\ell}^{|m|}(\cos\theta)$ is an associated Legendre polynomial [107, p. 155]. Alternative phase conventions for the spherical harmonics are discussed in [243, pp. 17 - 22].

The spherical harmonics $Y_{\ell}^{m}(\theta,\phi)$ are often called surface harmonics because the angles $\theta$ and $\phi$ characterize a point $\bm{r}/r$ on the surface of the three-dimensional unit sphere. In the literature, it is common to introduce the so-called regular and irregular solid harmonics

[TABLE]

The regular solid harmonic ${\mathcal{Y}}_{\ell}^{m}(\bm{r})$ is a homogeneous harmonic polynomial of degree $\ell$ in the Cartesian components of $\bm{r}=(x,y,z)$ [6, p. 71]:

[TABLE]

Moreover, $\mathcal{Y}_{\ell}^{m}(\bm{r})$ is for all $\bm{r}\in\mathbb{R}^{3}$ a solution of the homogeneous three-dimensional Laplace equation

[TABLE]

whereas the irregular solid harmonics are generalized solutions of the Poisson equation (8.17) of a unit multipole charge, yielding the spherical delta function (8.18).

The so-called Gaunt coefficient [244] is the integral of the product of three spherical harmonics over the surface of the unit sphere in ${\mathbb{R}}^{3}$:

[TABLE]

It follows from the orthonormality of the spherical harmonics that Gaunt coefficients linearize the product of two spherical harmonics:

[TABLE]

The symbol $\sum\!{}^{(2)}$ indicates that the summation proceeds in steps of two. The summation limits in (A) are given by [245, Eq. (3.1)]

[TABLE]

A compact review of the properties of Gaunt coefficients and additional references can be found in [3, Appendix C].

In this article, the following abbreviations are used:

[TABLE]

If the three orbital angular momentum quantum numbers $\ell_{1}$, $\ell_{2}$, and $\ell$ satisfy the summation limits (A.8), then these quantities are either positive integers or zero.