The 5-Electron case of Thomson’s Problem

∗
Richard Evan Schwartz
February 9, 2010

Abstract
We give a rigorous, computer-assisted proof that the triangular
bi-pyramid is the unique configuration of 5 points on the sphere that
globally minimizes the Coulomb (1/r) potential. We also prove the
same result for the (1/r 2 ) potential. The main mathematical contri-
bution of the paper is a fairly efficient energy estimate that works for
any number of points and any power law potential.

1 Introduction
1.1 Background
The problem of finding how electrons optimally distribute themselves on the
sphere is a well-known and difficult one. It is known as Thomson’s problem,
and dates from J. J. Thomson’s 1904 publication [T]. Thomson’s problem is
of interest not just to mathematicians, but also to physicists and chemists.
Here is a mathematical formulation. Let S 2 ⊂ R3 be the unit sphere.
Let P be a collection of n distinct point p1 , ..., pn ∈ S 2 . Let E : R+ → R+
be some function. the total energy to be the sum
X
E(P ) = E(kpi − pj k). (1)
i>j

Here k · k is the usual norm on R3 , and kpi − pj k is the distance from pi to pj .

The question is then: What does P look like if E(P ) is as small as possible?
∗
Supported by N.S.F. Research Grant DMS-0072607

1
 The question is perhaps too broad as stated, because the answer likely
depends on the function E. To consider a narrower question, one restricts
the class of functions in some way. For instance, a natural class of potential
functions is given by the power laws
1
E(r) = ; e ∈ (0, ∞). (2)
re
The case e = 1 is specially interesting to physicists. It is known as the
Coulomb potential . When energy is measured with respect to the Coulomb
potential, the points are naturally considered to be electrons.
There is a large literature on Thomson’s problem. One early work on
Thomson’s problem is [C]. The paper [SK] gives a nice survey in the two
dimensional case, with an emphasis on the case when n is large. The paper
[BBCGKS] gives a survey of results, both theoretical and experimental,
about highly symmetric configurations in higher dimensions. The fairly re-
cent paper [RZS] has some theoretical bounds for the logarithmic potential,
and also has a large amount of experimental information about configura-
tions minimizing the power laws on the 2-sphere. The website [CCD] has
a list of experimentally determined (candidate) minimizers for the Coulomb
potential for n = 2, ..., 972.
There are certain values of n where the minimal configuration is rigorously
known for all the power laws.
• When n = 2, the points of P are antipodal.

• When n = 3, the points of P make an equilateral triangle in an equator.

• When n = 4, the points of P are the vertices of a regular tetrahedron.

• When n = 6, the points of P are the vertices of a regular octahedron.

• When n = 12, the points of P are the vertices of a regular icosahedron.

The cases n = 2 is trivial to prove, and the case n = 3 is an easy exercise.
The cases n = 4, 6, 12 are all covered in [CK, Theorem 1.2], a much broader
result concerning a fairly general kind of energy potential and a class of point
configurations (not restricted to 2 dimensions) called sharp configurations.
For the specific cases cases n = 4, 6 see the older work [Y].
The case n = 5 is conspicuously absent from the above list of known
results. Everyone agrees, based on numerical simulation, that the triangular

2
bi-pyramid (TBP) seems to be the global minimizer for the Coulomb poten-
tial. In the TBP, two points are antipodal points on S 2 and the remaining
3 points form an equilateral triangle on the equator midway between the
two antipodal points. More generally, numerical experiments 1 suggest the
following.
• The TBP is a local minimizer for the power law potential with exponent
e if and only if e ∈ (0, e1 ), with e1 ≈ 21.147123.
• The TBP is a global minimizer for the power law potential with expo-
nent e if and only if e ∈ (0, e2 ), with e2 ≈ 15.040808.
For e > e2 , it seems that the global minimizer is a pryamid with square base.
The precise configuration depends on e in this range.
In spite of detailed experimental knowledge about the case n = 5, it
seems there has not ever been a proof that the TBP is global minimum for
any power law potential. In particular, this has not been proved for the
Coulomb potential. As far as we know, there are two rigorous results for the
case n = 5.
• The paper [DLT] contains a (traditional) proof that the TBP maxi-
mizes the geometric mean of the pairwise distances between the points.
This case corresponds to the logarithmic potential E(r) = − log(r).
• The paper [HS] contains a computer-aided proof that the TBP maxi-
mizes the potential for the exponent e = −1. That is, 5 points on the
sphere arrange themselves into a TBP so as to maximize the total sum
of the pairwise distances.

1.2 Results
It is the purpose of this paper to prove the following result.
Theorem 1.1 The TBP is the unique configuration of minimal energy with
respect to the Coulomb potential.
Our proof is computer aided, and similar in spirit to [HS]. While the
argument in [HS] is exactly tailored to understanding the sums of the dis-
tances, our method is rather insensitive to the precise power law being used.
Just to illustrate this fact, we prove
1
Lacking a handy reference for these experiments, we performed our own.

3
Theorem 1.2 The TBP is the unique configuration of minimal energy with
respect to the 1/r 2 potential.

We could certainly add other exponents to our list of results. However,

we currently have implemented the interval arithmetic in such a way that
it works just for the exponents e = 1, 2. We did this to avoid using the
pow function pow(a, b) = ab , which is not covered by the same guarantees in
the IEEE standards as are the basic arithmetic operations. See [I] and [I2].
Mainly what stops us is a sense of diminishing returns.
However, with a view towards an eventually broader application, we try
to state as many estimates as we can for the general power law. In every
situation, we try to point out the exact generality with which the construction
holds. For example, our main technical result, Theorem 5.1, holds for any
function E satisfies the conditions discussed in §5.2. Also, the version of
Theorem 5.1 we prove works for general n-point configurations.

1.3 Outline of the Proof

It has probably been clear for a long time that one could prove a result like
Theorem 1.1 using a computer program. The main difficulties are technical
rather than conceptual. The hard part is getting estimates that are sharp
enough so as to lead to a feasible calculation. Our paper doesn’t really have
any dramatic new ideas. It just organizes things well enough to get the job
done.
We begin by eliminating some obviously bad configurations from consider-
ation. For example we show that no two points in a minimizing configuration
lie within 1/2 units of each other. This result holds for any power law po-
tential. See Lemma 2.1. After eliminating these bad configurations, we are
left with a compact configuration space Ω of possible minimizers.
We use stereographic projection to transfer the configurations on S 2 to
configurations in C ∪ ∞. By keeping one point at ∞ and another on the
positive real axis, we end up with a natural description of Ω as the set
[0, 4] × [−2, 2]6 ⊂ R7 . Working in Ω, we define a natural way to subdivide a
rectangular solid subset of Ω into smaller rectangular solids.
We use a divide-and-conquer algorithm to show that any minimizing con-
figuration in Ω must lie very close to the TBP. Assuming that Q is a rect-
angular solid that does not lie too close to the point(s) of Ω representing the
TBP, we try to eliminate Q with a 3-step procedure.

4
 • We eliminate Q if we can see, based on the fact that the regular tetra-
hedron minimizes energy for 4 points, that no configuration in Q can
minimize energy. This is described in §3.4. This method of eliminating
Q cuts away a lot of the junk, so to speak, and focuses our attention
on the configurations where are fairly near the TBP.

• We eliminate Q if the configurations in it are redundant, in the sense

that some permutation of the vertices or obvious application of symme-
try changes these configurations to ones in a form that we deem more
standard. See §4.

• If the first two methods fail, we evaluate the energy E on the vertices
of Q and then apply an a priori estimate on how far E differs from a
linear function on Q. Our main result along these lines is Theorem 5.1.
This is our most powerful and general method of elimination, but we
only use it when all else fails.

If it is not possible to eliminate Q by any of our methods, we subdivide Q

into smaller rectangular solids and try again on each piece. This algorithm
is discussed, in outline, in §2.5. In the chapters following §2.5, we fill in the
details. Theorem 5.1 is the most subtle and important part of the paper. We
prove this result in §5-8.
The end result of our finite calculation is that the true minimizer of Q
must lie very close to the TBP. See Lemma 2.3. To finish the proof, we
use calculus to show that the TPB can be the only local minimum in the
region where we have confined the minimizer. We do this by showing that
the Hessian of E, the matrix of second partials, is positive definite throughout
the region of interest to us. See Lemma 2.4.

1.4 Computational Issues

We implement our computer program in Java, using interval arithmetic to
control the round-off errors. We will discuss this in §10. Our code only uses
the basic operations plus, minus, times, divide, and sqrt, and these operations
are performed in such a way as to avoid and arithmetic errors (such as taking
the square root of a negative number.) It takes about 6 hours for our program
to eliminate all configurations that are not within 2−14 , in the L∞ norm, of
a suitably normalized version of the TBP. See Lemma 2.3. Without the
interval arithmetic, the code runs in about an hour.

5
 Theorem 5.1 is the rate-determining step in our calculations. Numerical
experimentation suggests that our bound in Theorem 5.1 is off by a factor
somewhere between 2 and 4. In light of this fact, our calculation probably
ought to be about 10 times faster than it is. We can certainly improve
Theorem 5.1, but we haven’t been been able to think of a simple or dramatic
improvement.
A nice feature of our program is that we have embedded it in a graphical
user interface. The reader can watch the program in action and see how it
samples the configuration space. the reader can also manually construct a
rectangular solid subset of the configuration and then see a printout of all
the computational tests that are applied to it. This graphical aspect doesn’t
add anything to the formal proof, but it makes it less likely we have made a
gross computing error. We have tried to isolate the relatively small amount
of computer code that goes into the actual proof, so that it can be more
easily inspected.
The entire Java program is available from my website. See
[Link] The code
is fairly well documented, and we’re still working to improve the documen-
tation. Aside from graphical support files, all the files involved in the proof
have the Interval prefix. The directory contains a number of other files,
which support other features of the program. Even though the proof portion
of the code is done and working, the whole program is still somewhat in flux.
I plan to gradually improve the code as time passes, and update it as I go.

1.5 Acknowledgements
I would like to thank Henry Cohn for helpful conversations about placing elec-
trons on a sphere. Henry’s great colloquium talk at Brown university this fall
inspired me to work on this problem. I would also like to thank Jeff Hoffstein
and Jill Pipher for their interest and encouragement while I worked on this
problem. I would like to thank John Hughes for a very interesting discussion
about interval arithmetic. Finally, I would like to say that I learned how to do
interval arithmetic in Java by reading the source files from Tim Hickey’s im-
plementation, [Link]

6
2 Proof in Broad Strokes
2.1 Stereographic Projection
We find it convenient to work mainly with C ∪ ∞ rather than on S 2 . Our
reason for this is that a configuration space based on points in C has a
natural flat structure, and lends itself well to a nice subdivision scheme. The
subdivision scheme, which essentially amounts to cutting rectangular solids
into smaller rectangular solids, feeds into our divide-and-conquer algorithm.
All this is discussed in §2.5 below.
We map S 2 to C ∪ ∞ using stereographic projection:
x y
Σ(x, y, z) = +i . (3)
1−z 1−z
Σ is a conformal diffeomorphism which maps circles on S 2 to generalized
circles in C ∪ ∞. A generalized circle is either a circle or a straight line. We
have
Σ(0, 0, 1) = ∞. (4)
Thus, Σ maps a circle C ⊂ S 2 to a straight line if and only if (0, 0, 1) ∈ C.
The inverse map is given by
 2x 2y 2 
−1
Σ (x + iy) = , ,1 − . (5)
1 + x2 + y 2 1 + x2 + y 2 1 + x2 + y 2

Let k · k denote the usual norm on R3 . Here are two pieces of metric
information we will use later on:
2
kΣ−1 (z) − Σ−1 (∞)k = p (6)
1 + |z|2

dΣ−1 dΣ−1 2
= = . (7)
dx dy 1 + x2 + y 2
Equations 6 and 7 both have straightforward derivations, which we omit.
Slightly abusing notation, we define

E(z1 , z2 ) = E(Σ−1 (z1 ), Σ−1 (z2 )) (8)

Here E can be any energy potential. for points z1 , z2 ∈ C ∪ ∞. In this way,

we can talk about the energy of a configuration of points in C ∪ ∞.

7
2.2 The Triangular Bi-Pyramid
As in the introduction k · k denotes the usual norm on R3 . Unless stated
otherwise, all the distances we measure in R3 are taken with respect to the
Euclidean metric. What we say here works for any energy potential.
In the TBP, we have the following information.

• One pair of points is 2 units apart.

√
• 3 pairs of points are 3 units apart.
√
• 6 pairs of points are 2 units apart.

Accordingly, the energy of the TBP, with respect to E, is

√ √
ME = E(1/2) + 3E( 3) + 6E( 2). (9)

It is well known that the

p regular tetrahedron minimizes the energy for
4 points. All points are 8/3 units apart. For this reason, the regular
tetrahedron has energy p
TE = 6E( 8/3). (10)
We can use this fact to get some crude bounds on 5-point configurations.
Given a 5 point configuration P and a point p ∈ P , we define
X
E(P, p) = E(kp − qk). (11)
q∈P {p}

We have the immediate estimate

E(P ) ≥ E(P, p) + TE . (12)

Thus, if P is an energy minimizer we must have

E(P, p) ≤ ME − TE ; ∀p ∈ P. (13)

This is one of the criteria we will use to eliminate certain configurations from
consideration.
We can also use Equation 13 to give information about pairs of points
within a minimizing configuration. We will consider this in the next section.

8
2.3 Estimates for the Power Laws
The fairly weak results in this section are designed to give us some control on
the size of the configuration space we must consider. For any given exponent,
these results are just short calculations. It is only our desire to handle all
exponents at the same time that adds complexity to the proof.
Lemma 2.1 Let p, q1 , q2 be 3 points of an energy minimizer with respect to
a power law potential. Then kp − q1 k > 1/2.

