Thomson, One N at a Time

Recently, we published a solution of the Thomson problem for $N = 7$ electrons, showing that seven electrons on a sphere have the lowest energy when they form a regular pentagonal bipyramid: two poles, with a pentagon around the equator. It took ten Claude Sonnet 5.5 agents around 15 hours to come up with a 17,895-line proof of the statement in Lean. This post sheds light on the underlying method and suggests some values of $N$ to tackle next.

The short version: the proofs in this post bound the energy from below by sums that can’t be negative. Bounds built from pairs of points are exact only for the most balanced $N$, while bounds built from triangles of points, plus a case split, are what proved $N = 7$ and $8$. Our numerics suggest the same recipe should reach $N = 10$ then $9$ next, but that $N = 32$ will need a new idea.

1. The problem

Let $e_1, \dots, e_N$ be points on the unit sphere in $\mathbb{R}^3$.¹ They repel like electrons—each pair has electrical potential energy inversely proportional to their distance, which we denote $\operatorname{dist}(e_i, e_j)$. The Thomson problem asks for the configuration of points that minimizes the total electrical potential energy $E = \sum_{i<j} \operatorname{dist}(e_i, e_j)^{-1}$.

Before we prove anything, there’s a simple algorithm for finding good configurations: start from random points and let the charges push each other around until they settle. Adding a dissipative drag force, we see that potential energy decreases over time, guiding the system toward a local minimum.

But how can we prove that the resulting configuration is a global minimum? To do so, we need a global lower bound on the energy, achieved at our proposed configuration. Concretely, we’re looking for a constant $B$ such that $E \ge B$ for all configurations and $E = B$ for the proposed minimal configuration.

2. Some minimal configurations

Here are the answers for $N$ up to 12,² and whether and how each one has been proved. Three of them are regular polyhedra: the tetrahedron, the octahedron and the icosahedron.³ We lean on two throughout: $N = 7$, the subject of our previous post, and $N = 12$, the cleanest case where our method works perfectly.

3. A lower bound from polynomials

Most of the proofs above start from one idea, due to Delsarte and adapted to energy by Yudin. The plan is to find a sum, built from the angles between points, that can never be negative for any configuration, and then use it to bound the energy from below. We state the key fact first, check it by hand in two simple cases, turn it into a bound, and then ask when the bound is exact.

Setup: energy depends only on angles

Take two points $e_i,e_j$ on the unit sphere and let $t_{ij}=\langle e_i,e_j\rangle$ be the cosine of the angle between them. So $t_{ij}=1$ means the same point and $t_{ij}=-1$ means opposite points.

Lemma 1 · energy depends only on angles

The distance between two unit vectors is $\lVert e_i-e_j\rVert=\sqrt{2-2t_{ij}}$. So a pair costs $f(t_{ij})$, where $f(t)=1/\sqrt{2-2t}$, and the total energy is $E=\tfrac12\sum_{i\ne j}f(t_{ij})$.

Proof. Both vectors have length 1, so:

Recall that the electrical potential energy between $e_i$ and $e_j$ is the inverse of dist$(e_i,e_j)$, which is precisely $f(t_{ij})$. Summing over every pair once gives $E=\sum_{i\lt j}f(t_{ij})$. Summing over ordered pairs $(i,j)$ with $i\ne j$ counts each pair twice, because $t_{ij}=t_{ji}$, which gives the factor $\tfrac12$. $\square$

We use the ordered version because sums over all ordered pairs turn into squared lengths, as we’re about to see. From here on, $t_{ij}$ is the only thing we track about a pair of points, and $f$ is the energy curve that everything below is measured against.

3.1 Legendre polynomials

The polynomials $P_0,P_1,P_2,\dots$ are the Legendre polynomials. You get them by running the Gram–Schmidt orthogonalization process on $1,t,t^2,\dots$ over the interval $[-1,1]$ and rescaling so that $P_k(1)=1$. The first few are $P_0=1$, $P_1=t$, $P_2=\tfrac12(3t^2-1)$, $P_3=\tfrac12(5t^3-3t)$.

