A series of the form =0c(κ,)Mκ,+1/2(r0)Wκ,+1/2(r)P(cos(γ)) is evaluated explicitly where c(κ, ) are suitable complex coefficients, Mκ,μ and Wκ,μ are the Whittaker functions, P are the Legendre polynomials, r0 < r are radial variables, γ is an angle and κ is a complex parameter. The sum depends, as far as the radial variables and the angle are concerned, on their combinations r + r0 and (r2+r022rr0cos(γ))1/2. This addition formula generalizes in some respect Gegenbauer’s Addition Theorem and follows rather straightforwardly from some already known results, particularly from Hostler’s formula for Coulomb Green’s function. In addition, several complementary summation formulas are derived. They suggest that a further extension of this addition formula may be possible.

Green’s function of a Hamiltonian is an important object in quantum physics as it contains, in principal, all information about the respective physical system. Particularly the set of singular points of Green’s function coincides with the spectrum. Green’s function is in fact the integral kernel of the resolvent of the Hamiltonian which is regarded as an integral operator. Sometimes, however, the integral kernel should be interpreted in the distributional sense. Green’s function is also closely related to the heat kernel or to the propagator. Namely, Green’s function is the Laplace transform of the heat kernel.

Green’s function can be explicitly expressed in a compact form for some quantum systems which are usually distinguished by their symmetry properties. As a rule, such systems frequently enjoy rotational symmetry. If so, this also opens the way to an alternative construction of Green’s function based on the method of separation of variables. The problem then effectively reduces to a one-dimensional one. Finally one deals with a positive second-order ordinary differential operator of Sturm-Liouville type though on the half-line rather than on a bounded interval. This is a substantial simplification since the construction of Green’s function for a Sturm-Liouville operator is commonly known and, in fact, this is a text-book matter. Thus in distinguished solvable cases the radial part of Green’s function can be expressed in terms of appropriate special functions. The full Hamiltonian depending on both radial and angle variables is then expressed as a sum with the summation running over eigenmodes of the spherical part of the Hamiltonian. Equaling the compact form of Green’s function to the sum obtained via the method of separation of variables leads to an addition formula for the involved special functions.

This procedure can be successfully applied to the Hamiltonian of the hydrogen atom, resulting in an addition formula for the Whittaker functions Mκ,μ and Wκ,μ. The formula turns out to be a generalization of Gegenbauer’s Addition Theorem in some respect. To the best of author’s knowledge, this possibility has not been exploited yet and still remains overlooked. And this is despite the fact that a compact formula for Coulomb Green’s function has been derived by Hostler rather long time ago in Ref. 10. As a companion of the addition formula for the Whittaker functions we further derive another addition formula concerning the Laguerre polynomials. There already exists a well-known addition formula for the Laguerre polynomials but that one reported here is completely different.

From the mathematical point of view, the formula for the Whittaker functions cannot be considered fully satisfactory, however, as the resulting sum involves only the Whittaker functions with the parameter μ = 1/2. A more general formula for arbitrary parameters κ and μ seems to be lacking. Nevertheless here we present some partial results in this direction which indicate that the derived formula could be further generalized. In addition to the parameter κ, with μ being restricted to the values 1/2 modulo integers, the formula depends on the radial variables r0 and r and on an angle γ. In the particular case γ = π we show that there exists an addition formula admitting general values of both κ and μ. Further we derive a summation formula for the Whittaker functions Wκ,μ only and another one for the Whittaker functions Mκ,μ. Again, in both cases, the parameters κ and μ can take arbitrary values.

The paper is organized as follows. In Sec. II we summarize some known formulas and results which are essential for the solution of our problem. The main result of the paper, namely a derivation of an addition formula for the Whittaker functions, is the content of Sec. III. Section IV is devoted to an addition formula for the Laguerre polynomials. Finally, Sec. V contains some complementary results suggesting that further generalizations could be possible, as discussed above.

