Converting between the Cartesian tensor and spherical harmonic expansion of solutions to the Boltzmann equation

Abstract

The Vlasov–Fokker–Planck (VFP) equation, a variant of the Boltzmann equation, is frequently used to model the dynamics of laboratory and astrophysical plasmas. A common approach in solving the VFP equation in numerical and analytic studies is to expand the momentum part of the distribution function $f$ (i.e. the particle density in phase-space) in terms of known functions. The Cartesian tensor and spherical harmonic expansion have been widely used, leading to the question of how to convert between the different coefficients of the two expansions. This problem is also familiar in multipole expansions of an electrostatic (or gravitational) potential. The coefficients of the Cartesian tensor expansion of the potential are called (Cartesian) multipole moments and the ones of the spherical harmonic expansion are called spherical multipole moments. In this paper, we investigate the relation between the two kinds of multipole moments and provide a general formalism to convert between them. We subsequently apply this formalism to the coefficients of the expansions of the distribution function $f$. A free, open-source command-line tool which implements this formalism is provided.

1. Introduction

Plasma kinetic simulations are an essential tool in the study of weakly collisional/collisionless systems. Such conditions are realised in both astrophysical and laboratory plasmas. Collisions, either via binary Coulomb or frequent wave–particle interactions can be appropriately modelled with the aid of a Fokker–Planck approach (Chandrasekhar 1943). Thus treated, one arrives at the Vlasov–Fokker–Planck equation, which has the following form

(1.1)\begin{equation} \partial_{t} f + \boldsymbol{v}\boldsymbol{\cdot} \nabla_{\boldsymbol{r}} f + \boldsymbol{F} \boldsymbol{\cdot} \nabla_{\boldsymbol{p}}\, f = \left(\frac{\delta f}{\delta t}\right)_{\text{Sc}} . \end{equation}

Here, $f = f(t, \boldsymbol {r}, \boldsymbol {p})$ is the particle (number) density in the six-dimensional phase-space, or simply the distribution function. Here, $\boldsymbol {v}$ is the velocity of the particles with $\boldsymbol {F}$ the force acting on them. The derivatives are taken with respect to the $\boldsymbol {r}$- and $\boldsymbol {p}$-coordinates. The right-hand side of the equation represents the effect of the particles’ interactions; ‘Sc’ stands for ‘scattering’.

Several codes have been developed to solve the multi-dimensional Vlasov equation numerically (for example, von Alfthan et al. 2014; Juno et al. 2018). Simulations with this approach are expensive in multi-dimensions. When possible, it is pragmatic to exploit the fact that collisions tend to drive the distribution towards isotropy. A large body of work has exploited this fact to derive transport coefficients in ionized plasmas including Coulomb collisions, for example in the form of resistivity and conductivity tensors (e.g. Braginskii 1965). Reliable values for these coefficients are essential for accurate fluid models, necessary for fusion and other applications (see for example Thomas et al. (2012) for a review). In the case of small Knudsen numbers, i.e. when the Coulomb mean free path is short compared to the relevant macroscopic gradients of the system, it is common to adopt a perturbative approach, where the equilibrium solution is a local Maxwellian distribution (Spitzer & Härm 1953; Braginskii 1958). In the small Knudsen limit, the distribution is (following our notation below) well approximated by the leading-order terms in a Cartesian expansion $f(\boldsymbol {p}) \approx {\mathsf{F}}^{(0)}+ {\mathsf{F}}^{({1})}_i p^{i}/p$. The relevant transport coefficients can be determined by solving for $\boldsymbol{\mathsf{F}}^{({1})}(p)$ given $\boldsymbol{\mathsf{F}}^{(0)}(\boldsymbol {r}, p)$. In his seminal work, Braginskii (1958) derived approximate solutions using a truncated expansion in Laguerre (Sonine) polynomials. Other works applied a more direct numerical finite-difference approach (e.g. Epperlein & Haines 1986). The transport coefficients that result from the two approaches differ in certain regimes, which Epperlein & Haines (1986) attributed to the finite truncation of the Laguerre polynomial expansion. However, convergence is still only accurate to the expansion order in the Knudsen number, which for larger values necessitates more terms in the particle anisotropy.

To improve the angular space resolution, one can expand $f$ as a series of Laplace's spherical harmonics (Allis 1956; Bell et al. 2006; Reville & Bell 2013; Tzoufras et al. 2013)

(1.2)\begin{equation} f(t, \boldsymbol{r}, \boldsymbol{p}) = \sum^{\infty}_{l = 0} \sum^{l}_{m ={-}l} f^{m}_{l}(t,\boldsymbol{r}, p) Y^{m}_{l}(\theta, \varphi) . \end{equation}

Alternatively, one can expand $f$ in a series of Cartesian tensors (Johnston 1960; Shkarofsky, Johnston & Bachynski 1966; Epperlein & Haines 1986; Thomas et al. 2012)

(1.3)\begin{equation} f(t, \boldsymbol{r}, \boldsymbol{p}) = \sum^{\infty}_{l = 0} {\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}(t,\boldsymbol{r}, p) \frac{p^{i_{1}} \ldots p^{i_{l}}}{p^{l}} , \end{equation}

where $p^{1}, p^{2}$ and $p^{3}$ are the components of $\boldsymbol {p}$. The objects $\boldsymbol{\mathsf{F}}^{({l})}$ are called Cartesian tensors and its components ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ are the coefficients of the above expansion. Note that we are implicitly summing over all repeated indices. Moreover, the product of the components of $\boldsymbol {p}$ can be considered as a component of the (Cartesian) tensor $\otimes {\mathsf{p}}^{(l), i_{1} \ldots i_{l}} \equiv p^{i_{1}} \ldots p^{i_{l}}$. Hence, the summation over repeated indices may be seen as a contraction of two tensors (formerly called a Cartesian tensor scalar product).

Both expansions can be substituted into (1.1) and the coefficients $f^{m}_{l}$ and ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ can be computed numerically or, depending on the complexity of the problem, analytically. Since both expansions of $f$ represent the same distribution function, there must exist a relation between $f^{m}_{l}$ and ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$, see Courant & Hilbert (1953).

Johnston (1960) investigated this relation and derived a way to express the components of the Cartesian tensors ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ as a sum containing the $f^{m}_{l}$ for low values of $l$, namely $l \leq ~4$. We revisit the problem of relating $f^{m}_{l}$ and ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ and derive general formulae for converting between them (in both directions). We do not do this in the context of the distribution function $f$ and its expansions, but in the context of the multipole expansion of an electrostatic (or gravitational) potential. Also, the multipole expansion is either given as a Cartesian tensor or a spherical harmonic expansion. The expansion coefficients are called (Cartesian) multipole moments (i.e. monopole, dipole, quadrupole and high-order moments) or spherical multipole moments. The relation between these two kinds of moments is the same as the relation between ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ and $f^{m}_{l}$, but multipole moments play an important role in multiple branches of physics, such as general relativity (Thorne 1980), molecular physics (Cipriani & Silvi 1982), plasma physics (Johnston 1960; Epperlein & Haines 1986) and other areas like computational physics (Ludwig 1991).

Hence, in this paper, we derive a way to convert between the Cartesian multipole moments and the spherical multipole moments and apply it to convert between ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ and $f^{m}_{l}$. While the equivalence of the two multipole moments is long known (e.g. Courant & Hilbert 1953, p. 517), systematic methods that relate them are not readily found in the literature.

In § 2.1, we begin with deriving a definition of the components of the Cartesian multipole moments. We will see that harmonic and homogeneous polynomials are part of this definition.

A polynomial $p$ of degree $l$ is homogeneous if and only if

(1.4)\begin{equation} p(\lambda \boldsymbol{r}) = \lambda^{l} p(\boldsymbol{r}) \quad \text{for all } \lambda \in \mathbb{R} . \end{equation}

Moreover, a polynomial $p$ is harmonic if it is a solution to Laplace's equation, i.e.

(1.5)\begin{equation} {\rm \Delta} p = 0 . \end{equation}

We denote the space of homogeneous and harmonic polynomials of degree $l$ (in three variables) by $\mathcal {H}^{l}(\mathbb {R}^{3})$.

We call the polynomials, which are part of the Cartesian multipole moment's definition, multipole functions and we show that a subset of these polynomials is a basis of $\mathcal {H}^{l}(\mathbb {R}^{3})$. We call this subset multipole basis functions.

In § 2.2, we revisit the definition of the spherical multipole moments (prominently derived in Jackson 1998, p. 146). We will see that it contains (regular) solid harmonics. We prove that the solid harmonics are homogeneous and harmonic polynomials as well. Furthermore, we show that they are also a basis of $\mathcal {H}^{l}(\mathbb {R}^{3})$. We conclude that it is possible to find a basis transformation between the multipole basis functions and the solid harmonics. This basis transformation will allow us to express the spherical multipole moments in terms of the Cartesian multipole moments.

We can compute a preliminary basis transformation with the help of Efimov's definition of the Cartesian multipole moments (Efimov 1979), which we introduce in § 3. We prove that his definition and our definition are equivalent. In his definition, he uses a differential operator (Efimov's ladder operator) and, in conjunction with the ladder operator $L_{\pm }$ (known from the theory of angular momentum in Quantum Mechanics), Efimov derives the preliminary basis transformation, i.e. he derives a way to express the solid harmonics as a sum of the multipole functions.

In § 4, we transform Efimov's expression for the solid harmonics, taking into account that only a subset of the multipole functions is a basis of $\mathcal {H}^{l}(\mathbb {R}^{3})$, to become an actual basis transformation between the solid harmonics and the multipole basis functions.

In § 5, we write the basis transformation as a matrix equation and compute the inverse basis transformation by inverting the basis transformation matrix. We show that inverting the basis transformation matrix is equivalent to invert four upper triangular matrices, if real solid harmonics are used and if they and the multipole basis functions are ordered in a specific way. The inverse basis transformation gives an expression of the multipole basis functions in terms of the solid harmonics. We use the basis transformation matrix and its inverse to present the relation between the Cartesian multipole moments and the spherical multipole moments in the form of a matrix equation.

