How to prove Gegenbauer's addition theorem?

Question

I asked this at: https://math.stackexchange.com/questions/210153/, but didn't get any reply, so I am trying here, since I actually need this in physics anyway.

How can one prove the following identity: $$ V_k(r_1, r_2) = {2k+1\over 2 r_1 r_2}\int_{|r_1 - r_2|}^{r_1+r_2} e^{-{r\over D}} P_k\left(r_1^2 - r^2 + r_2^2 \over 2 r_1 r_2 \right) d r= $$ $$ =(2k+1) {I_{k+{1\over2}}\left(r_<\over D\right) \over \sqrt{r_<}} {K_{k+{1\over2}}\left(r_>\over D\right) \over \sqrt{r_>}} $$ ? Here $r_<=\min(r_1, r_2)$, $r_>=\max(r_1, r_2)$, $I_\nu$ and $K_\nu$ are the modified Bessel functions of the first and second kind, $P_k(x)$ are Legendre polynomials and $k=0, 1, 2, 3, \dots$.

I verified the identity numerically, so I know that it works, but I didn't figure out how one could prove it. I know it has something to do with Gegenbauer's addition theorem and the equation 10.23.8 at [1]. The $V_k(r_1, r_2)$ is just the coefficient in the Legendre expansion of: $$ {e^{{|{\bf r}_1 - {\bf r}_2|}\over D}\over |{\bf r}_1 - {\bf r}_2|} =\sum_{k=0}^\infty V_k(r_1, r_2) P_k(\cos\theta) $$ Note: when $D\to\infty$, then $V_k$ becomes ${r_<^k \over r_>^{k+1}}$, which follows from the properties of the Bessel functions and this is just the well known multipole expansion [2].

[1] http://dlmf.nist.gov/10.23#E8

[2] http://en.wikipedia.org/wiki/Legendre_polynomials#Applications_of_Legendre_polynomials_in_physics

Answer 1

EDIT (11/18/2012) As I failed to obtain a proof using my initial reasoning (see below, at the end), let me offer a complete proof using a different approach.

Let us transform the integral in the left-hand side first, using the following substitution: $x=\frac{r_1^2-r^2+r_2^2}{2 r_1 r_2}$: $\int^{r_1+r_2}_{|r_1-r_2|}\exp(-\frac{r}{d})P_k(\frac{r_1^2-r^2+r_2^2}{2 r_1 r_2})dr=r_1 r_2\int^{1}_{-1}\frac{\exp\left(-\frac{\sqrt{r_1^2+r_2^2-2 r_1 r_2 x}}{d}\right)}{\sqrt{r_1^2+r_2^2-2 r_1 r_2 x}}P_k(x)dx=$ $=x_1 x_2 d^2\int^{1}_{-1}\frac{\exp\left(-\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}\right)}{d\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}}P_k(x)dx$, where $x_1=\frac{r_1}{d},x_2=\frac{r_2}{d}.$

Let us transform the right-hand side now (I assume, without loss of generality, that $r_1\leq r_2$): $(2 k+1)\frac{I_{k+\frac{1}{2}}(\frac{r_1}{d})}{\sqrt{r_1}}\frac{K_{k+\frac{1}{2}}(\frac{r_2}{d})}{\sqrt{r_2}}=(2 k+1)\frac{1}{d}\frac{I_{k+\frac{1}{2}}(x_1)}{\sqrt{x_1}}\frac{K_{k+\frac{1}{2}}(x_2)}{\sqrt{x_2}}=(2 k+1)\frac{1}{d}i_k(x_1)k_k(x_2)$, where $i_k(x)=\sqrt{\frac{\pi}{2 x}}I_{k+\frac{1}{2}}(x)$ and $k_k(x)=\sqrt{\frac{2}{\pi x}}K_{k+\frac{1}{2}}(x)$ are the modified spherical Bessel functions of the first and the second kind, respectively ([1],[2]). Thus, the equality in the question is equivalent to the following one:

$\frac{1}{2}\int^{1}_{-1}\frac{\exp\left(-\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}\right)}{\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}}P_k(x)dx=i_k(x_1)k_k(x_2).$ (1)

It is well known that modified spherical Bessel functions of the first and the second kind $i_k(z)$ and $k_k(z)$ satisfy the following differential equation ([3], eq.10.47).2:

$z^2\frac{d^2w}{dz^2}+2 z \frac{dw}{dz}-(z^2+k(k+1))w=0.$ (2)

If equality (1) is correct, its left-hand side (let us denote it $Q$) must satisfy eq.(2) with respect to $x_1$ and with respect to $x_2$. Let us prove that it satisfies eq.(2) with respect to $x_1$ (the proof for $x_2$ is identical), i.e. that

