I detected that spherical_harmonic function switches phase of functions with odd values of m. Help file reference towards Merzbacher 9.64 (in fact this is 9.65), but results from Maxima doesn't correspond to sample formulas in 9.68 (and http://mathworld.wolfram.com/SphericalHarmonic.html\) - see Y1, +-1 and Y2, +- 1. Maxima gives formula with opposite phase (sign). this is probably due to use factor (-1)^m both in spherical_function and assoc_legendre_p (which gives proper results according to http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html\). Formula 9.64 in Merzbacher includes (-1)^m factor, but it is ommited in 9.59 - definition of associated Legendre polynomials.
The problem is described in http://mathworld.wolfram.com/Condon-ShortleyPhase.html. I suggest omitting (-1)^m factor in definition of spherical_harmonic and/or some description in help file for spherical_harmonic. Temporary formula
Y(l, m, theta, phi):=(-1)^m*spherical_harmonic(l, m, theta, phi);
works fine.
Thank you for the detailed analysis. I proposed using the definition in dlmf;
specifically http://dlmf.nist.gov/14.30#E3.