# Are solid spherical harmonics implemented in Mathematica?

In certain applications, solid spherical harmonics can be very useful. They are essentially the usual, 'surface' spherical harmonics, with the appropriate power of the radius inserted:

$$S_l^m(x,y,z)=r^l Y_l^m(\theta,\phi).$$

They are particularly useful because, for integer $l$ and $m$, they are homogeneous polynomials of degree $l$ in $x$, $y$ and $z$, and do not therefore require any ugly square roots or inverse trigonometric functions. The first few such functions are as follows: $$\begin{array}{ccc} \frac{1}{2 \sqrt{\pi }} & 0 & 0 \\ \frac{1}{2} \sqrt{\frac{3}{\pi }} z & -\frac{1}{2} \sqrt{\frac{3}{2 \pi }} (x+i y) & 0 \\ -\frac{1}{4} \sqrt{\frac{5}{\pi }} \left(x^2+y^2-2 z^2\right) & -\frac{1}{2} \sqrt{\frac{15}{2 \pi }} z (x+i y) & \frac{1}{4} \sqrt{\frac{15}{2 \pi }} (x+i y)^2 \\ \end{array}$$

More generally, though, it is useful to be able to implement a spherical harmonic when the natural input is the cartesian components of the argument.

Is this implemented in Mathematica? There is nothing very obvious at all in the documentation, though maybe it is hidden away in some third-party package. If there isn't a built-in function, is there some specific reason for that?

-
In addition to the answers below, see my answer on PhysicsSE about computing irregular real solid harmonic expansions of a current loop. – DumpsterDoofus Apr 11 '14 at 23:47
I chose to use the real spherical harmonics SphericalHarmonicYr or $Y_{L,m}$ rather than the standard built-in quantum spherical harmonics SphericalHarmonicY or $Y_L^m$ because they're useful for data which is purely real-valued. As an example of converting to Cartesian coordinates, one can do Table[FullSimplify[ SolidHarmonicRr[L, m, ##] & @@ CoordinateTransform["Cartesian" -> "Spherical", {x, y, z}], Assumptions -> Element[{x, y, z}, Reals]], {L, 0, 2}, {m, -L, L}]. – DumpsterDoofus Apr 11 '14 at 23:58
This yields $$\begin{array}{ccccc} 1 & \text{} & \text{} & \text{} & \text{} \\ y & z & x & \text{} & \text{} \\ \sqrt{3} x y & \sqrt{3} y z & -\frac{x^2}{2}-\frac{y^2}{2}+z^2 & \sqrt{3} x z & \frac{1}{2} \sqrt{3} (x-y) (x+y) \\ \end{array}$$ which are the real regular solid harmonics up to order $L=2$. – DumpsterDoofus Apr 11 '14 at 23:59
I'm not sure, but I suspect that the Cartesian definitions might become numerically useless for high enough values of $L$ because the polynomial coefficients become gigantic enough that catastrophic cancellation can occur. A similar problem happens with ChebyshevT when you first compute ChebyshevT[60,x] and then do %/.x->2.0/3 which gives an incorrect answer. So if you're doing numerics for high $L$ you may actually want to evaluate them in spherical coordinates (but again I haven't tested to see if this is necessary). – DumpsterDoofus Apr 12 '14 at 0:07

You could also do something like the following (taking advantage of TransformedField:

solidHarmonicS[l_?IntegerQ, m_?IntegerQ, x_, y_, z_] :=
Module[{r, θ, ϕ, xx, yy, zz},
FullSimplify@
Evaluate[
TransformedField["Spherical" -> "Cartesian",
r^l SphericalHarmonicY[l, m, θ, ϕ],
{r, θ, ϕ} -> {xx, yy, zz}]] /. {xx -> x, yy -> y, zz -> z}
]


$$\begin{array}{ccc} \frac{1}{2 \sqrt{\pi }} & 0 & 0 \\ \frac{1}{2} \sqrt{\frac{3}{\pi }} z & -\frac{1}{2} \sqrt{\frac{3}{2 \pi }} (x+i y) & 0 \\ -\frac{1}{4} \sqrt{\frac{5}{\pi }} \left(x^2+y^2-2 z^2\right) & -\frac{1}{2} \sqrt{\frac{15}{2 \pi }} z (x+i y) & \frac{1}{4} \sqrt{\frac{15}{2 \pi }} (x+i y)^2 \\ \end{array}$$

-

This can be implemented by appropriate replacement rules, though they are a bit fiddly to get just right (note in particular the exact form of the first one). The implementation below uses partial memoization as in this question, and it will only match the case when $l$ and $m$ are integers; otherwise, it will go through a direct SphericalHarmonicY evaluation.

SolidHarmonicS[l_, m_, x_, y_, z_] := If[
IntegerQ[l] && IntegerQ[m],
Block[{x1, y1, z1},
SolidHarmonicS[l, m, x1_, y1_, z1_] = FullSimplify[ Module[{θ, ϕ},
(x1^2 + y1^2 + z1^2)^(l/2)
SphericalHarmonicY[l, m, θ, ϕ] /. {
Power[E,Times[Complex[0, μ_], ϕ]] -> (x1 + I y1)^μ/(x1^2 + y1^2)^(μ/2),
Sin[θ] -> (x1^2 + y1^2)^(1/2)/(x1^2 + y1^2 + z1^2)^(1/2),
Sin[θ]^μ_ -> (x1^2 + y1^2)^(μ/2)/(x1^2 + y1^2 + z1^2)^(μ/2),
Cos[θ] -> z1/(x1^2 + y1^2 + z1^2)^(1/2),
Cos[θ]^λ_ -> z1^λ/(x1^2 + y1^2 + z1^2)^(λ/2)
}
]];
SolidHarmonicS[l, m, x, y, z]
],
(x^2 + y^2 + z^2)^(l/2)
SphericalHarmonicY[l, m, ArcCos[z/(x^2 + y^2 + z^2)^(l/2)],Arg[x + I y]]
]


For some sample output, see the table in the question, which is a direct copy of the output of

TeXForm[Grid[Table[
SolidHarmonicS[l, m, x, y, z]
, {l, 0, 2}, {m, 0, 2}]]]


I would of course be interested in seeing other implementations.

-