How to find expansion of slightly modified Coulomb's potential? From here I know that,
${{\frac {1}{|\mathbf {r}_1 -\mathbf {r}_2|}}=\sum _{\ell =0}^{\infty }{\frac {4\pi }{2\ell +1}}\sum _{m=-\ell }^{\ell }(-1)^{m}{\frac {r_1^{\ell }}{r_2^{\ell +1}}}Y_{\ell }^{-m}(\theta ,\varphi )Y_{\ell }^{m}(\theta ',\varphi ').}$
where $|r_1|<|r_2|.$ I was wondering if there exists a strategy to find expansions for general potentials, for instance
$$\frac{1}{|r_1-r_2|^{1+\alpha}}, \text{    }\alpha>0$$
or more generally for functions $\phi(|r_1-r_2|)?$
 A: The usual multipole expansion follows from the Legendre identity
$$
{\displaystyle {\frac {1}{\sqrt {1-2xt+t^{2}}}}=\sum _{n=0}^{\infty }P_{n}(x)t^{n}}
$$
The generalization to arbitrary powers requires the Gegenbauer identity
$$
{\frac {1}{(1-2xt+t^{2})^{\alpha }}}=\sum _{{n=0}}^{\infty }C_{n}^{{(\alpha )}}(x)t^{n}
$$
which yields
$$
{\displaystyle {\frac {1}{|\mathbf {x} -\mathbf {y} |^\alpha}}=\sum _{k=0}^{\infty }{\frac {|\mathbf {x} |^{k}}{|\mathbf {y} |^{k+\alpha}}}C_{k}^{(\alpha/2)}\biggl({\frac {\mathbf {x} \cdot \mathbf {y} }{|\mathbf {x} ||\mathbf {y} |}}\biggr)}
$$
You can write this in terms of spherical harmonics, if you wish.
For arbitrary potentials, you can use the fact that the spherical harmonics are a basis, and therefore you can expand
$$
\phi(\vec r_1,\vec r_2)=\sum_{\ell,m\\\ell',m'}c_{\ell,m;\ell',m'}(r_1,r_2) Y_\ell^m(\theta,\phi)Y_{\ell'}^{m'}(\theta',\phi')
$$
where the coefficients are
$$
c_{\ell,m;\ell',m'}(r_1,r_2)=\int \phi(\vec r_1,\vec r_2)(Y_\ell^m(\theta,\phi)Y_{\ell'}^{m'}(\theta',\phi'))^*
$$
Unless you carefully choose $\phi$ to be a special function (e.g., with some interesting symmetries), this expression cannot in general be simplified much further.
