Thomson Problem

The Nobel Laureate physicist J. J. Thomson posed a problem in 1904 of determining the minimum electrostatic energy (and configuration) of \(N_{e}\) electrons constrained on the surface of a unit sphere, each electron repelling another from Coulomb's law. In reality, Thomson's model of the atom became outdated and Earnshaw's theorem suggests that no system of discrete electric charges can be held in stable equilibrium under the influence of their electrical interaction alone (Aspden 1987). Regardless, the Thomson problem is still related to spherical codes, sphere packing, computational geometry (for spheres of higher dimensions), and especially to #7 of Smale's list of unsolved problems in mathematics, and more. The Thomson problem presents itself as a non-convex optimization problem, and remains to be generally unsolved for all values of \(N_{e}>1\).

To frame the optimization problem, we will arrive at the objective function as a function of \(N_{e}\). Arriving at Coulomb's law, using Gauss's law of Maxwell's equations, in SI-units, the net outflow of the electric field \(\vec{E}\) through a closed surface (sphere), \(S\), with volume \(V\), around a charge \(q\), where the field travels towards negative charges (and radially outward for a positive charge), is proportional to the enclosed charge, consisting of a charge density across the enclosed volume, \(\rho\left( \vec{r} \right)\). Using permittivity of free space (\(\epsilon_{0}\)) as the coefficient of proportion,

