This paper considers a recently introduced D-dimensional generalized Kuramoto model for many (N1) interacting agents, in which the agent states are D-dimensional unit vectors. It was previously shown that, for even (but not odd) D, similar to the original Kuramoto model (D=2), there exists a continuous dynamical phase transition from incoherence to coherence of the time asymptotic attracting state (time t) as the coupling parameter K increases through a critical value which we denote Kc(+)>0. We consider this transition from the point of view of the stability of an incoherent state, where an incoherent state is defined as one for which the N distribution function is time-independent and the macroscopic order parameter is zero. In contrast with D=2, for even D>2, there is an infinity of possible incoherent equilibria, each of which becomes unstable with increasing K at a different point K=Kc. Although there are incoherent equilibria for which Kc=Kc(+), there are also incoherent equilibria with a range of possible Kc values below Kc(+), (Kc(+)/2)Kc<Kc(+). How can the possible instability of incoherent states arising at K=Kc<Kc(+) be reconciled with the previous finding that, at large time (t), the state is always incoherent unless K>Kc(+)? We find, for a given incoherent equilibrium, that, if K is rapidly increased from K<Kc to Kc<K<Kc(+), due to the instability, a short, macroscopic burst of coherence is observed, in which the coherence initially grows exponentially, but then reaches a maximum, past which it decays back into incoherence. Furthermore, after this decay, we observe that the equilibrium has been reset to a new equilibrium whose Kc value exceeds that of the increased K. Thus, this process, which we call “Instability-Mediated Resetting,” leads to an increase in the effective Kc with continuously increasing K, until the equilibrium has been effectively set to one for which KcKc(+). Thus, instability-mediated resetting leads to a unique critical point of the t time asymptotic state (K=Kc(+)) in spite of the existence of an infinity of possible pretransition incoherent states.

The dynamical phase transition from incoherence to coherence for a recently proposed, higher-dimensional generalization of the Kuramoto model is examined from the point of view of the stability of the incoherent state. It is found that, due to the higher dimensionality, there is a continuum of different possible pretransition incoherent equilibrium states, each with distinct stability properties. This, in turn, leads to a novel phenomenon, which we call “Instability-Mediated Resetting,” which enables the existence of a unique critical transition point in spite of the infinite continuum of possible pretransition states. In general, these results provide an example illustrating that, for systems with a large number of entities described via a macroscopic variable(s), a degeneracy of microscopic states corresponding to the same macroscopic variable may occur and that signatures of such a degeneracy may be observable in the transient macroscopic system dynamics.

Motivated by a host of applications, much recent research has been focused on efforts aimed at understanding the behavior of large systems of many interacting dynamical agents. An important tool elucidating issues in this general area has been the study of simplified paradigmatic models. A prime example of such a model is the Kuramoto model1–4 

(1)

where N is the number of agents (i=1,2,,N), θi is an angle variable that specifies the state of agent i, the parameter K characterizes the coupling strength, and ωi is the natural frequency of agent i (θ˙i=ωi in the absence of coupling), where ωi is typically chosen randomly for each i from some distribution function g(ω) [g(ω)dω=1]. Because the parameter ωi characterizing the dynamics of each agent i is different for each agent, the agents are said to be heterogeneous. This model and its many generalizations have been used to study a wide variety of applications and phenomena. Examples include synchronously flashing fireflies,5 cellular clocks in the brain,6 Josephson junction circuits,7 pedestrian-induced oscillation of foot bridges,8 and motion direction alignment in large groups of agents (e.g., drones or flocking animals),9–11 among many others. In the first four of these examples, θi represents the phase angle of an oscillation experienced by agent i, while, in contrast, in the fifth example, θi specifies the direction in which agent i moves.

One aspect of the Kuramoto model and is previous generalizations is that the state of agent i is given by the single scalar angle variable θi(t). Recently, a generalization of these models has been introduced in which the state of the agent i is a D-dimensional unit vector, σi(t), thus allowing for more degrees of freedom in the dynamics of the individual agents. In this generalized model, the D-dimensional unit vector, σi(t), is taken to evolve according to the real equation12–14 

(2)

where the D-dimensional vector ρ(t) (to be specified subsequently) is a common field felt by all the agents, and Wi [analogous to ωi in Eq. (1)] is a D×D antisymmetric matrix (WiT=Wi), which we refer to as the rotation rate matrix. Note that for K=0, Eq. (1) becomes σ˙i=Wiσi, which represents a uniform rate of rotation of σi in the D-dimensional space, σi(t)=[exp(Wit)]σi(0), analogous to the action of the frequency ωi in D=2. Dotting Eq. (2) with σi, we obtain d|σi|2/dt=0, as required by our designation of σi as a unit vector. In general, depending on the situation to be modeled, ρ(t) can be chosen in different ways.12,15 In this paper, we focus on the simplest interesting choice

(3)

and we call |ρ(t)|, the “order parameter.” We note that, as shown in Ref. 12, Eqs. (2) and (3) reduce to Eq. (1) for D=2 with

thus justifying Eqs. (2) and (3) as a “generalization” of the Kuramoto model, Eq. (1), to higher dimensionality. One motivation for this generalization is the previously mentioned example of the application of Eq. (1) to model motion alignment in flocks: For D=2 [equivalent to the standard Kuramoto case, Eq. (1)], the direction of agent motion [characterized by the scalar angle θi or the unit vector (cosθisinθi)T] can be described for agents moving along a two-dimensional surface (like the surface of the Earth), while, if the agents are, e.g., moving in three dimensions (as for drones flying in the air), then the direction of an agent’s motion (σi for agent i) is necessarily given by a three-dimensional unit vector. In addition, Ref. 13 has considered the dynamics of the vector σi as characterizing the evolution of the opinions of an individual within a group of interacting individuals as the group evolves towards consensus. Another interesting point12 is that the inter-agent coupling for Eqs. (2) and (3) is the same as that for the classical, mean-field, zero-temperature, Heisenberg model for the evolution of N interacting spin states σi in the presence of frozen-in random site disorder (the terms Wiσi, with Wi randomly chosen).

Based on our previous work (see Ref. 12), we view Eqs. (2) and (3) as the simplest D-dimensional generalization of the Kuramoto model subject to the assumption that the state of any agent is a unit vector. See Sec. IV of Ref. 12 for a generalization, motivated by flocking drones, in which the agents are regarded as D-dimensional extended-body agents whose states of orientation are described by (D1) mutually perpendicular unit vectors. Although the model in Sec. IV of Ref. 12 is quite different from that considered here, Ref. 12 shows that it shares the same qualitative type of transition behavior as Eq. (2). Thus, we conjecture that the model we study in the this paper can provide a general guide to the possible behavior of other related systems.

Equation (2) with a zero rotation rate (Wi=0) or a uniform rotation rate (Wi=W) was introduced in Refs. 13 and 14. The generalization to heterogeneous rotation rates12 (the situation to be considered in this paper) makes Eq. (2) more similar to the original Kuramoto model and widens its range of applicability. In what follows, as in Ref. 12, we assume that the rotation rate matrix Wi is randomly generated for each i, by choosing each of its D(D1)/2 upper triangular matrix elements, wpq(i) (with p<q), independently from a zero-mean, Gaussian distribution function as described in Sec. II. Alternately, we can say that each of the Wi is randomly drawn from the ensemble of random antisymmetric matrices corresponding to the Gaussian distribution. It is important to note that this ensemble is invariant under rotations; i.e., the ensemble is unchanged when every matrix in the ensemble is subjected to the same rotation, WRW, for any orthogonal matrix R (e.g., Ref. 16).

