Generating harmonic polynomials in cartesian coordinates

Question

TLDR: Are these polynomials really harmonic polynomials, and how can I generate them?

---

Long version:

I want to describe an electrostatic potential $\Phi(x,y,z)$ over a source-free volume, by using a polynomial of the form

$$\begin{align} \Phi(x,y,z) &= a_{000} + a_{100} x + a_{010} y + \dots + a_{123}xy^2z^3 + \dots \\ &= \sum_{l,m,n \geq 0}^\infty a_{lmn}x^ly^mz^n \quad . \end{align} $$

The potential is known to me on a grid: I simulated it using a finite element simulation, and I want to find the coefficients $a_{lmn}$ so that I can give a smooth estimate of the potential around the origin. I verified that my potential $\Phi(x,y,z)$ is "benign" in the sense that I can neglect all terms higher than some order $N$. I can somewhat successfully find the coefficients $a_{lmn}$ with a crude fit, but of course the fit-functions $x^ly^mz^n$are not orthogonal in any way, which makes "finding the right order $N$" more important than in should be.

But since the volume that I am interested in is charge free, I know that the total polynomial must satisfy the Laplace equation $\Delta\Phi(x,y,z) = 0$. I want to use this fact to rewrite the polynomial as a sum of smaller polynomials that, on their own, satisfy the Laplace equation and are orthogonal to each other. I believe such polynomials are called harmonic polynomials, and the whole process reminds me of doing multipole expansions, only that instead of sitting far away from the charges that create the potential, I sit between them.

Using horror-inducing Mathematica expressions that transmogrify spherical harmonics to real-valued, cartesian harmonics, I found a set of polynomials $p_{n,m}(x,y,z)$ of degree $n$, with $-n \leq m \leq n$, which as far as I have tested,

• satisfy the Laplace equation,
• are orthogonal over the cube $x,y,z \in [-1;1]$
• are normalizied over the cube $x,y,z \in [-1;1]$

The first of these polynomials are (dropping all normalization factors to make them shorter)

$$\begin{align} p_{0,0}(x,y,z) &= 1\\ \\ p_{1,-1}(x,y,z) &= x\\ p_{1,0}(x,y,z) &= y\\ p_{1,1}(x,y,z) &= z\\ \\ p_{2,-2}(x,y,z) &= xy\\ p_{2,-1}(x,y,z) &= yz\\ p_{2,0}(x,y,z) &= x^2+y^2-2 z^2\\ p_{2,1}(x,y,z) &= xz\\ p_{2,2}(x,y,z) &= x^2-y^2\\ \\ p_{3,-3}(x,y,z) &= -3x^2y+y^3\\ p_{3,-2}(x,y,z) &= xyz\\ p_{3,-1}(x,y,z) &= y(x^2+y^2-4z^2)\\ p_{3,0}(x,y,z) &= z(-3(x^2+y^2)+2z^2)\\ p_{3,1}(x,y,z) &= x(x^2+y^2-4z^2)\\ p_{3,2}(x,y,z) &= (x^2-y^2)z\\ p_{3,3}(x,y,z) &= x^3-3xy^2 \end{align} $$

I vaguely remember encountering similar polynomials when doing multipole expansions of ion optic elements, such as magnetic fields used in particle accelerators. But the generating terms used in ion optics are typically two-dimensional.

I have tried finding how to generate such polynomials in the literature, but found nothing that did not use pretty heavy math. I am not entirely sure that "harmonic polynomials" is really the correct term form them. In any case, how can these polynomials be efficiently generated in a way that ensures that they are orthogonal?

Answer 1

Overview

These types of polynomials are used in quantum chemistry, potential theory, magnet shimming, and probably many other branches of science. One problem is that the nomenclature seems to be subtly different between the fields, so the following description is not definitive.

Let us call the solutions of the Laplace equation $\Delta\Phi = 0$, when they are stated in spherical coordinates, by the name of solid harmonics. There are two classes of solid harmonics: Regular solid harmonics, which vanish at the origin, and irregular solid harmonics, which have a singularity at the origin.

Solid harmonics can be parametrized using spherical harmonics $Y_l^m$. The parametrization of the regular solid harmonics is $r^l Y_l^m(\theta,\phi)$, and irregular solid harmonics can be stated as $r^{-(l+1)} Y_l^m(\theta,\phi)$, which immediately shows the singularity at the origin. These parametrizations of the regular and irregular solid harmonics have the benefit that they are orthonormal. (The corresponding scalar product is an integral over the unit sphere).

