
DYNAMICAL SYSTEMS IN PURE MATHEMATICS

A Dissertation

Presented to the Faculty of the Graduate School

of Cornell University

in Partial Fulfillment of the Requirements for the Degree of

Doctor of Philosophy

by

Max Lipton

May 2023


© 2023 Max Lipton

ALL RIGHTS RESERVED


DYNAMICAL SYSTEMS IN PURE MATHEMATICS

Max Lipton, Ph.D.

Cornell University 2023

The author’s research program involves several topics which differ at first glance.

However, they all share the common theme of exploring how geometry and topology

influences dynamical equilibria. The dissertation is broken into three parts: the hyper-

bolic geometry of higher-dimensional Kuramoto oscillators, electrostatic knot theory,

and minimal surfaces with Möbius energy on the boundary. Each part is further di-

vided into chapters adapting the author’s preprints and published papers, which have

appeared in journals in physics, applied mathematics, and pure mathematics.


BIOGRAPHICAL SKETCH

Max Lipton was born in Salem, Oregon on June 18, 1994 as the only child to Stephen

and Mythuan “Mitoo” Lipton. His father was a public criminal defense attorney, and

his mother was a refugee of the Vietnam War who later became a certified medical

technician. After graduating from South Salem High School, Max attended Willamette

University and earned a B.A., double majoring in mathematics and computer science.

During his junior year, Max spent a semester abroad in Russia, where he participated

in the Math in Moscow exchange program. Despite having visited every state in the

US as a teenager, except for Alaska and Hawaii, this trip was the first time he left

North America.

He is about to graduate with a Ph.D. in mathematics from Cornell University, ad-

vised by Steve Strogatz. After a “wrestling match” unlike anything his advisor has

ever seen in his thirty years of advising, Max was awarded an NSF Mathematical Sci-

ences Postdoctoral Research Fellowship at the Massachusetts Institute of Technology.

He is set to marry his beloved Courtney Elizabeth on October 21, 2023.

iii


To Courtney.

iv


ACKNOWLEDGEMENTS

There are so many people who have contributed to my mathematical development.

Everywhere I have been, I had the pleasure of meeting many excellent teachers, men-

tors, and friends who helped me in big ways and small. I owe all of them my deepest

gratitude. In particular, without the opportunity and support provided by the late

Bob Strichartz and his Analysis on Fractals REU, I would never have made it to Cor-

nell. I’m also incredibly grateful for the support of my thesis advisor Steve Strogatz,

who I am proud to also call my friend. I would also like to give special thanks to the

following people who, through one way or another, left a lasting impact on my in-

tellectual and personal development: Wayne Atwood, Warren Trotter, Inga Johnson,

Colin Starr, Fritz Ruehr, Charleen Gust, Gill Grindstaff, Aaron Calderon, Max Hall-

gren, Kelly Delp, Alex Townsend, Rennie Mirollo, Gokul Nair, and Xin Zhou. The

secret to life is enjoying the company of those you meet along the way.

I would also be remiss if I didn’t acknowledge the monumental support given

from my mom, dad, and Courtney. And last, but not least, I would like to give a very

special shout-out to the cats who have provided me with much love, support, and

mathematical insight: Latitude, Gratitude, Amplitude “Ampi”, Plus, Minus, Equals,

and Infinity-Algebra “Fifi.” (My parents came up with the names). Every meow is

special.

v


TABLE OF CONTENTS

Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi

1 Introduction 1

I The hyperbolic geometry of higher-dimensional Kuramoto os-
cillators 4

2 The Kuramoto model on a sphere 5
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.1 The Kuramoto Model on a Sphere . . . . . . . . . . . . . . . . . . 9
2.3 Reduction of system parameters . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Hyperbolic geometry and Möbius transformations . . . . . . . . . . . . . 13
2.5 Infinitesimal generators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.6 Comparison of Z versus W coordinates . . . . . . . . . . . . . . . . . . . 17
2.7 Continuum limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.8 Complex case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.9 Relation to previous research . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.10 An example: First-order linear order parameter gives gradient system . 26

2.10.1 Computer visualization . . . . . . . . . . . . . . . . . . . . . . . . 27
2.10.2 Existence of hyperbolic gradient . . . . . . . . . . . . . . . . . . . 29
2.10.3 Analysis of dynamics . . . . . . . . . . . . . . . . . . . . . . . . . 31

II Electrostatic knot theory 39

3 A lower bound on the critical points of the electric potential of a knot 40
3.1 Preliminaries and problem formulation . . . . . . . . . . . . . . . . . . . 40
3.2 Preliminary definitions and lemmas . . . . . . . . . . . . . . . . . . . . . 41
3.3 Proof of the tunnel number bound . . . . . . . . . . . . . . . . . . . . . . 45

3.3.1 Construction of the tunneling . . . . . . . . . . . . . . . . . . . . . 45
3.3.2 Constructing a smooth boundary around the tunneling . . . . . . 48
3.3.3 The deformation retraction to a handlebody . . . . . . . . . . . . 51

4 Further topological approaches to knotted electric charge distributions 53
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2 Bifurcations through degenerate critical points . . . . . . . . . . . . . . . 58
4.3 A partial sharpening of the tunnel number bound . . . . . . . . . . . . . 61

vi


4.4 The relationship between the critical set and the Morse code of equipo-
tential surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.5 Concluding Remarks and Future Directions . . . . . . . . . . . . . . . . . 70

5 Specific charge distributions and perspectives in physics 71
5.1 Motivations from physics . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.2 Charged loops in the plane: Specific examples . . . . . . . . . . . . . . . 72

5.2.1 Rectangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.2.2 Stadium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.2.3 Ellipse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5.3 Charged Knots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.3.2 Physical intuition about equipotential surfaces and equilibrium

points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.3 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.3.4 Results and conjectures for charged knots . . . . . . . . . . . . . 84

5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

III Plateau problems with Möbius energy on the boundary 90

6 Minimizers of the Möbius-Plateau energy 91
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.3 Energy Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

6.3.1 Lower Semicontinuity of EP . . . . . . . . . . . . . . . . . . . . . . 96
6.3.2 Equicontinuity at the Boundary . . . . . . . . . . . . . . . . . . . . 97
6.3.3 Convergence in the Interior . . . . . . . . . . . . . . . . . . . . . . 100

6.4 Möbius-Plateau Stationary Helix Pairs . . . . . . . . . . . . . . . . . . . . 101

7 Stationary helix pairs and complex asymptotics of the Möbius gradient 115
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7.2 Möbius-stationary symmetric double helices . . . . . . . . . . . . . . . . 118
7.3 Möbius-Plateau stationary symmetric screws . . . . . . . . . . . . . . . . 128

Bibliography 132

vii


CHAPTER 1

INTRODUCTION

When everything is said and done

Still looking for answers, if only one

Turn my back, the urge has gone

Left with no reason, we come undone
Rob Swire, The Island

Dynamical systems are differential equations. Differential equations are vector

fields. And vector fields are topology. This perspective is mathematically reduc-

tive, but it is an essential scientific heuristic for building bridges between disciplines

whose respective research communities are generally disjoint. My dissertation strad-

dles several different fields in pure and applied mathematics, but the common thread

is exploring how topology and geometry influence dynamical equilibria. Employing

techniques from areas including knot theory, calculus of variations, and dynamical

systems, I seek to build new bridges between disciplines. I also implement numerical

simulation and computer visualization to assist in the formulation of new conjectures

which are then proven rigorously. Each part and chapter stands independently, and is

adapted from the preprints and publications I’ve contributed to throughout my grad-

uate career. I have obtained what I hope are considered to be significant milestone

results, but by no means are any of these lines of inquiry complete. I believe all three

topics are manifestations of deep mathematical phenomena yet to be unraveled.

The first part is about the hyperbolic geometry of higher-dimensional Kuramoto

oscillators. Kuramoto dynamics describe systems of coupled particles on S 1 which

show long-term synchronization behavior under generic initial conditions. The re-

sults of Watanabe-Strogatz yield constants of motion which greatly reduce the sys-

1


tem’s dimension [128]. Later work by Marvel-Mirollo-Strogatz, Chen-Engelbrecht-

Mirollo, and Lohe [20, 74, 80] show these constants of motion have an interpretation

in terms of hyperbolic geometry. A paper I co-authored with Mirollo and Strogatz

deals with higher-dimensional oscillators living on S d rather than S 1 [69]. We used

the generalized Watanabe-Strogatz transform to reduce the dynamics to the Lie group

of hyperbolic isometries of the sphere’s Poincaré ball interior, whilst giving the condi-

tions for a linearly coupled Kuramoto system to synchronize.

The second part concerns electrostatic knot theory, which investigates the poten-

tial induced by a knotted curve with uniform electric charge distributed on it. The

overarching objective is to find relationships between classical knot invariants and the

dynamical properties of the electric field on the complement. One particular result

I proved so far is a lower bound on the potential’s critical set based on the tunnel

number [66] using methods from Morse theory and geometric topology. A recent pa-

per I co-authored with Strogatz and Townsend expands this result by analyzing some

explicitly given knot parametrizations [71].

The third and final part discusses Plateau problems with Möbius boundary en-

ergy, a class of problems in the calculus of variations seeking to minimize novel con-

vex combinations of two well-known energies associated with a curve: its Möbius

energy and the area of its minimally spanning surface. A seminal paper of Freedman-

He-Wang proves several nice properties of the Möbius energy, including regularity

of minimizers and invariance under Möbius transformations [34]. In a paper co-

authored with Nair [70], we proved the existence of minimizers of our novel Möbius-

Plateau energy. In future work, Nair and I hope to use this result to deduce boundary

self-avoidance in other geometric optimization settings, particularly the Euler-Plateau

problems from elasticity theory. To add concrete examples to our theory, I searched

2


for critical helicoidal strips and proved a characterization of these strips via hard esti-

mates of a certain oscillatory integral in the Euler-Lagrange equation [70] in addition

to complex asymptotic methods [67].

3


Part I

The hyperbolic geometry of

higher-dimensional Kuramoto

oscillators

4


CHAPTER 2

THE KURAMOTO MODEL ON A SPHERE

This chapter originally appeared as “The Kuramoto model on a sphere: Explaining

its low-dimensional dynamics with group theory and hyperbolic geometry” by

Lipton, Mirollo, and Strogatz in Chaos Vol. 31, Iss. 9 [69].

2.1 Introduction

In 1975, Kuramoto introduced a model for a large population of coupled oscilla-

tors with randomly distributed natural frequencies. [60] Kuramoto’s model displayed

many remarkable features: It was exactly solvable (at least in some sense), despite

being nonlinear and infinite-dimensional. [61] Its solution shed analytical light on a

phase transition to mutual synchronization that Winfree had previously discovered in

a similar but less convenient system of oscillators. [132, 133] Since then, the Kuramoto

model has been an object of fascination for nonlinear dynamicists, as well as a sim-

plified model for many real-world instances of coupled oscillators in physics, biology,

chemistry, and engineering. [120, 102, 119, 1, 30, 101, 108, 10]

From a mathematical standpoint, one of the most intriguing problems has been

to explain the tractability of the Kuramoto model. What symmetry or other hidden

structure accounts for its solvability?

The first clues came from work on an adjacent topic: series arrays of N identi-

cal overdamped Josephson junctions. The governing equations for those supercon-

ducting oscillators are closely related to the equations of the Kuramoto model [130,

131], and themselves displayed remarkable dynamical features, such as surprisingly

low-dimensional invariant tori[125, 122] and ubiquitous neutral stability of splay

5


states[87], despite the presence of damping and driving in the governing equations.

These features were explained in 1993 by the discovery of a certain change of variables,

now called the Watanabe-Strogatz transformation [127, 128], which showed that the

governing equations have N − 3 constants of motion for all N ≥ 3. Goebel [38] then

pointed out that the Watanabe-Strogatz transformation could be viewed as a time-

dependent version of a linear fractional transformation, a standard tool in complex

analysis. For more than a decade, however, these results did not attract much atten-

tion, perhaps because they were assumed to be restricted to problems about Josephson

junctions, and within that specialized setting, even further restricted to junctions that

were strictly identical.

A breakthrough occurred in 2008 with the work of Ott and Antonsen. [96, 97]

They found an astonishing way to capture the macroscopic dynamics of the infinite-N