$0=\left(x_1^2\frac{\partial^2}{\partial x_1^2}+2 x_1 \frac{\partial}{\partial x_1}-(x_1^2+k(k+1))\right)\int^{1}_{-1}\frac{\exp\left(-\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}\right)}{\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}}P_k(x)dx=$

$=\int^{1}_{-1}\left(x_1^2\frac{\partial^2}{\partial x_1^2}+2 x_1 \frac{\partial}{\partial x_1}-(x_1^2+k(k+1))\right)\frac{\exp\left(-\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}\right)}{\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}}P_k(x)dx=$

$=\int^{1}_{-1}\left(x_1^2\frac{\partial^2}{\partial x_1^2}+2 x_1 \frac{\partial}{\partial x_1}-x_1^2\right)f P_k(x)dx-k(k+1)\int^{1}_{-1}f P_k(x)dx$, (3)

where

$f=f(x_1,x_2,x)=\frac{\exp\left(-\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}\right)}{\sqrt{x_1^2+x_2^2-2 x_1 x_2 x}}.$

As $\frac{d}{dx}((1-x^2)\frac{d}{dx}P_k(x))+k(k+1)P_k(x)=0$ ([4]),

$-k(k+1)\int^{1}_{-1}f P_k(x)dx=\int^{1}_{-1}f d((1-x^2)\frac{d}{dx}P_k(x))=-\int^{1}_{-1}(1-x^2)(\frac{d}{dx}P_k(x))f_x dx=-\int^{1}_{-1}(1-x^2)f_x d P_k(x)=\int^{1}_{-1}((1-x^2)f_x)_x P_k(x)dx=\int^{1}_{-1}(-2 x f_x+(1-x^2)f_{xx})P_k(x)dx$,

where differentiation by parts was used twice, and subscript $x$ means differentiation with respect to $x$. Thus, to prove eq.(3), it is sufficient to prove that

$\left(x_1^2\frac{\partial^2}{\partial x_1^2}+2 x_1 \frac{\partial}{\partial x_1}-x_1^2\right)f-2 x f_x+(1-x^2)f_{xx}=0.$ (4)

Derivatives of $f$ with respect to $x$ and $x_1$ (we denote them $f_x$ and $f_1$, respectively) can be obtained using the chain rule, and comparing the results, we obtain $f_1=f_x\frac{2 x_1-2 x_2 x}{-2 x_1 x_2}=f_x\frac{x_2 x-x_1}{x_1 x_2}=f_x(\frac{x}{x_1}-\frac{1}{x_2})$. The second derivative of $f$ with respect to $x_1$ equals $f_{11}=\frac{\partial}{\partial x_1}(f_x(\frac{x}{x_1}-\frac{1}{x_2}))=f_{1 x}(\frac{x}{x_1}-\frac{1}{x_2})+f_x(-\frac{x}{x_1^2})=(\frac{\partial}{\partial x}(f_x(\frac{x}{x_1}-\frac{1}{x_2})))(\frac{x}{x_1}-\frac{1}{x_2})-\frac{x}{x_1^2}f_x$. Therefore, the left-hand side of eq.(4) equals $f_{xx}(x-\frac{x_1}{x_2})^2+f_x(x-\frac{x_1}{x_2})-x f_x+2 f_x(x-\frac{x_1}{x_2})-x_1^2 f-2 x f_x+(1-x^2)f_{xx}=f_{xx}((x-\frac{x_1}{x_2})^2+1-x^2)+f_x(x-\frac{x_1}{x_2}-x+2 x -2\frac{x_1}{x_2}-2 x)-x_1^2 f=f_{xx}(-\frac{2 x x_1}{x_2}+\frac{x_1^2}{x_2^2}+1)+f_x(-3\frac{x_1}{x_2})-x_1^2 f=f_{xx}\frac{x_1^2+x_2^2-2 x x_1 x_2}{x_2^2}+f_x(-3\frac{x_1}{x_2})-x_1^2 f$.(5)

