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

The collective behavior of nonlinear oscillator networks has been used to study systems ranging from biology to physics. In this context, the Kuramoto model (KM) is of great importance. However, it remains difficult to directly relate the structure of a specific network adjacency matrix to the resulting dynamics in nonlinear systems. Here, we use a complex-valued matrix formulation for the KM, whose argument has a correspondence with the original model. This allows us to obtain an analytical approach for the transient dynamics in individual simulations of Kuramoto oscillator networks. We then apply this to study phase synchronization, chimera states, and waves. Our approach gives us a new, geometric perspective of synchronization phenomena in terms of complex eigenmodes, which, in turn, offers a unified geometry for synchrony, chimera states, and waves in nonlinear oscillator networks.

Networks of nonlinear oscillators have attracted interest across physics, neuroscience, biology, and applied mathematics as models where order can arise without a central coordinator.^1,2 Oscillator networks have been applied to study the behavior of insects,^3,4 patterns of social behavior,^5,6 neural systems,^7,8 and physical systems.^9,10 In this domain, the Kuramoto oscillator has emerged as the central mathematical model for studying these synchronization phenomena.^11–15 While it is quite simple, the Kuramoto model has led to the discovery of new dynamical phenomena well beyond the initial study of nodes evolving to the same phase. Three major dynamical phenomena studied in this model are complete phase synchronization,¹² partially synchronized “chimera” states,^16–19 and traveling waves.²⁰ While much study in mathematics has focused on the existence and stability of these states, the link between network topology (i.e., the structure of connections) and transient network dynamics remains unclear. Here, we report a new mathematical approach to nonlinear oscillator networks that can provide insight into how the structure of an individual network drives synchronization phenomena for any single simulation. We find that this new approach can explain the major synchronization phenomena that have been found in the Kuramoto model, not only in terms of stable, steady state solutions but also in terms of transient behaviors in finite time.

The original Kuramoto model (KM) with phase-lag is given by

where $θi(t)∈R$ is the phase of the $i$th 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 alli∈[1,N]$⁠). The value of $ϕ$ transforms the interaction function from the standard attractive case (⁠$ϕ=0$⁠) to the cosine version with neutral coupling (⁠$ϕ=π/2$⁠).²¹ 

Following recent work,²² we can use a complex-valued, algebraic approach to the Kuramoto model. Starting with Eq. (1), we can change to a rotating coordinate frame²³ and set $ω=0$ without loss of generality. By subtracting an additional imaginary component in the interaction term

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 $xi=eiψi$⁠, we can write

which has the general solution

where $x→=(x1,…,xN)$ 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→∈CN$ has complex-valued elements $xi(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 $|xi(t)|$ in the complex-valued system can exhibit amplitude dynamics that increase in the complex plane. Surprisingly, however, when initialized with unit-modulus initial conditions $|xi(0)|=1for alli$⁠, with complex arguments $Arg[xi(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→(t0)$ at the beginning of each window, (ii) propagating the analytical expression by evaluating the matrix exponential in Eq. (5) for a short time window (down to 1 ms in the most challenging case of the chimera state), and then (iii) taking the argument of the result (⁠$Arg[(x→)i(t)]for alli∈[1,N]$⁠), as done when we compare the solution of the complex-valued system with the numerical solution of the original KM. With this approach, the trajectories in the complex-valued model match those in the original Kuramoto model in a variety of conditions and across realizations over random initial conditions (Figs. S6, S7, S8, and S9). A detailed comparison of dynamics in these two models is provided in the supplementary material (Sec. I B).

We, thus, compare the argument of our analytical expression [Eq. (5), hereafter denoted “analytical”] to computer simulations of the original nonlinear Kuramoto model [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 $|xi(0)|=1for alli$⁠, with arguments $Arg[xi(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 eigenvectors and eigenvalues of $K$⁠,

where $λk$ is the eigenvalue associated with the eigenvector $v→k$⁠. The projection on to the $k$th eigenmode is given by the inner product $μk(t)=⟨x→(t),v→k⟩$⁠.

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 $k$th 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).²⁴ 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 (⁠$θi(0)=Arg[xi(0)]∼U[−π,π]$⁠) to synchrony. This transition is measured by the Kuramoto order parameter $R(t)=N−1|∑j=1Nexp⁡[iθj(t)]|$ (where the phases are obtained from the original KM or complex-valued model), which also exhibits a precise correspondence between the original KM and analytical evaluation [Fig. 1(c)]. Using this analytical approach, we can then analyze the system in terms of the contribution of individual eigenmodes, represented in color in Fig. 1(d) (⁠$log⁡|μ|$⁠). 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 $v→1=(1,1,…,1)T$ (and $Arg[(v→1)i]=0for alli∈[1,N]$⁠), 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 do not synchronize.²⁵ 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–28 Here, we consider the case where connections follow a deterministic, distance-dependent 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).²⁴ 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³⁵ and power grids.³⁶ 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)], 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 phase-synchronization 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 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 alli∈[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,²⁸ by the iterative application of two simple operators.

Finally, because the Kuramoto model has been extensively studied both as a model for neural dynamics³⁷ 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⁴⁰ and artificial⁴¹ 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.³⁴ 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.