Kuramoto model, even when the oscillators’ frequencies were non-identical and ran-

domly distributed. First, they wrote down an ansatz — seemingly pulled out of thin

air — for the density function ρ(θ, ω, t) of oscillators having phase θ and intrinsic fre-

quency ω at time t. Their ansatz had the form of a time-dependent Poisson density (a

density better known for its role in the study of partial differential equations, specifi-

cally for the solution of Laplace’s equation on a disk, given the values of the unknown

function on the bounding circle). By making this ansatz of a Poisson density, Ott

and Antonsen reduced the infinite-N Kuramoto model, an integro-partial differential

equation, to an infinite set of coupled ordinary differential equations. Then, by further

assuming that the intrinsic frequencies of the oscillators were randomly distributed

according to a Lorentzian (Cauchy) distribution, Ott and Antonsen showed that the

order parameter dynamics of the Kuramoto model could be reduced tremendously,

all the way down to an ordinary differential equation for a single scalar variable, the

amplitude of the order parameter. [96] With this discovery, the floodgates were now

6


open. Almost immediately the Ott-Antonsen ansatz was used to solve many long-

standing problems about the Kuramoto model and its variants, as well as to generate

and solve many new problems. [101]

Still, a lot of old questions hung in the air. Both the Watanabe-Strogatz transfor-

mation and the Ott-Antonsen ansatz appeared somewhat unmotivated and almost

miraculous. Where did they come from, and why did they work? It also was not

clear whether they were connected or perhaps even equivalent. There were reasons to

doubt they were linked: the Watanabe-Strogatz transformation could be used for any

finite N ≥ 3, but seemed restricted to identical oscillators, whereas the Ott-Antonsen

ansatz allowed for non-identical oscillators but seemed restricted to the continuum

limit of infinite N. Also, why were linear fractional transformations and Poisson den-

sities — tools from other branches of mathematics — popping up in these studies of

dynamical systems?

Later work made sense of all of this. The Josephson arrays and the Kuramoto

model both turned out to have deep mathematical ties to group theory, hyperbolic ge-

ometry, and projective geometry, and both the Watanabe-Strogatz transformation and

the Ott-Antonsen ansatz were tapping into these structures. [101, 100, 80, 117, 20, 21]

For the Josephson arrays, the governing equations turned out to be generated by a

group action, specifically the action of the Möbius group of linear fractional trans-

formations of the unit disk to itself. Seen in this light, the constants of motion for the

Josephson arrays were cross-ratios, and the invariant tori were group orbits. The same

group-theoretic structure was found to underlie the Kuramoto model (in the special

case where all the oscillator frequencies are identical) as well as other sinusoidally

coupled systems of identical phase oscillators.[80, 19]

In the past few years, several researchers wondered how far this story could be

7


pushed. Are there quantum or higher-dimensional extensions of the Kuramoto model

that might show similar reducibility? A number of results along these lines have now

been found. [72, 73, 123, 23, 24, 42, 43, 74, 52, 18, 17, 75, 28, 44, 76, 41, 53, 27] In par-

ticular, several researchers have explored a generalization of the Kuramoto model in

which the oscillators move on spheres instead of the unit circle; the spheres could be

either the ordinary two-dimensional sphere or higher-dimensional spheres. In par-

ticular, a system of particles moving on the two-sphere has been used to model the

orientation dynamics of swarms of drones flying around in three dimensions. [93] As

we shall demonstrate, these higher-dimensional oscillator models exhibit a dynami-

cal reduction generalizing the reduction for the Kuramoto model. Consequently, the

computational effort required to simulate these systems can be substantially reduced.

A counterpart of the Ott-Antonsen ansatz has also been discovered for the continuum

version of the Kuramoto model for identical oscillators on the d-dimensional sphere,

and used to reduce its infinite-dimensional dynamics to a lower dimensional set of

ordinary differential equations (ODEs). [17] But as before, some of the results appear

disconnected and a bit miraculous.

Our goal in this paper is to show that hyperbolic geometry and group theory can

unify and clarify our understanding of the Kuramoto model on a sphere and make

all the latest results seem natural, just as they did before for the traditional Kuramoto

model. We focus exclusively on the case of identical oscillators, though we expect

that our methods will extend to systems with multiple populations of identical oscil-

lators, or a continuum of oscillator families with a distribution of natural frequencies.

Our approach explains the model’s reducibility for any finite number of oscillators,

as well as for the continuum limit, and it reveals why Poisson densities arise again in

this setting. There is a close connection to Laplace’s equation and harmonic analysis,

as we will see in Section V. We also find that complex analysis is not really essential,

8


which is just as well, since it does not generalize to the higher-dimensional spheres

being considered here. Instead, the proper mathematical setting is harmonic analysis

and hyperbolic geometry on higher-dimensional balls. Our work also allows us to

go beyond merely unifying existing results. For instance, by establishing that linearly

coupled systems of identical Kuramoto oscillators on a sphere have a hyperbolic gra-

dient structure, we can prove new global stability results about convergence to the

synchronized state, as described in Section VIII.

2.2 Preliminaries

2.2.1 The Kuramoto Model on a Sphere

In a pioneering paper, Lohe [72] observed that there are at least two natural gen-

eralizations of the Kuramoto model to higher dimensions. One of them replaces the

phases θ j of the original Kuramoto model[60, 61] with complex numbers exp(iθ j) on the

unit circle and then views those as equivalent to 2 × 2 rotation matrices parametrized

by a rotation angle θ j. From there, it is a natural step to consider other Lie groups of

matrices, many of which are non-Abelian.

Our concern in this paper, however, is with a different generalization of the Ku-

ramoto model. Instead of regarding oscillators as particles moving on the unit circle,

we think of them as particles moving on the unit sphere. The sphere could be the sur-

face of the ordinary unit ball in three dimensions, or some higher-dimensional sphere

S d−1 in Rd. When d = 2, the sphere reduces to the unit circle in the plane, and the

model reduces to the original Kuramoto model.

The governing equations for the Kuramoto model on a sphere are

ẋi = Aixi + Z − ⟨Z, xi⟩xi, i = 1, . . . ,N, (2.1)

9


where xi is a point on the unit sphere S d−1 ⊂ Rd, each Ai is an antisymmetric d×d matrix,

and Z ∈ Rd is a d-dimensional vector analogous to the complex order parameter for

the classic Kuramoto model. In Eq. (2.1), the matrix Ai and the vector Z are functions of

the configuration (x1, . . . , xN) of points on the sphere. Note that Z does not depend on

i; like the usual Kuramoto order parameter, it plays the role of a mean-field quantity

that couples all the “oscillators” xi together. The antisymmetric matrix Ai is the higher-

dimensional counterpart of an intrinsic frequency ωi in the original Kuramoto model.

A straightforward computation shows that the dot product between an oscillator’s

instantaneous position and instantaneous velocity satisfies ⟨xi, ẋi⟩ = 0, which proves

that oscillators that start on the unit sphere stay on it forever. The state space for this

system is the N-fold product X = (S d−1)N , which has dimension N(d − 1). Later we will

also consider the natural infinite-N analogue of (2.1), where a state is a probability

measure on S d−1.

In what follows, we allow Z to be any smooth function on the state space X, though

in examples we usually restrict to fairly simple functions, like a linear combination of

the form

Z =
N∑

i=1

aixi

where the ai are real constants.

2.3 Reduction of system parameters

There is a general technique for dimensional reduction of systems which we pause

to describe. Suppose we have a smooth manifold X, which we think of as a state space,

and a group G acting on X, where G is also a smooth manifold (in other words, G is a

Lie group). Then the group action induces a space of vector fields on X, the so-called

infinitesimal generators of the action.

10


To construct these generators, let γ(t) be a smooth curve in G with γ(0) = e, the

identity element in G. Then the derivative γ̇(0) = v, where v is a vector in the tangent

space TeG of G at e. This vector v is in turn associated very naturally with a corre-

sponding vector ṽ in the tangent space of X, as follows. For each x ∈ X, t 7→ γ(t)x is a

smooth curve in X, and its derivative at t = 0 defines a vector ṽx in the tangent space at

x. The vector field ṽ = (ṽx) is the infinitesimal generator corresponding to the element

v in the tangent space of G at e. At each point x ∈ X, the infinitesimal generators vx

span a linear subspace Vx of the tangent space TxX, which is exactly the set of vectors

tangent to the group orbit Gx at the element x.

Now suppose we have a vector field ξ on X, which defines a dynamical system

on the state space X. If ξx ∈ Vx at each point x ∈ X, then the flow corresponding to

the vector field ξ will be constrained to lie on the group orbits Gx. If the dimension

of G is less than the state space X, then this will give us a dimensional reduction

of the dynamics from dim X to dim G. Now suppose, as is the case in applications

of this methodology, the correspondence G → Gx is one-to-one for generic x ∈ X;

equivalently, the stabilizer subgroup Gx = {e} for generic x ∈ X. Then for each x ∈ X, the

flow on the group orbit Gx is equivalent to a flow on G, which is a lower-dimensional

space than the state space X.

Here is a familiar example: Consider the orthogonal group G = S O(d) consisting

of orientation-preserving linear isometries of Rd. (For an intuitive picture, think of

these isometries as rotations.) Then G acts on Rd, and the corresponding infinitesimal

generators are the linear vector fields νx = Ax, where A is any skew-symmetric matrix.

We can also let G act on the product space X = (Rd)N of N-tuples x = (xi), xi ∈ Rd, and

then the infinitesimal generators have the form (νx)i = Axi for some skew-symmetric

matrix A. We could also, if we like, restrict X to N-tuples (xi) with xi ∈ S d−1, the state

11


space for (2.1). Now suppose we had a dynamical system on X of the form ẋi = Axi,

where the skew-symmetric matrix A is a function of the configuration x = (xi); this is

just the special case of (2.1) with Z = 0. Then the dynamics on X reduces to dynamics

on G, which has dimension d(d − 1)/2, and for large N this is much smaller than the

dimension of X, which is N(d − 1). Basically the configuration (xi) of points on the

sphere S d−1 can only move collectively by a rotation of the sphere, so the dynamics

reduces to a dynamical system on S O(d).

We want to apply this methodology to the system (2.1). But since that system

generally has Z , 0, we need a different group action to make this strategy work.

Fortunately, vector fields of the form seen in the Kuramoto model,

ẋ = Ax + Z − ⟨Z, x⟩x, (2.2)

turn out to arise as the infinitesimal generators of the group action of a larger group

G acting on the sphere S d−1 and its interior, the unit ball Bd. This larger group is

the Möbius group of isometries of the hyperbolic geometry on Bd. It contains the

orthogonal group as a proper subgroup, but has bigger dimension d(d + 1)/2.

So for the Kuramoto system (2.1), if the matrices Ai all happen to be identical, we

can reduce the dynamics of the system to a much smaller system on this Möbius group

G, and use this reduction to understand the dynamics on the larger state space X. This

is the essence of the approach we take. Ultimately we apply it to prove a synchro-

nization theorem for the system (2.1) for the special order parameter Z =
∑N

i=1 aixi with

ai > 0. But first we need to show that vector fields on S d−1 of the form in (2.2) are

indeed the infinitesimal generators of the action of the larger Möbius group.

12


2.4 Hyperbolic geometry and Möbius transformations

In this paper, a Möbius transformation is a composition of Euclidean isometries and

spherical inversions of Rd mapping the unit ball homeomorphically to itself and pre-

serving orientation. This is a more restrictive definition than the commonly defined

Möbius transformations which in general do not need to preserve the unit ball.

As in the case d = 2, flows of the form (1) are intimately related to the natural

hyperbolic geometry on the unit ball Bd with boundary S d−1. This geometry has metric

ds =
2|dx|

1 − |x|2
,

where |dx| is the ordinary Euclidean metric. Isometries are assumed to be with respect

to this hyperbolic geometry, unless otherwise qualified as Euclidean. The metric ds has

constant (sectional) curvature equal to −1, and we can describe its isometries, which

generalize the Möbius transformations preserving the unit disc for d = 2. For d = 2,

let w ∈ B2 and consider the Möbius transformation

Mw(x) =
x − w

1 − wx
,

which preserves the unit disc B2 and its boundary S 1.

2.5 Infinitesimal generators

Having parametrized the Möbius transformations from a prior section, we are now

