# Expansion of $\frac{\text{erfc}({|\vec{r} - \vec{r'}|})}{|\vec{r} - \vec{r'}|}$ in spherical harmonics?

How can I derive the spherical harmonic expansion coefficients for the function $$\frac{\text{erfc}({|\mathbf{r} - \mathbf{r'}|})}{{|\mathbf{r} - \mathbf{r'}|}}$$ by expressing it as $$f(\theta, \phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} c_{lm} Y_{lm}(\theta, \phi)$$?

The coefficients are determined through $$c_{lm} = \int\int_{\text{sphere}} \frac{\text{erfc}({|\mathbf{r} - \mathbf{r'}|})}{{|\mathbf{r} - \mathbf{r'}|}} \cdot Y_{lm}^*(\theta, \phi) , \sin(\theta) , d\theta , d\phi$$

I am aware that the expression $$\frac{1}{{r - r'}}$$ can be expanded in spherical harmonics using the following formula:

$$\frac{1}{{\lvert \mathbf{r} - \mathbf{r}' \rvert}} = 4\pi \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{1}{{2l+1}} \frac{r_<^l}{{r_>^{l+1}}} Y_{lm}^*(\hat{\mathbf{r}}') Y_{lm}(\hat{\mathbf{r}})$$

Is it possible to perform a similar expansion for $$\frac{\text{erfc}({|\mathbf{r} - \mathbf{r'}|})}{{|\mathbf{r} - \mathbf{r'}|}}$$, where the numerator involves the complementary error function?

What i have done so far was, using the expansion of inverse distance between two vectors arises often in physics in Legendre polynomials:

$$\frac{1}{|\mathbf{r} - \mathbf{r'}|} = \frac{1}{r} \left(1 + \frac{{r'}^2}{r^2} - 2\frac{r'}{r} \cos \gamma\right)^{-1/2} = \frac{1}{r} \sum_{\ell=0}^{\infty} \left(\frac{r'}{r}\right)^\ell P_\ell(\cos \gamma)$$

I would actually assume that starting point of the expansion is something like:

$$\frac{\text{erfc}({\lvert \mathbf{r} - \mathbf{r}' \rvert})}{{\lvert \mathbf{r} - \mathbf{r}' \rvert}} = \text{erfc}({\lvert \mathbf{r} - \mathbf{r}' \rvert}) \, 4\pi \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{1}{{2l+1}} \frac{r_<^l}{{r_>^{l+1}}} Y_{lm}^*(\hat{\mathbf{r}}') Y_{lm}(\hat{\mathbf{r}})$$