The Kuramoto Model

The Kuramoto model is the canonical framework for studying synchronization in populations of coupled oscillators. Derived from weakly coupled nonlinear oscillators via averaging theory, it exhibits a continuous phase transition from incoherence to collective synchrony—a phenomenon with deep connections to urban rhythms, from traffic signal coordination to commuter behaviour.

1. Derivation from Weakly Coupled Oscillators

Consider $N$ nonlinear oscillators, each described by a limit cycle with natural frequency $\omega_i$. In the uncoupled case, each oscillator has a well-defined phase $\theta_i(t)$ satisfying $\dot{\theta}_i = \omega_i$.

When oscillators are weakly coupled (coupling strength $\varepsilon \ll 1$), the averaging theorem shows that the amplitude dynamics decouple from the phase dynamics on the slow timescale. The phase interaction function $\Gamma(\theta)$ is obtained by averaging the coupling over one period:

$$\Gamma(\theta_j - \theta_i) = \frac{1}{T_i} \oint Z_i(\phi) \cdot h_{ij}(\phi, \phi + \theta_j - \theta_i) \, d\phi$$

where $Z_i$ is the phase response curve and $h_{ij}$ is the coupling function. For many physical systems, the leading Fourier mode of $\Gamma$ is sinusoidal:

$$\Gamma(\Delta\theta) \approx \sin(\Delta\theta)$$

This yields the Kuramoto model as the universal phase-reduction of weakly coupled oscillators.

2. The Kuramoto Equations

The all-to-all coupled Kuramoto model with $N$ oscillators:

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

The parameters are:

$\theta_i \in [0, 2\pi)$: phase of oscillator $i$.
$\omega_i$: natural frequency, drawn from a distribution $g(\omega)$.
$K \geq 0$: coupling strength. The only control parameter.

Each oscillator is attracted toward the mean phase of all others, but resists being pulled away from its natural frequency. The competition between disorder (spread of $\omega_i$) and coupling (strength $K$) drives the phase transition.

3. The Order Parameter

Kuramoto's brilliant insight was to introduce a complex-valued order parameter that measures the collective coherence:

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

Here $r \in [0, 1]$ is the synchronization magnitude and $\psi$ is the mean phase. The interpretation:

$r = 0$: Phases are uniformly distributed (complete incoherence). Each oscillator runs independently.
$r = 1$: All phases are identical (perfect synchronization). The population oscillates in unison.
$0 < r < 1$: Partial synchronization. A fraction of oscillators lock together while others drift.

Using the order parameter, the Kuramoto equations can be rewritten in mean-field form:

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

Each oscillator interacts not with every other oscillator individually, but with the collective field $(r, \psi)$. The coupling strength felt by each oscillator is $Kr$: stronger coupling and greater coherence both amplify the synchronizing force. This self-reinforcing feedback is the engine of the phase transition.

4. Self-Consistency and Critical Coupling

In the thermodynamic limit $N \to \infty$, we seek a stationary solution where $r$ is constant. Setting $\psi = 0$ by symmetry, the locked oscillators satisfy $\dot{\theta}_i = 0$, giving:

$$\sin\theta_i = \frac{\omega_i}{Kr}$$

This has a solution only when $|\omega_i| \leq Kr$. Oscillators with natural frequencies within this window lock to the mean field; those outside drift incoherently. The self-consistency equation for $r$ is obtained by computing the order parameter from the locked population:

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

For a Lorentzian (Cauchy) frequency distribution:

$$g(\omega) = \frac{\gamma}{\pi(\omega^2 + \gamma^2)}$$

the self-consistency equation can be solved exactly. The critical coupling is:

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

The phase transition has the characteristic square-root scaling near onset:

$$r = \begin{cases} 0 & K < K_c \\ \sqrt{\frac{16(K - K_c)}{K_c^3 \pi g''(0) + 16K_c}} \approx \sqrt{\frac{K - K_c}{K_c/2}} & K > K_c \end{cases}$$

For the Lorentzian, this simplifies to $r = \sqrt{1 - K_c/K}$. This is a supercritical pitchfork bifurcation in $r$ as $K$ crosses $K_c$.

5. Simulating the Phase Transition

We simulate the Kuramoto model for increasing coupling strength, measuring the steady-state order parameter to trace the phase transition. We also visualize the oscillators on the unit circle at different coupling strengths.

Kuramoto Model: Phase Transition and Unit Circle Visualization

Python

script.py109 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

np.random.seed(42)

N = 200         # number of oscillators
gamma = 1.0     # Lorentzian half-width

# Natural frequencies from Lorentzian distribution
omega = gamma * np.tan(np.pi * (np.random.rand(N) - 0.5))
omega = np.clip(omega, -10, 10)  # clip extreme outliers

def kuramoto_rhs(t, theta, K):
    """Right-hand side of Kuramoto equations."""
    dtheta = np.zeros(N)
    sin_diff = np.sin(theta[np.newaxis, :] - theta[:, np.newaxis])
    coupling = (K / N) * np.sum(sin_diff, axis=1)
    return omega + coupling

def order_parameter(theta):
    """Compute magnitude of order parameter."""
    z = np.mean(np.exp(1j * theta))
    return np.abs(z), np.angle(z)

# --- Sweep coupling strength ---
K_values = np.linspace(0, 6, 40)
r_steady = np.zeros(len(K_values))
T_transient = 50.0
T_measure = 20.0

