Dynamics and Equilibria of N Point Charges on a 2D Ellipse or a 3D Ellipsoid

Barah Makhdum, Ali Nadim

Institute of Mathematical Sciences, Claremont Graduate University, Claremont, USA.

DOI: 10.4236/am.2023.144015                

Abstract

We consider the so-called Thomson problem which refers to finding the equilibrium distribution of a finite number of mutually repelling point charges on the surface of a sphere, but for the case where the sphere is replaced by a spheroid or ellipsoid. To get started, we first consider the problem in two dimensions, with point charges on circles (for which the equilibrium distribution is intuitively obvious) and ellipses. We then generalize the approach to the three-dimensional case of an ellipsoid. The method we use is to begin with a random distribution of charges on the surface and allow each point charge to move tangentially to the surface due to the sum of all Coulomb forces it feels from the other charges. Deriving the proper equations of motion requires using a projection operator to project the total force on each point charge onto the tangent plane of the surface. The position vectors then evolve and find their final equilibrium distribution naturally. For the case of ellipses and ellipsoids or spheroids, we find that multiple distinct equilibria are possible for certain numbers of charges, depending on the starting conditions. We characterize these based on their total potential energies. Some of the equilibria found turn out to represent local minima in the potential energy landscape, while others represent the global minimum. We devise a method based on comparing the moment-of-inertia tensors of the final configurations to distinguish them from one another.

Keywords

Thomson Problem, Point Charges, Ellipsoid, Constrained Optimization

1. Introduction

In this contribution, we consider the so-called Thomson Problem on an ellipsoidal surface. What is the Thomson Problem? Nobel laureate J. J. Thomson who was credited with the discovery of electrons was interested in the structure of the atom and considered the problem of distribution of point charges on the surface of a sphere [1] [2] [3] ; problems of this type have since been referred to as the Thomson Problem. Prominent mathematician and Fields Medal winner Steve Smale included this as one of 18 problems worthy of further study during the current century. Given that N electrons are constrained to lie on a unit sphere, the goal was to find the equilibrium configuration having the minimum electrostatic potential energy labeled $U\left(N\right)$ [4] - [8] . From Coulomb’s Law which states that the potential energy of a pair of charges is proportional to the reciprocal of their distance: $U_{ij}=e_i{e}_j/\left(4\pi {\epsilon }_o{r}_{ij}\right)$ , assuming identical charges and scaling the energy with the pre-factor, the energy of N point charges takes the form

$U\left(N\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\displaystyle \underset{j=i+1}{\overset{N}{\sum }}}\frac{1}{r_{ij}}.$ (1)

The global minimum of $U\left(N\right)$ over all possible collections of N distinct points is typically found by some numerical minimization algorithm. For instance, in the Thomson Problem for two electrons, the solution is obtained when they are as far apart as possible across a diameter of the sphere of unit radius, so that $r_{12}=2$ or $U\left(2\right)=1/2$ .

The minimum energy configurations are known for different N. For the first few cases, the configurations are as follows: When $N=2$ , the electrons will be at antipodal points, across a great diameter of the sphere; for $N=3$ , the electrons form the vertices of an equilateral triangle around a great circle; for N ranging from 4 through 7, the vertices respectively form a regular tetrahedron, a triangular bipyramid, a regular octahedron and a pentagonal bipyramid.

Much work has been done on the Thomson Problem as a testing ground for various optimization algorithms [9] - [16] . Also, it is interesting to note that some physical applications, such as an atom model used in Quantum Mechanics, appeared after Thomson. For instance, a recent reference [17] generalizes the classical Thomson problem to the quantum regime. Other examples include the structure of viruses [18] [19] , the study of topological defects [20] - [25] , and crystals on curved surfaces [24] . In addition, Millikan’s oil drop experiment which successfully measured the elementary charge of an electron could be thought of as being related to the Thomson Problem, since the few discrete charges on each oil drop likely distribute themselves on the surface of the spherical oil drop [26] . A notable recent work [27] generalizes the Thomson Problem to curves in two dimensions, such as an ellipse. They find the minimum energy state of N point charges in that case and show possible multiplicity in the solutions. We tackle this same problem and the corresponding 3D problem with an ellipsoid using a different methodology. Another recent reference [28] considers another variation on this problem.

The significance of the present contribution is that it enables one to generalize the Thomson Problem to surfaces other than a sphere, and it also provides the dynamics, i.e., time evolution of the positions, instead of focusing on just the final equilibrium configuration that minimizes the energy. In the process, other equilibrium states that represent local minima in energy, rather than the global minimum, are also identified and their multiplicities explored.

In the following sections, we formulate the equations for the time evolution of the positions of a finite number of point charges that are initially placed at random on a curve or surface of a given constant shape. We start in two dimensions with the case of a circle and an ellipse, followed by the three-dimensional case of an ellipsoid. In each case, the problem can also be regarded as a constrained optimization where we seek to find the minimum Coulombic energy

${\displaystyle {\sum }_{i=1}^N}{\displaystyle {\sum }_{j=i+1}^N}{\left(r_{ij}\right)}^{-1}$ , subject to the constraint that the particles must remain on the curve or surface of interest. While taking the gradient of the potential energy function with respect to the coordinates of the particles is straightforward, yielding the force on the particles, the constraint that the particles stay on the given surface makes it somewhat challenging to find the right evolution equations. We show below that there are two approaches to obtaining the appropriate evolution equations. Conceptually, they are as follows. Suppose the curve or surface that the particles must reside upon has unit normal $\hat{n}$ . If the force on a particular particle, obtained by taking the negative gradient of the potential energy with respect to the particle position $x$ , is denoted by $F$ , in the absence of any constraints, the particle might move according to the evolution equation $\stackrel{\dot{}}{x}=F$ , taking the particle velocity to be in the direction of the force exerted on it. However, in order to satisfy the constraint that the particle stays on the surface of interest, one would have to modify that equation to: $\stackrel{\dot{}}{x}=F-\hat{n}\hat{n}\cdot F=\left(I-\hat{n}\hat{n}\right)\cdot F$ . That is, we subtract from the force any component it might have had in the direction of the local normal, and keep only the components of the force tangent to the surface. Here, tensor $\left(I-\hat{n}\hat{n}\right)$ is a projection operator that projects vector force $F$ onto the surface. To obtain the equations of motion, one can either directly use this vector evolution equation, or take its normal and tangential components that read: $\hat{n}\cdot \stackrel{\dot{}}{x}=0$ and $\hat{t}\cdot \stackrel{\dot{}}{x}=\hat{t}\cdot F$ , where $\hat{t}$ is the unit tangent to the surface. For a curve in two dimensions, there is just one unit tangent at each location. For a surface in three dimensions, there are two such unit tangents and the component along each of them should be taken. We show and compare these two formulations explicitly for the case of an ellipse but will use the simpler one in three dimensions for the case of an ellipsoid. In each case, we evolve the positions of the particles until they naturally find their local equilibrium state.

