A crowd of clocks that falls into step

Two chaotic orbits made to agree coupled identical copies of one chaotic map and found a sharp threshold. Below a certain strength of coupling the copies wander apart, above it they agree to every digit, and the strength needed was read off the Lyapunov exponent — the rate at which the copies would have separated if left alone. It ended on the case it could not handle. Real populations are not identical. Fireflies flash at slightly different natural rates, the pacemaker cells of a heart have different intrinsic periods, and the pedestrians on a footbridge walk at different paces. For units that differ, agreement to every digit is impossible: two clocks running at different rates cannot show the same time for ever unless something holds one back and pushes the other on.

Yoshiki Kuramoto wrote down the simplest model of such a population in 1975, and it has become the standard example of how collective order switches on. Each of $N$ oscillators is described by a single phase $\theta_i$, an angle going round the circle. Left alone, oscillator $i$ advances at its own natural frequency $\omega_i$. Coupled, each is pulled towards every other with a strength proportional to the sine of their phase difference:

$$\frac{d\theta_i}{dt} = \omega_i + \frac{K}{N}\sum_{j=1}^{N} \sin(\theta_j - \theta_i).$$

The natural frequencies are drawn from a spread — a probability density $g$ that is symmetric and peaked at nought, which can always be arranged by watching from a frame that rotates at the mean frequency. The question is what happens as the coupling strength $K$ increases from nothing. This essay answers it with a simulation of a few thousand oscillators, and finds Kuramoto’s answer on the screen to three decimal places.

One number that measures agreement #

Kuramoto’s first move turns a sum over all pairs into something each oscillator can feel on its own. Draw each oscillator’s phase as an arrow of length one pointing at angle $\theta_j$, and take the average of all the arrows:

$$r\,e^{i\psi} = \frac{1}{N}\sum_{j=1}^{N} e^{i\theta_j}.$$

The length $r$ of that average is the order parameter. If the phases are spread evenly round the circle, the arrows cancel and $r$ is nearly nought; if every oscillator has the same phase, $r$ is exactly one. The angle $\psi$ is the phase of the crowd. Taking imaginary parts of the same identity after rotating by $\theta_i$ rewrites the equation of motion as

$$\frac{d\theta_i}{dt} = \omega_i + K r \sin(\psi - \theta_i).$$

Every oscillator now sees only one thing: the crowd’s average phase, pulling it with an effective strength $Kr$. This is a mean-field description, and it makes the feedback plain. The more ordered the crowd is, the larger $r$, the harder each oscillator is pulled towards $\psi$, and the more ordered the crowd becomes. Whether that loop can start from nothing is the whole question. When $r$ is small, the pull $Kr$ is small too, and an oscillator whose natural frequency is far from the centre will slip past the crowd’s phase over and over without ever being caught.

The hero figure shows both outcomes. The frequencies there are drawn from the Lorentzian density of width $\gamma = \tfrac12$, the bell with heavy tails that an average that never settles met as the Cauchy distribution — the one whose average of samples is as spread as one sample. It is the spread Kuramoto’s theory is exact for, which is why it is used in every figure but the last. At $K = 0.5$ the average arrow has length 0.07; at $K = 2.5$ it has length 0.78, and the warm points are the oscillators that have given up their own frequency to run with the crowd. To keep the experiment reproducible, the natural frequencies in every run are the evenly spaced quantiles of the chosen density rather than random draws, so a population of $N$ is a deterministic sample of the distribution, and only the starting phases are random.

The threshold Kuramoto computed #

Suppose the crowd has settled into a steady state with order $r$, rotating at the mean frequency, so that $\psi$ can be taken to be nought. Each oscillator obeys $d\theta/dt = \omega - Kr\sin\theta$. If $|\omega| \le Kr$ this has a fixed point, $\sin\theta = \omega/(Kr)$: the oscillator is locked, sitting at a fixed angle to the crowd and turning with it. If $|\omega| > Kr$ there is no fixed point, and the oscillator drifts round and round the circle, slowed down where the pull opposes it and hurried along where it helps.

