Are solid spherical harmonics implemented in Mathematica?

Question

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?

---

Edit: Since January 2021, the Wolfram Function Repository contains an implementation, which it calls SolidHarmonicR.

Answer 1

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}}

Answer 2

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}$$

Answer 3

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.