I have been trying to use the 3D fast multipole expansion,
$$ \Phi(P) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \frac{M_n^m}{r^{n+1}}Y_n^m (\theta, \phi) $$
to approximate the source potential in the far field,
$$ \Phi(P) = \frac{q}{r} $$
, where $ q $ is the source strength, while $ M_n^m $ are constants associated with the location and strength $( q )$ of the source particle only, and are independent of $ (\theta, \phi) $; $ Y_m^n $ is the spherical harmonics,
$$ Y_n^m = \sqrt{\frac{(n - |m|)!}{(n + |m|)!}} . P_n^{|m|}(\cos{\theta}) \exp(i m \phi) $$
, where $ P_n^{|m|} $ is the associated Lengendre function,
$$ P_n^m = \frac{(-1)^m}{2^n n!} (1 - x^2)^{\frac{m}{2}} \frac{\partial ^{m+1}}{ \partial x ^{m+n}} (x^2 - 1)^n $$
My goal is to calculate the spatial gradient, $ \nabla\Phi $, using the multipole expansion, which results respectively in the $ \bf{e} _{\theta}, \bf{e} _{\phi} $ directions,
$$ e _{\theta} \cdot \nabla \Phi = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \frac{M_n^m}{r^{n+2}} \sqrt{\frac{(n - |m|)!}{(n + |m|)!}} \frac{\partial P_n^{|m|}(\cos{\theta})}{\partial \theta} (- \sin \theta) \exp(i m \phi) $$
and
$$ e _{\phi} \cdot \nabla \Phi = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \frac{i m M_n^m}{r^{n+2} \sin \theta} \sqrt{\frac{(n - |m|)!}{(n + |m|)!}} P_n^{|m|}(\cos{\theta}) \exp(i m \phi) $$
The derivatives $ \frac{\partial P_n^{|m|}(\cos{\theta})}{\partial \theta} $ goes to infinity at $ \theta = \pm \pi $, which has been noted in here. As well, in $ e _{\phi} \cdot \nabla \Phi $ the $ \sin \theta $ in the denominator may cause an infinity problem at $ \theta = 0 $ and $ \pi $. So I'm wondering if there is any trick to avoid this happen in the numerical implementation?
Appreciated.