Proof: Let E(r) = 1/r e . Since E is decreasing, E(kp − rk) ≥ E(2) for all
p, r ∈ S 2 , because S 2 has chordal diameter 2. In light of Equation 12, this
lemma is true provided that
TE + 3E(2) + E(1/2) − ME > 0 (14)
When E is as above, our problem boils down to showing that
φ(e) = 21−e − 31−e/2 + 2e − 31−e/2 + 21−(3e)/2 31+e/2 > 0. (15)
This is an exercise in calculus. We compute
φ(0) = 0; φ′ (0) > 0; φ′′ (0) > 1. (16)
We compute φ′′′ (e) = A(e) − B(e), where
log(3) 1−e/2
A(e) = 3 log(2)2 2−2−e/2 + 2e log(2)3 + 3 > 0;
8
B(e) = 21−e log(2)3 + 2−2−(3e)/2 31+e/2 log(8/3)3 < 2. (17)
In short, φ′′′ (e) > −2. Taylor’s theorem with remainder now tells us that
e2 e3
φ(e) > − . (18)
2 3
Hence φ(e) > 0 for e ∈ (0, 3/2). A similar computation for φ′ shows that
φ′ (e) > −10. (19)
We compute that
φ(3/2 + j/10) > 1; j = 0, ..., 85. (20)
Combining the last two equations, we see that φ(e) > 0 for all e ∈ [3/2, 10].
Finally, for e > 10, the result is obvious. ♠

9
Lemma 2.2 Let p, q1 , q2 be 3 points of an energy minimizer with respect
to a power law potential. Assume that E(kp − q2 k) ≤ E(kp − q1 k). Then
kp − q1 k > 1/2 and kp − q2 k > 1.

Proof: We have the inequality

E(kp − q1 k) + E(kp − q2 k) ≤ ME − TE − 2E(2). (21)

Hence
ME − TE − 2E(2)
E(kp − q2 k) ≤ . (22)
2
Establishing this inequality boils down to showing that

φ(e) = 2 + 21−e − (3)(21−e/2 ) − 2−e − 31−e/2 + 21−(3e)/2 31+e/2 > 0. (23)

This time we compute

φ(0) = 0; φ′ (0) > 0; φ′′ (0) > 1/4. (24)

Examining the terms of φ′′′ , as in the previous lemma, we find that

φ′′′ (e) > −12. (25)

Taylor’s theorem now tells us that φ > 0 on [0, 1/16). A similar computation
shows
φ′ (e) > −10. (26)
We now compute that

φ(1/16 + j/2000) > 1/200; j = 0, ..., 1875. (27)

Combining the last two equations, we see that φ(e) > 0 for e ∈ [1/16, 1].
Next, we compute that

φ(1 + j/100) > 1/10; j = 0, ..., 900. (28)

This shows that φ(e) > 0 for e ∈ [1, 10]. For e > 10 the result is again
obvious.

10
2.4 Planar Configurations
We took the trouble to prove Lemmas 2.1 and 2.2 for any power law so that
what we say in this section works for any power law.

Basic Definition: Let P = {p0 , ..., p4 } be a configuration of 5 points on

S 2 . Z = Σ(P ) = {z0 , ..., z4 } be the corresponding configuration in C ∪ ∞.
We rotate S 2 so that

z4 = ∞; z0 ∈ R+ E(z0 , z4 ) = max E(zi , zj ). (29)

i<j

The TBP: Now we discuss what the TBP looks like with this normalization.
The TBP has two kinds of points. We say that the polar points are the two
antipodal points in the configuration. We say that the three remaining points
are equitorial . When z4 corresponds to a polar point, we have the following
configuration, which is unique up to the permutation of z1 , z2 , z3 :

z0 = 1; z1 = exp(−2πi/3); z2 = 0; z3 = exp(2πi/3);
(30)
When z4 corresponds to an equitorial point, we have
√ √
z0 = 1; z1 = −i 3/2; z2 = −1 z3 = i 3/2; (31)

Later on, we will use a symmetry argument to avoid having to deal with the
second of these configurations.

The Space of Minimizers: We check 3 facts using Equation 6.

√
1. If |z| = 1 then kΣ−1 (z) − Σ−1 (∞)k = 2.

2. If |z| ≥ 2 then kΣ−1 (z) − Σ−1 (∞)k < 1.

3. If |x| ≥ 4 then kΣ−1 (z) − Σ−1 (∞)k < 1/2.

Suppose now that Z is a minimizer for some power law potential. It

follows from Lemma 2.1 and Item 3 that |zj | < 3 for j = 1, 2, 3. If follows
from Item 1 and Lemma 4.2 that |z0 | ≥ 1. We conclude that z0 ∈ [1, 4]. It
follows Lemma 2.2 and Item 2 that |zj | ≤ 2 for j = 1, 2, 3. We conclude that
all the configurations we need to consider lie in the compact region compact
region Ω = [0, 4] × [−2, 2]6 ⊂ R7 . The map is given by just stringing out the

11
coordinates of the points z0 , z1 , z2 , z3 .

Dyadic Objects: Given squares Q1 , Q2 ⊂ C, we write Q1 → Q2 if Q2

is one of the 4 squares obtained by dividing Q1 in half along both directions.
We say that a square Q is a dyadic square if there is a finite

[−2, 2]2 → Q1 → . . . → Qn = Q. (32)

The sides of Q are necessarily parallel to the coordinate axes, and the vertices
have dyadic rational coordinates. We call the dyadic square normal if it does
not cross the coordinate axes. A single subdivision of [−2, 2] produces normal
dyadic squares.
We make the same definition for line segments as for squares., except that
the notation S1 → S2 means that S2 is one of the two segments obtained by
cutting S1 in half. We say that a dyadic segment is a line segment S such
that there is a finite chain

[0, 4] → S1 → . . . → Sn = S. (33)

We say that a dyadic box in Ω is a set of the form Q0 × Q1 × Q2 × Q3 , where

Q0 is a dyadic segment and Qj is a dyadic square for j = 1, 2, 3. The whole
space Ω itself is a dyadic box.

Subdivision: Let Bj = Qj0 × Qj1 × Qj2 × Qj3 be a dyadic box for j = 1, 2.

We write B1 →k B2 if

• Q1i = Q2i for i 6= k.

• Q1k → Q2k .

The kth subdivision of B1 is the union of all B2 such that B1 →k B2 . When

k = 0, this union consists of two dyadic boxes. When k = 1, 2, 3, the union
consists of 4 dyadic boxes. The set of all dyadic boxes in Ω forms a directed
tree. Each dyadic box points to 14 = 2 + 4 + 4 + 4 smaller dyadic boxes.
The divide-and-conquer algorithm will be a depth-first search through
this tree. It seems that the speed of the program depends a lot on how we
do the subdivisions, so we will explain this in detail in the next section.

12
2.5 The Divide and Conquer Algorithm
Our discussion applies to a general power law potential, but we only apply
the algorithm to the Coulomb potential.
Let ǫ > 0 be some small number. In this section, we explain in the ab-
stract how we show, with a finite calculation, that any winning configuration
lies within ǫ of one of the two TBP configurations described in §2.4. There
are 5 components to our program:

Confinement: We say that a dyadic box

Q = Q0 × Q1 × Q2 × Q3 . (34)

is ǫ-confined if, for each j = 0, 1, 2, 3, the dyadic square Bj is contained in

the open square of side length ǫ centered on the point zj from Equation 30.

Tetrahedral Eliminator: We eliminate Q is if satisfies the hypotheses

of Lemma 3.6. In this situation, every configuration in Q violates Equation
13. This function speeds up our calculations quite a bit.

Redundancy Eliminator: We will isolate a subset Ω′ ⊂ Ω with the fol-

lowing property: If Z ∈ Ω − Ω′ , then there is some Z ′ ⊂ Ω′ (obtained by
permuting the points of applying some obviously energy-decreasing move)
such that E(Z ′ ) ≤ E(Z). In §4.4 we will describe some simple tests we use
to show that Q ∈ Ω − Ω′ . We eliminate Q if it passes one of these tests.
One convenient property of Ω′ is that it contains the configuration in
Equation 30 but not the configuration in Equation 31. This means that we
can automatically eliminate configurations near the one in Equation 31. This
is why our notion of ǫ-confinement only mentions the configuration in Equa-
tion 30.

Energy Estimator: Now we come to the main point. Say that an en-
ergy estimator is a function Φ : S → R such with the following property.
For every configuration Z = {z0 , ..., z4 } with zj ∈ Qj , we have

E(Z) ≥ Φ(Q).

See Theorem 5.1 for the definition of our Energy Estimator. Theorem 5.1 is
our main technical result.

13
Depth First Search: Recall that Me is the energy of the TBP. Our pro-
gram maintains a list L of dyadic boxes. Initially, L just has the single box
Ω, the whole space. At a given stage of the program, the algorithm examines
the last box Q and eliminates it if one of three things happens:
1. Q is ǫ-confined.
2. Lemma 3.6 eliminates Q.
3. Onf the the tests in §4.4 eliminates Q.
4. Ψ(Q) > Me ,.
Otherwise, Q is eliminated and to L we append the dyadic boxes in the kth
subdivision of Q for some k ∈ {0, 1, 2, 3}. We will explain how k is deter-
mined momentarily. The algorithm halts if L is the empty list. In this case,
we have shown that, up to symmetries, any minimizer lies within ǫ of the
TBP.

Subdivision Rule: The error term in Theorem 5.1 has the form
3 X
X
ǫ(Qi , Qj ). (35)
i=0 j6=i

Here ǫ(Qi , Qj ) is a quantity that depends on the geometry of Qi and Qj .

Roughly, it varies quadratically with the side length of Qi . We write
X
ǫ(i) = ǫ(i, j). (36)
j6=i

Then, again, ǫ(i) depends roughly quadratically on the side length of Qi . We

find the index k that maximizes the function i → Q(i) and then we use the
kth subdivision rule. Thus, we subdivide in such a way as to try to make the
error estimate in Theorem 5.1 as small as possible. When we compare this
method with a more straightforward method of subdividing so as to keep the
squares all about the same size, our method leads to a vastly faster compu-
tation.

Remark: The main difficulty in our proof is choosing an energy estimator

that leads to a feasible calculation. We found the subtle Energy Estimator
from Theorem 5.1 after quite a bit of trial and error.

14
2.6 The Main Results
Let E denote the energy function on Ω. For any s > 0 let Ωs denote those
configurations {zk } so that zk is contained in a square of side length s centered
at the kth point of the TBP, normalized as in Equation 30. We shall be
interested in the cases when

s = 2−11 . (37)
Running our computer program, we prove the following result.

Lemma 2.3 (Confinement) Any Coulomb energy minimizer in Ω has the

same energy as some configuration in Ωs/4 .

Let He denote the Hessian of E, with respect to the potential E(r) = r −e .

In §9 we will prove the following result.

Lemma 2.4 For e = 1, 2 and for any W ∈ Ωs , the matrix He (W ) is positive

definite.

Now let Z1 be a 5-point minimizer for the Coulomb potential. By the

Confinement Lemma, we can find a new configuration Z2 ∈ Ωs with the
same energy. But E has positive definite Hessian throughout Ωs , and (by
symmetry) the gradiant ∇E vanishes at Z0 , the TBP. Restricting E to a
straight line segment connecting Z0 to Z2 we see that E is a convex function
with a local minimum at Z0 . Hence E(Z2 ) > E(Z0 ). This proves Theorem 1.1.

Remarks:
(i) Our program also establishes the Confinement Lemma for the function
E(r) = r −2 . The same argument as above now establishes Theorem 1.2.
(ii) We didn’t need to compute all the way down to Ωs/4 . We could have
stopped at Ωs in the Confinement Lemma. However, an earlier version of
this paper had a weaker result on the Hessian, and we did the extra comput-
ing to accomodate this. There doesn’t seem to be any reason to throw out
our stronger computaional result since we (or, rather, the computer) took
the trouble to get it.

15
3 Separation Estimates
3.1 Overview
Let Σ be stereographic projection. The sets of interest to us have the form

Q∗ = Σ−1 (Q). (38)

where Q is either a dyadic segment or a dyadic square or the point ∞. We

will call Q a dyadic planar set and Q∗ a dyadic spherical patch. A dyadic
spherical patch is either the point (0, 0, 1), or an arc of a great circle, or else
a subset of the sphere bounded by 4 arcs of circles.
Let (Q1 , Q2 ) where Qj is a dyadic planar set. In this chapter we are
interested in the following two quantities.

ψmax (Q1 , Q2 ) = max kp1 − p2 k; pj ∈ Q∗j . (39)

ψmin (Q1 , Q2 ) = min kp1 − p2 k; pj ∈ Hull(Q∗j ). (40)

Here Hull(Q∗j ) is the convex hull of Q∗j . For technical reasons the lower bound
we seek need to work for a slightly larger range of points.
In this chapter, we will estimate ψ(Q1 , Q2 ) in terms of quantities that can
be determined by a finite computation. At the end of the chapter, we will
give, as an application of these estimates, the definition of the Tetrahedral
Eliminator discussed in §2.5.
In the next definitions, we set z = x + iy. If Q is a dyadic square or
dyadic segment, we define

x = min |x|; x = max |x|; y = min |y|; y = max |y|;

z∈Q z∈Q z∈Q z∈Q

2(x − x) p
δ= ; τ= dim(Q). (41)
1 + x2 + y 2
These quantities depend on Q but we usually suppress them from our nota-
tion. Here x − x is just the sidelength of Q. The quantity δ is an estimate√
on
∗
the side length of Q . Note that τ = 1 if Q is a dyadic segment and τ = 2
if Q is a dyadic square.
When Q = {∞} we set δ = τ = 0 and we don’t define the other quantities
know Q∗ exactly in this case.

16
 Given a pair (Q1 , Q2 ) of dyadic objects, we

δ1 τ1 + δ2 τ2
D = kΣ−1 (z1 ) − Σ−1 (z2 )k; δ= . (42)
4
Here zj is the center of Qj and δj = δ(Qj ), as defined above. Here is our
main result.

Lemma 3.1 (Bound) The following is true relative to any pair of dyadic
patches.
√ √
ψmin ≥ D(1 − δ 2 /2) − ( 4 − D 2 )δ; ψmax ≤ D + ( 4 − D 2 )δ.