ready to derive the associated infinitesimal generators of the Möbius group action on

the ball Bd. We will show that they correspond to flows of the form

ẏ = Ay − ⟨Z, y⟩y +
1
2

(1 + |y|2)Z, (2.3)

where A is an antisymmetric d × d matrix and Z ∈ Rd is a vector. Note that, as adver-

tised, this flow extends to a flow on S d−1 of the Kuramoto form in (2.2), as we can see

13


by restricting (2.3) to vectors y on the unit sphere where |y| = 1.

To derive (2.3), we work separately with the boost and rotation components. Let

us start with the boost component. Replace w by tw and expand Mtw(x) to first order in

t:

Mtw(x) ≈
x − |x|2tw

1 − 2t⟨w, x⟩
− tw ≈ x + t

(
2⟨w, x⟩x − (1 + |x|2)w

)
.

The derivative of this expression at t = 0 (which is just the coefficient of t) gives us the

infinitesimal generator: it is an “infinitesimal boost” of the form (2.3) with Z = −2w

and A = 0. Next, recall that the infinitesimal generators corresponding to the rotation

components are flows of the form ẋ = Ax for an antisymmetric matrix A. Together

with the infinitesimal boosts we then get all flows of the form (2.3). The group G acts

on the space X in the natural way (component by component) and the infinitesimal

generators of this group action on X are flows of the form (1) with all Ai identical.

Therefore, by the general philosophy discussed earlier, the evolution of any initial

point p ∈ X under the system (1) with all Ai = A lies in the group orbit Gp.

The given Kuramoto system has N(d − 1) degrees of freedom, for some large N.

However, since the flow of the system is determined via an action of the d(d + 1)/2

dimensional Lie group G, we can alternatively study the auxiliary dynamical system

on G, which we call the reduced equations. By ignoring rotations, we can further

restrict our attention to a system on the d-dimensional quotient G/S O(d) � Bd. The

dimensional reduction not only makes the reduced equations easier to analyze than

the original Kuramoto system, but the reduced equations require fewer computational

resources to numerically integrate.

As before, assume all the terms Ai in (1) are equal, so that Ai = A for some skew-

symmetric matrix A. Fix a base point p = (p1, . . . , pN) ∈ X. Then if the points pi are

in sufficiently general position, every element in the G-orbit of p can be expressed

14


uniquely as gp for some g ∈ G, with parameters w, z and ζ. We wish to derive the

corresponding evolution equations for w, z and ζ. Let (xi(t)) be any solution to (1)

in the group orbit Gp; we do not require that the initial point (xi(0)) = p. Then for

i = 1, . . . ,N we have xi(t) = gt(pi) for a unique gt ∈ G, which determines the parameters

w, z, ζ as functions of t. Now consider the equation (2.3), with coefficients A and Z

evaluated at (xi(t)). This is a non-autonomous ODE on Bd, and its time-t flow must

be given by some g̃t ∈ G. This ODE has solutions xi(t) = gt(pi) = gt(g−1
0 (xi(0))), which

implies that g̃t = gtg−1
0 .

So for any y0 ∈ Bd,

y(t) = gt(g−1
0 (g0(y0))) = gt(y0) = ζMw(y0) = M−z(ζy0)

must satisfy the ODE (2.3) with A and Z evaluated at (xi(t)) at time t. In particular, if

we let y0 = 0, then y(t) = −ζw = z, so z satisfies the ODE (2.3).

Now expand y = ζMw(y0) = M−z(ζy0) to first order in y0, using the variables z and ζ:

y ≈ z + (1 − |z|2)ζy0,

so

ẏ ≈ ż − 2⟨ż, z⟩ζy0 + (1 − |z|2)ζ̇y0.

On the other hand, (2.3) gives

ẏ = Ay +
1
2

(
1 + |y|2

)
Z − ⟨Z, y⟩y

≈ Az +
1
2

(1 + |z|2)Z − ⟨Z, z⟩z

+ (1 − |z|2) (Aζy0 + ⟨z, ζy0⟩Z − ⟨Z, z⟩ζy0 − ⟨Z, ζy0⟩z) .

Setting y0 = 0 gives the ż equation

ż = Az +
1
2

(1 + |z|2)Z − ⟨Z, z⟩z (2.4)

15


as expected, and since ⟨Az, z⟩ = 0, this in turn implies that

⟨ż, z⟩ =
1
2

(1 − |z|2)⟨Z, z⟩.

Equating the y0 terms, factoring out 1 − |z|2 and canceling the common term ⟨Z, z⟩ζy0

gives

ζ̇y0 = Aζx0 + ⟨z, ζy0⟩Z − ⟨Z, ζy0⟩z.

Together, the last two terms above define a special type of antisymmetric operator of

ζy0: Given any y1, y2 ∈ Rd, define the antisymmetric operator α as

α(y1, y2)y = ⟨y1, y)y2 − ⟨y2, y)y1;

this operator has range = span(y1, y2) providing y1 and y2 are linearly independent;

otherwise α(y1, y2) = 0. Then for all y0 ∈ Rd,

ζ̇y0 = Aζy0 + α(z,Z)ζy0

and therefore

ζ̇ = (A + α(z,Z))ζ.

Differentiating z = −ζw gives

Az +
1
2

(1 + |z|2)Z − ⟨Z, z⟩z = −ζẇ − ζ̇w

so

ζẇ = (A + α(z,Z))z − Az −
1
2

(1 + |z|2)Z + ⟨Z, z⟩z

= Az + |z|2Z − ⟨Z, z⟩z − Az −
1
2

(1 + |z|2)Z + ⟨Z, z⟩z

= −
1
2

(1 − |z|2)Z;

hence

ẇ = −
1
2

(1 − |w|2)ζ−1Z. (2.5)

16


Summing up, the evolution equations for the (z, ζ) coordinate system on Gp are

ż = Az +
1
2

(1 + |z|2)Z − ⟨Z, z⟩z (2.6a)

ζ̇ = (A + α(z,Z))ζ, (2.6b)

with A and Z evaluated at M−z(ζp), and for the (w, ζ) coordinate system on Gp are

ẇ = −
1
2

(1 − |w|2)ζ−1Z (2.7a)

ζ̇ = (A − α(ζw,Z))ζ, (2.7b)

with A and Z evaluated at ζMw(p). Note that these equations generalize the evolution

equations for the parameters w and ζ given in Chen et al.[20] for the classic case d = 2.

2.6 Comparison of Z versus W coordinates

The ż equation (2.4) is an extension of the system equation (2.1) on S d−1. However,

for finite N, the ż equation does not uncouple from ζ, since Z is evaluated at M−z(ζp).

The exception to this is in the infinite-N limit: if the base point p is now the uniform

density on S d−1, then ζp = p (the uniform density is invariant under rotations) and

the density M−z(p) is a hyperbolic Poisson density on S d−1 whose centroid is a function

of z. In the case d = 2, this Poisson density has centroid z. Unfortunately this simple

relationship is false for d ≥ 3 (we will give more details on this in the next section).

The advantage of the ẇ equation (2.5) is that for an order parameter function of the

form

Z =
N∑

i=1

aixi,

with ai ∈ R, ζ drops out of the ẇ equation and we get the reduced equation

ẇ = −
1
2

(1 − |w|2)Z(Mw(p)).

17


The parameter w essentially defines the “phase relations” among the xi; two configu-

rations have the same w if and only if they are related by a rotation. So w is the key

parameter that determines whether the system is approaching synchrony or incoher-

ence.

The w variable also has a nice invariance under change of base points. Suppose

p′ = M(p) ∈ Gp; then we have coordinates w′, ζ′ associated to the base point p′. Any

q ∈ Gp has two expressions

q = ζMw(p) = ζ′Mw′ p′ = ζ′Mw′(M(p)).

Assuming the coordinates of p are in sufficiently general position, this implies ζMw =

ζ′Mw′ ◦ M, and hence

0 = ζMw(w) = ζ′Mw′(M(w)).

But the unique solution to Mw′(y) = 0 is w′, and hence w′ = M(w). In other words, the

coordinates w and w′ transform exactly as the base points p and p′.

2.7 Continuum limit

Next, we consider the dynamics of the Kuramoto model (2.1) in the limit N →

∞. We assume that the rotation terms Ai = A are constant across the population,

corresponding to identical “oscillators.”

Let us also assume that the order parameter Z is is proportional to the centroid of

the population:

Z =
K
N

N∑
i=1

xi

In the continuum limit, a state of the system is a probability measure ρ on S d−1, and

the order parameter becomes

Z = K
∫

S d−1
x dρ(x).

18


The measure ρ evolves according to the continuity equation (also known as the noise-

less Fokker-Planck equation) associated to the flow in (1). Naturally, this flow must

preserve group orbits under the action of G. Recall that if M ∈ G, then the measure

M∗ρ is defined by the adjunction formula∫
S d−1

f (x) d(M∗ρ)(x) =
∫

S d−1
f (M(x)) dρ(x).

In particular, we can consider the G-orbit of the uniform probability measure σ on

S d−1. This orbit is special; whereas a typical group orbit Gρ has dimension equal to

the dimension of G, namely d(d + 1)/2, the orbit Gσ has dimension only d. This is

because the stabilizer of σ is S O(d); any rotation fixes σ, whereas the boosts deform

σ. Hence the orbit Gσ has dimension d. Any element in Gσ can be written as (M−z)∗σ,

with z ∈ Bd. The evolution equation for z is (2.4), with

Z(z) = K
∫

S d−1
x d(M−z)∗σ(x) = K

∫
S d−1

M−z(x) dσ(x). (2.8)

In the case d = 2 with x = ζ ∈ S 1, we have

dσ(ζ) =
1

2πi
dζ
ζ
,

so the integral

Z(z) =
K

2πi

∫
S 1

ζ + z
1 + zζ

·
dζ
ζ
= K

ζ + z
1 + zζ

∣∣∣∣∣∣
ζ=0

= Kz

by the Cauchy integral formula. Therefore (2.4) simplifies to the equation

ż = iωz +
K
2

(1 − |z|2)z

when d = 2. Unfortunately, the formula Z(z) = Kz is not correct for d ≥ 3; though as we

shall see later, this formula is correct in higher dimensions for the complex hyperbolic

model in even dimensions, which we discuss in the next section. For d = 2 the two

geometries agree, which explains the coincidence for d = 2.

19


Any Riemannian manifold X has a Laplace-Beltrami operator ∆ associated to its

metric; functions f on X satisfying the equation ∆ f = 0 are called harmonic. For func-

tions on the ball Bd with the hyperbolic metric, this operator is

∆hyp = (1 − |x|2)2∆euc + 2(d − 2)(1 − |x|2)
d∑

i=1

xi
∂

∂xi
,

where

∆euc =

d∑
i=1

∂2

∂x2
i

is the standard Laplace operator (see Stoll[118], Chapter 3). We will call solutions

to the equation ∆hyp f = 0 hyperbolic harmonic functions; for d = 2 these coincide with

ordinary (Euclidean) harmonic functions. We can consider the hyperbolic analogue of

the classical Dirichlet problem: given a continuous function f on S d−1, extend f to a

hyperbolic harmonic function f̃ on Bd. Assuming this problem has a unique solution,

then for any rotation ζ ∈ S O(d) we must have f̃ ◦ ζ = f̃ ◦ ζ, since rotations preserve the

hyperbolic metric. If we average f ◦ ζ on S d−1 over all rotations ζ ∈ S O(d) we get the

constant function

fave =

∫
S d−1

f (x) dσ(x)

on S d−1, and any constant is hyperbolic harmonic on Bd. Therefore the average on Bd

of f̃ ◦ ζ = f̃ ◦ ζ over all ζ ∈ S O(d) must be the constant fave. But f̃ (ζ(0)) = f̃ (0) for all ζ,

so we must have

f̃ (0) =
∫

S d−1
f (x) dσ(x).

Now let z ∈ Bd; since M−z preserves the hyperbolic metric, we must have f̃ ◦ M−z =

f̃ ◦ M−z, which implies

f̃ (z) = f̃ ◦ M−z(0)

=

∫
S d−1

f (M−z(x)) dσ(x)

=

∫
S d−1

f (x) d((M−z)∗σ)(x).

20


As shown in Chapter 5 in Stoll [118], the measure (M−z)∗σ is given by the formula