We are interested in the case where N1, and, to facilitate analysis, we consider the limit N, for which we characterize the system state for dimensionality D by a distribution function F(W,σ,t) such that the fraction of the agents lying in the differential volume element dσdW centered at (σ,W) in σ-W space is F(W,σ,t)dσdW at time t. We define an incoherent equilibrium distribution to be such that F/t=0 and |ρ|=0, where, since we consider the limit N, Eq. (3) is replaced by

(4)

As shown in Sec. II, for D>2, there is an infinite continuum of equilibrium (i.e., time-independent) distribution functions F for which |ρ|=0. We can think of these distributions as defining a manifold M in the space of distribution functions.

Within this manifold, a given F is neutrally stable to a perturbation δF such that F+δF also lies in M. Section III is devoted to an analysis of the stability of the manifold M; i.e., what happens if δF, the perturbation to F, is transverse to M. Before discussing what we find in Sec. III for the case D>2, it is first useful to recall the well-known results for the original Kuramoto model, corresponding to D=2, as well as relevant results from Ref. 12 for D>2.

In the case of the original (D=2) Kuramoto model, Eq. (1), for N, one can consider a distribution function in (ω,θ), i.e., F(W,σ,t)f(ω,θ,t). In this D=2 case, in contrast to the D>2 generalized model, Eq. (2), there is only one |ρ|=0 equilibrium distribution function, namely, f=g(ω)/(2π). Furthermore, it has long been well-established for D=2, that, as K increases continuously from zero, the long-time (t) stable value of the order parameter |ρ| undergoes a continuous transition from incoherence (|ρ|=0) to partial coherence (0<|ρ|<1) as K passes a critical value that depends on g(ω), see the green curve marked by the star symbols in Fig. 1. We denote this critical value by Kc(+). This transition has been studied from two different points of view (see Refs. 1–4)

FIG. 1.

Dynamical phase transition for the generalized Kuramoto model for D=2 (green stars), 4 (orange triangles), 6 (magenta squares), and 8 (blue diamonds) dimensions. The |ρ| values indicated by the plotted markers are obtained by choosing the values of σi(t=0) and Wi for each of the N=105 agents randomly [where the probability distribution of σi(0) is isotropic in direction and that of Wi is as given in Sec. II] and then integrating Eq. (2) from each such initial condition until |ρ(t)| attains a steady value. These steady state values attained appeared to be independent of this choice of initial condition. The theoretical predictions from Ref. 12 for the critical coupling strength, Kc(+), above which stable |ρ|>0 steady states exist are indicated by correspondingly colored vertical arrows on the x-axis.

FIG. 1.

Dynamical phase transition for the generalized Kuramoto model for D=2 (green stars), 4 (orange triangles), 6 (magenta squares), and 8 (blue diamonds) dimensions. The |ρ| values indicated by the plotted markers are obtained by choosing the values of σi(t=0) and Wi for each of the N=105 agents randomly [where the probability distribution of σi(0) is isotropic in direction and that of Wi is as given in Sec. II] and then integrating Eq. (2) from each such initial condition until |ρ(t)| attains a steady value. These steady state values attained appeared to be independent of this choice of initial condition. The theoretical predictions from Ref. 12 for the critical coupling strength, Kc(+), above which stable |ρ|>0 steady states exist are indicated by correspondingly colored vertical arrows on the x-axis.

Close modal

Method (i): It is assumed that f reaches a steady-state distribution (f/t=0) and the resulting nonlinear equation for f is then analytically solved, yielding two possible solutions for the order parameter |ρ|; one has |ρ|=0 and corresponds to f=g(ω)/(2π); the other satisfies a transcendental equation for |ρ| as a function of K involving an integral of the ω-distribution function g. Taking g to be continuous, unimodal, symmetric, and peaked at ω=0, the transcendental root for |ρ| only exists for KKc(+)>0 and gives the |ρ|>0 branch in Fig. 1. In this approach, an analytical result for Kc(+) is obtained by taking the limit |ρ|0+ in the expression for the transcendental branch. This is the approach originally taken by Kuramoto, who then essentially assumed that the |ρ|=0 branch applies for KKc(+), and the |ρ|>0 branch applies for K>Kc(+).

Method (ii): Considering the |ρ|=0 equilibrium distribution, a linear stability analysis was applied,1–4,17 and it was found that the |ρ|=0 equilibrium distribution (which exists for allK) becomes unstable when K increases through a critical value which is the same as that found for Kc(+) by method (i).

Thus, the value of Kc(+) for the original Kuramoto problem can be obtained straightforwardly by following either method (i) or method (ii).

In Ref. 12 using Method (i), previously employed for the original Kuramoto problem, analysis giving the critical transition values for even D were obtained. These values are indicated by the vertical arrows in Fig. 1 and agree well with the plotted numerical results.

Parenthetically, we note that for odd D3, which is not considered in this paper, the transition is qualitatively different from that shown in Fig. 1. Namely, as shown in Ref. 12, when D is odd, as K increases from negative values through zero there is a discontinuous jump in the coherence |ρ|.

Motivated by the results in Fig. 1, in Sec. III, we report results of a stability analysis of the incoherent equilibria for even D greater than two. That is, we attempt an analysis similar to method (ii), previously applied to the original Kuramoto model. We find that the straightforward correspondence that applies for D=2 between the method (i) result for Kc(+) and the method (ii) stability result does not hold for D=4,6,8,, and that the apparent paradox presented by this finding is resolved by a novel phenomenon that we call Instability-Mediated Resetting (IMR).

Specifically, our stability analysis in Sec. III applied to the infinity of possible incoherent equilibrium distributions found in Sec. II, shows that different incoherent equilibria have different stability properties. Considering one such incoherent equilibrium, as K increases, the equilibrium will become unstable as K passes through some value Kc which depends on the specific incoherent equilibrium considered. There are thus many possible values of Kc, in fact we find a continuum of such Kc values spanning a range between (Kc(+)/2) and Kc(+).

These stability results for D=4,6, suggest the following question. How can instability of incoherent equilibrium distributions for K<Kc(+) be reconciled with the numerical results of Fig. 1 and the corresponding method (i) analytical results (the vertical arrows in Fig. 1)? The answer to this question is given in Sec. IV which reports the following results on the nonlinear evolution of the instability found in Sec. III: Considering an incoherent equilibrium which becomes unstable at K=Kc<Kc(+), if one starts with K<Kc and then rapidly increases K to lie in the range Kc<K<Kc(+), the order parameter |ρ| initially experiences growth consistent with the existence of instability. This growth, however, slows as |ρ| reaches a maximum, and subsequently decays back to zero. But, after this short-lived macroscopic burst, once |ρ| returns to essentially zero the resulting incoherent equilibrium is different from that which existed before the instability occurred, and this resulting new incoherent equilibrium loses stability only at a value of the coupling strength between the value that K has been increased to and Kc(+).