Remark: Our proof will yield the better bounds

√ √
ψmin ≥ D cos(δ) − ( 4 − D 2 ) sin(δ) ψmax ≤ D cos(δ) + ( 4 − D 2 ) sin(δ),

but for computational reasons we want to avoid the trig functions. We have
used rational replacements which are quite close in practice to the trig func-
tions.

There is one situation where we can get a completely sharp upper bound.
We define
′
ψmax (Q1 , Q2 ) = max kp1 − p2 k; pj ∈ Q∗j . (43)
This time we mazimize over pairs (p1 , p2 ), where pj is a vertex of Q∗j . A
′ ′
finite computation gives ψmax . We define ψmin similarly. Recall that a dyadic
square is normal if it doesn’t cross the coordinate axes.

Lemma 3.2 (Perfect Bound) If Q1 is normal and Q2 = {∞} then we

′ ′
have ψmax = ψmax and ψmin = ψmin .

We define Ψmin to be the best bound we can get from the two lemmas
above (or 0, if no lemma applies.) We define Ψmax to be the best bound we
can get from the two lemmas above (or 2, if no lemma applies.)

The rest of the chapter is devoted to proving Lemma 3.1 and Lemma 3.2.

17
3.2 Proof of Lemma 3.1
Lemma 3.3 Let A and B be arcs of the unit circle. Let ∆A and ∆B denote
the arc lengths of A and B Let DA and DB denote the distance between the
endpoinds of A and B respectively. Let
∆A − ∆B
δ= .
2
Then q 
2
DB = DA cos(δ) ± 4 − DA sin(δ).

Proof: We have the relations

DA = 2 sin(∆A /2); DB = 2 sin(∆B /2) = 2 sin(∆A /2 − δ). (44)

Hence, by the angle addition formula,

DB = DA cos(δ) − 2 cos(sin−1 (DA /2)) sin(δ). (45)

Noting that √
cos(sin−1 (x)) = 1 − x2 , (46)
we see that Equation 45 is the same as the result we want. ♠

Lemma 3.4 Let Q be a dyadic set. Every point of Q∗ is within δτ /2 units

of Σ−1 (zQ ).

Proof: When evaluated at all points of Q, the quantity in Equation 7 is

maximized at the point (x, y). Its value at that point is exacrly

2 δ
= .
1 + x2 + y 2 s

Therefore Σ−1 expands distances on Q by at most a factor of δ/s, and every

point of Q is within sτ /2 of z. From here the result is obvious. ♠

Now we give the main argument for Lemma 3.1. We consider the lower
bound first. Given a point p ∈ S 2 and some r > 0 we let H(p, ǫ) ⊂ S 2 denote

18
the convex hull of the set of points on S 2 that are within ǫ of p in terms of
arc length on S 2 . We call such a set a cap.
It follows from Lemma 3.4 and the convexity of caps that
Hull(Q∗j ) ⊂ Hj = H(pj , δj τj /2). (47)
Here pj = Σ−1 (zj ), where zj is the center of Qj .
Suppose first that H1 and H2 are disjoint. Let (q1 , q2 ) ∈ H1 × H2 be two
points which realize the minimum of kq1 − q2 k. We must have qj ∈ ∂Hj .
Also, the segment joining q1 to q2 must be perpendicular to both ∂H1 and
∂H2 . This situation leads to the result that q1 and q2 are contained on the
great circle C joining p1 and p2 .
We can now reduce everything to a problem in the plane. Let Π be the
plane containing C. We identify Π with C, so that C is the unit circle. Figure
3.1 shows the situation. We have drawn the case when both intersections lie
in a half-disk, but this feature is not a necessary part of the proof.

p2

H2
q2

q1

H1
p1

Figure 3.1: Two circular caps.

Let A be the short circular arc joining p1 and p2 , and let B be the short
circular arc joining q1 and q2 . Referring to Lemma 3.3, we have
DA = D(Q1 , Q2 ) = kp1 − p2 k; (48)
Here D = D(Q1 , Q2 ) is the quantity in the statement of the lemma. We also
have
ψmin ≤ DB = kq1 − q2 k. (49)

19
Finally, the length of the circular arc joining pj to qj is δj τj /2. Therefore
δ1 τ1 /2 + δ2 τ2 /2
δ(A, B) = ∆A − ∆B = = δ(Q1 , Q2 ). (50)
2
Applying Lemma 3.3, and using our three equations, we get
√ 
ψmin ≥ D cos(δ) − 4 − D 2 sin(δ) (51)

But cos(δ) ≥ 1 − δ 2 and sin(δ) ≤ δ. This gives us the lower bound from
Lemma 3.1 in the case the caps are disjoint.
When the caps intersect, we have ψmin = 0. We just want to show that
the lower bound in Lemma 3.1 is nonpositive. We want to reduce this case
to the previous one. Choose some small ǫ > 0 and consider smaller caps H1′
and H2′ that are separated by a distance of exactly ǫ. These two caps are
based on some number δ ′ < δ. Our argument in the previous case gives
√ 
ǫ > D(1 − (δ ′ )2 /2) − 4 − D2 δ′ (52)
√  √ 
′ 2 2 ′ 2 2
D(1 − (δ ) /2) − 4 − D δ > D(1 − δ /2) − 4 − D δ. (53)
Combining these two results, we see that the lower bound in Lemma 3.1 is
less than ǫ. But ǫ is arbitrary. Hence, the lower bound in Lemma 3.1 is
nonpositive in this case. This completes the proof of the lower bound.

The proof of the upper bound is similar. Suppose first that H1 ∩ H2 does
not contain a pair of antipodal points. Then, by Lagrange Multipliers, the
pair of points (q1 , q2 ) realizing the maximum lie on the great circle through
p1 and q2 . We then rotate as above and apply the same argument. This time
we get √ 
ψmax ≤ D cos(δ) + 4 − D 2 sin(δ) (54)
But cos(δ) ≤ 1 and sin(δ) ≤ δ. This gives us the lower bound from Lemma
3.1 in the case H1 ∪ H2 does not contain a pair of antipodal points.
If H1 ∪ H2 contains a pair of antipodal points then one of two things is
true. If p1 and p2 are antipodal points, then obviously the upper bound in
Lemma 3.1 gives a number greater than 2. If p1 and p2 are not antipodal we
can use the same shrinking trick that we used in the previous case to show
that the upper bound in Lemma 3.1 gives a number that is at least 2. In
either case, the upper bound still holds.
This completes the proof of Lemma 3.1.

20
3.3 Proof of Lemma 3.2
Lemma 3.5 Suppose that Q1 is a dyadic segment or square and Q2 = {∞}.
Let p2 = (0, 0, 1). Then, for point p1 ∈ Q∗1 we have ψmin
′ ′
≤ kp1 − p2 k ≤ ψmax .
In particular, the upper bound in Lemma 3.2 is true.

Proof: Consider the lower bound first. The point z1 ∈ Q1 which minimizes
kΣ−1 (z1 ) − Σ−1 (∞)k (55)
is the point furthest from the origin. But, since disks are convex, the point
of Q1 farthest from the origin is a vertex.
Now for the upper bound. The point of Q1 that maximizes Equation 55 is
the one closest to the origin. If a vertex of Q1 does not minimize the distance
to the origin, then (by Lagrange multipliers) some ray through the origin in-
tersects a side of Q1 at a right angle. Since the sides of Q1 are parallel to
the coordinate axes, this can only happen if Q1 crosses one of the coordinate
axes. But Q1 does not cross the coordinate axes. ♠

Lemma 3.5 immediately gives the upper bound. For the lower bound
(which makes a different statement) we need to deal with all the points in the
convex hull H1 = hull(Q∗1 ). Any point q ∈ H1 that minimizes kq − (0, 0, 1)k
must lie on ∂H1 . We claim that ∂H1 is the union of the following
• Q∗1 .
• The flat quadrilateral F that is the convex hull of the vertices of Q∗1 .
• The convex hulls H(Ej ) of the edges E1 , ..., E4 of Q∗1 .
Each of these sets is a subset of H1 and we easily check, in each case, that
every point in each set is a boundary point. Since the union of these sets is
a topological sphere, it must account for the entire boundary.
Our lemma above takes care of the case when q ∈ Q∗1 . If q is an interior
point of F , then the segment joining (0, 0, 1) to q is perpendicular to F . But
F is completely contained in a hemisphere that has (0, 0, 1) on its boundary,
so no such segment can exist.
Finally, suppose that q ∈ H(Ej ). Note that H(Ej ) is contained in a plane
Π which also contains (0, 0, 1). But, considering the distance minimizimation
problem in Π, we see that q cannot be an interior point of H(Ej ). (The
picture looks just like the one drawn in Figure 3.1.) But then either q ∈ Q∗1
or q ∈ F , the two cases we have already handled. This completes the proof.

21
3.4 The Tetrahedral Eliminator
Let Q be a dyadic box, as in the previous section. As in Equation 13, the
quantity ME is the energy of the TBP and the quantity TE is the energy
of the regular tetrahedron. Let Ψmax be the upper bound on the distance
between Q∗1 and Q∗1 as in §3.1.

Lemma 3.6 (Tetrahedral Eliminator) Let Q be a dyadic box. Suppose

that there is some index i such that
X
E(dij ) > ME − TE ; dij = Ψmax (Qi , Qj ). (56)
j6=i

Then no configuration in Q is a minimizer for the E-energy.

Proof: By construction kpi −pj k ≤ dij for all pi ∈ Q∗i and pj ∈ Q∗j . Therefore
X
E(P, pi ) > E(dij ) > ME − TE (57)
j6=i

for any configuration P corresponding to a point in Q. But then every such

configuration violates Equation 13 and cannot be a minimizer. ♠

3.5 Discussion
We can get a cheap Energy Estimator using the fact that
X
E(Z) ≥ E(dij ); dij = Ψmax (Qi , Qj ). (58)
i<j

However, experiments lead us to believe that the main calculation would

effectively take forever using this estimator.
Here is the problem: Even though E(dij ) gives a good estimate for the
minimum energy of a pair of points (pi , pj ) ∈ Q∗i × Q∗j , there is no guarantee
that we can find a single configuration {pj } such that each dij is nearly
realized by E(kpi − pj k). So, the minimum energy of a configuration we
can actually produce might be much higher than our estimate. That is, our
estimate might not be that good globally (for all 10 interactions) even though
it is good locally (for pairs of interactions.) We need to work harder to get a
globally good energy estimator. This is the content of Theorem 5.1.

22
4 The Redundancy Eliminator
The constructions in this chapter work for any power law, and most of the
constructions (just symmetry) work for any decreasing energy function.

4.1 Inversion
An inversion is a conformal involution ρ of C ∪ ∞ that fixes some circle C
pointwise. We call ρ stereo-isometric if Σ−1 (C) is a great circle of S 2 . This
condition means that stereographic projection conjugates ρ to an isometry of
S 2 . The purpose of this section is to prove the following lemma. The reader
might want to just note the result and skip the proof on the first reading.

Lemma 4.1 Let ρ be a stereo-isometric inversion fixing a circle centered on

a point of [1, ∞) ⊂ R. Let X ⊂ C be a strip bounded by horizontal lines,
√ lying below R. Let X− ⊂ X denote those points
one lying above R and one
z = x + iy with x < 1 − 2. Then ρ(X− ) ⊂ interior(X) − X− .

Proof: Being an inversion, the map ρ has 3 nice geometric properties.

• ρ interchanges the disk D bounded by C with its complement.

• ρ interchanges the center of C with ∞.

• ρ preserves each ray through the center of C.

We will use these properties in our proof.

Let z0 be the center of C. We√claim that C separates z0 from any point
on R that lies to the left of 1 − 2. Let p0 = Σ−1 (z0 ). As usual, (0, 0, 1) =
Σ−1 (∞). Let Π be the plane equidistant between p0 and (0, 0, 1). Then

C = Σ(Π ∩ S 2 ) (59)

Let H be the half-space bounded by Π that contains p0 .

The point p0 lies on the great circle connecting Σ−1 (1) = (1, 0, 0) to
(0, 0, 1). Since z0 ∈ [1, ∞), the point p0 lies on the short arc connecting
(1, 0, 0) to (0, 0, 1). But then
√  1 1 
q = Σ−1 (1 − 2) = − √ , 0, − √ 6∈ interior(H). (60)
2 2

23
The extreme cases occurs when p0 = (1, 0, 0). In this extreme case q ∈ Π∩S 2 .
As p0 moves towards (0, 0, 1), H ∩ S 2 moves away from q. Since H does not
contain (0, 0, 1) either, the entire arc connecting q to (0, 0, 1) is disjoint from
the interior of H. But this √ means that C separates z0 from all points of R
that lie to the left of 1 − 2. This proves our claim.
Let D be the disk bounded by C that contains z0 . From what we have
already
√ shown, the interior of D is disjoint from the vertical line L through
1 − 2. Therefore D ∩ X− = ∅. Since ρ swaps D and its complement, we
have

ρ(X− ) ⊂ D. (61)
Hence
ρ(X− ) ∩ X− = ∅. (62)
At the same time, ρ preserves all rays through the z0 , the center of D.
Noting again that z0 ∈ R, we see geometrically that ρ(w) is closer to R than
w for all w ∈ X − D. See Figure 4.1. In particular, we have

ρ(X− ) ⊂ interior(X). (63)

The first statement of the lemma now follows from Equations 62 and 63. The
second statement follows from the first statement and from the fact that ρ
is an involution. ♠

L
D

X_ z0 R
ρ( w)

w X

Figure 4.1: Inversion

24
4.2 A Lemma about Five Points
The result in this section is certainly well known. We suggest that perhaps
the reader just note the result and skip the proof on the first reading.

Lemma 4.2 Let P be any configuration

√ of 5 distinct points on S 2 . Then
some pair points of P are within 2 units of each other.

Proof: A quarter sphere B is a region on S 2 bounded by two semicircular

