Continuous versus Discontinuous Transitions in the $D$-Dimensional Generalized Kuramoto Model: Odd $D$ is Different

The Kuramoto model, originally proposed to model the dynamics of many interacting oscillators, has been used and generalized for a wide range of applications involving the collective behavior of large heterogeneous groups of dynamical units whose states are characterized by a scalar angle variable. One such application in which we are interested is the alignment of orientation vectors among members of a swarm. Despite being commonly used for this purpose, the Kuramoto model can only describe swarms in two dimensions, and hence the results obtained do not apply to the often relevant situation of swarms in three dimensions. Partly based on this motivation, as well as on relevance to the classical, mean-field, zero-temperature Heisenberg model with quenched site disorder, in this paper we study the Kuramoto model generalized to $D$ dimensions. We show that in the important case of three dimensions, as well as for any odd number of dimensions, the $D$-dimensional generalized Kuramoto model for heterogeneous units has dynamics that are remarkably different from the dynamics in two dimensions. In particular, for odd $D$ the transition to coherence occurs discontinuously as the interunit coupling constant $K$ is increased through zero, as opposed to the $D=2$ case (and, as we show, also the case of even $D$) for which the transition to coherence occurs continuously as $K$ increases through a positive critical value $K_c$. We also demonstrate the qualitative applicability of our results to related models constructed specifically to capture swarming and flocking dynamics in three dimensions.

Popular Summary

Flashing fireflies, pulsing cardiac cells, and oscillating neutrinos are all examples of collective behavior and synchronization in nature. In 1975, physicist Yoshiki Kuramoto proposed a mathematical model to describe synchronization in a general way. This model describes the interaction between a large number of components whose time-dependent state is given solely by a single variable—for example, herding animals moving on a 2D surface. There are, however, many applications, such as the orientation of animals flying in the air, for which a single variable is insufficient. Motivated by such situations, we generalize the Kuramoto model to states characterized by an arbitrary number of dimensions.

We observe that the three-dimensional case, as well as other higher odd-dimensional cases, behaves remarkably differently from the dynamics of the well-understood case of two dimensions. In that original 2D case, a swarm achieves some small amount of coherence when the interagent coupling strength is larger than some critical amount, but we find that for odd dimensions, the swarm stays coherent for any coupling strength.

In addition, we provide new analytical methods for studying our extended Kuramoto model and build new understanding for a large class of problems relevant for a variety of contexts. We also lay the groundwork for future important extensions to this problem, such as the incorporation of network-based interaction and the generalization of the many important Kuramoto-like systems to higher dimensions.

Figure 1

Phase transitions for the standard two-dimensional Kuramoto model from theory (see Ref. [31]), shown as the black dashed curve, and for the Kuramoto model generalized to three dimensions as calculated from the theory in Eq. (18), shown as the solid red curve. Note the continuous transition in the two-dimensional Kuramoto model at a critical coupling of $K_c>0$, and the discontinuous transition of the three-dimensional Kuramoto model at $K_c=0$. The blue dotted line represents the maximum possible value of coherence, corresponding to $|\mathit{\rho }|=1$.

Figure 2

Illustration showing that $[{\mathit{\sigma }}_j-({\mathit{\sigma }}_j·{\mathit{\sigma }}_i){\mathit{\sigma }}_i]={\widehat{\mathit{\theta }}}_i\sin ({\theta }_j-{\theta }_i)$, where ${\widehat{\mathit{\theta }}}_i$ is a unit vector in the direction of increasing ${\theta }_i$.

Figure 3

${\mathit{\sigma }}_i$ precesses around ${\widehat{\mathit{\omega }}}_i$ with an angular frequency of ${\omega }_i$ for $K=0$ with $D=3$.

Figure 4

(a) Phase transitions for the generalized Kuramoto model for $D=2$ (red plus signs), 4 (blue inverted triangles), 6 (green triangles), and 8 (magenta stars) dimensions, numerical observations for $N={10}^5$. (b) Phase transitions for the Kuramoto model generalized to $D=3$ (red plus signs), 5 (blue inverted triangles), 7 (green triangles), and 9 (magenta stars) dimensions, numerical observations for $N={10}^4$. $\mathrm{\Delta }=1$ in each. The theoretical predictions from Eqs. (42) and (44) for the critical coupling strength for even dimensions has been shown in correspondingly colored arrows on the $x$ axis in (a). For a discussion on the slight mismatch between the theory and the numerical results, see Sec. 3c. We expect this mismatch to decrease with increasing $N$. The theoretical estimates from Eq. (13) for the magnitude of the discontinuity, i.e., $|\mathit{\rho }|$ at $K\to 0^{+}$, are shown in correspondingly colored arrows on the $y$ axis in (b). Note the close match between the theoretical result for the discontinuity and the numerical observation for $|\mathit{\rho }|$ at $K=0.2$.