For the steady state to be consistent, the arrows of all these oscillators must average back to the same $r$ that produced them. The drifting oscillators contribute nothing. Each spends more time where it moves slowly, but the oscillator with frequency $-\omega$ spends its time at the mirror-image angles, and by the symmetry of $g$ their contributions cancel. The locked oscillators contribute $\cos\theta$ each, and substituting $\omega = Kr\sin\theta$ turns the sum into an integral:

$$r = Kr\int_{-\pi/2}^{\pi/2}\cos^2\theta\; g(Kr\sin\theta)\,d\theta.$$

The equation always has the solution $r = 0$ — the incoherent state. A second solution exists when, after dividing by $r$, the right-hand side can equal one. As $r$ shrinks towards nought the integral tends to $g(0)\,\pi/2$, so a nonzero $r$ first appears when

$$K_c = \frac{2}{\pi\,g(0)}.$$

That is Kuramoto’s threshold. It depends on the density of frequencies at the centre and on nothing else — not the tails, not the width in any other sense, not the number of oscillators. For the Lorentzian $g(0) = 1/(\pi\gamma)$, so $K_c = 2\gamma$, which is $1$ for the width used here. And for the Lorentzian the integral can be done in closed form, giving the whole curve: above the threshold, $r = \sqrt{1 - K_c/K}$.

The square root, measured #

The calculation above is for infinitely many oscillators, in a state assumed steady. A simulation has neither luxury. It has a finite number of oscillators, starts them at random phases, and lets them run. The figure shows the result for 400 and for 2,000 oscillators at fourteen couplings, each run integrated for 200 time units to let the transient die and averaged over the next 100.

The agreement is close enough to be startling for a formula derived under such idealised assumptions. At $K = 1.1$, just past the threshold, the curve predicts $\sqrt{1 - 1/1.1} = 0.3015$, and both populations measure 0.302. At $K = 2$ the prediction is $\sqrt{1/2} = 0.7071$ and the measurement for 2,000 oscillators is 0.707. At $K = 4$ the prediction is 0.866 and the measurement is 0.865. The corner at $K_c = 1$ is already sharp at these sizes: at $K = 1$ itself the order is 0.06, and one step of a tenth later it is 0.30.

Two features of the shape are worth naming. The first is that the onset is continuous. At the threshold $r$ starts from nought rather than jumping, so this is what physicists call a second-order or continuous phase transition, of the kind that turns a magnet’s scattered atomic spins into an aligned field as a metal cools. The second is the square root itself. Near the threshold $r$ is proportional to $\sqrt{K - K_c}$, so the curve leaves the axis vertically; a tiny increase in coupling past $K_c$ buys a disproportionate amount of order. This exponent of one half is the signature of a mean-field transition, and it is the same square root that governs the drift speed past a saddle-node, met again below, because it comes from the same place: a quantity that must balance its own square.

Below the threshold: not quite nothing #

The curve says $r = 0$ below the threshold. The simulation does not quite agree, and the disagreement is informative rather than a defect. Below $K_c$ the order for 400 oscillators sits at about 0.045, and for 2,000 at about 0.025. A finite population of arrows pointing in scattered directions does not cancel exactly.

The size of the residue follows the oldest law in probability. Adding $N$ arrows in independent random directions is a random walk in the plane, and the bell curve built from coin flips is the reason its length grows like $\sqrt N$. Dividing by $N$ to take the average leaves something of size $1/\sqrt N$. The figure measures the residual order at half the threshold coupling for populations from fifty to three thousand two hundred, and the fitted slope on logarithmic axes is $-0.45$, against the $-\tfrac12$ the random walk predicts. The small shortfall comes from the weak coupling, which is not nothing: at $K = 0.5$ the oscillators still nudge one another, and correlated nudges make the imbalance a little larger than independent arrows would. The average settles and the wobble does not drew the same two scalings in one picture: the average goes to its limit, and the fluctuation around it shrinks only as the square root.

This is the sense in which the sharp threshold is a statement about a limit. With fifty oscillators the residual order is 0.12, and a reader looking at the phases could not reliably say whether the population was below the threshold or a little above it. With three thousand two hundred the residue is 0.018, and the corner is unmistakable. Every real population is finite, so every real transition is slightly rounded. The theory describes what the rounding converges to.

How long order takes to grow #