\[ \begin{align} \vec{\nabla}_{r}\cdot\vec{E}\left(\vec{r}\right) = \frac{\rho\left(\vec{r}\right)}{\epsilon_{0}} \Rightarrow \iiint_{V} \left(\vec{\nabla}_{r}\cdot\vec{E}\left(\vec{r}'\right)\right) d^{3} \vec{r}’ = \frac{1}{\epsilon_{0}}\iiint_{V} \rho\left(\vec{r}'\right) d^{3} \vec{r}'\ = \frac{q}{\epsilon_{0}}.\label{eq:1} \end{align} \]

Using the divergence (or Gauss/Green-Ostrogradski) theorem, for a general field \(\vec{F}\),

\[ \begin{align} \iiint_{V} \left( \vec{\nabla} \cdot \vec{F}\right) dV = {\large\bigcirc}\kern-1.55em\iint_{S} \left(\vec{F}\cdot \hat{n}\right)dS.\label{eq:2} \end{align} \]

From eqs. \(\eqref{eq:1}\) and \(\eqref{eq:2}\), using spherical symmetry, and \(\vec{E}\) having a radial component pointing away from \(q\),

\[ \begin{align} {\large\bigcirc}\kern-1.55em\iint_{S} \left(\vec{E}\left(\vec{r}'\right)\cdot \hat{r}'\right) d^{2} \vec{r}’ = \frac{q}{\epsilon_{0}} \Rightarrow 4\pi r^{2} \vec{E}\left(\vec{r}\right) \cdot \hat{r}= \frac{q}{\epsilon_{0}} \Rightarrow \vec{E}\left(\vec{r}\right) = \frac{1}{4\pi\epsilon_{0}}\frac{q}{r^{2}} \hat{r} = k\frac{q}{r^{2}}\hat{r}\label{eq:3} \end{align} \]

where \(k = \frac{1}{4\pi\epsilon_{0}}\). The electrostatic potential energy \(U\left(\vec{r}\right)\) of two-point charge system \(\left(q, q'\right)\) is the work required to assemble the system of these two charges by bringing them close together, usually considered to be an infinite distance away (without undergoing any acceleration),

\[ \begin{align} U\left(\vec{r}\right) = \lim_{a\to\infty} \left(-\int_{a}^{r}\vec{F}\left(\vec{r}'\right)\cdot d\vec{r}'\right) = \lim_{a\to\infty} \left(-\int_{a}^{r}q'\vec{E}\left(\vec{r}'\right)\cdot d\vec{r}'\right) .\label{eq:4} \end{align} \]

From \(\eqref{eq:3}\) and \(\eqref{eq:4}\),

\[ \begin{align} U\left(\vec{r}\right) = \lim_{a\to\infty} \left(k \frac{qq’}{r’} \right)_{a}^{r} = k\frac{qq’}{r}. \label{eq:5} \end{align} \]

Likewise, the change in electrostatic potential energy of a system of \(N_{e}\) point charges is the work required to assemble the system of charges by bringing them close together from an infinite distance away, without undergoing acceleration. For a particle pair \(i\), \(j\), \(U_{ij} = k\frac{q_{i}q_{j}}{r_{ij}}\) where \(r_{ij} = \| \vec{r}_{i} - \vec{r}_{j} \|_{2}\). Without loss of generality, we shall ignore the magnitude of \(k\) and \(q_{i}\) \(\forall i\) \(\in\) \(\left(1,..,N_{e}\right)\), and arrive at the total electrostatic potential energy, \(\widetilde{U}\left(N_{e},\vec{r}_{1},…,\vec{r}_{N_{e}}\right)\), for the system of \(N_{e}\) point charges for a numerical minimization framework:

\[ \begin{align} \widetilde{U}\left(N_{e},\vec{r}_{1},…,\vec{r}_{N_{e}}\right) = \displaystyle\sum_{i<j<N_{e}} \frac{1}{r_{ij}}, i \in \left(1,..,j\right),j \in \left(1,..,N_{e}\right).\label{eq:6} \end{align} \]

Thus, from \(\eqref{eq:6}\) the objective function for the Thomson problem is:

\[ \begin{equation} \begin{aligned} & \underset{\vec{r}_{1}, …, \vec{r}_{N_{e}}}{\text{minimize}} & & \widetilde{U}\left(N_{e},\vec{r}_{1},…,\vec{r}_{N_{e}}\right) \\ & \text{subject to} & & \| \vec{r}_{i} \|_{2} = 1, \forall i \in \left(1,…,N_{e}\right). \end{aligned} \end{equation} \]

Minimizing the electrostatic energy across \(N_{e}\) electrons constrained on a sphere, presents a non-convex optimization problem where the number of local minima appears to grow exponentially with \(N_{e}\) (Laszlo, 2022). Here, I use a stochastic pattern search global optimizer (SP-GOPT Solver) for finding solutions for varying \(N_{e}\). I benchmark my results with those published already, and I attempt to pioneer for minimizing \(\widetilde{U}\left(N_{e},\vec{r}_{1},…,\vec{r}_{N_{e}}\right)\) for larger values of \(N_{e}\) than what I can publicly find (with the help of distributed optimization using MPI).

Numerical solutions (global minima) and configurations, of smallest electrostatic energies \(N_{e}>1\)

For \(N_{e} \le 100\), SP-GOPT is used to calculate \(\widetilde{U}_{\text{min}}\left(N_{e}\right)\), called \(\widetilde{U}_{\text{min}}^{\text{SP-GOPT}}\left(N_{e}\right)\), and is compared against already-published results, \(\widetilde{U}_{\text{min}}^{\text{benchmark}}\left(N_{e}\right)\), for benchmarking purposes. For the optimization routine, the initial evaluation is such that \(\left( \vec{r}_{1},…,\vec{r}_{N_{e}}\right)\) are on a single point, and a re-start (for premature terminations) uses the previous start output as the new initial set of \(\left( \vec{r}_{1},…,\vec{r}_{N_{e}}\right)\). The optimal coordinates, denoted as \(\left( \vec{r}_{1}^{*},…,\vec{r}_{N_{e}}^{*}\right)\), are subject to being rotated about a unit sphere. The results in green are smaller than published results.

For large \(N_{e}\), a distributed computing approach for SP-GOPT is used, to propose newly unpublished results for large \(N_{e}\).

This is a great problem for confronting the challenges of non-conex objective functions, not only seen in the field of optimization, but those found in real-world applications like those in Data Science, professional sports betting, quantitative finance, engineering, physics, and more.

Reference

1. Aspden, H. “Earnshaw's Theorem.” Amer. J. Phys. 55, 199-200, 1987.

2. Wikipedia contributors. (2022, October 22). “Thomson problem.” Wikipedia. https://en.wikipedia.org/wiki/Thomson_problem

3. Laszlo, H. “Numerical Solutions for the Thomson Problems.” Draft 0.4, January 20, 2022. https://www.hars.us/Papers/Numerical_Thomson.pdf

4. Wales, D.J., and Ulker, S., “Structure and dynamics of spherical crystals characterized for the Thomson problem.” Phys. Rev. B 74, 212101, 2006. https://journals.aps.org/prb/abstract/10.1103/PhysRevB.74.212101

5. Wales, D.J., McKay, H., and Altschuler, E., “Defect motifs for spherical topologies.” Phys. Rev. B 79, 224115, 2009. https://journals.aps.org/prb/abstract/10.1103/PhysRevB.79.224115