In fact, if the initial burst occurred due to a value of K roughly in the middle of the range Kc<K<Kc(+), the resulting equilibrium may be one which loses stability only at Kc(+) itself, i.e., upon further increase of K, |ρ| remains near zero until K increases past Kc(+). If K is suddenly increased through Kc(+), there is unstable growth of |ρ|, as for when K is increased suddenly through Kc, but now |ρ| asymptotically approaches a positive value consistent with Fig. 1 for K>Kc(+). The essential point is that the instability for Kc<K<Kc(+) resets the equilibrium to one which is stable for K<Kc(+) and becomes unstable only when K exceeds Kc(+), consistent with the plot (Fig. 1) of the t order parameter vs K. This is the IMR phenomenon previously referred to.

This paper focuses on the case of even dimensional generalizations of the Kuramoto model of the form Eq. (2). A main message of this paper is that, although the curves, |ρ(t)| versus K plotted in Fig. 1 for D=4,6,, are qualitatively similar to the curve for D=2, the transient dynamics of ρ(t) starting from a given incoherent distribution at t=0 are surprisingly different for even D4 as compared with D=2. We will demonstrate in Sec. II that for even D>2, in contrast to D=2, there is an infinite continuum of incoherent stable equilibria in the limit of N. In Sec. III, we will perform a linear stability analysis of these equilibria, and show that these equilibria have different critical coupling strengths, i.e., values of K beyond which the equilibria are unstable. Further, we also show that in a continuous range of K, each value of K corresponds to the critical coupling strength of some incoherent equilibrium. The upper limit of this range corresponds to earlier results for the critical coupling strength for the t macroscopic phase transition of the order parameter shown in Fig. 1. To reconcile these lower values of critical stability coupling strengths for incoherent equilibria, with the phase transition of Fig. 1, we will examine the dynamics of the incoherent equilibria beyond their critical coupling strengths. This examination results in the observation of short-lived macroscopic bursts of |ρ| which lead to the phenomenon of Instability-Mediated Resetting, which we demonstrate and describe in Sec. IV. We also discuss the effect of finite N on the evolution of these incoherent equilibria in Sec. IV.

We reiterate that in this paper, we will only consider the case of even D. For each W, there are D/2 two-dimensional invariant subspaces for the |ρ|=0 evolution equation

(5)

To see this, we define the rotation RD to be a D×D orthogonal matrix that puts W in block-diagonal form

(6)

with ωk real. Furthermore, we define Pk to be the projection operator that projects a D-vector onto the kth invariant subspace of W, i.e., RDTPkRD=Pk~, has all elements zero except for the (2k1)th and (2k)th elements on the diagonal which are set to 1. By construction

where 1 is the D-dimensional identity matrix. Setting σ~=RDTσ transforms the |ρ|=0 evolution equation (5) to

(7)

Thus, for each k,

(8)

is a constant of motion for the |ρ|=0 evolution equation dσ/dt=Wσ.

Since we are interested in the case where the number of agents, N, is large, N1, it is appropriate to simplify the analysis by considering the limit N, in which case the state of the system can be described by a distribution function, F(W,σ,t) satisfying

(9)

where S(vF) represents the divergence of the vector field vF on the spherical surface |σ|=1. Hence, any time independent distribution function, F0¯(W,σ), for the |ρ|=0 dynamics must satisfy

(10)

where S represents the gradient operator on the spherical surface |σ|=1. The first equality follows from the fact that S(Wσ)=0 for W an antisymmetric matrix. Since

(11)

by comparing Eqs. (10) and (11), we see that the most general solution for a time-independent distribution F0¯ is

(12)

where c denotes the (D/2)-vector (C1,,CD/2)T, i.e., F0¯ depends on W and the (D/2) constants of the motion. There are two constraints. The first one is that, since |σ|=1, we have that |c|=1. The second constraint is that

(13)

which is automatically satisfied if, as we henceforth assume, F0¯ is isotropic in the sense that

(14)

for any rotation matrix R. Thus

(15)

since the constants Ck are invariant to such rotations. Equation (13) for ρ0 then yields ρ0=Rρ0 for any rotation R, which then implies that the integral σF0¯dWdσ=0, as required by our definition of an incoherent distribution, Eq. (13).

In our work, we consider the case where the marginal distribution of W expressed in terms of the matrix elements

(16)

is Gaussian. That is,

(17)

where gM(w) is the Gaussian distribution

(18)

Thus,

(19)

Since Trace(WTW)=Trace(W2) is invariant to rotations of W (i.e., WRTWR) and dW=d(RW) [since det(R)=1], we see that G(W) as defined above is isotropic in the sense that

(20)

for any D×D rotation matrix R.

According to random matrix theory, the distribution of block frequencies ωk in Eq. (6) for such a Gaussian ensemble of even-dimensional random antisymmetric matrices with w2 set to 1 is16 

(21)

where κ is a constant chosen to ensure that the integral of the distribution g~(ω1,,ωD/2) is normalized to 1. Note that g~ is symmetric to interchanges of any two of its arguments.

As an aside, we also mention that using Eq. (6), W=RDW~RDT, an alternative representation of G(W)dW is

where μ is the Haar measure for D×D rotation matrices. (The Haar measure for rotation matrices essentially gives a formal rigorous specification of what we loosely refer to as isotropy.18 In what follows, we use our informal notion of “isotropy” and do not invoke Haar measures.)

Returning to the distribution function F0, we define F^0 by

(22)

where

(23)

Note that |σ|2=C1++CD/2=1. Clearly, even with G(W) specified as Gaussian, there is still an infinity of choices for F^0 and hence F0. These choices specify how σ is distributed over the D/2 subspaces of W that are invariant for the |ρ|=0 dynamics of σ.

We linearize Eq. (2) about distributions corresponding to incoherent equilibria, i.e., |ρ|=0, by setting σ=σ0+δσ and ρ=δρ for small perturbations δσ and δρ. This yields

(24)
(25)

Transforming Eq. (24) to the basis that block-diagonalizes W [as in Eq. (6)], we obtain

(26)

Thus, each two-dimensional subspace k will undergo independent rotation with frequencies corresponding to real ωk frequencies of W~. This gives the solution

(27)

where Q(t) is a block diagonal matrix with (D/2) blocks of dimensions 2×2 given by

(28)

with

(29)

for 1kD/2. We can equivalently represent Eq. (27) as

(30)

for each k, where xk(t) is the two-dimensional vector formed by the (2k1) and 2k components of σ~0.

Now, assuming that δρ(t)=estδρ(0), Eq. (25) yields

(31)

where σ0(τ)=RDQ(t)RDTσ0(0).

We note that the order parameter of the perturbed system, δρ, will be given by the average of δσ over each agent (corresponding to an average over all W). We also perform an ensemble average over all choices of initial conditions corresponding to a given incoherent equilibrium characterized by F^0(W,c). Thus,

(32)

where σ0(0) denotes an average over σ0(0) at fixed W and W denotes an average over W. We first average Eq. (31) over σ0(0)

(33)

We focus on the evaluation of the term

(34)
(35)
(36)

Note that σ0σ0T is a D×D matrix which can be constructed from (D/2)×(D/2) blocks of 2×2 matrices, where the block at index (k,l) will be xkxlT for 1k,lD/2. Defining xk(0)=(yk+,yk)T, we obtain from Eq. (27)