In § 6, we come back to our original problem, namely the Vlasov–Fokker–Planck equation and the expansions of the distribution function $f$. We apply the formalism, which we derived to convert between the Cartesian multipole moments and the spherical multipole moments, to solve it and reproduce the results found by Johnston (1960).

We conclude the paper with a summary of the formalism to convert between the Cartesian multipole moments and the spherical multipole moments and assess it critically. Furthermore, we provide a link to the code of a command-line tool, which implements it.

2. Definition of the multipole moments

In electrodynamics (or equivalently, on exchanging constants, in gravitational theory), the multipole expansion is an expansion of the solution to Poisson's equation, namely

(2.1)\begin{equation} 4 {\rm \pi}\epsilon_{0} \phi(\boldsymbol{r}) = \int \frac{\rho(\boldsymbol{r})}{|\boldsymbol{r} - \boldsymbol{r}'|} \,\mathrm{d}^{3}r' . \end{equation}

Either the denominator of the integrand is Taylor expanded in $r'/r$ or it is interpreted as the generating function of the Legendre polynomials. The former leads to the (Cartesian) multipole expansion with the Cartesian multipole moments and the latter to the spherical multipole expansion with the spherical multipole moments.

In the first part of this section, we look at the Cartesian multipole expansion and derive a definition of the Cartesian multipole moments. The second part dedicates itself to the spherical multipole expansion and investigates the definition of the spherical multipole moments.

2.1. Definition of the Cartesian multipole moments

In textbooks (e.g. Jackson 1998, p. 146), no general definition of the Cartesian multipole moments is given. Instead, the Taylor expansion of the denominator in the integrand of (2.1) is computed up to, for example, second order, i.e.

(2.2)\begin{equation} \frac{1}{|\boldsymbol{r} - \boldsymbol{r}'|} = \sum^{\infty}_{l = 0} \frac{({-}1)^{l}}{l!} (\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})^{l} \frac{1}{r} = \frac{1}{r} - x'_{i}\partial_{i}\frac{1}{r} + x'_{i}x'_{j} \partial_{i}\partial_{j}\frac{1}{r} - \cdots . \end{equation}

The potential then becomes

(2.3)\begin{align} 4{\rm \pi}\epsilon_{0} \phi(\boldsymbol{r}) = \int \rho(\boldsymbol{r}') \, \mathrm{d}^{3}r' \frac{1}{r} + \int \rho(\boldsymbol{r}')x'_{i} \, \mathrm{d}^{3}r' \frac{x^{i}}{r^{3}}+ \frac{1}{2} \int \rho(\boldsymbol{r}')x'_{i}x'_{j} \, \mathrm{d}^{3}r' \frac{(3x^{i}x^{j} - r^{2}\delta^{ij})}{r^{5}} + \cdots . \end{align}

The primed and unprimed components may be interchanged. As an example, consider the second-order term:

(2.4)\begin{equation} \frac{1}{2}\int \rho(\boldsymbol{r}') (3x'_{i}x'_{j} - r'^{2}\delta_{ij}) \, \mathrm{d}^{3} r' \frac{x^{i}x^{j}}{r^{5}} , \end{equation}

where we used that $r^{2}\delta ^{ij}x'_{i}x'_{j} = r'^{2}\delta _{ij} x^{i} x^{j}$. This transforms the potential into

(2.5)\begin{equation} 4{\rm \pi}\epsilon_{0} \phi(\boldsymbol{r}) = \frac{q}{r} + \frac{d_{i}}{r^{3}}x^{i} + \frac{{\mathsf{Q}}_{ij}}{r^{5}}x^{i}x^{j} + \cdots , \end{equation}

with the usual definitions of the monopole and the dipole. The components of the quadrupole moment are defined as

(2.6)\begin{equation} {\mathsf{Q}}_{ij} \equiv \int \rho(\boldsymbol{r}') (3x'_{i}x'_{j} - r'^{2}\delta_{ij}) \,\mathrm{d}^{3}r' . \end{equation}

Hence, a general method to define the components of the Cartesian multipole moment of an arbitrary order $l$ consists in taking $l$ partial derivatives of $1/r$ and interchanging the components of $\boldsymbol {r}$ and $\boldsymbol {r}'$, as in the above example. However, note that the interchange did not change the (functional) form of the numerator of $\partial _{i}\partial _{j}1/r$, it merely replaced the components of $\boldsymbol {r}$ with the components of $\boldsymbol {r}'$. This implies that we could alternatively have obtained the expression in parentheses in the quadrupole moment's definition by computing $\partial '_{i}\partial '_{j} 1/r'$ and ‘removing’ the denominator. This is the idea behind our definition of the components of the Cartesian multipole moment of rank $l$, which is

(2.7)\begin{equation} {\mathsf{Q}}^{({l})}_{i_1 \ldots i_{l}} \equiv \int \rho(\boldsymbol{r}') \mathcal{K}\left[({-}1)^{l}\partial'_{i_{1}} \ldots \partial'_{i_{l}} \frac{1}{r'} \right] \, \mathrm{d}^{3}r' . \end{equation}

Here, the calligraphic $\mathcal {K}$ denotes the Kelvin transform (Axler, Bourdon & Wade 2001, p. 61), which is defined as

(2.8)\begin{equation} \mathcal{K}[f](\boldsymbol{r}) \equiv \frac{1}{r}f(\tilde{\boldsymbol{r}}) \quad\text{with } \tilde{\boldsymbol{r}} \equiv \boldsymbol{r}/r^{2} . \end{equation}

Applying the definition to the $l = 2$ case recovers the quadrupole moment and demonstrates in which sense the Kelvin transform removes the denominator. With this definition at hand, the (Cartesian) multipole expansion is

(2.9)\begin{equation} 4 {\rm \pi}\epsilon_{0} \phi(\boldsymbol{r}) = \sum^{\infty}_{l = 0} \frac{1}{l!r^{2l + 1}} {\mathsf{Q}}^{({l})}_{i_{1} \ldots i_{l}}x^{i_{1}} \ldots x^{i_{l}} . \end{equation}

Of course, the first three multipole moments are the monopole ${\mathsf{Q}}^{({0})} = q$, the dipole moment ${\mathsf{Q}}^{({1})}_{i} = d_{i}$ and the quadrupole moment ${\mathsf{Q}}^{({2})}_{ij} = {\mathsf{Q}}_{ij}$.

We show now that the definition in (2.7) gives the known Cartesian multipole moments for an arbitrary rank $l$. The definition is correct if for all $l$, the interchange of $\boldsymbol {r}$ and $\boldsymbol {r}'$ does not change the form of the numerator of the $l$th partial derivative of $1/r$ and if the Kelvin transform removes its denominator.

A sufficient condition for the Kelvin transform to remove the denominator is that the numerator of the $l$th partial derivative of $1/r$ is an homogeneous polynomial (see (1.4)). That this is a sufficient condition can be seen by noting that

(2.10)\begin{equation} ({-}1)^{l} \partial_{i_{1}} \ldots \partial_{i_{l}} \frac{1}{r} = \frac{{\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r})}{r^{2l + 1}}, \end{equation}

where ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is an as yet undetermined function (the above equality can be proven by induction). The indices reflect that the function is different for different values of the indices. We call these functions multipole functions. Moreover, Axler, while proving Lemma 5.15 in Axler et al. (2001, p. 86), shows that ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is an homogeneous polynomial of degree $l$. Taking the Kelvin transform then yields

(2.11)\begin{equation} \mathcal{K}\left[({-}1)^{l}\partial_{i_{1}} \ldots \partial_{i_{l}} \frac{1}{r}\right] = r^{2l + 1} \frac{{\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r})}{r^{2l + 1}} = {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r}) , \end{equation}

which shows that it suffices that ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is an homogeneous polynomial to remove the denominator.

In the example of the second-order term, the interchange of $\boldsymbol {r}$ and $\boldsymbol {r}'$ did not change the (functional) form of ${\mathsf{M}}^{({2})}_{ij}$ because this form was ‘somehow’ special. We already know that ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is an homogeneous polynomial of degree $l$. We can further restrict its form. Remember that ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is the numerator of $g(\boldsymbol {r}) \equiv \partial _{i_{1}} \ldots \partial _{i_{l}} 1/r$. The following two considerations fix its form.

First, the form of the function $g$ does not change if it is defined in a rotated coordinate system. Assume that the relation between the rotated coordinates $\hat {\boldsymbol {r}}$ and the original coordinates $\boldsymbol {r}$ is $\hat {\boldsymbol {r}} = \boldsymbol{\mathsf{A}}\boldsymbol {r}$. The former statement then says that $\hat {g}(\hat {\boldsymbol {r}}) = g(\hat {\boldsymbol {r}})$, where $\hat {g}(\hat {\boldsymbol {r}}) \equiv \hat {\partial }_{i_{1}} \ldots \hat {\partial }_{i_{l}} 1/\hat {r}$. A comparison of $g$'s definition and $\hat {g}$ proves it. This restricts the form of the numerator of $g$, as the following example for $l = 2$ illustrates:

(2.12)\begin{equation} \hat{g}(\hat{\boldsymbol{r}}) = \hat{\partial}_{i}\hat{\partial_{j}}\frac{1}{\hat{r}} = ({\mathsf{A}})_{i}^{k}\partial_{k} ({\mathsf{A}})_{j}^{l}\partial_{l} \frac{1}{r}= ({\mathsf{A}})_{i}^{k}({\mathsf{A}})_{j}^{l} \frac{3x_{k}x_{l} - r^{2}\delta_{kl}}{r^{5}} = \frac{3\hat{x}_{i}\hat{x}_{j} - \hat{r}^{2}\delta_{ij}}{\hat{r}^{5}} = g(\hat{\boldsymbol{r}}) , \end{equation}

where we used that $\hat {r} = r$ and $\hat {\boldsymbol {\nabla }} = \boldsymbol{\mathsf{A}} \boldsymbol {\nabla }$. Note the rotation matrices must either transform a component of $\boldsymbol {r}$ or they must be multiplied with each other, i.e. $({\mathsf{A}})_{i}^{k}({\mathsf{A}})_{j}^{l} \delta _{kl}= ({\mathsf{A}})_{i}^{k}({\mathsf{A}}^{T})_{kj} = \delta _{ij}$. Hence, all the terms of ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ must be products of components of $\boldsymbol {r}$ and Kronecker deltas to ‘cope’ with the $l$ rotation matrices which appear as a consequence of $\hat {g}(\hat {\boldsymbol {r}}) = g(\hat {\boldsymbol {r}})$.