d((M−z)∗σ)(x) = Phyp(z, x) dσ(x),

with hyperbolic Poisson kernel function

Phyp(z, x) =
(
1 − |z|2

|z − x|2

)d−1

. (2.9)

Thus the solution to the hyperbolic Dirichlet problem with boundary function f on

S d−1 is given by the hyperbolic Poisson integral

f̃ (z) =
∫

S d−1
Phyp(z, x) f (x) dσ(x), z ∈ Bd.

The orbit Gσ consists of all hyperbolic Poisson measures P(z, x) dσ(x), parametrized by

z ∈ Bd. By contrast, the Euclidean Poisson kernel function is

Peuc(z, x) =
1 − |z|2

|z − x|d
,

so the hyperbolic Poisson measures agree with the Euclidean Poisson measures only

if d = 2.

Now we can calculate the expression Z(z) in the general case d ≥ 2. We see from

(2.8) that Z(z) is the hyperbolic Poisson integral of the function Kx on S d−1. The func-

tion Kx is (Euclidean) harmonic and homogeneous of degree 1 on Rd; following the

recipe in Chapter 5 in Stoll [118], we see that its extension from S d−1 to a hyperbolic

harmonic function on Bd is given by

Z(z) = K
F(1, 1 − d/2; 1 + d/2; |z|2)
F(1, 1 − d/2; 1 + d/2; 1)

z, (2.10)

where F is the hypergeometric function

F(a, b; c; t) =
∞∑

k=0

(a)k(b)k

(c)k

tk

k!
,

with (a)0 = 1 and (a)k = a(a + 1) · · · (a + k − 1) for k ≥ 1. Notice that if a or b = 0, then

F(a, b; c; t) = 1; this gives Z(z) = Kz for d = 2, as expected.

21


2.8 Complex case

There is an alternative generalization of Kuramoto networks to higher-

dimensional oscillators when d = 2m is even. Then Rd = Cm, and we can study systems

of the form

ẋ j = A jx + Z − ⟨x j,Z⟩x j, i = 1, . . . ,N, (2.11)

where now xi is a point on the unit sphere S 2m−1 ⊂ Cm, Ai is an anti-Hermitian m × m

complex matrix, Z ∈ Cm and ⟨, ⟩ denotes the complex-valued Hermitian inner product.

These systems are the same as the real case when d = 2,m = 1 but are different for

m ≥ 2. To see this, suppose

Ax + Y − ⟨x,Y⟩R x = Bx + Z − ⟨x,Z⟩C x

for all x ∈ S 2m−1 ⊂ Cm = Rd, where A is antisymmetric, B is anti-Hermitian, Y,Z ∈ Cm

and we use the subscripts R and C to distinguish the real and complex inner products.

Then

(A − B)x = Z − Y +
(
⟨x,Y⟩R − ⟨x,Z⟩C

)
x

and so (A − B)(−x) = (A − B)x for all x ∈ S 2m−1, which implies A = B. This implies

Y − Z =
(
⟨x,Y⟩R − ⟨x,Z⟩C

)
x

for all x ∈ S 2m−1, hence Y − Z ∈ spanC(x) for all x ∈ Cm; if m ≥ 2, this implies Y = Z. But

then we have

⟨x,Y⟩R = ⟨x,Y⟩C

for all x ∈ Cm, which can only hold if Y = 0. Hence for m ≥ 2, the only flows simulta-

neously of the form (1) and (2.11) have Z = 0 and A anti-Hermitian.

Flows of the form (2.11) are related to the complex hyperbolic geometry on the

complex unit ball Bm with the Bergman metric (see Rudin[110], Chapter 1). The

22


orientation-preserving isometries of this metric are generated by unitary transforma-

tions ζ ∈ U(m) and boost transformations of the form

Mw(x) =

√
1 − |w|2 x +

(
⟨x,w⟩

1+
√

1−|w|2
− 1

)
w

1 − ⟨x,w⟩

=

x − w + ⟨x,w⟩w−|w|
2 x

1+
√

1−|w|2

1 − ⟨x,w⟩
.

Notice that when m = 1, this reduces to the standard complex Möbius map Mw. As

in the real case M0 is the identity, M−1
w = M−w, Mw(w) = 0 and Mw(0) = −w. Any

orientation-preserving isometry of Bd can be expressed uniquely in the form

g(x) = ζMw(x) = M−z(ξx),

where w, z ∈ Bm but now ζ, ξ ∈ U(m), the complex unitary group. Linearizing at x = 0

gives

g(x) ≈ ζ

−w − ⟨x,w⟩w +
√

1 − |w|2 x +
⟨x,w⟩w

1 +
√

1 − |w|2


≈ ζ

−w +
√

1 − |w|2 x −

√
1 − |w|2⟨x,w⟩w

1 +
√

1 − |w|2


≈ z +

√
1 − |z|2 ξx −

√
1 − |z|2⟨ξx, z⟩z

1 +
√

1 − |z|2

which implies z = −ζw (hence |z| = |w|) and ξ = ζ, as before. The corresponding

infinitesimal transformations are given by flows on the complex unit ball Bm of the

form

ẏ = Ay + Z − ⟨y,Z⟩y, (2.12)

with A anti-Hermitian m × m and Z ∈ Cm. This flow extends to a flow on S 2m−1 of

the form in (2.11). Note the absence of the quadratic term |y|2Z here. To derive these

infinitesimal transformations, we can apply our prior power series expansion, noting

that |x| = 1 to obtain

Mtw(x) ≈
x − tw

1 − t⟨x,w⟩
≈ x + 2 (⟨x,w⟩x − w) t.

23


Figure 2.1: A first-order linear Kuramoto system on the two-dimensional sphere S 2

with equal weights ai = 1/N, and randomly chosen initial conditions. The states
shown are at t = 0, t = 10, and t = 40 respectively. This simulation was written in
Python and visualized with Plotly.

So the infinitesimal generator is an “infinitesimal boost” of the form (2.12) with Z =

−1
2w and A = 0. The infinitesimal generators corresponding to the rotation components

are flows of the form ẋ = Ax with A anti-Hermitian; together with the infinitesimal

boosts we get all flows of the form (2.12).

We mention in passing that the complex model (2.11) has a natural quantum net-

work analogue, studied by Lohe [73], in which the complex vector xi is replaced by a

normalized wave function |ψi⟩. We expect that our reduction techniques extend to this

infinite-dimensional quantum model.

2.9 Relation to previous research

Many of the results above can be found in some form in the work of earlier au-

thors. [72, 123, 23, 43, 74, 52, 18, 17, 75, 41, 53] Three papers in particular overlap

considerably with the present work.

Tanaka [123] demonstrates in his 2014 paper that the dynamics of (2.1) can be re-

duced using Möbius transformations that fix the unit ball, similar to what Marvel,

Mirollo, and Strogatz found[80] for the traditional Kuramoto model. Tanaka writes

his Möbius transformations differently from ours, but he uses the same group of trans-

24


formations and he also gets reduced equations for his Möbius parameters. Tanaka’s

equation (10b) looks similar to our ż equation (2.4), except without the |z|2 term, which

is puzzling. He does not mention the reduction down to dimension d in the finite-N

case that we get with the ẇ equation (2.5). Tanaka also notes that the complex case

when d = 2m is different, and generalizes the Ott-Antonsen residue calculation to this

case, which is the highlight of his paper. In the real case, Tanaka’s equation (15) is sim-

ilar to our equation (2.10), though we were not able to show that the two expressions

are equivalent. Finally, Tanaka also presents a generalization of the Ott-Antonsen

reduction[96] for the complex version of the system.

Lohe [74] also looks at the same system as (1) (see his equation (22)) and he derives

a similar reduction as ours by using Möbius transformations for the finite-N model.

His transformation (30) on S d−1 is our Mw (with v = w) and his equation (31) is the

same as our ż equation (2.4). He also has something that looks like the ẇ equation (2.5),

which he says is independent of the rest of the reduced system for (in our notation) an

order parameter function of the form

Z =
1
N

N∑
i=1

λiQixi,

where Qi ∈ O(d) and λi ∈ R. But such a Z does not satisfy the identity ζZ(p) = Z(ζp) for

all rotations ζ, unless Qi = ±I, so we do not see how the ζ term cancels.

Lohe’s map M in his equation (55) (ignoring the R factor) agrees with our map M−v

on the sphere S d−1, but not on the ball Bd. So it is not a Möbius transformation of the

type we are using. For example, M(−v) = v whereas M−v(−v) = 0. We are not sure

why Lohe [74] prefers these maps over the boosts; he claims that M preserves cross-

ratios, but we do not see why this is advantageous. His map F in equation (63) (again

ignoring the R factor) is exactly our M−v.

Chandra, Girvan, and Ott [17] concentrate on the infinite-N or continuum limit

25


Figure 2.2: A first-order linear Kuramoto system on S 2 with weights distributed ac-
cording to a Riemann sum which approximates the integral of a normal distribution,
and randomly chosen initial conditions. Pink particles contribute to the order parame-
ter with greater weights than the blue particles do. The states shown are at t = 0, t = 10,
and t = 40 respectively.

system, and derive a dynamical reduction for a special class of probability densities

on S d−1, generalizing the Poisson densities used in the Ott-Antonsen reduction. They

proceed directly to the infinite-N version of (1). They make a very clever guess (their

equation (7)) of the form of the special densities that generalize the Poisson densities

for d = 2, and then calculate the exponent in the denominator of their expression,

getting exactly the hyperbolic Poisson kernel densities in (2.9) above. Their equation

(15) is exactly the same as our ż equation (2.4) in the infinite-N limit. The integral in

their equation (19) can be evaluated, as shown above in (2.10).

2.10 An example: First-order linear order parameter gives gradient

system

We conclude with an analysis of the system (1) with a weighted order parameter

Z =
N∑

i=1

aixi, (2.13)

where the ai are real constants.

26


Figure 2.3: A first-order linear Kuramoto system on S 2 with a majority cluster, where
one particle is chosen to have a weight which exceeds the combined weights of all
other particles (or equivalently, where all the particles have equal weight but a ma-
jority of them cluster into a single point and therefore act if they were a single giant
particle; hence the name “majority cluster”). The states shown are at t = 0, t = −10,
and t = −40 respectively; we have chosen to depict time running backward to high-
light that the backwards-time limit tends toward an antipodal configuration. In this
simulation, one particle, depicted in pink, was chosen to have a weight of 0.6, and the
remaining 99 particles, depicted as blue, were chosen to have equal weights of 0.4/99.

2.10.1 Computer visualization

We have been discussing a system of particles on a sphere that coalesce (in other

words, they spontaneously synchronize) when certain conditions are met. To visualize

this coalescence, we implemented the Runge-Kutta algorithm to numerically solve the

Kuramoto model on the 2-dimensional sphere in 3-space, corresponding to the case

d = 3 in (2.1). For simplicity, we simulated N = 100 particles with equal weights (ai =

1/N in the order parameter Z) and set the rotation term A to zero for all the particles, a

choice that is tantamount to ignoring the rotational influence, or equivalently, rotating

the frame of reference along with the entire system as it evolves.

In the simulation shown in Fig. 2.2, randomly chosen points on the sphere were

used as initial conditions. As time increases, one can see that the particles coalesce to

a limit point, mimicking the spontaneous synchronization that is well known for the

traditional Kuramoto model (d = 2) when the oscillators are identical.

Later in this section, we will prove this synchronization behavior holds more gen-

27


erally for Kuramoto models on the sphere having weighted order parameters of the

form (2.13), provided the weights ai are all positive and sum to 1 and no individual

weight exceeds 1/2. A partial result in this direction was obtained by Choi et al. [24]

These authors prove that in the case where all ai are equal and positive, initial con-

ditions satisfying ⟨xi(0),Z(0)⟩ > 0 for all i will synchronize. Geometrically, this condi-

tion is equivalent to requiring that the oscillators all lie on the hemisphere given by

⟨xi, y⟩ > 0 for some vector y , 0. More generally, similar partial synchronization results

are obtained by Ha et al. [41] for the case

Z = (aI +W)
N∑

i=1

xi,

with W skew-symmetric, and by Ha and Park [44] for the complex system (2.11) with

Z =
∑

xi.

Figure 2.2 shows a simulation in which we weighted each particle according to the