A steady state says nothing about how it is reached. The next figure follows the order parameter in time for 2,000 oscillators started at random phases, at three couplings: below the threshold, a quarter above it, and two and a half times it.

Below the threshold nothing grows; the order flickers at the level of the finite-size noise for the whole run. Well above the threshold, at $K = 2.5$, the order climbs to its steady value of 0.77 within about five time units. Just above it, at $K = 1.25$, the climb takes more than twenty, and the order settles at 0.45 against the predicted $\sqrt{1 - 1/1.25} = 0.447$.

The slowness near the threshold has a precise source. The incoherent state is a solution at every coupling; what changes at $K_c$ is whether it is stable. Strogatz and Mirollo showed in 1991 that for the Lorentzian a small ripple of order in the incoherent state grows at the rate $K/2 - \gamma$, which is negative below the threshold, nought at it, and $0.125$ at $K = 1.25$. A growth rate of one eighth means the ripple takes eight time units to grow by a factor of $e$, and the random starting order of $1/\sqrt{2000}$ needs several such factors to reach 0.45. The figure’s delay is that arithmetic, and it is the dynamical face of the threshold: near a continuous transition the system is only barely persuaded to leave its old state, and it leaves slowly. Physicists call this critical slowing down, and it is the same phenomenon as the long bottleneck that the window that opens with a stutter found just before a periodic window opens.

Who locks and who drifts #

The self-consistency argument split the population cleanly. Oscillators whose natural frequency lies inside the band $|\omega| \le Kr$ lock to the crowd; those outside it drift. The simulation can check this oscillator by oscillator, by measuring each one’s average frequency over a long run.

At $K = 2$ the order is $r = 0.707$, so $Kr = 1.414$. Every oscillator with natural frequency inside $\pm 1.414$ has average frequency nought to within measurement: all of them, without exception, have abandoned their own rhythm and run with the crowd. Outside the band each oscillator keeps a frequency of its own, but a reduced one, and it lies exactly on the dashed curve $\pm\sqrt{\omega^2 - (Kr)^2}$. That curve is the frequency of an angle obeying $d\phi/dt = \omega - Kr\sin\phi$, which spends a long time crawling past the place where the pull nearly cancels its natural speed. It is the bottleneck that how a lock comes apart found at the edge of an Arnold tongue, where a locked orbit disappears in a saddle-node bifurcation and the drift that replaces it starts with zero average speed and grows as a square root. Every drifting oscillator here is sitting outside its own tongue, and the crowd’s pull $Kr$ is the width of the tongue.

The locked share can be predicted as well. It is the fraction of the Lorentzian within $Kr$ of nought, which is $(2/\pi)\arctan(Kr/\gamma)$ — 78.4% at these numbers. The simulation counts 784 of 1,000 oscillators locked, 78.4%. The drifters are the remaining fifth, spread over the tails that the heavy-tailed Lorentzian makes long; a tighter spread would lock a larger share at the same coupling.

A different spread, the same threshold #

The Lorentzian is chosen because it makes everything exact. It is not chosen because the threshold depends on it. Kuramoto’s formula $K_c = 2/(\pi g(0))$ uses only the height of the density at the centre, and the last figure tests that by changing the spread to the ordinary normal distribution with standard deviation one.

For the normal density $g(0) = 1/\sqrt{2\pi}$, so the threshold is $2\sqrt{2\pi}/\pi$, about 1.596. The simulation finds the order flat at about 0.03 up to $K = 1.5$, then 0.11 at $K = 1.6$, then 0.43 at $K = 1.7$. The transition sits where the formula put it, and what happens after it is different. The normal density has a sharper peak relative to its tails than the Lorentzian, so a small increase in the pull beyond the threshold captures a large band of oscillators at once. Near the threshold the order still grows as a square root of $K - K_c$, but with a coefficient set by how sharply the bell is curved at its peak, and for the normal spread that coefficient is large — the curve is nearly a cliff. Further from the threshold no closed form is known, and the curve is found by solving the self-consistency integral numerically.

The comparison makes the meaning of the threshold clear. The oscillators that matter at the onset are the ones whose natural frequency is already almost the crowd’s; they are the first to be caught, and the density of them is $g(0)$. A spread with more of them in the middle starts its cluster at a weaker pull. The tails decide how far the cluster can spread afterwards, and so the shape of the curve above the threshold, but they have no say in where it begins.

