Exploring the phase diagrams of multidimensional Kuramoto models

In two dimensions we can use the dimensional reduction approach of Ott and Antonsen for Lorentzian distributions of natural frequencies [45]. In this case one can derive differential equations for the module and phase of the order parameter [21] and we will use these equations to construct the full 2D phase diagram analytically. Taking the limit $N\rightarrow\infty$ we define the function $f(\omega,\theta,t)$ describing the density of oscillators with natural frequency $\omega$ at position $\theta$ in time $t$. Since the total number of oscillators is conserved, $f$ satisfies a continuity equation with velocity field

where $\vec{\sigma}=(\cos\theta,\sin\theta)$, $\vec{q}\equiv q(\cos\xi,\sin\xi)$ and $\hat{\theta}=(-\sin\theta,\cos\theta)$. The continuity equation reads

The Ott-Antonsen ansatz consists in expanding $f$ in Fourier series and choose the coefficients so that all of them depend on a single complex parameter $\nu(t)$ as

Inserting (18) and (16) into (17) we obtain the following differential equation for $\nu(t)$:

where we defined the complex number $u=qe^{i\xi}$. Setting $z=pe^{i\psi}$ and using the definition $\vec{q}=\mathbf{K}\vec{p}$ (see Eq.(7)), we find

The last step is to relate the ansatz parameter $\nu$ with the complex order parameter $z$. To do that we note that in the limit of infinitely many oscillators, Eq.(2) becomes

Using the Fourier series (18) we see that only the term proportional to $e^{-i\theta}$ contributes to the integral and we are left with

This equation can be integrated exactly for

[45]. In this case we can write $\nu=\rho e^{i\Phi}$ and we obtain

The integral can now be performed in the complex $\omega$-plane using a closed path from $-R$ to $+R$ over the real line and closing back from $R$ to $-R$ with a half circle in the positive complex half plane, taking $R\rightarrow\infty$. From Eq.(19) we see that $\Phi$ should be proportional to $\omega$ for large $\omega$, implying that the integral over the half circle goes to zero. Using the residues theorem at the pole $\omega=\omega_{0}+i\Delta$ we obtain

Calculating Eq.(19) at $\omega_{0}+i\Delta$ allows us to replace $\nu$ by $z$, resulting in

Finally, separating real and imaginary parts we obtain equations for the module and phase of the order parameter $z=pe^{i\psi}$:

and

where we have defined $\theta=2\psi+\beta$.

Non-trivial stationary solutions of Eqs.(27) and (28) are given by

and

where $a=2\omega_{0}-(1+p^{2})K\sin\alpha$ and $b=J(1+p^{2})$. These solutions are real provided

and $|a|<|b|$, or

Therefore, to find $p$ and $\psi$ for each pair $(\omega_{0},\Delta)$ we need to solve Eqs.(29) and (30) and check that conditions (31) and (32) are satisfied.

The trivial (asynchronous) state $p=0$ is always a solution of Eq.(27). Although the phase of $z$ for $p=0$ is mostly irrelevant, it plays a role in analysis of its stability. At $p=0$ the linearized version of Eq.(27) is

whose solution is

Therefore, if $\theta$ oscillates (for $\omega_{0}$ outside the interval in Eq.(32)) $p=0$ is stable for $\Delta>K\cos\alpha/2\equiv\Delta_{0}$. When $\theta$ converges to a constant, stability requires $\Delta>K\cos\alpha/2-J\cos\theta/2$. Since the linearized equation for $\theta$ at $p=0$ is $\delta\dot{\theta}=J\cos\theta\delta\theta$, for $J>0$, the trivial solution is stable if $\cos\theta<0$ ($\pi 2/<\theta<3\pi/2$) and for $J<0$ if $\cos\theta>0$ ($-\pi/2<\theta<\pi/2$).

For $\Delta>\Delta_{0}$ the line separating the async from the static sync region is given in parametric form by $(\omega(\theta),\Delta(\theta))$ with $\omega(\theta)=(K\sin\alpha-|J|\sin\theta)/2$, $\Delta(\theta)=(K\cos\alpha-|J|\cos\theta)/2$, for $\theta\in(\pi/2,3\pi/2)$.

For $\Delta<\Delta_{0}$ the solution $p=0$ is unstable and $p$ either converges to a constant (static solutions) or it oscillates (active solutions). The boundary between the two kinds of behavior is obtained setting $\sin\theta=\pm 1$ and $\cos\theta=0$. From Eq.(29) we find

Setting $a=b$ we obtain

where $\omega_{max}\equiv(K\sin\alpha+|J|)/2$. Equating these two relations gives the boundary curve