Second, the exchange of any two of the indices in $g$'s definition does not change its numerator, because it is possible to exchange the partial derivatives with each other (Schwarz's theorem). Hence, the object ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is symmetric in all its indices. This implies that a sum of the mentioned products of components of $\boldsymbol {r}$ and Kronecker deltas is needed. For example, consider the $l = 3$ case,

(2.13)\begin{equation} {\mathsf{M}}^{({3})}_{ijk} = c_{3,0} x_{i}x_{j}x_{k} + c_{3,1}r^{2}(x_{i}\delta_{jk} + x_{j}\delta_{ik} + x_{k}\delta_{ij}), \end{equation}

where $c_{3,0}, c_{3,1} \in \mathbb {R}$ are coefficients. Note that the sum is over pairs of indices and that this guarantees that any two of the indices can be exchanged without changing ${\mathsf{M}}^{({3})}_{ijk}$. The factor $r^{2}$ reflects that ${\mathsf{M}}^{({3})}_{ijk}$ must be an homogeneous polynomial of degree three.

In summary, the fact that ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is an homogeneous polynomial of degree $l$ that can be invariantly defined in a rotated coordinate system and that is symmetric in its indices implies that its functional form must be

(2.14)\begin{equation} {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}} = \sum^{\lfloor {l}/{2} \rfloor}_{k=0} c_{l,k} r^{2k} R_{l,2k}, \end{equation}

where $\lfloor x \rfloor$ is the floor function. We also introduce $R_{l,2k} \equiv {P}(\delta _{i_{1}i_{2}} \ldots \delta _{i_{2k-1}i_{2k}} x_{i_{2k + 1}} \ldots x_{i_{l}})$, where ${P}$ is an operator which produces the sum over the pairs of indices needed to assure symmetry (cf. with the previous example). This operator hides a lot of the complexity, but is explained in detail in Appendix A.2. A closed-form expression for the coefficients $c_{l,k}$ is derived in Appendix A.3. Similar arguments are used by Efimov to fix the functional form of the multipole functions (Efimov 1979, p. 426).

This is the ‘somehow’ special (functional) form of the multipole functions ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ which permits that the interchange of $\boldsymbol {r}$ and $\boldsymbol {r}'$ does not change the form but only replaces the components of $\boldsymbol {r}$ with the components of $\boldsymbol {r}'$. To see this, note that, in the light of (2.10) and (2.14), a term in the Taylor expansion of (2.2) of arbitrary order $l$ contains

(2.15)\begin{equation} r^{2k} x'^{i_{1}} \ldots x'^{i_{l}} {P}(\delta_{i_{1}i_{2}} \ldots \delta_{i_{2k-1}i_{2k}} x_{i_{2k+1}} \ldots x_{i_{l}}). \end{equation}

Using that the exponent of $r^{2k}$ is even and remembering that $r^{2}\delta ^{ij}x'_{i}x'_{j} = r'^{2}\delta _{ij}x^{i}x^{j}$, we can replace the components of $\boldsymbol {r}$ and $\boldsymbol {r}'$.

We conclude that the definition of the Cartesian multipole moments introduced in (2.7) is correct. The $\boldsymbol{\mathsf{Q}}^{({l})}$ are tensors, i.e. they transform as required under coordinate transformations. In Appendix A.1, we use our definition to prove this.

Before we close this section, we emphasise that the multipole functions ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ are harmonic functions, i.e.

(2.16)\begin{equation} {\rm \Delta} {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}} = 0 . \end{equation}

This is also proven in (the already quoted) Lemma 5.15 in Axler et al. (2001, p. 86).

In what follows, it is convenient to introduce a different notation for the indices of $\boldsymbol{\mathsf{M}}^{({l})}$, namely

(2.17)\begin{equation} {\mathsf{M}}^{({l})}_{pqr} \equiv \mathcal{K}\left[({-}1)^{l} \partial^{p}_{x} \partial^{q}_{y} \partial^{r}_{z} \frac{1}{r}\right] \quad \text{with } p + q + r = l. \end{equation}

This notation implicitly takes advantage of the symmetry. For example, ${\mathsf{M}}^{({3})}_{112} = {\mathsf{M}}^{({3})}_{121}$ in the ${\mathsf{M}}^{({l})}_{i_1\dots i_l}$ notation and both components are ${\mathsf{M}}^{({3})}_{210}$ in the newly introduced $pqr$ notation.

Last, but not least, Axler shows in Theorem 5.25 in Axler et al. (2001, p. 92) that

(2.18)\begin{equation} \mathcal{H}^{l}(\mathbb{R}^{3}) = \operatorname{span}\{{\mathsf{M}}^{({l})}_{pqr} \mid p \leq 1 \text{ and } p + q + r = l\} . \end{equation}

In words, the functions ${\mathsf{M}}^{({l})}_{pqr}$ are homogeneous and harmonic polynomials of degree $l$, and a subset of them is a basis of the space of these polynomials. We call this subset the multipole basis functions.

2.2. Definition of the spherical multipole moments

The denominator of the integrand in (2.1) can be interpreted as the generating function of the Legendre polynomials, i.e.

(2.19)\begin{equation} \frac{1}{|\boldsymbol{r} - \boldsymbol{r}'|} = \sum^{\infty}_{l = 0} \frac{r'^{l}}{r^{l+1}} P_{l}(\cos\gamma) . \end{equation}

Here, $\gamma$ is the angle between $\boldsymbol {r}$ and $\boldsymbol {r}'$ (Jackson 1998, p. 102, (3.38)). Together with the addition theorem for Laplace's spherical harmonics (cf. Jackson 1998, p. 110, (3.62)), which is

(2.20)\begin{equation} P_{l}(\cos\gamma) = \frac{4{\rm \pi}}{2l + 1} \sum^{l}_{m ={-}l} {Y^{m}_{l}}^{*}(\theta', \varphi') Y^{m}_{l}(\theta, \varphi) , \end{equation}

we obtain the spherical multipole expansion

(2.21)\begin{equation} 4 {\rm \pi}\epsilon_{0} \phi(\boldsymbol{r}) = \sum^{\infty}_{l = 0} \sum^{l}_{m ={-}l} \frac{4 {\rm \pi}}{2l + 1} \frac{q^{m}_{l}}{r^{2l + 1}} r^{l}Y^{m}_{l}(\theta, \varphi) . \end{equation}

The function $r^{l}Y^{m}_{l}(\theta, \varphi )$ is called a solid harmonic of degree $l$ and order $m$. Furthermore, we call the coefficients $q^{m}_{l}$ spherical multipole moments (Jackson 1998, p. 145). They are defined as