First let us recall the definition of the Whittaker functions and summarize several useful formulas related to them.1,9 The Whittaker functions are defined in terms of the confluent hypergeometric functions,
Mκ,μ(r)er/2r1/2+μF11μκ+12;2μ+1;r,
(1)
Wκ,μ(r)er/2r1/2+μUμκ+12,2μ+1,r.
(2)
For the derivatives we shall use the notation
Wν,1/2(x)Wν,1/2(x)x,Mν,1/2(x)Mν,1/2(x)x.
The modified Bessel functions Kν, Iν are particular cases of the Whittaker functions,’
Kν(z)=π2zW0,ν(2z),Iν(z)=122νΓ(ν+1)2zM0,ν(2z).
(3)
Furthermore,
M0,1/2(z)=2sinhz2,W0,1/2(z)=ez/2.
(4)
For nZ+, the functions Wn+(α+1)/2,α/2 and Mn+(α+1)/2,α/2 are linearly dependent,
Mn+(α+1)/2,α/2(z)=(1)nΓ(α+1)Γ(n+α+1)Wn+(α+1)/2,α/2(z).
Moreover, the generalized (associated) Laguerre polynomials are related to the Whittaker functions,
Lnα(z)=(1)nn!z(α+1)/2Wn+(α+1)/2,α/2(z).
Therefore, for nZ+,
Wn,1/2(z)=(1)n+1n!Mn,1/2(z)
(5)
and, for n ≥ 1,
Mn,1/2(z)=1nez/2zLn11(z).
(6)
Regarding the asymptotic forms, we have
er/2Wκ,μ(r)=rκ1+μκ+12μ+κ12r+O1r2,asr
(7)
(see, for example, Ref. 2 and equation 13.5.2 in Ref. 1). Furthermore (Ref. 4, Eqs. 13.20.1, 13.20.2), when μ in C, Re(μ) > 0, and κC is fixed,
Mκ,μ(z)=zμ+1/21+O(μ1)
(8)
uniformly for z in a bounded region in C, and
Wκ,μ(x)=Γ(κ+μ)πx41/2μ1+O(μ1)
(9)
uniformly for bounded positive values of x [one can also consult (Ref. 13, Chap. 7, Sec. 11.1) or Ref. 5].
Next let us shortly recall, without going into all details, a standard construction of Green’s function of a Sturm-Liouville operator for it is essential for our purposes. Thus we consider a second-order ordinary differential operator
Lf(x)p(x)f(x)+q(x)f(x)
on an interval which can be bounded or unbounded. Here p(x) > 0 and q(x) ≥ 0 are sufficiently regular functions. At the finite endpoints one imposes mixed boundary conditions or, if the interval is unbounded, one requires functions from the domain of L to be square integrable on a neighborhood of infinity. One assumes that L with properly chosen boundary conditions is positive definite. To describe Green’s function of L one finds two nontrivial solutions v0, v1 of the differential equation
pvjpvj+qvj=0,j=0,1,
on the given interval such that v0 satisfies the boundary condition at the left endpoint (or minus infinity) only while v1 satisfies the boundary condition at the right endpoint (or plus infinity) only. Then L−1 is an integral operator with the integral kernel
G(x,y)=1pwϑ(yx)v0(x)v1(y)+ϑ(xy)v0(y)v1(x)
where wv0v1v1v0 is the Wronskian of v0 and v1. Note that p(x)w(x) is in fact a constant function. Here and in the sequel ϑ denotes the Heaviside step function.
It may be instructive to illustrate the procedure leading to an addition formula, as described in Sec. I, on the well-known example of the operator −∇2 + k2, k > 0, in R2. Naturally, the partial differential operator is expressed in polar coordinates. Using the method of separation of variables one finds that Green’s function can be written in the form
G(r,φ;r0,φ0)=12πrr0f0(r,r0)+2n=1fn(r,r0)cos(n(φφ0)).
(10)
The functions fn(r, r0) can be obtained as solutions of the Sturm-Liouville problem in the radial variable, as described above, and we get
fn(r,r0)=12k(2n)!Γn+12×ϑ(rr0)W0,n(2kr)M0,n(2kr0)+ϑ(r0r)M0,n(2kr)W0,n(2kr0).
The RHS can be also expressed in terms of the modified Bessel functions, see Ref. 3. At the same time, a compact formula for Green’s function is well known,
G(r,r0)=12πK0(k|rr0|).
We can let φ0 = 0 and, by comparison, we obtain an addition formula for the modified Bessel functions,
I0(kr0)K0(kr)+2n=1In(kr0)Kn(kr)cos(nφ)=K0(kR)for 0r0<r,
(11)
where
R=R(φ)r2+r022rr0cos(φ).
(12)
As a matter of fact, formula (11) is a corollary of substantially more general Graf’s Addition Theorem (Ref. 1, Eq. 9.1.79).
Let us rewrite just to have a comparison to Hostler’s result (20) which is mentioned below. From Ref. 10 we deduce that
I0(v)K0(u)+2n=1In(v)Kn(u)=K0(uv)for0v<u.
Let (r = |r|, r0=r0)
xr+r0+|rr0|,yr+r0|rr0|.
(13)
Then
2πG(r,r0)=I0ky2K0kx2+2n=1Inky2Knkx2.
One can proceed very analogously in case of the operator −∇2 + k2, k > 0, in R3. Doing so one obtains a particular case of Gegenbauer’s Addition Theorem. We are not going to discuss this case, however. Instead, in Sec. III, we will focus on the operator −∇2g/|r| + k2, g > 0 and k > 0, in R3. Nevertheless let us recall what Gegenbauer’s theorem claims if specialized to the modified Bessel functions. For 0 ≤ r0 < r and νR,
2νΓ(ν)(rr0)νn=0(ν+n)Kν+n(r)Iν+n(r0)Cn(ν)(cos(γ))=Kν(R)Rν
(14)
with R defined in Ref. 11, see (Ref. 16, Sec. II.4) [and also Ref. 1 (Eq. 9.1.80)]. Here Cn(ν)(z) are the Gegenbauer polynomials.
Another addition formula which is crucial for our purposes is Spherical Harmonic Addition Theorem, also called Legendre Addition Theorem. Recall that the spherical harmonics are defined for Z+ (the set of non-negative integers), mZ, |m| ≤ ,
Ym(θ,φ)(2+1)(m)!4π(+m)!Pm(cos(θ))eimφ.
(15)
Here θ ∈ [ 0, π ], φ ∈ [ 0, 2π ] are coordinates on the unit sphere S2 and Pm(z) is the associated Legendre polynomial. We have
Pm(z)=(1)m(m)!(+m)!Pm(z).
(16)
and
Ym(θ,φ)=(1)mYm(θ,φ)̄.
The spherical harmonics {Ym} form an orthonormal basis in L2(S2, dΩ) and
2Ym(θ,φ)=(+1)r2Ym(θ,φ).
Spherical Harmonic Addition Theorem tells us that, for every Z+,
m=Ym(θ,φ)Ym(θ0,φ0)̄=m=(1)mYm(θ,φ)Ym(θ0,φ0)=2+14πP(cos(γ))
(17)
where
cos(γ)cos(θ)cos(θ0)+sin(θ)sin(θ0)cos(φφ0),
that is cos(γ) = n · n0, nsin(θ)cos(φ),sin(θ)sin(φ),cos(θ) and n0 is defined similarly. P(z)P0(z) are the Legendre polynomials. Referring to (15) and (16), Eq. (17) means that
P(cos(γ))=P(cos(θ))P(cos(θ0))+2m=1(m)!(+m)!Pm(cos(θ))Pm(cos(θ0))cos(m(φφ0)).
This is a classical result with a long history,8 and a dozen different proofs of it have been provided, some of them quite intricate.12 A straightforward derivation, as encountered in physical literature, is based on symmetry considerations and can be briefly rephrased as follows. Observe that
P̃(n,n0)m=Ym(θ,φ)Ym(θ0,φ0)̄
is the integral kernel of the orthogonal projection in L2(S2, dΩ) onto the eigenspace of minus the Laplace-Beltrami operator on S2 (denoted as ΔS2) corresponding to the eigenvalue ( + 1). Thus we have (the differential operator acts in variables θ and φ)
ΔS2P̃(n,n0)=(+1)P̃(n,n0).
(18)
Owing to rotational symmetry P̃(n,n0) should depend on the distance of the points n, n0S2 only which in turn is a function of n · n0 = cos(γ). Writing P̃(n,n0)=f(nn0) Eq. (18) reduces to the ordinary second-order differential equation
(1z2)f(z)+2zf(z)(+1)f(z)=0.
A general solution has the form f(z) = c1P(z) + c2Q(z). Here Q(z) is the Legendre function of the second kind which is another independent solutions of the differential equation and is known to be singular for z = ±1. Therefore P̃(n,n0)=c1P(nn0). The multiplicative constant is easily found to be c1 = (2 + 1)/(4π) when letting n = n0 = (0, 0, 1) and taking into account that P(1) = 1 and Pm(1)=0 for m ≠ 0.
Using spherical coordinates r > 0, θ ∈ [0, π], φ ∈ [0, 2π], we denote
r=rsin(θ)cos(φ),rsin(θ)sin(φ),rcos(φ).
We write the Hamiltonian of the hydrogen atom in a dimensionless form as H = −∇2g/r, g > 0, and we wish to apply the procedure leading to an addition formula, as described in Sec. I, to the operator
H+k2=1r2rr2r1r21sin(θ)θsin(θ)θ+1sin2(θ)2φ2gr+k2,k>0.
While the continuous spectrum of H coincides with the positive real half-line, the discrete spectrum consists of eigenvalues En = −g2/4n2, nN, the multiplicity of En equals n2. The corresponding normalized eigenfunctions are
ψn,,m(r,θ,φ)=g3/2n+2(n1)!2(n+)!(gr)expgr2nLn12+1grnYm(θ,φ),
where nN, Z+ and mZ are the principal, the azimuthal and the magnetic quantum number, respectively, and |m| ≤ n − 1.
Application of the method of separation of variables to this operator again leads to the Sturm-Liouville problem in the radial variable whose solution is rather straightforward, as outlined in Sec. II, and has already been described in the literature.14,15 We have
G(r,θ,φ,r0,θ0,φ0)=1rr0=0m=12k(2+1)!Γ+1g2k×ϑ(rr0)Mg/(2k),+1/2(2kr0)Wg/(2k),+1/2(2kr)+ϑ(r0r)Mg/(2k),+1/2(2kr)Wg/(2k),+1/2(2kr0)Ym(θ,φ)Ym(θ0,φ0)̄.
With the aid of16 this equation can be further simplified,
G(r,θ,φ,r0,θ0,φ0)=18πkrr0=01(2)!Γ+1g2k×ϑ(rr0)Mg/(2k),+1/2(2kr0)Wg/(2k),+1/2(2kr)+ϑ(r0r)Mg/(2k),+1/2(2kr)Wg/(2k),+1/2(2kr0)P(cos(γ)).
(19)
Notably, there exists a remarkable compact formula for Green’s function due to Hostler,10,
G(r,θ,φ,r0,θ0,φ0)=14πRΓ1g2k×Mg/(2k),1/2(ky)Wg/(2k),1/2(kx)Mg/(2k),1/2(ky)Wg/(2k),1/2(kx)
(20)
where in this case we have
R|rr0|=r2+r022rr0sin(θ)sin(θ0)cos(φφ0)+cos(θ)cos(θ0),
and
xr+r0+|rr0|,yr+r0|rr0|
(21)
(formally the same equations as in Ref. 12 but now the dimension is 3 rather than 2).
Comparing (19) to (20) while still using notation (21) we have, for 0 ≤ r0 < r,
1krr0=01(2)!Γ+1g2kMg/(2k),+1/2(2kr0)Wg/(2k),+1/2(2kr)P(cos(γ))=Γ1g2k2RMg/(2k),1/2(ky)Wg/(2k),1/2(kx)Mg/(2k),1/2(ky)Wg/(2k),1/2(kx).
After substitution g = 2 and rescaling rr/(2k), r0r0/(2k) we get an addition formula for the Whittaker functions. Moreover, the parameter κ can be extended to complex values by analyticity.