great arcs, meeting at antipodal points, such that the interior angles at the
intersection points are π/2. The axis of B is the circular arc, contained in
B, that connects the midpoints of the two bounding arcs. We will use the
following easy-to-prove principle: If p is any point on the axis of B, then the
closed hemisphere centered at p contains B. Call this Property X .
The lemma has a trivial proof if P contains two antipodal points, so we
assume this is not the case. Let Hj be the closed hemisphere centered at pj
of P . Another way to state this lemma is that there are indices i 6= j such
that pi ∈ Hj . Assume this is false, for the sake of contradiction.
Let N denote the hemisphere centered at (0, 0, 1), and let S denote the
opposite hemisphere. We normalize so that p0 = (0, 0, 1). By assumption,
p1 , ..., p4 ∈ S − N. Since P does not contain antipodal points, P does not
contain (0, 0, −1). But then, for j = 1, 2, 3, 4, there is a unique point p∗j on
the equator N ∩ S which is is as close as possible to pj .
Let Hj∗ denote the hemisphere centered centered at p∗j . By construction,
Hj∗ ∩ S is a quarter sphere and pj is a point on its axis. By Property X,
we have Hj∗ ∩ S ⊂ Hj . But then, by assumption, pi 6∈ Hj∗ if i 6= j. Let
π : S → R2 be the projection map π(x, y, z) = (x, y). Then π(S) is the unit
disk and π(S ∩ Hj∗) is a half disk. Finally, π(pi) and π(p∗i ) lie on the same
ray through the origin. For this reason, pi ∈ Hj∗ if and only if p∗i ∈ Hj∗ . Since
pi 6∈ Hj∗ , we conclude that p∗i 6∈ Hj∗ .
Now we know that p∗i 6∈ Hj∗ if i 6= j. Hence, the distance from p∗i to p∗j
is greater than π/2. But then we have 4 points p∗1 , ..., p∗4 , contained on the
equator, each of which is more than π/2 from any of the other points. This
is a contradiction. ♠

25
4.3 The Set of Good Configurations
Recall that Ω is our configuration space. Our construction here works for
any energy function E. In this section we define the set Ω′ ⊂ Ω, as discussed
in §2.5 in connection with the Redundancy Eliminator.
Let Z = {z0 , ..., z4 } be a configuration of points in C ∪ ∞, normalized as
in §2.4. We write zj = xi + yi . Also, we set pj = Σ−1 (zj ). Let Ω′ be the set
of configurations Z having all the following properties.
1. kpi − pj k ≤ kp0 − p4 k for all indices i 6= j. Also x0 ≥ 1.
2. y1 ≤ 0 ≤ y2 ≤ y3 .
√
3. If y1 < 0 < y3 then x2 ≥ 1 − 2.
Lemma 4.3 For any Z ∈ Ω, there exists Z ′ ∈ Ω′ such that E(Z ′ ) ≤ E(Z).

Proof: Permuting the points, we can the first half of Property 1. √ By Lemma
4.2 and the first half of Property 1, we must have kp0 , p4 k ≤ 2. But this
fails if z0 ∈ [0, 1). Hence x0 ≥ 1. In short, the first half of Property 1 implies
the second half.
Reflecting in R, we can arrange that y2 ≥ 0. Permuting again, we can
retain Property 1 and arrange that y1 ≤ y2 ≤ y3 . To get Property 2, we just
have to deal with the situation when 0 < y1 . Let Z ′ be the new configuration
obtained when we replace z1 with the conjugate z 1 and otherwise keep the
points the same. Let P and P ′ be the corresponding configurations in S 2 .
The points p0 , ..., p5 are all contained in the hemisphere bounded by the great
circle
C = Σ−1 (R ∪ ∞). (64)
The configuration P ′ is obtained from P simply by reflecting p1 across C.
But, as can be seen from e.g. the Pythagorean theorem, we then have kp′1 −
pj k ≥ kp1 − pj k for j 6 3. The other distances do not change. Since E
is a monotone decreasing function, we have E(Z ′ ) ≤ E(Z). This gives us
Property 2.
It remains
√ to deal with Property 3. Suppose that y1 < 0 < y3 and
x2 < 1 − 2. Let ρ be the stereo-isometric inversion that swaps z0 and z4 .
Let Z ′ = ρ(Z). Let zj′ = ρ(zj ). We set zj′ = x′j + iyj′ . By construction Z ′ has
property 1. By symmetry,
y1′ < 0 < y3′ ; y2′ ≥ 0. (65)

26
There are two cases to consider.

Case 1: Suppose first that

y2′ ≤ y3′ (66)
′
Equations 65 and 66 combine to say that Z has Property 2 as well. Let
X be the strip bounded by the horizontal lines through z1 and z3 . Then X
satisfies the hypotheses of Lemma 4.1. Moreover, z2 ∈ X− , the half-strip
√
from Lemma 4.1. Then z2′ ∈ X − X− by Lemma 4.1. Hence x′2 > 1 − 2.
This shows that Z ′ has Property 3. Hence Z ′ ∈ Ω′ .

Case 2: Suppose that y2′ > y3′ . Then we let Z ′′ be the configuration ob-
tained from Z ′ by swapping z2′ and z3′ . Then Z ′′ satisfies Properties 1 and
2. The argument in Case 1 again shows that z2′ ∈ X − X− . In particular,
z2′ ∈ X. By construction,
z3 6∈ interior(X). (67)
By Lemma 4.1, we have
ρ(X− ) ⊂ interior(X). (68)
Combining these last two equations, we see that
z3 6∈ ρ(X). (69)
Since ρ is an involution, this last equation gives us
z3′ 6∈ X− . (70)
Summarizing the situation, we now know that
z2′ ∈ X; z3′ 6∈ X− . (71)
Since 0 < y3′ < y2′ and z2′ ∈ X, we have
z3′ ∈ X. (72)
Combining the last two equations, we have
√
z3′ ∈ X − X ′ ; =⇒ x′′2 = x′3 ≥ 1 − 2. (73)
The first statement implies the second. The second statement shows that Z ′′
has Property 3 as well. Finally, E(Z ′′ ) = E(Z). ♠

27
4.4 The Main Construction
Now we define the Redundancy Eliminator discussed in §2.5. The redun-
dancy eliminator performs 4 tests, one per property discussed above. Let
Q = Q0 × Q1 × Q2 × Q3 be a dyadic box. As usual, a configuration Z ∈ Q
defines points p0 , ..., p4 , with pj ∈ Q∗j .

Property 1: Let Ψmax and Ψmin be the separation functions from §3.1.
We eliminate Q if
Ψmax (Qi , Qj ) < Ψmin (Q0 , Q4 ) (74)
for some pair of indices i < j such that (i, j) 6= (0, 4). We also eliminate Q
if x0 < 1.

Property 2: Let xj , etc be as in §3.1. We eliminate Q if any of the following

is true.

• y1 ≥ y2 ;

• y2 ≥ y3 .

• y 2 ≤ 0.

• y 1 ≥ 0.

In the first 3 cases, the fact that we are using a weak inequality rather than
a strict inequality means that sometimes we eliminate some configurations
that lie in ∂Ω′ . Since we want to consider every configuration in Ω′ , we need
to justify this. We give the justification in the next section.

Property 3: We eliminate Q if all of the following happen.

• y 1 < 0.

• y 3 > 0.
√
• x2 < 1 − 2.

This time there are no boundary cases to worry about.

28
4.5 Discussion of Boundary Cases
Here we discuss the use of weak inequalities in connection with Property 2 in
the previous section. First of all, the reason why we want to do this is that it
speeds up the computation. For example, were we to keep dyadic boxes with
y 2 = 0 we would need to consider many more dyadic boxes that are near the
TBP.
The justification for why we can use weak inequalities is that all the
configurations in ∂Ω′ that we eliminate are actually counted twice, and we
do not eliminate the relevant dyadic boxes both times.
To clarify the situation, we make the interpretation that our dyadic
squares Q2 and Q3 are missing their top boundaries and Q1 is missing its
bottom boundary. With this convention, the divide and conquer algorithm
from §2.5 examines every point in the configuration space except those for
which |yj | = 2 for some j. These configurations are not minimizers. See the
discussion in §2.4. The main point here is that the union of “quarter-open”
dyadic squares in the subdivision of a “quarter open” dyadic square is still
equal to the original “quarter open” dyadic square. In the case of Q2 and
Q3 , the bottom edges fill in for the top edges. In the case of Q1 the top edges
fill in for the bottom ones.
With this interpretation, we can simply eliminate the dyadic box Q men-
tioned above, because it contains no configurations in Ω′ . Adopting this
convention has no effect whatsoever on our program. It is simply a question
of how we interpret the output of the program.

29
5 The Energy Estimator
5.1 Preliminaries
We fix some value of n. Say that a block Q is a collection Q0 , ..., Qn of
dyadic objects, with Qn = {∞} and all the other objects either segments or
squares. Say that a collection of points z0 , ..., zn of points is dominated by
the block if zk ∈ Qk for all k. In our application, we will take n = 4. In this
case, the set of configurations dominated by a block is precisely a dyadic box.

Remark: Though we always take Qk to be a dyadic segment or square

for k = 0, ..., n − 1, the interested reader will note that everything we do
works the same way with the milder constraint that each Qk is either a seg-
ment in [0, ∞) or a square with sides parallel to the coordinate axes that
does not cross the coordinate axes.

We say that {zk } is a vertex configuration if zk is a vertex of Qk for

k = 0, ..., n. There is a finite list of vertex configurations. Given an energy
function E, we define
E(Q) = min E(Z). E ′ (Q) = min E(Z ′ ). (75)
Z Z′

The first minimum is taken over all configurations Z dominated by Q and

the second minimum is taken over all vertex configurations Z ′ . The purpose
of our main result in this chapter is to bound the quantity
ERR(Q) = E ′ (Q) − E(Q) (76)
in terms of Q.
We are interested mainly in the case when the function E is a power law,
but we will state things more generally just to clarify the formulas we get.
We suppose that E is convex decreasing and satisfies the following property.
If 0 < r2 ≤ r1 then there are constants c ≥ 1 and h ≥ 0 (depending on
everything) such that
E ′ (r2 ) E ′ (r1 )
=c ; E ′′ (r2 ) = (c + h)E ′′ (r1 ). (77)
r2 r1
This condition could be stated more simply, but we have stated exactly the
version we will use in Lemma 6.4, the one place where we use it. Other than
this one place, the rest of our argument works for any convex decreasing
function. The power law functions all satisfy these conditions with h = 0.

30
5.2 The Main Result
b we define
We defined Ψmin in §3.1. Given two dyadic objects Q and Q,
b
R = Ψmin (Q, Q). (78)

When R = 0, our bound gives ∞, a result that holds no matter what. So,
without loss of generality, we treat the case when R > 0. In this case, the
two sets Q∗ and Qb∗ are contained in disjoint convex sets.
We define
b = max(0, Λ1) + Λ2 .
ǫ(Q, Q) (79)
When Q is a dyadic segment,

R ′  1 R2  E ′ (R)
Λ1 = E (R) + − E ′′ (R); Λ2 = − (80)
32 8 32 8
When Q is a dyadic square,

R ′  1 R2  E ′ (R) p p 
Λ1 = E (R) + − E ′′ (R); Λ2 = − 1 + x2 + 1 + y 2 .
16 4 16 7.98
(81)
Let δi be as in Equation 41, relative to Qi . Recall that δi is a good estimate
for the side length of of the spherical patch Q∗i . Also, let ǫij = ǫ(Qi , Qj ).

Theorem 5.1 (Energy Estimator) Let Q = (Q0 , ..., Qn ) be any block.

Then
n−1 X
X
ERR(Q) ≤ ǫij δi2 .
i=0 j6=i

Again, we remark that the case n = 4 is the case of interest to us. For
reference, we work out the power-law case E(r) = r −e explicitly. When Q is
a dyadic segment,

e(e + 1) e(e + 2) e
Λ1 = − ; Λ2 = . (82)
8Re+2 32Re 8Re+1
When Q is a dyadic square,

e(e + 1) e(e + 2) e p p 
2 2
Λ1 = − ; Λ2 = 1 + x + 1 + y (83)
4Re+2 16Re 7.98Re+1

31
5.3 The Beginning of the Proof
Recall that a partition of unity is a collection of functions that sum to 1 at
every point. A point weighting assigns a partition of unity

λab : Q → [0, 1]; a, b ∈ {0, 1} (84)

to each dyadic square Q and a partition of unity

λa : Q → [0, 1]; a ∈ {0, 1} (85)

to each dyadic segment. We have a specific point weighting in mind, and we

will define it in §8.
We let the vertices of a dyadic square Q be Qab for a, b ∈ {0, 1}. Here Q00
is the lower left vertex and the remaining vertices, traced in counterclockwise
order, are Q10 , Q11 , Q01 . We let the vertices of a dyadic segment Q be Q0
and Q1 , with Q0 on the left.
Given two points z, w ∈ C ∪ ∞, we define
1
f (z, w) = . (86)
kΣ−1 (z) − Σ−1 (w)ke

b as in the previous section.

Let X denote the set of disjoint pairs (Q, Q),
Let ǫ : X → R be as in the previous section. We will spend the next 3
chapters proving the following result.

Lemma 5.2 There exists a point weighting such that the following is true
b ∈ X.
for all (Q, Q)

• When Q is a segment
X 
b δ(Q)2 ;
λa (z)f (Qa , w) −f (z, w) < ǫ(Q, Q) b
∀(z, w) ∈ Q× Q.
a

• When Q is a square
X 
b δ(Q)2 ;
λab (z)f (Qab , w) − f (z, w) < ǫ(Q, Q) b
∀(z, w) ∈ Q × Q.
a,b

32
Proof of Theorem 5.1: For ease of exposition, we will assume that the
dyadic objects Q0 , ..., Qn−1 are all dyadic squares. The case when there are
some segments involved presents only notational complications. Consider
a configuration z0 , ..., zn , domainated by Q, that realizes E(Q). For each
i = 0, ..., n, define

