In the implementation of multipole_psi, the radial prefactor of the spherical multipole expansion seems to fall off as $\frac{1}{r}$.
Starting from the standard spherical multipole expansion
$\phi(\mathbf{r})=\sum_{\ell,m}\frac{4\pi}{2\ell+1}\frac{q_{\ell,m}}{r^{2\ell+1}}r^{\ell}Y_{\ell,m}$,
where $q_{\ell,m}$ are the corresponding multipole moments, one can recognize $Y_{\ell,m}$ as Laplace's spherical harmonics, do some math, and in the end should end up at the definition of the multipole moments that go into $\Psi$ (computed in multipole_psi).
When I do the derivation, however, I arrive at a slightly different expression than what is implemented in multipole_psi. Is it possible that instead of having a prefactor of $\frac{4\pi}{(2\ell+1)\cdot r^\ell}$, it needs to be $\frac{4\pi}{(2\ell+1)\cdot r^{(\ell+1)}}$?
The $\frac{1}{r^{\ell+1}}$ is just the simplification of $\frac{r^\ell}{r^{2\ell+1}}$, isn't it? I am wondering where in the derivation one loses the '+1' in the exponent of $r$. If the code currently uses $\frac{1}{r^\ell}$, the monopole term would be constant instead of scaling with $\frac{1}{r}$, which seems inconsistent.
A clarification would be much appreciated!