Theorem 1.
For 0 ≤ r0 < r, κC\N and γR,
1rr0=0Γ(+1κ)Γ(1κ)(2)!Mκ,+1/2(r0)Wκ,+1/2(r)P(cos(γ))=1RMκ,1/2y2Wκ,1/2x2Mκ,1/2y2Wκ,1/2x2
(22)
where R = R(γ), see Ref. 11, and
x=r+r0+R,y=r+r0R.
(23)

One has to exclude the values κN. As a matter of fact, the equation holds for these values, too, but it should be achieved in the limit after the singular terms in the equation have been combined. For instance, for κ = 1 we have
1rr0=1(1)!(2)!M1,+1/2(r0)W1,+1/2(r)Pcos(γ)=κ1RMκ,1/2y2Wκ,1/2x2Mκ,1/2y2Wκ,1/2x21rr0Mκ,1/2(r0)Wκ,1/2(r)κ=1.

Remark 2.
Let us check two particular cases. For γ = 0 (hence R = rr0, x = 2r, y = 2r0) we have
1rr0=0Γ(+1κ)(2)!Mκ,+1/2(r0)Wκ,+1/2(r)=Γ(1κ)rr0Mκ,1/2(r0)Wκ,1/2(r)Mκ,1/2(r0)Wκ,1/2(r),
(24)
and for γ = π [hence R = r + r0, x = 2(r + r0), y = 0] we have
1rr0=0(1)Γ(+1κ)(2)!Mκ,+1/2(r0)Wκ,+1/2(r)=Γ(1κ)r+r0Wκ,1/2(r+r0).
(25)