Bi = {z0 } × {zi−1 } × Qi × ... × Qn . (87)

Just to be explicit, the two extreme cases are

B0 = Q0 × ... × Qn ; Bn = {z0 } × ... × {zn }. (88)

(Recall that Qn = {zn }.) The set Bi is a (2n − 2i)-dimensional rectangular

solid. Let E 1 (i) denote the minimum energy taken over all configurations of
Bi . Let E 2 (i) denote the minimum energy over all vertex configurations of
Bi . The conclusion of this lemma is exactly
n−1 X
X
E 2 (0) − E 1 (0) ≤ ǫij δ(Qi )2 (89)
i=0 j6=i

We have E 2 (n) = E 1 (n), because Bn is a single point. Hence, we have the

telescoping sum
n−1 
X   
E 2 (0) − E 1 (0) = E 2 (i) − E 1 (i) − E 2 (i + 1) − E 1 (i + 1) (90)
i=0

Our choice of {zi } as a minimizing configuration implies that E 1 (i) is

independent of i Combining this with Equation 90, we get
n−1
X
E 2 (0) − E 1 (0) = E 2 (i) − E 2 (i + 1). (91)
i=0

To establish Equation 89, it suffices to prove

X
E 2 (i) − E 2 (i + 1) ≤ ǫij δi2 ; ∀i = 0, ..., n − 1. (92)
j6=i

For the remainder of the proof, we fix some i ∈ {0, ..., n − 1} once and for
all. We can find vertices vi+1 , ..., vn of Qi+1 , ..., Qn respectively such that

E(z1 , ..., zi , vi+1 , ..., vn ) = E 2 (i + 1). (93)

33
Here E(...) is the energy of the configuration. Note that Qn is a singleton,
so that our “choice” of vn is forced.
For any choice of a, b ∈ {0, 1}, let Qiab be the corresponding vertex of Qi .
Let
λab = λab (zi ) (94)
be as guaranteed by Lemma 5.2. The function λab depends on i, but we
suppress this from our notation.
Since E 2 (i) is the minimum energy of any vertex configuration of Bi ,

E(z1 , ..., zi−1 , Qiab , vi+1 , ..., vn ) ≥ E 2 (i). (95)

Hence,
E 2 (i) − E 2 (i + 1) ≤ Aab ;
Aab = E(z1 , ..., zi−1 , Qiab , vi+1 , ..., vn ) − E(z1 , ..., zi , vi+1 , ..., vn ). (96)
Setting wj = zj for j < i and wj = vj for j > i, we have

E 2 (i) − E 2 (i + 1) ≤
X
λab Aab =1
a,b
X X 
λab f (Qiab wj ) − f (zi , wj ) =2
a,b j6=i
XX  
λab f (Qiab , wj ) − f (zi , wj ) =3
j6=i a,b
XX 
λab f (Qiab , wj ) − f (zi , wj ) ≤
j6=i a,b
X
ǫij δi2 . (97)
j6=i

P The first inequality comes from Equation 96 and from the fact that
λab = 1. Equality 1 comes from the cancellation of all terms not in-
volving the ith index when we subtract the two sums for Aab . Equality 2
comes fromP switching the order of summation. Equality 3 comes from the
fact that ab λab = 1. The last inequality follows from Lemma 5.2. This
establishes Equation 92, which is all we need to prove Theorem 5.1. ♠

34
6 An Estimate for Line Segments
6.1 The Main Estimate
In this chapter we prove an estimate that relates directly to the Λ1 term in
our Energy Estimator. The reader might want to simply note the main result
in this chapter on the first reading, and then come back to the proof later
on. Let E be as in the previous chapter.
Let p ∈ S 2 be some point. In terms of Lemma 5.2, we think of p as being
some point in the spherical patch (Q) b ∗ . Let A′ ⊂ R3 − {p} be a segment
whose endpoints lie in S 2 . In terms of Lemma 5.2, we think of A′ as joining
two boundary points of the spherical patch Q∗ . We define

F (q) = E(kp − qk) (98)

Here F depends on p, but we suppress this from our notation.

Define
A′x = (1 − x)A′0 + xA′1 ; x ∈ [0, 1]. (99)
Then x → A′x is a constant speed parameterization of A′ .

Lemma 6.1 (Segment Estimate) Suppose that p is at least R units from

every point of A′ . Let δ be the length of A′ .

Xδ 2
(1 − x)F (A′0 ) + xF (A′1 ) − F (A′x ) ≤ ,
8
where 
R ′ R2  ′′
X = E (R) + 1 − E (R).
4 4
Remark: We wish to point out one unfortunate feature of our notation.
The quantities E ′ and E ′′ are derivatives of E whereas the quantity A′ is
simply a chord of S 2 . We make this notation because, in the next chapter,
we will consider an arc A of S 2 and the chord A′ that joins the endpoints of
A. In some sense, the chord A′ is a linear approximation to the arc A and
the derivative E ′ is a linear approximation to the function E.

35
6.2 Strategy of the Proof
Lemma 6.1 really just involves the single variable function φ(x) = F (A′x )
defined for x ∈ [0, 1]. In §6.3 we prove the following easy estimate.

Lemma 6.2
H
(1 − x)φ(0) + xφ(1) − φ(x) ≤ ; H = sup φ′′ (x).
8 x∈[0,1]

Here φ′′ (x) is the second derivative with respect to x.

Recall that δ is the length of S. Let s denote the arc-length parameter
along the segment S. We set things up so that s = 0 corresponds to S0 and
s = δ corresponds to S1 . In general, the parameter s ∈ S corresponds to
x = s/δ ∈ [0, 1]. In §6.5 we establish the following result.

Lemma 6.3 Suppose we have the hypotheses in Lemma 6.1. Then,

d2 φ R ′  R2  ′′
≤ E (R) + 1 − E (R).
ds2 4 4
By the Chain Rule, we have

d2 φ
φ′′ (x) = δ 2 . (100)
ds2
at corresponding points.
Lemma 6.1 follows from Lemma 6.2, Lemma 6.3, and Equation 100.

6.3 Proof of Lemma 6.2

We will prove Lemma 6.2 for a smooth function h : [0, 1] → R. We want to
show that
H
(1 − x)h(0) + xh(1) − h(x) ≤ ; H = sup h′′ (x). (101)
8 x∈[0,1]

Let c be any nonzero constant and let L be any linear function. Equation
101 holds for the function h if and only if it holds for ch. Likewise, Equation
101 holds for h − L if and only if Equation 101 holds for h. Using these two

36
symmetries, it suffices to prove Equation 101 in the case when h(0) = h(1) =
0 and H = 1. In this case, Equation 101 simplifies to
1
−h(x) ≤ . (102)
8
Let a ∈ [0, 1] be a point where −h attains its maximum. That is, h attains
its minimum at a. Replacing h by the function x → h(1 − x) if necessary, we
can suppose without loss of generality a ≥ 1/2.
We have h′ (a) = 0 and, by the Fundamental Theorem of Calculus,
Z b
′
h (x) = h′′ (t)dt ≤ x − a (103)
a

for all x ∈ [a, 1]. But then

Z 1 Z 1
′ (1 − a)2 1
−h(a) = h(1) − h(a) = h (x)dx ≤ x − a dx = ≤ (104)
a a 2 8
This completes the proof. ♠

6.4 Comparison Lemmas

Let p and f be as above. Say that a flag is a pair (q, L) where L is a line
and q ∈ L is a point, and q 6= p. We define the following quantities.
• r(q, L) is the distance from p to q.
• d(q, L) is the distance from p to L.
• θ(q, L) is the small angle between pq and L.
• D(q, L) = F ′′ (q), the second derivative w.r.t. arc length on L.
We want to compare flats (q1 , L1 ) and (q2 , L2 ). We set θ1 = θ(q1 , L1 ), etc.
Our main result is Corollary 6.6, proved at the end. This result will help us
establish Lemma 6.3.
For each j we can rotate so that everything takes place in C, and pj = irj
and Lj = R, and qj = xj . Here x2j + d2j = rj2 . Suppressing the index, we have
d2 √ ′′ x2 E ′ (r) d2
D= E( d 2 + x2 ) = E (r) + × 2. (105)
dx2 x r2 r r
′ ′′
Here E and E are derivatives taken with respect to r.

37
Lemma 6.4 If θ2 = θ1 and r2 ≤ r1 and D1 > 0 then D2 ≥ D1 .

Proof: In our situation, there is some constant ρ ∈ (0, 1] so that

r2 x2 d2
= = = ρ. (106)
r1 x2 d1
It follows from Equation 77 that there are constants c ≥ 1 and h ≥ 0 such
that
E ′ (r2 ) E ′ (r1 )
=c ; E ′′ (r2 ) = (c + h)E ′′ (r1 ). (107)
r2 r1
It now follows from Equation 105 that

D2 − cD1 = hE ′′ (r1 ) ≥ 0. (108)

Hence D2 ≥ cD1 for some c ≥ 1. The lemma follows immediately. ♠

Lemma 6.5 If θ2 ≤ θ1 and r2 = r1 then D2 ≥ D1 .

The quantity d(q, L) is monotone increasing with θ(q, L). Thus, we have
d2 ≤ d1 . This time, we have

r := r2 = r1 ; d2 ≤ d1 ; x2 ≥ x1 . (109)

We again have Equation 105. Note that E ′′ (r) > 0 and E ′ (r) < 0. When
we change from the index j = 1 to the index j = 2, we do not decrease
the positive coefficient of E ′′ (r) in the first term and we do not increase the
positive coefficient of E ′ (r)/r in the second term. Hence, D2 ≥ D1 . ♠

The previous two results combine to prove the following corollary.

Corollary 6.6 If θ2 ≤ θ1 and r2 ≤ r1 and D1 ≥ 0, then D2 ≥ D1 .

6.5 Proof of Lemma 6.3

We continue using the notation from above. There is a unique plane Π such
that p ∪ L ⊂ Π. We rotate the picture so that Π is the xy-plane. Let
C = S 2 ∩ Π. Then, as shown in Figure 6.1,

38
 1. C is a circle whose radius is at most 1;

2. Every point of S is at least R units from p.

θ
S q

Figure 6.1: The chord and the circle

We are interested in bounding the quantity D(q, L), where L is the line
containing the chord S. The chord S is subject to the two constraints men-
tioned above. If C has radius less than 1, we replace C by a unit radius
circle C ′ , and S by a larger segment S ′ such that the pair (C ′ , S ′ ) satisfies
the same constraints, and the flag (q, L′ ) is the same as the flag (q, L). Figure
2.2 shows the construction.

39
 p

C’ C
L
q
S
S’

Figure 6.2: Expanding the circle

So, without loss of generality, we can assume that C is the unit circle.
Also, it suffices to consider the case when D(q, L) ≥ 0. Let q1 = q. Let r1
be the distance from p to q1 . Note that r1 ≥ R. Let q2 denote a point on C
that is exactly R = r2 units from p. Let L2 be the line tangent to C at q2 .
We want to apply Corollary 6.6 to the flats (q1 , L1 ) and (q2 , L2 ). We already
know that r2 ≤ r1 .

Lemma 6.7 θ2 ≤ θ1 .

Proof: See Figure 6.3. Let q3 be the endpoint of S such that the small
angle θ1 subtends the arc of C between p and q3 , as shown in Figure 6.3.
Let L3 be the line tangent to C at q3 . Let θ3 = θ(q3 , L3 ). The angle θ1 is
half the length of the two thick arcs in Figure 6.3 whereas the angle θ3 is
half thelength of the thick arc joining p to q3 . Hence θ3 ≤ θ1 . But the angle
θ3 decreases as we move q3 towards p along C. Therefore θ2 ≤ θ3 . Putting
these two inequalities together, we find that θ2 ≤ θ1 . ♠

40
 L2
p
θ2

q2 L3
C
θ3 L1

q3
θ1
q1

Figure 6.3: Comparing two flags.

Corollary 6.6 now says that D(q2 , L2 ) ≥ D(q1 , L1 ). To finish our proof,
it remains only to compute the quantity D(q2 , L2 ).
We know that r2 = R. Some elementary geometry shows that
R2
d2 = ; x2 = r 2 − d2 (110)
2
Plugging Equation 110 into Equation 105 and simplifying, we get
R  R2  ′′
D(q2 , L2 ) = E ′ (R) + 1 − E (R), (111)
4 4
the bound from Lemma 6.3. Now we know that
R  R2  ′′
D(q, L) ≤ E ′ (R) + 1 − E (R), (112)
4 4
for all flags (q, L) such that kp − qk ≥ R. But D(q, L) is just another name
for the quantity d2 φ/ds2 featured in Lemma 6.3. This completes the proof
of Lemma 6.3.

41
7 Parametrizing Arcs and Segments
7.1 Overview
Let S ⊂ C be a line segment. In this chapter we define a certain parametriza-
tion of S by a parameter x ∈ [0, 1]. Once we define our parametrization, we
will state and prove several geometric results about it. As with the last chap-
ter, the reader might want to just note the results on the first reading and
then come back later for the proofs.
Let A be a circular arc on S 2 , and let A′ be the chord that joins the
endpoints A0 and A1 of A. Throughout the chapter, we assume that A is
contained in a semicircle. We let x → A′x be the affine map from [0, 1] to
A′ . Let C be the circle containing A, and let c ∈ C be the point which is
diametrically opposed to the midpoint of A. We define Ax so that the three
points c, A′x , Ax are always collinear. Here is our first result.
Lemma 7.1 The function f (x) = kAx − A′x k attatains its maximum at x =
1/2.
Now we return to our main task of parametrizing a segment S ⊂ C. Let
S ∗ = Σ−1 (S). The method above gives us a parametrization of S ∗ . Now we
define
Sx = Σ(Sx∗ ). (113)
As usual Σ denotes stereographic projection.
Suppose Q is a normal dyadic square. This means that Q does not cross
the coordinate axes, and the side length of Q is at most 1. We consider the
case when Q is contained in the positive quadrant. The other cases have
symmetric treatments. Let Q0 and Q1 be the left and right edges of Q. For
any x ∈ [0, 1], let Qx denote the segment connecting (Q0 )x and (Q1 )x . The
main result in this chapter gives estimates on the size and shape of the image
of Q∗x = Σ−1 (Qx ).
Lemma 7.2 Q∗x has arc length at most
δ × 1.0013.
and is contained in a circle of radius at least
1
p .
1 + y2