For each $k$ there is a map $\varphi_k$ from the sphere into a space of dimension $2k+1$ with $P_k(\langle x,y\rangle)=\langle\varphi_k(x),\varphi_k(y)\rangle$. This is the addition theorem for spherical harmonics. Granting it, the bound takes one line. Add $P_k(t_{ij})$ over every ordered pair, including $i=j$, and call the total $S_k$:

A single term can be negative: $P_1(t)=t$ is negative for an obtuse pair. Only the sum over the whole configuration is guaranteed to be non-negative.

But what is this mysterious map $\varphi_k$? The next two sections answer that question for $k=1$ and $k=2$.

3.2 The centre of mass

Let’s start with $k=1$, where $P_1(t)=t$. Add $P_1(t_{ij})$ over every ordered pair, including $i=j$, and call the total $S_1$:

The sum of all the cosines is the squared length of $e_1+\dots+e_N$, which is $N$ times the centre of mass. A squared length can’t be negative, and it is zero exactly when the centre of mass sits at the origin.

The $N$ terms with $i=j$ each equal $t_{ii}=1$, so they contribute $N$. Moving them across gives $\sum_{i\ne j}t_{ij}\ge -N$. In words, the average cosine over the $N(N-1)$ ordered pairs of distinct points is at least $-1/(N-1)$,⁴ so the points can’t all be far apart from each other. For $N=4$ the floor is $-1/3$, and the regular tetrahedron puts every one of its pairs exactly on it.

There is another way to read $S_1=0$ that carries over to larger $k$. Rotate the whole configuration about the origin and the vector $\sum_i e_i$ rotates with it. The only vector that stays put under every rotation is $0$. So $S_1=0$ says the average position of the points doesn’t change when you rotate the configuration.

3.3 The same trick with a matrix

At $k=2$, $P_2(t)=\tfrac12(3t^2-1)$ involves $t^2$, and a squared inner product is a quadratic form in an outer product⁵: $t_{ij}^2=e_i^{\mathsf T}(e_je_j^{\mathsf T})e_i$. So where $k=1$ added up the vectors $e_i$, $k=2$ adds up their outer products into the $3\times3$ matrix $M=\sum_{i=1}^{N}e_ie_i^{\mathsf T}$, whose entries are $M_{ab}=\sum_{i=1}^{N}(e_i)_a(e_i)_b$. Its diagonal entry $M_{aa}$ is the total squared extent of the configuration along axis $a$, and its trace is $\operatorname{tr}M=\sum_i\lVert e_i\rVert^2=N$.

Step 1. Write $\sum t_{ij}^2$ in terms of $M$:

Line 3 uses the identity $v^{\mathsf T}Bv=\operatorname{tr}(B\,vv^{\mathsf T})$.

Step 2. Put that into $S_2$. There are $N^2$ ordered pairs, and each contributes $-\tfrac12$:

Step 3. Complete the square. The Frobenius norm of a matrix, $\lVert A\rVert_F^2=\operatorname{tr}(A^{\mathsf T}A)$, is the sum of its squared entries. Using that $M$ is symmetric, $\operatorname{tr}M=N$ and $\operatorname{tr}I=3$:

Step 4. Compare Steps 2 and 3:

This has the same shape as $k=1$: a sum over all pairs turned into a squared length. It also has the same rotation reading. Rotating the configuration by a rotation $R$ turns $M$ into $RMR^{\mathsf T}$, and the only matrices that survive every rotation unchanged are multiples of the identity. So $S_2=0$ means $M=\tfrac N3I$: the configuration extends equally far along every axis, with no preferred direction at second order.

