I am working with
$Y_{l,m}(\theta-\theta', \phi -\phi')$
and I believe there is a nice way to write that as a product of Spherical Harmonics, but I cannot derive it or find it anywhere. Is it possible to write a Spherical Harmonics of a sum/difference as the product of Spherical Harmonics in a simple form?
EDIT:
Specifcally, I am interested in the way that
$$ \int_{-\infty}^{\infty} (\sum\limits_{l=0}^{\infty} \sum\limits_{m=-l}^{+l} f_{lm}(r-r')Y_{lm}(\theta-\theta',\phi-\phi')) (\sum\limits_{l=0}^{\infty} \sum\limits_{m=-l}^{+l} g_{lm}(r') Y_{lm}(\theta',\phi'))) dV' $$
would behave. In order to simplify this, it is necessary to break up the $Y_{lm}$ in the first expansion into part which are only evaluated at $(\theta', \phi')$ in order to rely on the orthonormality of Spherical Harmonics.
To answer your first question: there is a paper by Rami Mehrem http://arxiv.org/abs/0909.0494 that gives this in equation 5.12. How you get it: write the Y_lms in terms of solid harmonics R_lm, and then use the solid harmonic addition theorem. http://en.wikipedia.org/wiki/Solid_harmonics#Addition_theorems
An unfortunate aspect of this is that I think it does mix the angular and radial dependences so the factorization you have above won't be clean anymore.
I think a better way to do this problem is to see that if you pull all the sums outside the integral, you are trying to convolve two functions with spherical harmonic angular dependence. If you move to Fourier space, this becomes a product by the Convolution Theorem. I wrote a blog post on how to do this problem:
http://humblecosmologist.blogspot.com/2015/03/convolution-of-two-functions-with.html
Note that one does need some conditions on f_lm and g_lm in your notation above, or K_lm and Q_l'm' in mine in the blog post, to guarantee that the radial integrals converge.