2. Point Charges on a Circle

We begin with the simple case of just two point charges. We take the first one to have position vector $x_1$ and the second $x_2$ . Both are arbitrarily placed around a circle of scaled radius 1. The angle between the x-axis and vector $x_1$ is called ${\theta }_1$ and that with $x_2$ is called ${\theta }_2$ as illustrated in Figure 1. There is a force on ion 1 by ion 2 which we call $F_{21}$ ; that on ion 2 by ion 1 is $F_{12}$ where

Figure 1. A schematic of two particles on a circle.

$F_{12}=-F_{21}$ .

The position vectors of the two particles are given by: $x_1=\hat{i}\cos \left({\theta }_1\right)+\hat{j}\sin \left({\theta }_1\right)$ and $x_2=\hat{i}\cos \left({\theta }_2\right)+\hat{j}\sin \left({\theta }_2\right)$ . As such the vector $r_{12}$ connecting $x_1$ to $x_2$ is given by $r_{12}=\hat{i}\left[\cos \left({\theta }_2\right)-\cos \left({\theta }_1\right)\right]+\hat{j}\left[\sin \left({\theta }_2\right)-\sin \left({\theta }_1\right)\right]$ , whose magnitude is $\Vert r_{12}\Vert =\sqrt{{\left(\cos \left({\theta }_2\right)-\cos \left({\theta }_1\right)\right)}^2+{\left(\sin \left({\theta }_2\right)-\sin \left({\theta }_1\right)\right)}^2}$ . The scaled Coulomb force on ion 2 by ion 1 is $F_{12}=r_{12}/{\Vert r_{12}\Vert }^3$ which reduces to

$F_{12}=\frac{\hat{i}\left(\cos \left({\theta }_2\right)-\cos \left({\theta }_1\right)\right)+\hat{j}\left(\sin \left({\theta }_2\right)-\sin \left({\theta }_1\right)\right)}{{\left[{\left(\cos \left({\theta }_2\right)-\cos \left({\theta }_1\right)\right)}^2+{\left(\sin \left({\theta }_2\right)-\sin \left({\theta }_1\right)\right)}^2\right]}^{3/2}}\mathrm{.}$ (2)

In this work, we aim to find the equilibrium distribution of point charges on a surface by allowing an initially random set of points on the surface to relax to their equilibrium under the action of the sum of forces acting on each of them. In particular, if the position vector of one of the charges is $x\left(t\right)$ at time t, we allow it to move with a velocity $\stackrel{\dot{}}{x}$ proportional to the force acting on it. By scaling the time variable, the proportionality constant can always be taken to be unity so that $\stackrel{\dot{}}{x}=F$ with $F$ representing the sum of all forces on the point charge of interest. However, since we wish to restrict the points to remain on the surface of interest, we need to modify that equation to the form: $\stackrel{\dot{}}{x}=\left(I-\hat{n}\hat{n}\right)\cdot F$ . Here $\hat{n}$ is the unit normal to the surface at the point $x$ , $I$ is the unit (identity) tensor, and the projection operator $\left(I-\hat{n}\hat{n}\right)$ , when dotted with any vector, removes any component of that vector in the direction normal to the surface and leaves only the component tangential to the surface. As such, the velocity of the point charge is purely tangential to the surface as well and the evolution equation maintains the particles on the desired surface.

With those in mind, the equations of motion for the two ions can be written as

${\stackrel{\dot{}}{x}}_1=\left(I-{\hat{n}}_1{\hat{n}}_1\right)\cdot F_{21}={\hat{t}}_1{\hat{t}}_1\cdot F_{21}\mathrm{,}$ (3)

${\stackrel{\dot{}}{x}}_2=\left(I-{\hat{n}}_2{\hat{n}}_2\right)\cdot F_{12}={\hat{t}}_2{\hat{t}}_2\cdot F_{12}\mathrm{,}$ (4)

where ${\hat{n}}_1$ and ${\hat{t}}_1$ are the normal and tangential unit vectors to the curve at $x_1$ , namely:

${\hat{n}}_1={{\hat{e}}_r|}_{{\theta }_1}=\cos \left({\theta }_1\right)\hat{i}+\sin \left({\theta }_1\right)\hat{j}\mathrm{,}{\hat{t}}_1={{\hat{e}}_{\theta }|}_{{\theta }_1}=\cos \left({\theta }_1\right)\hat{j}-\sin \left({\theta }_1\right)\hat{i}\mathrm{,}$

with similar expressions at ${\theta }_2$ . The unit tensor $I$ in two dimensions is given by $I=\hat{i}\hat{i}+\hat{j}\hat{j}=\hat{n}\hat{n}+\hat{t}\hat{t}$ .

Using the above results, the tangential component of the force on ion 2, for instance, simplifies to:

${\hat{t}}_2\cdot F_{12}=\frac{\sin \left({\theta }_1-{\theta }_2\right)}{{\left[{\left(\cos \left({\theta }_2\right)-\cos \left({\theta }_1\right)\right)}^2+{\left(\sin \left({\theta }_2\right)-\sin \left({\theta }_1\right)\right)}^2\right]}^{3/2}}\equiv F\left({\theta }_1\mathrm{,}{\theta }_2\right)\mathrm{.}$ (5)

Here, we have defined a scalar function of two variables $F\left({\theta }_1\mathrm{,}{\theta }_2\right)$ which will make it easier to write the equations of motion for point charges on a circle. It represents the tangential component (in the direction of increasing $\theta$ ) of the force by particle 1 (first argument of the function) on particle 2 (second argument).