(2.22)\begin{equation} q^{m}_{l} \equiv \int \rho(\boldsymbol{r}') r'^{l}{Y^{m}_{l}}^{*}(\theta', \varphi') \, \mathrm{d}^{3}r' . \end{equation}

A (Laplace's) spherical harmonic of degree $l$ and order $m$ is (cf. Jackson 1998, p. 108, (3.53))

(2.23)\begin{equation} Y^{m}_{l}(\theta, \varphi) \equiv \sqrt{\frac{2l + 1}{4{\rm \pi}} \frac{(l - m)!}{(l + m)!}} P^{m}_{l}(\cos\theta)\, {\rm e}^{\mathrm{i} m\varphi} \equiv N^{m}_{l} P^{m}_{l}(\cos\theta)\, {\rm e}^{\mathrm{i} m\varphi}. \end{equation}

Here, $N^{m}_{l}$ is the normalisation and the functions $P^{m}_{l}$ are the associated Legendre polynomials.

The solid harmonics are homogeneous polynomials of degree $l$. This can be seen by using the following definition of the associated Legendre polynomials (cf. Jackson 1998, p. 108, (3.49)):

(2.24)\begin{equation} P^{m}_{l}(\cos\theta) = ({-}1)^{m}\sin^{m}\theta \frac{\mathrm{d}^{m}}{\mathrm{d} \cos^{m}\theta} P_{l}(\cos\theta) . \end{equation}

Then, the solid harmonics with order greater than or equal to zero ($m \geq 0$) are

(2.25)\begin{align} r^{l}Y^{m}_{l}(\theta, \varphi) & = ({-}1)^{m} N^{m}_{l} r^{l-m} \frac{\mathrm{d}^{m}}{\mathrm{d}\cos^{m}\theta} P_{l}(\cos\theta) (r\sin\theta\cos\varphi + \mathrm{i} r\sin\theta\sin\varphi)^{m} \nonumber\\ & = ({-}1)^{m} N^{m}_{l}\sum^{\lfloor({l-m})/{2}\rfloor}_{k = 0} a_{k}r^{2k}z^{l - m - 2k} (x + \mathrm{i} y)^{m} \equiv p(\boldsymbol{r}). \end{align}

The $a_{k} \in \mathbb {R}$ are coefficients, which collect numerical factors. In the second equation, we expressed the spherical coordinates in Cartesian coordinates and exploited that the $m$th derivative of a degree $l$ polynomial yields a degree $l-m$ polynomial. The Legendre polynomials, in particular, consist of either odd or even powers of $\cos \theta$. It can be checked that the definition of an homogeneous polynomial of degree $l$, which was given in (1.4), applies to $p$ as defined in (2.25). The argument holds true for negative $m$ as well. In this case, $(x + \mathrm {i} y)^{m}$ becomes $(x - \mathrm {i} y)^{m}$. The definition of an homogeneous polynomial still applies.

Furthermore, the solid harmonics (as their name suggests) are harmonic functions. In spherical coordinates the Laplace operator can be split into a radial part and an angular part, i.e. $\varDelta = \varDelta _{r} + \varDelta _{\theta, \varphi }/r^{2}$. Laplace's spherical harmonics are eigenfunctions of the angular part with eigenvalue $-l(l+1)$ (Landau & Lifshitz 1977, p. 91, (28.7)). Hence,

(2.26)\begin{equation} {\rm \Delta} r^{l}Y^{m}_{l} = Y^{m}_{l}\left(\partial^{2}_{r}r^{l} + \frac{2}{r}\partial_{r}r^{l}\right) + r^{l-2}\varDelta_{\theta, \varphi}Y^{m}_{l} = l(l+1) r^{l-2}Y^{m}_{l} - l(l+1)r^{l-2}Y^{m}_{l} = 0 . \end{equation}

We conclude this section by proving

(2.27)\begin{equation} \mathcal{H}^{l}(\mathbb{R}^{3}) = \operatorname{span}\{r^{l}Y^{m}_{l} \mid{-}l \leq m \leq l \} , \end{equation}

i.e. the solid harmonics are also a basis of the space of homogeneous and harmonic polynomials.

We proceed in two steps. First, we determine the dimension of $\mathcal {H}^{l}(\mathbb {R}^{3})$ and, second, we show that the solid harmonics are independent. The dimension of $\mathcal {H}^{l}(\mathbb {R}^{3})$ is $2~l + 1$. This can be derived in different ways, for example, as in Müller (1966, p. 11, (11)) or, alternatively, as done by Axler in Proposition 5.8 in Axler et al. (2001, p. 78). Note that there are $2l + 1$ solid harmonics of degree $l$. The linear independence follows from the orthogonality of Laplace's spherical harmonics, namely (cf. Jackson 1998, p. 108, (3.55))

(2.28)\begin{equation} \int^{2{\rm \pi}}_{0}\int^{\rm \pi}_{0} {Y^{m'}_{l'}}^{*}Y^{m}_{l} \sin\theta \,\mathrm{d}\theta \,\mathrm{d}\varphi = \delta_{l'l}\delta_{m'm}. \end{equation}

Hence, (2.27) is proven.

We gave explicit definitions of the Cartesian multipole moments and the spherical multipole moments. Both definitions contained homogeneous and harmonic polynomials: the multipole functions ${\mathsf{M}}^{({l})}_{pqr}$ in the former case, solid harmonics $r^{l}Y^{m}_{l}$ in the latter. Moreover, a subset of the functions ${\mathsf{M}}^{({l})}_{pqr}$ (namely, the multipole basis functions) are a basis of $\mathcal {H}^{l}(\mathbb {R}^{3})$. The same is true of the solid harmonics. Hence, the relation between the Cartesian multipole moments and spherical multipole moments can be formalised as a basis transformation between the solid harmonics and the multipole basis functions. In the next section, we introduce Efimov's ladder operator to derive the basis transformation.

3. Efimov's ladder operator

3.1. An alternative definition of the Cartesian multipole moments

In their paper on multipole moments, Efimov derives an alternative definition of the Cartesian multipole moments (Efimov 1979). To derive it, they multiply the Taylor expansion given in (2.2) by $r$, which yields

(3.1)\begin{equation} r \frac{1}{|\boldsymbol{r} - \boldsymbol{r}'|} = \sum^{\infty}_{l = 0} r \frac{(-\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})^{l}}{l!}\frac{1}{r} = \sum^{\infty}_{l = 0} \frac{{\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r})}{l! r^{2l}} x'_{i_{1}} \ldots x'_{i_{l}} . \end{equation}

Taking advantage of the fact that the functions ${\mathsf{M}}^{({l})}_{i_{1} \ldots {i_{l}}}$ are homogeneous polynomials of degree $l$ and using the definition of $\tilde {\boldsymbol {r}}$ (see (2.8)) leads to the equivalent equation

(3.2)\begin{equation} \sum^{\infty}_{l = 0} \frac{1}{\tilde{r}} \frac{(-\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})^{l}}{l!}\tilde{r} = \sum^{\infty}_{l = 0} \frac{1}{l!}{\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\tilde{\boldsymbol{r}}) x'_{i_{1}} \ldots x'_{i_{l}} . \end{equation}

Since every term of the sums on both sides of the equation is an homogeneous polynomial of the same degree, one finds

(3.3)\begin{equation} ({-}1)^{l} \frac{1}{\tilde{r}}(\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})^{l}\tilde{r} = {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\tilde{\boldsymbol{r}}) x'^{i_{1}} \ldots x'^{i_{l}} . \end{equation}

Note the derivatives on the left-hand side of the equation are taken with respect to $\boldsymbol {r}$ (and not $\tilde {\boldsymbol {r}}$). We may rewrite the above to get a convenient grouping:

(3.4)\begin{equation} ({-}1)^{l} \overbrace{\frac{1}{\tilde{r}}(\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})\tilde{r}\underbrace{\frac{1}{\tilde{r}}(\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})\tilde{r} \ldots \frac{1}{\tilde{r}}(\boldsymbol{r}' \boldsymbol{\cdot} \boldsymbol{\nabla})\tilde{r}}_{{\equiv} f(\tilde{\boldsymbol{r}})}}^{l \text{ times}} = ({-}1)^{l+1}x'^{i} (2\tilde{x}_{i}\tilde{x}^{j}\tilde{\partial}_{j} - \tilde{r}^{2}\tilde{\partial}_{i} + \tilde{x}_{i})f(\tilde{\boldsymbol{r}}) . \end{equation}

The expression in parentheses is Efimov's ladder operator. We introduce the notation

(3.5)\begin{equation} {{D}}_{i} \equiv (2x_{i}x^{j}\partial_{j} - r^{2}\partial_{i} + x_{i}), \end{equation}

which on applying $l$ times turns (3.3) into

(3.6)\begin{equation} (\boldsymbol{r}' \boldsymbol{\cdot} \tilde{\boldsymbol{{{D}}}})^{l}\boldsymbol{1} = {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\tilde{\boldsymbol{r}})x'^{i_{1}} \ldots x'^{i_{l}}. \end{equation}

Here, $\mathbf {1}(\boldsymbol {r}) = 1$ is a constant function. This implies an alternative definition of the multipole functions

(3.7)\begin{equation} {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r}) \equiv {{D}}_{i_{1}} \ldots {{D}}_{i_{l}} \boldsymbol{1} . \end{equation}

Replacing the Kelvin transform in (2.7) with Efimov's ladder operator yields an alternative definition for the Cartesian multipole moments as well.

We explicitly prove that the two definitions of the multipole functions are equivalent, i.e.

(3.8)\begin{equation} \mathcal{K}\left[({-}1)^{l}{\partial_{i_{1}} \ldots \partial_{i_{l}}\frac{1}{r}}\right] = {{D}}_{i_{1}} \ldots {{D}}_{i_{l}} \boldsymbol{1}. \end{equation}

The first step consists in showing that

(3.9)\begin{equation} \mathcal{K}\left[({-}1)^{l}\partial_{i_{1}} \ldots \partial_{i_{l}} \frac{1}{r}\right] =\prod^{l}_{k = 1} ((2l - (2k - 1)) x_{i_{k}} - r^{2}\partial_{i_{k}})\boldsymbol{1} . \end{equation}

This can be proven by induction. A quick computation shows that the statement holds true for the base case $l = 1$. From (2.10), and noting that $\mathcal {K}[\mathcal {K}[f]] = f$, it follows that

(3.10)\begin{equation} ({-}1)^{l}\partial_{i_{1}} \ldots \partial_{i_{l}} \frac{1}{r} = \mathcal{K}[ {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}](\boldsymbol{r}) = \frac{{\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r})}{r^{2l + 1}} . \end{equation}

Taking the derivative with respect to $x_{j}$ yields

(3.11)\begin{equation} \partial_{j}\partial_{i_{1}} \ldots \partial_{i_{l}} \frac{1}{r} = \frac{({-}1)^{l + 1}}{r^{2l + 3}}((2l + 1) x_{j} - r^{2}\partial_{j}) {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r}) . \end{equation}

Note that acting with the expression in parentheses on ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ gives an homogeneous polynomial of degree $l + 1$. We again take the Kelvin transform, use the induction hypothesis and see that

(3.12)\begin{align} \mathcal{K}\left[({-}1)^{l+1}\partial_{j}\partial_{i_{1}} \ldots \partial_{i_{l}} \frac{1}{r}\right] & = ((2l + 1) x_{j} - r^{2}\partial_{j}) {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}(\boldsymbol{r})\nonumber\\ & =((2l + 1) x_{j} - r^{2}\partial_{j}) \prod^{l}_{k = 1} ((2l - (2k - 1)) x_{i_{k}} - r^{2}\partial_{i_{k}})\boldsymbol{1}. \end{align}

Relabelling the indices ($j \rightarrow i_{1}$ and $i_{k} \rightarrow i_{k + 1}$) and shifting the index of the product by one, yields the conclusion.

In a second step, we prove that each factor in the product in (3.9) is equivalent to Efimov's ladder operator. We choose an arbitrary factor with index $k$ equal to $n$ and express the right-hand side of (3.9) as

(3.13)\begin{equation} ((2l - 1)x_{i_{1}} - r^{2}\partial_{i_{1}}) \times \cdots \times (2(l - n)x_{i_{n}} - r^{2}\partial_{i_{n}} + x_{i_{n}}) {\mathsf{M}}^{({l-n})}_{i_{n+1}\ldots i_{l}}(\boldsymbol{r}), \end{equation}

where we used that

(3.14)\begin{align} \prod^{l}_{k = n + 1}((2l - (2k - 1)) x_{i_{k}} - r^{2}\partial_{i_{k}}) \boldsymbol{1}& = \prod^{l - n}_{k = 1}((2(l-n) - (2k - 1)) x_{i_{k + n}} - r^{2}\partial_{i_{k + n}})\boldsymbol{1} \nonumber\\ & = {\mathsf{M}}^{({l-n})}_{i_{n+1} \ldots i_{l}}(\boldsymbol{r}) . \end{align}

To show that the factor in (3.13) corresponding to the index $k = n$ and ${{D}}_{i_{n}}$ are equivalent, we have to prove that

(3.15)\begin{equation} 2(l - n)x_{i_{n}}{\mathsf{M}}^{({l-n})}_{i_{n+1} \ldots i_{l}} = 2x_{i_{n}}x^{m}\partial_{m}{\mathsf{M}}^{({l-n})}_{i_{n+1} \ldots i_{l}} . \end{equation}

