This paper considers a recently introduced $D$-dimensional generalized Kuramoto model for many $(N\u226b1)$ 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\u2192\u221e$) 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\u2192\u221e$ 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)\u2264Kc<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\u2192\u221e)$, 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 $Kc\u2248Kc(+)$. Thus, instability-mediated resetting leads to a unique critical point of the $t\u2192\u221e$ 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.

## I. INTRODUCTION

### A. Background

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 model^{1–4}

where $N$ is the number of agents ($i=1,2,\u2026,N$), $\theta i$ is an angle variable that specifies the state of agent $i$, the parameter $K$ characterizes the coupling strength, and $\omega i$ is the natural frequency of agent $i$ ($\theta \u02d9i=\omega i$ in the absence of coupling), where $\omega i$ is typically chosen randomly for each $i$ from some distribution function $g(\omega )$ [$\u222bg(\omega )d\omega =1$]. Because the parameter $\omega 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, $\theta i$ represents the phase angle of an oscillation experienced by agent $i$, while, in contrast, in the fifth example, $\theta 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 $\theta 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, $\sigma 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, $\sigma i(t)$, is taken to evolve according to the real equation^{12–14}

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

and we call $|\rho (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 $\theta i$ or the unit vector $(cos\u2061\theta isin\u2061\theta 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 ($\sigma i$ for agent $i$) is necessarily given by a three-dimensional unit vector. In addition, Ref. 13 has considered the dynamics of the vector $\sigma 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 point^{12} 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 $\sigma i$ in the presence of frozen-in random site disorder (the terms $Wi\sigma 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 $(D\u22121)$ 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.

### B. The rotation rates $Wi$

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 rates^{12} (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(D\u22121)/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, $W\u2192RW$, for any orthogonal matrix $R$ (e.g., Ref. 16).

### C. The $N\u2192\u221e$ limit and the multiplicity of incoherent equilibria

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

As shown in Sec. II, for $D>2,$ there is an infinite continuum of equilibrium (i.e., time-independent) distribution functions $F$ for which $|\rho |=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 $\delta F$ such that $F+\delta F$ also lies in $M$. Section III is devoted to an analysis of the stability of the manifold $M$; i.e., what happens if $\delta 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$.

### D. The dynamical phase transition

In the case of the original ($D=2$) Kuramoto model, Eq. (1), for $N\u2192\u221e,$ one can consider a distribution function in $(\omega ,\theta )$, i.e., $F(W,\sigma ,t)\u2192f(\omega ,\theta ,t)$. In this $D=2$ case, in contrast to the $D>2$ generalized model, Eq. (2), there is only one $|\rho |=0$ equilibrium distribution function, namely, $f=g(\omega )/(2\pi )$. Furthermore, it has long been well-established for $D=2$, that, as $K$ increases continuously from zero, the long-time ($t\u2192\u221e$) stable value of the order parameter $|\rho |$ undergoes a continuous transition from incoherence ($|\rho |=0$) to partial coherence ($0<|\rho |<1$) as $K$ passes a critical value that depends on $g(\omega )$, 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)

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

Method (ii): Considering the $|\rho |=0$ equilibrium distribution, a linear stability analysis was applied,^{1–4,17} and it was found that the $|\rho |=0$ equilibrium distribution (which exists for *all* $K$) 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 $D\u22653$, 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 $|\rho |$.

### E. Linear stability of the incoherent equilibria

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,\u2026$, 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(+)$.

### F. Instability-mediated resetting (IMR)

These stability results for $D=4,6,\u2026$ 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 $|\rho |$ initially experiences growth consistent with the existence of instability. This growth, however, slows as $|\rho |$ reaches a maximum, and subsequently decays back to zero. But, after this short-lived macroscopic burst, once $|\rho |$ 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$, $|\rho |$ remains near zero until $K$ increases past $Kc(+)$. If $K$ is suddenly increased through $Kc(+)$, there is unstable growth of $|\rho |$, as for when $K$ is increased suddenly through $Kc$, but now $|\rho |$ 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\u2192\u221e$ order parameter vs $K$. This is the IMR phenomenon previously referred to.

### G. Main points of this paper

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, $|\rho (t\u2192\u221e)|$ versus $K$ plotted in Fig. 1 for $D=4,6,\u2026$, are qualitatively similar to the curve for $D=2$, the transient dynamics of $\rho (t)$ starting from a given incoherent distribution at $t=0$ are surprisingly different for even $D\u22654$ 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\u2192\u221e$. 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\u2192\u221e$ 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 $|\rho |$ 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.

## II. INCOHERENT EQUILIBRIA

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 $|\rho |=0$ evolution equation

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

with $\omega 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 $(2k\u22121)th$ and $(2k)th$ elements on the diagonal which are set to $1$. By construction

where $1$ is the $D$-dimensional identity matrix. Setting $\sigma ~=RDT\sigma $ transforms the $|\rho |=0$ evolution equation (5) to

Thus, for each $k$,

is a constant of motion for the $|\rho |=0$ evolution equation $d\sigma /dt=W\sigma $.

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

where $\u2207S\u22c5(vF)$ represents the divergence of the vector field $vF$ on the spherical surface $|\sigma |=1$. Hence, any time independent distribution function, $F0\xaf(W,\sigma )$, for the $|\rho |=0$ dynamics must satisfy

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

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

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

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

for any rotation matrix $R$. Thus

since the constants $Ck$ are invariant to such rotations. Equation (13) for $\rho 0$ then yields $\rho 0=R\rho 0$ for any rotation $R$, which then implies that the integral $\u222b\sigma F0\xafdWd\sigma =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

is Gaussian. That is,

where $gM(w)$ is the Gaussian distribution

Thus,

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

for any $D\xd7D$ rotation matrix $R$.

According to random matrix theory, the distribution of block frequencies $\omega k$ in Eq. (6) for such a Gaussian ensemble of even-dimensional random antisymmetric matrices with $\u27e8w2\u27e9$ set to 1 is^{16}

where $\kappa $ is a constant chosen to ensure that the integral of the distribution $g~(\omega 1,\u2026,\omega 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 $\mu $ is the Haar measure for $D\xd7D$ 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

where

Note that $|\sigma |2=C1+\cdots +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 $\sigma $ is distributed over the $D/2$ subspaces of $W$ that are invariant for the $|\rho |=0$ dynamics of $\sigma $.

## III. STABILITY OF INCOHERENT EQUILIBRIA

We linearize Eq. (2) about distributions corresponding to incoherent equilibria, i.e., $|\rho |=0$, by setting $\sigma =\sigma 0+\delta \sigma $ and $\rho =\delta \rho $ for small perturbations $\delta \sigma $ and $\delta \rho $. This yields

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

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

with

for $1\u2264k\u2264D/2$. We can equivalently represent Eq. (27) as

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

Now, assuming that $\delta \rho (t)=est\delta \rho (0)$, Eq. (25) yields

where $\sigma 0(\tau )=RDQ(t)RDT\sigma 0(0)$.

We note that the order parameter of the perturbed system, $\delta \rho $, will be given by the average of $\delta \sigma $ 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,

where $\u27e8\u2219\u27e9\sigma 0(0)$ denotes an average over $\sigma 0(0)$ at fixed $W$ and $\u27e8\u2219\u27e9W$ denotes an average over $W$. We first average Eq. (31) over $\sigma 0(0)$

We focus on the evaluation of the term

Note that $\sigma 0\sigma 0T$ is a $D\xd7D$ matrix which can be constructed from $(D/2)\xd7(D/2)$ blocks of $2\xd72$ matrices, where the block at index $(k,l)$ will be $xkxlT$ for $1\u2264k,l\u2264D/2$. Defining $xk(0)=(yk+,yk\u2212)T$, we obtain from Eq. (27)

Since $Ck=(yk+)2+(yk\u2212)2$, we write

Thus,

We interpret the average to be performed in Eq. (36) as an average over $\theta k$ and $Ck$ for each $k$, with the differential element $d\sigma $ transforming to $\u220fkCkdCkd\theta k$.

Noting that $\u27e8xk\u27e9$ averaged over $\theta k$ is zero, we see that $\u27e8xkxlT\u27e9$ can only be nonzero if $k=l$. Further, in averaging $xkxkT$, the diagonal terms corresponding to $Ckcos2\u2061(\omega k\tau \u2212\theta k)$ and $Cksin2\u2061(\omega k\tau \u2212\theta k)$ will yield $(Ck/2)$ when averaged over $\theta k$, and the cross terms corresponding to $Cksin\u2061(\omega k\tau \u2212\theta k)cos\u2061(\omega k\tau \u2212\theta k)$ will yield zero. Thus, we obtain

where $12$ represents the $2\xd72$ identity matrix. Note that the average over $\theta k$ removes all $\tau $ dependence in Eq. (36). Performing the average over $Ck$, we obtain

where

with the domain $\Gamma $ corresponding to the set of all $c$ such that $0\u2264Ck\u22641$ for all $k$, and $\u2211kCk=1$.

Thus, the quantity $\u27e8\sigma 0(\tau )\sigma 0(\tau )T\u27e9\sigma 0(0)$ in Eq. (35) becomes

where $C\xaf(W)$ is the $D$-dimensional diagonal matrix,

or

Integrating over $\tau $, we obtain

Using the change of basis Eq. (6),

Since $RD(W)=RD(\u2212W)$, $G(W)=G(\u2212W)$, and by Eq. (15)$C\xaf(W)=C\xaf(\u2212W)$, we can replace the $(s1\u2212W~)\u22121$ term in Eq. (44) by

Noting that

the quantity in Eq. (45) becomes

which when inserted into Eq. (44) yields

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., ${\omega 1,\u2026,\omega D/2}$. Noting that $Trace(V)=Trace(RVRT)$ for any rotation $R$ and that the average $\u27e8RVRT\u27e9R$ over an isotropic ensemble of rotations $R$ must, by the isotropy, be a scalar multiple of the $D\xd7D$ identity matrix, we obtain

Using Eqs. (48) and (21), we find that, for $\delta \rho (0)\u22600$, Eq. (47) yields the scalar equation

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\xafk$, we write $C\xafk\u2261C\xafk(W)$ [see Eq. (42)], we can more specifically write it as a function only of the rotation invariant block frequencies ${\omega 1,\u2026,\omega D/2}$ characterizing $W$

Due to the isotropy of the ensemble of matrices $W$, the function $C\xafk(\omega 1,\u2026,\omega 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

To obtain $Kc$, the critical coupling constant at instability onset, we consider the limit $Re(s)\u21920$ 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+iq\u2192iq$ from $p>0$. Note that

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

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\xaf1$ 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

and

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

For a given distribution $g~(\omega 1,\u2026,\omega D/2)$ and a given $C\xaf1(\omega 1,\u2026,\omega D/2)$, Eq. (56) can be solved to obtain a set of solutions for $q$, which we denote as $Q$. Note that $q=0\u2208Q$. The $q$ dependence of $Kc$ indicates that for each value of $q\u2208Q$, 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,

For notational convenience, we define

Recalling Eq. (42), we see that $C\xaf1$ is the expected value of the fraction of $|\sigma i|2$ lying in the first invariant subspace of $W$. Hence, for $D\u22654$,

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\xaf1=1$. Thus, Eq. (54) shows that $Kc(q)$ must lie in the range

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

Using the above inequality, we make the following observations:

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

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 $|\rho |>0$.

For an arbitrary choice of $C\xafk$ 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\xaf1$ 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\xaf1=1$, resulting in a single critical coupling constant $Kc=Kc(+)=2/[\pi h(0)]=2/[\pi 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* $\sigma $: For each $Wi$, the corresponding unit vector $\sigma i$ is chosen randomly with uniform probability in all directions. Thus, the expected value of the magnitude squared of the projection $\sigma iPk\sigma i$ onto subspace $k$ [see Eq. (8)] is the same for all of the $D/2$ subspaces, and, since $|\sigma i|2=1$, this expected value is $(2/D)$, i.e.,

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=0\u2208Q$, thus

giving $Kc(u)=4/[3\pi 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(\u2212)$. 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

For $D=4,$ this corresponds to

Note that for this distribution $C\xaf1(0,\omega 2)=0$ for all $\omega 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/[\pi h(0)]=Kc(\u2212)$.

*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 $\sigma 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

As earlier, the integration of Eq. (54) with the above $C\xaf1$ 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/[\pi h(0)]=Kc(+)$. (An analogous construction of setting $Cmin=1$ for each agent does not work to construct a maximally stable distribution in $D\u22656$. 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 $D\u22656$. 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,\sigma )$ are initialized (e.g., when the uniform $\sigma $ 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 $D\u22654$, as shown via the arrows marked in Fig. 1. Thus, for $D\u22654,$ we obtain

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(\u2212)$, $Kc(u)$, and $Kc(+)$ for the cases of even $D\u22648$ in Table I.

. | $h(0)$ . | $Kc(\u2212)$ . | $Kc(u)$ . | $Kc(+)$ . |
---|---|---|---|---|

$D=2$ | $(2\pi )\u22121/2$ | N/A | 1.596 | 1.596 |

$D=4$ | $(3/4)\xd7(2\pi )\u22121/2$ | 1.064 | 1.418 | 2.128 |

$D=6$ | $(5/8)\xd7(2\pi )\u22121/2$ | 1.277 | 1.532 | 2.553 |

$D=8$ | $(35/64)\xd7(2\pi )\u22121/2$ | 1.459 | 1.667 | 2.918 |

. | $h(0)$ . | $Kc(\u2212)$ . | $Kc(u)$ . | $Kc(+)$ . |
---|---|---|---|---|

$D=2$ | $(2\pi )\u22121/2$ | N/A | 1.596 | 1.596 |

$D=4$ | $(3/4)\xd7(2\pi )\u22121/2$ | 1.064 | 1.418 | 2.128 |

$D=6$ | $(5/8)\xd7(2\pi )\u22121/2$ | 1.277 | 1.532 | 2.553 |

$D=8$ | $(35/64)\xd7(2\pi )\u22121/2$ | 1.459 | 1.667 | 2.918 |

In order to demonstrate that any $Kc$ value between $Kc(\u2212)$ 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 $\sigma 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(\u2212)$ (for the minimally stable case), with

corresponding to

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

## IV. MACROSCOPIC BURSTS AND INSTABILITY-MEDIATED RESETTING

In Sec. III, we considered $N\u2192\u221e$ and showed that, for even $D\u22654$, 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(\u2212)$ to $Kc(+)=2Kc(\u2212)$, where $Kc$ depends on $F$. By definition, for any $K<Kc(\u2212)$, 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, $|\rho |$, characterizing the coherence of the agent population for $t\u2192\u221e$? In particular, in Fig. 1, how is the time-asymptotic value for $|\rho |$ maintained at zero for $Kc(\u2212)<K<Kc(+)$ despite multiple incoherent equilibria losing their stabilities at critical coupling strengths $Kc<K$?

### A. Macroscopic bursts of coherence

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(\u2212)$. We then consider a sudden increase in $K$ to a value satisfying $Kc(\u2212)<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(\u2212)\u22481.064$ (see Table I). In Fig. 3(a), we plot two quantities—in the orange solid curve we present $|\rho (t)|$, and in the blue dashed curve we show the average value of $Cmin$ over all agents, $\u27e8Cmin\u27e9$. Note the rapid rise and fall of $|\rho |$ which is accompanied by a change in value of $\u27e8Cmin\u27e9$. This rapid change in $\u27e8Cmin\u27e9$ indicates the evolution of the system away from the initialized incoherent distribution (constructed to have$\u27e8Cmin\u27e9\u22480$) to a different different incoherent distribution with a larger value of $\u27e8Cmin\u27e9$.

To explain the origin and consequences of this short-lived macroscopic burst of $|\rho |$ 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(\u2212)\u2264Kc<Kc(+)$ (in the numerical example presented above, $F$ was constructed to be a minimally stable distribution with $Kc=Kc(\u2212)$). Denote a distribution of agents for a system initialized close to this incoherent steady-state distribution by $F+\delta F$, for some perturbation $\delta 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 $\delta F$, the distribution $F+\delta 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+\delta F$, at a rate governed by Eq. (49), with the perturbation $\delta F$ increasing as $\delta 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 $|\rho |$ described earlier [orange curve in Fig. 3(a)]. Note, however, that for $K\u2264Kc(+)$, the analysis in Ref. 12 shows that there are no time-independent attractors with $|\rho |>0$, and, further, our numerical experiments indicate that there are no $|\rho |>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 $|\rho |$ 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 $F\u2208M$ to the distribution $F1\u2208M$ with $F1\u2260F$ is not distinguishable solely by observation of the time-asymptotic values of $\rho $, 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 $\rho $ in the form of a rapid short-lived burst of $|\rho |$ away from its steady state value near zero.

### B. Instability-mediated resetting

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 $F\u2208M$ to $F1\u2208M$, 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 $|\rho |=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(+)\u22482.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), $|\rho |$ and $\u27e8Cmin\u27e9$ do not change significantly, whereas in Fig. 3(c) for $K=2.0,$ we see a characteristic short burst of $|\rho (t)|$, accompanied by a change in $\u27e8Cmin\u27e9$. 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 $|\rho |$ 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 $|\rho |$ and hence the maximum value is near zero; for $K>Kc(1)$, the burst in $|\rho |$ results in a large value of this maximum, and this transition from zero indicates a value of $Kc(1)\u22481.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 $\u27e8Cmin\u27e9=0$), the uniform $\sigma $ distribution (setup as described in Sec. III), and the maximally stable distribution (constructed with $\u27e8Cmin\u27e9=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 $|\rho |$ attained during the evolution $t\u22650$. 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 $|\rho (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 $|\rho |$ will be apparent with a larger maximum value of $|\rho |$. Thus, we expect the onset of such transient bursts for the three cases at the theoretically described values $Kc(\u2212)$, $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 $|\rho |>0$ as opposed to the rapid rise and fall described earlier.)

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 $F1\u2208M$.

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 $|\rho (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$|\rho |$ as a function of $K$ and time. We observe $|\rho |$ to show a small fluctuating increase at $K\u2248Kc(\u2212)$, which is sustained until $K\u2248Kc(+)$, after which $|\rho |$ approaches the steady state value of $|\rho |>0$ shown in Fig. 5. If we consider successively slower rates of increase of $K$, the resulting plot of $|\rho |$ as a function of $K$ displays smaller fluctuations from $|\rho |=0$ sustained through $Kc(\u2212)<K<Kc(+)$. In the limit of an infinitely slow rate of increase of $K$, we expect that $|\rho |$ will remain at zero for all $K<Kc(+)$, reproducing Fig. 1.

### C. Resetting due to finite size

So far, we have restricted our analysis to the $N\u2192\u221e$ 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\u2192\u221e$ as small, we can regard the system with $N$ finite but large as being akin to the$N\u2192\u221e$ limit system with small added perturbations. Thus, we might suspect finite, large $N$ to induce a slow secular evolution of the $N\u2192\u221e$ 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 $\u27e8Cmin\u27e9$ as a function of time for evolution of $N=103$ agents initialized to the minimally stable distribution with $\u27e8Cmin\u27e9=0$, evolved with $K=0.9<Kc(\u2212)$. Note that $\u27e8Cmin\u27e9$ undergoes slow growth and eventually asymptotes to a large value of $\u27e8Cmin\u27e9\u22480.7$ at very long times. After this long time, the large value of $\u27e8Cmin\u27e9$ 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(\u2212)<K<Kc(+)$, we do not observe any characteristic burst in the value of $|\rho |$, 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\u2192\u221e$. 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 $\u27e8Cmin\u27e9=0$ and evolved at $K=0.9<Kc(\u2212)$ takes about ten times longer time to reach an asymptotic state for $\u27e8Cmin\u27e9$ (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\u2192\u221e$ limit) and wait for sufficiently long time, then one could continuously increase $K$ without the incoherent distribution becoming unstable until $K$ reaches $Kc(+)$.

## V. CONCLUSIONS

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 $D\u22654,$ 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 $|\rho |=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(\u2212)<Kc<Kc(+)=2Kc(\u2212)$, 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 $|\rho |=0$ manifold $M$ of steady-state distributions as $K$ is increased in the range $[Kc(\u2212),Kc(+)]$, which leave their signatures as short-lived macroscopic bursts in the value of $|\rho |$ (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 $|\rho (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 $|\rho |$ 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.

## ACKNOWLEDGMENTS

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.

## REFERENCES

*D*-dimensional generalized Kuramoto model: Odd

*D*is different,”