Since the particles move tangent to the circle, we have ${\stackrel{\dot{}}{x}}_1={\hat{t}}_1{\stackrel{\dot{}}{\theta }}_1$ and ${\stackrel{\dot{}}{x}}_2={\hat{t}}_2{\stackrel{\dot{}}{\theta }}_2$ . The above equations of motion therefore simplify to a system of two equations:

${\stackrel{\dot{}}{\theta }}_1\mathrm{<mo>\; =\; </mo>}F\left({\theta }_2\mathrm{,}{\theta }_1\right)\mathrm{,}$ (6)

${\stackrel{\dot{}}{\theta }}_2=F\left({\theta }_1,{\theta }_2\right).$ (7)

It is also straightforward to obtain the equations of motion for three or more particles by summing the pairwise forces. For instance, for the case of three particles we will have:

${\stackrel{\dot{}}{\theta }}_1=F\left({\theta }_2,{\theta }_1\right)+F\left({\theta }_3,{\theta }_1\right),$ (8)

${\stackrel{\dot{}}{\theta }}_2=F\left({\theta }_1,{\theta }_2\right)+F\left({\theta }_3,{\theta }_2\right),$ (9)

${\stackrel{\dot{}}{\theta }}_3=F\left({\theta }_1,{\theta }_3\right)+F\left({\theta }_2,{\theta }_3\right).$ (10)

More generally, with N particles numbered from 1 to N, the equation of motion for the i-th particle is given by:

${\stackrel{\dot{}}{\theta }}_i={\displaystyle \underset{j=1,j\ne i}{\overset{N}{\sum }}}F\left({\theta }_j,{\theta }_i\right).$ (11)

We use MATLAB to integrate these systems of equations for two or more particles starting with random initial positions of the point charges. As expected intuitively, in all cases the N point charges move in response to the repulsive forces from all their neighbors until they reach an equilibrium state in which they are uniformly spaced on the circle. To show just one illustration of this, in Figure 2 we show the dynamics of 5 point charges on a circle. To be able to track the dynamics in time, we use polar plots in which the radial coordinate represents

Figure 2. Trajectories of five particles on a circle. (a) Long-time equilibrium; (b) Initial transient.

the time variable rather than the radius of the circle which is unity. The left panel in the figure shows the dynamics over a longer time of up to 50. It is clear that the five charges have reached and stay at the vertices of a regular pentagon. The right panel, on the other hand, shows the shorter time dynamics up to a time of 0.2, during which the five point charges gradually move apart along the unit circle due to the repulsive forces between them, creating curved paths in our polar plot.

3. Point Charges on an Ellipse

While the equilibrium position of repelling point charges on a circle can be trivially guessed to be an arrangement that distributes the charges along the circle with equal arc segments between them, the situation for point charges on an elliptical curve is not so straightforward. In order to study this problem, let’s consider an ellipse which is described by the equation

$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1.$ (12)

Without loss of generality, we can take either parameter a or b as our length scale (in effect taking it to be unity) and interpret the other as the dimensionless ratio of the semi-axis lengths of the ellipse. In the following, we take $b=1$ and allow a to be a free parameter. When $a=1$ we recover the circle shape; when $a>1$ we will have an ellipse whose length is in the horizontal direction is greater than its height in the vertical direction.

While it is possible to parameterize the position of a point charge on the ellipse by one angle variable (i.e., one degree of freedom), it turns out to be easier to use the Cartesian coordinates of each point (two degrees of freedom) to track the particles. The two coordinates are constrained to lie on the curve described by Equation (12), so in actuality, there is still one degree of freedom per particle. If we start with two particles as before, with their position vectors and velocities given by $x_i\left(t\right)=\hat{i}{x}_i\left(t\right)+\hat{j}{y}_i\left(t\right)$ and ${\stackrel{\dot{}}{x}}_i=\hat{i}{\stackrel{\dot{}}{x}}_i+\hat{j}{\stackrel{\dot{}}{y}}_i$ the Coulomb force on ion 1 by ion 2 is:

$f_{21}=\frac{x_1-x_2}{{\left|x_1-x_2\right|}^3}=\frac{\hat{i}\left(x_1-x_2\right)+\hat{j}\left(y_1-y_2\right)}{{\left[{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2\right]}^{3/2}}\mathrm{.}$ (13)

The normal and tangent vectors (not normalized) at any point $\left(x\mathrm{,}y\right)$ on the ellipse (with $b=1$ ) are given by:

$n=\hat{i}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(x/a^2\right)+\hat{j}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}y\mathrm{,}t=-\hat{i}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}y+\hat{j}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(x/a^2\right)\mathrm{.}$ (14)

If needed, these can be normalized by dividing each by $\sqrt{\left(x^2/a^4\right)+y^2}$ .

The equations of motion for the point charges can now be derived as follows. One approach is to take the equation of motion for each particle, say ion 1, which reads ${\stackrel{\dot{}}{x}}_1=\left(I-{\hat{n}}_1{\hat{n}}_1\right)\cdot F_{21}$ , and obtain one equation by taking the tangential component of both side by dot multiplying with vector $t$ ; noting that $t\cdot \left(I-{\hat{n}}_1{\hat{n}}_1\right)=t$ , this yields

$t\cdot {\stackrel{\dot{}}{x}}_1=t\cdot F_{21}$ (15)

which reduces to