Remark 3.
From (25) one can derive a summation formula for Whittaker functions which seems to be also new. For κC\N and zC,
1z=0(1)Γ(+1κ)Γ(1κ)(2)!Mκ,+1/2(z)=ez/2
(26)
(where the singularity at z = 0 on the LHS is removable). It can be proven by exploring the asymptotic behavior of both sides of (25), as r. We have, in virtue of (7),
1r+r0Wκ,1/2(r+r0)=er/2r0/2r1+κ1+κκ2(1κ)r0r+O1r2
and
1rWκ,+1/2(r)=er/2r1+κ1+(+1κ)(+κ)r+O1r2.
Thus, writing z instead of r0, we see that (26) holds for z > 0. But the asymptotic behavior of Mκ,+1/2(z) for large, as recalled in Eq. (8), which is locally uniform in z implies that the LHS is an entire function of z. Since the same is true for the RHS we conclude that (26) must hold for all zC.

Let us point out a relation of Theorem 1 to Gegenbauer’s Addition Theorem. Confining ourselves to the value κ = 0 (corresponding to g = 0) in (22) we obtain a simplified equation
=02+1rr0K+1/2r2I+1/2r02P(cos(γ))=1ReR/2=1πRK1/2R2.
(27)
This is a particular case of Gegenbauer’s Addition Theorem, however, see Eq. (14) with ν = 1/2. To derive (27) from (22) we have used (3) and (4) and also the equation
M0,1/2y2W0,1/2x2M0,1/2y2W0,1/2x2=eR/2=RπK1/2R2.
Moreover, note that Cn(1/2)(z)=Pn(z).
As for Eq. (26), letting κ = 0 we get
π2r=0(1)(2+1)I+1/2(z)=ez.
This is a particular case of the identity
2νΓ(ν)=0(+ν)C(ν)(γ)I+ν(z)=zνeγz,
which is well known, see [Ref. 7, Sec. 7.15(1)]; note that C(1/2)(1)=(1).

Theorem 4.
For nN, 0 ≤ r0 < r and γR,
=0n1(2+1)(n1)!(n+)!(rr0)Ln12+1(r)Ln12+1(r0)P(cos(γ))=12RxLn11x2Lny2yLn11y2Lnx2
(28)
where x, y are defined in (23), with R = R(γ), see Eq. (12).

Remark.
Note that formula (28) is completely different from the well known addition formula for the Laguerre polynomials which claims that (Ref. 1, Eq. 22.12.6)
j=0nLjα(u)Lnjβ(v)=Lnα+β+1(u+v).