Setting $a=-b$ gives the other boundary curve

where $\omega_{min}\equiv(K\sin\alpha-|J|)/2$. The value of $p$ along the curve is given br Eq.(29) with $\cos\theta=0$, and it approaches 1 as $\Delta\rightarrow 0$. The phase is $\psi=3\pi/4-\beta$ on the upper curve and $\psi=\pi/4-\beta$ on the lower curve. Between these two curves the solution is static and outside this interval the phase $\theta$, and therefore $\psi$, oscillates, leading to an oscillation of $p(t)$, corresponding to active states.

For $J=0$ the static sync region collapses and the active regions become regions of rotation, where the module of $p$ remains constant whereas the phase $\psi$ rotates with angular velocity $\omega_{0}-K\sin\alpha+\Delta\tan\alpha$. A line of static sync appears for $\omega_{0}=K\sin\alpha-\Delta\tan\alpha$. The full diagram is illustrated in Figure 1.

Although in $D=2$ we have the exact phase diagrams for the Lorentz distributions, simulations are important for two reasons: first, we don’t have the diagrams for other distributions; second, since we need to automatize the simulations for other distributions and higher dimensions, we need to make sure we can extract the different phases from the simulated diagrams. Therefore, we first simulate the diagrams for the Lorentz distribution

and see how they compare with the exact solutions. We also simulated the equations for the Gaussian

and uniform distributions,

Note that $\Delta$ is a measure of the width of the distributions, but has different meanings in each case. For the Gaussian distribution $\Delta^{2}$ is the variance; for the uniform distribution the variance is $\Delta^{2}/12$ whereas the Lorentz has infinite variance.

For each distribution we simulated the dynamics with three coupling matrices (see Eq.(9)):

and

These choices correspond to coupling matrices with real eigenvalues, leading to static states for $\Delta=\omega_{0}=0$ ($\mathbf{K}_{S}$), purely complex eigenvalues, leading to rotations as in the Kuramoto-Sakaguchi model ($\mathbf{K}_{R}$) and complex eigenvalues corresponding to active states ($\mathbf{K}_{A}$).

Fig. 2 displays results for the Lorentz distribution, which we simulate only for $\mathbf{K}_{S}$ (top panels) and $\mathbf{K}_{R}$ (lower panels), as our intention is to validate the numerical results comparing with the exact diagrams. Continuous black lines show the boundary curves as in Fig.1 and they divide the plane int o three (top) or two (bottom) regions. In both cases the right most region corresponds to non-synchronized states. The top panels show static synchronization, with $\vec{p}$ constant, in the lower left corner, as can be seen by the low values of both $\delta p$ and $\delta_{max}$. The upper left region, on the other hand, shows active states, with $\langle p\rangle$ constant but rotating and oscillating $\vec{p}(t)$, as indicated by significant values of $\delta p$ and large values of $\delta_{max}$. For the bottom panels, corresponding to the Kuramoto-Sakaguchi model, $\vec{p}$ rotates but keeps its module constant ($\delta p\approx 0$).

Fig. 3 shows similar plots for the Gaussian distribution, Eq.(40). The top two panels are qualitatively similar to the panels in Fig.2, although not identical: critical transition values of $\Delta$ and $\omega_{0}$ are different and transitions are much sharper, as natural frequencies far from $\omega_{0}$ are much less likely to be sampled. Interestingly, for the case of $\mathbf{K}_{A}$, the line of fixed $\vec{p}$ immersed in the rotating zone (see Fig.1(b) ) is enlarged into an area (red stripe in the $\delta_{mas})$, showing the great sensitivity of phase diagram with the coupling matrix.

Finally, Fig. 4 shows results for the uniform distribution, Eq.(41) and the same coupling matrices, Eqs.(42) -(44). Except for the case of $\mathbf{K}_{R}$ (middle row) which is qualitatively similar to the cases of Lorentz and Gaussian distributions, new regions develop in this case. For $\mathbf{K}_{S}$ part of the region corresponding to non-synchronized states for Gaussian and Lorentz distributions becomes partially synchronized with static states (first row) and for $\mathbf{K}_{A}$ the active states near the line of static states develop larger oscillations (yellow area in the middle plot in the third row). Enlargement of the line of static sync states is also observed in this case.

To explore the phase diagrams in 3D we set the coupling matrix as in Eq.(11) with $a=1$ and $\alpha=0.5$ and consider two values of the parameter $b$. For $b=0.5$ we define

and for $b=1$