42
7.2 Proof of Lemma 7.1
Our result is scale-invariant. It suffices to prove the result when A is an arc
of the unit circle, as shown in Figure 7.1. The arc cy is evidently shorter
than the diameter cx. On the other hand, the arc cz is evidently longer than
the arc cw. Hence the arc yz is shorter than the arc wx. This is what we
wanted to prove.

z y

c w x

Figure 7.1: The relevant points

7.3 A Lemma about Slopes

The rest of the chapter is devoted to proving Lemma 7.2. Here we reduce
Lemma 7.2 to a more geometric statement.

Lemma 7.3 Qx has slope in (0, 0.051) for all Q and x ∈ (0, 1).

Proof of Lemma 7.2: Let S = Qx . We deal first with the arc length of
S ∗ . This is really the same argument as in Lemma 3.4. When evaluated at
all points of S, the quantity in Equation 7 is at most δ/s. Since S has slope
in (0, 0.051) and Q has side length s, the segment S has length at most
p
s × 1 + (0.051)2 < 1.0013. (114)

43
The arc length estimate on S ∗ follows by integration.
Since the line L through S has positive slope, it intersects the imaginary
axis in a point of the form iy where y < y. But then, according to Equation
6, there are two points on (L ∪ ∞)∗ which are at least
1
p
1 + y2

apart. ♠

The rest of the chapter is devoted to proving Lemma 7.3.

7.4 Mobius Geometry

Let the vertices of Q be Q00 , Q10 , Q11 , Q01 , starting in the bottom left corner
and going counterclockwise around. Associated to Q is an auxilliary map
φ : [0, 1] → [0, 1], defined as follows.

• The point (1 − x)Q00 + xQ01 is (Q0 )s for some s.

• The point (Q1 )s is (1 − y)Q10 + yQ11 for some y = φ(x).

In other words, we consider the natural map (Q0 )s → (Q1 )s but we precom-
pose and postcompose with affine maps to make the domain and range equal
to [0, 1]. Lemma 7.2 is equivalent to the statement that

0 < φ(x) − x < 0.051; ∀x ∈ (0, 1). (115)

Given 4 distinct points A, B, C, D ∈ Rn , we have the cross ratio

kA − Ck kB − Dk
χ(A, B, C, D) = . (116)
kA − Bk kC − Dk

We say that a homeomorphism from one curve to another is Mobius if the

map preserves cross ratios. In case the curves are line segments in the plane,
a Mobius map between them is the restriction of a linear fractional transfor-
mation of C ∪ ∞.

Lemma 7.4 φ is Mobius.

44
Proof: Since similarities are Mobius transformations, it suffices to prove
that the map (Q0 )s → (Q1 )s is a Mobius transformation. Stereographic
projection is well known to be a Mobius map from any line segment in C to
the corresponding arc on S 2 . Referring to the construction in the beginning
of §7.1, the map from Ax to A′x is just the composition of affine maps with
(one dimensional) stereographic projection. Hence, this map is also Mobius.
Let A0 = Q∗0 and A′0 be the chord connecting the endpoints of A0 . Like-
wise define A1 and A′1 . The map of interest to us is the composition

Q0 → A0 → A′0 → A′1 → A1 → Q1 . (117)

The outer maps are Mobius, from what we have already said, and the middle
map is affine. ♠

A Mobius map from [0, 1] to [0, 1] is completely determined by its deriva-

tive at either endpoint. Thus, we can understand φ by computing or esti-
mating φ′ (0). In our next result, we think of φ′ (0) as a function of the choice
of dyadic square Q. As in Equation 7, define
2
g(x, y) = . (118)
1 + x2 + y 2
Define
g(Q00 )g(Q11 )
GQ = . (119)
g(Q10 )g(Q10 )

Lemma 7.5 Relative to Q, we have

p
φ′ (0) = GQ .

In particular, φ is the identity iff GQ = 1.

Proof: If follows from symmetry that

φ′ (0)φ′(1) = 1. (120)

Therefore s
φ′ (0)
φ′ (0) = . (121)
φ′ (1)

45
Looking at the composition in Equation 117, all the maps except the outer
two have the same derivative at either endpoint. For the affine map in the
middle, this is obvious: the derivative is constant. In the case of the map
Aj → A′j this follows from the fact that we are projecting from a point cj
that is symmetrically located with respect to A′j and Aj . Call this property
of the derivatives the symmetry property.
By Equation 7, the quantity g(Q00 ) is the norm of the derivative of Σ−1
at Q00 . The other quantities g(Qij ) have similar interpretations. It therefore
follows from the symmetry property and the Chain Rule that
φ′ (0) g(Q00 )/g(Q01 )
′
= . (122)
φ (1) g(Q10 )/g(Q11 )
This Lemma now follows from Equations 121 and 122. ♠

7.5 The End of the Proof

For the purposes of doing calculus, we define φ and G relative to any square
that is contained in the positive quadrant. We remind the reader that we
only consider squares that have side length at most 1.

Lemma 7.6 It never happens that GQ = 1.

Proof: When Q00 = (x, y) and Q has side length r, we compute

(1 + r 2 + 2rx + x2 + y 2 )(1 + r 2 + 2ry + x2 + y 2)

GQ = (123)
(1 + x2 + y 2)(1 + 2r 2 + x2 + y 2 + 2rx + 2ry)
Every factor in Equation 123 is a polynomial with only constant coefficients
and at least one constant term. Hence this expression never vanishes. ♠

When Q = [0, 1]2 , we compute that GQ = 4/3. Hence φ′ (0) > 1 relative
to this choice of Q. It tollows from continuity that φ′ (0) > 1 relative to any
square in the positive quadrant. But this means that φ is increasing. (Here,
of course, we are crucially using the fact that φ is a Mobius map from [0, 1]
to [0, 1].) Hence φ(x) > x. This proves that the slope of the segment Qx is
positive for all Q and all x > 0. This proves half of Lemma 7.3. Now we
turn to the other half.

46
Lemma 7.7 The quantity GQ is maximized when Q has side length 1 and
√
3−1
Q00 = (ξ, ξ); ξ= .
2
Proof: Let ψ(x, y, r) be the function in Equation 123. We compute symbol-
ically that
dψ ∆(x, y, r)
= , (124)
dr (1 + x + y ) (1 + 2r 2 + x2 + y 2 + 2rx + 2ry)2
2 2 2

Where ∆ is a polynomial with entirely positive terms, at least one of which

involves only r. (The polynomial is rather long and unenlightening.) From
this we conclude that dψ/dr > 0. Hence, the maximum value of GQ must
occur when r = 1.
We compute that ψ(ξ, ξ, 1) = 3/2. So, we just have to show that the
function h(x, y) = 3/2 − ψ(x, y, 1) is non-negative on the positive quadrant.
We compute that h is a rational function. The denominator is a polynomial
with only positive terms, and the numerator N equals
1−2x+4x2 +2x3 +x4 −2y −8xy +2x2 y +4y 2 +2xy 2 +2x2 y 2 +2y 3 +y 4 (125)
In fact, N is non-negative on the entire plane. To see this, we make the
change of variables
x = ξ + u; y = ξ + v. (126)
With this change of variables, we find that N − (2u − 2v)2 is a polynomial
involving only positive terms. ♠

In light of the previous result, the quantity φ′ (0) is maximized for the
special square Q0 from Lemma 7.7. Since G = 3/2 in this case, we have
p
φ′ (0) = 3/2. (127)
But this equation pins down φ uniquely, and we observe that the map
√
3 2x
φ(x) = √ √ √ . (128)
2 3 + 3 x − 2 3x
has the same derivative. Hence, this is the correct formula for φ. A bit of
calculus now shows that
φ(x) − x < 0.051; ∀x ∈ [0, 1]. (129)
This completes the proof.

47
8 Proof of Lemma 5.2
8.1 The Geometry of Circles
We need one more result about circles.

Lemma 8.1 Suppose that A is an arc of a circle C. Let d be the arc-length

of A. Let A′ be the segment connecting the endpoints of A. Let r be the radius
of C. Let µ be the smallest constant such that every point of A is within µ
of some point of A′ . Then
d2
µ< .
8r

Proof: Let’s first consider the case r = 1. We rotate so that C is the unit
circle, and A is the arc bounded by the points exp(−iθ) and exp(iθ). Here
θ ∈ (0, π). Then
d = 2θ; µ = 1 − cos(θ). (130)
The claim of this lemma boils down to the statement that
θ2
> 2, (131)
1 − cos(θ)

This is equivalent to the statement that

φ(θ) = θ2 + 2 cos(θ) − 2 > 0. (132)

We have φ(0) = 0 and