terms in a Riemann sum approximating the integral of the normal probability distri-

bution. The pink particles, which have higher weights, exert greater influence over

the final synchronization location of the particles, but there is still synchronization.

When time runs backwards, almost all initial conditions of the particles tend to-

wards a limiting configuration where their centroid is at the origin. The exception

is when we have a majority cluster, depicted in Fig. 2.3, where one particle has a

weight which exceeds the weight of all other particles. When this condition holds, it

is impossible to arrange the particles so their weighted centroid is at the origin, so the

backwards time limit will tend towards an antipodal configuration, where all parti-

cles not in the majority cluster will coalesce around the antipode of the cluster. This

is the configuration which minimizes the magnitude of the weighted centroid. We do

not include the proof that the backwards-time limit is antipodal in this case, but it is a

straightforward generalization of the result which we do prove below.

28


2.10.2 Existence of hyperbolic gradient

As mentioned above, if Z has the form in (2.13), then the ẇ equation in (2.5) reduces

to

ẇ = −
1
2

(1 − |w|2)Z(Mw(p)), (2.14)

independent of the parameter ζ. We will show that this is a gradient flow on the

unit ball Bd with respect to the hyperbolic metric. In the presence of a Riemannian

metric we can associate a 1-form to any vector field, and the vector field is gradient if

and only if the associated 1-form is exact; since the unit ball is simply connected, this

holds if and only if the associated 1-form is closed. For the Euclidean metric on Bd (or

any open subset of Rd) and standard coordinates w1, . . . ,wd, the 1-form associated to

the vector field with components f1, . . . , fd is

ω = f1 dw1 + · · · + fd dwd.

If we scale the Euclidean metric by a positive smooth function φ, then the associated

1-form with respect to the metric ds = φ|dw| is then

ω = φ2( f1 dw1 + · · · + fd dwd).

Therefore the gradient of a function Φwith respect to this scaled metric is given by

∇Φ = φ−2∇eucΦ,

where ∇euc denotes the ordinary Euclidean gradient operator. We have φ(w) = 2(1 −

|w|2)−1 for the hyperbolic metric, so the hyperbolic gradient operator on Bd is given by

∇hypΦ(w) =
1
4

(1 − |w|2)2∇eucΦ(w).

Now let’s consider the vector field V defined by (2.14). By linearity, it suffices to

treat the case Z = xi, and we can take i = 1 without loss of generality. Then the

associated 1-form is

29


ω =
4

(1 − |w|2)2

(
−

1
2

(1 − |w|2)
)

·

d∑
j=1

(
(1 − |w|2)(p1, j − w j)

|p1 − w|2
− w j

)
dw j

= −2
d∑

j=1

(
p1, j − w j

|p1 − w|2
−

w j

1 − |w|2

)
dw j

where p1, j denotes the jth component of the point p1 ∈ S d−1. Let E j denote the coeffi-

cient of dw j in parentheses above; then

dω = −2
d∑

j,k=1

∂E j

∂wk
dwk ∧ dw j.

Applying the chain and quotient rules gives

∂E j

∂wk
=

2(p1, j − w j)(p1, k − wk)
|p1 − w|4

+
2w jwk

(1 − |w|2)2

for j , k, which is symmetric in j and k; hence the sum above for dω simplifies to

dω = 0. Thus ω is closed and we see that the flow (2.14) is gradient for any order

parameter function of the form (2.13).

Next, we show that the hyperbolic potential for V , up to an additive constant, is

given by

Φ(w) =
N∑

i=1

ai log
1 − |w|2

|w − pi|
2 =

1
d − 1

N∑
i=1

ai log Phyp(w, pi). (2.15)

Here we follow the convention that the potential decreases along trajectories, so we are

asserting that ∇hypΦ = −V . To derive this, use the identity ∇euc|w − w0|
2 = 2(w − w0), for

30


any constant vector w0 ∈ Rd. Then

∇eucΦ(w) =
N∑

i=1

ai

(
−

2w
1 − |w|2

−
2(w − pi)
|w − pi|

2

)

=
2

1 − |w|2

N∑
i=1

ai

(
(1 − |w|2)(pi − w)
|w − pi|

2 − w
)

=
2

1 − |w|2

N∑
i=1

aiMw(pi) =
2

1 − |w|2
Z(Mw(p)).

Hence we see that

∇hypΦ(w) =
1
2

(1 − |w|2)Z(Mw(p)) = −V(w),

as desired.

2.10.3 Analysis of dynamics

We can use the existence of the potential Φ(w) for the flow on Bd to prove a global

synchrony result for the system (2.1) when the coefficients ai in the order parameter Z

are all positive. Specifically, we assume that 0 < ai < 1/2 for all i, and
∑N

i=1 ai = 1. We

also assume N ≥ 3 and all the rotation terms Ai in (2.1) are equal. Under these condi-

tions, almost all trajectories for (1) converge in forward time to the (d−1)-dimensional

diagonal manifold ∆ ⊂ X as t → ∞, meaning that the system self-synchronizes. In

contrast, in backwards time the system tends to an incoherent state having zero order

parameter: as t → −∞ almost all trajectories for (1) converge to the codimension-d

subspace Σ ⊂ X consisting of states with Z(p) = 0.

The proof is modeled after Theorem 1 in Chen et al. [21] and will be based on two

preliminary lemmas. In each of these lemmas we assume the conditions on the ai

above, and that the base point p = (pi) for the flow (2.14) has all distinct coordinates.

We begin with a general observation about gradient flows in the ball Bd: if w0 ∈ Bd

is any initial condition and w∗ ∈ Bd is in the forward limit set Ω+(w0), then w∗ is a fixed

31


point for the flow. To see this, let Φ be a potential for the flow, and suppose w(tn) →

w∗ ∈ Bd for some sequence tn → ∞. Since the potential decreases along trajectories,

lim
t→∞
Φ(w(t)) = lim

n→∞
Φ(w(tn)) = Φ(w∗).

Let Ft denote the time-t flow map. If w∗ is not a fixed point, then for any s > 0,

lim
t→∞
Φ(w(t)) = lim

n→∞
Φ(w(tn + s))

= lim
n→∞
Φ(Fs(w(tn)))

= Φ(Fs(w∗))

< Φ(w∗),

which is a contradiction, so w∗ must be a fixed point. (Compact limit sets are con-

nected, so Ω+(w0) cannot consist of two or more but finitely many fixed points; how-

ever it is possible that forward or backward limits sets for gradient flows consist of a

continuum of fixed points. We will see that this is not the case for our system on Bd.)

Lemma 2.10.1. Any fixed point for the flow (2.14) in Bd is repelling.

Proof. Suppose w∗ ∈ Bd is a fixed point for (2.14). As discussed above, an advantage

of using the w-parameter is the equivariance with respect to change of base point p.

Consequently we can assume w∗ = 0 without loss of generality, so Z(p) =
N∑

i=1
ai pi = 0.

To first order in w,

Mw(pi) =
pi − w

1 − 2⟨w, pi⟩
− w

= (pi − w)
(
1 + 2⟨w, pi⟩

)
− w

= pi − 2w + 2⟨w, pi⟩pi.

32


The linearization of (2.14) at the fixed point w∗ = 0 is

ẇ = −
1
2

N∑
i=1

ai

(
pi − 2w + 2⟨w, pi⟩pi

)
= w −

N∑
i=1

ai⟨w, pi⟩pi.

We claim that the linear map

Tw =
N∑

i=1

ai⟨w, pi⟩pi

has ||T || < 1; to see this, suppose |w| = 1. Then |⟨w, pi⟩pi| ≤ 1 and Tw is a convex

combination of the vectors ⟨w, pi⟩pi. We can only obtain |Tw| = 1 if all terms ⟨w, pi⟩pi =

u with |u| = 1, which implies all pi = ±u, and this cannot happen if at least three of

the pi are distinct. Hence ||T || < 1 and so the eigenvalues µi of T satisfy |µi| < 1. The

eigenvalues for the ẇ linearization are λi = 1 − µi, so wee see that Re λi > 0 for all i,

establishing that the fixed point w∗ is repelling. □

Lemma 2.10.2. lim
|w|→1
Φ(w) = −∞.

Proof. It suffices to show that

lim
n→∞
Φ(wn) = −∞

for any sequence wn ∈ Bd with wn → x ∈ S d−1. The result is clear if x , pi: as n → ∞

the terms |wn − pi| in the potential (3.1) are bounded away from 0, and 1− |wn|
2 → 0. So

let’s say that wn → p1. We rewrite Φ(wn) as

33


Φ(wn) = log(1 − |wn|
2) − 2a1 log |wn − p1|

− 2
N∑

i=2

ai log |wn − pi|

= log(1 − |wn|) − 2a1 log |wn − p1|

+ log(1 + |wn|) − 2
N∑

i=2

ai log |wn − pi|.

The latter two terms above have finite limit as n → ∞, so we focus on the first two

terms. We have 1 − |wn| ≤ |wn − p1|, so

log(1 − |wn|) − 2a1 log |wn − p1| ≤ (1 − 2a1) log |wn − p1| → −∞

as n → ∞, which proves our result. Notice that we need the assumption ai < 1/2 for

this argument. □

Theorem 2.10.3. Under the conditions above, almost all trajectories for (2.1) converge to ∆

as t → ∞ and to Σ as t → −∞.

Proof. Let p = (p1, . . . , pN) ∈ X be any point with all distinct coordinates. The points

on Gp are parametrized by w ∈ Bd and ζ ∈ S O(d), and the dynamics for these pa-

rameters are given by (2.5). We begin with the dynamics as t → −∞. Let w(t) be a

trajectory for (2.14) with initial condition w0 ∈ Bd, and consider the backward time

limit set Ω−(w0); this is a nonempty, compact, connected subset of Bd. The potential Φ

is decreasing along all trajectories w(t), hence bounded below as t → −∞, so Lemma

2 implies that the limit set Ω−(w0) must be contained in the interior Bd. We know that

any w∗ ∈ Ω−(w0) is a fixed point for the flow. By Lemma 1, w∗ is repelling and so any

trajectory w(t) which comes sufficiently close to w∗ must have w(t) → w∗ as t → −∞;

therefore Ω−(w0) = {w∗}. This proves the existence of fixed points for (2.14), and that

every trajectory w(t) converges to a fixed point as t → −∞. If the flow had multiple

34


fixed points, we would obtain a partition of Bd into the disjoint open basins of repul-

sion of the fixed points, violating connectedness of the ball. Therefore (2.14) has a

unique fixed point w∗, and w(t) → w∗ as t → −∞ for all trajectories. The fixed point w∗

has Z(Mw∗(p)) = 0, so all trajectories in Gp converge to Σ as t → −∞.

In forward time, the limit set Ω+(w0) for any w0 , w∗ must be completely contained

in the boundary S d−1, since the unique fixed point w∗ ∈ Bd is repelling. Suppose we

remove the factor (1/2)(1 − |w|2) in the flow (2.14); the scaled vector field on Bd given

by

ẇ = −
N∑

i=1

aiMw(pi) = w −
N∑

i=1

ai

(
(1 − |w|2)(pi − w)
|pi − w|2

)
(2.16)

has the same trajectories as the original flow, just with different time parametrizations.

Observe that this scaled vector field extends smoothly to Rd − {pi}, and coincides with

the radial vector field x at any x ∈ S d−1 with x , pi. Therefore there is a unique

trajectory passing through each point x ∈ S d−1, flowing from the interior to the exterior

of the sphere, as long as x , pi. Consequently the original flow (2.14) has a unique

trajectory w(t) in Bd with w(t) → x as t → ∞, as long as x , pi. This also shows

that there is a neighborhood U of S d−1 − {pi} such that if w(t0) ∈ U for some t0, then

w(t) → x , pi for some x ∈ S d−1. So if Ω+(w0) contains some x , pi, then the trajectory

w(t) of w0 must enter the neighborhood U, and therefore w(t)→ x ∈ S d−1 as t → ∞.

Since limit sets are connected, the only other possibility is Ω+(w0) = {pi} for some

i; equivalently, w(t) → pi. We will show that there is a unique trajectory with this

behavior for each pi. Assuming this, we see that with N + 1 exceptions, any trajectory

w(t) converges to a point x ∈ S d−1 with x , pi (the exceptions are the N trajectories

converging to the base point coordinates pi, and the fixed point trajectory w∗). The