$-y_1{\stackrel{\dot{}}{x}}_1+\left(x_1/a^2\right){\stackrel{\dot{}}{y}}_1=\frac{-y_1\left(x_1-x_2\right)+\left(x_1/a^2\right)\left(y_1-y_2\right)}{{\left[{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2\right]}^{3/2}}\mathrm{.}$ (16)

A second equation involving ${\stackrel{\dot{}}{x}}_1$ and ${\stackrel{\dot{}}{y}}_1$ can be obtained by taking the normal component of the equation of motion which reads, $n_1\cdot {\stackrel{\dot{}}{x}}_1=0$ ; this is the same as differentiating Equation (12) for the ellipse with respect to time, yielding

$\frac{x_1^2}{a^2}+y_1^2=0\Rightarrow \frac{x_1}{a^2}{\stackrel{\dot{}}{x}}_1+y_1{\stackrel{\dot{}}{y}}_1=0.$ (17)

A similar pair of equations can be obtained governing the evolution of $x_2\left(t\right)$ and $y_2\left(t\right)$ , together resulting in a system of 4 coupled equations for the motion of the two ions.

Written in matrix form, these 4 equations are given by

$\left[\begin{array}{cccc}-y_1& x_1/a^2& 0& 0\\ x_1/a^2& y_1& 0& 0\\ 0& 0& -y_2& x_2/a^2\\ 0& 0& x_2/a^2& y_2\end{array}\right]\left[\begin{array}{c}{\stackrel{\dot{}}{x}}_1\\ {\stackrel{\dot{}}{y}}_1\\ {\stackrel{\dot{}}{x}}_2\\ {\stackrel{\dot{}}{y}}_2\end{array}\right]=\left[\begin{array}{c}\frac{-y_1\left(x_1-x_2\right)+\left(x_1/a^2\right)\left(y_1-y_2\right)}{{\left[{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2\right]}^{3/2}}\\ 0\\ \frac{-y_2\mathrm{<mo\; stretchy="false">\; (\; </mo>}{x}_2-x_1\mathrm{<mo\; stretchy="false">\; )\; </mo>}+\left(x_2/a^2\right)\left(y_2-y_1\right)}{{\left[{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2\right]}^{3/2}}\\ 0\end{array}\right]\mathrm{.}$ (18)

The approach outlined above derives a pair of equations for the evolution of the Cartesian coordinates of each point charge, one equating the tangential velocity of the particle to the tangential component of the force, and the other ensuring that the particle stays on the ellipse as it moves. The resulting system of ODEs has a coefficient matrix on the left-hand side, multiplying the column vector of velocities.

Alternatively, we can use the two components of the vector equation of motion directly, without first taking its dot product with the local unit normal and tangent. In other words, for two particles, we start with ${\stackrel{\dot{}}{x}}_1=\left(I-{\hat{n}}_1{\hat{n}}_1\right)\cdot F_{21}$ for particle one. The right-hand side can be shown to be

$F_{21}-{\hat{n}}_1\left({\hat{n}}_1\cdot F_{21}\right)=\hat{i}f\left(x_2\mathrm{,}{y}_2\mathrm{<mo>\; ;\; </mo>}{x}_1\mathrm{,}{y}_1\right)+\hat{j}g\left(x_2\mathrm{,}{y}_2\mathrm{<mo>\; ;\; </mo>}{x}_1\mathrm{,}{y}_1\right)\mathrm{.}$ (19)

where

$f\left(r_2\mathrm{,}{r}_1\right)=\frac{a^2{y}_1\left(x_1\left(\left(a^2-1\right)y_1+y_2\right)-a^2{x}_2{y}_1\right)}{\left(a^4{y}_1^2+x_1^2\right){\left({\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2\right)}^{3/2}}\mathrm{,}$ (20)

$g\left(r_2\mathrm{,}{r}_1\right)=\frac{x_1\left(a^2{x}_2{y}_1-x_1\left(\left(a^2-1\right)y_1+y_2\right)\right)}{\left(a^4{y}_1^2+x_1^2\right){\left({\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2\right)}^{3/2}}\mathrm{.}$ (21)

Here, we have defined the functions $f\left(x_2\mathrm{,}{y}_2\mathrm{<mo>\; ;\; </mo>}{x}_1\mathrm{,}{y}_1\right)$ and $g\left(x_2\mathrm{,}{y}_2\mathrm{<mo>\; ;\; </mo>}{x}_1\mathrm{,}{y}_1\right)$ whose arguments we have abbreviated as $f\left(r_2\mathrm{,}{r}_1\right)$ and $g\left(r_2\mathrm{,}{r}_1\right)$ . With their help, the equations of motion for the two particles give rise to the following system of four coupled first-order differential equations:

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_1=f\left(r_2\mathrm{,}{r}_1\right)\\ {\stackrel{\dot{}}{y}}_1=g\left(r_2\mathrm{,}{r}_1\right)\\ {\stackrel{\dot{}}{x}}_2=f\left(r_1\mathrm{,}{r}_2\right)\\ {\stackrel{\dot{}}{y}}_2=g\left(r_1\mathrm{,}{r}_2\right)\end{array}$ (22)

Equation (22) can be obtained from Equation (18) by inverting the coefficient matrix that appears on the left-hand side of the latter. For the case with N point charges, we would have 2N equations. In that case, the equations of motion for particle i take the form:

$(\begin{array}{l}{\stackrel{\dot{}}{x}}_i={\displaystyle {\sum }_{j=1,j\ne i}^{N}f\left(r_j\mathrm{,}{r}_i\right)}\\ {\stackrel{\dot{}}{y}}_i={\displaystyle {\sum }_{j=1,j\ne i}^{N}g\left(r_j\mathrm{,}{r}_i\right)}\end{array}$ (23)

Both approaches yield the same results for the dynamics of point particles on the ellipse.

While tracking the positions of the particles on the curve or surface, it is also helpful to be able to monitor the total potential energy in the system and verify that it reaches a minimum as the system approaches its equilibrium. The electrostatic energy of two point charges is given by $q^2/\left(4\pi {\epsilon }_{o}r\right)$ where r is the distance between the point charges of charge q each and ${\epsilon }_o$ is the permittivity of vacuum. We can scale the energy with the factor out front so that it is simply given by 1/r for two particles. In free space, the minimum of that energy occurs when the two particles are infinitely far apart. When constrained to remain on a circle or sphere, they minimum occurs when they are across a diameter. For a system of N particles, the scaled total potential energy is the sum of the pairwise energies and would be given by Equation (1) where

$r_{ij}={\left[{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2\right]}^{1/2}$ .

Results

To illustrate a typical time evolution, for the case of two point charges and with $a=3$ , in Figure 3 we show successive positions of the two particles (distinguished by their colors) at discrete time steps when their random initial positions start near each other close to the top portion of the ellipse. As expected, they move along the ellipse and after some time, they position themselves near the left and right ends of the ellipse, as far away from each other as possible.

Next, we focus on the final equilibrium configuration of a few point charges as their number is varied. We choose $a=3$ so the ellipse is elongated in the horizontal direction. In Figure 4, we display the equilibrium configurations of 3 to 7 point charges, each starting with random initial positions and evolving to their final rest state which is being displayed. As we see, in many of the cases, the equilibrium configuration is not unique, though some of the configurations represent local minima in energy, rather than being the global minimum energy state of the system.

Figure 4(a) shows the result for $N=3$ . At equilibrium, one particle is at the top of the ellipse, while the other two are near the left and right end points, but slightly below, forming an isosceles triangle. An identical configuration which is flipped upside down relative to this one is also obtained, depending on what initial conditions we begin with. Of course, the final energy minimum would be the same for these two configurations, given by $\text{PE}=0.79$ .

Figure 4(b) and Figure 4(c) show two different final equilibrium configurations of 4 point charges on this ellipse. In the latter configuration which has reflection symmetry about the vertical axis, two of the particles are on the top

Figure 3. Dynamics of two particles on an ellipse in the (x, y)-plane.

Figure 4. Equilibrium configurations of N particles on an ellipse in the (x, y)-plane.

segment of the ellipse, closer to the middle, while the other two are near the right and left ends, but below. Connecting them would yield a regular trapezoid shape. The other configuration does not have the reflection symmetry. In that case connecting the four points produces a rhomboid shape. The potential energy for first case is 1.90 and for the second it is 2.02.

Figure 4(d) and Figure 4(e) show two different final equilibrium configurations of 5 point charges. While both configurations have the reflection symmetry about the vertical, they are different in nature. In both cases, two of the particles are on one side (top or bottom) of the ellipse, and three on the other. But for the three particles on one side, in one case they are toward the middle, and in the other, one is in the middle and two are near the ends. Of course, in all these cases, the flipped versions of these configurations are also possible but they are essentially one of the two configurations displayed ( $\text{PE}=3.9$ and $\text{PE}=3.53$ ).

Figure 4(f) and Figure 4(g) display the equilibrium configurations of 6 point charges. As before, these two configurations are essentially different from each other, while flipped versions of each are also possible ( $\text{PE}=5.66$ and $\text{PE}=5.84$ ).

Finally, Figure 4(h) shows the equilibrium configuration of 7 point charges, three of which end up on one side of the ellipse (top or bottom) and found on the other side. In the side that has four, two are near the ends and two are closer to the middle ( $\text{PE}=8.39$ ).

4. Point Charges on an Ellipsoid

We now apply the same procedures that were described above to the problem of point charges on a three-dimensional ellipsoid. This surface in this case is described by

$\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1.$ (24)

Later, without loss of generality, we choose our length scale so that one of the lengths is unity (i.e., $c=1$ ) and allow the other two parameters a and b to vary. When $a=b$ , the ellipsoid reduces to an axisymmetric spheroid, a special case that we will also consider. When $a=b=1$ , this is a sphere and the original Thomson Problem is recovered.

Proceeding in a manner similar to the ellipse case in two dimensions, we start with the position vectors and their derivatives in three dimensions: $x_i\left(t\right)=\hat{i}{x}_i\left(t\right)+\hat{j}{y}_i\left(t\right)+\hat{k}{z}_i\left(t\right)$ and ${\stackrel{\dot{}}{x}}_i=\hat{i}{\stackrel{\dot{}}{x}}_i+\hat{j}{\stackrel{\dot{}}{y}}_i+\hat{k}{\stackrel{\dot{}}{z}}_i$ . The Coulomb force on ion 1 by ion 2 is:

$F_{21}=\frac{x_1-x_2}{{\left|x_1-x_2\right|}^3}=\frac{\hat{i}\left(x_1-x_2\right)+\hat{j}\left(y_1-y_2\right)+\hat{k}\left(z_1-z_2\right)}{{\left[{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2\right]}^{3/2}}\mathrm{.}$ (25)

Also, the unit normal at any point on the ellipsoid is given by

$\hat{n}=\frac{\hat{i}\left(\frac{x}{a^2}\right)+\hat{j}\left(\frac{y}{b^2}\right)+\hat{k}\left(\frac{z}{c^2}\right)}{{\left[{\left(\frac{x}{a^2}\right)}^2+{\left(\frac{y}{b^2}\right)}^2+{\left(\frac{z}{c^2}\right)}^2\right]}^{1/2}}\mathrm{.}$ (26)

After some careful algebra, the expression for the tangential force by particle 2 on particle 1 becomes:

$F_{21}-{\hat{n}}_1\left({\hat{n}}_1\cdot F_{21}\right)=\hat{i}f\left(r_2\mathrm{,}{r}_1\right)+\hat{j}g\left(r_2\mathrm{,}{r}_1\right)+\hat{k}h\left(r_2\mathrm{,}{r}_1\right)$ (27)

where

$f\left(r_2\mathrm{,}{r}_1\right)=\frac{a^2{x}_1\left(b^4{z}_1\left(z_1\left(a^2-c^2\right)+c^2{z}_2\right)+c^4{y}_1^2\left(a^2-b^2\right)+b^2{c}^4{y}_1{y}_2\right)-a^4{x}_2\left(b^4{z}_1^2+c^4{y}_1^2\right)}{{\left({\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2\right)}^{3/2}\left(a^4\left(b^4{z}_1^2+c^4{y}_1^2\right)+b^4{c}^4{x}_1^2\right)}$

(28)

$g\left(r_2\mathrm{,}{r}_1\right)=\frac{b^2\left(a^4{z}_1\left(y_1\left(z_1\left(b^2-c^2\right)+c^2{z}_2\right)-b^2{y}_2{z}_1\right)-c^4{x}_1^2\left(y_1\left(a^2-b^2\right)+b^2{y}_2\right)+a^2{c}^4{x}_1{x}_2{y}_1\right)}{{\left({\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2\right)}^{3/2}\left(a^4\left(b^4{z}_1^2+c^4{y}_1^2\right)+b^4{c}^4{x}_1^2\right)}$

(29)

$h\left(r_2\mathrm{,}{r}_1\right)=\frac{c^2\left(a^4{y}_1\left(b^2{y}_2{z}_1-y_1\left(z_1\left(b^2-c^2\right)+c^2{z}_2\right)\right)-b^4{x}_1^2\left(z_1\left(a^2-c^2\right)+c^2{z}_2\right)+a^2{b}^4{x}_1{x}_2{z}_1\right)}{{\left({\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2\right)}^{3/2}\left(a^4\left(b^4{z}_1^2+c^4{y}_1^2\right)+b^4{c}^4{x}_1^2\right)}$

(30)

The six equations of motion for the coordinates of the two particles are thus:

$(\begin{array}{l}{\stackrel{\dot{}}{x}}_1=f\left(r_2\mathrm{,}{r}_1\right)\hfill \\ {\stackrel{\dot{}}{y}}_1=g\left(r_2\mathrm{,}{r}_1\right)\hfill \\ {\stackrel{\dot{}}{z}}_1=h\left(r_2\mathrm{,}{r}_1\right)\hfill \\ {\stackrel{\dot{}}{x}}_2=f\left(r_1\mathrm{,}{r}_2\right)\hfill \\ {\stackrel{\dot{}}{y}}_2=g\left(r_1\mathrm{,}{r}_2\right)\hfill \\ {\stackrel{\dot{}}{z}}_2=h\left(r_1\mathrm{,}{r}_2\right)\hfill \end{array}$ (31)

For the general case with N particles, we have:

$(\begin{array}{l}{\stackrel{\dot{}}{x}}_i={\displaystyle {\sum }_{j=1,j\ne i}^{N}f\left(r_j\mathrm{,}{r}_i\right)}\\ {\stackrel{\dot{}}{y}}_i={\displaystyle {\sum }_{j=1,j\ne i}^{N}g\left(r_j\mathrm{,}{r}_i\right)}\\ {\stackrel{\dot{}}{z}}_i={\displaystyle {\sum }_{j=1,j\ne i}^{N}h\left(r_j\mathrm{,}{r}_i\right)}\end{array}$ (32)

with i ranging from 1 to N.

The potential energy for the three-dimensional configurations is again given by Equation (1) with $r_{ij}={\left[{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2+{\left(z_i-z_j\right)}^2\right]}^{1/2}$ .

Results

If we set $a=b=c=1$ , the ellipsoid becomes a sphere of radius unity. This is the original Thomson Problem [29] . For all cases we examined, varying the number of particles, we recovered the known results of the Thomson Problem. For instance, with five particles we obtain the triangular bipyramid analyzed by [30] in detail.

The prolate spheroid is the case where the ellipsoid has $a=b<1$ with $c=1$ . The spheroid is elongated along the z-axis while its horizontal cross sections are circles. We take $a=b=1/3$ and show some of the results as the number of point charges varies.

Figure 5(a) and Figure 5(b) shows two different equilibrium configurations of four particles on the prolate spheroid. Similar to the case of four particles on the ellipse in 2D, these four particles can find two different equilibrium configurations. In that equilibrium, the particles lie in the same plane and connecting them produces a rhomboid or regular trapezoid, like the 2D result shown earlier. For the case of five particles, we found one equilibrium configuration that is shown in Figure 5(c). These particles do not lie in a plane.

Figure 5. Equilibrium configurations on a prolate spheroid in the (x, y, z)-space.

Next we consider an oblate spheroid. This is the case when $a=b>1$ . In this case, the spheroid would appear flattened along the z-axis. In Figure 6, we have shown two equilibrium configurations for the cases $N=3$ and $N=5$ . In the top two panels, we see that at equilibrium three particles can arrange themselves around the equator of the spheroid forming an equilateral triangle, or one particle could find itself at one of the poles, with the other two on the other side of the spheroid but close to the equator, forming an isosceles triangle. Similarly, for five particles, one may find all five equally spaced around the equator, or have a particle at one of the poles and four close to the equator on the other side, forming a pyramid with a square base. The triangular bipyramid that was the shape for five particles on the sphere does not form in this case.

Finally, we consider the more general triaxial ellipsoid in which all three length parameters are different. We still set $c=1$ and take $a=7$ and $b=1/2$ as a test case. With $N=5$ point charges, we find that they arrange themselves in a plane this time, and within that plane, their configuration is similar to one on the ellipse, with two of them on one side and close to the local pole, while the other three are on the other side, with two close to the ends and one at the pole. This is shown in Figure 7.

5. Characterizing the Multiplicity of Equilibrium Configurations

As we saw in the cases of point charges on an ellipse or ellipsoid, there were some cases that had more than one final equilibrium configuration. The potential

Figure 6. Equilibrium configurations on an oblate spheroid in the (x, y, z)-space.

energies of those final configurations could be evaluated to see which one represents the true global minimum and which others represent local minima that are still at equilibrium but did not achieve the global minimum configuration. It is thus apparent that the equilibrium configurations of N particles are not necessarily unique, having different energy levels. This is also observed in references [12] [31] [32] .

To see this more explicitly, let us consider the case of four point charges on an ellipse with parameter values $a=3$ and $b=1$ . Starting with $N=4$ and running the program over a hundred times with random starting positions of the four charges, we plotted the value of the Coulombic energy as a function of time all on the same graph, producing the one in Figure 8. If we focus on the tail end of these energy curves, we notice that they are all asymptotically approaching one of two possible values. The inset in the plot zooms in on the tail end to show the two asymptotes more clearly. The final equilibrium energy levels are nearly 2.02 or 1.90, with a clear gap between the two. So the random starting configurations all converge to one of two possible states whose final energy values are different but close.

It is not always easy to compare the energy levels of the multiple equilibrium configurations since they may be very close to one another, and depending on

Figure 7. Five particles on a triaxial ellipsoid in the (x, y, z)-space.

Figure 8. Time variation of the potential energy of four particles on an ellipse with different initial configurations. The inset shows a blow-up of the tails of the curves.

the tolerance chosen for the numerical integration of the system of differential equations, the differences in energy may be obscured by the accuracy in the numerical solution. As such, we sought other methods to be able to compare two or more final equilibrium configuration or arrangements of the N point charges. One method that proved successful was based on the concept of moment-of-inertial from classical mechanics.

Once the N particles have found their final equilibrium state, in order to calculate the moment-of-inertia tensor (matrix) of that configuration we proceed as follows. First we find the centroid or center-of-mass of the N particles, defined by

$\langle r\rangle =\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{r}_i\mathrm{.}$ (33)

The x-, y- and z-coordinates of the centroid are just the arithmetic means of those of the N particles. The moment-of-inertia tensor is then defined by

$M=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\left(r_i-\langle r\rangle \right)\left(r_i-\langle r\rangle \right)\mathrm{.}$ (34)

In 3D, $M$ is equivalent to a symmetric 3-by-3 matrix. Its first diagonal component, for instance, is given by

$M_{11}=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\left(x_i-\langle x\rangle \right)}^2,$

while the entries $M_{12}$ and $M_{21}$ are given by

$M_{12}=M_{21}=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\left(x_i-\langle x\rangle \right)\left(y_i-\langle y\rangle \right).$

In 2D, the moment-of-inertial matrix is a symmetric 2-by-2 matrix. Given the final configuration of the N particles in either 2D or 3D, we can calculate the moment-of-inertia tensor corresponding to that configuration. If a given arrangement of particles is flipped or rotated as a rigid body, while the entries in the matrix will change, the principal moments of inertia which are the eigenvalues of that symmetric matrix do not change. Thus two configurations that are the same to within a solid-body rotation, or reflection across the vertical or horizontal line, yield identical eigenvalues, while configurations that are genuinely different give us different eigenvalues.

When we carry out this process for the two different 4-particle configurations on the ellipse with $a=3$ , we find the principal moments of inertia to be (4.78, 0.46) for one case, and (5.11, 0.29) in the other, better distinguishing one from the other. We have also applied the same procedure to 3D configurations on an ellipsoid and have been able to distinguish multiple equilibria from one another by comparing their three principal moments of inertia.

6. Conclusions

In this paper, we have provided an approach to finding the equilibrium configuration of N identical point charges that are constrained to reside on a particular surface. In 2D, we considered an ellipse and in 3D an ellipsoid. The latter generalizes the Thomson Problem from a sphere to the more general ellipsoid or spheroid shapes.

Our approach is to follow the trajectories of the N points as each moves under the action of the repulsive forces from all the other points until they find their equilibrium, while maintaining the constraint that they should remain on the surface. This necessitates removing any component of the net force in the direction of the unit normal to the surface while writing the equations of motion. We presented two different approaches to achieve this. While the two formulations provide systems of differential equations that look different, they can be reduced to one another and the resulting dynamics are the same.

By applying this method to cases of a few point charges on both ellipses and ellipsoids, we were able to obtain the equilibrium configurations. For the case of a circle, as expected the final equilibrium is one where the particles are equally spaced along the boundary of the circle. For the case of a sphere, the known configurations for the Thomson Problem were recovered. But for the ellipse in 2D and the ellipsoid in 3D more interesting and complicated results emerged. In both cases, it was possible to find multiple different equilibrium configurations for certain values of N, the number of charges, depending on the starting random positions. We devised a method to compare the configurations based on both their final Coulombic energy and based on the moment-of-inertia tensors associated with their final state. When multiple equilibria are possible, some have the lowest energy, representing the global minimum of the Coulombic energy. Others have slightly higher energy but represent a local minimum in configuration space. Finding the moment-of-inertia tensor and comparing the principal moments of inertia (eigenvalues of that symmetric tensor) can be used to distinguish the final equilibrium configurations, even when their energy levels are relatively close.

One of the advantages of the approach introduced in this paper is having an explicit system of evolution equations that can begin with any given, possibly random, set of initial positions of the particles and allow the system to evolve on its own to obtain the final equilibrium configuration. With reasonable numbers of point charges, the computations are fast and reach the final steady state quickly. Smooth surface shapes such as spheroids and ellipsoids can be explored easily. Some limitations of the approach are that as the number of particles approaches many hundreds or thousands, the computations become slow and the system becomes stiff. In addition, only smooth surface geometries can be explored which precludes surfaces that have kinks or cusps. One of the promising directions for future research is to allow the particles to move into the bulk when their number density on the surface becomes too large. In that case, lower overall energy may be achieved when some of the particles move inside the domain rather than being forced to stay on the boundary.

Acknowledgements

Barah Makhdum is grateful for financial support for her graduate studies from her country Saudi Arabia and King Salman bin Abdulaziz Al Saud.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Thomson, J.J. (1913) The Structure of the Atom. Academie Royale de Belgique, Brussels.
[2]	Thomson, J.J. (1904) XXIV. On the Structure of the Atom: An Investigation of the Stability and Periods of Oscillation of a Number of Corpuscles Arranged at Equal Intervals around the Circumference of a Circle; with Application of the Results to the Theory of Atomic Structure. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 7, 237-265. https://doi.org/10.1080/14786440409463107
[3]	Fadeev, S.A., Dedok, V.A. and Bondarenko, A.N. (2022) Polynomial Classification Algorithm for Solutions of the Thomson Problem. Journal of Applied and Industrial Mathematics, 16, 189-202. https://doi.org/10.1134/S1990478922020028
[4]	Cohn, H. (1956) Stability Configurations of Electrons on a Sphere. Mathematics of Computation, 10, 117-120. https://doi.org/10.1090/S0025-5718-1956-0081133-0
[5]	Ballinger, B., Blekherman, G., Cohn, H., Giansiracusa, N., Kelly, E. and Schürmann, A. (2009) Experimental Study of Energy-Minimizing Point Configurations on Spheres. Experimental Mathematics, 18, 257-283. https://doi.org/10.1080/10586458.2009.10129052
[6]	Cohn, H. and Kumar, A. (2007) Universally Optimal Distribution of Points on Spheres. Journal of the American Mathematical Society, 20, 99-148. https://doi.org/10.1090/S0894-0347-06-00546-7
[7]	Wales, D.J. and Ulker, S. (2006) Structure and Dynamics of Spherical Crystals Characterized for the Thomson Problem. Physical Review B, 74, Article ID: 212101. https://doi.org/10.1103/PhysRevB.74.212101
[8]	Calkin, M.G., Kiang, D. and Tindall, D.A. (1986) Minimum Energy Configurations. Nature, 319, 454-454. https://doi.org/10.1038/319454b0
[9]	Erber, T. and Hockney, G.M. (1991) Equilibrium Configurations of N Equal Charges on a Sphere. Journal of Physics A: Mathematical and General, 24, L1369-L1377. https://doi.org/10.1088/0305-4470/24/23/008
[10]	Altschuler, E.L., Williams, T.J., Ratner, E.R., Tipton, R., Stong, R., Dowla, F. and Wooten, F. (1997) Possible Global Minimum Lattice Configurations for Thomson’s Problem of Charges on a Sphere. Physical Review Letters, 78, 2681-2685. https://doi.org/10.1103/PhysRevLett.78.2681
[11]	Wales, D.J., McKay, H. and Altschuler, E.L. (2009) Defect Motifs for Spherical Topologies. Physical Review B, 79, Article ID: 224115. https://doi.org/10.1103/PhysRevB.79.224115
[12]	Bondarenko, A.N., Bugueva, T.V. and Kozinkin, L.A. (2016) Numerical Study of the Structure of Metastable Configurations for the Thomson Problem. Russian Physics Journal, 59, 121-129. https://doi.org/10.1007/s11182-016-0746-3
[13]	Birtea, P. and Comănescu, D. (2017) Newton Algorithm on Constraint Manifolds and the 5-Electron Thomson Problem. Journal of Optimization Theory and Applications, 173, 563-583. https://doi.org/10.1007/s10957-016-1049-0
[14]	Morris, J.R., Deaven, D.M. and Ho, K.M. (1996) Genetic-Algorithm Energy Minimization for Point Charges on a Sphere. Physical Review B, 53, R1740-R1743. https://doi.org/10.1103/PhysRevB.53.R1740
[15]	Von Brecht, J.H., Uminsky, D., Kolokolnikov, T. and Bertozzi, A.L. (2012) Predicting Pattern Formation in Particle Interactions. Mathematical Models and Methods in Applied Sciences, 22, Article ID: 1140002. https://doi.org/10.1142/S0218202511400021
[16]	Ashby, N. and Brittin, W.E. (1986) Thomson’s Problem. American Journal of Physics, 54, 776-777. https://doi.org/10.1119/1.14440
[17]	Yang, L. and Yao, Z. (2018) Two and Three Electrons on a Sphere: A Generalized Thomson Problem. Physical Review B, 97, Article ID: 235431. https://doi.org/10.1103/PhysRevB.97.235431
[18]	Zandi, R., Reguera, D., Bruinsma, R.F., Gelbart, W.M. and Rudnick, J. (2004) Origin of Icosahedral Symmetry in Viruses. Proceedings of the National Academy of Sciences of the United States of America, 101, 15556-15560. https://doi.org/10.1073/pnas.0405844101
[19]	Longuet-Higgins, M.S.(2009) Snub Polyhedra and Organic Growth. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 465, 477-491. https://doi.org/10.1098/rspa.2008.0219
[20]	Pérez-Garrido, A., Dodgson, M.J.W. and Moore, M.A. (1997) Influence of Dislocations in Thomson’s Problem. Physical Review B, 56, 3640-3643. https://doi.org/10.1103/PhysRevB.56.3640
[21]	Bowick, M., Cacciuto, A., Nelson, D.R. and Travesset, A. (2002) Crystalline Order on a Sphere and the Generalized Thomson Problem. Physical Review Letters, 89, Article ID: 185502. https://doi.org/10.1103/PhysRevLett.89.185502
[22]	Bausch, A.R., Bowick, M.J., Cacciuto, A., Dinsmore, A.D., Hsu, M.F., Nelson, D.R., Nikolaides, M.G., Travesset, A. and Weitz, F.A. (2003) Grain Boundary Scars and Spherical Crystallography. Science, 299, 1716-1718. https://doi.org/10.1126/science.1081160
[23]	Backofen, R., Voigt, A. and Witkowski, T. (2010) Particles on Curved Surfaces: A Dynamic Approach by a Phase-Field-Crystal Model. Physical Review E, 81, Article ID: 025701. https://doi.org/10.1103/PhysRevE.81.025701
[24]	Irvine, W.T.M., Vitelli, V. and Chaikin, P.M. (2010) Pleats in Crystals on Curved Surfaces. Nature, 468, 947-951. https://doi.org/10.1038/nature09620
[25]	Azadi, A. and Grason, G.M. (2014) Emergent Structure of Multidislocation Ground States in Curved Crystals. Physical Review Letters, 112, Article ID: 225502. https://doi.org/10.1103/PhysRevLett.112.225502
[26]	Jones, R.C. (1995) The Millikan Oil-Drop Experiment: Making It Worthwhile. American Journal of Physics, 63, 970-977. https://doi.org/10.1119/1.18001
[27]	Amore, P. and Jacobo, M. (2019) Thomson Problem in One Dimension: Minimal Energy Configurations of N Charges on a Curve. Physica A: Statistical Mechanics and Its Applications, 519, 256-266. https://doi.org/10.1016/j.physa.2018.12.040
[28]	Leng, T. and Wu, Y.C. (2023) A Reverse Thomson Problem on the Unit Circle. Proceedings of the American Mathematical Society, 151, 327-337. https://doi.org/10.1090/proc/16110
[29]	Whyte, L.L. (1952) Unique Arrangements of Points on a Sphere. The American Mathematical Monthly, 59, 606-611. https://doi.org/10.1080/00029890.1952.11988207
[30]	Schwartz, R.E. (2010) The 5 Electron Case of Thomson’s Problem. ArXiv: 1001.3702
[31]	Calef, M., Griffiths, W. and Schulz, A. (2015) Estimating the Number of Stable Configurations for the Generalized Thomson Problem. Journal of Statistical Physics, 160, 239-253. https://doi.org/10.1007/s10955-015-1245-6
[32]	Saff, E.B. and Kuijlaars, A.B.J. (1997) Distributing Many Points on a Sphere. The Mathematical Intelligencer, 19, 5-11. https://doi.org/10.1007/BF03024331