Figure 5

Orientations of each of the $N=5\times {10}^3$ agents at time $T=0$, $T=500$, and $T=50000$. We visualize the orientations by plotting the end points of the orientation vectors on the sphere $\mathcal{S}$. The sphere is then mapped onto a rectangle using an area-preserving transformation. We choose the $z$ axis along $\mathit{\rho }$, and arbitrarily choose mutually orthogonal $x$ and $y$ axes. $\theta$ then represents the angle measured from the $z$ axis ($\cos \theta =\widehat{\mathit{\rho }}·\mathit{\sigma }$), and $\varphi$ represents the azimuthal angle measured anticlockwise from the $x$ axis. The agent state vectors are initialized with a uniform distribution on the sphere and evolved with a coupling constant $K=0.1$. Note how all of the agents tend to uniformly distribute themselves on one hemisphere.

Figure 6

Phase transition for Kuramoto model generalized to three dimensions. The solid (dashed) black curve represents the derived stable (unstable) fixed points for the order parameter, and the red plus sign markers represent numerical results from a simulation with $N={10}^4$ agents and $\mathrm{\Delta }=1$.

Figure 7

$N=1000$ agents were simulated with a coupling strength of $K=-2$. Histograms of $z_i={\mathit{\sigma }}_i·{\widehat{\mathit{\omega }}}_i$ have been plotted at $T=0$ (in blue, corresponding to the initial condition having ${\mathit{\sigma }}_i$ uniformly spread on $\mathcal{S}$) and after $T=1.25\times {10}^5$ time units (in red). Note how the distributions concentrate at 1 and $-1$ for large $T$. In the insets, we show plots of $z_i$ as a function of time for 50 randomly chosen agents.

Figure 8

Illustration of an extended-body agent. Unlike the agents in the generalized Kuramoto model Eq. (2), we assume that the state of an extended-body agent cannot be described by a single unit vector $\mathit{\sigma }$. Rather, the pair of vectors $\mathit{\sigma }$ and $\mathit{\eta }$ together describe the orientation and state of the agent. The direction of agent velocity is assumed to be along the direction $\mathit{\sigma }$ as earlier. The unit vector $\mathit{\nu }$ is defined as $\mathit{\nu }=\mathit{\sigma }\times \mathit{\eta }$.

Figure 9

Phase transition for interacting three-dimensional extended-body agents described by Eqs. (56)–(58). The dynamics of individual agents in this system have been constructed to satisfy constraints imposed by extended-body dynamics, and are not equivalent to the dynamics of the generalized Kuramoto model described in Sec. 3. Despite this, we continue to observe a discontinuous jump in the asymptotic steady-state value of $|\mathit{\rho }|$ as $K$ is increased through 0.

Figure 10

The shaded regions (blue, green, and orange) correspond to the domain $\mathrm{\Gamma }$ in which ${\sum }_{k}1/{\mu }_k^2>1$ for the case of $D=4$ ($\mathrm{\Lambda }=2$) in the $\{{\mu }_1,{\mu }_2\}$ space. The subdomain ${\mathrm{\Gamma }}_0$, shown in blue, is the part of $\mathrm{\Gamma }$ inside the circle of radius $L$, and the subdomains ${\mathrm{\Gamma }}_i$ are the parts of the domain $\mathrm{\Gamma }$ that lie outside ${\mathrm{\Gamma }}_0$ which do not contain the ${\mu }_i$ axis (${\mathrm{\Gamma }}_1$ is shown in orange and ${\mathrm{\Gamma }}_2$ in green). The width of the strips in $\mathrm{\Gamma }$ far away from the origin is 1; hence, the volume of the subdomain ${\mathrm{\Gamma }}_0$ will scale as $\mathcal{O}(L^{\mathrm{\Lambda }-1})$ for large $L$.