The issue you have seems very similar to mine, which is answered in https://scicomp.stackexchange.com/questions/45262/m2m-in-python-using-scipy-sph-harm-y

The TLDR:
You must include a factor $ \sqrt{\frac{4 \pi}{2n +1}}$ in ALL the function (P2M, M2M and potential evaluation) to keep the normalization consistent.