This equality holds because the function $u \equiv {\mathsf{M}}^{({l-n})}_{i_{n+1} \ldots i_{l}}$ is an homogeneous polynomial of degree $l-n$. Note that every homogeneous polynomial of an arbitrary degree $l$ can be written as a linear combination of monomials $x^{j_{1}} \ldots x^{j_{l}}$, i.e. $p(\boldsymbol {r}) = \alpha _{j_{1} \ldots j_{l}}x^{j_{1}} \ldots x^{j_{l}}$. Here, $\alpha$ is a collection of coefficients. There are $(l + 2)(l + 1)/2$ different monomials. If the object ${\boldsymbol{\mathsf{\alpha}}}$ is symmetric in all its indices, it contains an equal amount of independent coefficients. Thus, there exists a symmetric ${\boldsymbol{\mathsf{\alpha}}}$ such that $u(\boldsymbol {r}) = \alpha _{j_{1} \ldots j_{l-n}} x^{j_{1}} \ldots x^{j_{l-n}}$. This implies that

(3.16)\begin{align} x^{m}\partial_{m}u(\boldsymbol{r}) & = x^{m}\alpha_{j_{1} \ldots j_{l-n}}\partial_{m}(x^{j_{1}} \ldots x^{j_{l-n}}) = (l-n) x^{m}\alpha_{j_{1} \ldots j_{l-n}}\delta^{j_{1}}_{m} \ldots x^{j_{l-n}} \nonumber\\ & = (l-n) \alpha_{j_{1} \ldots j_{l-n}}x^{j_{1}} \ldots x^{j_{l-n}} = (l-n) u(\boldsymbol{r}). \end{align}

In the second line, we used that ${\boldsymbol{\mathsf{\alpha}}}$ is symmetric. Since $n$ was arbitrary, we conclude that all factors in the product in (3.9) are equivalent to a corresponding component of the ladder operator $\boldsymbol {{{D}}}$.

Before we end this section, we use the alternative definition of the Cartesian multipole moments (and the multipole functions) to learn more about them.

We first note that Efimov's ladder operators commute, i.e.

(3.17)\begin{equation} [{{D}}_{i}, {{D}}_{j}] \equiv {{D}}_{i} {{D}}_{j} - {{D}}_{j}{{D}}_{i}= 0 . \end{equation}

Note that this must be the case, because ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ must be symmetric. Moreover,

(3.18)\begin{equation} {{D}}_{i}{{D}}_{i} = r^{4}\varDelta . \end{equation}

An immediate consequence of these two equations is that

(3.19)\begin{equation} {\mathsf{M}}^{({l})}_{i_{1} \ldots i_{k} \ldots i_{k} \ldots i_{l}} = 0 . \end{equation}

We say that the object ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ is traceless, i.e. the contraction of two arbitrary indices is zero (cf. Efimov 1979, (1.5)–(1.7)). This property is carried over to the Cartesian multipole moments. Thus, the Cartesian multipole moments $\boldsymbol{\mathsf{Q}}^{({l})}$ are symmetric and traceless tensors of rank $l$.

We pointed out that an homogeneous polynomial of degree $l$ can be written as a linear combination of monomials with the help of a symmetric collection of coefficients ${\boldsymbol{\mathsf{\alpha}}}$. It can be shown that if the polynomial is harmonic as well, ${\boldsymbol{\mathsf{\alpha}}}$ must be traceless (Jeevanjee 2011, ex. 3.26). This implies that not only the multipole functions ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$ are homogeneous and harmonic polynomials, but also the functions in the Cartesian multipole expansion in (2.9), namely

(3.20)\begin{equation} {\mathsf{Q}}^{({l})}_{i_{1} \ldots i_{l}} x^{i_{1}} \ldots x^{i_{l}} , \end{equation}

are homogeneous and harmonic polynomials of degree $l$.

3.2. Preliminary basis transformation

With the tools derived up to now, we are in a position to express the solid harmonics as a sum of the multipole functions ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$.

We start again with the denominator of the integrand in (2.1). We fix $\boldsymbol {r}'$ to be the unit vector pointing into the $z$-direction. This implies

(3.21)\begin{equation} \frac{1}{|\boldsymbol{r} - \boldsymbol{e}_{z}|} = \sum^{\infty}_{l = 0} \frac{(-\boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{\nabla})^{l}}{l!}\frac{1}{r} = \sum^{\infty}_{l = 0} \frac{1}{r^{l+1}} P_{l}(\cos\theta) . \end{equation}

We used the Taylor expansion given in (2.2) for the left-hand side and (2.19) for the right-hand side of the last equation. Since $\boldsymbol {r}' = \boldsymbol {e}_{z}$, the angle $\gamma$ is the polar angle $\theta$. We take the Kelvin transform of this equation, i.e.

(3.22)\begin{equation} \sum^{\infty}_{l = 0} \frac{1}{l!}\mathcal{K}\left[({-}1)^{l}\partial^{l}_{z}\frac{1}{r}\right] = \sum^{\infty}_{l = 0} \mathcal{K}\left[\frac{1}{r^{2l+1}} r^{l}P_{l}(\cos\theta)\right] . \end{equation}

Note that the Kelvin transform's definition implies that it is linear, namely $\mathcal {K}[f + \lambda g] = \mathcal {K}[f] + \lambda \mathcal {K}[g]$. We now use that $N^{0}_{l} r^{l} P_{l} = r^{l}Y^{0}_{l}$. Furthermore, since a solid harmonic of degree $l$ is an homogeneous polynomial of degree $l$, the Kelvin transform removes the denominator of the right-hand side of the last equation (cf. § 2.1). Together with the equivalence of the Kelvin transform of the partial derivatives of $1/r$ and Efimov's ladder operator, we obtain the equivalent statement

(3.23)\begin{equation} \sum^{\infty}_{l = 0} r^{l}P_{l}(\cos\theta) = \sum^{\infty}_{l = 0}\frac{1}{l!}(\boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{{{D}}})^{l}\boldsymbol{1}. \end{equation}

Since every term on both sides of the equation is a polynomial of the same degree (namely $l$), the equation implies that

(3.24)\begin{equation} r^{l}Y^{0}_{l} = \frac{1}{l!}\sqrt{\frac{2l + 1}{4{\rm \pi}}}(\boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{{{D}}})^{l}\boldsymbol{1} = \frac{1}{l!}\sqrt{\frac{2l + 1}{4{\rm \pi}}} {\mathsf{M}}^{({l})}_{00l}. \end{equation}

Compare with Efimov (1979, p. 428, (2.9)) and note that the notation introduced in (2.17) is used for the indices of the multipole function.

To express a solid harmonic of order $m$ greater than zero as a sum of multipole functions, we use that Laplace's spherical harmonics $Y^{m}_{l}$ appear as the common eigenfunctions of the square of the angular momentum operator (known from Quantum Mechanics) and one of its components. The angular momentum operator is defined as (cf. Landau & Lifshitz 1977, p. 83, (26.2))

(3.25)\begin{equation} \textbf{L} \equiv{-}\mathrm{i} \boldsymbol{r} \times \boldsymbol{\nabla}. \end{equation}

We are particularly interested in the angular momentum ladder operator

(3.26)\begin{equation} \mathrm{L}_{{\pm}} = \mathrm{L}_{x} \pm \mathrm{i} \mathrm{L}_{y} ={-}\mathrm{i} (y \partial_{z} - z \partial_{y} \pm \mathrm{i} (z \partial_{x} - x \partial_{z})) = {\rm e}^{{\pm} \mathrm{i} \varphi}({\pm} \partial_{\theta} + \mathrm{i} \cot\theta \partial_{\phi}) . \end{equation}

The last line shows the ladder operator in spherical coordinates (Landau & Lifshitz 1977, p. 85, (26.15)). The ladder operators are used to increase (or decrease) the order of a spherical harmonic, e.g.

(3.27)\begin{equation} Y^{m}_{l} = ((l + m)(l - m + 1))^{{-}1/2}\mathrm{L}_{+}Y^{m-1}_{l}. \end{equation}

The numerical factor is an implication of (27.12) in Landau & Lifshitz (1977, p. 89). We can apply this operator $m$ times to (3.24) and we get

(3.28)\begin{equation} r^{l}Y^{m}_{l} = \sqrt{\frac{(l-m)!}{(l+m)!}}r^{l}\mathrm{L}^{m}_{+}Y^{0}_{l} = \frac{1}{l!} \sqrt{\frac{2l + 1}{4{\rm \pi}}\frac{(l-m)!}{(l+m)!}} \mathrm{L}^{m}_{+}{{D}}^{l}_{z}\boldsymbol{1} , \end{equation}

where we used the fact that $\mathrm {L}_{+}$ does not contain a derivative with respect to $r$, as can be seen in (3.26). We can compute the commutator of $\mathrm {L}_{+}$ and ${{D}}_{z}$. This yields

(3.29)\begin{equation} [L_{+}, {{D}}_{z}] ={-}({{D}}_{x} + \mathrm{i} {{D}}_{y}) . \end{equation}

This computation is done in detail in Appendix B.

One can prove by induction that

(3.30)\begin{equation} \mathrm{L}^{m}_{+}{{D}}^{l}_{z} \boldsymbol{1} = ({-}1)^{m} \frac{l!}{(l-m)!} {{D}}^{l-m}_{z} ({{D}}_{x} + \mathrm{i} {{D}}_{y})^{m} \boldsymbol{1} , \end{equation}

which, when substituted into the previously derived expression for $r^{l}Y^{m}_{l}$, yields the final result of this section: the preliminary basis transformation, i.e.

(3.31)\begin{align} r^{l} Y^{m}_{l} & = a^{m}_{l} {{D}}^{l-m}_{z} ({{D}}_{x} + \mathrm{i} {{D}}_{y})^{m}\boldsymbol{1} = a^{m}_{l} \sum^{m}_{p = 0} \mathrm{i}^{m-p}\binom{m}{p} {{D}}^{p}_{x} {{D}}^{m-p}_{y} {{D}}^{l-m}_{z} \boldsymbol{1} \nonumber\\ & = a^{m}_{l} \sum^{m}_{p = 0} \mathrm{i}^{m-p}\binom{m}{p} {\mathsf{M}}^{({l})}_{p(m-p)(l-m)} . \end{align}