Two examples we’ll keep using. For the icosahedron ($N=12$), symmetry forces $M=4I=\tfrac{12}{3}I$, so $S_2=0$. For the pentagonal bipyramid ($N=7$), the ring contributes $\tfrac52$ along each of the two horizontal axes and the two poles contribute $2$ along the vertical one, so $M=\operatorname{diag}(\tfrac52,\tfrac52,2)$. That isn’t $\tfrac73I$, and

In summary, we find familiar objects in $\varphi_1(x)=x$ and $\varphi_2(x)=\sqrt{3/2}\,\bigl(xx^{\mathsf T}-\tfrac13I\bigr)$. For larger $k$, the coordinates of $\varphi_k$ are the degree-$k$ spherical harmonics, the sphere’s version of Fourier modes. If you’re familiar with kernels from machine learning: $P_k(\langle x,y\rangle)$ is a positive-definite kernel and $\varphi_k$ is its feature map.

Why $k=0$ is special. The first row is the odd one out. The map $\varphi_0(x)=1$ gives every point the same one-dimensional feature, so $S_0=(1+1+\dots+1)^2=N^2$ for every configuration. It never changes and is never zero. Think of $P_0$ as the average level of a signal, like the constant offset of a sound wave, and of $P_k$ for $k\ge1$ as the wiggles on top. The wiggle sums $S_k$ depend on the configuration, are never negative, and are zero exactly when the configuration is balanced at that level. This difference is why the constant term $a_0$ gets special treatment in the next subsection.

A configuration with $S_1=\dots=S_s=0$ is called a spherical $s$-design: averages of polynomials up to degree $s$ over its points match the averages over the whole sphere. Which $S_k$ vanish for a given configuration will decide whether this method gives the exact answer. We’ll compare $N=12$ and $N=7$ in 3.5.

3.4 The polynomial program

Now combine the $P_k$ into a polynomial $h(t)=\sum_{k=0}^{d}a_kP_k(t)$. Two rules:

• $a_k\ge0$ for every $k\ge1$. This makes $h$ a positive combination of things whose sums are non-negative.
• $h(t)\le f(t)$ for all $t\in[-1,1)$. This keeps $h$ below the energy curve.

Then for any configuration:

The last line doesn’t mention the configuration at all, so it is a lower bound for every configuration at once. This is where $P_0$ is treated differently. The term $a_0S_0=a_0N^2$ is known exactly, so we keep it. The terms $a_kS_k$ with $k\ge1$ are unknown but non-negative, so we throw them away. That is why only $a_1,a_2,\dots$ need a sign condition: we never discard $a_0$‘s term, so there is no inequality for a sign to protect.

Every $P_k(1)=1$, so $h(1)=a_0+a_1+\dots+a_d$, and the bound simplifies to

Each unit of $a_0$ earns $N(N-1)/2$ and each unit of $a_k$ with $k\ge1$ costs $N/2$. So you want $h$ as high as the curve $f$ allows, using as few wiggles as possible. Finding the best such $h$ is a linear program in the coefficients, and the answer is usually called the LP bound.

Polynomial program

Variables: the coefficients $a_0,a_1,\dots,a_d$.

Constraints: $a_k\ge0$ for $k\ge1$, and $\sum_ka_kP_k(t)\le f(t)$ everywhere on $[-1,1)$.⁶

Maximise: $\tfrac12\bigl[N(N-1)a_0-N\sum_{k\ge1}a_k\bigr]$.

Output: a proof. Any feasible $h$ gives a valid lower bound for every configuration, and checking it doesn’t require trusting the software that found it.

3.5 When is the bound exact?

For the bound to equal the energy of a particular configuration, both inequalities we used must be equalities for it:

• Touching: $h(t)=f(t)$ at every angle $t$ that occurs in the configuration, so that $h$ meets $f$ from below at those points.
• Balance: $a_kS_k=0$ for every $k\ge1$, so $S_k=0$ at every level where the polynomial has $a_k\gt 0$.

