Multipole Moments Make Life Easier

The Multipole Expansion in Python

Physics

Photo by Hans-Jurgen Mager on Unsplash

In electrodynamics as well as in gravitational physics, we are often faced with the problem to calculate some potential (electric or gravitational) based on some distribution of charge or matter. The solution to this problem is easily stated, but hard to compute. But there is an easier approach. If we are only interested in how the potential looks from far away, we can generate an approximate solution much easier. In fact, we can make this approximation as accurate as we want. We are talking about the multipole expansion.

Suppose we have some electric charge distribution $\rho$, which is a function of the space variables $x, y, z$. We know that a charge distribution causes a potential which can in principle be calculated from

$$ \Phi({\bf r}) = \frac{1}{4\pi}\int \frac{\rho({\bf r}')}{|{\bf r}-{\bf r'}|}d^3r' $$

where the integration runs over all space. In practice, this integral is hard to compute, even numerically. The problem is, that we would have to compute an integral for each and every point in space! But we can resort to a simple idea, which makes life much easier.

Marcus Guy, CC BY-SA 3.0 https://creativecommons.org/licenses/by-sa/3.0, via Wikimedia Commons

Consider the same charge distribution, but imagine we are very very far away from these charges. Then we don't really see the detailed structure, and effectively, we only see kind of a point charge. This is the monopole. If we want to be a bit more accurate, we can imagine, that there is some additional dipole (two point charges) that give some improvement to our approximation. And if that is still not accurate enough, we could add a configuration of quadrupoles (four point charges) and another improvent. So basically we have

$$ \Phi({\bf r}) \approx \Phi_{mono}({\bf r}) + \Phi_{dipole}({\bf r}) + \Phi_{quadrupole}({\bf r}) + \ldots $$

The thing is, each of these terms can be computed much more easily than the original integral. And there is a rigorous way how to generate as many correction terms as we want. Let's see how.

The crucial point to know is that you can expand the function $1/|({\bf r})-({\bf r}')|$ in terms of spherical harmonics. This is called Laplace expansion.

$$ \frac{1}{|{\bf r}-{\bf r}'|} = \sum_{l=0}^\infty \frac{4\pi}{2l+1}\sum_{m=-l}^l \frac{r_<^{l}}{r_>^{l+1}} Y^*_{lm}(\theta', \varphi') Y_{lm}(\theta, \varphi) $$

where $\theta, \varphi$ and $\theta', \varphi'$ are the spherical coordinates of the points ${\bf r}$ and ${\bf r}'$, and $Y_{lm}$ are the spherical harmonics.

If you need to brush up your knowledge about spherical harmonics, you may want to have a look at this article .

The fraction with the $r$'s means

$$ r_< = \min(r, r')\\r_> = \max(r, r') $$

I must admit, that this expression looks much more complicated than the original problem. But it isn't!

In our case, where we are looking from far away $r >> r'$, we have

$$ r_< = r'\\ r_> = r $$

So if we insert the Laplace expansion in our integral, we get

$$ \Phi({\bf r}) = \frac{1}{4\pi}\sum_{l=0}^\infty \sum_{m=-l}^l \frac{4\pi}{(2l+1)r^{l+1}}Y_{lm}(\theta, \varphi) \int r'^l \rho({\bf r'} )Y^*_{lm} d^3r' $$

Now we simply define

$$ Q_{lm} = \sqrt\frac{4\pi}{2l+1}\int r'^l \rho({\bf r'} )Y^*_{lm} d^3r' $$

which are called multipole moments, and we have our multipole expansion

$$ \Phi({\bf r}) = \frac{1}{4\pi} \sum_{l=0}^\infty \sum_{m=-l}^l \sqrt\frac{4\pi}{2l+1}\frac{Q_{lm}Y_{lm}(\theta, \varphi)}{(2l+1)r^{l+1}} $$

So we have effectively replaced our original problem of calculating an integral for every point in space by calculating a different integral for every $(l, m)$ and then summing up. For a give level of accuracy, we just have to compute a small number of these integrals, and this is the big advantage of the multipole method.

Now let's see how we can use multipole expansions in Python. There is a convenient package, aptly called multipoles, that greatly simplify our job. Simply install it using pip install multipoles. As an example, we will calculate the multipole for a charge distribution

$$ \rho({\bf r}) = e^{-|{\bf r}-{\bf r_A}|/\sigma_A^2} - e^{-|{\bf r}-{\bf r_B}|/\sigma_B^2} $$


! pip install multipoles


from multipoles import MultipoleExpansion
import numpy as np

# First we set up our grid, a cube of length 10 centered at the origin:

npoints = 101
edge = 10
x, y, z = [np.linspace(-edge/2., edge/2., npoints)]*3
XYZ = np.meshgrid(x, y, z, indexing='ij')


# We model our smeared out charges as gaussian functions:

def gaussian(XYZ, xyz0, sigma):
    g = np.ones_like(XYZ[0])
    for k in range(3):
        g *= np.exp(-(XYZ[k] - xyz0[k])**2 / sigma**2)
    g *= (sigma**2*np.pi)**-1.5
    return g

sigma = 1.5   # the width of our gaussians

# Initialize the charge density rho, which is a 3D numpy array:
rho = gaussian(XYZ, (0, 0, 1), sigma) - gaussian(XYZ, (0, 0, -1), sigma)


# Prepare the charge distribution dict for the MultipoleExpansion object:

charge_dist = {
    'discrete': False,     # we have a continuous charge distribution here
    'rho': rho,
    'xyz': XYZ
}

# The rest is the same as for the discrete case:

l_max = 2   # where to stop the infinite multipole sum; here we expand up to the quadrupole (l=2)

Phi = MultipoleExpansion(charge_dist, l_max)

The multipole moments can then be obtained by the property multipole_moments


Phi.multipole_moments


{(0, 0): (4.169275454499798e-17+0j),
 (1, 0): (1.9992843202867914+0j),
 (1, 1): (-4.1021078134943617e-19+1.5215054472701573e-17j),
 (1, -1): (4.1021078134943617e-19+1.5215054472701573e-17j),
 (2, 0): (1.1877987280187917e-16+0j),
 (2, 1): (5.936664102910997e-17+2.0706668214237114j),
 (2, -1): (-5.936664102910997e-17+2.0706668214237114j),
 (2, 2): (4.43852837391959e-19-8.631695303251424e-19j),
 (2, -2): (4.43852837391959e-19+8.631695303251424e-19j)}

We don't have to sum up all the terms manually, the package does the job for us. Just evaluate the object at the point of interest. For example to evaluate the expansion at ${\bf r } = (30, 40, 50)$, we can say


Phi(30, 40, 50)


0.00028847990030335093

In conclusion, the multipole expansion serves as a compelling technique that can significantly simplify the calculations of electric or gravitational potential. Rather than performing a complex integration for each point in space, the multipole expansion provides an efficient workaround, reducing the computations to a smaller number of distinct integrals for each (l, m) pair and summing up.

By using Python and the multipoles package, we can conveniently perform these calculations, modeling different charge or matter distributions, and even control the level of accuracy we desire. The example shared in this article offers a practical demonstration of how the multipole expansion technique can be applied using Python, significantly easing the process of scientific computation in physics.