Proof.
Recall that
ResΓ(z);z=n=(1)nn!,nZ+,
whence
ResΓ1g/(2z);z=En=(1)ng22n!n2,nN.
If we substitute k=z in Hostler’s formula (20) then the residue of Green’s function at z = En equals
(1)ng28πn!n2RMn,1/2gy2nWn,1/2gx2nMn,1/2gy2nWn,1/2gx2n.
(29)
The residue also equals minus the projection Pn onto the eigenspace corresponding to eigenvalue En. Pn is an integral operator with the integral kernel
Pn(r,θ,φ,r0,θ0,φ0)==0n1m=ψn,,m(r,θ,φ)ψn,,m(r0,θ0,φ0)̄.
(30)
Hence, for nN,
(1)n+1n(n1)!gRMn,1/2gy2nWn,1/2gx2nMn,1/2gy2nWn,1/2gx2n=expg(r+r0)2n=0n1(2+1)(n1)!(n+)!g2rr0n2×Ln12+1grnLn12+1gr0nP(cos(γ)).
(31)
In regard of (5) and (6), we derive
Mn,1/2(z)=ddz1nez/2ezzLn11(z)=12Mn,1/2(z)+1nez/2ddzezzLn11(z)=12Mn,1/2(z)+ez/2Ln(z).
Here we have used the formula [Ref. 11, Eq. (9.12.8)]
ddzezzαLn1α(z)=nezzα1Lnα1(z).
Furthermore,
(1)n+1g28πRn!n2Mn,1/2gy2nWn,1/2gx2nMn,1/2gy2nWn,1/2gx2n=g28πRn2Mn,1/2gx2nMn,1/2gy2nMn,1/2gy2nMn,1/2gx2n=g316πRn4expg2n(r+r0)xLn11gx2nLngy2nyLn11gy2nLngx2n.
(32)
After rescaling r → (n/g)r, r0 → (n/g)r0, we get from (31) and (32)
=0n1(2+1)(n1)!(n+)!(rr0)Ln12+1(r)Ln12+1(r0)P(cos(γ))=(1)n+1(n1)!Rexpr+r02Mn,1/2y2Wn,1/2x2Mn,1/2y2Wn,1/2x2=12RxLn11x2Lny2yLn11y2Lnx2.
This concludes the proof.□
Let us check two particular cases of (28). For γ = 0 we have
=0n1(2+1)(n1)!(n+)!(rr0)Ln12+1(r)Ln12+1(r0)=1rr0rLn11(r)Ln(r0)r0Ln11(r0)Ln(r)
and for γ = π we get
=0n1(1)(2+1)(n1)!(n+)!(rr0)Ln12+1(r)Ln12+1(r0)=Ln11(r+r0),
Note that Ln(0) = 1 [and Ln1(0)=n+1].

After the shift nn + 1 one observes, in these two particular cases, that both sides are symmetric polynomials in r and r0, and therefore r, r0 can be replaced by arbitrary complex variables.

Corollary 5.
For every nZ+ and all u,vC,
=0n(2+1)(n)!(n++1)!(uv)Ln2+1(u)Ln2+1(v)=1uvuLn11(u)Ln(v)vLn11(v)Ln(u)
(33)
and
=0n(1)(2+1)(n)!(n++1)!(uv)Ln2+1(u)Ln2+1(v)=Ln1(u+v).
(34)

As a short digression let us note that the projection Pn, as introduced in Eq. (30) in the proof, has already been discussed in the literature .3 Using expression for the residue of Green’s function, see (29) and (32), we get, for n ≥ 1,
Pn(r,θ,φ,r0,θ0,φ0)=g316πRn4expg2n(r+r0)×xLn11gx2nLngy2nyLn11gy2nLngx2n.
The diagonal in fact does not depend on angles and equals
Pn(r,θ,φ,r,θ,φ)=g38πn4expgrn×LngrnLn11grngrnLngrnLn22grn+grnLn11grn2.
Pn(r,θ,φ,r0,θ0,φ0) is called the density function in Ref. 3, and
Dn(r)4πr2Pn(r,θ,φ,r,θ,φ)
is called the radial distribution function. It holds true that
4π0Pn(r,θ,φ,r,θ,φ)r2dr=n2
meaning that
0exp(r)Ln(r)Ln11(r)rLn(r)Ln22(r)+rLn11(r)2r2dr=2n3.

The following proposition presents a summation formula for the Whittaker functions Wκ,μ and is a straightforward corollary of Theorem 8 below which in turn generalizes the addition formula (25). Nevertheless this proposition should be proven independently because, conversely, it is used in the Proof of Theorem 8.

Proposition 6.
For every nN, r > 0, κ,μC,μ ≠ −1, −2, −3, …,
=0n(1)n2μ+2(2μ+)n+1Wκ,μ+(r)=(1)nrn/2Wκn/2,μ+n/2(r).
(35)