Here, $m \geq 0$ and we used the notation of (2.17) for the indices of the multipole functions. The factor is

(3.32)\begin{equation} a^{m}_{l} \equiv ({-}1)^{m} \sqrt{\frac{2l + 1}{4{\rm \pi}}\frac{1}{(l+m)! (l-m)!}} . \end{equation}

We called this expression for the solid harmonics a preliminary basis transformation, because the sum on the right-hand side of (3.31) contains multipole functions with $p > 1$, which are not multipole basis functions. The solid harmonics with negative order $m$ are obtained with the help of $Y^{-m}_{l} = (-1)^{m}{Y^{m}_{l}}^{*}$ (Jackson 1998, p. 146, (4.3), (4.7)). Efimov derived (3.31) with similar arguments (see Efimov 1979, p. 429, (2.10)).

The preliminary basis transformation already allows to express spherical multipole moments in terms of Cartesian multipole moments. We can take the complex conjugate of (3.31) and plug it into the definition of the spherical multipole moments (see (2.22)) and obtain

(3.33)\begin{equation} q^{m}_{l} = a^{m}_{l} \sum^{m}_{p = 0} (-\mathrm{i})^{m-p}\binom{m}{p}{\mathsf{Q}}^{({l})}_{p(m-p)(l-m)} . \end{equation}

In his textbook ‘Classical Electrodynamics’, Jackson computes the spherical multipole moments with degree $l = 0, 1$ and $2$ (Jackson 1998, p. 146). The above formula allows a computation for arbitrary $l$ and $m$. For example,

(3.34)\begin{align} q^{3}_{3} & = \frac{1}{24}\sqrt{\frac{7}{5{\rm \pi}}}( \mathrm{i} {\mathsf{Q}}^{({3})}_{030} - 3 {\mathsf{Q}}^{({3})}_{120} - 3\mathrm{i} {\mathsf{Q}}^{({3})}_{210} + {\mathsf{Q}}^{({3})}_{300}) \nonumber\\ & = \frac{1}{24}\sqrt{\frac{7}{5{\rm \pi}}}( \mathrm{i} {\mathsf{Q}}^{({3})}_{222} - 3 {\mathsf{Q}}^{({3})}_{122} - 3\mathrm{i} {\mathsf{Q}}^{({3})}_{112} + {\mathsf{Q}}^{({3})}_{111}) . \end{align}

In the last equation, we used the conventional tensor indices, namely ${\mathsf{Q}}^{({3})}_{ijk}$.

4. Basis transformation

As already pointed out, we called the expression derived in (3.31) a preliminary basis transformation, because the sum on the right-hand side contained multipole functions with index $p > 1$. However, only the multipole functions with $p \leq ~1$ are a basis of the space of homogeneous and harmonic polynomials of degree $l$ (see (2.18)). To turn (3.31) into an actual basis transformation, we have to express all ${\mathsf{M}}^{({l})}_{pqr}$ with $p > 1$ as a linear combination of the multipole basis functions.

We can accomplish this by exploiting that the object $\boldsymbol{\mathsf{M}}^{({l})}$ is symmetric and traceless (cf. (3.19)). In a first step, we use the notation for the indices of $\boldsymbol{\mathsf{M}}^{({l})}$, which we introduced in (2.17), to express the latter property of $\boldsymbol{\mathsf{M}}^{({l})}$ in the following way:

(4.1)\begin{equation} {\mathsf{M}}^{({l})}_{pqr} + {\mathsf{M}}^{({l})}_{(p-2)(q+2)r} + {\mathsf{M}}^{({l})}_{(p-2)q(r+2)} = 0, \end{equation}

with $p > 1$ and, by definition, $p + q + r = l$. This becomes plausible when an example is considered, e.g. $l = 5$ and $p = 2, q = 1$ and $r = 2$,

(4.2)\begin{equation} 0 = {\mathsf{M}}^{({5})}_{212} + {\mathsf{M}}^{({5})}_{032} + {\mathsf{M}}^{({5})}_{014} = {\mathsf{M}}^{({5})}_{11233} + {\mathsf{M}}^{({5})}_{22233} + {\mathsf{M}}^{({5})}_{33233} = {\mathsf{M}}^{({5})}_{ii233} . \end{equation}

Note that we used the symmetry of ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$. An immediate consequence of (4.1) is that we can express a multipole function with index $p > 1$ as a sum of multipole functions with index $p - 2$. These, in turn, can again be expressed as a sum of multipole functions with index $p - 4$. Depending on whether $p$ is even or odd, this chain ends when $p = 0$ or $p = 1$. This computation can be represented in the form of Pascal's triangle (see figure 1). Note the alternating minus sign. Or, more compactly, in the following formula:

(4.3)\begin{equation} {\mathsf{M}}^{({l})}_{pqr} = ({-}1)^{\lfloor p/2 \rfloor} \sum^{\lfloor p/2 \rfloor}_{k = 0} \left(\begin{array}{c}\left\lfloor\dfrac{p}{2}\right\rfloor \\ k \end{array}\right) \begin{cases} {\mathsf{M}}^{({l})}_{0(q + p -2k)(r+2k)} & p \text{ is even} , \\ {\mathsf{M}}^{({l})}_{1(q + p -(2k +1))(r + 2k)} & p \text{ is odd}. \end{cases} \end{equation}

Figure 1. Dependence of the multipole functions with $p > 1$ on functions with smaller $p$ given by Pascal's triangle.

We can split the sum in (3.31) into even and odd $p$ values and plug in the expression for the multipole functions, which we just derived. This yields

(4.4)\begin{align} r^{l}Y^{m}_{l} & = a^{m}_{l} \left[ \sum^{a}_{n=0} \mathrm{i}^{m - 2n}\binom{m}{2n} {\mathsf{M}}^{({l})}_{2n(m - 2n)(l - m)} \right. \nonumber\\ & \quad \left. + \sum^{b}_{n=0} \mathrm{i}^{m - (2n + 1)}\binom{m}{2n + 1}{\mathsf{M}}^{({l})}_{(2n + 1)(m - (2n + 1))(l - m)} \right] \nonumber\\ & = a^{m}_{l} \left[ \sum^{a}_{n=0}\sum^{n}_{k = 0} \mathrm{i}^{m}\binom{m}{2n} \binom{n}{k} {\mathsf{M}}^{({l})}_{0(m - 2k)(l - m +2k)} \right. \nonumber\\ & \quad \left. - \sum^{b}_{n=0}\sum^{n}_{k = 0} \mathrm{i}^{m}\binom{m}{2n + 1} \binom{n}{k}{\mathsf{M}}^{({l})}_{1(m -(2k +1))(l - m + 2k)} \right]. \end{align}

The limits of the sums are

(4.5a,b)\begin{equation} a = \lfloor m/2 \rfloor \quad \text{and}\quad b = \lfloor (m - 1)/2 \rfloor . \end{equation}

In a last step, we note that we can exchange the two sums. For example, for even $p$, we find

(4.6)\begin{equation} \sum^{a}_{n=0}\sum^{n}_{k = 0} \mathrm{i}^{m} \binom{m}{2n} \binom{n}{k} {\mathsf{M}}^{({l})}_{0(m - 2k)(l - m +2k)} = \sum^{a}_{k = 0} \left(\mathrm{i}^{m} \sum^{a}_{n=k}\binom{m}{2n} \binom{n}{k}\right) {\mathsf{M}}^{({l})}_{0(m - 2k)(l - m + 2k)} . \end{equation}

The same holds true for odd $p$. We arrive at the final result of this section, namely

(4.7)\begin{equation} r^{l}Y^{m}_{l} = \sum^{a}_{k = 0} \beta^{m}_{l,0,k} {\mathsf{M}}^{({l})}_{0(m - 2k)(l - m +2k)} + \sum^{b}_{k = 0} \beta^{m}_{l,1,k} {\mathsf{M}}^{({l})}_{1(m - (2k + 1))(l - m + 2k)} . \end{equation}

For $m \geq ~0$ and with the coefficients

(4.8a,b)\begin{equation} \beta^{m}_{l,0,k} \equiv \mathrm{i}^{m} a^{m}_{l} \sum^{a}_{n=k}\binom{m}{2n} \binom{n}{k} \quad\text{and}\quad \beta^{m}_{l,1,k} \equiv{-}\mathrm{i}^{m + 1} a^{m}_{l} \sum^{b}_{n=k}\binom{m}{2n + 1} \binom{n}{k} , \end{equation}

the solid harmonics with $m < 0$ can be computed with the formula $Y^{-m}_{l} = (-1)^{m}{Y^{m}_{l}}^{*}$.

We thus have a way to express an arbitrary solid harmonic of degree $l$ and order $m$ as a linear combination of the multipole basis functions with degree $l$. We call the coefficients $\beta ^{m}_{l,0,k}$ and $\beta ^{m}_{l,1,k}$ the basis transformation between the solid harmonics and the multipole basis functions in the space of homogeneous and harmonic polynomials of degree $l$.

Last, but not least, the complex conjugate of the basis transformation relates the spherical multipole moments with the Cartesian multipole moments. A practical implication of this section's considerations is that it is only necessary to compute $2l + 1$ components of the Cartesian multipole moments (those corresponding to the multipole basis functions given in (2.18)) instead of the $(l + 2)(l + 1)/2$ components of the symmetric tensor $\boldsymbol{\mathsf{Q}}^{({l})}$. A way to reconstruct the dependent components of $\boldsymbol{\mathsf{M}}^{({l})}$ (or $\boldsymbol{\mathsf{Q}}^{({l})}$) will be found in the next section while deriving the inverse basis transformation.

5. Inverse transformation

We derive the inverse of the basis transformation given in (4.7) by, first, collecting the coefficients $\beta ^{m}_{l,0,k}$ and $\beta ^{m}_{l,1,k}$ in a matrix, the basis transformation matrix $\boldsymbol{\mathsf{B}}$, and second, by inverting this matrix.