(37)

Since Ck=(yk+)2+(yk)2, we write

(38)

Thus,

(39)

We interpret the average to be performed in Eq. (36) as an average over θk and Ck for each k, with the differential element dσ transforming to kCkdCkdθk.

Noting that xk averaged over θk is zero, we see that xkxlT can only be nonzero if k=l. Further, in averaging xkxkT, the diagonal terms corresponding to Ckcos2(ωkτθk) and Cksin2(ωkτθk) will yield (Ck/2) when averaged over θk, and the cross terms corresponding to Cksin(ωkτθk)cos(ωkτθk) will yield zero. Thus, we obtain

(40)

where 12 represents the 2×2 identity matrix. Note that the average over θk removes all τ dependence in Eq. (36). Performing the average over Ck, we obtain

(41)

where

(42)

with the domain Γ corresponding to the set of all c such that 0Ck1 for all k, and kCk=1.

Thus, the quantity σ0(τ)σ0(τ)Tσ0(0) in Eq. (35) becomes

(43)

where C¯(W) is the D-dimensional diagonal matrix,

Now performing the average over W as prescribed in Eq. (32), we obtain from Eqs. (33) and (43)

or

Integrating over τ, we obtain

Using the change of basis Eq. (6),

(44)

Since RD(W)=RD(W), G(W)=G(W), and by Eq. (15)C¯(W)=C¯(W), we can replace the (s1W~)1 term in Eq. (44) by

(45)

Noting that

the quantity in Eq. (45) becomes

(46)

which when inserted into Eq. (44) yields

(47)

where

Noting that G(W) is isotropic in the sense of Eq. (20), we can average RDVRDT (equivalently V) over an isotropic ensemble of rotations and replace dWG(W) by the distribution of the rotation invariant quantities characterizing W, i.e., {ω1,,ωD/2}. Noting that Trace(V)=Trace(RVRT) for any rotation R and that the average RVRTR over an isotropic ensemble of rotations R must, by the isotropy, be a scalar multiple of the D×D identity matrix, we obtain

(48)

Using Eqs. (48) and (21), we find that, for δρ(0)0, Eq. (47) yields the scalar equation

(49)

where after averaging over the ensemble of rotations, we have replaced G(W)dW in Eq. (47) by

with g~ being the distribution of block frequencies [Eq. (21)] corresponding to the distribution G(W). Note that, by the invariance of F^0(W,c) with respect to rotations of W, although in our definition of C¯k, we write C¯kC¯k(W) [see Eq. (42)], we can more specifically write it as a function only of the rotation invariant block frequencies {ω1,,ωD/2} characterizing W

Due to the isotropy of the ensemble of matrices W, the function C¯k(ω1,,ωD/2) will be invariant to any swapping of indices, i.e.,

for all k. Since g~ is also invariant to swapping of its arguments [see Eq. (21)], we obtain

(50)

To obtain Kc, the critical coupling constant at instability onset, we consider the limit Re(s)0 from Re(s)>0. Denoting the real and imaginary parts of s by p and q, respectively, we hence consider the limit of s=p+iqiq from p>0. Note that

(51)
(52)

where δ(x) represents the Dirac delta function at x and PV represents the Cauchy Principal value of the integral over ω1. Thus, we find from Eq. (49) that

(53)

where Kc(q) is the critical coupling strength at which a small perturbation to the distribution F0 begins to have an unstable mode growing as est with Im(s)=q. Given our choice of an isotropic ensemble of rotation matrices W, the functions g~ and C¯1 must be even functions in each of their arguments. Further, since Kc(q) represents a coupling strength, it must be real. Thus, from the real and imaginary parts of Eq. (53), we obtain

(54)

and

(55)

Using Eq. (54), the above expression reduces to

(56)

For a given distribution g~(ω1,,ωD/2) and a given C¯1(ω1,,ωD/2), Eq. (56) can be solved to obtain a set of solutions for q, which we denote as Q. Note that q=0Q. The q dependence of Kc indicates that for each value of qQ, there exists a mode of instability that arises at the corresponding value of Kc(q). However, the critical coupling strength Kc of a distribution F0 is the smallest value of K for which there is an unstable mode. Thus,

(57)

For notational convenience, we define

(58)

Recalling Eq. (42), we see that C¯1 is the expected value of the fraction of |σi|2 lying in the first invariant subspace of W. Hence, for D4,

for all realizations of W. For the case of D=2 (i.e., the standard Kuramoto model), there is only a single frequency associated with W, and hence C¯1=1. Thus, Eq. (54) shows that Kc(q) must lie in the range

(59)

Following the form of g~ given in Eq. (21), we observe that h(q) is maximized at q=0Q. Thus, minimizing each of the three terms in the above inequality over qQ,

(60)

Using the above inequality, we make the following observations:

  • For all incoherent equilibria, the corresponding Kc is greater than Kc(). Thus, any incoherent equilibrium will be stable for coupling strengths K<Kc().

  • There does not exist any incoherent equilibrium distribution whose Kc is greater than Kc(+). Thus, all incoherent equilibria become unstable for coupling strengths K>Kc(+). This is consistent with Fig. 1, where we see that for K>Kc(+), the system attains an equilibria with |ρ|>0.

  • For an arbitrary choice of C¯k it is not necessary that Kc(q) will be minimized at q=0. However, for several of the examples we consider below we will consider simple choices for C¯1 such that the minima will occur at Kc(0).

  • The inequality in Eq. (60) does not have an explicit D dependence. However, as noted above for D=2, C¯1=1, resulting in a single critical coupling constant Kc=Kc(+)=2/[πh(0)]=2/[πg~(0)].

In the subsequent discussion, we will consider the special case of D=4 and give examples of distributions and their corresponding critical coupling strengths for the onset of instability.

Uniformσ: For each Wi, the corresponding unit vector σi is chosen randomly with uniform probability in all directions. Thus, the expected value of the magnitude squared of the projection σiPkσi onto subspace k [see Eq. (8)] is the same for all of the D/2 subspaces, and, since |σi|2=1, this expected value is (2/D), i.e.,

(61)

The uniform distribution is of particular interest because of its ease of implementation in computer simulations and because of the intuitive naturalness of this choice. From Eq. (54), we obtain

and since h(q) is minimized at q=0Q, thus

(62)

giving Kc(u)=4/[3πh(0)] for D=4.

Minimally Stable Distribution: We define a minimally stable distribution to be one whose critical coupling constant for the onset of instability corresponds to the lower bound of Eq. (60), i.e., Kc=Kc(). To construct such a distribution, we initialize each agent arbitrarily but restricted to the subspace that is orthogonal to the subspace corresponding to the smallest absolute value of the frequency, i.e., for each agent, we set Cmin=0, where

(63)

For D=4, this corresponds to

(64)

Note that for this distribution C¯1(0,ω2)=0 for all ω2. To see why this results in a minimally stable distribution, we compute the integral in Eq. (54) and observe that Kc(q) for this distribution is minimized at q=0 [see Fig. 2; for this minimally stable distribution Kc(q) has been labelled as Kc(min)(q), shown in purple]. This gives Kc=Kc(0)=1/[πh(0)]=Kc().

FIG. 2.

