Geometry unites synchrony, chimeras, and waves in nonlinear oscillator networks

One of the simplest mathematical models in the study of nonlinear systems is the Kuramoto model, which describes synchronization in systems from swarms of insects to superconductors. We have recently found a connection between the original, real-valued nonlinear Kuramoto model and a corresponding complex-valued system that permits describing the system in terms of a linear operator and iterative update rule. We now use this description to investigate three major synchronization phenomena in Kuramoto networks (phase synchronization, chimera states, and traveling waves), not only in terms of steady state solutions but also in terms of transient dynamics and individual simulations. These results provide new mathematical insight into how sophisticated behaviors arise from connection patterns in nonlinear networked systems.


ARTICLE
scitation.org/journal/cha stable, steady state solutions but also in terms of transient behaviors in finite time. The original Kuramoto model (KM) with phase-lag is given bẏ where θ i (t) ∈ R is the phase of the ith oscillator at time t, ω i is its natural frequency, N is the number of oscillators, scales the coupling strength, A is the adjacency matrix, and φ is the phase lag. We focus here on the case where all oscillators have the same natural frequency (ω i = ω for all i ∈ [1, N]). The value of φ transforms the interaction function from the standard attractive case (φ = 0) to the cosine version with neutral coupling (φ = π/2). 21 Following recent work, 22 we can use a complex-valued, algebraic approach to the Kuramoto model. Starting with Eq. (1), we can change to a rotating coordinate frame 23 and set ω = 0 without loss of generality. By subtracting an additional imaginary component in the interaction terṁ (2) we arrive at a complex-valued system in a new variable ψ i ∈ C that provides an analytical, algebraic approach to the dynamics of the original nonlinear KM. Specifically, while the addition of the imaginary component results in a modified system, using the analytical solution of Eq. (2) (detailed below) to propagate across multiple time windows produces trajectories that match those in the original KM for long timescales (Figs. S1, S2, and S3 in the supplementary material). It is important to note that this specific form of the additional imaginary component is necessary (cf. supplementary material, Sec. I B).
We now show that Eq.
(2) has a closed-form, analytical solution. To do this, we can now multiply by i and use Euler's formula to obtain (cf. Sec. I B in the supplementary material) Letting x i = e iψ i , we can writė which has the general solution where x = (x 1 , . . . , x N ) and where K now collates the network topology A, the coupling strength , and the phase lag φ. We now have two distinct dynamical systems, the original KM and the complex-valued system with the explicit solution in Eq. (5). In the complex-valued system, x ∈ C N has complex-valued elements x i (t) ∈ C whose argument we compare with the numerical solution of the original KM θ i (t) ∈ R (video 1 in the supplementary material). In more detail, let ψ = ψ re + i ψ im be the decomposition of ψ into the real and imaginary parts. Then, we have From this, we can see that the argument of x(t) is the real part of ψ. With this, we can show that the real and imaginary parts of the complex-valued model follow the dynamics (see Sec. I B in the supplementary material) Note that Eqs. (8) and (9) show that the complex-valued system is not identical to the original KM since moduli |x i (t)| in the complexvalued system can exhibit amplitude dynamics that increase in the complex plane. Surprisingly, however, when initialized with unitmodulus initial conditions |x i (0)| = 1 for all i, with complex arguments Arg[x i (0)] that match the initial phases θ i (0) in the original KM, the trajectories in the original and complex-valued KM correspond for a non-trivial window of time (Fig. S1). In this work, we introduce an approach to evaluate Eq. (5) in short time windows, with the goal of utilizing the analytical form to generate insight into the underlying mechanism of synchronization, chimeras, and traveling waves in nonlinear oscillator networks. The approach involves (i) taking a unit-modulus state vector x(t 0 ) at the beginning of each window, (ii) propagating the analytical expression by evaluating the matrix exponential in Eq. (1), "original KM"] throughout the rest of this work. The numerical simulations were performed using Euler's method, and no significant differences were observed using a different integration technique. Importantly, we note that initial conditions (at t = 0) in the complex-valued Kuramoto model always have unit modulus |x i (0)| = 1 for all i, with arguments Arg[x i (0)] equal to the initial angles θ i (0) in the original KM.
In general (see supplementary material Sec. I B for specific conditions), we can write the temporal dynamics of x in terms of the Chaos ARTICLE scitation.org/journal/cha eigenvectors and eigenvalues of K, where λ k is the eigenvalue associated with the eigenvector v k . The projection on to the kth eigenmode is given by the inner product This expression now allows us to understand the three major synchronization phenomena previously discovered in the Kuramoto model from a comprehensive mathematical perspective, in terms of the geometry of µ k representing the contribution of kth eigenmode. Furthermore, when A follows the same connectivity rule for all nodes (i.e., where connections can be fully specified by a single vector cyclically permuted across all nodes, resulting in a circulant matrix), the circulant diagonalization theorem (CDT) allows obtaining the eigenspectrum analytically, with eigenvectors following an ordering of Fourier modes from low to high spatial frequencies (supplementary material, Sec. I B, Figs. S10 and S13). 24 Using Eq. (10), we can then analyze the system in terms of the eigenmode contributions µ k .
We first studied the emergence of phase synchronized states in Kuramoto oscillators with attractive coupling (φ = 0), where the adjacency matrix A is given by a complete graph on N = 50 nodes. Figure 1 shows the precise correspondence between (a) the original Kuramoto and (b) the analytical evaluation during the transition from random initial conditions ( to synchrony. This transition is measured by the Kuramoto order parameter Prior to the transition to synchrony, the eigenmode contributions are almost uniform. After the transition, however, the contributions shift away from uniformity, and the first eigenmode becomes dominant (yellow line at µ 1 ).
The analytical approach we have introduced provides an opportunity to understand this transition geometrically. Specifically, we can study how the first eigenmode, which represents the synchronized state, behaves in relation to the other modes, which represent more complicated configurations of phase. As the system transitions to synchronization, the contribution of the first eigenmode increases in magnitude (modulus) in the complex plane, while the other eigenmode contributions collapse to the origin (video 2 in the supplementary material). Because it is the argument of Eq. (5) that determines the phase dynamics in the analytical approach, the contribution of the first eigenmode, which is associated with the eigenvector , explicitly brings the network toward the synchronized state. This analysis provides a novel geometric insight into the transition to synchrony and how the pattern of connections in a network determines the transition from an incoherent to a coherent state during an individual simulation. Furthermore, when φ = π/2, interactions in the network are no longer attractive, and in this case, Kuramoto networks are plotted in color-code (dark tones indicate values close to −π and light tones to π). As time evolves, the spatiotemporal dynamics become coherent and phase synchronization appears. A quantitative comparison between the two models is provided in Fig. S2. The level of synchronization, measured by the Kuramoto order parameter (c), starts at a low value with the random initial conditions and quickly approaches unity, which indicates the phase synchronized state of the system. The eigenmode contribution (log |µ|) is plotted in color (d). The contribution of the first eigenmode, which represents zero phase difference across oscillators, dominates when the phase synchronized state appears.
do not synchronize. 25 Consistent with this, our analysis shows that the eigenmode contributions remain uniform across time ( Fig. S5 and video 3 in the supplementary material). These results demonstrate that the geometric view provided by the eigenmodes can provide insight into the transient dynamics of synchronized and incoherent states.
We now use this approach to understand two more sophisticated dynamical phenomena that have been discovered in the Kuramoto model: partially synchronized chimera states and traveling waves. Chimera states, which are a sophisticated mixture of synchronized and non-synchronized clusters in oscillator networks, are known to arise in models with distance-dependent (non-local) connections and a constant phase lag (φ). 17,[26][27][28] Here, we consider the case where connections follow a deterministic, distancedependent power rule that specifies a real-valued connection weight (supplementary material, Sec. I E 2). Figure 2 depicts the spatiotemporal dynamics for the original KM and our analytical solution for the cases φ = 1.15 (a) and φ = 1.30 (b). In Fig. 2(a), one can see a transient chimera, where the system transitions from incoherence, to a chimera state, and finally to phase synchronization. In Fig. 2(b), the system transitions from incoherence to a chimera state that continues for all times we simulated. Importantly, in both cases, the spatiotemporal dynamics produced by our analytical evaluation depicts a good correspondence with the chimera states observed in the numerical simulation (videos 4 and 5 in the supplementary material). A quantitative comparison between the analytical and original KM results in this case is provided in Fig. S3. We now use the geometric approach to the Kuramoto dynamics to derive insight into the mechanism underlying chimera states in oscillator networks. Specifically, in networks where the same connectivity rule is applied to each node, the eigenmodes higher than µ 1 take the form of traveling waves from low to high spatial scales in their arguments (Fig. S13). 24 Analysis of the eigenmode contributions (log |µ|), plotted as solid lines representing an average over 100 simulations (with shaded regions representing the standard deviation over simulations), reveals that the chimera is created by an interplay of the synchronizing mode µ 1 and a set of modes representing waves traveling in opposite directions with progressively higher spatial frequencies (t = 10 s, blue line, Fig. 3, right panel). This is in contrast to the completely synchronous state, where the contribution of the first eigenmode and the rest differ by more than 12 orders of magnitude 10 s into the simulation (burgundy line, Fig. 3, middle panel), and in contrast to the incoherent case (φ = π/2), where the contributions across eigenmodes remain uniform (yellow line, Fig. 3, right panel). These results allow us to understand the emergence of the chimera as an interplay between specific types of modes in the Kuramoto system. Furthermore, analysis of the aggregate connectivity matrix K illustrates that the underlying mechanism for this interplay of modes is a rotation of the eigenvalues in the complex plane, which reduces the difference between the real part of the eigenvalue associated with the synchronized state and the real part of the rest (Figs. S10 and S11).
To illustrate the insight into the chimera state this geometric approach can provide, we can now rewrite the analytical solution using a subset of contributing modes and compare this truncated approximation to the numerical simulation [ Fig. 4(a), see also Fig. S12]. With only the first 10 modes in the solution ({µ 1 , µ 2 , . . . , µ 10 }), the synchronization in the system is overestimated ( Fig. 4(b)). Moreover, with only the last 10 modes in the solution ({µ 216 , µ 217 , . . . , µ 225 }), the synchronization is underestimated [ Fig. 4(d)], leaving the system with no signature of the chimera. With the first and last 10 modes ({µ 1 , µ 2 , . . . , µ 10 , µ 216 , µ 217 , . . . , µ 225 }), however, the dynamics in the numerical simulation are recovered, and the main structures defining the chimera state are observed in the analytical approximation [ Fig. 4(c)]. These results demonstrate this analytical approach can capture highly non-trivial dynamical phenomena such as chimera states that emerge from the network structure and nonlinear dynamics of the Kuramoto model. Finally, we considered traveling waves in the Kuramoto model. While, when φ = 0, the synchronized state occurs starting from a broad range of random initial conditions (as in Fig. 1), Kuramoto networks can also exhibit traveling wave states that exist for arbitrary lengths of time [original KM and analytical, Fig. 5(a) and 5(b). A quantitative comparison between these two results is provided in Fig. S4]. While the stability of such traveling wave, or "twisted," states has been the subject of much investigation, 29-32 our geometric approach provides insight into transient wave dynamics, which have recently been observed to play important roles in several fields, from neural systems 33,34 to ecology 35 and power grids. 36 Analysis of eigenmode contributions during a traveling wave on a ring graph reveals that these states exist as contributions from a single eigenmode representing a more structured configuration of phase than µ 1 (Figs. S14 and S15).
If the network receives a perturbation with sufficient magnitude, however, the traveling wave state can transition to synchrony [ Fig. 5(c)]. In this case, our analytical approach also captures the timecourse of these transient dynamics well (see video 6 and Fig. S15 in the supplementary material). In this case, the transition to synchrony has a specific signature: while the eigenmode contribution is localized to µ 2 prior to the perturbation, the contributions become spread out during the transition before collapsing to µ 1 [ Fig. 5(d)],

FIG. 5.
Original KM (a) and analytical evaluation (b) for an example of a wave state on a ring graph (N = 100, k = 1, see Sec. I E 2 in the supplementary material), resulting from a special set of initial conditions, a "twisted" state [see Eq. (27) in the supplementary material]. In this context, traveling wave states constitute the contribution of a single eigenmode to the dynamics (2 nd eigenmode, in this case) (cf. Fig. S13). The Kuramoto order parameter (c) in the traveling wave state remains at R = 0. However, when a finite perturbation is applied to the system at t = 2 s, the network transitions to a phase synchronized state (R = 1). The relative instantaneous frequencies following the perturbation exhibit a non-trivial self-organized state during the transition to synchrony (c, inset). This transition is captured by the eigenmode contributions (d), where the 2 nd eigenmode gives way to the 1 st eigenmode when phase synchronization is reached.
since the system transitions from a wave state to phase synchronization [where R = 1 in Fig. 5(c)]. Furthermore, after the perturbation is applied, the instantaneous frequency is no longer the same across oscillators, which results in a non-trivial, self-organized "falcon" shape in the instantaneous frequencies leading to phasesynchronization after the transition [inset, Fig. 5(c)]. These results demonstrate that the geometric view can provide novel insight into transient dynamics in Kuramoto systems following perturbation.
The connection between network structure and resulting nonlinear dynamics is a central question in modern physics. In this work, we have studied an analytical approach to the Kuramoto model, a canonical model for synchronization dynamics in nature. The geometric view enabled by our analytical approach provides insight into three major synchronization phenomena in the Kuramoto model (phase synchronization, partially synchronized chimera states, and traveling waves). The key novel feature of our analytical approach is that it is valid at finite scales and for individual realizations of the model, which can provide more detailed insight into specific, moment-by-moment dynamics in the system than statistical or asymptotic theoretical approaches. Here, the insight provided by this approach allows us to explain the fully synchronized state as a dominant first eigenmode of the complex system (Fig. 1), partially synchronized chimera states as an interplay between modes representing a set of waves traveling in opposite directions (Figs. 2-4), and traveling waves as localizations in higher eigenmodes (Fig. 5). This unifying insight into three main Chaos ARTICLE scitation.org/journal/cha dynamical phenomena that have been discovered in the Kuramoto model demonstrates the utility of the non-asymptotic, algebraic approach to network structure and nonlinear dynamics in this work.
Our approach involves a complex-valued system that admits an explicit analytical solution and whose trajectories correspond to the original, nonlinear KM for a non-trivial window of time. Motivated to see how precisely the trajectories between the two models could match, here we developed an approach to evaluate the analytical expression in Eq. (5) in short time windows. Again, briefly, we start with unit-modulus complex-valued initial conditions, corresponding to the initial phases of the original KM, at the start of each window. We then evaluate the matrix exponential for a short time window (down to 1 ms in the most challenging case of the chimera state). Finally, we take the argument of each element in the resulting state vector (Arg[( x) i (t)] for all i ∈ [1, N]), as always done when comparing the solution of the complex-valued system with the numerical solution of the original KM.
This approach to evaluating the complex-valued system represents an initial step toward an operator-based approach to nonlinear dynamics. Specifically, while it is possible to evaluate the complex-valued system directly, by applying the matrix exponential to produce a solution that agrees with the original KM for a short, but non-trivial, length of time, here we find this windowed approach produces trajectories that precisely match those in the original, nonlinear KM across many different dynamical regimes, including synchronization, chimeras, and traveling waves. The approach taken in this work can be viewed as the combination of two operators: the matrix exponential and Arg. Importantly, while this approach is more complicated than a standard evaluation of the solution to a linear system by using the matrix exponential, this approach allows us to describe the microscopic evolution of the Kuramoto system in terms of a linear operator acting on an instantaneous state. This, in turn, permits analytical insight into the mechanisms underlying states such as chimeras in terms of the spectrum of the adjacency matrix. This approach thus allows us to connect the nonlinear dynamics generated in an individual simulation of a network of Kuramoto oscillators to the specific adjacency matrix for that system. In addition to investigating the mechanisms of these dynamical phenomena further using spectral graph theory, we aim to understand more fully how the nonlinear dynamics in the original KM can be described to the high precision necessary to match a chimera state, which is known to be a chaotic transient, 28 by the iterative application of two simple operators.
Finally, because the Kuramoto model has been extensively studied both as a model for neural dynamics 37 and as a fundamental model for computation, 33,38,39 these results open up several possibilities for studying the connections between network structure, nonlinear dynamics, and computation. The importance of a recurrent network structure is becoming increasingly recognized both in biological 40 and artificial 41 visual processing, and, in recent work, we introduced a theoretical framework for studying how recurrent connections and traveling waves shape computation in the visual system. 34 Understanding how precise network structure can support specific computations in the brain and artificial learning systems through these analytical approaches may lead to a more comprehensive mathematical understanding of neural computation in future work.

SUPPLEMENTARY MATERIAL
See the supplementary material for additional mathematical details about the methodology, models and networks, complementary figures, description of the parameters used in the analyses, and description of the supplementary movies.