Introduction

Synchronization — the spontaneous emergence of coordinated timing among coupled oscillators — appears across vastly different physical and biological systems. Fireflies in South-East Asia flash in near-perfect unison[buck1988synchronous]. Power-grid generators must lock to the same frequency or risk cascading failure[dorfler2013synchronization]. Cardiac pacemaker cells synchronize to produce a reliable heartbeat[winfree1967biological]. And in classical ballet, a corps de ballet of twelve or more dancers must synchronize both timing and spatial position to achieve the aesthetic ideal of a unified visual line.

The Kuramoto model[kuramoto1975self] is the canonical mathematical framework for studying phase synchronization in populations of coupled oscillators. Each oscillator $i$ has a natural frequency $\omega_i$ and a phase $\theta_i$ that evolves according to $$\dot{\theta}_i = \omega_i + \frac{K}{|{\mathcal{N}$$i}|} \sum{j \in \mathcal{N}_i} \sin(\theta_j - \theta_i), where $K \geq 0$ is the coupling strength and $\mathcal{N}$iNi​ is the set of neighbors of oscillator $i$. When $K$ exceeds a critical value $K_c$, the system undergoes a continuous phase transition from incoherence to partial (and eventually full) synchronization. The global synchrony is quantified by the Kuramoto order parameter $$r(t) = \frac{1}{N} \left| \sum$${j=1}^{N} e^{i\theta_j(t)} \right|, where $r = 0$ indicates fully desynchronized phases and $r = 1$ indicates perfect unison.

The vast majority of Kuramoto studies assume homogeneous, all-to-all coupling. Two extensions are, however, directly relevant to real ensemble performance. First, spatial embedding: dancers occupy fixed positions on a stage, and visual or auditory coupling naturally diminishes with distance. Second, \emph{hierarchical structure}: professional ballet ensembles are organized around a principal dancer whose movement is followed by soloists, who in turn anchor the corps. Whether this hierarchy lowers the coupling threshold for synchronization — and if so, by how much — is an open empirical question with direct implications for ensemble pedagogy.

We make three contributions:

• A spatially-embedded Kuramoto simulator with four configurable network topologies (all-to-all, nearest-$k$, hierarchical, ring) and four domain presets (ballet-corps, fireflies, drum-circle, power-grid), integrated with 4th-order Runge-Kutta (RK4) for numerical accuracy near the phase transition.
• A multi-evaluator panel combining four complementary synchronization metrics with three aggregation strategies, enabling robust $K_c$ detection across topology types.
• An agent-executable skill (SKILL.md) encoding the full 1,440-simulation experiment pipeline, reproducible by any AI coding agent without internet access or GPU.

Model

Spatially-Embedded Kuramoto Dynamics

We simulate $N$ dancer-agents on a 2D stage of side $L = 10.0$. Each agent $i$ is characterized by its phase $\theta_i \in [0, 2\pi)$, natural frequency $\omega_i ~ \mathcal{N}(\omega_0, \sigma^2)$ (drawn once at initialization and held fixed), and a fixed spatial position $(x_i, y_i)$ on the stage. Agents do not move; the "dance" is the phase oscillation representing one full movement cycle (e.g., a repeated ballet step). Synchronization means all agents arrive at the same phase simultaneously.

Dynamics follow Eq. \eqref{eq:kuramoto} and are numerically integrated via RK4: $$k_1 = \Delta t \cdot f(\theta(t)),$$

$$k_2 = \Delta t \cdot f(\theta(t) + k_1/2),$$

$$k_3 = \Delta t \cdot f(\theta(t) + k_2/2),$$

$$k_4 = \Delta t \cdot f(\theta(t) + k_3),$$

$$\theta(t + \Delta t) = \theta(t) + \frac{k_1 + 2k_2 + 2k_3 + k_4}{6},$$ with $\Delta t = 0.01$ and $T = 1,000$ timesteps ($t_{\max} = 10$ time units). RK4 is preferred over Euler integration because first-order errors can shift the apparent $K_c$ when the order-parameter gradient is steep near the transition.