How many coefficients does touching take? Meeting $f$ from below at an interior angle forces the two curves to have the same height and the same slope there, which is two conditions. At the endpoint $t=-1$ only the height has to match, which is one. The icosahedron has three angles besides $t=1$, namely $-1$ and $\pm1/\sqrt5$, so touching costs $1+2+2=5$ conditions. A degree-4 polynomial has exactly five coefficients, $a_0,\dots,a_4$, so the touching polynomial is fully determined. Solving that system gives non-negative $a_k$, and the resulting $h$ stays below $f$ everywhere. The icosahedron is also balanced at every level up to $k=5$, so nothing the polynomial uses costs anything.

In general, $m$ touching angles mean up to $2m$ conditions, hence a polynomial of degree up to $2m-1$, hence balance needed up to level $2m-1$. This is why Cohn and Kumar’s sharp configurations are the ones with $m$ distinct angles that are also $(2m-1)$-designs. Every one of them gets an exact bound. On the sphere, those are exactly $N=2,3,4,6,12$.

Both conditions hold for the icosahedron, so the bound equals its energy exactly. Call the difference between the true minimum energy and the best bound the gap; for the icosahedron it is zero.

The pentagonal bipyramid fails both conditions. Its angles are $-1$, $\cos144^\circ$, $0$ and $\cos72^\circ$. No polynomial that obeys the rules can touch $f$ at all four of them and also avoid paying for imbalance, so the best one compromises. Its bound is about 14.40 against the true energy of about 14.45, a shortfall of 0.36%. Of the 0.052 gap, 0.027 comes from $h$ sitting below $f$ at the real angles, and 0.025 comes from the balance payments $a_kS_k$, mostly at $k=2$.

3.6 The adversary program

The polynomial program never looks at the points themselves. It sees only the list of $N(N-1)$ ordered-pair angles and the facts $S_k\ge0$. That suggests a mirror-image question. Suppose an adversary wants to beat our bound by proposing a list of angles, with no points behind it, that merely passes every balance test $S_k\ge0$. How low can the energy of such a list go?

Describe a list by a distribution $\mu$ on $[-1,1)$ that puts one unit of mass at the angle of each ordered pair. The $N$ terms with $i=j$ contribute $N$ to every $S_k$, so the test $S_k\ge0$ reads $\int P_k\,d\mu\ge-N$.

Adversary program

Variables: a distribution $\mu\ge0$ of angles on $[-1,1)$, with total mass $N(N-1)$.

Constraints: $\int P_k\,d\mu\ge-N$ for every $k=1,\dots,d$, with the same degree $d$ as the polynomial program.

Minimise: the energy $\tfrac12\int f\,d\mu$.

Output: the lowest energy any angle list can have while passing every test. If no real configuration has that list, call it a fake.

These two programs are duals of each other. Every variable in one is a constraint in the other:

If you write the adversary program with a Lagrange multiplier for each of its constraints, those multipliers are exactly $a_0$ and $a_1,a_2,\dots$, and the polynomial program falls out.

Duality has two consequences. First, every polynomial gives a number below every adversary’s energy. That is the derivation from 3.4, applied to the adversary’s list instead of a real configuration’s. Second, the best polynomial and the best adversary meet.⁷ At the meeting point the two exactness conditions from 3.5 become optimality conditions: the adversary puts mass only at angles where $h$ touches $f$, and $a_k\gt 0$ only at levels where the adversary’s test is tight, which means $S_k=0$.

A real configuration is itself an adversary, because its own list of angles passes every test. So the bound is exact precisely when the real minimal configuration is also the best adversary. In that case there is no gap, and no fake beats the real thing.

That is the gap in concrete terms. For $N=7$, the 0.052 is the energy the adversary saves with a list of angles that no seven points can have. Every improvement from here on is a new test that fakes fail and real configurations pass.

4. Three-point bounds and cuts