Kc(q) vs q for the case of the minimally stable distribution [shown in green, labelled Kc(min)(q)] corresponding to Eq. (64) and the maximally stable distribution [shown in purple, labelled Kc(max)(q)] corresponding to Eq. (65) for D=4. Note that Kc(q) is always bounded by 1/[πh(q)] (shown as the red dashed curve) and 2/[πh(q)] (shown as the orange dashed curve) as indicated in Eq. (59). The critical coupling strength for the onset of instability, Kc for a given distribution is given by the minimum value attained by Kc(q), which for Kc(min)(q) and Kc(max)(q) is at q=0. (Kc(max)(q) appears to be approximately minimized at q2.12, corresponding to a value of Kc(q)2.141. The true minima however is at q=0, corresponding to Kc(q)2.128).

FIG. 2.

Kc(q) vs q for the case of the minimally stable distribution [shown in green, labelled Kc(min)(q)] corresponding to Eq. (64) and the maximally stable distribution [shown in purple, labelled Kc(max)(q)] corresponding to Eq. (65) for D=4. Note that Kc(q) is always bounded by 1/[πh(q)] (shown as the red dashed curve) and 2/[πh(q)] (shown as the orange dashed curve) as indicated in Eq. (59). The critical coupling strength for the onset of instability, Kc for a given distribution is given by the minimum value attained by Kc(q), which for Kc(min)(q) and Kc(max)(q) is at q=0. (Kc(max)(q) appears to be approximately minimized at q2.12, corresponding to a value of Kc(q)2.141. The true minima however is at q=0, corresponding to Kc(q)2.128).

Close modal

Maximally Stable Distribution: We define a maximally stable distribution to be one whose critical coupling constant for the onset of instability corresponds to the upper bound of Eq. (60), i.e., Kc=Kc(+). In D=4, such a distribution can be set up similar to the case of the minimally stable distribution, by choosing the σi to lie entirely in the subspace corresponding to the smallest absolute value of the frequency, i.e., by setting Cmin=1 for each agent. This corresponds to

(65)

As earlier, the integration of Eq. (54) with the above C¯1 results in an expression for Kc(q) which is again minimized at q=0 [see Fig. 2; for this maximally stable distribution Kc(q) has been labelled as Kc(max)(q), shown in green]. This gives Kc=Kc(0)=2/[πh(0)]=Kc(+). (An analogous construction of setting Cmin=1 for each agent does not work to construct a maximally stable distribution in D6. While this implies that Cmin=1 for each agent is not always a maximally stable distribution, it does not imply that there is no such distribution in D6. We leave the construction of such a distribution to future work.)

In addition to yielding an upper bound on Kc, maximally stable distributions are of particular interest because they surprisingly tend to arise naturally in our numerical simulations performed on necessarily finite system size, even when other equilibrium distributions F0(W,σ) are initialized (e.g., when the uniform σ distribution is initialized); see Sec. IV C. Note that it is not necessary for a maximally stable distribution to have Cmin=1 for each agent; for example, the maximally stable distributions attained due to the long-time limit of finite-N-effects as shown in Fig. 6 do not have Cmin=1 for each agent.

The largest possible value of the critical coupling constant, Kc(+), beyond which no stable incoherent equilibria exist, corresponds to the calculation of Kc performed in Ref. 12 for D4, as shown via the arrows marked in Fig. 1. Thus, for D4, we obtain

(66)

In particular, for the choice of the distribution of rotations matrices Eq. (19), Ref. 12, presents an expression for h(0), which we use to give values of Kc(), Kc(u), and Kc(+) for the cases of even D8 in Table I.

TABLE I.

Expressions for h(0) and numerical values of Kc(), Kc(u), and Kc(+) for D=2, 4, 6, and 8. The expression for h(0) is obtained from Ref. 12, and the values of the various critical coupling strengths are obtained from Eq. (66).

h(0)Kc()Kc(u)Kc(+)
D=2 (2π)1/2 N/A 1.596 1.596 
D=4 (3/4)×(2π)1/2 1.064 1.418 2.128 
D=6 (5/8)×(2π)1/2 1.277 1.532 2.553 
D=8 (35/64)×(2π)1/2 1.459 1.667 2.918 
h(0)Kc()Kc(u)Kc(+)
D=2 (2π)1/2 N/A 1.596 1.596 
D=4 (3/4)×(2π)1/2 1.064 1.418 2.128 
D=6 (5/8)×(2π)1/2 1.277 1.532 2.553 
D=8 (35/64)×(2π)1/2 1.459 1.667 2.918 

In order to demonstrate that any Kc value between Kc() and Kc(+) can occur depending on the equilibrium, we consider a particular simple example: For every Wi, in our randomly chosen W-ensemble, we determine σi according to either one of the three protocols specified above with probabilities p(u) (for the uniform case), p(+) (for the maximally stable case), or p() (for the minimally stable case), with

(67)

Using the expected value interpretation of C¯1, we thus obtain from Eqs. (61), (64), and (65)

corresponding to

(68)

Hence, for D=4, by choosing values of p(u,+,), we can construct a distribution to have any given value of Kc between Kc() and Kc(+). (We expect similar constructions to exist for all even D4.) Furthermore, for any given Kc, in the range equations (60), (67), and (68) represent only two constraints on the three parameters, p(u), p(), and p(+). Thus, for each value of Kc in the range Kc()<Kc<Kc(+), there are an infinity of possible choices for p(u), p(+), and p() (i.e., an infinite number of distribution functions) satisfying Eq. (68).

In Sec. III, we considered N and showed that, for even D4, the stability of incoherent equilibria depends on their associated equilibrium distribution function, F. In particular, we observe a range of critical parameter values Kc for instability onset from Kc() to Kc(+)=2Kc(), where Kc depends on F. By definition, for any K<Kc(), all incoherent distributions are stable and for K>Kc(+), all incoherent distributions are unstable.

We now return to the central question posed in Sec. I, namely, how can we reconcile the loss of the stability of an incoherent distribution at a critical coupling value of Kc<Kc(+) with Fig. 1, which shows the attracting value of the magnitude of the order parameter, |ρ|, characterizing the coherence of the agent population for t? In particular, in Fig. 1, how is the time-asymptotic value for |ρ| maintained at zero for Kc()<K<Kc(+) despite multiple incoherent equilibria losing their stabilities at critical coupling strengths Kc<K?

To examine the aforementioned question, we first consider the following setup: We initialize the system to a minimally stable incoherent equilibrium distribution by setting Cmin=0 for each agent (see Sec. III). This initial setup will be invariant to evolution with a coupling strength of K<Kc(). We then consider a sudden increase in K to a value satisfying Kc()<K<Kc(+).

The dynamics observed following this change of K is represented in Fig. 3(a), for a numerical simulation of N=106 agents in D=4 dimensions, initialized to the minimally stable distribution with Cmin=0, and then numerically integrated according to Eq. (2) with a coupling strength of K=1.4>Kc()1.064 (see Table I). In Fig. 3(a), we plot two quantities—in the orange solid curve we present |ρ(t)|, and in the blue dashed curve we show the average value of Cmin over all agents, Cmin. Note the rapid rise and fall of |ρ| which is accompanied by a change in value of Cmin. This rapid change in Cmin indicates the evolution of the system away from the initialized incoherent distribution (constructed to haveCmin0) to a different different incoherent distribution with a larger value of Cmin.