Proof.
We shall proceed by induction in n. For n = 0 the equation is trivial. For n = 1 this is a well known identity [for instance, this is a combination of Eqs. 9.234 ad(1) and ad(2) in Ref. 9]
Wκ,μ+1(r)Wκ,μ(r)=2μ+1rWκ1/2,μ+1/2(r).
(36)
Suppose n > 0. Set, for k = 0, 1, …, n,
A(k)=0nk1(1)n2μ+2(2μ+)n+1Wκ,μ+(r)+(1)n+kn1kWκ,μ+nk(r)(2μ+nk)n1r=nkn1(1)n12μ+2+1(2μ++1)nWκ1/2,μ++1/2(r).
In particular,
A(0)==0n1(1)n2μ+2(2μ+)n+1Wκ,μ+(r)+(1)nWκ,μ+n(r)(2μ+n)n==0n(1)n2μ+2(2μ+)n+1Wκ,μ+(r).
Thus A(0) coincides with the LHS of Eq. (35).
We claim that A(k + 1) = A(k) for k = 0, 1, …, n − 1. Suppose 0 ≤ k < n. Then
A(k+1)A(k)=(1)n+knk+12μ+2n2k2(2μ+nk1)n+1Wκ,μ+nk1(r)+(1)n+k+1n1k+1Wκ,μ+nk1(r)(2μ+nk1)n+(1)n+k+1n1kWκ,μ+nk(r)(2μ+nk)n+1r(1)n+kn1k2μ+2n2k1(2μ+nk)nWκ1/2,μ+nk1/2(r).
With the aid of (36) one finds that this expression equals
(1)n+k+1n1k+1(2μ+2nk1)Wκ,μ+nk1(r)(2μ+nk1)n+1(1)n+k+1nk+1(2μ+2n2k2)Wκ,μ+nk1(r)(2μ+nk1)n+1+(1)n+k+1n1k(2μ+nk1)Wκ,μ+nk1(r)(2μ+nk1)n+1.
Now it is elementary to see that the expression actually equals 0, and therefore A(k + 1) − A(k) = 0.
In A(n) we apply substitutions κ′ = κ − 1/2, μ′ = μ + 1/2, and obtain, by the induction hypothesis,
A(n)=1r=0n1(1)n12μ+2(2μ+)nWκ,μ+(r)=1r(1)n1r(n1)/2Wκ(n1)/2,μ+(n1)/2(r)=(1)nrn/2Wκn/2,μ+n/2(r).
Since A(0) = A(n), the identity follows.□

Remark 7.
From (35) it follows that, for nZ+,
=0n(1)n2μ+2(2μ+)n+1=δn,0.

Indeed, one just has to recall (7) and compare the asymptotic expansions of both sides of (35) as r.

Here is a generalization of the addition formula for the Whittaker functions (25).

Theorem 8.
For 0 ≤ r0 < r, κ,μC, Re μ > 0,
1(rr0)μ+1/2=0(1)(μκ+1/2)(+2μ)!Mκ,+μ(r0)Wκ,+μ(r)=1(r+r0)μ+1/2Wκ,μ(r+r0).
(37)

Remark.

The equation for r0 = 0 should be understood as a limiting case of (37), and it is trivial.

Remark 9.
Regarding the convergence of the series, it is guaranteed by the assumption 0 ≤ r0 < r. Let us shortly analyze this point. Put
tr+r0rr0μ+1/2(μκ+1/2)Mκ,+μ(r0)Wκ,+μ(r)(+2μ)!Wκ,μ(r+r0),
so that (37) becomes =0(1)t=1. Referring to (8) and (9), we replace the Whittaker functions in the numerator by their leading asymptotic terms and get
t=22μ1(r+r0)μ+1/2πr2μWκ,μ(r+r0)(μκ+1/2)Γ(+κ+μ)!(+2μ)4r0r1+O(1).
Next we do the same for the factorial and the Pochhammer symbols in the expression while using Stirling’s asymptotic formula. We find that for large ,
t=(r+r0)μ+1/2Γ(μκ+1/2)r2μWκ,μ(r+r0)2μ1r0r1+O(1).
This obviously guarantees the convergence.

But on the other hand, the series turns out to be numerically quite unstable for large values of μ. In such a case we can be dealing with an alternating series with many summands in its beginning attaining huge values. Then significant cancellations of the terms necessarily happen. From the numerical point of view this is a troublesome situation. As an example let us consider the case with μ = 20, κ = 1, r0 = 1, and r = 2. The Computer Algebra System Mathematica, as of version 14.0.0, gives the values t0 = 1.072 39 × 107, t145 = 3214.65, and is not capable to compute t for higher indices. Nonetheless replacing the involved Whittaker functions by their leading asymptotic terms for μ large one finds that the first index for which t attains a value smaller than 0.1 is = 168. This shows that this concrete series starts to rapidly converge to its final sum only for very large summation indices.

Proof.
In view of (1) and (2) the equation can be rewritten as
=0(1)(μκ+1/2)(+2μ)!(rr0)×F11+μκ+12;2+2μ+1;r0U+μκ+12,2+2μ+1,r=Uμκ+12,2μ+1,r+r0.
(38)
We have (Ref. 1, Eq. 13.4.21)
jrjU(a,b,r)=(1)j(a)jU(a+j,b+j,r),jZ+.
whence
Ua,b,r+r0=n=0(1)n(a)nn!U(a+n,b+n,r)r0n.
Using, in (38), this power expansion as well as the power expansion of F11 we get
=0(1)(1/2+μκ)(l+2μ)!j=0(+μκ+1/2)j(2l+2μ+1)jj!×U+μκ+12,2+2μ+1,rrr0+j=n=0(1)n(μκ+1/2)nn!Uμκ+n+12,2μ+n+1,rr0n.
Comparing the coefficients at the same powers of r0 we obtain an equivalent countable system of equations, with nZ+,
=0n(1)n(2+2μ)(+2μ)n+1rUμ+κ+12,2μ+2+1,r=(1)nUμκ+n+12,2μ+n+1,r.
Expressing reversely the confluent hypergeometric function U in terms of the Whittaker function W we get
=0n(1)n2μ+2(2μ+)n+1Wκ,μ+(r)=(1)nrn/2Wκn/2,μ+n/2(r).
This identity has been proven in Proposition 6.□