Section 3 ended with an adversary who beats the polynomial bound by proposing lists of angles that no real configuration could have. The polynomial program only sees pairs, so the adversary never has to make its angles fit together into triangles. Three points carry more structure than their three pair angles taken separately: given two of the angles, the third can’t be arbitrary. The three-point bound of Bachoc and Vallentin, adapted to energy by Cohn and Woo, adds some of those tests. Section 3 built the two-point bound from a squared length. This section builds the three-point bound from a squared length too, but one that lives at each point of the configuration separately.

This is the most technical part of the post. If you only want the story, skip to 4.4: the one thing to take from 4.1 to 4.3 is that triangles give the bound new tests, which shrink the room fakes have.

4.1 Looking from one point

Fix one point $x$ of the configuration and put it at the north pole. Every other point $y$ has a height $u=\langle x,y\rangle$, the cosine of its angle from $x$, and a bearing $\theta_y$, its longitude as seen from $x$. The picture below lets you move two such points around.

Writing out coordinates, two other points $y,z$ with heights $u,v$ and bearings $\theta_y,\theta_z$ have these coordinates and this inner product:

Up to rotating or reflecting the sphere, a triple of points $x,y,z$ is determined by its three mutual inner products, $u=\langle x,y\rangle$, $v=\langle x,z\rangle$ and $t=\langle y,z\rangle$. We call $(u,v,t)$ the triple’s triangle type. It records the shape of the triangle the three points form, and the order in which we list the three numbers doesn’t matter. The three numbers can’t be chosen freely, though. Given $u$ and $v$, the formula above forces $t$ into the interval $uv\pm\sqrt{(1-u^2)(1-v^2)}$, and the bearing difference $\sigma=\theta_y-\theta_z$ decides where in that interval it lands. So the triangle types that can occur fill a three-dimensional region, not the whole cube $[-1,1]^3$. A list of pair angles never says which pairs share a point, so the adversary of 3.6 never had to fit its angles into real triangles. The three-point bound closes part of that loophole.

4.2 Fourier modes seen from one point

In section 3, positivity came from a sum of vectors: $S_k$ was the squared length of $\sum_i\varphi_k(e_i)$. Here we do the same thing around each base point $x$. Seen from $x$, every other point is just a height and a bearing, so the points form a pattern around a circle.

A pattern of bearings $\theta_y$ has Fourier coefficients. The $k$-th one, $\sum_ye^{ik\theta_y}$, is a complex number, which is a vector in the plane. For $k=1$ it adds up the points’ directions, and for $k=2$ it first doubles every bearing, and so on. Its squared length is non-negative, and expands as $\bigl|\sum_ye^{ik\theta_y}\bigr|^2=\sum_{y,z}\cos\bigl(k(\theta_y-\theta_z)\bigr)\ge0$. This is the circle version of the squared-length trick from 3.2.

So far every point’s bearing counts the same, whatever its height. But the three-point bound needs to see bearing and height together, because a triple’s inner products $(u,v,t)$ involve both. The squared-length argument gives us a way in, since it works for any weights. Give each point $y$ a weight $W(u_y)$ that depends only on its height $u_y$, and form the complex number $\sum_yW(u_y)\,e^{ik\theta_y}$. Its squared length is non-negative whatever $W$ is. Every choice of $W$ gives its own non-negative quantity, and we want a whole family of them.

We build the family from polynomials, for two reasons:

1. everything in the end has to be a polynomial in the inner products, because a semidefinite program can check polynomial inequalities.
2. polynomials of higher degree can approximate richer weightings, so the degree of the certificate controls how many weightings it can use.

Concretely, take the weights $W_a(u)=p_a(u)\,(1-u^2)^{k/2}$ for $a=0,1,2,\dots$. Here $p_a$ is the $a$-th Legendre polynomial (any polynomial basis would do), and the extra factor $(1-u^2)^{k/2}$ gives no weight to a point at either pole, $u=\pm1$, where its bearing means nothing. It also keeps the algebra polynomial, as we’ll see. For $k\ge1$, the weight with $a=0$ is positive away from the poles and largest at the equator. The weight with $a=1$ counts points in the northern half positively and points in the southern half negatively. The weight with $a=2$ gives points near the equator the opposite sign from points nearer the poles. For each level $k$ and weight number $a$, attach to $x$ the complex number

