12.3. Multipole and basis-function expansions¶
In the previous section, we demonstrated that it is in fact quite straightforward to compute the gravitational potential, field, and rotation curve for any spheroidal or ellipsoidal mass distribution. However, as we saw in Section 12.1, many elliptical galaxies and likely all dark-matter halos have shapes and/or orientations that change with radius. To model such systems, we need more general methods for solving the Poisson equation than we have discussed so far. We already discussed one such solution in Chapter 7.3.5 that works best for disky mass distributions, but it involves difficult-to-work with Bessel functions. As such, this solution is not that useful in practice except for models with additional symmetries or constraints that simplify the expression substantially. In this section, we will look at a few different ways to solve the Poisson equation in full generality that are also useful in practice. For full numerical N-body work, we will discuss the solution of the Poisson equation in Section 12.4.1 below.
Because the Poisson equation (2.2) is linear, a general approach—which we have already taken many times!—is to break up the density into a (possibly infinite) sum of additive contributions, solve for the potential of each contribution, and add up all of the resulting potentials to obtain the full gravitational potential. We followed this approach in the previous section to determine the potential of an oblate spheroid by decomposing it into an infinite sequence of thin oblate shells. But the density constituents do not have to be simple shells, we simply want the density decomposition to be general and easy to compute, and for the Poisson equation to be easily solvable for each density contribution. For example, in Chapter 7.3.2, we determined the gravitational potential of an axisymmetric, razor-thin disk by breaking the surface density up into contributions proportional to \(k\,J_0(kR)\), which is general because of the Fourier-Bessel theorem.
In all of the methods that we discuss in this section, we will obtain the gravitational potential \(\Phi(\vec{x})\) by expanding the density \(\rho(\vec{x})\) into a complete set of basis functions \(\rho_{\vec{n}}(\vec{x})\) as \begin{equation}\label{eq-dens-basis-expansion} \rho(\vec{x}) = \sum_{\vec{n}} A_{\vec{n}}\,\rho_{\vec{n}}(\vec{x})\,, \end{equation} where each basis function has a solution \(\Phi_{\vec{n}}(\vec{x})\) of the Poisson equation \begin{equation} \nabla^2 \Phi_{\vec{n}}(\vec{x}) = 4\pi\,G\,\rho_{\vec{n}}(\vec{x})\,. \end{equation} Using the linearity of the Poisson equation, we then obtain the gravitational potential corresponding to \(\rho(\vec{x})\) as \begin{equation} \Phi(\vec{x}) = \sum_{\vec{n}} A_{\vec{n}}\,\Phi_{\vec{n}}(\vec{x})\,. \end{equation} Note that while we write a summation sign here, one or more of the dimensions may be integrated over instead of summed if a continuous basis-function set is used. For example, in the multipole-expansion method that we discuss in Section 12.3.1 below, the radial part of the density and potential is an integral over infinitesimally-thin shells.
The Poisson equation is a partial differential equation that is particularly suited to solving by a basis function approach, because:
The Laplace operator \(\nabla^2\) is Hermitian (see below for proof); this implies that eigenfunctions of the operator—functions \(f(\vec{x})\) for which \(\nabla^2 f(\vec{x}) \propto f(\vec{x})\) analogous to eigenvectors in linear algebra—are orthogonal (see proof below). Typically, the eigenfunctions are a manifestly complete basis set for the space of functions under consideration—smooth functions that are well-behaved on the boundary of the volume being considered. In the methods described in this section, completeness is guaranteed by Sturm-Liouville theory.
The Poisson equation can be solved using the separation-of-variables technique in a variety of different coordinate systems (e.g., Cartesian, cylindrical, spherical, oblate/prolate spheroidal, ellipsoidal; see Section 12.2). This provides a convenient way to obtain eigenfunctions of the Laplace operator (or other orthogonal sets of complete basis functions).
To prove the statements in the first bullet point above, we define the inner product \(\langle f,g\rangle\) of two functions \(f\) and \(g\) on \(\mathbb{R}^k\) as \(\langle f,g\rangle = \int\mathrm{d} \vec{x}\, f^*(\vec{x})\,g(\vec{x})\), where \(f^*(\vec{x})\) denotes the complex conjugate of \(f(\vec{x})\). We can then show that the Laplacian is Hermitian (at least if the functions and their gradients go to zero at infinity), that is, that \(\langle f,\nabla^2 g\rangle = \langle \nabla^2 f, g\rangle\), because \begin{align} \langle f,\nabla^2 g\rangle & = \phantom{-}\int\mathrm{d} \vec{x}\, f^*(\vec{x})\,\nabla^2 g(\vec{x})= \int\mathrm{d}\vec{x}\, f^*(\vec{x})\,\mathrm{div} \,\vec{\nabla} g(\vec{x})\nonumber\\ & = -\int\mathrm{d} \vec{x}\, \vec{\nabla} f^*(\vec{x})\,\vec{\nabla} g(\vec{x})\\ & = \phantom{-}\int\mathrm{d} \vec{x}\, \mathrm{div}\vec{\nabla} f^*(\vec{x})\,g(\vec{x}) = \int\mathrm{d} \vec{x} \,[\nabla^2 f(\vec{x})]^*\,g(\vec{x})= \langle \nabla^2 f, g\rangle\nonumber\,, \end{align} where we have used the definition of the Laplacian \(\nabla^2 = \mathrm{div}\,\vec{\nabla}\) twice and used multi-dimensional integration by parts to go to the third and fourth lines. Good sets of basis functions can then be found as eigenfunctions of the Laplacian, that is, functions \(f_{\vec{n}}(\vec{x})\) that satisfy \(\nabla^2 f_{\vec{n}} = \lambda_{\vec{n}}\,f_{\vec{n}}\) (this equation is also known as the Helmholtz equation). Because the Laplacian is Hermitian, these basis functions are orthogonal, since then \begin{align} \langle f_{\vec{n}},g_{\vec{n}'}\rangle & = \frac{1}{\lambda_{\vec{n}'} }\,\langle f_{\vec{n}},\nabla^2 g_{\vec{n}'}\rangle = \frac{1}{\lambda_{\vec{n}'} }\,\langle \nabla^2 f_{\vec{n}}, g_{\vec{n}'}\rangle= \frac{\lambda_{\vec{n}}}{\lambda_{\vec{n}'} }\,\langle f_{\vec{n}}, g_{\vec{n}'}\rangle\,, \end{align} and therefore \(\langle f_{\vec{n}},g_{\vec{n}'}\rangle =0\) if \(\lambda_{\vec{n}} \neq\lambda_{\vec{n}'}\); if an eigenvalue has a multiplicity larger than one, then the eigenfunctions for that eigenvalue can be orthogonalized using Gram-Schmidt orthogonalization (see Equation B.68 in Appendix B.3.6). We can always define the eigenfunctions such that \(\langle f_{\vec{n}},f_{\vec{n}}\rangle = 1\) and the eigenfunctions are then orthonormal; let’s denote them as \(\rho_{\vec{n}}\) in that case. Eigenfunctions provide a convenient basis set, because the density can then be decomposed in the form of Equation (12.52) by computing \begin{equation} A_{\vec{n}} = \langle \rho_{\vec{n}},\rho \rangle = \int \mathrm{d}^3\vec{x}\,\rho^*_{\vec{n}}(\vec{x})\,\rho(\vec{x})\,. \end{equation} To demonstrate that the set of eigenfunctions forms a complete basis function set requires details about the ordinary differential equations obtained using the separation-of-variables approach. Often, these differential equations will be of the Sturm-Liouville type and we can apply Sturm-Liouville theory (Courant & Hilbert 1953; Zettl 2010): this theory states that the eigenfunctions of Sturm-Liouville-type differential equations form a complete basis.
12.3.1. Multipole expansions¶
As a first application of the basis-function-expansion approach, we work out the standard multipole-expansion method of decomposing any density field into an infinite sequence of spherical shells with non-uniform densities. It’s obvious that any density distribution can be decomposed in this way and for elliptical galaxies and dark matter halos that are not too far from spherical it works quite well. Similar to how we proceeded in the Section 12.2, we first solve for the potential of a single thin, non-uniform shell and then compute the full potential by integrating over all shells. Everywhere outside of the shell, the Poisson equation reduces to the Laplace equation \(\nabla^2 \Phi = 0\) and we can therefore again obtain the full potential by solving the Laplace equation inside and outside of the shell in question and matching them at the radius of the shell using the appropriate boundary conditions.
The thin spherical shell with non-uniform density has a density in spherical coordinates \(\rho(r,\phi,\theta) = \delta(r-a)\,\sigma(\phi,\theta)\), where \(\sigma(\phi,\theta)\) is the surface density of the shell. To solve the Poisson equation for this density, we can follow the approach introduced above, but on the sphere rather than in three dimensions: we decompose \(\sigma(\phi,\theta)\) into basis functions, which we will suggestively denote as \(Y_l^m(\theta,\phi)\) and then solve the Poisson equation \(\nabla^2 \Phi = 4\pi G\,\sigma_{lm}\,\delta(r-a)\,Y_l^m(\theta,\phi)\) for each basis-function contribution, where \(\sigma_{lm}\) are expansion coefficients. To obtain a suitable set of basis functions \(Y_l^m(\theta,\phi)\) on the sphere, we look for eigenfunctions of the Laplace operator on the sphere using the separation-of-variables approach, that is functions \(\Phi(\phi,\theta) = F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\) that satisfy (using the Laplacian from Equation A.8 without the radial first term and with \(r=1\) in the other terms) \begin{equation}\label{eq-poisson-separate-spherical-zerothstep} \nabla^2 \Phi = \frac{1}{\sin\theta}\,\frac{\partial }{\partial \theta}\left(\sin\theta\,\frac{\partial \Phi}{\partial \theta}\right) +\frac{1}{\sin^2\theta}\,\frac{\partial^2 \Phi}{\partial \phi^2} = \lambda_{\vec{n}}\Phi\,, \end{equation} where \(\vec{n}\) is a vector indexing the eigenfunctions that we will determine as part of the solution. Substituting \(\Phi(\phi,\theta) = F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\) and bringing all of the \(\phi\) dependence to the right-hand side of the equation, we find \begin{equation}\label{eq-poisson-separate-spherical-firststep} \frac{\sin\theta}{T_{\vec{n}}(\theta)}\,\frac{\mathrm{d}}{\mathrm{d} \theta}\left(\sin\theta\,\frac{\mathrm{d} T_{\vec{n}}(\theta)}{\mathrm{d} \theta}\right) - \lambda_{\vec{n}}\,\sin^2\theta = - \frac{1}{F_{\vec{n}}(\phi)}\,\frac{\mathrm{d}^2 F_{\vec{n}}(\phi)}{\mathrm{d} \phi^2}\,. \end{equation} Because the left-hand side only depends on \(\theta\) and the right-hand side only depends on \(\phi\), both sides must equal a constant, say \(m^2\). The right-hand side then leads to the following equation \begin{equation}\label{eq-triaxialgrav-sphereFphi-diffeq} \frac{\mathrm{d}^2 F_{\vec{n}}(\phi)}{\mathrm{d} \phi^2} = -m^2\,F_{\vec{n}}(\phi)\,, \end{equation} with solutions \begin{equation}\label{eq-triaxialgrav-sphereFphi-solution} F_{\vec{n}}(\phi) = F_m\,e^{im\phi}\,. \end{equation} Because \(F_{\vec{n}}(\phi+2\pi)= F_{\vec{n}}(\phi)\), we have that \(e^{i2\pi\,m} = 1\), and therefore \(m\) needs to be an integer: \(m = 0,\pm1,\pm2,\ldots\).
The left-hand side of Equation (12.59) becomes \begin{equation}\label{eq-triaxialgrav-sphereTtheta-diffeq} \frac{m^2}{\sin^2\theta}\,T_{\vec{n}}(\theta)-\frac{1}{\sin\theta}\,\frac{\mathrm{d}}{\mathrm{d} \theta}\left(\sin\theta\,\frac{\mathrm{d} T_{\vec{n}}(\theta)}{\mathrm{d} \theta}\right) = -\lambda_{\vec{n}}\,T_{\vec{n}}(\theta)\,. \end{equation} If we write this equation in terms of \(x = \cos \theta\) and set \(\lambda_{\vec{n}} = -l(l+1)\), we obtain the Legendre differential equation (Equation B.86). As discussed in Appendix B.3.6, the solutions are given by the associated Legendre polynomials of the first kind \(P_l^m(x)\) and those of the second kind \(Q_l^m(x)\). For a physically-reasonable solution, we require that the solution remains finite at the poles \(\theta = 0,\pi\) or \(x=\pm1\). All Legendre functions of the second kind and all Legendre functions of the first kind with non-integer \(l\) diverge at at least one of \(x=-1\) or \(x=+1\). Therefore, physically-reasonable solutions are given by \(P_l^m(x)\), where \(l\) is an integer. These solutions further only exist if \(l \geq |m|\), which further restricts the solutions \(F_{\vec{n}}(\phi)\) to those satisfying this constraint. Finally, the Legendre differential equation is invariant under the transformation \(l \rightarrow -l-1\); therefore we only need to consider non-negative \(l\) to build a basis set (\(P_{-l}^m[x] = P_{l-1}^m[x]\)). Thus, we have that \(T_{\vec{n}}(\theta) = P_l^m(\cos \theta)\) with \(l\geq 0\) and \(|m| \leq l\).
This set of solutions \(T_{\vec{n}}(\theta)\,F_{\vec{n}}(\phi)\) often gets written in terms of a single two-dimensional function \(Y_l^m(\theta,\phi)\) known as spherical harmonics defined as \begin{equation}\label{eq-triaxialgrav-spherharm-def} Y_l^m(\theta,\phi) \equiv \sqrt{\frac{2l+1}{4\pi}\,\frac{(l-m)!}{(l+m)!}}\,P_l^m(\cos \theta)\,e^{im\phi}\,. \end{equation} The normalization of these is chosen such that these functions are orthonormal, that is \begin{equation}\label{eq-triaxialgrav-spherharm-ortho} \int_0^\pi\mathrm{d}\theta\,\sin\theta\,\int_0^{2\pi}\mathrm{d}\phi\,Y_l^{m *}(\theta,\phi)\,Y_{l'}^{m'}(\theta,\phi) = \delta_{mm'}\,\delta_{ll'}\,, \end{equation} and as we have shown, they are eigenfunctions of the two-dimensional, spherical Laplacian that satisfy \begin{equation} \nabla^2 Y_l^m(\theta,\phi) = -l(l+1)\,Y_l^m(\theta,\phi)\,. \end{equation} Because of the orthonormality of spherical harmonics, we can then easily decompose a shell’s non-uniform surface density \(\sigma(\phi,\theta)\) as \begin{equation}\label{eq-triaxialgrav-shell-decompose-1} \sigma(\phi,\theta) = \sum_{l=0}^{\infty}\sum_{m=-l}^l\,\sigma_{lm}\,Y_l^m(\theta,\phi)\,, \end{equation} where the coefficients \(\sigma_{lm}\) are given by \begin{equation}\label{eq-triaxialgrav-shell-decompose-2} \sigma_{lm} = \int_0^\pi\mathrm{d}\theta\,\sin\theta\,\int_0^{2\pi}\mathrm{d}\phi\,Y_l^{m *}(\theta,\phi)\,\sigma(\phi,\theta)\,. \end{equation} Because the differential equations (12.60) and (12.62) are of the Sturm-Liouville type, their set of solutions forms a complete basis for all smooth functions on the sphere and we can therefore decompose any smooth surface density in terms of spherical harmonics.
Now that we have decomposed the surface density of the non-uniform shell, we can return to solving the three-dimensional Poisson equation and we solve it for each \((l,m)\) component of the shell separately before summing over all \((l,m)\) contributions. That is, we solve \begin{equation} \nabla^2 \Phi = 4\pi G\,\sigma_{lm}\,\delta(r-a)\,Y_l^m(\theta,\phi)\,. \end{equation} Because the \(Y_l^m(\theta,\phi)\) are eigenfunctions of the Laplacian that satisfy \(\nabla^2 Y_l^m(\theta,\phi) = -l(l+1) Y_l^m(\theta,\phi)\), we can find the solution of this equation as \(\Phi = R(r)\,Y_l^m(\theta,\phi)\), where \(R(r)\) satisfies (using the Laplacian in spherical coordinates from Equation A.8) \begin{equation} \frac{\mathrm{d}}{\mathrm{d} r}\left(r^2\,\frac{\mathrm{d} R(r)}{\mathrm{d} r}\right) -l(l+1)\,R(r) = 4\pi G\,\sigma_{lm}\,r^2\,\delta(r-a)\,. \end{equation} We can solve this equation by solving it at \(r < a\) and \(r > a\) and applying the appropriate boundary condition at \(r=a\). At \(r \neq a\), we have that \begin{equation} \frac{\mathrm{d}}{\mathrm{d} r}\left(r^2\,\frac{\mathrm{d} R(r)}{\mathrm{d} r}\right) -l(l+1)\,R(r) = 0\,, \end{equation} which has the general solution \begin{equation} R(r) = A\,r^l + B\,r^{-(l+1)}\,. \end{equation} Thus, inside the shell, we have that \(R(r) = A_{\mathrm{in}}\,r^l + B_{\mathrm{in}}\,r^{-(l+1)}\), while outside the shell \(R(r) = A_{\mathrm{out}}\,r^l + B_{\mathrm{out}}\,r^{-(l+1)}\). We set \(A_{\mathrm{out}} = 0\) such that \(R(r) \rightarrow 0\) (and thus \(\Phi(r,\phi,\theta) \rightarrow 0\)) when \(r \rightarrow \infty\). The radial component of the gravitational field is \(-Y_l^m(\theta,\phi)\,\mathrm{d} R / \mathrm{d} r\) and Gauss’s law (Chapter 2.2) then states that \begin{equation}\label{eq-triaxialgrav-multipole-gauss} 4\pi GM(< r) = \int_0^\pi\mathrm{d}\theta\,\sin\theta\,\int_0^{2\pi}\mathrm{d}\phi\,Y_l^{m}(\theta,\phi)\,\left[Al\,r^{l+1} - B(l+1)\,r^{-l}\right]\,, \end{equation} where \(M(<r)\) is the enclosed mass within a radius \(r\). When \(r < a\), \(M(<r) = 0\) and the right-hand side of this equation should therefore vanish. The first term in the integrand always does, because \(l\,\int_0^\pi\mathrm{d}\theta\,\sin\theta\,\int_0^{2\pi}\mathrm{d}\phi\,Y_l^{m}=0\) for all \((l,m)\). The second term doesn’t, because it diverges as \(r^{-l}\) as \(r \rightarrow 0\). Therefore, we must set \(B_{\mathrm{in}}=0\) to satisfy Gauss’s law. Because no work is done when passing through an infinitesimally-thin shell, the potential at \(r=a\) needs to be continuous and therefore \(A_{\mathrm{in}}\,a^l = B_{\mathrm{out}}\,a^{-(l+1)}\) or \(B_{\mathrm{out}} = A_{\mathrm{in}}\,a^{2l+1}\). Finally, to determine \(A_{\mathrm{in}}\) we integrate the Poisson equation over an infinitesimal volume that encompasses a part of the spherical shell; similar to the derivation of Equation (7.9) that relates the vertical field and the surface-density for a razor-thin disk, we have in this case that \begin{equation} \left({\partial \Phi_{\mathrm{out}} \over \partial r}\right)-\left({\partial \Phi_{\mathrm{in}} \over \partial r}\right) = 4\pi G\,\sigma_{lm}\,Y_l^m(\theta,\phi)\,. \end{equation} Multiplying this equation by \(Y_l^{m *}(\theta,\phi)\) and integrating over \((\phi,\theta)\) using \(\Phi_{\mathrm{in/out}} = R(r)_{\mathrm{in/out}}\,Y_l^m(\theta,\phi)\) and the orthonormality condition from Equation (12.64), we find that \begin{equation} -B_{\mathrm{out}}(l+1)\,a^{-l-2}-A_{\mathrm{in}}l\,a^{l-1} = 4\pi G\,\sigma_{lm}\,, \end{equation} and using that \(B_{\mathrm{out}} = A_{\mathrm{in}}\,a^{2l+1}\), we find that \(A_{\mathrm{in}} = -4\pi G \sigma_{lm} a^{1-l}/(2l+1)\). Thus, the potential of a spherical shell with radius \(a\) and with non-uniform surface density \(\sigma_{lm}\,Y_l^{m}(\theta,\phi)\) is given by \begin{equation} \Phi(r,\phi,\theta) = -{4\pi G\,a\,\sigma_{lm}\over 2l+1}\,Y_l^{m}(\theta,\phi)\,\begin{cases} \left(r/a\right)^l& \quad r \leq a\\\left(a/r\right)^{l+1}& \quad r > a\end{cases}\,. \end{equation} For a uniform shell, we recover Newton’s shell theorems from Chapter 2.2: \(\sigma_{00}\,Y_0^{0}(\theta,\phi) = M/(4\pi\,a^2)\) in this case and therefore (a) the potential within the shell is constant and equal to \(-4\pi G\,a\,M/(4\pi\,a^2) = -GM/a\) and (b) outside of the shell, the potential is \(-4\pi G\,a\,M/(4\pi\,a^2)\,(a/r) = -GM/r\).
Using the decomposition from Equation (12.66), we then immediately arrive at the gravitational potential of a spherical shell with an arbitrary surface density \(\sigma(\phi,\theta)\) \begin{equation}\label{eq-triaxialgrav-shell-potential} \Phi(r,\phi,\theta) = -4\pi G\,a\,\begin{cases} \sum_{l=0}^{\infty}{ \left(r/a\right)^{l\phantom{+1}}\sum_{m=-l}^{l}{ {\sigma_{lm}\over 2l+1}\,Y_l^{m}(\theta,\phi)}}& \quad r \leq a\\\sum_{l=0}^{\infty}{ \left(a/r\right)^{l+1}\sum_{m=-l}^{l}{ \,{\sigma_{lm}\over 2l+1}\,Y_l^{m}(\theta,\phi)}}& \quad r > a\end{cases}\,. \end{equation} At large distances from the shell, the monopole (\(l=0\)) always dominates and higher-order multipoles decay faster than lower-order ones. This makes physical sense because far away from a shell, its detailed structure should become unimportant and the potential should approach that of a point mass.
For an extended body, we can decompose the density \(\rho(r,\phi,\theta)\) into a set of shells as (using \(a\) to denote the \(r\) coordinate) \begin{equation}\label{eq-triaxialgrav-body-decompose-1} \rho(a,\phi,\theta) = \sum_{l=0}^{\infty}\sum_{m=-l}^l\,\rho_{lm}(a)\,Y_l^m(\theta,\phi)\,, \end{equation} where the coefficients \(\rho_{lm}(a)\) are given by \begin{equation}\label{eq-triaxialgrav-body-decompose-2} \rho_{lm}(a) = \int_0^\pi\mathrm{d}\theta\,\sin\theta\,\int_0^{2\pi}\mathrm{d}\phi\,Y_l^{m *}(\theta,\phi)\,\rho(a,\phi,\theta)\,. \end{equation} This corresponds to an infinite series of shells with surface density expressed in spherical harmonics \(\sigma_{lm} = \rho_{lm}\mathrm{d}a\) and we can thus obtain the full potential by integrating over all shells. This gives \begin{equation}\label{eq-triaxialgrav-multipole-potential} \Phi(r,\phi,\theta) = -4\pi G\,\sum_{l=0}^{\infty}{\sum_{m=-l}^{l}{{Y_l^{m}(\theta,\phi)\over 2l+1}\,\left[{1\over r^{l+1}}\,\int_0^r\mathrm{d}a\,a^{l+2}\,\rho_{lm}(a)+r^l\,\int_r^\infty\mathrm{d}a\,a^{1-l}\,\rho_{lm}(a)\right]}}\,. \end{equation} For a spherical density we have that \(\rho_{00}(a)\,Y^0_0(\theta,\phi) = \rho(a)\) and therefore the expression simplifies to that for the potential of a spherical density (Equation 2.21). Equation (12.79) is the multipole expansion of the gravitational potential, with the \((l,m) = (0,0)\) term the monopole, the \(l=1\) terms are dipole terms, the \(l=2\) terms are quadrupole terms, \(l=3\) terms are octopole terms, \(l=4\) terms are hexadecapole terms, and so on. If we put everything together, computing the gravitational potential through this multipole expansion is quite complicated, because it involves a triple integral and a double sum. However, the situation is not as bad as it sounds because (i) we generally do not need to go to very high multipole orders (large \(l\) and \(|m|\)) to get a close match to the density, (ii) for tasks like orbit integration, the multipole coefficient functions \(\rho_{lm}(a)\) only need to be computed once, and (iii) because spherical harmonic decompositions are so common in many areas of science, fast numerical codes exist to efficiently perform the integrals over \((\phi,\theta)\) and sums over \((l,m)\).
We have worked out the multipole expansion here in terms of spherical harmonics and an expansion in terms of non-uniform, spherical shells. However, we could have similarly expanded the density in terms of a different set of orthogonal functions, for example, as a set of non-uniform spheroidal or ellipsoidal shells (using, e.g., oblate or prolate spheroidal harmonics). For galaxies or dark-matter halos that are close to having a constant shape or orientation, this could create a more compact representation of the gravitational potential as fewer multipole terms would be needed to attain a given accuracy in the density representation. But in practice these methods are not used in galactic astrophysics, largely because the generalizations of the spherical harmonics that are required are not a standard part of an astrophysicist’s toolbox. These methods are used in situations where very high-precision in the gravitational field is required, for example, for the Earth’s gravitational field, and planets and asteroids in the solar system that one wants to, for example, land spacecraft on (e.g., Hofmann-Wellenhof & Moritz 2005; Maus 2010; Fukushima 2014).
12.3.2. Basis-function expansions¶
In the multipole-expansion method from the previous section, we expand densities and potentials using a countable set of basis-functions in two dimensions combined with a continuous set in the third dimension. However, for smooth density distributions, we can replace the continuous set with a countable set of complete basis functions in the third dimension as well and obtain the potential as a three-dimensional sum. This is what we will refer to as a basis-function expansion, even though the multipole expansion described above also uses basis functions.
Following the approach we have followed in the previous section to derive spherical harmonics as a complete, orthonormal basis set of smooth functions on the sphere, the obvious way to create a three-dimensional set of basis functions is to look for eigenfunctions of the three-dimensional Laplace operator using the separation of variables approach. Similar to what we did above, we look for functions functions \(\Phi(r,\phi,\theta) = R_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\) that satisfy \begin{equation}\label{eq-triaxialgrav-3dlaplace-eigenfunctions} \nabla^2 \left[R_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\right] = \lambda_{\vec{n}}\,R_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\,, \end{equation} where \(\vec{n}\) is again a vector indexing the eigenfunctions that we will determine as part of the solution and where \(\lambda_{\vec{n}}\) are to-be-determined eigenvalues. Because it will be useful later on, to work out this equation we generalize it to \(\Phi(r,\phi,\theta) = U_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\) and work out \(\nabla^2 \left[U_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\right] = \lambda_{\vec{n}}\,R_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\); eigenfunctions are then obtained by additionally requiring that \(U_{\vec{n}}(r) \equiv R_{\vec{n}}(r)\). Using the Laplacian in spherical coordinates from Equation (A.8), we can work out Equation (12.80) as \begin{align} \frac{F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)}{r^2}\,&\frac{\mathrm{d}} {\mathrm{d} r}\left(r^2\,\frac{\mathrm{d} U_{\vec{n}}(r)}{\mathrm{d} r}\right) + \frac{U_{\vec{n}}(r)\,F_{\vec{n}}(\phi)}{r^2\sin\theta}\,\frac{\mathrm{d}}{\mathrm{d} \theta}\left( \sin\theta\,\frac{\mathrm{d} T_{\vec{n}}(\theta)}{\mathrm{d} \theta}\right) \nonumber\\ &\quad +\frac{U_{\vec{n}}(r)\,T_{\vec{n}}(\theta)}{r^2\sin^2\theta}\,\frac{\mathrm{d}^2 F_{\vec{n}}(\phi)}{\mathrm{d} \phi^2} = \lambda_{\vec{n}}\,R_{\vec{n}}(r)\,F_{\vec{n}}(\phi)\,T_{\vec{n}}(\theta)\,. \end{align} This equation is similar to Equation (12.58) above and we can follow the usual approach of bringing all dependence on a single coordinate dimension on one side and the remaining coordinate dimensions to the other side. Similar to how we handled Equation (12.58), we can bring all \(\phi\) dependence to the right-hand side and all \((r,\theta)\) dependence to the left hand-side and both sides have to then equal a constant that we again call \(m^2\). This again leads to Equation (12.60) for \(F(\phi)\) with the solution of Equation (12.61). Following the same approach of separating coordinate dimensions for the \((r,\theta)\) part of the separated equation leads to \begin{align}\label{eq-poisson-separate-spherical-secondstep} \frac{1}{U_{\vec{n}}(r)}\,\frac{\mathrm{d}} {\mathrm{d} r}\left(r^2\,\frac{\mathrm{d} U_{\vec{n}}(r)}{\mathrm{d} r}\right) - \, \lambda_{\vec{n}}\,r^2\,\frac{R_{\vec{n}}(r)}{U_{\vec{n}}(r)} = \frac{m^2}{\sin^2\theta}-\frac{1}{\sin\theta}\,\frac{1}{T_{\vec{n}}(\theta)}\,\frac{\mathrm{d}}{\mathrm{d} \theta}\left(\sin\theta\,\frac{\mathrm{d} T_{\vec{n}}(\theta)}{\mathrm{d} \theta}\right) \,. \end{align} Again, because the left-hand side only depends on \(r\) and the right-hand side only depends on \(\theta\), both sides need to be equal to a constant, say \(l(l+1)\). The right-hand side then reduces to the same differential equation as Equation (12.62) from the previous section if we set \(\lambda_{\vec{n}} = -l(l+1)\) in that equation and the solution is therefore again given in terms of associated Legendre polynomials. The combination \(T_{\vec{n}}(\theta)\,F_{\vec{n}}(\phi)\) is therefore again equal to the spherical harmonics \(Y_l^m(\theta,\phi)\) and two dimensions of the eigenvalue vector \(\lambda_{\vec{n}}\) are therefore given by the integers \((l,m)\) that satisfy \(l \leq 0\) and \(|m| \leq l\). Finally, we turn to the left-hand side of Equation (12.82), which gives an ordinary differential equation for the radial part of the basis \begin{equation}\label{eq-poisson-separate-spherical-radial-basis} r^2\,\frac{\mathrm{d}^2 U_{\vec{n}}(r)}{\mathrm{d} r^2} + 2\,r\,\frac{\mathrm{d} U_{\vec{n}}(r)}{\mathrm{d} r}+\left(-\lambda_{\vec{n}}\,{R_{\vec{n}}(r) \over U_{\vec{n}}(r)}\,r^2-l\,[l+1]\right)\,U_{\vec{n}}(r) = 0\,. \end{equation} If we require \(U_{\vec{n}}(r) \equiv R_{\vec{n}}(r)\) such that we find eigenfunctions, this equation is a form of the defining equation (B.57) of spherical Bessel functions (see Appendix B.3.5) and the basis functions can therefore be written in terms of spherical Bessel functions (e.g., Weinberg 1999). While spherical Bessel functions provide a complete basis (their defining differential equation is of the Sturm-Liouville type), these are oscillating functions that on their own do not well represent a galactic mass distribution.
If we were to use Bessel functions, we would require a large number of expansion coefficients to represent the density. Without requiring that \(U_{\vec{n}}(r)\) is an eigenfunction, however, we have much greater flexibility in finding a complete set of functions that solves Equation (12.83). A general approach is to find solutions where the lowest-order basis-function is a realistic galactic potential–density pair \([\Phi_0(r),\rho_0(r)]\), positing a general solution of the form \begin{align} \Phi_{nlm}(r,\phi,\theta) & = \phantom{\lambda_{nl}\,}\sqrt{4\pi}\,\Phi_0(r)\,U_{nl}(r)\,Y_l^m(\theta,\phi)\,,\\ \rho_{nlm}(r,\phi,\theta) & = \lambda_{nl}\,\sqrt{4\pi}\,\rho_0(r)\,U_{nl}(r)\,Y_l^m(\theta,\phi)\, \end{align} (e.g., Hernquist & Ostriker 1992; Weinberg 1999), where the factor of \(\sqrt{4\pi}\) is introduced to cancel \(Y_0^0 = 1/\sqrt{4\pi}\) in order to make the lowest-order function the \([\Phi_0(r),\rho_0(r)]\) pair. Plugging the radial part of this form into Equation (12.83), we get \begin{equation}\label{eq-poisson-separate-spherical-radial} r^2\,\frac{\mathrm{d}^2 [\Phi_0(r)\,U_{nl}(r)]}{\mathrm{d} r^2} +2r\,\frac{\mathrm{d} [\Phi_0(r)\,U_{nl}(r)]}{\mathrm{d} r} -l(l+1)\,\Phi_0(r)\,U_{nl}(r) = r^2\,\lambda_{nl}\,\rho_0(r)\,U_{nl}(r)\,. \end{equation} For any choice of \([\Phi_0(r),\rho_0(r)]\), this equation can be written in Sturm-Liouville form and the solutions \(U_{nl}(r)\) therefore always provide a complete, orthogonal basis for the space of functions under consideration. Because of this orthogonality, decomposing any density \(\rho(r,\phi,\theta)\) into \(\rho_{nlm}(r,\phi,\theta)\) components is straightforward \begin{equation} \rho(r,\phi,\theta) = \sum_{n=0}^\infty{\sum_{l=0}^{\infty}{\sum_{m=-l}^{l}{A_{nlm}\,\rho_{nlm}(r,\phi,\theta)}}}\,, \end{equation} where the expansion coefficients are given by \begin{equation}\label{eq-triaxialgrav-bfe-coeffs} A_{nlm} = {\int \mathrm{d}r\,\mathrm{d}\phi\,\mathrm{d}\theta\,r^2\,\sin\theta\,\rho(r,\phi,\theta)\,\Phi_{nlm}^*(r,\phi,\theta) \over \int \mathrm{d}r\,\mathrm{d}\phi\,\mathrm{d}\theta\,r^2\,\sin\theta\,\rho_{nlm}(r,\phi,\theta)\,\Phi_{nlm}^*(r,\phi,\theta)}\,, \end{equation} where \(\Phi_{nlm}^*\) is the complex conjugate of \(\Phi_{nlm}\) (in practice, by combining terms with opposite-sign \(m\), we can compute two sets of real coefficients by replacing \(e^{im\phi}\) with either of \(\cos[m\phi]\) and \(\sin[m\phi]\)). The gravitational potential is then given by \begin{equation} \Phi(r,\phi,\theta) = \sum_{n=0}^\infty{\sum_{l=0}^{\infty}{\sum_{m=-l}^{l}{A_{nlm}\,\Phi_{nlm}(r,\phi,\theta)}}}\,, \end{equation} A general approach to coming up with complete sets of basis functions that can be used to represent an arbitrary density is then to choose the \([\Phi_0(r),\rho_0(r)]\) pair to closely resemble to overall radial profile of the density one wants to model and to then solve Equation (12.86). This approach combines the advantages of using eigenfunctions with the flexibility to choose physically-meaningful radial basis functions and if the lowest-order model is a close match to the full density, few terms are necessary in the expansion to accurately represent the density. Equation (12.86), however, in general cannot be solved in terms of known functions and for most \([\Phi_0(r),\rho_0(r)]\) pairs one therefore has to solve it numerically (e.g., for the NFW profile; Dai et al. 2018). Because Sturm-Liouville problems are so common, software exists to compute such numerical solutions efficiently (see, e.g., Weinberg 1999 for an application of this). From the definition, we know that \(U_{00}(r) = 1\) and \(\lambda_{00} =1\).
For computational convenience and speed, it is useful to have analytical basis-function sets and Equation (12.86) can in fact be solved for two important cases: (a) if \([\Phi_0(r),\rho_0(r)]\) is the Plummer potential–density pair (see Chapter 2.4.3; Clutton-Brock 1973) or (b) when \([\Phi_0(r),\rho_0(r)]\) is the Hernquist potential–density pair (see Chapter 2.4.6; Hernquist & Ostriker 1992). These provide good sets of basis functions to represent the potential of globular clusters (for the Plummer-based set) and dark-matter halos and elliptical galaxies (for the Hernquist-based set). Because it is widely used in galactic dynamics and to further illustrate the approach, we briefly discuss the Hernquist & Ostriker (1992) method based on the Hernquist potential–density pair, which the authors refer to as the self-consistent field method.
For convenience, we work with a unit mass Hernquist profile with scale radius \(a=1\); because a non-unit mass is a simple rescaling of the density and potential and because a non-zero scale radius simply scales all positions as \(r/a\), it is straightforward to apply the formulae below in these cases. The Hernquist potential–density pair is then derived from Equations (2.49), (2.51), and (2.53) \begin{align} \Phi_0(r) & = -{1\over 1+r}\,;\quad \rho_0(r) = {1\over 2\pi}\,{1\over r\,(1+r)^3}\,. \end{align} The main task at hand is then to solve Equation (12.86). To make headway in this, Hernquist & Ostriker (1992) posited that the \(U_{nl}(r)\) have the following form \begin{equation}\label{eq-triaxialgrav-scf-Uansatz} U_{nl}(r) = {r^l\over (1+r)^{2l}}\,W_{nl}(\xi)\,. \end{equation} This form is inspired by the multipole-expansion that we discussed above (e.g., Equation 12.76), where the potential is \(\propto r^l\) at small \(r\) and \(\propto 1/r^{l+1}\) at large \(r\). If \(W_{0l}(\xi)=1\), then the lowest-order radial basis functions have this behavior at \(r \ll 1\) and \(r \gg 1\), respectively. The unknown functions \(W_{nl}(\xi)\) then represent deviations from this basic behavior through the higher-order radial basis functions. We write \(W_{nl}\) as a function of a coordinate \(\xi(r)\) that is determined by the requirement that the differential equation for \(W_{nl}\) that results from plugging Equation (12.91) into Equation (12.86) is of known form. Hernquist & Ostriker (1992) show that it is useful to define \(\xi(r)\) as \begin{equation} \xi = {r-1\over r+1}\,, \end{equation} because then (12.86) becomes \begin{equation}\label{eq-triaxialgrav-scf-define} (1-\xi^2)\,{\mathrm{d}^2 W_{nl}\over \mathrm{d}\xi^2}-4\,(l+1)\,\xi\,{\mathrm{d}W_{nl}\over \mathrm{d}\xi}+2\,[K_{nl}-(l+1)(2l+1)]\,W_{nl} =0\,, \end{equation} where we have defined \(K_{nl} = \lambda_{nl}/[2\pi]\). This may not seem simple, but it is equivalent to the equation defining the Gegenbauer (or ultraspherical) polynomials \(C_n^{(\alpha)}(\xi)\) if \(\alpha = 2l+3/2\) and \(K_{nl}\) is given by (see Equation B.81) \begin{equation} K_{nl} = {1\over 2}\,n\,(n+4l+3)+(l+1)(2l+1)\,. \end{equation} Thus, the full orthogonal basis-function set is \begin{align} \Phi_{nlm}(r,\phi,\theta) & = \!\!\!\!\!\phantom{{K_{nl}\over 2\pi}\,}-\sqrt{4\pi}\,{r^l\over (1+r)^{2l+1}}\,C_n^{(2l+3/2)}(\xi)\,Y_l^m(\theta,\phi)\,,\\ \rho_{nlm}(r,\phi,\theta) & = \phantom{-}\!\!{K_{nl}\over 2\pi}\,\sqrt{4\pi}\,{r^{l-1}\over (1+r)^{2l+3}}\,C_n^{(2l+3/2)}(\xi)\,Y_l^m(\theta,\phi)\,. \end{align} The expansion coefficients can be computed using Equation (12.88).
The approach we have followed so far in this section works well for densities that are not too far from spherical (like elliptical galaxies and dark matter halos). But it does not work well for disks, which are so far from spherical that they would require very high order in \(l\) to be well-described by an expansion in terms of spherical harmonics. Thus, to represent disk potentials through basis-function expansions, it is necessary to come up with different sets of basis functions. Indeed, many of the first applications of basis-function expansions in galactic dynamics were for disks (e.g., Clutton-Brock 1972b; Kalnajs 1976). These first applications were for the gravitational field within a razor-thin disk, where the separation-of-variables approach to the Helmholtz equation in two dimensions leads to basis functions that are combinations of Bessel functions for the radial component and \(e^{im\phi}\) for the azimuthal component (this can be shown by following a similar approach as we took in Chapter 7.3.5). Bessel functions are again not realistic galactic components, so we can follow a similar approach as we did in the three-dimensional case above of finding radial potential-density pairs for which the lowest-order component is a realistic disk density and that form a complete basis.
The situation is trickier for thick disks (that is, not razor-thin in this context). For those we could again look for a complete set of basis functions by solving the Helmholtz equation using the separation-of-variables approach in cylindrical coordinates. This leads to basis functions that are combinations of (i) Bessel functions for the radial component, (ii) \(e^{im\phi}\) for the azimuthal component, and (iii) \(e^{ihz}\) for the vertical component. However, the vertical profile of disks is not well described as \(e^{ihz}\) and, thus, we now have the unfortunate situation where two of the dimensions of our basis functions do not describe realistic galactic density profiles. The approach that we followed above of replacing the offending radial eigenfunctions with a complete set of potential-density pairs is now much harder to follow, because we need to replace both the radial and vertical eigenfunctions. Earn (1996) describes an approach of this form, but it leads to a complete set of basis functions that is not orthogonal. In that case, we cannot compute expansion coefficients using an expression like Equation (12.88), but instead we need to solve for the coefficients \(A_{\vec{n}}\) by solving the equation \begin{equation} \int\mathrm{d}\vec{x}\,\Phi^*_\vec{n}(\vec{x})\,\rho(\vec{x}) = \sum_{\vec{m}}{M_{\vec{n}\vec{m}}\,A_{\vec{n}}}\,, \end{equation} where the matrix \(M_{\vec{n}\vec{m}}\) is given by \begin{equation} M_{\vec{n}\vec{m}} = \int\mathrm{d}\vec{x}\,\Phi^*_{\vec{n}}(\vec{x})\,\rho_{\vec{m}}(\vec{x})\,. \end{equation} Alternatively, we can orthogonalize a non-orthogonal set using Gram-Schmidt orthogonalization, but in the end the computational cost is similar. When using a basis-function expansion method, solving for the expansion coefficients is generally not the most time-consuming step, so the fact that we now have to solve a matrix equation rather than directly computing the coefficients is not too burdensome. Aside from this complication, the basis-function expansion in terms of disky basis functions is not that different from that which we described in detail above and we refer the reader to Earn (1996) for more details. Robijn & Earn (1996) and Brown & Papaloizou (1998) give similar, related approaches.
Using a simple trick, the basis-function approach based on spherical harmonics discussed in this section can also be applied to derive the potential of general disk mass distributions, by using the basis-function expansion to describe the difference between the actual density and a good disky approximation to it. If the approximation is good enough, the difference is much closer to spherical than the original density and the difference can therefore efficiently be represented using either the multipole expansion or a basis-function expansion based on spherical harmonics. See Kuijken & Dubinski (1995) and Dehnen (1998) for more details.
.
.
.