Understanding the Coulomb Problem

In both gravitational physics and electrodynamics, a frequent challenge arises known as the Coulomb problem. Essentially, the question is this: given a certain distribution of matter in space, such as stars in a galaxy or galaxies within a supercluster, what is the corresponding Newtonian gravitational potential? Alternatively, if you have a specific charge distribution, what would be the electrostatic potential? Mathematically, this scenario can be expressed through a Poisson equation, characterized by special boundary conditions.

In contrast to typical scenarios with the Poisson equation, the Coulomb problem stipulates boundary conditions that extend to infinity.

This problem is indeed quite prevalent, yet surprisingly, there are limited resources online that provide guidance on numerical solutions. Most discussions focus on finite domains; however, our interest lies in infinite domains, which we will explore here.

An Ineffective Analytic Solution

While the Coulomb problem may appear straightforward at first glance, it is deceptively complex. Although one can quickly derive an analytic solution as taught in electrodynamics courses, this solution proves to be largely impractical, as we will soon demonstrate.

To derive this analytic solution, consider the charge located in an infinitesimal volume element (d^3r):

This infinitesimal charge behaves like a point charge and, as learned in high school, produces a Coulomb potential that varies inversely with distance. Therefore, for this volume element, we can express it as:

To find the total potential, we can "sum" the contributions from all other volume elements. The linearity of the Poisson equation allows us to achieve this:

This analytic solution inherently satisfies the boundary condition at infinity. However, as mentioned, it is not very useful. The underlying reason is that for most interesting scenarios, the integral cannot be solved analytically. Furthermore, numerically, we face challenges because we essentially need to compute an integral for every point (r) in space. Thus, we must seek alternative approaches.

Numerical Solutions

When tackling the standard Poisson equation, we usually discretize our finite domain by employing a grid. We will follow this method here, even though our boundary conditions suggest an infinitely large domain. Initially, we need a strategy to solve the Poisson equation for a finite domain with Dirichlet boundary conditions. We will utilize the Jacobi method due to its simplicity and stability. Although not the fastest method available, it can later be replaced by more advanced techniques like successive over-relaxation or multigrids.

When we discretize the Laplace operator, we define the grid point ((i, j, k)) as follows:

This approach employs second-order finite differences, where (h) represents the grid spacing. Note that this reflects the principle that the Laplacian is the difference between the local point and the average of its neighbors, as discussed in previous literature.

Inserting this into the Poisson equation and solving for (Phi) at ((i, j, k)) yields:

This sets the stage for implementation in Python, prior to addressing how to manage the boundary conditions at infinity.

import numpy as np

def jacobi_step(f, rho, h):

"""Performs a single Jacobi step on the array f and the charge distribution rho."""

f_new = f.copy()

f_new[1:-1, 1:-1, 1:-1] = (

f[2:, 1:-1, 1:-1] + f[:-2, 1:-1, 1:-1] +

f[1:-1, 2:, 1:-1] + f[1:-1, :-2, 1:-1] +

f[1:-1, 1:-1, 2:] + f[1:-1, 1:-1, :-2]

) + 4 * np.pi * rho[1:-1, 1:-1, 1:-1] * h ** 2

f_new[1:-1, 1:-1, 1:-1] /= 6

return f_new

Addressing Boundary Conditions

Now, what about the boundary conditions at infinity? Since we cannot have an infinite domain on a computer, we need a method to approximate the potential at the boundary. If our grid is sufficiently large, we can apply the concept of multipoles.

From a distance, if you are sufficiently far from the charge distribution, it resembles a point charge. As you approach, you may perceive a dipole, and with closer examination, you can identify more contributions from complex arrangements of virtual point charges. This concept of multipoles will aid us in calculating the potential at the boundary.

For clarity, we will utilize the method of manufactured solutions to generate a scenario that allows for straightforward comparisons, as we have the exact solution by definition.

This manufactured solution represents a smoothed Coulomb potential, where the smoothing parameter (epsilon) mitigates the singularity at (r=0). By calculating the Laplacian of this function, we derive:

As anticipated, our smoothed Coulomb potential is generated by a corresponding smeared-out point charge, modeled as a Gaussian with width (epsilon).

Next, we will implement this in Python using the multipoles package:

from multipoles import MultipoleExpansion

# Set up the grid x = y = z = np.linspace(-10, 10, 51) X, Y, Z = np.meshgrid(x, y, z, indexing="ij") R = np.sqrt(X ** 2 + Y ** 2 + Z ** 2)

# Set up the charge distribution rho = np.pi ** -1.5 * np.exp(-R ** 2)

# Calculate the multipole expansion, up to the octopole (l_max=8) mpe = MultipoleExpansion({'discrete': False, 'rho': rho, 'xyz': (X, Y, Z)}, l_max=8)

# Evaluate the multipole expansion on the boundary of our grid: boundary = np.ones_like(rho, dtype=bool) boundary[1:-1, 1:-1, 1:-1] = False

phi = np.zeros_like(rho) phi[boundary] = mpe[boundary]

Now that we've established our boundary, we merely need to compute the interior using the Jacobi method, which we will tackle next.

Bringing It All Together

Having developed our Jacobi step routine and set the boundary conditions via the multipole expansion, we can now integrate everything and iterate until convergence:

eps = 1.E-6 max_iter = 10000 h = x[1] - x[0]

for it in range(max_iter):

phi_new = jacobi_step(phi.copy(), rho, h)

diff = np.max(np.abs(phi_new - phi))

if diff < eps:

break

phi = phi_new

it += 1

And that's the complete process!

If you found this article helpful, consider becoming a Medium member for unlimited access to articles. Signing up through this link also supports me as a writer.