corresponding trajectory in Gp has coordinates

ζ(t)Mw(t)(pi) = ζ(t)
(
(1 − |w(t)|2)(pi − w(t))

|pi − w(t)|2
− w(t)

)
.

35


We have |w(t)| → 1 and |pi − w(t)| is bounded away from 0 as t → ∞, so Mw(t)(pi) → −x

for each i and therefore the trajectory ζ(t)Mw(t)(p) in Gp converges to ∆ as t → ∞.

This analysis breaks down at x = pi because the scaled vector field above does

not have a unique limit as w → pi; rather, its limit depends on the direction of the

approach. To see this, write w = p1 − ru, where 0 < r < 1 and |u| = 1 (with this

convention, u = p1 corresponds to w approaching p1 radially). Then |p1 − w| = r and

|w|2 = 1 − 2r⟨p1, u⟩ + r2

so
(1 − |w|2)(p1 − w)
|p1 − w|2

=
(2r⟨p1, u⟩ − r2)ru

r2 = (2⟨p1, u⟩ − r) u.

As r → 0, the magnitude of this term is 2⟨p1, u⟩, which depends on the angle of ap-

proach given by u (note that ⟨p1, u⟩ > 0 because u points outwards at p1).

To complete the proof, we will examine the scaled system (2.16) using the polar

representation (r, u), and show that the polar system has the unique fixed point r∗ =

0, u∗ = p1, which has a unique attracting trajectory because it is a saddle with a (d − 1)-

dimensional unstable manifold.

We see that the scaled system has

ẇ = p1 − ru − a1 (2⟨p1, u⟩ − r) u + O(r)

= p1 − 2a1⟨p1, u⟩u + O(r),

where the O(r) term is a smooth function of r and u for |r| < ε = min |pi − p1|, i ≥ 2.

This condition insures that |pi − w| ≥ |pi − p1| − |r| > 0, so the i ≥ 2 terms in the scaled

ẇ equation are all smooth functions of r and u. And we can allow r < 0 here, even

though it is not relevant to the ẇ system. Now r2 = |w − p1|
2, so

rṙ = ⟨w − p1, ẇ⟩ = −r⟨u, ẇ⟩,

36


which gives

ṙ = −(1 − 2a1)⟨p1, u⟩ + O(r).

Differentiating ru = p1 − w gives

ru̇ = −ṙu − ẇ

= (1 − 2a1)⟨p1, u⟩u −
(
p1 − 2a1⟨p1, u⟩u

)
+ O(r)

= ⟨p1, u⟩u − p1 + O(r).

Hence the scaled system in polar form can be written

rṙ = −(1 − 2a1)r ⟨p1, u⟩ + O(r2),

ru̇ = ⟨p1, u⟩ u − p1 + O(r).

We emphasize that the O(r) and O(r2) terms are smooth functions of r, u as long as

|r| < ε. We consider the “semi-scaled” polar system

ṙ = −(1 − 2a1)r ⟨p1, u⟩ + O(r2), (2.17a)

u̇ = ⟨p1, u⟩ u − p1 + O(r), (2.17b)

which has the same trajectories as the original system, just with different time

parametrizations. The advantage of this modified system is that the equations are

smooth on (−ε, ε) × S d−1.

Observe that the system (2.17) has {0} × S d−1 invariant, and has fixed point (r∗, u∗) =

(0, p1). The fixed point p1 is repelling on the invariant manifold {0} × S d−1; to see this,

observe that

⟨p1, u⟩˙= ⟨p1, u⟩2 − 1

when r = 0. In fact, if we assign the coordinate θ on any great circle joining p1 and

−p1 on S d−1 so that u = eiθ and p1 = 1, then the system reduces to θ̇ = sin θ. We also

37


see that the ṙ equation linearized at (0, p1) is ṙ = −(1 − 2a1)r, so the linearization of

(2.17) has the single negative eigenvalue −(1 − 2a1) and d − 1 positive eigenvalues +1.

Therefore (0, p1) is a saddle with a one-dimensional stable manifold, and hence has a

unique trajectory (r(t), u(t))→ (0, p1) with r(t) > 0.

Now suppose we have a trajectory w(t) → p1 in our original system (2.14). The

corresponding trajectory for (2.17) will have r(t) → 0; we cannot achieve r(t) = 0 in

finite time because the manifold {r = 0} is invariant for (2.17). We must prove that

u(t)→ p1, so that this trajectory is in fact the saddle stable manifold. Observe that

⟨p1, u⟩˙= ⟨p1, u⟩2 − 1 + O(r).

Also note that ⟨p1, u(t)⟩ > 0 since |w(t)| < 1. Let 0 < c < 1; then for some T ≥ 0, t ≥ T

implies O(r(t)) ≤ (1 − c2)/2. Now suppose 0 < ⟨p1, u(t0)⟩ < c for some t0 ≥ T ; then

0 < ⟨p1, u(t)⟩c for all t ≥ t0. This is because the function t 7→ ⟨p1, u(t)⟩ is decreasing if

0 < ⟨p1, u(t)⟩ < c:

⟨p1, u(t)⟩˙≤ c2 − 1 +
1
2

(1 − c2) = −
1
2

(1 − c2)

as long as 0 < ⟨p1, u(t)⟩ < c. But this also implies that eventually ⟨p1, u(t)⟩ < 0, which

is a contradiction. Hence we must have ⟨p1, u(t)⟩ ≥ c for all t ≥ T , which proves that

u(t)→ p1. □

38


Part II

Electrostatic knot theory

39


CHAPTER 3

A LOWER BOUND ON THE CRITICAL POINTS OF THE ELECTRIC

POTENTIAL OF A KNOT

This chapter originally appeared as ”A lower bound on critical points of the electric

potential of a knot” in Journal of Knot Theory and its Ramifications, Vol. 30. [66]

3.1 Preliminaries and problem formulation

Our novel problem of interest is to analyze the zeros of the electric field around a

charged, knotted wire fixed in place. Let K ⊆ R3 ⊆ S 3 be a smooth knot parametrized

by the curve r(t), t ∈ [0, 2π] with r(0) = r(2π). We will take the convention that S 3 is

the union of R3 and a single compactifying point at infinity. Suppose K is endowed

with a uniform charge distribution. With a choice of units, the electric potential be-

tween a point k ∈ K and a point charge at x at a distance R from k is proportional to

R−1. It therefore makes sense to define the electric potential Φ : S 3 − K → R, on the

complement of K by the line integral

Φ(x) =
∫

k∈K

dk
|x − k|

=

∫ 2π

0

|r′(t)|
|x − r(t)|

dt, x ∈ R3 − K. (3.1)

We set Φ(∞) = 0 to ensure smoothness. By differentiating under the integral sign

with respect to x, we can see the electric potential is smooth and harmonic. The elec-

tric field is defined by E = −∇Φ. We want to describe the critical points of the potential

(equivalently, the zeros of the electric field) and their behaviors. These represent equi-

librium points where a charged particle at rest will continue to experience no electric

force from the charge distribution. Some conventions use the negative of the poten-

tial, so that the electric field points towards the knot, but it is more convenient for our

purposes to work with a nonnegative potential function.

40


Define the knot invariant cp(K) to be the smallest number of critical points of

the electric potential among all parametrizations in the knot isotopy class [K]. All

parametrizations have a critical point at infinity, which we include in the count. We

will now assume r is a parametrization which yields a critical set of minimal size. We

obtain a lower bound for the number of critical points of the electric potential based on

a well known topological invariant called the tunnel number t(K). The tunnel number

was originally introduced by Clark [25], and remains an active topic of knot theory

research [5].

We now come to the main result of this article. Proving this theorem uses Morse

theory and stable manifold theory.

Theorem 3.1.1. For all knots K, cp(K) ≥ 2t(K) + 2.

3.2 Preliminary definitions and lemmas

The theorems of Morse theory require us to work on a compact manifold, so in the

sequel we will define the knot complement of K in S 3 to be the complement of an open

tubular neighborhood of K of sufficiently small radius. The theorems from Morse

theory we invoke will still hold on this compact manifold with boundary because the

potential function is proper, with gradient transversely intersecting the boundary. We

will still denote our domain by S 3 − K.

Many of the following definitions and results are standard and can be found in

Nicolaescu [88] and Burde [15]. A critical point of a smooth real valued function f on

a manifold M is a point p such that the differential d f (p) is zero. The critical set of f is

the set of critical points, and is denoted Crit( f ). We say f is Morse if its critical points

are nondegenerate, which means the Hessian matrix H( f ) of second partial derivatives

is nonsingular. The index of p ∈ Crit( f ) is the number of negative eigenvalues of H( f ),

41


which is invariant under the choice of local coordinates. We denote the set of critical

points of index i by Criti( f ). If f is fixed, we denote the index of p by λ(p) and we

denote the number of critical points of index i by mi. In the space of all smooth real

valued maps, under a suitable function space topology, the set of Morse functions

is dense. Therefore, we may assume the electric potential Φ is Morse by adding a

perturbation if necessary.

We write WS (p) and WU(p) to denote the stable and unstable manifolds of p ∈

Crit( f ) respectively. Recall that the stable manifold (resp. unstable manifold) of p is

the set of all points which flow to p along the gradient vector field ∇ f as time tends to

infinity (resp. negative infinity). If f is Morse, the dimension of WS (p) is λ(p) and the

dimension of WU(p) is dim M − λ(p).

The Morse Reconstruction Theorem states the domain of a Morse function on a

compact manifold can be expressed as a cell complex by attaching closed discs with

dimensions given by the indices of the critical points. The attaching maps are obtained

by a process known as surgery, but we will not need to discuss the attaching maps in

any further detail.

The Morse inequalities state for a fixed Morse function f on a manifold M, mi ≥ bi,

where bi = dim Hi(M) is the ith Betti number of our domain manifold. We will need

the stronger result which states

dim M∑
i=0

(−1)imi =

dim M∑
i=0

(−1)ibi = χ(M), (3.2)

with χ(M) denoting the Euler characteristic of M. See Nicolaescu, Corollary 2.3.3 [88].

Before turning to the proof of Theorem 3.1.1, we need to prove a few preliminary

lemmas.

42


Lemma 3.2.1. For all knots K, Hi(S 3 − K) = Z, for i = 0, 1 and Hi(S 3 − K) = 0 for i ≥ 2.

Proof. See Rolfsen, Proposition 3.A.3 [109]. □

From the homology of the knot complement, along with (3.2), we can deduce the

following lemma.

Lemma 3.2.2. For all knots K, and with Φ defined above, m1 − m2 = 1.

Proof. Equation (3.2) states the Euler characteristic equals the alternating sum of the

mi’s. In other words, χ(S 3 − K) =
3∑

i=0
(−1)imi. From Lemma 3.2.2, we can see that χ(S 3 −

K) = 1 − 1 = 0. Since Φ is harmonic, every critical point has index 1 or 2, save for

the point at infinity, which has index 0. We can conclude m0 − m1 + m2 − m3 = 0, or

equivalently, m1 − m2 = 1 as desired. □

Remark 3.2.3. As m2 ≥ 0, m1 = 1 + m2 ≥ 1. That is, there is always a critical point of index

1.

The next set of definitions and results are standard in 3-manifold topology, and fur-

ther exposition can be found in [113] and [111]. A handlebody is a topological space

homotopic to the three dimensional ball with solid handles attached (by “attaching

handles” we mean there are copies of D2 × [0, 1] where the boundary discs D2 × 0 and

D2 × 1 are embedded on the boundary of the three-ball). Given a knot K ⊆ S 3, the

tunnel number t(K) is the least number of arcs we must add to K such that the com-

plement in S 3 is a handlebody. A collection of arcs with this property is known as a

tunneling.

In the proof of Theorem 3.1.1, we shall use the above definition of the tunnel num-

ber, but there is an equivalent definition that is more visually intuitive. A Heegaard

43


splitting of a three-manifold M is an embedding of a closed, compact, and orientable

surface H such that the interior and exterior of H in M are both handlebodies. We say

the genus of H is the genus of the splitting.

Theorem 3.2.4. Let H be an unknotted embedding of a genus g surface in S 3. That is, let H

be the boundary of a tubular neighborhood around a wedge sum of g unknotted circles. Then