φ′ (θ) = 2(θ − sin(θ) > 0. (133)

This proves what we need.

If r 6= 1, we let T be a dilation that scales distances by a factor of 1/r.
The arc T (A) has length d/r and T (C) has radius 1. We apply our result to
the pair (T (A), T (C)) and find that every point of T (A) is within

(d/r)2
8
of T (C). Applying T −1 , we get the desired result for the pair (A, C). ♠

48
8.2 Lemma 5.2 for dyadic segments
8.2.1 Defining the Weighting
Suppose that (Q, Q)b is a reasonable pair, and Q is a dyadic line segment.
The construction in the previous chapter gives us a parametrization x → Qx .
We take Q0 to be the left endpoint and Q1 to be the right endpoint. The
corresponding endpoints of Q∗ are Q∗0 and Q∗1 . We define our weighting as
follows. Letting z = Qx , we define

λ0 (z) = 1 − x; λ1 (z) = x; (134)

8.2.2 Setting up the Calculation

To bring our notation in line with Lemma 6.1, we define

A = A(Q∗0 , Q∗1 ) = Q∗ ; A′ = A′ (Q∗0 , Q∗1 ). (135)

b Define
Let (z, w) ∈ Q × Q.

p = Σ−1 (w); q = Σ−1 (z). (136)

Setting x = λ1 (z), we have

q = Ax . (137)
We also define
q ′ = A′x (138)
The idea of our proof is to estimate things with q ′ in place of q, and then to
estimate the error we get when replacing q ′ by q.

8.2.3 Using Lemma 6.1

Let F be as in Equation 98. We have
X
λu (z)f (Qu , w) = (1 − x)F (A′0 ) + xF (A′1 );
u

F (A′x ) = E(kq ′ − pk). (139)

Lemma 6.1 now tells us that
X
λu (z)f (Qu , w) − E(kp − q ′ k) ≤ max(0, Λ1 ) δ 2 (140)
u

49
8.2.4 The Easy Case
Now we need to see what happens when we replace q ′ by q. Define

r = kp − qk r ′ = kp − q ′ k. (141)

If r ′ ≥ r then (since E is decreasing)

E(kp − qk) ≥ E(kp − q ′ k). (142)

In this case, our proof is done: Equations 142 and 140 combine to give a
tighter bound than what Lemma 5.2 gives.

8.2.5 The Hard Case

Now suppose r ′ < r. Since E is convex, E ′ is monotone decreasing. Therefore

E(r ′ ) − E(r) ≤ E ′ (r ′)kq − q ′ k ≤ E ′ (R)kq − q ′ k. (143)

Here R ≥ r ′ is as in Lemma 5.2.

The same proof as in Lemma 7.2 shows that A has arc length at most
δ. Also A is contained in a great circle – i.e. a circle of radius 1. Lemma
8.1 now says that every point of A is within δ 2 /8 of A′ . But the point of A′
closest to A1/2 is A′1/2 . Therefore

δ2
kq − q ′ k = kAx − A′x k ≤∗ kA1/2 − A′1/2 k ≤ . (144)
8
The starred inequality is Lemma 7.1. Combining Equations 143 and 144, we
find that
δ2
E(r ′ ) − E(r) ≤ E ′ (R) = Λ2 δ 2 . (145)
8
Note that
f (z, w) = E(kp − qk). (146)
Hence, by Equation 146,

E(kp − q ′ k) − f (z, w) = E(r ′) − E(r) ≤ Λ2 δ 2 . (147)

Adding Equations 140 and 147, we get the bound in Lemma 5.2. This com-
pletes the proof in case Q is a dyadic segment.

50
8.3 The Weighting for Dyadic Squares
Let Q be a dyadic square and let Q∗ = Σ−1 (Q). Let Qab be the vertices of
Q, as in §7.4. Let Q∗ab be the corresponding vertex of Q∗ . As we make our
construction, the reader should picture the letter ‘H’, with the horizontal bar
very slightly slanted. We will use coordinates (h, v) ∈ [0, 1]2. The v vari-
able moves along vertical segments and the h variable moves along (roughly)
horizontal segments.
Let Q0 be the left edge of Q. The two endpoints of Q0 are Q00 and Q01 .
Likewise, let Q1 be the right edge of Q. The two endpoints of Q10 are Q11 .
Using the parametrization from the previous chapter, we define

Q0v = (Q0 )v ∈ Q0 ; Q1v = (Q1 )v ∈ Q1 . (148)

Next, we define Qv to be the segment joining Q0v to Q1v . Note that Qv is

nearly horizontal but not exactly horizontal. See Lemma 7.3. Finally, we
define
Qhv = (Qv )h , (149)
again using the parametrization discussed in the previous sections. We define
our weighting as follows. Assuming that z = Qvh ,

λ00 (z) = (1 − h)(1 − v); λ10 (z) = (h)(1 − v);

λ01 (z) = (1 − h)(t); λ11 (z) = (h)(v). (150)

In the previous section, we defined a weighting for dyadic segments. We
could equally well define this weighting for any segment, since we have already
explained how to parametrize any such segment. There are three ways in
which our weighting here is compatible with the weighting for segments.

• With respect to the segment Q0 , the weighting for the point S0v is
λ0 = 1 − v and λ1 = v.

• With respect to the segment Q1 , the weighting for the point S1v is
λ0 = 1 − v and λ1 = v.

• With respect to the segment Qv , the weighting for the point Qhv is
λ0 = 1 − h and λ1 = h.

51
8.4 The End of The Proof
For ease of notation, we will assume that our dyadic square lies in the positive
quadrant. The other cases are similar, and indeed follow from symmetry. We
gather 4 pieces of information.

1. The circles containing the arcs Q∗0 and Q∗1 have radius at least
1
√ .
1+x
This follows from Equation 6 and from the fact that the line extending
a vertical edge of Q comes within x of the origin.

2. The arcs Q∗0 and Q∗1 have length at most δ. This follows from the same
argument as in Lemma 7.2.

3. For any v ∈ [0, 1], the line extending the segment Qv comes within y
of the origin. Hence, the circle containing Q∗v has radius at least
1
√ .
1+y
See Lemma 7.2.

4. For any v ∈ [0, 1], the arc Q∗v has length at most (1.0013) × δ. See
Lemma 7.2.

Now we are ready for the main argument. The basic idea is to make
repeated appeals to the segment case of Lemma 5.2 and then to suitably
average the result.
We define
Λ1x = Λ1y = Λ1 /2;
E ′ (R) p E ′ (R) p
Λ2x = − 1 + x2 ; Λ2y = − 1 + y2. (151)
8 7.98
We have
Λ1 = Λ1x + Λ1y ; Λ2 = Λ2x + Λ2y . (152)
b ∗ be some point. Let w = Σ(p). Let f be as in Lemma 5.2.
Let p ∈ (Q)
That is
f (z, w) = E(kΣ−1 (z) − Σ−1 (w)k). (153)

52
 Applying the segment case of Lemma 5.2 to the arc Q0 , we find that
 
Z1 := (1−v)f (Q00 , w)+vf (Q01, w)−f (Q0v , w) ≤ max Λ1x +Λ2x δ 2 (154)

The proof is exactly the same as in the previous section, except for the one
point that the circular arc Q∗0 lies not necessarily in a great circle but rather
a circle whose radius is bounded by Item 1 above. The designation Z1 is for
algebraic purposes which will become clear momentarily. Similarly
 
Z2 := (1−v)f (Q10 , w)+vf (Q11, w)−f (Q1v , w) ≤ max Λ1x +Λ2x δ 2 (155)

Finally, an argument just like the one given for dyadic segments also works
for the segment Qs . The only property we used about dyadic segments is
that they don’t cross the coordinate axes, and Qs has this property. Using
Items 3 and 4 above in place of Items 1 and 2, the same argument gives

 
Z3 := (1−h)f (Q0v , w)+hf (Q1v , w)−f (Qhv , w) ≤ max Λ1y +Λ2y δ 2 (156)

The key point here is that

7.98 < 8/1.0013.
The number 1.0013 comes up in Item 4 above.
Concentrating on the left hand sides of Equations 154, 155, and 156, we
have
(1 − h)Z1 + hZ2 + Z3 =
(1 − h)(1 − v)f (Q00 , w) + (1 − h)(v)f (Q01 , w)+
(h)(1 − v)f (Q10 , w) + (h)(v)f (Q11 , w) − f (Qhv , w) =
X 
λab (z)f (Qab , w) − f (z, w); z = Qhv . (157)
ab
Concentrating on the right hand sides, we see that the quantity in Equation
157 is at most
     
(1 − h) max Λ1x + Λ2x δ 2 + (h) max Λ1x + Λ2x δ 2 + max Λ1y + Λ2y δ 2 =
 
max(Λ1 , 0) + Λ2 δ 2 . (158)
This proves Lemma 5.2.

53
9 The Hessian and its Variation
9.1 Main Part of the Proof
In this chapter we prove Lemma 2.4. Let He denote the Hessian of E, the en-
ergy function, relative to the function E(r) = r −e . Let Z be the configuration
corresponding to the TBP, normalized as in Equation 30.
Lemma 9.1 The lowest eigenvalue of He exceeds 1/10 for e = 1, 2.

Proof: Let M be either of the two matrices. Let I7 be the 7 × 7 identity ma-
trix. Using a modified version of the Cholesky Decomposition, as discussed
in [Wa, p 84], we write
1
(M − I7 ) = LDLt . (159)
10
where L is lower triangular, D is diagonal matrix with all diagonal entries
1
positive, and Lt is the transpose of L. This suffices to show that M − 10 I7 is
positive definite. Hence, the lowest eigenvalue of M is at least 1/10. ♠

Given a square matrix M, we define

X 1/2
kMk2 = Mij2 . (160)
i,j

We mention a familiar and useful property of this norm.

Lemma 9.2 For any unit vector v ∈ R7 we have kM(v) · vk ≤ kMk2 .

Proof: We have
X X 1/2
kM(v) · vk2 = Mij vi vj ≤∗ kMk2 (vi vj )2 =
2
r  X 
X
kMk2 vi2 vj2 = kMk2 kvk2 = kMk2 . (161)

The starred inequality is the Cauchy-Schwarz inequality. ♠

We write M ′ ≺ M if |Mij′ | ≤ Mij for all indices. Recall that s = 2−11 , as

in Lemmas 2.3 and 2.4. Let Dk denote the partial derivative with respect to
the kth direction in R7 .

54
Lemma 9.3 (Variation Bound) Let Ψ1,e , ..., Ψ7,e be the smallest non-negative
matrices such that Dk He (Z) ≺ Ψk,e for all k and all Z ∈ Ωs . Let
7
X
Ψ(e) = Ψk,e .
2
k=1

Then Ψ(1) < 345 and Ψ(2) < 140 and, supe Ψ(e) < 463.

Corollary 9.4 Let e = 1 or e = 2. Let Z be the center of Ωs (i.e., the TBP)

and let W ∈ Ωs be any other point. Then
1
kHe (Z) − He (W )k2 < ; ∀e ∈ (0, ∞).
10

Proof: Let H = He . Say that a special path in Ω is a 7-segment polygonal

path γ such that the ith segment is parallel to the ith coordinate direction.
We can connect Z = Z0 to W = Z7 by a special path γ, all of whose
segments have length at most 2−12 . Let Z0 , ..., Z7 be the vertices of γ. Let
∆ij = H(Zi) − H(Zj ). Integrating along the kth segment, we get the bound
∆k,k−1 ≺ 2−12 Ψk . Hence
P7
k=1 Ψk
7
X
2 345
kH(Z) − H(W )k2 = k∆70 k2 = ∆k,k−1 ≤ < .
k=1
2 212 4096
(162)
This last quantity is less than 1/10. ♠

Proof of Lemma 2.4: Let e = 1 or e = 2. Let H = He . Let W ∈ Ωs be

some point. We want to show that H(W ) is positive definite. Let A = H(Z)
and B = H(W ). By the preceding Corollary, we have B = A + ∆ with
k∆k2 < 1/10.
Now let v ∈ R7 be any unit vector. We have
1
Bv · v = Av · v + ∆v · v ≥ Av · v − |∆v · v| > − k∆k2 > 0. (163)
10
Hence B is positive definite. The fact that Av · v > 1/10 comes from the fact
that the lowest eigenvalue of the symmetric matrix A is greater than 1/10. ♠

55
9.2 A Stronger Bound
Consider the function
φ(x) = xe/2 (164)
and the interval  
1 1 −9
I= + +2 . (165)
4 2
Define
dk φ
ck (e) = sup (x) (166)
x∈I dxk
Setting ck = ck (e), we define
q
Υ(e) = 19336c21 + 19036c1c2 + 4922c22 + 1474c1 c3 + 772c2 c3 + 31c23 . (167)

Our next result uses the notation from the Variation Bound.

Lemma 9.5 (Variation Bound II) Ψ(e) < Υ(e) for all e ∈ (0, ∞).

We find easily that

(c1 (1), c2 (1), c3 (1)) = (1, 2, 12). (168)

Plugging this into Equation 167, we get Ψ(1) < 345. Similarly, we have

(c1 (2), c2 (2), c3 (2)) = (1, 0, 0). (169)

This yields Ψ(2) < 140.

Some elementary calculus shows that

sup c1 (e) < 1.1; sup c2 (e) < 3.2; sup c3 (e) < 16. (170)
e∈(0,∞) e e∈(0,∞)

We omit the details, because we don’t use the general bounds anywhere in
our main proof. When we plug in these bounds we find that Ψ(e) < 463 for
all e. Hence, the Variation Bound II implies the Variation Bound.

Remark: The Variation Bound II is generally much better than the Varia-
tion Bound. For instance, Υ(e) → 0 as e → 0 or as e → ∞.

56
9.3 Proof of the Variation Lemma II
Our proof usually suppresses the dependence on the exponent e.
A point (z0 , z1 , z2 , z3 ) ∈ Ωs is such that each zm lies within a square
∆m+1 of side length 2−11 about one of the points of the TBP configuration,
normalized as in Equation 30. Let R1 be the union of these 4 squares,
∆1 , ..., ∆4 . Let R2 ⊂ C 2 denote the set of points z1 , z2 which arise as a
disjoint pair of finite points of a configuration of Ωs . Here R2 consists of
12 components, all of the form ∆i × ∆j . The components ∆i × ∆j and
∆j × ∆i are partners. We let ∆5 , ..., ∆10 be 6 components, no two of which
are partners.
Recall that Σ−1 is inverse stereographic projection. See Equation 5. Set-
ting z = x + iy, define
1 1 + x2 + y 2
F (z) = = . (171)
kΣ−1 (z) − Σ−1 (∞)k2 4

1 (1 + x21 + y12)(1 + x22 + y 2 )

G(z1 , z2 ) = = . (172)
kΣ−1 (z1 ) − Σ−1 (z2 )k2 4(x1 − x2 )2 + (y1 − y2 )2
Let φ be as in Equation 164. According as m ≤ 4 or m ≥ 5, define
bm = φ ◦ Um ;
U Um = F or G (173)

Define

Φ(i, j, k, m) = sup |Di Dj Dk Um |; Φ(i, bm |. (174)

b j, k, m) = sup |Di Dj Dk U
∆m ∆m

Here i, j, k ∈ {1, ..., 7} and m ∈ {1, ..., 10} and Dk is the kth partial deriva-
tive.
We have
10
X
E(z0 , z1 , z2 , z3 ) = bm ,
U (175)
m=1

where the arguments of the functions on the right hand side are suitably
chosen sub-lists of (z0 , z1 , z2 , z3 ). Therefore
10
X
|Dk Di Dj (E)| ≤ Ψk (i, j) := b j, k, m).
Φ(i, (176)
m=1

57
With this definition of Ψk , we have Dk H ≺ Ψk . Summing over k, we have
7
X X
Ψk ≤ b j, k, m) .
Φ(i, (177)
2 2
k=1 k,m

Say that i ∈ {1, ..., 7} is related to an index m ∈ {1, ...., 10} if a change
in the coordinate xi moves one of the points in the argument of U bm . For
instance i = 2 is related to m = 5 because changing x2 changes the location
of the point z1 = x1 + ix2 , and the argument of U5 is the point (z0 , z1 ) ∈ ∆5 .
b j, k, m) = 0 unless i, j, k are all related to m.
Clearly Φ(i,
When all indices are related to m, we use the chain rule
bm = φ′ × Di Dj Dk Um + φ′′′ × Di Um × Dj Um × Dk Um +
Di D j Dk U

φ′′ ×Di Dj Um ×Dk Um +φ′′ ×Dj Dk Um ×Di Um +φ′′ ×Dk Di Um ×Dj Um . (178)
Given the definitions of the regions R1 and R2 we have Um (∆m ) ⊂ I,
where I is the interval from Equation 165. Combining this fact with Equation
166 and Equation 178, we have

b j, k, m) ≤ c1 Φ(i, j, k, m) + c3 Φ(i, m)Φ(j, m)Φ(k, m)+