$w_a(x)=\sum_{y}p_a(u_y)\,(1-u_y^2)^{k/2}\,e^{ik\theta_y}.$

The sum runs over every point of the configuration, and $x$ itself drops out for $k\ge1$ because its height is 1. For $k=1$ and $a=0$, $w_0(x)$ is exactly the horizontal part of $\sum_yy$, the vector from 3.2.

For example, consider our solution for N = 7, the pentagonal bipyramid, and stand at one of its poles. The picture shows what the other points look like from there, and the arrows that add up to $w_0(x)$ at each level $k$.

Collect the numbers $w_0(x),w_1(x),\dots$ and form the matrix of their inner products. Here it is for three weights, with $\bar w$ the complex conjugate:

A matrix of inner products is always positive semidefinite, written $\succeq0$, which means $c^{\mathsf T}Mc\ge0$ for every real vector $c$. The reason is the same as for a squared length being non-negative. For our matrix, $c^{\mathsf T}M_k(x)\,c=\bigl|\sum_ac_aw_a(x)\bigr|^2\ge0$. Multiplying out, the bearings drop out: each entry of $M_k(x)$ is a sum over pairs of points $y,z$ of a polynomial in the three inner products $\langle x,y\rangle,\langle x,z\rangle,\langle y,z\rangle$. Collect these polynomials into a matrix $Y_k$, and add up over every choice of base point $x$:

The sum runs over all ordered triples, repeats included. A sum of positive semidefinite matrices is positive semidefinite, so $M_k\succeq0$ for every configuration and every level $k$.⁸ This is the three-point version of $S_k\ge0$.

In the polynomial program each $S_k$ came with a weight $a_k\ge0$; here each matrix $M_k$ gets a matrix weight $F_k\succeq0$, and $\operatorname{tr}(F_kM_k)$ plays the role of $a_kS_k$; setting every $F_k=0$ gives back the polynomial program.

4.3 The three-point program

The facts $M_k\succeq0$ do the same job as $S_k\ge0$ did in 3.4. For any configuration, $\operatorname{tr}(F_kM_k)\ge0$, because the trace of a product of two positive semidefinite matrices is non-negative. So these terms can be added to an inequality for free and then thrown away, as before.

What ties the certificate to the energy is a per-triple inequality. A triple of points has an energy share, $\tfrac13\bigl(f(u)+f(v)+f(t)\bigr)$, the average of its three pair energies, where $f(t)=1/\sqrt{2-2t}$ is the energy curve from Lemma 1. Each pair of points lies in $N-2$ triples, so adding the shares of all $\binom{N}{3}$ triples gives $\tfrac{N-2}{3}E$.

On the certificate’s side, each triple gets a matching term $R(u,v,t)$, assembled from pieces we have already met. The formula below is schematic: the repeated-point terms are real bookkeeping that we don’t spell out.

• The first line is the matrices $Y_k$ of 4.2, weighted by $F_k$, just as $a_k$ weighted $P_k$ in 3.4.
• The second line is the wiggles of the polynomial from 3.4, with $g(t)=\sum_{j\ge1}a_jP_j(t)$. The constant $a_0$ of 3.4 is absorbed into the bound $e$.
• The third line is the same bookkeeping as the $i=j$ terms in 3.4, for triples in which two points coincide.

Add $R$ up over the $\binom{N}{3}$ distinct triples of a configuration. The repeated-point terms fill in the triples with repeated points that $M_k$ and $S_j$ include and (T) does not, and the pieces become a positive multiple of $\sum_k\operatorname{tr}(F_kM_k)+\sum_{j\ge1}a_jS_j$, which can’t be negative. So for a number $e$ that we hope is a lower bound, the requirement is