H defines a Heegaard splitting.

Theorem 3.2.5. Genus g Heegaard splittings of S 3 are unique up to isotopy.

The previous result is known as Waldhausen’s Theorem. Clearly, a tunneling of a

knot defines a Heegaard splitting. Therefore, we can view the tunnel number of K as

the least number of arcs we must add to K such that it is isotopic to a wedge sum of

unknotted circles. We can immediately deduce that tunnelings and the tunnel number

always exist.

Lemma 3.2.6. Every smooth knot K has a tunnel number.

Proof. Take a diagram of K with finitely many crossings. Over each crossing, introduce

an arc connecting the top and bottom strands. Collapse each arc so that the top and

bottom strands intersect, and project onto the diagram’s plane so that we are left with

a wedge sum of say, g circles. By Theorem 3.2.4, it follows that the complement in S 3

is also a handlebody with g handles.

□

Remark 3.2.7. We just proved the tunnel number is bounded above by the crossing number,

the least number of crossings needed in a knot diagram of K.

As some elementary examples, the tunnel number of the unknot is zero, and the

tunnel number of the trefoil is one. Indeed, the tunnel number need not be the crossing

44


number. For example, torus knots have tunnel number 1, yet can have arbitrarily high

crossing number. See Clark [25].

3.3 Proof of the tunnel number bound

We now come to the proof of our main result. To prove the result, we construct

a tunneling with m2 arcs. This proves cp(K) ≥ m2 ≥ t(K). Then, by applying Lemma

3.2.2, we get that cp(K) = m0 + m1 + m2 + m3 ≥ 1 + (t(K) + 1) + t(K) + 0 = 2t(K) + 2 as

desired.

3.3.1 Construction of the tunneling

We will apply the Morse Rearrangement Lemma to allow us to make some ad-

ditional convenient assumptions about Φ without losing generality. A proof can be

found in Nicolaescu, Chapter 2.4 [88]. The theorem states we can find a smooth

Φ̂ : S 3 − K → R satisfying the following properties:

• Φ and Φ̂ share the same critical points, and each critical point has the same index.

• Suppose p and q are distinct critical points. Then Φ̂(p) , Φ̂(q). If Φ̂(p) < Φ̂(q),

then λ(p) ≤ λ(q).

• Inside of a neighborhood of each critical point, the gradient flows for Φ and Φ̂

are identical.

• If γ is an integral curve of ∇Φ̂, then Φ(γ(t)) is strictly increasing in t. A vector

field with this property is called gradient-like with respect to Φ.

• For p, q ∈ Crit(Φ̂), WS (p) and WU(q) intersect transversely. Morse functions with

this property are called Morse-Smale.

45


This theorem allows us to perturb the values of the critical points so we can reorder

them ascending by index without affecting the topological data encoded by the origi-

nal potential. At this point, we are not necessarily working with the physical potential

whose formula is given in (3.1), but for simplicity we will still refer to the perturbation

as Φ.

Our rearrangement restricts the limiting behavior of trajectories. Let γ : (−∞,∞)→

S 3 − K be a trajectory for Φ. When t → −∞, Φ(γ(t)) strictly decreases, but it is bounded

below by 0, by assumption. Should lim
t→−∞
Φ(γ(t)) = 0, then lim

t→−∞
γ(t) = ∞, the point of

infinity on S 3, because it is the only point in the knot complement with zero potential.

Otherwise, lim
t→−∞

γ(t) is a critical point. Similarly, we can deduce that for t → ∞, either

γ(t) tends to a critical point or K, since Φ(γ(t)) is strictly increasing.

Should both ends of γ be critical points, then we know the index of the critical

point at t = −∞ is less than or equal to that of the critical point at t = ∞.

Consider the critical points of index 2. For p2
1, . . . , p2

m2
∈ Crit2(Φ), let Γi = WU(p2

i ).

See Fig. 3.1 for a schematic of the trefoil case. Notice each Γi is a union of two tra-

jectories leaving p2
i , since the unstable manifold has dimension 1. Since Φ strictly

increases along trajectories, we have that the trajectories will either tend to K or to

another critical point of index 2 as t → ∞. However, should either end of Γi tend to a

critical point q ∈ Crit2(Φ), then the corresponding trajectory will be a submanifold of

WS (q). However, WS (q) and Γi are 2 and 1-submanifolds respectively, and for them to

intersect transversely as per our assumption, their intersection cannot be more than

0 dimensions. Therefore, we conclude both ends of Γi eventually reach the tubular

neighborhood of the knot. Our tunneling is only concerned with the arc outside of the

tubular neighborhood, so we can assume Γi is defined only on a compact interval.

46


Figure 3.1: The arcs we add to K are the unstable manifolds associated to critical points
of index 2. Note that this diagram does not necessarily depict the specific situation
accurately for the trefoil.

Now consider the critical points of index 1. For p1
1, . . . , p1

m1
∈ Crit1(Φ), let Θ j =

WS (p1
j). Analogous to before, each Θ j is a union of two trajectories tending to p1

j . Sim-

ilar reasoning will tell us that both ends tend to the point at infinity. Indeed, each Θ j is

a union of two trajectories tending towards a critical point p1
j of index 1. As t → −∞,

Φ will strictly decrease along these trajectories, so we know that the negative infinite

limits of these trajectories must either be another critical point of index 1, or the point

at infinity. By the transversality assumption of the stable and unstable manifolds, we

cannot have that the endpoints ofΘ j are critical points of index 1. Therefore, both ends

are at the point at infinity. We may view the union of all the Θ j arcs as a wedge sum of

circles at the point of infinity, which we denote
m1∨
j=1
Θ j.

Notice
m1∨
j=1
Θ j is homotopy equivalent to a handlebody. Using a standard maneuver

from differential topology, we will flow along the (negative) gradient to perform a

deformation retraction from S 3 − (K ∪ Γ1 ∪ · · · ∪ Γm2) to
m1∨
j=1
Θ j. We will opt to work with

smooth tubular neighborhoods of both of these spaces, but there is a crucial technical

lemma we must prove before proceeding.

47


Figure 3.2: A sketch depicting the tubular neighborhood around a regular point of
Γi in specially chosen local coordinates which makes the gradient parallel. This is a
projection to the xz-plane.

3.3.2 Constructing a smooth boundary around the tunneling

In this subsection, we prove the following lemma.

Lemma 3.3.1. There are tubular neighborhoods of both K ∪ Γ1 ∪ · · · ∪ Γm2 and
m1∨
j=1
Θ j with

smooth boundary such that the gradient vector field points inwards and outwards respectively.

Proof. It is a standard fact from electrostatics that neighborhoods of K and the point

at infinity exist such that the gradient points inwards and outwards respectively. This

proof will make use of the tubular flow lemma, which states for every regular point of

S 3 − K, there are local coordinates such that the gradient flow takes the constant form

∇(x0,y0,z0)Φ = (1, 0, 0) = ∂
∂x . Take local coordinates centered at a regular point of Γi so that

the portion of the arc inside our coordinate chart is mapped to the path of unit speed

along the x-axis, γ(t) = (t, 0, 0).

48


Figure 3.3: A sketch depicting the tubular neighborhood around a critical point on Γi,
in specially chosen local coordinates such that Φ is a quadratic with signature (1, 2).
This is a projection to the xz-plane.

Around a segment of the x-axis, we can choose our tubular neighborhood of Γi to be

the interior of the slanted tube defined by y2+ z2 = ( 1
2 x+1)2, where we possibly restrict

our local coordinates so −1 < x < 1. See Fig. 3.2. Notice the gradient points inwards

from the boundary. If we expand around a regular point on a Θ j arc, we can reflect

the tube in the x direction so the gradient points outwards. Note that the properties of

a vector field pointing inwards and outwards from a boundary are invariant under a

change of coordinate charts in an orientation-preserving atlas.

Now suppose we want to take local coordinates around a critical point p2
i of

Γi, which has index 2. By Lemma 2.2 of [82], we can find local coordinates cen-

tered at the critical point such that Γi corresponds to the x-axis, and Φ takes the

form Φ(x0, y0, z0) = 1
2 (x2

0 − y2
0 − z2

0) + c, for some constant c. Thus, ∇Φ takes the form

∇Φ(x0, y0, z0) = (x0,−y0,−z0). Consider the tube around the x-axis defined by y2+ z2 = 1.

49


See Fig. 3.2. Using an abuse of notation, we will refer to this tube as ∂Γi. We will show

the gradient points inside the tube, as Figure 3.3 shows. At a point (x0, y0, z0) ∈ ∂Γi, the

tangent plane is spanned by (1, 0, 0) and (0,−z0, y0). As the point varies on the tube, this

oriented basis of the tangent plane varies smoothly, thus defining an orientation of the

tube. The triple of tangent vectors {(1, 0, 0), (0,−z0, y0),∇Φ(x0, y0, z0)} therefore defines

a smooth choice of orientation for the three dimensional ambient tangent spaces sur-

rounding the tube. Therefore, if the gradient points inside the tube at one point of the

tube, it does so throughout the whole tube. For example, at the point p = (0, 1, 0) ∈ ∂Γi,

∇Φ(p) = (0,−1, 0) clearly points inside the tube. The existence of a surface surrounding

the arc with which the gradient points inwards is a significant topological obstruction

to proving our main theorem. The construction would not be possible if the critical

point had index 1, or if we were asked to place a tube around another axis. The case

for a critical point of Θi is analogous, and we get the result that the gradient points

outwards.

By compactness of the arcs Γi and Θ j we only need to construct finitely many tubes

around arc segments to encapsulate the entire arc. When the tubes cover the same

part of the arc, we may have to shrink the radius of one of the tubes so the boundaries

will intersect, but the resulting tube will still have the gradient pointing in the correct

direction. By taking the final boundary to be the points of minimal radial distance

to the arc, we obtain a connected boundary to the tubular neighborhood that is only

piecewise smooth. Likewise, the tubular neighborhoods around the Γi arcs intersect

the tube around K, and the gradient points inwards on the boundary of the union.

This construction is sufficient to prove the lemma. To get a smooth boundary out of

the piecewise smooth boundary, one could either make a density argument in a space

of manifolds [88], justify why the theorems work in the piecewise smooth case [82],

or use mollifiers to smooth out the kinks [116]. For the sake of space, we omit the

50


technical details.

For reasons that will be clear when we perform the deformation retraction, we will

want to include a ball around the point at infinity in the tubular neighborhood. In the

standard R3 coordinates, this is the complement of a large open ball. We can assume

this neighborhood around ∞ contains Φ−1([0, ε]) for some ε > 0. On the boundary

sphere, we can again use mollifiers to connect the tube smoothly whilst preserving

their orientations against the gradient flow. □

3.3.3 The deformation retraction to a handlebody

The final step is to use the flow of E to perform the deformation retraction. We will

use the closed tubular neighborhoods described in Lemma 3.3.1 as an alternative to

just the knot with arcs attached and the wedge of circles. Let A be the aforementioned

tubular neighborhood of the knot with the Γi arcs attached, and let B denote the afore-

mentioned tubular neighborhood of the Θ j arcs connected to a ball around the point

at infinity.

Let E(x, t) be the flow of the negative gradient. The point E(x, t) refers to the loca-

tion of the path at time t starting from the unique integral curve starting at x. First, we

prove every point in S 3 − A will eventually flow to B. Suppose x ∈ S 3 − A. There are

three possibilities for the limit of the negative gradient flow of x: it could flow to K,

it could flow to a critical point of index 1 or 2, or it could flow to the point at infinity.

Since the negative gradient points outwards from the boundary of A, x cannot flow to

K or a critical point of index 2. Therefore, x either flows to a critical point of index 1, or

to the point at infinity, which means that x eventually flows to B. Furthermore, since

the negative gradient points into B, once x enters B, it will never leave. This fact still

holds even in the vacuous case where m2 = 0 and therefore A = K.

51


By smoothness of the boundary of B, the function which assigns each x ∈ S 3 − A

the minimum time t such that E(x, t) ∈ B is smooth. Call this function C(x). Notice

C assigns 0 to each point already in B. By compactness of the domain, the func-

tion reaches a maximum value Cmax. Define a homotopy H on (S 3 − A) × [0,∞) by

H(x, t) = E(x,min(t,C(x))), which we can see is continuous. Also notice that for x ∈ B,