Finally we prove a summation formula for the Whittaker functions Mκ,μ.

Proposition 10.
For zC, κC, μ > 0 and c ∈ [−1, 1],
=0(μκ+1/2)(2μ)2zF11+μκ+12;2+2μ+1;zC(μ)(c)=F11μκ+12;μ+12;1+c2z,
(39)
or, if rewritten in terms of the Whittaker functions, with γ ∈ [−π, π],
zμ1/2=0(μκ+1/2)(2μ)2Mκ,+μ(z)C(μ)(cos(γ))=ez/2F11μκ+12;μ+12;cos2γ2z.
(40)

Remark 11.

Here we tacitly assume that zμ−1/2Mκ,+μ(z) is understood as an entire function – first defined on the positive half-line and then admitting an unambiguous continuation to the entire complex plane as an analytic function.

The particular case γ = 0 gives
=0(μκ+1/2)(+2μ)!Mκ,+μ(z)=ez/2F11μκ+12;μ+12;z,
and the particular case γ = π gives
zμ1/2=0(1)(μκ+1/2)(l+2μ)!Mκ,+μ(z)=ez/2.
(41)
Note that
C(μ)(1)=(2μ)!,C(μ)(1)=(1)(2μ)!.
We remark that (41) is a generalization of (26).

Remark 12.
Convergence of the series in (40) can be justified by exploring the asymptotic behavior of the summands. Regarding the Gegenbauer polynomials, an elaborate asymptotic expansion for large orders is derived in a recent paper.6 For our purposes a comparatively simple approach is sufficient. From equation (3.30) in Ref. 2 one infers that
C(ν)(cos(θ))C(ν)(1)=2Γν+12πΓ(ν)×0π/2cos2ν1(φ)cosarccoscos(θ)1sin2(θ)cos2(φ)(1sin2(θ)cos2(φ))/2dφ
for ν > 0, 0 ≤ θπ and Z+. It follows that
|C(ν)(cos(θ))|C(ν)(1)2Γν+12πΓ(ν)0π/2cos2ν1(φ)dφ=1.
Thus for c ∈ [−1, 1] we have
|C(μ)(c)|C(μ)(1)=(2μ)!=2μ1Γ(2μ)1+O(1)as.
Furthermore, for sufficiently large, surely for any such that + 2μ + 1 > |μκ + 1/2|, we can estimate
n0,|(+μκ+1/2)n|(2+2μ+1)n1,
and therefore
F11+μκ+12;2+2μ+1;zK(κ,μ)e|z|
where K(κμ) is a constant independent of and z although it may depend on κ and μ. The convergence now becomes obvious.

These estimates even show that the LHS of (39) is an entire function.

Proof.
Using the power series expansion for the hypergeometric series we have to show that
=0j=0(μκ+1/2)(2μ)(l+2μ)(+μκ+1/2)j(2+2μ+1)jj!C(μ)(c)z+j=n=0(μκ+1/2)n(μ+1/2)nn!1+c2nzn.
Comparing coefficients at the same powers of z we get an equivalent countable set of equations, with nZ+,
=0n(μκ+1/2)(2μ)2(+μκ+1/2)n(2+2μ+1)n(n)!C(μ)(c)=(μκ+1/2)n(μ+1/2)nn!1+c2n.
These equations can be simplified using straightforward manipulations,
=0n(2μ+2)(2μ)n++1(n)!C(μ)(c)=2n(μ)n(2μ)2nn!(1+c)n.
Recall that [Ref. 11, Eq. (9.8.19)]
C(μ)(c)=(2μ)!2F1,2μ+;μ+12;1c2.
After the substitution c = 1 − 2w we get the equation
=0n(2μ+2)(2μ)n++1(n)!(2μ)!2F1,2μ+;μ+12;w=22n(μ)n(2μ)2nn!(1w)n.
And after the power series expansion of the hypergeometric function and interchanging the order of summation we obtain
j=0n=jn(2μ+2)(2μ)n++1(n)!(2μ)!()j(2μ+)j(μ+1/2)jj!wj=22n(μ)n(2μ)2nn!(1w)n.
Comparing coefficients at the same powers of w leads to the equations
=jn(2μ+2)(2μ)n++1(n)!(2μ)!()j(2μ+)j(μ+1/2)j=(1)j(nj)!22n(μ)n(2μ)2n,0jn,
and this can be further rewritten,
=jnnjj2μ+2(2μ++j)nj+1=(μ+1/2)j(μ+1/2)n,0jn.
Shifting the summation index + j and applying the substitutions n = N + j, ν = μ + j, we obtain an identity [Eq. (42)] which is proven in Lemma 13 below thus concluding this proof.□

Lemma 13.
For NZ+ and νC, ν ≠ 0, −1, −2, …,
=0NN2ν+2(2ν+)N+1=1(ν+1/2)N.
(42)