FIG. 3.

Representative plots demonstrating the short-lived macroscopic burst of coherence and the resulting IMR. (a) The magnitude of the order parameter (orange solid curve) and Cmin (blue dashed curve) as a function of time for a system setup with a minimally stable distribution corresponding to Cmin=0 and evolved with K=1.4. Note the sharp rise and fall of |ρ|, i.e., the macroscopic burst of coherence, accompanied by the increase of the value of Cmin (i.e., IMR). This results in an increase of the critical coupling constant for instability onset of the new incoherent distribution. Panels (b) and (c) show the order parameter evolution beginning with the distribution function at the last time-step of (a) but with K increased to K=1.6 and 2.0, respectively. The presence of a macroscopic burst of |ρ| in (c) and not in (b) indicates that Kc has been reset to a value between 1.6 and 2.0. In panel (d) |ρ|max indicates the largest value of |ρ| for systems initialized similar to (b) or (c) following a discontinuous increase of the coupling constant to a value K plotted on the horizontal axis. |ρ|max is macroscopically observable (i.e., distinguishable from finite-N-induced fluctuations) for bursts of |ρ|, and approximately 0 for incoherent steady-state distributions without any such burst. Hence (d) indicates that, by the end of the simulation in panel (a), due to IMR the critical coupling strength has been reset to Kc1.75. See text for more details.

FIG. 3.

Representative plots demonstrating the short-lived macroscopic burst of coherence and the resulting IMR. (a) The magnitude of the order parameter (orange solid curve) and Cmin (blue dashed curve) as a function of time for a system setup with a minimally stable distribution corresponding to Cmin=0 and evolved with K=1.4. Note the sharp rise and fall of |ρ|, i.e., the macroscopic burst of coherence, accompanied by the increase of the value of Cmin (i.e., IMR). This results in an increase of the critical coupling constant for instability onset of the new incoherent distribution. Panels (b) and (c) show the order parameter evolution beginning with the distribution function at the last time-step of (a) but with K increased to K=1.6 and 2.0, respectively. The presence of a macroscopic burst of |ρ| in (c) and not in (b) indicates that Kc has been reset to a value between 1.6 and 2.0. In panel (d) |ρ|max indicates the largest value of |ρ| for systems initialized similar to (b) or (c) following a discontinuous increase of the coupling constant to a value K plotted on the horizontal axis. |ρ|max is macroscopically observable (i.e., distinguishable from finite-N-induced fluctuations) for bursts of |ρ|, and approximately 0 for incoherent steady-state distributions without any such burst. Hence (d) indicates that, by the end of the simulation in panel (a), due to IMR the critical coupling strength has been reset to Kc1.75. See text for more details.

Close modal

To explain the origin and consequences of this short-lived macroscopic burst of |ρ| we describe the evolution of the system in the space of distribution functions. Let us consider a given incoherent steady-state distribution, corresponding to a distribution F and having a corresponding critical coupling stability strength Kc for Kc()Kc<Kc(+) (in the numerical example presented above, F was constructed to be a minimally stable distribution with Kc=Kc()). Denote a distribution of agents for a system initialized close to this incoherent steady-state distribution by F+δF, for some perturbation δF. We then examine the expected dynamics for the evolution of the system under the dynamics of Eq. (2) for a coupling strength K abruptly increased from K<Kc to Kc<K<Kc(+).

For almost every perturbation δF, the distribution F+δF will no longer lie in the manifold of incoherent distributions M. Since the initially chosen incoherent distribution is unstable at the increased value of K, for small t, the system will rapidly evolve away from the initial distribution, F+δF, at a rate governed by Eq. (49), with the perturbation δF increasing as δFest, Re(s)>0. This corresponds to increasing distance away from the manifold of incoherent distributions, M, and hence appears as the sharp increase in |ρ| described earlier [orange curve in Fig. 3(a)]. Note, however, that for KKc(+), the analysis in Ref. 12 shows that there are no time-independent attractors with |ρ|>0, and, further, our numerical experiments indicate that there are no |ρ|>0 time-dependent attractors (e.g., periodic or chaotic). Hence, the distribution function must evolve to a stable steady-state distribution function on the manifold M. Thus, in the space of distribution functions, the evolution of the system will follow a trajectory that begins near the initial incoherent steady-state distribution in M, moves away from M, and is then attracted back towards M, but to a different incoherent steady-state distribution (corresponding to some distribution F1) that is stable for the chosen coupling strength K. Thus, observing this system at large finite N via the order parameter demonstrates an initially small value of |ρ| near zero, which rapidly rises to a large (macroscopic) value, and then falls back to a small value near zero as depicted in the representative illustration [Fig. 3(a)].

This transition from the distribution FM to the distribution F1M with F1F is not distinguishable solely by observation of the time-asymptotic values of ρ, since both F and F1 correspond to incoherent steady-state distributions. However, a signature of this transition is displayed in the transient dynamics of the macroscopic observable ρ in the form of a rapid short-lived burst of |ρ| away from its steady state value near zero.

An important expected consequence of the above described behavior is an “Instability-Mediated Resetting” of the system stability properties, which we define and describe as follows: The critical coupling constant of F1, denoted Kc(1), is necessarily greater than K. Hence, due to the evolution of the system from FM to F1M, the critical coupling strength of the system has been reset from Kc<K to Kc(1)>K. This change in critical coupling strength without change in the time-asymptotic macroscopic steady-state of the system (i.e., the system is on the manifold M corresponding to |ρ|=0 at the initial distribution and at the asymptotic final distribution) is what we call Instability-Mediated Resetting. To demonstrate this change in critical coupling strength, we choose the resulting distribution at the end of the aforementioned simulation [corresponding to time t=500 in Fig. 3(a)] as the initial distribution for the following two situations: (i) evolution with K=1.6<Kc(+)2.128, corresponding to Fig. 3(b), and (ii) evolution with K=2.0<Kc(+), corresponding to Fig. 3(c). Note that in Fig. 3(b), |ρ| and Cmin do not change significantly, whereas in Fig. 3(c) for K=2.0, we see a characteristic short burst of |ρ(t)|, accompanied by a change in Cmin. Thus, we infer that 1.6<Kc(1)<2.0, hence indicating this instability-mediated resetting of the critical coupling constant for instability. To more precisely pin down the value of Kc(1), we evolve the system for a range of values of K with each evolution having the initial condition described earlier. In Fig. 3(d), we plot the maximum value of |ρ| attained during the evolution as a function of K. We interpret Fig. 3(d) as follows: for all values of K<Kc(1), there is no burst in |ρ| and hence the maximum value is near zero; for K>Kc(1), the burst in |ρ| results in a large value of this maximum, and this transition from zero indicates a value of Kc(1)1.75.