Network Topologies

Four topologies define the neighbor set $\mathcal{N}_i$ in Eq. \eqref{eq:kuramoto}:

• All-to-all. Every dancer sees every other ($|\mathcal{N}_i| = N - 1$). Corresponds to open formations with full mutual visibility.
• Nearest-$k$. Each dancer couples to the $k = 4$ spatially nearest neighbors. Models tight formations where local visual cues dominate.
• Hierarchical. One principal dancer couples to all soloists; each soloist couples to the principal and its assigned corps subset; each corps member couples only to its soloist. This three-tier tree reflects actual ballet company structure.
• Ring. Each dancer couples to two circular neighbors. Models circular formations common in folk and processional dance.

Domain Presets

Four presets fix $(N, \sigma, \text{topology})$ to represent real-world synchronization contexts beyond ballet:

\caption{Domain presets. $\sigma$ controls frequency heterogeneity; larger $\sigma$ requires higher $K$ to achieve synchrony.}

Analytical Baseline

For all-to-all coupling with Gaussian natural frequencies, the critical coupling is known analytically[strogatz2000kuramoto]: $$K_c^{\text{all-to-all}} = \frac{2}{\pi g(\omega_0)} = \frac{2\sigma\sqrt{2\pi}}{\pi} \approx 1.596 \sigma,$$ where $g(\omega)$ is the Gaussian density evaluated at the mean frequency. For $\sigma = 0.5$, this gives $K_c \approx 0.798$. No closed-form expression exists for the other topologies; their empirical $K_c$ values constitute the primary scientific output.

Above $K_c$, the order parameter grows as $$r \propto (K - K_c)^\beta,$$ with $\beta = 0.5$ for the mean-field (all-to-all) universality class[acebron2005kuramoto]. Deviations from $\beta = 0.5$ for non-mean-field topologies would indicate a different universality class.

Evaluator Panel

Synchronization Metrics

Four evaluators independently score the phase history on $[0, 1]$ (0 = incoherent, 1 = perfect sync):

Evaluator 1: Kuramoto Order Parameter (KOP). Computes $\bar{r} = \langle r(t) \rangle$ over the final 20% of timesteps. Theoretically grounded and directly comparable to Eq. \eqref{eq:kc_analytic}. Weakness: ignores spatial arrangement of dancers.

Evaluator 2: Spatial Alignment (SA). Computes the Pearson correlation $\rho$ between pairwise spatial distances and pairwise circular phase distances (using $\min(|\theta_i - \theta_j|, 2\pi - |\theta_i - \theta_j|)$) over the final 20% of timesteps. Sync score $= \max(0, 1 - \rho)$, clipped to $[0,1]$. Positive $\rho$ means nearby dancers are out of phase (undesirable). Captures spatial coherence relevant to ballet aesthetics; less informative for all-to-all topology.

Evaluator 3: Velocity Synchrony (VS). Computes the variance of instantaneous angular velocities $\dot{\theta}_i$ across agents in the final 20% of timesteps, normalized by $\sigma^2$ (the variance at $K = 0$): sync score $= \max(0, 1 - \text{Var}(\dot{\theta}) / \sigma^2)$. Detects frequency-locking even without phase alignment; a frequency-synchronized but phase-offset ensemble still scores high.

Evaluator 4: Pairwise Entrainment (PE). For each connected pair $(i, j)$, measures the variance of the phase difference $|\theta_i(t) - \theta_j(t)|$ over the final 20% of timesteps. A pair is entrained if this variance is below 0.1. Sync score $=$ fraction of connected pairs that are entrained. Provides fine-grained connection-level evidence; most computationally demanding.

Panel Aggregation

Three aggregation strategies are reported for every experimental condition:

• Majority vote. Synchronized if $\geq 3$ of 4 evaluators score $> 0.5$.
• Weighted average. Scores weighted by empirical reliability, estimated on calibration runs at $K = 0$ (ground-truth incoherence) and large $K$ (ground-truth synchrony).
• Unanimous. All 4 evaluators must exceed 0.5 (most conservative; minimizes false positives).

Results

Phase Transition Curves

Figure (generated by run.py) shows the order parameter $r(K)$ for each topology at $N = 12$, $\sigma = 0.5$, averaged over 3 seeds. All topologies display an S-shaped transition from incoherence ($r \approx 0$) to synchrony ($r \approx 1$). Table reports the empirical $K_c$ estimated via two independent methods: sigmoid fit to $r(K)$ and susceptibility peak.

\caption{Critical coupling strength $K_c$ by topology ($N=12$, $\sigma=0.5$, mean $\pm$ std over 3 seeds; 95% bootstrap CI). Analytical all-to-all value: $K_c \approx 0.798$.}

Susceptibility and Critical Exponent

The susceptibility $\chi(K) = N \cdot \text{Var}_{\text{seeds}}[r(K)]$ peaks sharply at $K_c$ and provides a higher-resolution estimate than the sigmoid fit alone. Table reports agreement between the two estimation methods; discrepancies larger than 0.05 indicate finite-size effects.

Above $K_c$, we fit the scaling relation Eq. \eqref{eq:critical_exponent} via log-log regression of $r$ versus $(K - K_c)$ for $K \in [K_c, K_c + 0.5]$. The estimated critical exponents are:

\caption{Critical exponent $\beta$ by topology. Mean-field theory predicts $\beta = 0.5$ for all-to-all coupling. Values are from log-log regression ($R^2$ shown).}

A deviation of $\beta$ from 0.5 for the hierarchical or ring topologies would indicate a different universality class than classical mean-field Kuramoto — a physically significant finding.

Finite-Size Scaling

We exploit the three group sizes ($N \in {6, 12, 24}$) to perform finite-size scaling. The critical coupling shifts with system size according to $$K_c(N) = K_c(\infty) + a \cdot N^{-\nu},$$ where $K_c(\infty)$ is the thermodynamic-limit critical coupling and $\nu$ is the finite-size scaling exponent. We fit Eq. \eqref{eq:fss} by nonlinear least squares to the three $(N, K_c(N))$ data points per topology. With only three size points the fit is under-determined; we obtain $\nu = 1.0$ for all topologies, suggesting that larger system sizes are needed to resolve finite-size scaling exponents reliably.

Heterogeneity Effect

Comparing the two frequency-spread levels ($\sigma = 0.3$ vs. $\sigma = 0.8$) across all topologies and group sizes confirms the theoretical prediction (Eq. \eqref{eq:kc_analytic}): higher heterogeneity requires stronger coupling for synchronization. The empirical ratio $K_c(\sigma{=}0.8) / K_c(\sigma{=}0.3)$ is approximately $2.7$ for all-to-all (matching the prediction); the analytically predicted ratio for all-to-all coupling is $0.8 / 0.3 \approx 2.67$.

Evaluator Agreement

Table reports pairwise evaluator agreement (fraction of 1,440 simulations where both evaluators classify sync/no-sync identically at the majority-vote threshold) and the F1 score of each panel aggregation method against the KOP ground truth ($r > 0.5$).

Pairwise evaluator agreement and panel method F1 scores.

The Velocity Synchrony evaluator is expected to detect the transition earlier (at lower $K$) than Pairwise Entrainment because frequency-locking precedes phase-locking. Evaluator agreement near $K_c$ is expected to be lowest, since that is precisely where the system is most ambiguously poised between the two phases.

Discussion

Dance pedagogy implications. If the hierarchical topology achieves synchrony at a lower $K_c$ than flat topologies, it provides a quantitative justification for the traditional ballet corps structure: the principal-soloist-corps hierarchy is not merely an aesthetic convention but an efficient synchronization architecture. Ensemble directors could use the empirical $K_c$ estimates to calibrate rehearsal intensity — lower coupling overhead means faster synchrony with less practice time. The finite-size scaling result quantifies how smaller ensembles (e.g., pas de quatre) differ from full corps synchronization dynamics.