print("Sweeping coupling strength K...")
for idx, K in enumerate(K_values):
    theta0 = 2 * np.pi * np.random.rand(N)
    # Integrate through transient
    sol = solve_ivp(lambda t, y: kuramoto_rhs(t, y, K),
                    [0, T_transient + T_measure], theta0,
                    t_eval=np.linspace(T_transient, T_transient + T_measure, 50),
                    method='RK45', max_step=0.5)
    # Average r over measurement window
    r_vals = [order_parameter(sol.y[:, k])[0] for k in range(sol.y.shape[1])]
    r_steady[idx] = np.mean(r_vals)

K_c = 2 * gamma  # theoretical critical coupling
r_theory = np.where(K_values > K_c, np.sqrt(1 - K_c / np.maximum(K_values, K_c + 1e-10)), 0)

print(f"Theoretical K_c = {K_c:.2f}")
print(f"Measured r at K=0: {r_steady[0]:.3f}")
print(f"Measured r at K=K_c: {r_steady[np.argmin(np.abs(K_values - K_c))]:.3f}")
print(f"Measured r at K=2*K_c: {r_steady[np.argmin(np.abs(K_values - 2*K_c))]:.3f}")

# --- Snapshots on the unit circle ---
K_snapshots = [0.5, K_c, 3.0, 5.0]
snapshot_thetas = []
for K in K_snapshots:
    theta0 = 2 * np.pi * np.random.rand(N)
    sol = solve_ivp(lambda t, y: kuramoto_rhs(t, y, K),
                    [0, T_transient + 10], theta0,
                    t_eval=[T_transient + 10],
                    method='RK45', max_step=0.5)
    snapshot_thetas.append(sol.y[:, -1])

# --- Plot ---
fig = plt.figure(figsize=(14, 10))

# Phase transition curve
ax1 = fig.add_subplot(2, 2, (1, 2))
ax1.plot(K_values, r_steady, 'o', color='#10b981', markersize=5, alpha=0.8, label='Simulation')
ax1.plot(K_values, r_theory, '-', color='#f59e0b', linewidth=2, label=r'Theory: $r = \sqrt{1 - K_c/K}$')
ax1.axvline(x=K_c, color='#ef4444', linestyle='--', linewidth=1.5, label=f'$K_c = {K_c:.1f}$')
ax1.set_xlabel('Coupling strength K', fontsize=13, color='white')
ax1.set_ylabel('Order parameter r', fontsize=13, color='white')
ax1.set_title('Kuramoto Phase Transition', fontsize=15, color='white', fontweight='bold')
ax1.legend(fontsize=11)
ax1.set_facecolor('#0a0a0a')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.2)
ax1.set_xlim(0, 6)
ax1.set_ylim(-0.05, 1.05)

# Unit circle snapshots
for idx, (K, theta) in enumerate(zip(K_snapshots, snapshot_thetas)):
    ax = fig.add_subplot(2, 4, 5 + idx, projection='polar')
    r_val, psi_val = order_parameter(theta)

# Plot individual oscillators
    ax.scatter(theta % (2*np.pi), np.ones(N), c='#10b981', s=8, alpha=0.5)

# Plot order parameter
    ax.annotate('', xy=(psi_val, r_val), xytext=(0, 0),
                arrowprops=dict(arrowstyle='->', color='#ef4444', lw=2.5))

ax.set_title(f'K={K:.1f}, r={r_val:.2f}', color='white', fontsize=11, pad=12)
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='gray', labelsize=7)
    ax.set_rticks([0.5, 1.0])
    ax.set_rlabel_position(45)
    ax.grid(True, alpha=0.2)

fig.patch.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print("\nPlot saved: Kuramoto phase transition with unit circle snapshots")
print(f"\nN = {N} oscillators, Lorentzian distribution with gamma = {gamma}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Stability of the Incoherent and Synchronized States

The incoherent state $r = 0$ is described by a uniform density on the circle. Its linear stability can be analysed using the Ott-Antonsen ansatz, which reduces the infinite-dimensional dynamics to a single ODE for the order parameter:

$$\frac{dz}{dt} = -\gamma z + \frac{K}{2}(z - z|z|^2)$$

where $z = re^{i\psi}$ is the complex order parameter. Linearising around $z = 0$:

$$\frac{dz}{dt} \approx \left(\frac{K}{2} - \gamma\right) z$$

The incoherent state loses stability when the growth rate becomes positive, confirming:

$$K_c = 2\gamma$$

For $K > K_c$, the cubic term stabilises the order parameter at the fixed point $|z|^2 = 1 - 2\gamma/K$, confirming the square-root scaling $r = \sqrt{1 - K_c/K}$.

7. Universality and Critical Exponents

The Kuramoto transition belongs to a universality class determined by the frequency distribution $g(\omega)$ near $\omega = 0$:

Smooth, unimodal $g$ with $g''(0) < 0$: Supercritical transition with $r \sim (K - K_c)^{1/2}$. Mean-field exponent $\beta = 1/2$.
Flat-top distribution $g''(0) = 0$: Higher-order transition with $r \sim (K - K_c)^{1/4}$.
Bimodal $g$: Can produce explosive (first-order) synchronization with hysteresis.

Key Takeaway

The Kuramoto model captures a genuine phase transition: below the critical coupling, no macroscopic order exists despite microscopic interactions; above it, a macroscopic fraction of oscillators spontaneously synchronizes. The order parameter $r$ serves as the magnetisation analogue, and the critical coupling $K_c = 2/(\pi g(0))$ plays the role of the critical temperature.

← Part 14 Overview Network Synchronization →