Φ(i,
c2 Φ(i, j, m)Φ(k, m) + c2 Φ(j, k, m)Φ(i, m) + c2 Φ(k, i, m)Φ(j, m). (179)
Let δm be the center of ∆m . Define

a(i, j, k, m) = |Di Dj Dk Um (δm )|. (180)

We define a(i, j; m) and a(i; m) similarly. Define

 1   1 1 1
v = ((v(i), v(i, j), v(i, j, k)) = , 0, 0 or , , , (181)
1000 200 50 10
according as m ≤ 4 or m ≥ 5. Let b(i; m) = a(i; m) + v(i), etc.

Lemma 9.6 for all (i, j, k, m) we have

Φ(i, m) ≤ b(i, m); Φ(i, j, m) ≤ b(i, j, m); Φ(i, j, k, m) ≤ b(i, j, k, m).

58
 Combining Equation 177, Equation 179, and Lemma 9.6, we have
 
c1 b(i, j, k, m)+
 
 
 c2 (b(i, j, m)b(k, m)+ 
 
X X  
 
Ψk ≺ δ(i, j, k, m) 
 c 2 (b(j, k, m)b(i, m)+ 
 (182)
k k,m  
 
 c2 (b(k, i, m)b(j, m)+ 
 
 
c3 b(i, m)b(j, m)b(k, m)

Here δ(i, j, k, m) = 1 if i, j, k are all related to m, and otherwise 0. When we

compute the norm of the right hand side in Mathematica, we get the square
root of a polynomial in c1 , c2 , c3 . Rounding the coefficients of this polynomial
up to integers, we get Υ. This completes the proof of the Variation Bound
II, modulo the proof of Lemma 9.6.

9.4 Proof of Lemma 9.6

When m ≤ 4, we are dealing with the function F from Equation 171. We
have Di F = xi /2 and all higher derivatives are constant. Lemma 9.6 in this
trivial case now follows from the fact that radius(∆m ) < 1/500.
For the nontrivial cases, we fix some value of m ∈ {5, ..., 10} and consider
the cube ∆ = ∆m ⊂ C 2 . This cube has side-length 2−11 . The center point
is δ = δm . Define

ak = |Dk G(δ)|; Φk = max |Dk G(z1 , z2 )|. (183)

(z1 ,z2 )∈∆

The maxima here are taken over all kth partial derivatives Dk . An exact,
finite calculation, done for each of the 6 choices of m, yields

a1 ≤ 1; a2 ≤ 4; a3 ≤ 18; a4 ≤ 96; a5 ≤ 600. (184)

Now we give an estimate on Φ6 . Given a polynomial P ∈ R[x1 , y1 , x2 , y2],

let |P | denote the sum of the absolute values of the coefficients of P . We
have the easy upper bound

|P (z1 , z2 )| ≤ |P | max(|z1 |, |z2 |)d ; d = degree(P ). (185)

59
Here we have set zj = xj + iyj . Calculating symbolically we find that each
6th partial derivative D6 has the following structure. There is a polynomial
P , depending on the choice of derivative, such that
P
D6 G = ; |P | ≤ 4519440; deg(P ) ≤ 10. (186)
kz1 − z2 k7

We (easily) have

max(kz1 k, kz2 k, kz1 − z2 k−1 ) < 1 + 2−8 ; ∀(z1 , z2 ) ∈ ∆. (187)

Hence
Φ6 ≤ 4519440 × (1 + 2−8 )17 < 5000000. (188)

Lemma 9.7 Φk < ak + 2−10 × Φk+1 .

Proof: Let (w1 , w2 ) be a point of ∆. We can connect (w1 , w2 ) to a point

(z1 , z2 ) ∈ δ2 by a 4-segment polygonal path γ such that the ith segment has
length at most 2−12 and moves in the ith coordinate direction. This lemma
now follows from integration. ♠

Now we apply Lemma 9.7 in an iterative way.

Φ5 ≤ 600 + 5000000 × 2−10 < 5483 (189)

Φ4 ≤ 96 + 5483 × 2−10 < 102. (190)

1
|Φ3 − a3 | < 102 × 2−10 < ; Φ3 < 18.1 (191)
10
1
|Φ2 − a2 | < 18.1 × 2−10 < ; Φ2 < 4.02. (192)
50
1
|Φ1 − a1 | ≤ 4.02 × 2−10 < . (193)
200
The R.H.S. of Equation 191 comes from the L.H.S. and Equation 184. Like-
wise, the R.H.S. of Equation 192 comes from the L.H.S. and Equation 184.
This completes the proof of Lemma 9.6.

60
10 Computational Issues
10.1 General Remarks
Our calculations in this paper are of two kinds. The material in §9 makes
some exact calculations in Mathematica [W] and the rest of the paper makes
calculations in Java. For the Mathematica calculations, we need to take
derivatives of rational
√ functions or their square roots, and evaluate them
at elements of Q[ 3]. We also need to simplify and group terms of some
polynomials. Everything is manipulated exactly. The user who downloads
our Java program will also find our Mathematica files in the same directory.
The bulk of the calculations are done in Java, and these are what we
discuss below. We take the IEEE-754 standards for binary floating point
arithmetic [I] as our reference for the Java calculations. This 1985 document
has recently been superceded by a 2008 publication [I2]. We will stick to
the 1985 publication for three reasons. First, the 1985 version is shorter and
simpler. Second, the portion relevant to our computation has not changed in
any significant way. Third, we think that some of the computers running our
code will have been manufactured between 1985 and 2008, thus conforming
to the older standard.

10.2 Doubles
With a view towards explaining interval arithmetic, we first describe the way
that Java represents real numbers. Our Java code represents real numbers
by doubles, in a way that is an insignificant modification of the scheme dis-
cussed in [I, §3.2.2]. To see what our program does, read the documentation
for the longBitsToDouble method in the Double class, on the website
[Link] According to this documentation –
and experiments verify that it works this way on our computer – our program
represents a double by a 64 bit binary string, where

• The first bit is called s.

• The next 11 bits are a binary expansion of an integer e.

• If e = 0, the last 52 bits are the binary expansion of an integer m.

• If e 6= 0, the last 52 bits are binary expansion of the number m − 252 .

61
The real number represented by the double is

(−1)s × 2e−1075 × m. (194)

Example: The double representing −317 is stored as the 64 bit string

1/10000000111/001111010 . . . 0.

The slashes are put in to emphasize the breaks. Here s = 1 and (since e 6= 0)

e = 210 + 8 + 4 + 1 = 1031.

There are 44 zeros at the end of the word, and

m = 244 (0 + 0 + 32 + 16 + 8 + 4 + 0 + 1) + 252 = 244 (317).

Hence
(−1)s × 2e−1075 × m = −317.
Now we come to the main point of our discussion above. The non-negative
doubles have a lexicographic ordering, and this ordering coincides with the
usual ordering of the real numbers they represent. The lexicographic ordering
for the non-positive doubles is the reverse of the usual ordering of the real
numbers they represent. To increment x+ of a positive double x is the very
next double in the ordering. This amounts to treating the last 63 bits of
the string as an integer (written in binary) and adding 1 to it. With this
interpretation, we have x+ = x + 1. We also have the decrement x− = x − 1.
Similar operations are defined on the non-positive doubles. These operations
are not defined on the largest and smallest doubles, but our program never
encounters (or comes anywhere near) these.

10.3 The Basic Operations

Let D be the set of all doubles. Let

R0 = {x ∈ R| |x| ≤ 230 } (195)

Our choice of 230 is an arbitrary but convenient cutoff. Let D 0 denote the
set of doubles representing reals in R0 .
According to the discussion in [I, 3.2.2, 4.1, 5.6], there is a map R0 → D0
which maps each x ∈ R0 to some [x] ∈ D0 which is closest to x. In case

62
there are several equally close choices, the computer chooses one according
to the method in [I, §4.1]. This “nearest point projection” exists on a subset
of R that is much larger than R0 , but we only need to consider R0 . We also
have the inclusion r : D0 → R0 , which maps a double to the real that it
represents.
Our calculations use the 5 functions
√
+; −; ×; ÷; . (196)
These operations act on R0 in the usual way. Operations with the same name
act on D 0 . Regarding the first 5 basic operations, [I, §5] states that each of
the operations shall be performed as if it first produced an intermediate result
correct to infinite precision and with unbounded range, and then coerced this
intermediate result to fit into the destination’s format. Thus, for doubles x
and y.
√ hp i
x= r(|x|) ; x ∗ y = [r(x) ∗ r(y)]; ∗ ∈ {+, −, ×, ÷}. (197)

The operations on the left hand side represent operations on doubles and the
operations on the right hand side represent operations on reals.

Remark: Exceptions to Equation 197 can arise if we divide by 0 or a number

too close to zero. This will produce an overflow error – too large a number
to be accurately represented by a double. To see that this never happens, we
check that we never divide by a number smaller than 2−11 . Moreover, as we
discuss at the end of this chapter, we never take the square root of a number
less than 2−5 .

10.4 Interval Arithmetic

An interval is a pair I = (x, y) of doubles with x ≤ y. Say that I bounds
z ∈ R0 if x ≤ [z] ≤ y. This is true if and only if x ≤ z ≤ y. Define
[x, y]o = [x− , y+ ]. (198)
This operationis well defined for doubles in D 0 . We are essentially rounding
out the endpoints of the interval. Let I0 and I1 denote the left and right
endpoints of I. Letting I and J be intervals, we define
√ √ √
I ∗ J = (min Ii ∗ Ij , max Ii ∗ Ij )o ; I = ( I 0 , I 1 )0 (199)
ij ij

63
That is, we perform the operations on all the endpoints, order the results,
and then round outward. Given Equation 197, we have the following facts.

1. If I bounds x and J bounds y then I ∗ J bounds x ∗ y.

√ √
2. If I > 0 bounds x > 0 then I bounds x.

After we have defined intervals and their basic operations, we define other
Java objects based on intervals. Namely, interval versions of complex num-
bers, vectors in R3 , dyadic segments and squares, and dyadic boxes. For
instance, the interval version of a complex number is an object of the form
X + iY , where X and Y are intervals. The algebra of these interval objects
– e.g. complex addition or the dot product – is formally the same as the
corresponding algebra on the usual objects. At every step of the calculation,
the real version of the object is bounded, component by component, by the
interval version. If some particular interval object passes our test, it means
that all the real objects bounded by it also pass.

10.5 Implementation Details

We implement the interval arithmetic in a way that tries to minimize the
time we spend using it. We run the floating point algorithm until we notice
that a box has passed one of the floating point tests. Then, we re-perform
the interval arithmetic test on the interval arithmetic version of the box.

• If the box passes the interval arithmetic test, we eliminate it, and switch
back to the floating point algorithm.

• If the box fails the interval arithmetic test, we just act as if it failed
the corresponding floating point test, and resume the algorithm. We
call this case a mismatch. It is harmless.

Even though the mismatches are harmless, we prefer to have few mis-
matches, so that run time for the interval arithmetic version of the code is
easier to predict from the run time of the floating point version. There are
4 kinds of potential mismatches, corresponding to the 4 kinds of tests we
perform, as discussed in §2.5. The two mismatches that actually arise (or,
rather arose) in abundance are the Energy Estimator mismatches and the
Redundancy Eliminator mismatches. Once we make our modifications, the

64
code runs completely without mismatches for the 1/r potential and with only
a 4 mismatches for the 1/r 2 potential.
To prevent the Energy Eliminator mismatches, we arrange the floating
point eliminator so that a dyadic box passes only if the minimum energy is
above E0 +2−40 , where E0 is the energy of the T.B.P. In this way, the interval
versions of the boxes that are actually checked have a bit of a cushion. The
reason why this fudge factors does not cause our code to halt is that the
inequality E(X) − E0 < 2−40 only occurs well inside the L∞ neighborhood of
size 2−14 about the TBP.
Adding the fudge factor of 2−40 makes the floating point test harder to
pass, but does not change the validity of the program. Logically, the float-
ing point calculation simply finds a candidate partition of the configuration
space, in which each box in the partition passes one of our tests. Adding the
fudge factor does change the final partition a bit, but it doesn’t destroy the
basic property of the partition.
Mismatches occur for the Redundancy Eliminator because we sometimes
eliminate configurations where there is an exact equality of coordinates. This
equality will fail for the corresponding intervals, because of a tiny overlap.
To get around this problem, we associate to each dyadic object a Gaussian
integer that represents 225 z, where z is the center of the object. When we
subdivide a dyadic object, we perform the arithmetic on the Gaussian integer
exactly, so as to compute the exact value of the centers of the subdivided
objects. We never subdivide more than 24 times – and in fact the maximum
is about 17 – so we never arrive at a situation where 225 z is not an integer.
When it comes time to compare the various coordinates of a dyadic object,
we actually compare 225 times those coordinates, so that we are working en-
tirely with integers. This lets us make exact comparisons.

There is one more fine point√we would like to mention. We want to avoid
taking expressions of the form√ I, where I is an interval that is too close to
0. If 0 ∈ I, then I0 < 0 and I 0 causes an arithmetic error. The only time
this issue comes up is in the bounds from Lemma 3.1. The quantity 4 − D 2
might be close to 0 if D is close to 2. To avoid this problem, we define a new
function√σ, which has the property σ(I) = 2−5 if I1 < 2−10 and otherwise
σ(I) = I. (It never happens that I1 > 2−10 and I0 < 0.) When computing
√
the interval version of the bound in Lemma 3.1, we use σ in place of .
This causes no problem with the proof, because the bound is not as strong
√
when we use σ in place of .

65
11 References
[BBCGKS] Brandon Ballinger, Grigoriy Blekherman, Henry Cohn, Noah
Giansiracusa, Elizabeth Kelly, Achill Schurmann,
Experimental Study of Energy-Minimizing Point Configurations on Spheres,
arXiv: math/0611451v3, 7 Oct 2008

[C] Harvey Cohn, Stability Configurations of Electrons on a Sphere, Math-

ematical Tables and Other Aids to Computation, Vol 10, No 55, July 1956,
pp 117-120.

[CK] Henry Cohn and Abhinav Kumar, Universally Optimal Distributions

of Points on Spheres, J.A.M.S. 20 (2007) 99-147

[CCD] online website:

[Link] wales/CCD/Thomson/[Link]

[DLT] P. D. Dragnev, D. A. Legg, and D. W. Townsend, Discrete Loga-

rithmic Energy on the Sphere, Pacific Journal of Mathematics, Volume 207,
Number 2 (2002) pp 345–357

[HS], Xiaorong Hou and Junwei Shao, Spherical Distribution of 5 Points

with Maximal Distance Sum, arXiv:0906.0937v1 [[Link]] 4 Jun 2009

[I] IEEE Standard for Binary Floating-Point Arithmetic (IEEE Std 754-
1985) Institute of Electrical and Electronics Engineers, July 26, 1985

[I2] IEEE Standard for Floating-Point Arithmetic (IEEE Std 754-2008) In-
stitute of Electrical and Electronics Engineers, August 29, 2008.

[RSZ] E. A. Rakhmanoff, E. B. Saff, and Y. M. Zhou, Electrons on the

Sphere, Computational Methods and Function Theory, R. M. Ali, St. Ruscheweyh,
and E. B. Saff, Eds. (1995) pp 111-127

[SK] E. B. Saff and A. B. J. Kuijlaars, Distributing many points on a Sphere,

Math. Intelligencer, Volume 19, Number 1, December 1997 pp 5-11

66
[T] J. J. Thomson, On the Structure of the Atom: an Investigation of the
Stability of the 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. Philosophical magazine, Series 6,
Volume 7, Number 39, pp 237-265, March 1904.

[W] S. Wolfram, The Mathematica Book , 4th ed. Wolfram Media/Cambridge

University Press, Champaign/Cambridge (1999)

[Wa] D. Watkins, Fundamentals of Matrix Calculations, John Wiley and

Sons (2002)

[Y], V. A. Yudin, Minimum potential energy of a point system of charges

(Russian) Diskret. Mat. 4 (1992), 115-121, translation in Discrete Math
Appl. 3 (1993) 75-81

67