How difficult it is to invert a matrix depends on its structure, e.g. if it is a diagonal matrix, computing the inverse consists in taking the reciprocal of its elements. The structure of the matrix $\boldsymbol{\mathsf{B}}$ simplifies if we use real spherical harmonics instead of Laplace's spherical harmonics.

The real spherical harmonics are defined as

(5.1)\begin{equation} \begin{pmatrix} Y_{lm0} \\ Y_{l00} \\ Y_{lm1} \end{pmatrix} \equiv \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 0 & ({-}1)^{m} \\ 0 & \sqrt{2} & 0 \\ -\mathrm{i} & 0 & ({-}1)^{m}\mathrm{i} \end{pmatrix} \begin{pmatrix} Y^{m}_{l} \\ Y^{0}_{l} \\ Y^{{-}m}_{l} \end{pmatrix} . \end{equation}

Henceforth, we denote the above matrix by $\boldsymbol{\mathsf{S}}$. Note that $\boldsymbol{\mathsf{S}}$ is unitary, i.e. $\boldsymbol{\mathsf{S}}^{{\dagger} }\boldsymbol{\mathsf{S}} = \boldsymbol{\mathsf{1}}$. An explicit expression for the real spherical harmonics is

(5.2)\begin{equation} Y_{lms} = N_{lm} P^{m}_{l}(\cos\theta) (\delta_{s0}\cos{m \varphi} + \delta_{s1}\sin{m \varphi}), \end{equation}

with $m \geq 0$. Furthermore, the above definition implies the normalisation $N_{lm} \equiv \sqrt {2 - \delta _{m0}} N^{m}_{l}$.

Applying the unitary transformation to (4.7) results in the following pattern:

(5.3)\begin{equation} \left.\begin{array}{ll@{}} s = 0: & \\ \text{even } m: & \displaystyle r^{l}Y_{lm0} = \sqrt{2} \sum^{a}_{k = 0} \beta^{m}_{l,0,k} {\mathsf{M}}^{({l})}_{0(m - 2k)(l - m + 2k)} \\ \text{odd } m: & \displaystyle r^{l}Y_{lm0} = \sqrt{2} \sum^{b}_{k = 0} \beta^{m}_{l,1,k} {\mathsf{M}}^{({l})}_{1(m - (2k + 1))(l - m + 2k))}\\ s = 1: & \\ \text{even } m: & \displaystyle r^{l}Y_{lm1} ={-}\sqrt{2}\mathrm{i} \sum^{b}_{k = 0} \beta^{m}_{l,1,k} {\mathsf{M}}^{({l})}_{1(m - (2k + 1))(l - m + 2k)} \\ \text{odd } m: & \displaystyle r^{l}Y_{lm1} ={-}\sqrt{2}\mathrm{i} \sum^{a}_{k = 0} \beta^{m}_{l,0,k} {\mathsf{M}}^{({l})}_{0(m - 2k)(l - m + 2k)} \end{array}\right\}. \end{equation}

This pattern has two implications for the structure of the basis transformation matrix $\boldsymbol{\mathsf{B}}$. First, the unitary transformation brings it into a form in which it consists of four blocks which correspond to the above four cases. Second, since the limits of the sums $a$ and $b$ decrease with decreasing $m$, the corresponding real solid harmonic depends on less multipole basis functions. For example, let us consider the case $s = 0$ and even $m$ for $l = 5$. The even values of $m$ are zero, two and four. The corresponding limits $a$ are zero, one and two,respectively. Hence, the real solid harmonic of degree five and order zero depends on one multipole basis function, the one with order two on two and the last one with order four on three multipole basis functions. This pattern translates into a triangular matrix block. Hence, the basis transformation matrix can be transformed into a matrix, which consists of four triangular matrix blocks.

Whether the triangular matrices show up also depends on the ordering of the solid harmonics and the multipole basis functions. For the solid harmonics, we pick the ordering, which is implicitly suggested by (5.1), namely

(5.4)\begin{equation} r^{l}\boldsymbol{Y}^{l} \equiv r^{l} \begin{pmatrix} Y^{l}_{l} & Y^{l - 1}_{l} & \cdots & Y^{{-}l + 1}_{l} & Y^{{-}l}_{l} \end{pmatrix}^{T} . \end{equation}

Additionally, a corresponding ordering for the multipole basis functions, for example, is

(5.5)\begin{equation} \boldsymbol{M}^{l} \equiv \begin{pmatrix} {\mathsf{M}}^{({l})}_{0l0} & {\mathsf{M}}^{({l})}_{0(l-1)1} & \cdots & {\mathsf{M}}^{({l})}_{1(l-2)1} & {\mathsf{M}}^{({l})}_{1(l-1)0} \end{pmatrix}^{T} . \end{equation}

Assembling the matrix $\boldsymbol{\mathsf{B}}$ in accordance with these orderings brings (4.7) into the compact form:

(5.6)\begin{equation} r^{l}\boldsymbol{Y}^{l}(\theta, \varphi) = \boldsymbol{\mathsf{B}}\boldsymbol{M}^{l}(\boldsymbol{r}). \end{equation}

This equation can be multiplied with the unitary matrix $\boldsymbol{\mathsf{S}}$. On the left-hand side, the real solid harmonics will show up and on the right-hand side, a matrix with the four blocks, which we found in (5.3). If we now re-order the real solid harmonics (and the multipole basis functions) such that the harmonics with even order $m$ (and odd order $m$) are grouped together, four triangular matrices result. If we additionally reverse the order of the real solid harmonics with $s = 1$, we get four upper triangular matrices. Let $\boldsymbol{\mathsf{P}}$ denote the corresponding permutation matrix, then

(5.7)\begin{equation} r^{l}\boldsymbol{\mathsf{PS}} \boldsymbol{Y}^{l} = \boldsymbol{\mathsf{P}} \boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{P}}^{T} \boldsymbol{\mathsf{P}} \boldsymbol{M}^{l} . \end{equation}

An actual computation for even $l$ shows the following matrix structure:

(5.8)\begin{equation} \tilde{\boldsymbol{\mathsf{B}}} \equiv \boldsymbol{\mathsf{P}} \boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{P}}^{T} = \begin{pmatrix} \boldsymbol{\mathsf{U}}_{1} & & & \\ & & & \boldsymbol{\mathsf{U}}_{2}\\ & & \boldsymbol{\mathsf{U}}_{3} & \\ & \boldsymbol{\mathsf{U}}_{4} & & \end{pmatrix} . \end{equation}

For odd $l$, the four upper triangular matrices are at other positions in the matrix.

The inversion of this matrix is equivalent to inverting four upper triangular matrices. This can, for example, be done by applying repeatedly back-substitution with unit vectors as the right-hand sides (Press et al. 2007, § 2.2.1). We define the inverse basis transformation matrix as

(5.9)\begin{equation} \boldsymbol{\mathsf{A}} \equiv \boldsymbol{\mathsf{B}}^{{-}1} = \boldsymbol{\mathsf{P}}^{T} \tilde{\boldsymbol{\mathsf{B}}}^{{-}1} \boldsymbol{\mathsf{P}} \boldsymbol{\mathsf{S}} . \end{equation}

Thus, $\boldsymbol {M}^{l} = r^{l} \boldsymbol{\mathsf{A}} \boldsymbol {Y}^{l}$ and with the help of a suitable mapping between the matrix indices and the $p,q,r$ indices, we get

(5.10a,b)\begin{equation} {\mathsf{M}}^{({l})}_{0qr} = \sum^{l}_{m ={-}l} \alpha^{m}_{l,0,q,r} r^{l}Y^{m}_{l} \quad\text{and}\quad {\mathsf{M}}^{({l})}_{1qr} = \sum^{l}_{m ={-}l} \alpha^{m}_{l,1,q,r} r^{l}Y^{m}_{l} . \end{equation}

We call the $\alpha ^{m}_{l,0,q,r}$ and $\alpha ^{m}_{l,1,q,r}$ the inverse basis transformation.

The dependent multipole functions (i.e. the functions with index $p > 1$) can be expressed as a linear combination of the multipole basis functions (see (4.3)). Since the multipole basis functions are $\boldsymbol {M}^{l} = r^{l}\boldsymbol{\mathsf{A}}\boldsymbol {Y}^{l}$, it is possible to obtain an expression for the dependent multipole functions in terms of the solid harmonics by adding the rows of the inverse basis transformation matrix while multiplying them by the factor given in (4.3).

The relation between the (independent) Cartesian multipole moments and the spherical multipole moments is given by the complex conjugate of the inverse basis transformation (matrix), i.e.

(5.11a,b)\begin{equation} \boldsymbol{q}^{l} = \boldsymbol{\mathsf{B}}^{*} \boldsymbol{Q}^{l} \quad\text{and}\quad \boldsymbol{Q}^{l} =\boldsymbol{\mathsf{A}}^{*} \boldsymbol{q}^{l} . \end{equation}

Last, but not least, instead of inverting the basis transformation matrix $\boldsymbol{\mathsf{B}}$, it is possible to directly compute the inverse basis transformation matrix $\boldsymbol{\mathsf{A}}$. A derivation of the formulae for the $\alpha$s is provided in Appendix C.

6. The Vlasov–Fokker–Planck equation

In this section, we come back to the original problem, which stood at the beginning of our investigation of the relation between the Cartesian multipole moments and the spherical multipole moments: the two expansions of the distribution function $f$ and their relation.

In the first step, we show that their relation is the same as the relation between the Cartesian multipole moments and spherical multipole moments. To this end, we investigate the differences between (2.9) and (2.21) and the expansions of $f$ given in (1.3) and (1.2). The expansions of $f$ are done in momentum space (and not in configuration space), the coefficients (previously ${\mathsf{Q}}^{({l})}_{i_{1} \ldots i_{l}}$ and $q^{m}_{l}$ and now ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ and $f^{m}_{l}$) are functions, there are no numerical factors and, most importantly, the magnitude of $\boldsymbol {p}$ (previously $\boldsymbol {r}$) was set to one. This last difference accounts for the factor $p^{-l}$ in (1.3), since $\boldsymbol {p}/p$ restricts the components of $\boldsymbol {p}$ to the unit sphere in momentum space. The dependence on the magnitude of $p$ is now included in the coefficients $f^{m}_{l}$ and not explicitly given as $1/r^{2~l + 1}$.

