EDITAR (11/18/2012) Como yo no pudo obtener una prueba usando mi razonamiento inicial (ver más abajo, al final), me permiten ofrecer una completa prueba utilizando un enfoque diferente.
Transformemos la integral en el lado izquierdo primero, utilizando la siguiente sustitución: $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$ donde $x_1=\frac{r_1}{d},x_2=\frac{r_2}{d}.$
Transformemos el lado derecho ahora (es de suponer, sin pérdida de generalidad, que $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)$, donde $i_k(x)=\sqrt{\frac{\pi}{2 x}}I_{k+\frac{1}{2}}(x)$ $k_k(x)=\sqrt{\frac{2}{\pi x}}K_{k+\frac{1}{2}}(x)$ son modificados esférica funciones de Bessel de primera y segunda clase, respectivamente ([1],[2]). Por lo tanto, la igualdad en la pregunta es equivalente a la siguiente:
$\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)
Es bien sabido que modificó el esférico funciones de Bessel de primera y segunda especie $i_k(z)$ $k_k(z)$ satisfacer la siguiente ecuación diferencial ([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)
Si la igualdad (1) es correcta, su mano izquierda, nos vamos a denotar es $Q$) deben satisfacer la ecuación.(2) con respecto a $x_1$ y con respecto a $x_2$. Vamos a demostrar que satisface la ecuación.(2) con respecto a $x_1$ (la prueba de $x_2$ es idéntico), es decir, que
$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)
donde
$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}}.$
Como $\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$,
donde la diferenciación por partes se ha utilizado dos veces, y el subíndice $x$ significa que la diferenciación con respecto a $x$. Por lo tanto, para demostrar eq.(3), es suficiente para demostrar que
$\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)
Los derivados de la $f$ con respecto al $x$ $x_1$ (se denota ellos $f_x$$f_1$, respectivamente) puede ser obtenido usando la regla de la cadena, y la comparación de los resultados, obtenemos $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})$. La segunda derivada de $f$ con respecto al $x_1$ es igual a $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$. Por lo tanto, el lado izquierdo de la ecuación.(4) es igual a $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)
Si denotamos $z=\sqrt{x_1^2+x_2^2-2 x x_1 x_2}$,$f=\frac{\exp(-z)}{z}$, y $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})$. Por lo tanto, el lado derecho de la ecuación.(5) es igual a $\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.$ por lo Tanto, $Q$, el lado izquierdo de la igualdad (1), satisface la ecuación.(2) con respecto a $x_1$ y con respecto a $x_2$. Eso significa que $Q$ es igual a un producto de una combinación lineal de modificación esférica funciones de Bessel de primera y el segundo tipo de argumento $x_1$ y una combinación lineal de modificación esférica funciones de Bessel de primera y el segundo tipo de argumento $x_2$. Como $Q\rightarrow 0$ $x_1\rightarrow 0$ (EDICIÓN 11/19/2012: o $Q$ tiende a un valor finito si $k=0$), y $Q\rightarrow 0$ $x_2\rightarrow \infty$ (utilizamos aquí nuestra suposición de que $r_1\leq r_2$ y, por lo tanto, $x_1\leq x_2$), obtenemos que $Q=\alpha i_k(x_1)k_k(x_2)$ donde $\alpha=const$. Para determinar el $\alpha$ (y esto no parece fácil), supongamos que $x_1=x_2=y\rightarrow 0$. A continuación, $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}}.$ $(1-2 x t +t^2)^{-\frac{1}{2}}=\sum_{n=0}^{\infty}P_n(x)t^n$ ([5], eq. (38)), para $t=1$ obtenemos $\frac{1}{\sqrt{2-2 x}}=\sum_{n=0}^{\infty}P_n(x).$ $\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}$. El uso de las relaciones anteriores entre modificado esféricas de Bessel y funciones modificadas de Bessel funciones, el asintótica expresiones para las funciones de Bessel modificadas para los pequeños argumentos ([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$), y la propiedad de la $\Gamma$-función: $\Gamma(z+1)=z \Gamma(z)$, obtenemos que $\alpha=1$, y por lo tanto la igualdad en la pregunta está totalmente comprobado.
[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
FINAL DE EDICIÓN
Yo no tengo nada cerca a una prueba, solo algunas consideraciones, por lo que vale. Si la fórmula es correcta y, digamos, $r_1>r_2$, entonces la integral de la $I(r_1,r_2)$ sobre el lado izquierdo es igual a un producto de una función de $r_1$ y una función de $r_2$, por lo que la mezcla de derivados $\frac{\partial \ln(I)}{\partial r_1 \partial r_2}$ debe desaparecer. Puede ser aconsejable para tratar de demostrar que. Tomando derivados de $I$ rendimientos el integrando en la integración límites, donde es fácil calcular el polinomio de Legendre (supongo que es +1 o -1) y la integral de la derivada de el integrando, donde las relaciones recursivas para los polinomios de Legendre podría ser útil. Si la afirmación anterior es probado, sería suficiente para demostrar el enunciado de la pregunta para algún valor específico de, digamos, $r_2$, tal como cero (en ese caso deberían ser los límites de toma) o 1. A continuación, la primera fórmula de http://en.wikipedia.org/wiki/Plane_wave_expansion podría ser útil.