H(x, t) = x for all t. Running this homotopy on the time interval [0,Cmax] completes the

deformation retraction. This completes the proof of our main theorem.

52


CHAPTER 4

FURTHER TOPOLOGICAL APPROACHES TO KNOTTED ELECTRIC CHARGE

DISTRIBUTIONS

This chapter orignally appeared as “Topological approaches to knotted electric

charge distributions” in Partial Differential Equations and Applications Vol. 4, Iss. 2

[68].

4.1 Introduction

Imagine a closed knotted wire in space fixed in place, with uniform electric charge

distributed on it. Our novel problem of interest is to understand the properties of the

knot’s electric field, a topic with applications to the study of knotted DNA molecules

in viruses [95], the design of electrical circuits [35], and broader questions in physical

knot theory [55]. In our work, the knot is fixed in space. By contrast, others have

studied how knotted curves deform, as discussed in the extensive literature on the

energy methods and flows of knots [12, 34, 89, 137].

The electric field around a charged knot is the negative gradient of a potential

function, so the zeroes of the electric field correspond to critical points of the potential.

Our primary objects of study are the structurally stable critical points with respect to

small perturbations of the charge distribution. Physically, these are the points where

all of the electric forces cancel. Thus some basic questions arise: How many critical

points are there? What are their local behaviors? And how does the critical set relate

to topological properties of the knot?

In this paper, we present some results which address these broad questions

through the use of techniques from geometric topology, and offer some rigorous

53


proofs for some of the results communicated in [71]. Whilst the methods we apply

are standard in low-dimensional topology, research in electrostatics and physical knot

theory rarely implement these methods. We use a topological counting argument to

show that the size of a knot’s structurally stable critical set sharpens from a previously

proven lower bound based on the knot’s tunnel number [65]. Next, we calculate the

effect of the Morse surgery on their genera. A critical point of Morse index 1 corre-

sponds to an increase in genus as we raise the level of the potential, whilst a critical

point of index 2 corresponds to a decrease in genus. The tuple of equipotential sur-

face genera is a Scharlemann-Thompson handlebody width decomposition of the knot

complement [112], but our tuple differs from the classical thin decomposition in that

we exclusively consider those arising from potentials, and we do not lexicographically

reorder the tuples.

Let K ⊆ R3 ⊆ S 3 be a smooth, tame knot parametrized by the curve r(t), t ∈ [0, 2π]

with r(0) = r(2π). We will take the convention that S 3 is the union of R3 and a single

compactifying point at infinity. Suppose K is endowed with a uniform charge distri-

bution. With a choice of units, the electric potential between a point k ∈ K and a point

charge at x at a distance R from k is proportional to R−1 by Coulomb’s Law. It therefore

makes sense to define the electric potential Φ : S 3 − K → R, on the complement of K

by the line integral

Φ(x) =
∫

k∈K

dk
|x − k|

=

∫ 2π

0

|r′(t)|
|x − r(t)|

dt, x ∈ R3 − K. (4.1)

We set Φ(∞) = 0 to ensure smoothness. By differentiating under the integral sign

with respect to x, we can see the electric potential is smooth and harmonic. The elec-

tric field is defined by E = −∇Φ. The critical points of Φ represent equilibrium points

where a charged test particle at rest will continue to experience no electric force from

the charge distribution. Some conventions define the potential to be negative, in order

54


for the electric field to point towards the knot, but it is more convenient for our pur-

poses to work with a nonnegative potential function. Typically charged wires are pre-

sented in physics literature as the integral over a solid, thin torus, where our idealized

uniform charge distribution over a curve is approximated by a nonuniform charge dis-

tribution over a solid, which depends greatly on the geometry of the parametrization.

For background on the physical interpretation of the potential, see [39]. However, our

work is concerned with critical points which occur far away from the knot, so this

distinction does not alter our results.

The electric potential is difficult to analyze in general, as integral formulas can

be highly nonlinear and nonsymmetric, but the potential and electric field inte-

grals can be numerically approximated via Gaussian quadrature, and the solutions

to E(x) = 0 can be determined by the multivariable Newton method. Initial nu-

merical approximations of the critical sets for various knots were implemented by

Townsend and Lipton [64], showing that there are isotopic parametrizations of the

unknots which have critical sets of differing size. In particular, the trivial unknot em-

bedding r(t) = (cos(t), sin(t), 0) has a single critical point in the origin, whilst the (5, 1)

unknot with parametrization r(t) = (cos(t + 2) cos(5t), cos(t + 2) sin(5t),− sin(t)) has three

critical points, as depicted in Figure 1.

Figure 4.1: A top-down projection of two knots isotopic to the unknot with differing
critical sets.

It is unsurprising that the critical set is not invariant under isotopy, given the de-

pendence ofΦ on r, but a 2021 result of the author [65] proved there is indeed a restric-

tion on the critical set based on a knot invariant called the tunnel number. Let cp(K)

55


denote the smallest size of structurally stable critical sets of K over all parametriza-

tions, which includes the critical point at infinity in the count. Let t(K) denote the

tunnel number, which is the smallest number of arcs one must attach to K such that

the complement is a handlebody. By Waldhausen’s theorem, t(K) is equivalently the

smallest number of arcs one must attach such that a tubular neighborhood of the union

is isotopic to an unknotted embedding of a handlebody. For example, it is trivial

to see t(unknot) = 0, whilst a geometric argument depicted in Fig. 2 of [84] shows

t(trefoil) = 1. The aforementioned result of the author is stated as follows.

Theorem 4.1.1. For all knots K, cp(K) ≥ 2t(K) + 2.

Theorem 4.1.1 was proven using techniques in geometric topology, Morse theory,

and stable manifold theory. The proof does not apply Morse reconstruction to Φ itself,

but rather, it applies the Morse Rearrangement Lemma (see Lemma 2.4.12 of [88]) to

Φ, which allows us to alter the global, long-term dynamical behavior of the gradient to

conform to certain convenient properties, whilst leaving the local behavior around the

critical points unchanged. Critical sets of harmonic functions can be pathological in

general, but the conditions we impose avoid these edge cases, and so the critical sets

under consideration will always be finite [114] and lie in the convex hull of the charge

distribution [126], which we assume is bounded. The inclusion of the structurally sta-

ble requirement in our definition of cp(K) is important, because as we will discuss

in Section 2, generic knot isotopies (under the C∞ topology) can induce bifurcations

which create or annihilate pairs of critical points, with a structurally unstable degen-

erate critical point of high multiplicity occurring during the transition. Demanding

that the potential is Morse avoids undercounting the critical points.

In Section 3, we prove a partial sharpening of Theorem 4.1.1. The tunnel number is

a poor measure of a knot’s complexity. For instance, all torus knots have tunnel num-

56


ber 1, even though they can have arbitrarily high crossing number. Likewise, prime

knots with up to 10 crossings can only have tunnel numbers 1 or 2 [85]. However, we

can show the bound sharpens when the knot is sufficiently close to a planar curve,

defined by the notion of a δ-lifting which is given later on. We then show that given a

planar knot diagram with c crossings, there is a sufficiently small δ > 0 such that the

size of the critical set for its δ-liftings is at least 2c+2. This result is proven by applying

the Poincare-Hopf Index Theorem to the electric field of the projected knot to prove

the existence of a critical point in each connected component in the projection’s planar

complement, of which there are c+ 2, including the critical point at infinity which cor-

responds to the unbounded component. We then use the Morse inequalities to show

that lifting each crossing from the projection necessitates the creation of another criti-

cal point due to the change in the topology of the knot complement, which yields our

final count of 2c + 2.

In Section 4, we prove a result describing the topology of the equipotential surfaces

in terms of the critical set. The critical values of Φ partition the positive real numbers

into finitely many intervals, and a key idea from Morse theory states that the level

sets of two values from the same interval are homeomorphic. This leads us to the

central object of the section: the Morse code of a knot, which is a tuple of numbers

listing the genera of the equipotential surfaces corresponding to each of the intervals.

Another key idea from Morse theory states the passage from one equipotential surface

to another is determined, up to homotopy equivalence, by the attachment of a disc of

appropriate dimension. Even though the attachment maps can be quite pathological

in general, we prove that the potential yields particularly simple surgeries, showing

that a critical point of index 1 corresponds to an increase in the subsequent term of the

Morse code, whilst a critical point of index 2 corresponds to a decrease.

57


4.2 Bifurcations through degenerate critical points

A strategy for finding more precise statements than Theorem 3.1.1 is to bifurcate

the electric potential’s critical set via knot isotopies into a set with a sharper lower

bound. One such class of isotopies proposed by Lipton-Strogatz-Townsend is the

“flattenings,” where one of the coordinate functions is scaled by the parameter s, as s

ranges from 1 to 0. The final state is a self-intersecting planar curve, but we opt to con-

tinue referring to this process as an ambient isotopy in an abuse of notation. In [71], an

article written for a general audience of physicists, the authors implement numerical

simulations of charged knots yielding critical sets whose sizes exceed the proven the-

oretical lower bound. In this section, we discuss the general fashion with which these

critical sets transition from one to another as the electric potential bifurcates. Essen-

tially, as the bifurcation is through harmonic functions, this constrains the transitions

to be the creation or destruction of pairs, modulo structurally unstable symmetries

arising from the particular parametrization and isotopy.

A knot isotopy induces a path in C∞(M) via the corresponding potentials, but the

path could pass through non-Morse functions. The proofs of our theorems in the

subsequent sections make liberal use of the Morse property, and break down when

Φ is not Morse. However, using Cerf theory, we can extract a generic property of

knot isotopies in which the path through non-Morse potentials occurs only finitely

many times, and their degenerate critical points can be classified based on partial

derivatives of the isotopy. We will see that these degenerate critical points correspond

to saddle-node bifurcations. It should be emphasized that this result does not describe

the dynamics of a moving charged knot, as we do not take into account the induced

magnetic forces from the moving charges. Rather, we are describing a bifurcation

among a smooth family of possible fixed, rigid knot configurations.

58


We now state a part of the Cerf Structure Theorem in terms of the homotopies of

Morse functions. There is a more general structure theorem given by a stratification

[88, 33] of the class of isotopies via the coordinate expansion of general degenerate

critical points. A full exposition can be found in Ch. 23 of Freed [33].

Theorem 4.2.1. Let F̃ : M × [0, 1] → R be a one parameter family of smooth functions on

an n−manifold M indexed by the variable s. We can replace F̃ with a substitute function

(which we still call F̃) in a dense subspace of C∞(M × [0, 1]) which satisfies the following.

Each Fs = F̃(−, s) is Morse except for finitely many s where there could exist degenerate

critical points. For each such s0, and for each degenerate critical point p of Fs0 , there exist

local coordinates (x1, . . . , xn, s) of M × [0, 1] centered at (p, s0) such that F̃(x1, . . . , xn, s) =

x3
1 + ε1sx1 + ε2x2

2 + · · · + εnx2
n + C, where each εi ∈ {±1}, whose values depend on the partial

derivatives of F̃, and C is the degenerate critical value.

ε1 ε2 ε3 Bifurcation for increasing s
1 1 1 Impossible
1 1 −1 Destruction of critical pair
1 −1 1 Destruction of critical pair
1 −1 −1 Impossible
−1 1 1 Impossible
−1 1 −1 Creation of critical pair
−1 −1 1 Creation of critical pair
−1 −1 −1 Impossible

Figure 4.2: The generic classification of bifurcations of critical points of the potential
induced by knot isotopy.

We will apply this theorem to the one parameter family of electric potentials in-

duced by a knot isotopy. Suppose K0 and K1 are two smooth isotopic knots where the

isotopy is parametrized by r : [0, 2π] × [0, 1] → R3. For fixed s0 ∈ [0, 1], we will say

rs0 = r(−, s0) is the parametrization of the knot Ks0 , whose electric potential is Φs0 . We

will now apply Theorem 4.2.1 to the one parameter familyΦs, and if necessary, replace

59


Φs with another nearby generic family so that the property mentioned in the theorem

holds.

We will label the coordinates from Theorem 4.2.1 as (x, y, z, s) so we have that the

local coordinate formula of Φ at a specific s0 where Φs0 has a degenerate critical point

is Φ(x, y, z, s) = x3 + ε1sx + ε2y2 + ε3z2 +C.

Remark 4.2