The fact that the coefficients ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ and $f^{m}_{l}$ are functions does not change the relation between them. Actually, the object ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ is constructed (for example, it is chosen to be symmetric in its indices and it should be chosen to be traceless as well) such that the relation between it and the $f^{m}_{l}$s is the same as between the Cartesian multipole moments and the spherical multipole moments. The fact that $p$ (or $r$) was set to one also does not change the relation. This can be seen by setting $r$ to one and restricting the magnitude of $\boldsymbol {r}$ to one in (5.6). However, that there are no factors means that $4{\rm \pi} /(2~l + 1)$ was included in the definition of $f^{m}_{l}$ (see (2.22)) and $1/l!$ is included in the definition of ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$ (cf. with (2.7)). Thus, the relation is given by (5.11a,b), i.e.

(6.1)\begin{equation} \boldsymbol{f}^{l} = \boldsymbol{\mathsf{D}}_{1}\boldsymbol{\mathsf{B}}^{*}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}^{l} , \end{equation}

where we replaced $\boldsymbol {q}^{l}$ with $\boldsymbol {f}^{l}$, $\boldsymbol {Q}^{l}$ with $\boldsymbol {F}^{l}$, and multiplied with $\boldsymbol{\mathsf{D}}_{1} \equiv 4{\rm \pi} /(2l + 1) {1}$ and $\boldsymbol{\mathsf{D}}_{2} \equiv l! {1}$ to include the factors in the definitions of the coefficients. The inverse transformation is obtained by computing the inverse of the above matrix, as explained in § 5.

Johnston (1960) derives this inverse basis transformation for $l = 2$ and $l = 3$. However, when expanding $f$, he uses the real spherical harmonics without normalisation $N_{lm}$ and without including $(-1)^{m}$ in the definition of the associated Legendre polynomials. Thus, he implicitly starts with the following equation:

(6.2)\begin{align} f(t,\boldsymbol{r}, \boldsymbol{p}) & = \sum^{\infty}_{l = 0} \sum^{l}_{ m = 0} \sum^{1}_{s = 0} ({-}1)^{m}\frac{2}{1 + \delta_{0m}} \frac{(l-m)!}{(l + m)!} \frac{\bar{f}_{lms}(t, \boldsymbol{r}, p)}{N_{lm}}Y_{lms}(\theta, \varphi)\nonumber\\ & \equiv \sum^{\infty}_{l = 0} \sum^{l}_{ m = 0} \sum^{1}_{s = 0} f_{lms}(t, \boldsymbol{r}, p)Y_{lms}(\theta, \varphi). \end{align}

The fraction with factorials appears because the numerical factor in the addition theorem changes if (real) spherical harmonics without normalisation are used to derive it (cf. with (2.20)). Moreover, $1/N_{lm}$ removes the normalisation. Additionally, in the last definition of (6.2), all these factors are included in the coefficients $f_{lms}$.

Thus, we should be able to reproduce his results by inverting the matrix on the right-hand side of

(6.3)\begin{equation} \boldsymbol{f}^{l}_{{J}} = \boldsymbol{\mathsf{C}} \boldsymbol{\mathsf{D}}_{3} \boldsymbol{\mathsf{N}}^{{-}1} \boldsymbol{\mathsf{S}}^{*} \boldsymbol{\mathsf{B}}^{*}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}^{l} , \end{equation}

which yields the following matrix:

(6.4)\begin{equation} \boldsymbol{\mathsf{A}}_{{J}} = \boldsymbol{\mathsf{D}}^{{-}1}_{2} \boldsymbol{\mathsf{A}}^{*}\boldsymbol{\mathsf{S}}^{T}\boldsymbol{\mathsf{D}}^{{-}1}_{3}\boldsymbol{\mathsf{C}} , \end{equation}

where $\boldsymbol{\mathsf{A}}^{*} = {\boldsymbol{\mathsf{B}}^{*}}^{-1}$ and where we introduced the following definitions:

(6.5)\begin{equation} \left.\begin{array}{c@{}} \boldsymbol{f}^{l}_{{J}} \equiv (f_{ll0} f_{l(l-1)0} \ldots f_{l(l-1)1} f_{ll1} ) \\ \boldsymbol{\mathsf{N}} \equiv \operatorname{diag}(N_{ll} N_{l(l-1)}\ldots N_{l(l-1)} N_{ll})\\ \boldsymbol{\mathsf{C}} \equiv \operatorname{diag}(({-}1)^{l} ({-}1)^{l-1} \ldots ({-}1)^{l-1} ({-}1)^{l}) \quad \text{and}\\ \boldsymbol{\mathsf{D}}_{3} \equiv \operatorname{diag}\left(\dfrac{2}{(2l)!} \dfrac{2}{(2l - 1)!} \ldots 1 \ldots \dfrac{2}{(2l - 1)!} \dfrac{2}{(2l)!}\right) \end{array}\right\}. \end{equation}

Furthermore, $\boldsymbol{\mathsf{S}}$ is the unitary transformation between Laplace's spherical harmonics and the real spherical harmonics (see (5.2)). The complex conjugate is necessary, because we are using the complex conjugate of (5.6) (namely, we are using ${Y^{m}_{l}}^{*}$ and comparing with (5.11a,b)). The dependent components of ${\mathsf{F}}^{l}_{i_{1} \ldots i_{l}}$ can be computed as described at the end of § 5. The entries of the matrix $\boldsymbol{\mathsf{A}}_{{J}}$ can be compared with (9a–d) in the paper by Johnston (1960). Moreover, the entries of $\boldsymbol{\mathsf{A}}^{*}$ can be directly computed with formulae derived in Appendix C.

Last, but not least, the Vlasov–Fokker–Planck equation and the corresponding expansions of $f$ can be considered as an example partial differential equation (PDE). We note that whenever (Laplace's) spherical harmonics are used to expand a solution to a PDE, it is possible to work with the multipole functions instead (or vice versa). If the multipole functions are used, it is always enough to compute the $2l + 1$ independent components, namely the multipole basis functions.

7. Conclusion

We started off with the aim to find a way to convert between the coefficients appearing in the Cartesian tensor expansion and the spherical harmonic expansion of the distribution function $f$. A problem which was prominently investigated by Johnston (1960).

We pointed out that the relation between the Cartesian tensors $\boldsymbol{\mathsf{F}}^{({l})}$ and the $f^{m}_{l}$ is the same as the relation between the Cartesian multipole moments $\boldsymbol{\mathsf{Q}}^{({l})}$ and the spherical multipole moments $q^{m}_{l}$. Since multipole moments are ubiquitous in physics, we revisited the problem of converting between coefficients of expansions in the context of the multipole expansion.

In the first step, we introduced a new definition of the Cartesian multipole moments and showed that it reproduces the known multipole moments. We saw that this definition contained the multipole functions ${\mathsf{M}}^{({l})}_{i_{1} \ldots i_{l}}$, which are proven to be homogeneous and harmonic polynomials of degree $l$. Furthermore, we pointed out that a subset of the multipole functions is a basis of the space of homogeneous and harmonic polynomials $\mathcal {H}^{l}(\mathbb {R}^{3})$. We called this subset the multipole basis functions.

Subsequently, we took a look at the definition of the spherical multipole moments and noticed that it contains the complex conjugate of a solid harmonic of degree $l$ and order $m$. We showed that these functions are also harmonic and homogeneous polynomials of degree $l$ and that they are a basis of $\mathcal {H}^{l}(\mathbb {R}^{3})$ as well.

This led us to the new conclusion that the relation between the Cartesian multipole moments and the spherical multipole moments can be formalised as a basis transformation between the solid harmonics of degree $l$ and the multipole basis functions of degree $l$.

We derived closed-form formulae to compute the basis transformation, i.e. we derived a way to express the solid harmonics of degree $l$ and order $m \geq 0$ as a sum of multipole basis functions. We inverted the basis transformation matrix and simplified the computation of the inverse matrix considerably. Additionally, we also derived closed-form formulae for the inverse basis transformation (see Appendix C). We, thus, provided a systematic and easy-to-compute method to convert between the Cartesian multipole moments and the spherical multipole moments.

A key practical implication of our approach is the insight, that it is enough to compute $2\,l + 1$ components of the Cartesian multipole moments, namely the components corresponding to the multipole basis functions. The dependent components can be computed as explained at the end of § 5.

Eventually, we came back to our original problem and applied our knowledge about the relation between the Cartesian multipole moments and the spherical multipole moments in the context of the Cartesian tensor and the spherical harmonic expansion of $f$. We used the basis transformation matrix $\boldsymbol{\mathsf{B}}$ to express the coefficients $\,f^{m}_{l}$ as a linear combination of the independent components of the corresponding Cartesian tensor ${\mathsf{F}}^{({l})}_{i_{1} \ldots i_{l}}$. We explained how our results reproduce the findings in Johnston (1960).

We developed a free and open-source command-line tool, called multipole-conv (Schween & Reville 2022), which is able to compute the basis transformation matrix and its inverse. We hope that this tool will be useful for all who work with spherical harmonics or multipole moments.

Last, but not least, we would like to assess our approach critically. First, the elements in the basis transformation matrix $\boldsymbol{\mathsf{B}}$ become very small even for small $l$, say $l \approx 15$, and the elements of the inverse basis transformation matrix $\boldsymbol{\mathsf{A}}$ very large. This leads to numerical errors in the computation of $\boldsymbol{\mathsf{A}}$. We suspect that this happens because Laplace's spherical harmonics are normalised whereas the multipole basis functions are not. Second, the mathematically inclined reader may know that Laplace's spherical harmonics are an irreducible representation of $SO(3)$. Since both definitions of the multipole basis functions, i.e. the definition with the Kelvin transform (see (2.17) and, in particular, (3.9)) and the definition with Efimov's ladder operator (cf. with (3.7)) generalise to arbitrary dimensions, it could be that there also exists a basis transformation between these generalised multipole basis functions and irreducible representations of $SO(n)$.