Infinite integrals in calculation of dipole potential The Situation
Assume that a point dipole is positioned at $z=D$ and is pointing in an arbitrary direction so that its dipole moment is $\textbf{p} = p \cos \alpha \, \mathbf{\hat{r}} + p \sin \alpha \, \mathbf{\hat{\theta}}$ ($\alpha$ is the angle that the dipole moment vector makes with the $z$ axis, and we choose the $x$ axis so that the dipole moment has 0 component in the $y$ direction).
Now, I want to expand the potential of this dipole in whole space in terms of spherical harmonics (so that I may later solve a boundary value problem). I noticed that I may do so by using the multipole expansion and the fact that the charge density of a point dipole is equal to:
$$ \rho(\textbf{x}) = - \textbf{p} \cdot \nabla \delta(\textbf{x}-\textbf{x'}) $$
The Problem
In my case this should evaluate to:
$$ \rho(\textbf{x}) = -\frac{p \cos\alpha}{r^3} \left( \delta'(r-D)r - 2\delta(r-D) \right) \delta(\cos\theta-1) \delta(\varphi) + \frac{p\sin\alpha}{r^3} \sin\theta \delta(r-D)\delta'(\cos\theta - 1)\delta(\varphi) $$
The multipole expansion in spherical harmonics is (I'm interested only in potential at points where $r<D$):
$$ \phi(\textbf{x}) = \frac{1}{\varepsilon_0} \sum_{l=0}^{\infty} \sum_{m=-l}^l \frac{r^l}{2l+1} Y_{lm}(\theta, \varphi) \int \frac{\rho(\textbf{x'})}{(r')^{l+1}} \text{d}^3x'$$
Now that I calculate the integrals, I encounter a problem with the integral of the term with $\delta'(\cos\theta'-1)$, namely I have the integral which has the form:
$$ \int_0^{\pi} \text{d}\theta' P_l^m(\cos\theta')\sin^2(\theta')\delta'(\cos\theta'-1) $$
I tried to deal with it in the same manner as with the others, namely by substitution ($x=\cos\theta'$). However, as I do it, I get the following:
$$ \int_{-2}^0 \text{d}x \, P_l^m(x+1) \sqrt{-x^2-2x} \, \delta'(x) $$
which seems to be divergent. I also tried using the formula $x\delta'(x)=-\delta(x)$, or writing down the delta in terms of $\theta$ alone from the beginning. Unfortunately, I end up with a seemingly divergent integral either way. However, the potential of a point dipole should be finite, since the well-known formula for it is:
$$ \phi = \frac{1}{4\pi\varepsilon_0} \frac{\mathbf{p} \cdot \mathbf{x}}{x^3} $$
Can anyone tell me what am I doing wrong in this case? I tried to search the Internet for any clues but I have found nothing and am absolutely lost.
 A: Instead of projecting on the spherical harmonics, a faster way is to convert the formula:
$$
\rho=-p\cdot\nabla \delta
$$
in spherical coordinates by derivating with respect to the position of the monopole.
Let’s first do the spherical harmonic expansion of a monopole. When it is at coordinate $(r_0,\theta_0,\phi_0)$, then:
$$
\rho(r,\theta,\phi)= \sum_{l\in\mathbb N,|m|\leq l}\frac{\delta(r-r_0)}{r_0^2}Y_l^m(\theta_0,\phi_0)Y_l^m(\theta,\phi)
$$
Using (changed the potential’s name to avoid ambiguity)
$$
V(r,\theta,\phi)=\sum V^m_l(r)Y_l^m(\theta,\phi)
$$
You have:
$$
-{V^m_l}’’(r)-\frac{2}{r}{V^m_l}’(r)+\frac{l(l+1)}{r^2}V^m_l(r) = \frac{\delta(r-r_0)}{r_0^2}Y_l^m(\theta_0,\phi_0)
$$
which you solve with piecewise power laws:
$$
V_l^m(r)=\begin{cases}
\frac{1}{(2l+1)r_0}\left(\frac{r}{r_0}\right)^l Y_l^m(\theta_0,\phi_0) & r\leq r_0 \\
\frac{1}{(2l+1)r_0}\left(\frac{r}{r_0}\right)^{-l-1} Y_l^m(\theta_0,\phi_0) & r\geq r_0
\end{cases}
$$
Writing $p=p_re_r+p_\theta e_\theta+p_\phi e_\phi$ with the basis vectors $e_r,e_\theta,e_\phi$ taken at position $r_0,\theta_0,\phi_0)$, you get in general:
$$
-p\cdot \nabla f(\vec x-\vec x_0)= p_r\partial_{r_0}f+\frac{p_\theta}{r_0}\partial_{\theta_0}f+\frac{p_\phi}{r_0\sin\theta_0}\partial_{\phi_0}f
$$
Your new charge density:
$$
\rho= \sum \left[-p_r\left(\frac{\delta’(r-r_0)}{r_0^2}+2\frac{\delta(r-r_0)}{r_0^3} \right)Y_l^m(\theta_0,\phi_0)+\frac{\delta(r-r_0)}{r_0^2}\left(\frac{p_\theta}{r_0}\partial_{\theta_0}Y_l^m(\theta_0,\phi_0)+\frac{p_\phi}{r_0\sin\theta_0}\partial_{\phi_0}Y_l^m(\theta_0,\phi_0)\right)\right]Y_l^m
$$
and similarly you can apply the same trick to the potential:
$$
V= \sum \left[p_r \partial_{r_0}A_l^m(r,r_0) Y_l^m(\theta_0,\phi_0)+ A_l^m(r,r_0)\left(\frac{p_\theta}{r_0}\partial_{\theta_0}Y_l^m(\theta_0,\phi_0)+\frac{p_\phi}{r_0\sin\theta_0}\partial_{\phi_0}Y_l^m(\theta_0,\phi_0)\right)\right]Y_l^m
$$
with:
$$
A_l^m(r,r_0) = \begin{cases}
\frac{1}{(2l+1)r_0}\left(\frac{r}{r_0}\right)^l & r\leq r_0 \\
\frac{1}{(2l+1)r_0}\left(\frac{r}{r_0}\right)^{-l-1} & r\geq r_0
\end{cases} \\
\partial_{r_0}A_l^m(r,r_0)
= \begin{cases}
-\frac{l+1}{(2l+1)r_0^2}\left(\frac{r}{r_0}\right)^l & r\leq r_0 \\
\frac{l}{(2l+1)r_0^2}\left(\frac{r}{r_0}\right)^{-l-1} & r\geq r_0
\end{cases}
$$
Hope this helps.