Sum (T) over the $\binom{N}{3}$ triples of a configuration. The left side adds up to $\tfrac{N-2}{3}(E-e)$, the right side to something non-negative, and so $E\ge e$.

There is one inequality (T) for every possible triangle type, which is infinitely many, but each is linear in the unknowns. Because the unknowns include positive semidefinite matrices, this is a semidefinite program (SDP) and not a linear one. To solve it in practice, sample a grid of triangle types, solve the resulting problem, find the triple where (T) fails worst, add it to the list, and repeat. The final certificate is then checked over the whole region in exact arithmetic. So the three-point bound keeps every test from section 3 and adds the matrix tests on top, which means the best three-point bound is never weaker than the best polynomial bound.

As in 3.5, the bound is exact for a configuration when two things hold. Inequality (T) must be an equality at every triangle type the configuration has, and $\operatorname{tr}(F_kM_k)=0$ at every level the certificate uses. The triangle types of the minimal configuration play the role the angles played before. N = 7 and N = 8 have 5 each, and N = 9 and N = 10 have 12 each. The adversary of 3.6 now proposes a list of triangles instead of a list of angles. It has to pass the three-point tests as well, so it has fewer fakes available and the bound can only go up.

4.4 Cutting the space

The second tool is a case split. Suppose we only care about configurations in which every pair has its inner product in some range, such as “no pair closer than θ”. Then (T) only has to hold for triples whose three inner products all lie in that range. That is fewer constraints, so more certificates are allowed and the best bound can only go up. In the adversary’s picture the cut adds a constraint instead: its list must stay inside the range, which removes the fakes that lived outside it.

A cut also has to be paid for. A separate argument must show that any configuration breaking the cut has energy above the minimum. The proof then has two cases: either some pair breaks the cut, or none does and the restricted bound applies.

The two proofs look different, but they run the same algorithm.

Branch and bound over inner products

1. Branch. Cut configuration space by inequalities on inner products.
2. Bound. On each piece, solve the three-point SDP with the domain restricted to that piece.
3. Close. Every piece either has a bound above the minimum energy, or is a small neighbourhood of the minimal configuration, where a local argument applies.

For N = 7, Hung chose cuts on the most-antipodal pair, because the minimal configuration has an antipodal pair, its poles. For N = 8, Tooby-Smith and Zughaid chose a cut on the closest pair. Nothing in the method fixes the choice.

The local step at the end depends on how flat the energy is near the minimal configuration. At N = 7 the smallest nonzero eigenvalue of the Hessian (modulo rotations) is 0.0092, which makes the last step delicate. At N = 9 and 10 it is 0.144 and 0.096, 10 to 16 times larger.⁹

5. What the numbers say

Which N are close to provable? A natural guess is that the fewer distinct inner products a minimal configuration has, the fewer touching points $h$ needs, and so the smaller the gap. Comparing the two bounds across many N tests it.

The guess fails for the LP bound and mostly holds for the three-point bound.

• In the N we tested, LP gaps track N, not symmetry. N = 5 has only 3 distinct inner products but the worst LP gap of all. Past the sharp cases, the relative LP gap shrinks fairly steadily with N, from about 10⁻³ near N = 10 to 2×10⁻⁴ at N = 32. N = 32’s small LP gap reflects its size, not its symmetry.
• Three-point gaps track symmetry. The three-point bound is numerically tight at N = 5, in line with de Laat’s result for a related bound. At the symmetric N = 9 and 10 it shrinks the LP gap by a factor of 13 to 25. At the low-symmetry N = 11 and 13 it only roughly halves it, leaving relative gaps 10 to 25 times larger than at 9 and 10. N = 14 sits in between.

Working hypothesis

The LP bound only sees how many pairs sit at each angle, and its error depends mostly on N. The three-point bound also sees how triangles fit together, and that is where symmetry pays off: a minimal configuration with few distinct inner products and triangle types leaves fewer ways for a fake to undercut it. Cuts then remove the fakes that remain.

