# 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

Using this formula for the spherical harmonic function, and making a few simplifications, here is a direct implementation of the solid spherical harmonic function:

dpower[x_, y_] := Piecewise[{{1, y == 0}}, x^y]

SolidHarmonicS[λ_Integer?NonNegative, μ_Integer, x_, y_, z_] /; Abs[μ] <= λ :=
With[{s = Sign[μ], am = Abs[μ]},
(-1)^((1 - s) am/2) Sqrt[((2 λ + 1) (λ - am)!)/(4 π (λ + am)!)]
dpower[x + I s y, am] Sum[(-1)^((λ + am)/2 - k) (λ + am + 2 k - 1)!!
dpower[z, 2 k] dpower[x^2 + y^2 + z^2, (λ - am)/2 - k]/
((2 k)! (λ - am - 2 k)!!),
{k, Mod[λ - am, 2]/2, (λ - am)/2}]]


I chose to use the direct sum instead of using the hypergeometric representation to avoid having to do a polynomial division, which can be troublesome for zero arguments.

A test:

Table[SolidHarmonicS[λ, μ, x, y, z], {λ, 0, 2}, {μ, -λ, λ}] // Simplify
{{1/(2 Sqrt[π])},
{1/2 Sqrt[3/(2 π)] (x - I y), 1/2 Sqrt[3/π] z, -(1/2) Sqrt[3/(2 π)] (x + I y)},
{1/4 Sqrt[15/(2 π)] (x - I y)^2, 1/2 Sqrt[15/(2 π)] (x - I y) z,
-(1/4) Sqrt[5/π] (x^2 + y^2 - 2 z^2), -(1/2) Sqrt[15/(2 π)] (x + I y) z,
1/4 Sqrt[15/(2 π)] (x + I y)^2}}

• Oh, of course. This would be a bad idea at large $\lambda,\mu$ (and it's probably not how SphericalHarmonicY is implemented) but it should work fine for the applications I envisaged. (And, if it breaks numerically at large parameters for whoever is pushing it there, then that's well deserved and they should seriously think about coming back to the $Y_{lm}$ fold.) – Emilio Pisanty Aug 19 '16 at 17:49
• For Legendre functions of large degree and/or order, I would probably use an asymptotic expansion instead. – J. M. is away Aug 19 '16 at 19:29
• @Emilio, it's okay; I appreciate that you're reporting this. Now if only I can find a bit of spare time to dig deeper... – J. M. is away Dec 16 '16 at 2:11
• Hi! another minor bug: it works a bit wonkily at $\mathbf r=0$. This comes from the factor of z^(λ - Abs[μ] - 2 k), which shouldn't count because the Pochhammer symbol is zero when λ - Abs[μ] - 2 k is negative. – Emilio Pisanty Jan 26 '17 at 18:45
• @J.M. Oh, for sure, no pressure =). I'm just writing them down here instead of waiting to report bug #2 until after you've fixed bug #1 ;-). I've implemented a version here but I don't like the maze of Ifs it's become, and I didn't want to jump in and impose it on your answer. – Emilio Pisanty Jan 26 '17 at 18:51

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.