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.