We use a similar setup to verify Eq. (66). We consider three series of numerical simulations, corresponding to initial conditions of the minimally stable distribution (constructed with Cmin=0), the uniform σ distribution (setup as described in Sec. III), and the maximally stable distribution (constructed with Cmin=1). For each initial condition, we evolve the system with an abrupt increase from a coupling strength less than 0.5 at t=0 to a given value of K and note the maximum value of |ρ| attained during the evolution t0. This is then repeated for the same initial condition with a different value of K, over a range of values for K. The results are then plotted for this maximum attained value of |ρ(t)| as a function of K. As earlier, for K below the corresponding Kc, this maximum value will be approximately zero, and for K above Kc, the rapid macroscopic burst of |ρ| will be apparent with a larger maximum value of |ρ|. Thus, we expect the onset of such transient bursts for the three cases at the theoretically described values Kc(), Kc(u), and Kc(+), respectively, according to Eq. (66). These values have been marked with the vertical dashed lines in Fig. 4. Note the close agreement between these theoretically predicted values and the numerically observed burst onset. We expect improving agreement with increasing N. (Note that for K>Kc(+), the maximum attained value corresponds to the stable distribution with |ρ|>0 as opposed to the rapid rise and fall described earlier.)

FIG. 4.

Transitions demonstrating the results derived in Eq. (66) for a system with N=106. For the minimally stable distribution (red circles), the uniform σ distribution (blue triangles), and the maximally stable distribution (green stars), the system is evolved for various values of K. The maximum value attained by |ρ(t)| over a short evolution is shown as a function of K. For incoherent steady-state distributions that undergo stable evolution at a given value of K, |ρ|max is approximately zero, whereas instability of incoherent steady-state distributions results in a short-lived burst of coherence, resulting in a larger value of |ρ|max. The theoretical predictions for the transitions to instability are shown in the respective colors using vertical dashed lines, and agree well with the numerical results. (Note that for K>Kc(+), |ρ|max corresponds to the stable state of |ρ|>0 shown in Fig. 1 as opposed to the peak value during these short bursts.)

FIG. 4.

Transitions demonstrating the results derived in Eq. (66) for a system with N=106. For the minimally stable distribution (red circles), the uniform σ distribution (blue triangles), and the maximally stable distribution (green stars), the system is evolved for various values of K. The maximum value attained by |ρ(t)| over a short evolution is shown as a function of K. For incoherent steady-state distributions that undergo stable evolution at a given value of K, |ρ|max is approximately zero, whereas instability of incoherent steady-state distributions results in a short-lived burst of coherence, resulting in a larger value of |ρ|max. The theoretical predictions for the transitions to instability are shown in the respective colors using vertical dashed lines, and agree well with the numerical results. (Note that for K>Kc(+), |ρ|max corresponds to the stable state of |ρ|>0 shown in Fig. 1 as opposed to the peak value during these short bursts.)

Close modal

In each of the above cases, for a system initialized to a distribution F, with a corresponding critical instability coupling strength of Kc, we examined the case of an abrupt increase in K from a value of K<Kc to a value K>Kc. The distribution F remains invariant to evolution for K<Kc, and then, after the abrupt increase, there is an initial repulsion away from the state with distribution F, followed by an attraction back towards an invariant state with distribution F1M.

While the system state is away from M and is being attracted towards F1, if the value of K is altered again to one greater than Kc(1), then the system will again be repelled away from M. As the system is then attracted towards another distinct distribution F2 (with a critical coupling strength of Kc(2)>Kc(1)), the coupling strength can be varied again to K>Kc(2), resulting in additional delay in the attraction towards M. In this fashion, if we consider a slowly increasing coupling strength, then we can delay this attraction towards the manifold M for large amounts of time, resulting in |ρ(t)|>0 for extended periods of time without any such steady state existing at the corresponding coupling strength. This phenomenon is demonstrated in Fig. 5, where we consider a linearly increasing coupling constant K, and plot|ρ| as a function of K and time. We observe |ρ| to show a small fluctuating increase at KKc(), which is sustained until KKc(+), after which |ρ| approaches the steady state value of |ρ|>0 shown in Fig. 5. If we consider successively slower rates of increase of K, the resulting plot of |ρ| as a function of K displays smaller fluctuations from |ρ|=0 sustained through Kc()<K<Kc(+). In the limit of an infinitely slow rate of increase of K, we expect that |ρ| will remain at zero for all K<Kc(+), reproducing Fig. 1.

FIG. 5.

Evolution of |ρ(t)| for a system having N=106 initialized at a minimally stable incoherent steady-state distribution with slow temporally linear increase in K shown in orange, and a sliding average shown in red. The temporal increase of K is linear in time and is indicated by the horizontal axis at the top of the figure panel. The vertical dashed lines correspond to Kc() and Kc(+). For KKc(), the initialized steady-state distribution is stable and hence |ρ| maintains a value close to zero. For Kc()<K<Kc(+), the system demonstrates enhanced fluctuations of |ρ| about increased, nonzero values that are apparently sustained by the continuous increase of K. For KKc(+), no incoherent distribution is stable, and |ρ| attains a larger value similar to Fig. 1.

FIG. 5.

Evolution of |ρ(t)| for a system having N=106 initialized at a minimally stable incoherent steady-state distribution with slow temporally linear increase in K shown in orange, and a sliding average shown in red. The temporal increase of K is linear in time and is indicated by the horizontal axis at the top of the figure panel. The vertical dashed lines correspond to Kc() and Kc(+). For KKc(), the initialized steady-state distribution is stable and hence |ρ| maintains a value close to zero. For Kc()<K<Kc(+), the system demonstrates enhanced fluctuations of |ρ| about increased, nonzero values that are apparently sustained by the continuous increase of K. For KKc(+), no incoherent distribution is stable, and |ρ| attains a larger value similar to Fig. 1.

Close modal
FIG. 6.

Slow finite-N-induced evolution of the incoherent steady-state distributions for (a) N=103 and (b) N=104. Note the significantly longer timescales shown here as compared with the Fig. 3. |ρ| is shown as the orange curve, and Cmin is shown as the blue dashed curve. Further, note the larger timescale for (b) as compared with (a). In both cases, the system was initialized to a minimally stable incoherent steady-state distribution with Cmin=0. |ρ| remains approximately zero, indicating that the system remains in M, but the gradual increase in Cmin indicates the change in distribution on M. The final distribution achieved after long time evolution is a maximally stable distribution.

FIG. 6.

Slow finite-N-induced evolution of the incoherent steady-state distributions for (a) N=103 and (b) N=104. Note the significantly longer timescales shown here as compared with the Fig. 3. |ρ| is shown as the orange curve, and Cmin is shown as the blue dashed curve. Further, note the larger timescale for (b) as compared with (a). In both cases, the system was initialized to a minimally stable incoherent steady-state distribution with Cmin=0. |ρ| remains approximately zero, indicating that the system remains in M, but the gradual increase in Cmin indicates the change in distribution on M. The final distribution achieved after long time evolution is a maximally stable distribution.

Close modal