Both regular and irregular solid harmonics are used for spherical multipole expansions: Think of it as a Taylor expansion, or Fourier expansion, only that the base function are not $x^n$ or $\sin(n \omega t)$, but instead solid harmonics. This kind of expansion is very useful when trying to approximate some function that (at least locally) is known to satisfy the Laplace equation. When expanding it with solid harmonics, then even though you might stop at some multipole moment, you can be certain that the expansion still satisfies the Laplace equation.

When the charges are near the origin, and you want to describe the field (caused by these charges) far away from the origin, you would typically use the irregular solid harmonics as base functions: Their poles are on the origin, but that's OK, because all you care about is that the approximation far away from the origin is a good one. This is the "standard" multipole expansion.

For situations where you want to approximate the field near the origin, and the field is caused by charges far away, you would use the regular solid harmonics as base functions. This is called the interior multipole expansion, and it is used in beam optics and ion traps. Now the important part: When the regular solid harmonics are expressed in cartesian coordinates, the results turn out to be surprisingly simple polynomials in $x,y,z$. (Maybe it's not so surprising, after all, you could always do multidimensional Taylor expansion around the origin, and restrict some of the coefficients, so that the total polynomial always satisfies the Laplace equation.)

While solid harmonics are typically complex functions, a simple linear combination removes all instances of $i$, and you are left with polynomials in $x,y,z$ where all coefficients are real.

Generating Formula

There is a closed formula for generating real solid harmonics in cartesian coordinates. The resulting polynomials should, as far as I understand it, be orthogonal and normal over the unit sphere. Those that I checked turned out to be orthogonal over the cube $|r| \leq 1$, but not normalized! The following Mathematica code implements the generating formula.

gamma[l_, m_, k_] :=
  (-1)^k 2^(-l) Binomial[l, k] Binomial[2 l - 2 k, l] (l - 2 k)!/(l - 2 k - m)!;
pi[l_, m_] :=
  Sum[gamma[l, m, k] (x^2 + y^2 + z^2)^k z^(l - 2 k - m), {k, 0, Floor[(l - m)/2]}];
a[m_] :=
  Sum[Binomial[m, p] x^p y^(m - p) Cos[(m - p) Pi/2], {p, 0, m}];
b[m_] :=
  Sum[Binomial[m, p] x^p y^(m - p) Sin[(m - p) Pi/2], {p, 0, m}];
c[l_, m_] :=
  Sqrt[(2 - KroneckerDelta[m, 0]) (l - m)!/(l + m)!] pi[l, m] a[m];
s[l_, m_] :=
  Sqrt[2 (l - m)!/(l + m)!] pi[l, m] b[m];
genfun[l_, m_] := If[m >= 0, If[m <= l, c[l, m], 0], If[m >= -l, s[l, -m], 0]] // Simplify;
result =
 Table[Table[genfun[l, m], {m, {0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5}}], {l, 0, 5}] // TableForm

You can use

Map[Simplify[Laplacian[#, {x, y, z}]] &, result, {2}]

to ensure that all these polynomials satisfy the Laplace equation.

Alternative I

A different approach is to start from the two-dimensional multipole expansion of the field, and then generalize to the three-dimensional harmonic polynomials. This can be achieved by multiplying the two-dimensional polynomial with the third variable to any power, and then applying a special operator that I can only describe as funky. This procedure is outlined in the book Fundamental World of Quantum Physics. The corresponding Mathematica code is lengthy, but it can be evaluated rather quickly.

Alternative II

A third approach is to use Mathematica to build the solid harmonics from the spherical harmonics, then transforming them to cartesian coordinates, and finally combining them in a way that leaves only real-valued coefficients. This code is short (it is my original approach), but takes a while to execute. Also, Mathematica fails to simplify the polynomials of degree > 10 in a sensible way.

solidHarmonicS[l_?IntegerQ, m_?IntegerQ] := 
 Module[{r, \[Theta], \[Phi], xx, yy, zz}, 
  FullSimplify@
    Evaluate[
     TransformedField["Spherical" -> "Cartesian", 
      r^l SphericalHarmonicY[l, 
        m, \[Theta], \[Phi]], {r, \[Theta], \[Phi]} -> {xx, yy, 
        zz}]] /. {xx -> x, yy -> y, zz -> z}]
realSolidHarmonics[l_?IntegerQ, m_?IntegerQ] := 
 If[m == 0, solidHarmonicS[l, m], 
   If[m > 0, 
    1/Sqrt[2] (solidHarmonicS[l, m] + (-1)^m solidHarmonicS[l, -m]), 
    1/(I Sqrt[2]) (solidHarmonicS[l, 
        m] - (-1)^m solidHarmonicS[l, -m])]] // Simplify
result = Table[
   Table[realSolidHarmonics[i, 
     j], {j, {0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5}}], {i, 0, 5}] // 
  TableForm