In the first case the dominant eigenvalues have complex eigenvectors in the $\hat{e}_{1}$-$\hat{e}_{2}$ plane, corresponding to rotations for $\Delta=\omega_{0}=0$. In the second case the dominant eigenvector is real in the $\hat{e}_{3}$ direction leading to static synchronization. For each case we use three different types of natural frequencies distributions, described either in terms of the matrix $\mathbf{W}_{i}$ entries as in Eq.(11) or in terms of the associated vector $\vec{\omega}_{i}$, Eq.(13), as follows:

Figs. 5 to 7 show results of numerical simulations in each case. As in the 2D case we show the time-average value of the order parameter after the transient, $\langle p\rangle$, its mean square deviation $\delta p$, and $\delta_{max}$, the maximum between mean square deviation of the components of $\vec{p}$, indicating if $\vec{p}$ is rotating ($\delta_{max}>0$) or not ($\delta_{max}=0$). Capital letters indicate the asymptotic state of the system as in $D=2$. In all cases, when $\Delta$ and $\omega_{0}$ are sufficiently small the oscillators synchronize and rotate (case 1, first line in all figures) or converge to a static configuration (case 2, second line of figures). However, as $\Delta$ and $\omega_{0}$ increase, the behavior of the system depends significantly on the distribution of natural frequencies.

For $g_{ang}(\vec{\omega})$, synchronization decreases as $\Delta$ and $\omega_{0}$ increase. For the coupling matrix $\mathbf{K}_{3R}$ a phase transition from rotation (R) to static sync (S) and back to rotation (R) is observed as $\Delta$ and $\omega_{0}$ increase. A thin region of active states is also noted for small $\Delta$ and large $\omega_{0}$. For $\mathbf{K}_{3S}$ the diagram is dominated by a large area of static sync (S), although a similar region of active states is observed, next to rotations (R). Notice that, since the direction of the vectors $\vec{\omega}_{i}$ is uniformly sampled, the average value of these vectors is zero for $g_{ang}$, independent of the value of $\omega_{0}$. This makes the phases observed at $\omega_{0}=\Delta=0$ to extend over large regions of the diagram as compared to the other two distributions in Figs. 6 and 7.

The phase diagrams for $g_{gauss}$ and $g_{uni}$ are similar, but very different from that of $g_{ang}$. Synchronization is much facilitated in these cases, as we don’t see regions of disordered motion in this range of parameters. Moreover, states with nearly complete sync are possible even for large $\Delta$ in the static phase. The phase of active states (A) is also much larger than in Fig. 5 and occurs for small values of $\omega_{0}$ and large values of $\Delta$. Pure rotations tend to be suppressed for $g_{gauss}$, occurring in small regions for both $\mathbf{K}_{3R}$ and $\mathbf{K}_{3S}$ and in a larger region $g_{uni}$, especially for $\mathbf{K}_{3S}$. Finally we note that the uniform distribution, Fig. 7, produces sharper transitions between the different phases.

In this section we illustrate the phase diagrams in $D=4$ with two instances of the coupling matrix Eq. (14). In both cases we fixed $\alpha=0$, $\beta=0.5$ and $a_{2}=0.8$. Similar to the 3D system we choose the remaining parameters $a_{1}$ and $b$ as follows: first we set $a_{1}=2.5$ and $b=0.5$, so that the leading eigenvalue of $\mathbf{K}$ is real in the direction of $\hat{e}_{1}$:

For the second case we set $a_{1}=0.5$ and $b=2.5$, with a complex leading eigenvalue in the $\hat{e}_{3}$-$\hat{e}_{4}$ plane:

In $D=4$ there are six independent entries for each matrix of natural frequencies $\mathbf{W}_{i}$ and many ways to choose these values. For each choice of coupling matrix above we performed simulations for the following three distributions:

Fig. 8 shows results for $g_{ang4}$. Similar to the 3D case with $g_{ang}$, the average of the natural frequencies is zero for all $\omega_{0}$ and $\Delta$, since the entries of $\mathbf{W}_{i}$ are uniformly distributed in all directions, with only its average intensity controlled by $\omega_{0}$. The basic phases S and R dominate much of the phase diagrams for $\mathbf{K}_{4S}$ and $\mathbf{K}_{4R}$ respectively.

For $g_{gauss4}$ and $g_{uni4}$ we obtain qualitatively similar diagrams (Figs. 9 and 10), where rotations are suppressed for $\mathbf{K}_{4S}$ and active states are suppressed for $\mathbf{K}_{4R}$. The phase diagrams are very different from their 3D counterparts exhibited in Figs. 6 and 7, but somewhat similar to the 2D case, Figs. 3 and 4, highlighting the role of dimensional parity (even or odd) in the dynamics and equilibrium properties of the model [37]. Synchronization happens only for limited values of $\Delta$, especially for the Gaussian case.