If we denote $z=\sqrt{x_1^2+x_2^2-2 x x_1 x_2}$, then $f=\frac{\exp(-z)}{z}$, and $f_x=f_z z_x$, $f_z=\frac{-\exp(-z)z-\exp(-z)}{z^2}=-(1+\frac{1}{z})f$, $z_x=\frac{1}{2 z}(-2 x_1 x_2)$, $f_x=-(1+\frac{1}{z})f\frac{1}{2 z}(-2 x_1 x_2)=x_1 x_2 f(\frac{1}{z}+\frac{1}{z^2})$, $f_{xx}=x_1 x_2 f_x(\frac{1}{z}+\frac{1}{z^2})+x_1 x_2 f(-\frac{1}{z^2}-\frac{2}{z^3})\frac{1}{2 z}(-2 x_1 x_2)$, $\frac{z^2 f_{xx}}{x_2^2}=\frac{x_1}{x_2}f_x(z+1)+x_1^2 f(\frac{1}{z}+\frac{2}{z^2})$. Thus, the right-hand side of eq.(5) equals $\frac{x_1}{x_2}f_x(z+1)+x_1^2 f(\frac{1}{z}+\frac{2}{z^2})+f_x(-3\frac{x_1}{x_2})-x_1^2 f=f_x(\frac{x_1}{x_2}(z+1)-3\frac{x_1}{x_2})+x_1^2 f(\frac{1}{z}+\frac{2}{z^2}-1)=x_1 x_2 f(\frac{1}{z}+\frac{1}{z^2})(\frac{x_1}{x_2}(z+1)-3\frac{x_1}{x_2})+x_1^2 f(\frac{1}{z}+\frac{2}{z^2}-1)=x_1^2 f((\frac{1}{z}+\frac{1}{z^2})(z+1-3)+\frac{1}{z}+\frac{2}{z^2}-1)=0.$ Therefore, $Q$, the left-hand side of equality (1), does satisfy eq.(2) with respect to $x_1$ and with respect to $x_2$. That means that $Q$ equals a product of a linear combination of modified spherical Bessel functions of the first and the second kind of argument $x_1$ and a linear combination of modified spherical Bessel functions of the first and the second kind of argument $x_2$. As $Q\rightarrow 0$ as $x_1\rightarrow 0$ (EDIT 11/19/2012: or $Q$ tends to a finite value if $k=0$), and $Q\rightarrow 0$ as $x_2\rightarrow \infty$ (we use here our assumption that $r_1\leq r_2$ and, therefore, $x_1\leq x_2$), we obtain that $Q=\alpha i_k(x_1)k_k(x_2)$, where $\alpha=const$. To determine $\alpha$ (and this does not look easy), let us assume that $x_1=x_2=y\rightarrow 0$. Then $Q=\frac{1}{2}\int^{1}_{-1}\frac{\exp(-y\sqrt{2-2 x})}{y\sqrt{2-2 x}}P_k(x)dx\approx \frac{1}{2 y}\int^{1}_{-1}\frac{P_k(x)dx}{\sqrt{2-2 x}}.$ As $(1-2 x t +t^2)^{-\frac{1}{2}}=\sum_{n=0}^{\infty}P_n(x)t^n$ ([5], eq. (38)), for $t=1$ we obtain $\frac{1}{\sqrt{2-2 x}}=\sum_{n=0}^{\infty}P_n(x).$ As $\int^{1}_{-1}P_n(x)P_k(x)dx=\frac{2}{2 k+1}\delta_{nk}$ ([5], eq.(28)), $Q\approx \frac{1}{2 y}\int^{1}_{-1}(\sum_{n=0}^{\infty}P_n(x))P_k(x)dx=\frac{1}{2 y}\frac{2}{2 k+1}$. Using the relations above between modified spherical Bessel functions and modified Bessel functions, the asymptotic expressions for modified Bessel functions for small arguments ([6]): $I_{\beta}(x)\approx\frac{1}{\Gamma(\beta+1)}(\frac{y}{2})^{\beta}$, $K_{\beta}(y)\approx\frac{\Gamma(\beta)}{2}(\frac{2}{y})^{\beta}$ ($\beta>0$), and the property of the $\Gamma$-function: $\Gamma(z+1)=z \Gamma(z)$, we obtain that $\alpha=1$, and thus the equality in the question is completely proven.

[1] http://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheFirstKind.html

[2] http://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html

[3] http://dlmf.nist.gov/10.47

[4] http://en.wikipedia.org/wiki/Legendre_polynomials

[5] http://mathworld.wolfram.com/LegendrePolynomial.html

[6] http://en.wikipedia.org/wiki/Bessel_function

END EDIT

I don't have anything close to a complete proof, just some considerations, for what it's worth. If the formula is correct and, say, $r_1>r_2$, then the integral $I(r_1,r_2)$ on the left-hand-side equals a product of a function of $r_1$ and a function of $r_2$, so the mixed derivative $\frac{\partial \ln(I)}{\partial r_1 \partial r_2}$ should vanish. It may be advisable to try to prove that. Taking derivatives of $I$ yields the integrand at the integration limits, where it is easy to calculate the Legendre polynomial (I guess it's either +1 or -1) and an integral of the derivative of the integrand, where recursive relations for the Legendre polynomials might be useful. If the statement above is proven, it would be enough to prove the statement in the question for some specific value of, say, $r_2$, such as zero (in that case limits should be taken) or 1. Then the first formula of http://en.wikipedia.org/wiki/Plane_wave_expansion might be useful.