Proof.
For NZ+, consider the expression
=0NN(2ν+2)(ν+1/2)N(2ν+)N+1
(43)
which is regarded as a function of νC. This is a meromorphic function on C, and its limit, as ν, equals
=0NN12N=1.
(44)
The poles are located at the points ν = 0, −1/2, −1, …, − (2N − 1)/2, − N, and all of them are of first order.
Owing to the factor (ν + 1/2)N the singularities are removable for
ν=1/2,3/2,,(2N1)/2.
Hence it suffices to check the poles at ν = 0, −1, …, − N. The residue at a pole ν = −t, 0 ≤ tN, equals
(t+1/2)N=max{2tN,0}min{2t,N}N2t+2(2t+)2t(1)N2t+=(t+1/2)N=max{2tN,0}min{2t,N}(1)N2t+2(2t)!(N2t+)!.
Further we omit the nonzero factor 2 (−t + 1/2)N. Shifting the summation index t + , for tN/2 we get
l=2tNN(1)Nt+(2t)!(N2t+)!=(1)tN!=tNNt(1)Nt+Nt=0.
Just note that the summands are odd in . For tN/2 we similarly get
=02t(1)Nt+(2t)!(N2t+)!=(1)tN!=tt(1)Nt+Nt=0.
Hence all singularities are removable and therefore expression (43) equals 1 identically from (44) and Liouville’s Theorem.□

The author acknowledges partial support by European Regional Development Fund Project “Center for Advanced Applied Science” Grant No. CZ.02.1.01/0.0/0.0/16_019/0000778. The author is indebted to the reviewer for numerous comments helpful in improving the quality of the paper, and in particular for Remark 9.

The author has no conflicts to disclose.

Pavel Šťovíček: Investigation (lead).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1.
M.
Abramowitz
and
I. A.
Stegun
,
Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables
(
Dover Publications
,
New York
,
1972
).
2.
R.
Askey
,
Orthogonal Polynomials and Special Functions
(
SIAM
,
Philadelphia
,
1975
).
3.
S. M.
Blinder
, “
Generalized Unsöld theorem and radial distribution function for hydrogenic orbitals
,”
J. Math. Chem.
14
,
319
324
(
1993
).
4.
NIST Digital Library of Mathematical Functions, edited by
F. W. J.
Olver
,
A. B.
Olde Daalhuis
,
D. W.
Lozier
,
B. I.
Schneider
,
R. F.
Boisvert
,
C. W.
Clark
,
B. R.
Miller
,
B. V.
Saunders
,
H. S.
Cohl
, and
M. A.
McClain
, https://dlmf.nist.gov/, Release 1.1.11 of 15 September 2023.
5.
T. M.
Dunster
, “
Uniform asymptotic expansions for the Whittaker functions Mκ,μ(z) and Wκ,μ(z) with μ large
,”
Proc. R. Soc., Ser. A
477
,
0360
(
2021
).
6.
T. M.
Dunster
, “
Uniform asymptotic expansions for Gegenbauer polynomials and related functions via differential equations having a simple pole
,”
Constr. Approx.
(published online) (
2023
).
7.
A.
Erdélyi
,
W.
Magnus
,
F.
Oberhettinger
, and
F. G.
Tricomi
,
Higher Transcendental Functions
(
McGraw-Hill Book Company, Inc.
,
New York; Toronto, London
,
1953
),
Vol. II
.
8.
N. M.
Ferrers
,
An Elementary Treatise on Spherical Harmonics and Subjects Connected with Them
(
Macmillan
,
London
,
1877
).
9.
I. S.
Gradshteyn
and
I. M.
Ryzhik
, in
Table of Integrals, Series, and Products
, edited by
A.
Jeffrey
and
D.
Zwillinger
(
Academic Press
,
Amsterdam
,
2007
).
10.
L. C.
Hostler
and
R. H.
Pratt
, “
Coulomb Green’s function in closed form
,”
Phys. Rev. Lett.
10
,
469
470
(
1963
).
11.
R.
Koekoek
,
P. A.
Lesky
, and
R. F.
Swarttouw
,
Hypergeometric Orthogonal Polynomials and Their Q-Analogues
(
Springer
,
Berlin
,
2010
).
12.
K.
Maleček
and
Z.
Nádeník
, “
On the inductive proof of Legendre addition theorem
,”
Stud. Geophys. Geod.
45
,
1
11
(
2001
).
13.
F. W. J.
Olver
and
A. K.
Peters
,
Asymptotics and Special Functions
(
Wellesley
,
1997
).
14.
A.
Rauh
, “
On the singularities of the nonrelativistic Coulomb Green’s function
,”
Adv. Stud. Theor. Phys.
13
,
175
187
(
2019
).
15.
S. S.
Vasan
and
M.
Seetharaman
, “
The Coulomb Green’s function revisited
,”
Pramana
45
,
165
174
(
1995
).
16.
G. N.
Watson
,
A Treatise on the Theory of Bessel Functions
, 2nd ed. (
Cambridge University Press
,
Cambridge
,
1944
).
Published open access through an agreement with Czech Technical University in Prague Faculty of Nuclear Sciences and Physical Engineering