6. From one N to all N

A general proof for every N is out of reach. Nobody even knows the minimal configurations for large N except numerically, and even finding near-optimal configurations efficiently, for the closely related log energy, is Smale’s 7th problem. What does look reachable is making each N routine.

• The bounds can always be made exact in principle. De Laat (2020) built a hierarchy of $k$-point bounds and proved it becomes exact by level $N$, though turning that into a finite, exactly checkable certificate is a separate problem. The practical question is whether low levels plus a few cuts suffice, and the data so far say yes for small symmetric N.
• Automate the branching. Read off which fakes the SDP can’t exclude, choose the cut that removes them, and repeat, generating exact certificates for each piece. Then a new N costs compute, not an 18,000-line hand-built proof.
• Reuse structure across N. N = 10, 12 and 14 are the same shape with ring size 4, 5 and 6. A proof with the ring size as a parameter could cover the family.

6.1 Next targets

6.2 Challenge: N = 32

Its conjectured minimal configuration, an icosahedron nested with its dual dodecahedron, is about as symmetric as 32 points can be: the full icosahedral group of 120 symmetries, with no free parameters.

The three-point bound barely improves on the polynomial bound there, and no cut we know of changes that. So the method above stops short, and something new is needed: a stronger bound, a better case split, or an argument we haven’t thought of.

6.3 Prove it to me

Our expectation is that N = 10 falls next, then N = 9, each with the machinery above: a three-point certificate, a few cuts, and an exact check of every piece. Apart from the theorems we cite, the bounds and gaps in this post come from floating-point solvers, and the gap between “our solver says it is nearly exact” and “Lean accepts it” is the whole difficulty.

Finally, our challenge to you, dear reader: produce an exact, machine-checkable certificate for N = 32 and we’ll write it up and credit you. We’re just as excited for N = 10 or N = 9, if you get there first - the exact Lean statements for all three are on Prove2Me, where any proof can be submitted and checked.

A failure is worth reporting too. If the method above breaks on N = 9 or 10, tell us where. The sharpest version of that failure, a fake that no cut removes, would be a real result in its own right.

Bibliography & Footnotes

1. N. N. Andreev, An extremal property of the icosahedron, East J. Approx. 2 (1996), 459–462.
2. C. Bachoc, F. Vallentin, New upper bounds for kissing numbers from semidefinite programming, J. Amer. Math. Soc. 21 (2008). arXiv:math/0608426
3. H. Cohn, A. Kumar, Universally optimal distribution of points on spheres, J. Amer. Math. Soc. 20 (2007). arXiv:math/0607446
4. H. Cohn, J. Woo, Three-point bounds for energy minimization, J. Amer. Math. Soc. 25 (2012). arXiv:1103.0485
5. P. Delsarte, J. M. Goethals, J. J. Seidel, Spherical codes and designs, Geom. Dedicata 6 (1977), 363–388. doi:10.1007/BF03187604
6. R. E. Schwartz, The five-electron case of Thomson’s problem, Exp. Math. 22 (2013). arXiv:1001.3702
7. D. de Laat, Moment methods in energy minimization, Trans. Amer. Math. Soc. (2020). arXiv:1610.04905
8. Kryvonos, Liehr, Taylor, Energy minimization for eight points on the sphere (2026). arXiv:2609.22077
9. J. Tooby-Smith, A. Zughaid, Thomson-N-8-Warrant (2026). github.com/jstoobysmith/Thomson-N-8-Warrant
10. V. A. Yudin, Minimum potential energy of a point system of charges, Diskret. Mat. 4 (1992), 115–121 (in Russian); English translation in Discrete Math. Appl. 3 (1993), 75–81.
11. Thomson N = 7 Lean proof package, snapshot 2026-09-27. github.com/huwngtran/thomson-n7-lean