From Huygens’ clocks to the Millennium Bridge #

Christiaan Huygens noticed in 1665 that two pendulum clocks hung from the same beam came to swing in exact opposition, and stayed so whatever he did to disturb them. He had found mutual synchronisation through a shared support, and the phenomenon was rediscovered in organ pipes, electrical generators and radio circuits over the next three centuries. The population version began with the biologist Arthur Winfree, who in 1967 proposed that large groups of biological oscillators with different natural rates would show a threshold for collective synchrony, and found it in simulation.

Kuramoto simplified Winfree’s model to the sine coupling above in 1975 and solved for the threshold by the self-consistency argument, in a few pages of a conference proceedings. The argument was famously incomplete. It assumes a steady state and shows only that one exists; it says nothing about whether the incoherent state is unstable below and above the threshold, or whether the partly locked state is what the system actually approaches. Steven Strogatz’s survey of 2000, From Kuramoto to Crawford, tells the story of a quarter-century of attempts to fill those gaps, including the discovery that the incoherent state is neutrally stable below the threshold in a linearised sense, with perturbations decaying only through a mechanism borrowed from plasma physics. In 2008 Edward Ott and Thomas Antonsen found that for the Lorentzian the whole infinite population obeys one ordinary differential equation for the order parameter, which made the square-root curve and the growth rate exact consequences rather than steady-state guesses. Hayato Chiba’s proof of the conjectured bifurcation for general symmetric spreads followed in the 2010s.

The applications came alongside. In June 2000 the Millennium Bridge in London swayed sideways under its opening-day crowd, and the analysis that followed found a Kuramoto-like mechanism: each walker adjusted their step to the swaying, which fed the sway, and above a critical number of pedestrians the crowd locked into step. The bridge was fitted with dampers and has not swayed since. The same equations are now applied to power grids, whose generators must stay in phase with one another, to the circadian pacemaker in the brain, and to arrays of lasers and Josephson junctions.

Still open: a finite crowd at the edge #

The infinite-population theory is now complete for the Lorentzian and largely complete in general. What remains open is the finite population near the threshold, which is the only kind of population that exists. The finite-size figure above measured the residue well below the threshold, where it is simply the random walk’s $1/\sqrt N$. At the threshold itself the fluctuations are larger and decay more slowly with $N$, and how much more slowly is not settled. Numerical studies of the exponent disagree depending on whether the natural frequencies are drawn at random or spaced evenly as they are here — random draws add a second source of finite-size noise, the irregular density of frequencies near the centre — and none of the measured exponents has been derived.

A second question is where synchrony goes when the network is not complete. Kuramoto’s oscillators each feel all the others equally, which is what makes the mean field exact. In a ring, a grid or a social network, each oscillator feels only its neighbours, and the threshold depends on the network’s shape in ways that have been computed for a few families and are open in general. In one dimension with nearest-neighbour coupling, it is known that complete locking of an infinite chain needs a coupling growing with its length, so a long chain never fully synchronises — the same obstruction that two chaotic orbits made to agree found for rings of chaotic maps, arriving now from the side of oscillators that differ rather than units that are chaotic.

Agreement without sameness #

The threshold for two chaotic orbits was set by chaos: identical units, each amplifying its own differences at a rate the Lyapunov exponent measured, and a coupling that had to beat that rate. Here the units are as tame as a unit can be — each one alone simply goes round the circle at a steady speed — and the threshold is set by diversity. The quantity the coupling has to beat is the density of natural frequencies near the centre, and what it wins is not agreement to every digit but a cluster that runs together while the outliers continue on their own.

The two thresholds share a structure. Both compare a pull towards agreement against a tendency to come apart, and both switch on when one exceeds the other. And the self-consistency loop that makes Kuramoto’s transition work — order produces a stronger pull, and a stronger pull produces order — is the reason a few thousand pendulum-like equations, each with its own frequency, behave as one object with a single number describing it. Below the threshold that number is nought up to the noise of a sine seen as a turning circle averaged a thousand times; above it the number grows as a square root, and the simulation finds the root to the third decimal.