So far, we have restricted our analysis to the N limit, wherein several incoherent equilibria in the manifold M can simultaneously be stable to perturbations orthogonal to M and are neutrally stable to perturbations in M. Hence, a small perturbation within M can move an incoherent equilibrium in M to another nearby incoherent equilibrium in M, and many such small perturbations can cumulatively cause a large change from an initial incoherent distribution. As we have observed earlier, transient dynamics away from M appear to shift the critical coupling strength for loss of stability towards Kc(+). Thus, we suspect that perturbations away from M are biased towards maximally stable distributions. Since, in practice, N is always finite it is of interest to consider the effect of finite N. Viewing the difference between N finite but large and N as small, we can regard the system with N finite but large as being akin to theN limit system with small added perturbations. Thus, we might suspect finite, large N to induce a slow secular evolution of the N incoherent equilibria towards a maximally stable distribution within M. In particular, we observe that for large-but-finite N, a system initialized at any stable incoherent equilibrium undergoes slow evolution to a equilibrium corresponding to a maximally stable distribution. We demonstrate this effect in Fig. 6(a), where we plot Cmin as a function of time for evolution of N=103 agents initialized to the minimally stable distribution with Cmin=0, evolved with K=0.9<Kc(). Note that Cmin undergoes slow growth and eventually asymptotes to a large value of Cmin0.7 at very long times. After this long time, the large value of Cmin indicates that a large fraction of agents have moved to the subspace corresponding to the lowest frequency of rotation, similar to our setup of the maximally stable distribution in Sec. III. From this state if we consider sudden changes in K to values in the range of Kc()<K<Kc(+), we do not observe any characteristic burst in the value of |ρ|, indicating that our distribution was indeed a maximally stable distribution.

Since this evolution towards a maximally stable distribution appears to be mediated by finite-N effects, we expect this evolution to become progressively slower as N increases, with stationarity of the incoherent distributions restored as N. This picture is confirmed numerically in Fig. 6, where it can be clearly seen that for a larger value of N=104 initialized as earlier with Cmin=0 and evolved at K=0.9<Kc() takes about ten times longer time to reach an asymptotic state for Cmin (Note the different scales on the x-axes of the plots). Thus, for N large but finite, if one were to initialize an incoherent equilibrium distribution with K<Kc (where Kc is calculated in the N limit) and wait for sufficiently long time, then one could continuously increase K without the incoherent distribution becoming unstable until K reaches Kc(+).

In this paper, we look at a D-dimensional generalization of the Kuramoto model.12 Unlike the case of the standard (D=2) Kuramoto model, we have shown that for even D4, there are an infinite number of time-independent distributions of agents (defining the manifold M) that correspond to the completely incoherent distribution (i.e., having |ρ|=0) in the infinite system-size limit (Sec. II). We then proved that these distributions demonstrate different stabilities, with each distribution being stable for coupling strengths below a critical coupling strength Kc corresponding to the given distribution (Sec. III). Further, for each value of Kc within a range Kc()<Kc<Kc(+)=2Kc(), there exists an infinite number of distributions that become unstable as K is increased through Kc. In Sec. IV, we show that these properties result in transitions within the |ρ|=0 manifold M of steady-state distributions as K is increased in the range [Kc(),Kc(+)], which leave their signatures as short-lived macroscopic bursts in the value of |ρ| (Fig. 3). These transitions imply a change in the microscopic state of the system, with the distribution after a transient having a significantly larger critical coupling strength for instability due to an Instability-Mediated Resetting of the distribution function [Fig. 3(d)]. While for all K<Kc(+), the only stable steady-state distributions are on M, we demonstrate (Fig. 5) that considering a linearly increasing K results in a small positive fluctuating value of |ρ(t)| (and hence indicating evolution not on M), which can be sustained for long periods of time as K is linearly increased through the range Kc<K<Kc(+) (where Kc refers to the originally initialized distribution); also, these fluctuations in |ρ| become smaller as the rate of increase of K with time becomes slower. Since there are a multitude of stable distributions on M, with neutral stability to perturbations in M, noise can cause slow evolution of distributions in M. We observe such slow evolution due to noise induced by finite-N effects (Fig. 6), which evolves the system towards a maximally stable distribution.

We thank Michelle Girvan and Thomas M. Antonsen for useful discussion. We also thank the referees for their useful comments. This work was supported by ONR Grant No. N000141512134 and by AFOSR Grant No. FA9550-15-1-0171.

1.
Y.
Kuramoto
, “
Self-entrainment of a population of coupled non-linear oscillators
,” in
International Symposium on Mathematical Problems in Theoretical Physics
(
Springer
,
1975
), pp.
420
422
.
2.
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
,”
Rev. Mod. Phys.
77
,
137
(
2005
).
3.
S. H.
Strogatz
, “
From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators
,”
Physica D
143
,
1
20
(
2000
).
4.
E.
Ott
,
Chaos in Dynamical Systems
(
Cambridge University Press
,
2002
).
5.
J.
Buck
and
E.
Buck
, “
Synchronous fireflies
,”
Sci. Am.
234
,
74
85
(
1976
).
6.
Z.
Lu
,
K.
Klein-Cardeña
,
S.
Lee
,
T. M.
Antonsen
,
M.
Girvan
, and
E.
Ott
, “
Resynchronization of circadian oscillators and the east-west asymmetry of jet-lag
,”
Chaos
26
,
094811
(
2016
).
7.
K.
Wiesenfeld
,
P.
Colet
, and
S. H.
Strogatz
, “
Frequency locking in Josephson arrays: Connection with the Kuramoto model
,”
Phys. Rev. E
57
,
1563
(
1998
).
8.
B.
Eckhardt
,
E.
Ott
,
S. H.
Strogatz
,
D. M.
Abrams
, and
A.
McRobie
, “
Modeling walker synchronization on the millennium bridge
,”
Phys. Rev. E
75
,
021110
(
2007
).
9.
N.
Moshtagh
and
A.
Jadbabaie
, “
Distributed geodesic control laws for flocking of nonholonomic agents
,”
IEEE Trans. Automat. Contr.
52
,
681
686
(
2007
).
10.
J.
Zhu
,
J.
Lu
, and
X.
Yu
, “
Flocking of multi-agent non-holonomic systems with proximity graphs
,”
IEEE Trans. Circuits Syst. I Regular Papers
60
,
199
210
(
2013
).
11.
W.
Wang
and
J.-J. E.
Slotine
, “
On partial contraction analysis for coupled nonlinear oscillators
,”
Biol. Cybern.
92
,
38
53
(
2005
).
12.
S.
Chandra
,
M.
Girvan
, and
E.
Ott
, “Continuous versus discontinuous transitions in the D-dimensional generalized Kuramoto model: Odd D is different,”
Phys. Rev. X
9
,
011002
(
2019
); e-print arXiv:1806.01314.
13.
R.
Olfati-Saber
, “
Swarms on sphere: A programmable swarm with synchronous behaviors like oscillator networks
,” in
Proceedings of the 45th IEEE Conference on Decision and Control
(
IEEE
,
2006
), pp.
5060
5066
.
14.
J.
Zhu
, “
Synchronization of Kuramoto model in a high-dimensional linear space
,”
Phys. Lett. A
377
,
2939
2943
(
2013
).
15.
E.
Ott
and
T. M.
Antonsen
, “
Low dimensional behavior of large systems of globally coupled oscillators
,”
Chaos
18
,
37113
(
2008
).
16.
M.
Mehta
and
N.
Rosenzweig
, “
Distribution laws for the roots of a random antisymmetric Hermitian matrix
,”
Nucl. Phys. A
109
,
449
456
(
1968
).
17.
S. H.
Strogatz
and
R. E.
Mirollo
, “
Stability of incoherence in a population of coupled oscillators
,”
J. Stat. Phys.
63
,
613
635
(
1991
).
18.
J.
Faraut
,
Analysis on Lie Groups: An Introduction, Cambridge Studies in Advanced Mathematics
(
Cambridge University Press
,
2008
).