11

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.

Emilio Pisanty
  • 10,255
  • 1
  • 36
  • 69
  • 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
  • 1
    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
  • 1
    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 00:07

3 Answers3

12

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

mmal
  • 3,508
  • 2
  • 18
  • 38
chuy
  • 11,205
  • 28
  • 48
10

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}}
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
  • 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.'s missing motivation Aug 19 '16 at 19:29
  • @Jens, to help me debug, can you give a specific $\lambda,\mu$ with missing $z$ factors? I tried my best to simplify the formula from the Wolfram Functions site, but it is possible I may have accidentally lobbed off something important in the process. – J. M.'s missing motivation Aug 21 '16 at 03:21
  • @J.M. Sorry, it was just that the definition didn't get copied correctly! Everything is fine. I had to delete a line break that crept in before the third line of code. I had already upvoted your answer because I assumed it was some trivial omission... – Jens Aug 21 '16 at 04:49
  • Hi J.M., I just realized that this wasn't working for negative $\mu$, and I took the liberty to edit in the fix =). – Emilio Pisanty Dec 10 '16 at 00:04
  • @J.M. Another real-world bug - asking for e.g. SolidHarmonicS[1,1,1,0,0] will call 0^0 and therefore produce an error. With this one I'm less sure of the best way to fix it, since anything that redefines Power[0,0] can be dangerous in general-use code like this one. – Emilio Pisanty Dec 14 '16 at 14:48
  • @Emilio, I'll think about how to fix the thing later; I'm juggling other stuff at the moment. – J. M.'s missing motivation Dec 14 '16 at 15:24
  • @J.M. Sure, no worries. I implemented a fix on my machine, but I think it's best for you to decide how to edit the version on this answer - and on your own time and terms, of course ;-). – Emilio Pisanty Dec 14 '16 at 16:59
  • @J.M. Another small one - there is a sign discrepancy for odd, negative $m$ with respect to the standard normalization. Sorry to keep bugging - it's because I think it's the best solution here and I'd like to use this answer as the go-to place for this code ;-). I'm also sorry for the delay since you answered, this is the sort of thing that only comes out in the wash when you actually use things in anger. – Emilio Pisanty Dec 15 '16 at 20:34
  • 1
    @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.'s missing motivation Dec 16 '16 at 02:11
  • 1
    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
  • I promise @Emilio, when I get a big chunk of free time, I'll work on this... – J. M.'s missing motivation Jan 26 '17 at 18:47
  • 1
    @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
  • Quick update: there is now an implementation in the Wolfram Function Repository. – Emilio Pisanty Feb 16 '24 at 11:07
1

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.

Emilio Pisanty
  • 10,255
  • 1
  • 36
  • 69