Generalization beyond dance. The four domain presets demonstrate the model's applicability across domains: fireflies ($N=50$, $\sigma=1.0$, nearest-$k$) tests whether biological flash synchronization falls in the same universality class as ballet; power-grid ($N=10$, $\sigma=0.2$, ring) probes whether the low-heterogeneity regime produces sharper transitions, relevant to frequency-control engineering. The multi-evaluator panel generalizes readily: any of the four metrics can serve as the primary detection criterion in domains where one measure is more interpretable (e.g., pairwise entrainment is directly meaningful for power-grid stability monitoring).

Limitations. This framework treats dancer positions as fixed throughout the simulation, whereas real choreography involves continuous spatial repositioning. The Kuramoto model uses simple sinusoidal coupling; actual inter-dancer cueing involves visual, auditory, and haptic channels with non-trivial delay and gain structures. The natural frequency $\omega_i$ is drawn once and held fixed, whereas real dancers continuously adapt their tempo. Finally, our "hierarchical" topology assumes a clean tree structure; real ballet hierarchies may have cross-cutting couplings (e.g., two soloists watching each other) not captured here.

Future work. Natural extensions include: (i) mobile agents with choreographed spatial trajectories, (ii) adaptive frequencies capturing deliberate tempo modulation, (iii) time-delayed coupling modeling the finite reaction time of human dancers, (iv) application to 3D drone swarms where spatial embedding in three dimensions introduces new topological effects, and (v) fitting $\sigma$ and $K$ parameters to motion-capture data from real ballet performances to test whether the Kuramoto model quantitatively predicts observed synchronization statistics.

\section*{Reproducibility}

All code and results are packaged as an executable SKILL.md. Any AI coding agent can reproduce all 1,440 simulations by following the skill steps: set up the Python virtual environment, run unit tests with pytest, execute run.py, and validate outputs with validate.py. No API keys, GPU, or external data downloads are required; the framework is a pure Python simulation relying on numpy, scipy, and matplotlib with pinned versions.

References

• [kuramoto1975self] Y. Kuramoto, "Self-entrainment of a population of coupled non-linear oscillators," in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, vol. 39, pp. 420--422. Springer, Berlin, 1975.

• [strogatz2000kuramoto] S. H. Strogatz, "From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators," Physica D: Nonlinear Phenomena, vol. 143, no. 1--4, pp. 1--20, 2000.

• [acebron2005kuramoto] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, "The Kuramoto model: a simple paradigm for synchronization phenomena," Reviews of Modern Physics, vol. 77, no. 1, pp. 137--185, 2005.

• [buck1988synchronous] J. Buck, "Synchronous rhythmic flashing of fireflies. II," The Quarterly Review of Biology, vol. 63, no. 3, pp. 265--289, 1988.

• [dorfler2013synchronization] F. D"{o}rfler and F. Bullo, "Synchronization in complex networks of phase oscillators: a survey," Automatica, vol. 50, no. 6, pp. 1539--1564, 2014.

• [winfree1967biological] A. T. Winfree, "Biological rhythms and the behavior of populations of coupled oscillators," Journal of Theoretical Biology, vol. 16, no. 1, pp. 15--42, 1967.

• [rodrigues2016kuramoto] F. A. Rodrigues, T. K. DM. Peron, P. Ji, and J. Kurths, "The Kuramoto model in complex networks," Physics Reports, vol. 610, pp. 1--98, 2016.

• [breakspear2010generative] M. Breakspear, S. Heitmann, and A. Daffertshofer, "Generative models of cortical oscillations: neurobiological implications of the Kuramoto model," Frontiers in Human Neuroscience, vol. 4, p. 190, 2010.