Observing Microscopic Transitions from Macroscopic Bursts: Instability-Mediated Resetting in the Incoherent Regime of the $D$-dimensional Generalized Kuramoto Model

TL;DR

This paper investigates the stability of incoherent states in a generalized high-dimensional Kuramoto model, revealing a phenomenon called 'Instability-Mediated Resetting' that influences the system's transition to coherence.

Contribution

It introduces the concept of 'Instability-Mediated Resetting' in high-dimensional Kuramoto models, explaining how incoherent equilibria evolve with increasing coupling strength.

Findings

Incoherent equilibria exist for all coupling values within a certain range.

Rapid increases in coupling cause transient coherence bursts that reset the system's stability.

The effective critical coupling approaches the known transition point, unifying the behavior across dimensions.

Abstract

We consider a recently introduced $D$-dimensional generalized Kuramoto model for many $(N\gg 1)$ interacting agents in which the agent states are $D$-dimensional unit vectors. It was previously shown that, for even $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 as the coupling parameter $K$ increases through a critical value $K_c^{(+)}>0$. We consider this transition from the point of view of the stability of an incoherent state, i.e., where the $N\to\infty$ 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=K_c$. We show that there are incoherent equilibria for all $K_c$ within the range $(K_c^{(+)}/2)\leq K_c \leq K_c^{(+)}$. How can the possible instability of incoherent states arising at $K=K_c<K_c^{(+)}$ be reconciled with the previous finding that, at large time, the state is always incoherent unless $K>K_c^{(+)}$? We find, for a given incoherent equilibrium, that, if $K$ is rapidly increased from $K<K_c$ to $K_c<K<K_c^{(+)}$, due to the instability, a short, macroscopic burst of coherence is observed, which initially grows exponentially, but then reaches a maximum, past which it decays back into incoherence. After this decay, we observe that the equilibrium has been reset to one whose $K_c$ value exceeds that of the increased $K$. Thus this process, which we call `Instability-Mediated Resetting,' leads to an increase in the effective $K_c$ with continuously increasing $K$, until the equilibrium has been effectively set to one for which for which $K_c\approx K_c^{(+)}$, leading to a unique critical point of the time asymptotic state.

Tables

Table 1. Table 1: Expressions for h ​ ( 0 ) ℎ 0 h(0) and numerical values of K c ( − ) superscript subscript 𝐾 𝑐 K_{c}^{(-)} , K c ( u ) superscript subscript 𝐾 𝑐 𝑢 K_{c}^{(u)} and K c ( + ) superscript subscript 𝐾 𝑐 K_{c}^{(+)} for D = 2 𝐷 2 D=2 , 4 4 4 , 6 6 6 and 8 8 8 . The expression for h ​ ( 0 ) ℎ 0 h(0) is obtained from Ref. Chandra, Girvan, and Ott, 2018 , and the values of the various critical coupling strengths are obtained from Eq. ( 63 )

	$h\left(0\right)$	$K_c^{(-)}$	$K_c^{\left(u\right)}$	$K_c^{(+)}$
$D=2$	${\left(2\pi \right)}^{-1/2}$	N/A	1.596	1.596
$D=4$	$\left(3/4\right)\times {\left(2\pi \right)}^{-1/2}$	1.064	1.418	2.128
$D=6$	$\left(5/8\right)\times {\left(2\pi \right)}^{-1/2}$	1.277	1.532	2.553
$D=8$	$\left(35/64\right)\times {\left(2\pi \right)}^{-1/2}$	1.459	1.667	2.918

Equations

$$d{\theta }_i/dt={\omega }_i+\frac{K}{N}\sum _{j=1}^N\sin ({\theta }_j-{\theta }_i),$$

$$d{\mathbf{\sigma }}_i/dt=K[\mathbf{\rho }-(\mathbf{\rho }\cdot {\mathbf{\sigma }}_i){\mathbf{\sigma }}_i]+{\mathbf{W}}_i{\mathbf{\sigma }}_i,$$

$$\mathbf{\rho }(t)=\frac{1}{N}\sum _{i=1}^N{\mathbf{\sigma }}_i(t),$$

$${\mathbf{\sigma }}_i=\left(\begin{array}{c}\cos {\theta }_i\\ \sin {\theta }_i\end{array}\right)\text{ and }{\mathbf{W}}_i=\left(\begin{array}{cc}0& -{\omega }_i\\ {\omega }_i& 0\end{array}\right),$$

$$\mathbf{\rho }(t)=\int \int F(\mathbf{W},\mathbf{\sigma },t)\mathbf{\sigma }d\mathbf{\sigma }d\mathbf{W}.$$

$$d\mathbf{\sigma }/dt=\mathbf{W}\mathbf{\sigma }.$$

$${\mathbf{R}}_D^T\mathbf{W}{\mathbf{R}}_D=\tilde{\mathbf{W}}=\left(\begin{array}{ccccccc}0& {\omega }_1& & & & & \\ -{\omega }_1& 0& & & & & \\ & & 0& {\omega }_2& & & \\ & & -{\omega }_2& 0& & & \\ & & & & \ddots & & \\ & & & & & 0& {\omega }_{D/2}\\ & & & & & -{\omega }_{D/2}& 0\end{array}\right),$$

$$1=\sum _{k=1}^{D/2}\widetilde{{\mathbf{P}}_k}=\sum _{k=1}^{D/2}{\mathbf{P}}_k,$$

$$d\tilde{\mathbf{\sigma }}/dt=\tilde{\mathbf{W}}\tilde{\mathbf{\sigma }}.$$

$$C_k=(\widetilde{{\mathbf{P}}_k}\tilde{\mathbf{\sigma }}{)}^T(\widetilde{{\mathbf{P}}_k}\tilde{\mathbf{\sigma }})={\mathbf{\sigma }}^T{\mathbf{P}}_k^T{\mathbf{P}}_k\mathbf{\sigma }={\mathbf{\sigma }}^T{\mathbf{P}}_k\mathbf{\sigma }$$

$$\partial F/\partial t+{\nabla }_S\cdot (\mathbf{v}F)=0,\mathbf{v}=K[1_D-\mathbf{\sigma }{\mathbf{\sigma }}^T]\mathbf{\rho }+\mathbf{W}\mathbf{\sigma },$$

$${\nabla }_S\cdot [(\mathbf{W}\mathbf{\sigma })\overline{F_0}]=(\mathbf{W}\mathbf{\sigma })\cdot {\nabla }_S\overline{F_0}=0,$$

$$0=d{C}_k/dt={\nabla }_S{C}_k(\mathbf{\sigma })\cdot d\mathbf{\sigma }/dt=(\mathbf{W}\mathbf{\sigma })\cdot {\nabla }_S{C}_k(\mathbf{\sigma }),$$

$$\overline{F_0}(\mathbf{W},\mathbf{\sigma })=F_0(\mathbf{W};C_1,\dots ,C_{D/2})=F_0(\mathbf{W},\mathbf{c}),$$

$$0={\mathbf{\rho }}_0=\int {\int }_{|\mathbf{\sigma }|=1}\mathbf{\sigma }\overline{F_0}(\mathbf{W},\mathbf{\sigma })d\mathbf{W}d\mathbf{\sigma },$$

$$\overline{F_0}({\mathbf{R}}^T\mathbf{WR},\mathbf{R}\mathbf{\sigma })=\overline{F_0}(\mathbf{W},\mathbf{\sigma })$$

$$F_0({\mathbf{R}}^T\mathbf{WR},\mathbf{c})=F_0(\mathbf{W},\mathbf{c}),$$

$$w_{ii}=0,w_{ij}=-w_{ji}$$

$$G(\mathbf{W})d\mathbf{W}=\left[{\int }_{|\mathbf{\sigma }|=1}{F}_0(\mathbf{W},\mathbf{c})d\mathbf{\sigma }\right]d\mathbf{W}=\prod _{j=1}^D\prod _{i>j}^D{g}_M(w_{ij})d{w}_{ij},$$

$$g_M(w)=\frac{1}{\sqrt{2\pi ⟨w^2⟩}}{e}^{-\frac{w^2}{2⟨w^2⟩}},⟨w^2⟩={\int }_{-\infty }^{\infty }{w}^2{g}_M(w)dw.$$

$$G(\mathbf{W})=(2\pi ⟨w^2⟩{)}^{-D(D-1)/4}\exp \left[-\mathrm{Trace}({\mathbf{W}}^T\mathbf{W})/(4⟨w^2⟩)\right].$$

$$G(\mathbf{W})=G({\mathbf{R}}^T\mathbf{WR})$$

$$\begin{gathered} \tilde{g}({\omega }_1,\dots ,{\omega }_{D/2}) \\ =\kappa \prod _{1\le j\le k\le D/2}({\omega }_j^2-{\omega }_k^2{)}^2\exp \left(-\sum _{i=1}^{D/2}{\omega }_i^2/2\right), \end{gathered}$$

$$\tilde{g}({\omega }_1,\dots ,{\omega }_{D/2})d{\omega }_1\dots d{\omega }_{D/2}d\mu ({\mathbf{R}}_D),$$

$$F_0(\mathbf{W},\mathbf{c})=G(\mathbf{W}){\hat{F}}_0(\mathbf{W},\mathbf{c}),$$

$$\int {\hat{F}}_0(\mathbf{W},\mathbf{c})\delta (|\mathbf{\sigma }|-1)d\mathbf{\sigma }=1.$$

$$d{\mathbf{\sigma }}_0/dt$$

$$d\delta \mathbf{\sigma }/dt$$

$$d{\tilde{\mathbf{\sigma }}}_0/dt=\tilde{\mathbf{W}}{\tilde{\mathbf{\sigma }}}_0.$$

$${\tilde{\mathbf{\sigma }}}_0(t)=\mathbf{Q}(t){\tilde{\mathbf{\sigma }}}_0(0),$$

$$\mathbf{Q}(t)=\left(\begin{array}{ccccc}{\mathbf{Q}}_1(t)& & & & \\ & \ddots & & & \\ & & {\mathbf{Q}}_k(t)& & \\ & & & \ddots & \\ & & & & {\mathbf{Q}}_{D/2}(t)\end{array}\right),$$

Full text

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

I.1 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 modelKuramoto (1975); Acebrón et al. (2005); Strogatz (2000); Ott (2002),

[TABLE]

where $N$ is the number of agents ($i=1,2,\ldots,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$ ($\dot{\theta}_{i}=\omega_{i}$ in the absence of coupling), where $\omega_{i}$ is typically chosen randomly for each $i$ from some distribution function $g(\omega)$ ($\int g(\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 firefliesBuck and Buck (1976), cellular clocks in the brainLu et al. (2016), Josephson junction circuitsWiesenfeld, Colet, and Strogatz (1998), pedestrian-induced oscillation of foot bridgesEckhardt et al. (2007), and motion direction alignment in large groups of agents (e.g., drones or flocking animals)Moshtagh and Jadbabaie (2007); Zhu, Lu, and Yu (2013); Wang and Slotine (2005), 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, $\bm{\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, $\bm{\sigma}_{i}(t)$, is taken to evolve according to the real equationChandra, Girvan, and Ott (2018); Olfati-Saber (2006); Zhu (2013),

[TABLE]

where the $D$-dimensional vector $\bm{\rho}(t)$ (to be specified subsequently) is a common field felt by all the agents, and $\mathbf{W}_{i}$ (analogous to $\omega_{i}$ in Eq. (1)) is a $D\times D$ antisymmetric matrix ($\mathbf{W}_{i}^{T}=-\mathbf{W}_{i}$) which we refer to as the rotation rate matrix. Note that for $K=0$ Eq. (1) becomes $\dot{\bm{\sigma}}_{i}=\mathbf{W}_{i}\bm{\sigma}_{i}$ which represents a uniform rate of rotation of $\bm{\sigma}_{i}$ in $D$-dimensional space, $\bm{\sigma}_{i}(t)=[\exp(\mathbf{W}_{i}t)]\bm{\sigma}_{i}(0)$, analogous to the action of the frequency $\omega_{i}$ in $D=2$. Dotting Eq. (2) with $\bm{\sigma}_{i}$, we obtain $d|\bm{\sigma}_{i}|^{2}/dt=0$, as required by our designation of $\bm{\sigma}_{i}$ as a unit vector. In general, depending on the situation to be modeled, $\bm{\rho}(t)$ can be chosen in different waysOtt and Antonsen (2008); Chandra, Girvan, and Ott (2018). In this paper we focus on the simplest interesting choice,

[TABLE]

and we call $|\bm{\rho}(t)|$, the ‘order parameter.’ We note that, as shown in Ref.Chandra, Girvan, and Ott, 2018, Eqs. (2) and (3) reduce to Eq. (1) for $D=2$ with

[TABLE]

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\theta_{i}\;\sin\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 ($\bm{\sigma}_{i}$ for agent $i$) is necessarily given by a three-dimensional unit vector. In addition, Ref. Olfati-Saber, 2006 has considered the dynamics of the vector $\bm{\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 pointChandra, Girvan, and Ott (2018) 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 $\bm{\sigma}_{i}$ in the presence of frozen-in random site disorder (the terms $\mathbf{W}_{i}\bm{\sigma}_{i}$, with $\mathbf{W}_{i}$ randomly chosen).

Based on our previous work (see Ref. Chandra, Girvan, and Ott, 2018) 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. Chandra, Girvan, and Ott, 2018 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-1)$ mutually perpendicular unit vectors. Although the model in Sec. IV of Ref. Chandra, Girvan, and Ott, 2018 is quite different from that considered here, Ref. Chandra, Girvan, and Ott, 2018 shows that it shares the same qualitative type of transition behavior as Eq. (2). Thus we conjecture that the model we study in the present paper can provide a general guide to the possible behavior of other related systems.

I.2 The Rotation Rates $\mathbf{W}_{i}$

Equation (2) with zero rotation rate ($\mathbf{W}_{i}=0$) or a uniform rotation rate ($\mathbf{W}_{i}=\mathbf{W}$) was introduced in Refs. Olfati-Saber, 2006; Zhu, 2013. The generalization to heterogeneous rotation ratesChandra, Girvan, and Ott (2018) (the situation to be considered in the present paper) makes Eq. (2) more similar to the original Kuramoto model and widens its range of applicability. In what follows, as in Ref.Chandra, Girvan, and Ott, 2018, we assume that the rotation rate matrix $\mathbf{W}_{i}$ is randomly generated for each $i$, by choosing each of its $D(D-1)/2$ upper triangular matrix elements, $w_{pq}^{(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 $\mathbf{W}_{i}$ 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, $\mathbf{W}\to\mathbf{R}\mathbf{W}$, for any orthogonal matrix $\mathbf{R}$ (e.g., Ref.Mehta and Rosenzweig, 1968).

I.3 The $N\to\infty$ limit and the multiplicity of incoherent equilibria

We are interested in the case where $N\gg 1$, and, to facilitate analysis, we consider the limit $N\to\infty$, for which we characterize the system state for dimensionality $D$ by a distribution function $F(\mathbf{W},\bm{\sigma},t)$ such that the fraction of the agents lying in the differential volume element $d\bm{\sigma}d\mathbf{W}$ centered at $(\bm{\sigma},\mathbf{W})$ in $\bm{\sigma}$-$\mathbf{W}$ space is $F(\mathbf{W},\bm{\sigma},t)d\bm{\sigma}d\mathbf{W}$ at time $t$. Throughout this paper we will use the term state to denote this system state characterized by $F(\mathbf{W},\bm{\sigma},t)$. We define an incoherent equilibrium state to be such that $\partial F/\partial t=0$ and $|\bm{\rho}|=0$, where, since we consider the limit $N\to\infty$, Eq. (3) is replaced by

[TABLE]

As shown in Sec.II, for $D>2$ there is an infinite continuum of equilibrium (i.e., time-independent) distribution functions $F$ for which $|\bm{\rho}|=0$. We can think of these distributions as defining a manifold $\mathcal{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 $\mathcal{M}$. Section III is devoted to an analysis of the stability of the manifold $\mathcal{M}$; i.e., what happens if $\delta F$, the perturbation to $F$, is transverse to $\mathcal{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.Chandra, Girvan, and Ott, 2018 for $D>2$.

I.4 The Dynamical phase transition

In the case of the original ($D=2$) Kuramoto model, Eq. (1), for $N\to\infty$ one can consider a distribution function in $(\omega,\theta)$; i.e., $F(\mathbf{W},\bm{\sigma},t)\to f(\omega,\theta,t)$. In this $D=2$ case, in contrast to the $D>2$ generalized model, Eq. (2), there is only one $|\bm{\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\to\infty$) stable value of the order parameter $|\bm{\rho}|$ undergoes a continuous transition from incoherence ($|\bm{\rho}|=0$) to partial coherence ($0<|\bm{\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 $K_{c}^{(+)}$. This transition has been studied from two different points of view (see Refs.Kuramoto, 1975; Acebrón et al., 2005; Strogatz, 2000; Ott, 2002):

Method (i): It is assumed that $f$ reaches a steady state ($\partial f/\partial t=0$) and the resulting nonlinear equation for $f$ is then analytically solved, yielding two possible solutions for the order parameter $|\bm{\rho}|$; one has $|\bm{\rho}|=0$ and corresponds to $f=g(\omega)/(2\pi)$; the other satisfies a transcendental equation for $|\bm{\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 $|\bm{\rho}|$ only exists for $K\geq K_{c}^{(+)}>0$ and gives the $|\bm{\rho}|>0$ branch in Fig. 1. In this approach, an analytical result for $K_{c}^{(+)}$ is obtained by taking the limit $|\bm{\rho}|\to 0^{+}$ in the expression for the transcendental branch. This is the approach originally taken by Kuramoto, who then essentially assumed that the $|\bm{\rho}|=0$ branch applies for $K\leq K_{c}^{(+)}$, and the $|\bm{\rho}|>0$ branch applies for $K>K_{c}^{(+)}$.

Method (ii): Considering the $|\bm{\rho}|=0$ equilibrium state, a linear stability analysis was appliedKuramoto (1975); Acebrón et al. (2005); Strogatz (2000); Ott (2002); Strogatz and Mirollo (1991), and it was found that the $|\bm{\rho}|=0$ equilibrium state (which exists for all $K$) becomes unstable when $K$ increases through a critical value which is the same as that found for $K_{c}^{(+)}$ by method (i).

Thus the value of $K_{c}^{(+)}$ for the original Kuramoto problem can be obtained straightforwardly by following either method (i) or method (ii).

In Ref.Chandra, Girvan, and Ott, 2018 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\geq 3$, which is not considered in this paper, the transition is qualitatively different from that shown in Fig. 1. Namely, as shown in Ref. Chandra, Girvan, and Ott, 2018, when $D$ is odd, as $K$ increases from negative values through zero there is a discontinuous jump in the coherence $|\bm{\rho}|$.

I.5 Linear Stability of the incoherent state

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 $K_{c}^{(+)}$ and the method (ii) stability result does not hold for $D=4,6,8\ldots$, 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 states 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 $K_{c}$ which depends on the specific incoherent equilibrium considered. There are thus many possible values of $K_{c}$, in fact we find a continuum of such $K_{c}$ values spanning a range between $(K_{c}^{(+)}/2)$ and $K_{c}^{(+)}$.

I.6 Instability-Mediated Resetting (IMR)

These stability results for $D=4,6,\ldots$ suggest the following question. How can instability of incoherent equilibrium states for $K<K_{c}^{(+)}$ 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=K_{c}<K_{c}^{(+)}$, if one starts with $K<K_{c}$ and then rapidly increases $K$ to lie in the range $K_{c}<K<K_{c}^{(+)}$, the order parameter $|\bm{\rho}|$ initially experiences growth consistent with the existence of instability. This growth, however, slows as $|\bm{\rho}|$ reaches a maximum, and subsequently decays back to zero. But, after this short-lived macroscopic burst, once $|\bm{\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 $K_{c}^{(+)}$.

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

I.7 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, $|\bm{\rho}(t\to\infty)|$ versus $K$ plotted in Fig. 1 for $D=4,6,\ldots$, are qualitatively similar to the curve for $D=2$, the transient dynamics of $\bm{\rho}(t)$ starting from a given incoherent distribution at $t=0$ are surprisingly different for even $D\geq 4$ 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\to\infty$. 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\to\infty$ 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 $|\bm{\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 $\mathbf{W}$ there are $D/2$ two-dimensional invariant subspaces for the $|\bm{\rho}|=0$ evolution equation

[TABLE]

To see this, we define the rotation $\mathbf{R}_{D}$ to be a $D\times D$ orthogonal matrix that puts $\mathbf{W}$ in block-diagonal form,

[TABLE]

with $\omega_{k}$ real. Furthermore, we define $\mathbf{P}_{k}$ to be the projection operator that projects a $D$-vector onto the $k^{\text{th}}$ invariant subspace of $\mathbf{W}$, i.e., $\mathbf{R}_{D}^{T}\mathbf{P}_{k}\mathbf{R}_{D}=\widetilde{\mathbf{P}_{k}}$ has all elements zero except for the $(2k-1)^{\text{th}}$ and $(2k)^{\text{th}}$ elements on the diagonal which are set to $1$. By construction

[TABLE]

where $\mathbb{1}$ is the $D$-dimensional identity matrix. Setting $\widetilde{\bm{\sigma}}=\mathbf{R}_{D}^{T}\bm{\sigma}$ transforms the $|\bm{\rho}|=0$ evolution equation Eq. (5) to

[TABLE]

Thus, for each $k$,

[TABLE]

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

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

[TABLE]

where $\bm{\nabla}_{S}\cdot(\mathbf{v}F)$ represents the divergence of the vector field $\mathbf{v}F$ on the spherical surface $|\bm{\sigma}|=1$. Hence any time independent distribution function, $\overline{F_{0}}(\mathbf{W},\bm{\sigma})$, for the $|\bm{\rho}|=0$ dynamics must satisfy

[TABLE]

where $\bm{\nabla}_{S}$ represents the gradient operator on the spherical surface $|\bm{\sigma}|=1$. The first equality follows from the fact that $\bm{\nabla}_{S}\cdot(\mathbf{W}\bm{\sigma})=0$ for $\mathbf{W}$ an antisymmetric matrix. Since

[TABLE]

by comparing Eqs. (10) and (11), we see that the most general solution for a time-independent distribution $\overline{F_{0}}$ is

[TABLE]

where $\mathbf{c}$ denotes the $(D/2)$-vector $(C_{1},\ldots,C_{D/2})^{T}$, i.e., $\overline{F_{0}}$ depends on $\mathbf{W}$ and the $(D/2)$ constants of the motion. There are two constraints. The first one is that, since $|\bm{\sigma}|=1$, we have that $|\mathbf{c}|=1$. The second constraint is that

[TABLE]

which is automatically satisfied if, as we henceforth assume, $\overline{F_{0}}$ is isotropic in the sense that

[TABLE]

for any rotation matrix $\mathbf{R}$. Thus

[TABLE]

since the constants $C_{k}$ are invariant to such rotations. Equation (13) for $\bm{\rho}_{0}$ then yields $\bm{\rho}_{0}=\mathbf{R}\bm{\rho}_{0}$ for any rotation $\mathbf{R}$, which then implies that the integral $\int\bm{\sigma}\overline{F_{0}}d\mathbf{W}d\bm{\sigma}=0$, as required by our definition of an incoherent state, Eq. (13).

In our work we consider the case where the marginal distribution of $\mathbf{W}$ expressed in terms of the matrix elements

[TABLE]

is Gaussian. That is,

[TABLE]

where $g_{M}(w)$ is the Gaussian distribution

[TABLE]

Thus

[TABLE]

Since $\operatorname{Trace}(\mathbf{W}^{T}\mathbf{W})=-\operatorname{Trace}(\mathbf{W}^{2})$ is invariant to rotations of $\mathbf{W}$ (i.e., $\mathbf{W}\to\mathbf{R}^{T}\mathbf{W}\mathbf{R}$) and $d\mathbf{W}=d(\mathbf{R}\mathbf{W})$ (since $\det(\mathbf{R})$=1), we see that $G(\mathbf{W})$ as defined above is isotropic in the sense that

[TABLE]

for any $D\times D$ rotation matrix $\mathbf{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 $\langle w^{2}\rangle$ set to 1 isMehta and Rosenzweig (1968)

[TABLE]

where $\kappa$ is a constant chosen to ensure that the integral of the distribution $\widetilde{g}(\omega_{1},\ldots,\omega_{D/2})$ is normalized to 1. Note that $\widetilde{g}$ is symmetric to interchanges of any two of its arguments.

As an aside, we also mention that using Eq. (6), $\mathbf{W}=\mathbf{R}_{D}\widetilde{\mathbf{W}}\mathbf{R}_{D}^{T}$, an alternative representation of $G(\mathbf{W})d\mathbf{W}$ is

[TABLE]

where $\mu$ is the Haar measure for $D\times D$ rotation matrices. (The Haar measure for rotation matrices essentially gives a formal rigorous specification of what we loosely refer to as isotropyFaraut (2008). In what follows we use our informal notion of ‘isotropy’ and do not invoke Haar measures.)

Returning to the distribution function $F_{0}$, we define $\hat{F}_{0}$ by

[TABLE]

where

[TABLE]

Note that $|\bm{\sigma}|^{2}=C_{1}+\ldots+C_{D/2}=1$. Clearly, even with $G(\mathbf{W})$ specified as Gaussian, there is still an infinity of choices for $\hat{F}_{0}$ and hence $F_{0}$. These choices specify how $\bm{\sigma}$ is distributed over the $D/2$ subspaces of $\mathbf{W}$ that are invariant for the $|\bm{\rho}|=0$ dynamics of $\bm{\sigma}$.

III Stability of incoherent equilibria

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

[TABLE]

Transforming Eq. (24) to the basis that block-diagonalizes $\mathbf{W}$ (as in Eq. (6)), we obtain

[TABLE]

Thus each two-dimensional subspace $k$ will undergo independent rotation with frequencies corresponding to real $\omega_{k}$ frequencies of $\widetilde{\mathbf{W}}$. This gives the solution

[TABLE]

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

[TABLE]

with

[TABLE]

for $1\leq k\leq D/2$. We can equivalently represent Eq. (27) as

[TABLE]

for each $k$, where $\mathbf{x}_{k}(t)$ is the two-dimensional vector formed by the $(2k-1)$ and $2k$ components of $\widetilde{\bm{\sigma}}_{0}$.

Now, assuming that $\delta\bm{\rho}(t)=e^{st}\delta\bm{\rho}(0)$, Eq. (25) yields

[TABLE]

where $\bm{\sigma}_{0}(\tau)=\mathbf{R}_{D}\mathbf{Q}(t)\mathbf{R}_{D}^{T}\bm{\sigma}_{0}(0)$.

We note that the order parameter of the perturbed system, $\delta\bm{\rho}$, will be given by the average of $\delta\bm{\sigma}$ over each agent (corresponding to an average over all $\mathbf{W}$). We also perform an ensemble average over all choices of initial conditions corresponding to a given incoherent equilibrium characterized by $\hat{F}_{0}(\mathbf{W},\mathbf{c})$. Thus

[TABLE]

where \langle\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}\rangle_{\bm{\sigma}_{0}(0)} denotes an average over $\bm{\sigma}_{0}(0)$ at fixed $\mathbf{W}$, and \langle\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}\rangle_{\mathbf{W}} denotes an average over $\mathbf{W}$. We first average Eq. (31) over $\bm{\sigma}_{0}(0)$:

[TABLE]

We focus on the evaluation of the term

[TABLE]

Note that $\bm{\sigma}_{0}\bm{\sigma}_{0}^{T}$ is a $D\times D$ matrix which can be constructed from $(D/2)\times(D/2)$ blocks of $2\times 2$ matrices, where the block at index $(k,l)$ will be $\mathbf{x}_{k}\mathbf{x}_{l}^{T}$ for $1\leq k,l\leq D/2$. Defining $\mathbf{x}_{k}(0)=(y_{k}^{+},y_{k}^{-})^{T}$, we obtain from Eq. (27)

[TABLE]

Since $C_{k}=(y_{k}^{+})^{2}+(y_{k}^{-})^{2}$, we write

[TABLE]

Thus

[TABLE]

We interpret the average to be performed in Eq. (36) as an average over $\theta_{k}$ and $\sqrt{C_{k}}$ for each $k$, with the differential element $d\bm{\sigma}$ transforming to $\prod_{k}\sqrt{C_{k}}d\sqrt{C_{k}}d\theta_{k}$.

Noting that $\langle\mathbf{x}_{k}\rangle$ averaged over $\theta_{k}$ is zero, we see that $\langle\mathbf{x}_{k}\mathbf{x}_{l}^{T}\rangle$ can only be nonzero if $k=l$. Further, in averaging $\mathbf{x}_{k}\mathbf{x}_{k}^{T}$, the diagonal terms corresponding to $C_{k}\cos^{2}(\omega_{k}\tau-\theta_{k})$ and $C_{k}\sin^{2}(\omega_{k}\tau-\theta_{k})$ will yield $(C_{k}/2)$ when averaged over $\theta_{k}$, and the cross terms corresponding to $C_{k}\sin(\omega_{k}\tau-\theta_{k})\cos(\omega_{k}\tau-\theta_{k})$ will yield zero. Thus, we obtain

[TABLE]

where $\mathbb{1}_{2}$ represents the $2\times 2$ identity matrix. Note that the average over $\theta_{k}$ removes all $\tau$ dependence in Eq. (36). Performing the average over $C_{k}$, we obtain

[TABLE]

where

[TABLE]

with the domain $\Gamma$ corresponding to the set of all $\mathbf{c}$ such that $0\leq C_{k}\leq 1$ for all $k$, and $\sum_{k}C_{k}=1$.

Thus the quantity $\langle\bm{\sigma}_{0}(\tau)\bm{\sigma}_{0}(\tau)^{T}\rangle_{\bm{\sigma}_{0}(0)}$ in Eq. (35) becomes

[TABLE]

where $\bar{\mathbf{C}}(\mathbf{W})$ is the $D$-dimensional diagonal matrix,

[TABLE]

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

[TABLE]

or

[TABLE]

Integrating over $\tau$, we obtain

[TABLE]

Using the change of basis Eq. (6),

[TABLE]

Since $\mathbf{R}_{D}(\mathbf{W})=\mathbf{R}_{D}(-\mathbf{W})$, $G(\mathbf{W})=G(-\mathbf{W})$, and by Eq. (15) $\bar{\mathbf{C}}(\mathbf{W})=\bar{\mathbf{C}}(-\mathbf{W})$, we can replace the $(s\mathbb{1}-\widetilde{\mathbf{W}})^{-1}$ term in Eq. (44) by

[TABLE]

Noting that

[TABLE]

the quantity in Eq. (45) becomes

[TABLE]

which when inserted into Eq. (44), yields

[TABLE]

where

[TABLE]

Noting that $G(\mathbf{W})$ is isotropic in the sense of Eq. (20), we can average $\mathbf{R}_{D}\mathbf{V}\mathbf{R}_{D}^{T}$ (equivalently $\mathbf{V}$) over an isotropic ensemble of rotations and replace $d\mathbf{W}G(\mathbf{W})$ by the distribution of the rotation invariant quantities characterizing $\mathbf{W}$, i.e., $\{\omega_{1},\ldots,\omega_{D/2}\}$. Noting that $\operatorname{Trace}(\mathbf{V})=\operatorname{Trace}(\mathbf{R}\mathbf{V}\mathbf{R}^{T})$ for any rotation $\mathbf{R}$ and that the average $\langle\mathbf{R}\mathbf{V}\mathbf{R}^{T}\rangle_{\mathbf{R}}$ over an isotropic ensemble of rotations $\mathbf{R}$ must, by the isotropy, be a scalar multiple of the $D\times D$ identity matrix, we obtain

[TABLE]

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

[TABLE]

where after averaging over the ensemble of rotations, we have replaced $G(\mathbf{W})d\mathbf{W}$ in Eq. (47) by

[TABLE]

with $\widetilde{g}$ being the distribution of block frequencies (Eq.(21)) corresponding to the distribution $G(\mathbf{W})$. Note that, by the invariance of $\hat{F}_{0}(\mathbf{W},\mathbf{c})$ with respect to rotations of $\mathbf{W}$, although in our definition of $\bar{C}_{k}$ we write $\bar{C}_{k}\equiv\bar{C}_{k}(\mathbf{W})$ (see Eq. (42)), we can more specifically write it as a function only of the rotation invariant block frequencies $\{\omega_{1},\ldots,\omega_{D/2}\}$ characterizing $\mathbf{W}$:

[TABLE]

Due to the isotropy of the ensemble of matrices $\mathbf{W}$, the function $\bar{C}_{k}(\omega_{1},\ldots,\omega_{D/2})$ will be invariant to any swapping of indices, i.e.,

[TABLE]

for all $k$. Since $\widetilde{g}$ is also invariant to swapping of its arguments (see Eq. (21)), we obtain

[TABLE]

To obtain $K_{c}$, the critical coupling constant at instability onset, we consider the limit $\operatorname{Re}(s)\to 0$ from $\operatorname{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\to iq$ from $p>0$. Note that

[TABLE]

where $\delta(x)$ represents the Dirac delta function at $x$. Thus we find from Eq. (49) that

[TABLE]

where $K_{c}(q)$ is the critical coupling strength at which a small perturbation to the distribution $F_{0}$ begins to have an unstable mode growing as $e^{st}$ with $\operatorname{Im}(s)=q$. Given that our choice of an isotropic ensemble of rotation matrices $\mathbf{W}$, the functions $\widetilde{g}$ and $\bar{C}_{1}$ must be even functions in each of their arguments. Thus,

[TABLE]

The $q$ dependence of $K_{c}$ indicates that for each value of $q$ there exists a mode of instability that arises at the corresponding value of $K_{c}(q)$. However, the critical coupling strength $K_{c}$ of a distribution $F_{0}$ is the smallest value of $K$ for which there is an unstable mode. Thus

[TABLE]

For notational convenience we define

[TABLE]

Recalling Eq. (42), we see that $\bar{C}_{1}$ is the expected value of the fraction of $|\bm{\sigma}_{i}|^{2}$ lying in the first invariant subspace of $\mathbf{W}$. Hence, for $D\geq 4$,

[TABLE]

for all realizations of $\mathbf{W}$. For the case of $D=2$ (i.e., the standard Kuramoto model) there is only a single frequency associated with $\mathbf{W}$, and hence $\bar{C}_{1}=1$. Thus, Eq. (53) shows that $K_{c}(q)$ must lie in the range

[TABLE]

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

[TABLE]

Using the above inequality we make the following observations:

• •

For all incoherent equilibria, the corresponding $K_{c}$ is greater than $K_{c}^{(-)}$. Thus any incoherent equilibrium will be stable for coupling strengths $K<K_{c}^{(-)}$

• •

There does not exist any incoherent equilibrium distribution whose $K_{c}$ is greater than $K_{c}^{(+)}$. Thus, all incoherent equilibria become unstable for coupling strengths $K>K_{c}^{(+)}$. This is consistent with Fig. 1, where we see that for $K>K_{c}^{(+)}$ the system attains an equilibria with $|\bm{\rho}|>0$.

• •

For an arbitrary choice of $\bar{C}_{k}$ it is not necessary that $K_{c}(q)$ will be minimized at $q=0$. However, for several of the examples we consider below we will consider simple choices for $\bar{C}_{1}$ such that the minima will occur at $K_{c}(0)$.

• •

The inequality in Eq. (57) does not have an explicit $D$ dependence. However, as noted above for $D=2$, $\bar{C}_{1}=1$, resulting in a single critical coupling constant $K_{c}=K_{c}^{(+)}=2/(\pi h(0))=2/(\pi\widetilde{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 $\bm{\sigma}$: For each $\mathbf{W}_{i}$, the corresponding unit vector $\bm{\sigma}_{i}$ is chosen randomly with uniform probability in all directions. Thus the expected value of the magnitude squared of the projection $\bm{\sigma}_{i}\mathbf{P}_{k}\bm{\sigma}_{i}$ onto subspace $k$ (see Eq. (8)) is the same for all of the $D/2$ subspaces, and, since $|\bm{\sigma}_{i}|^{2}=1$, this expected value is $(2/D)$, i.e.,

[TABLE]

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. (53) we obtain

[TABLE]

and hence

[TABLE]

giving $K_{c}^{(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. (57), i.e., $K_{c}=K_{c}^{(-)}$. 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 $C_{min}=0$ where

[TABLE]

For $D=4$ this corresponds to

[TABLE]

Note that for this distribution $\bar{C}_{1}(0,\omega_{2})=0$ for all $\omega_{2}$. To see why this results in a minimally stable distribution we compute the integral in Eq. (53) and observe that $K_{c}(q)$ for this distribution is minimized at $q=0$ (see Fig. 2; For this minimally stable distribution $K_{c}(q)$ has been labelled as $K_{c}^{(min)}(q)$, shown in purple). This gives $K_{c}=K_{c}(0)=1/(\pi h(0))=K_{c}^{(-)}$.

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. (57), i.e., $K_{c}=K_{c}^{(+)}$. In $D=4$, such a distribution can be set up similar to the case of the minimally stable distribution, by choosing the $\bm{\sigma}_{i}$ to lie entirely in the subspace corresponding to the smallest absolute value of the frequency, i.e., by setting $C_{min}=1$ for each agent111An analogous construction of setting $C_{min}=1$ for each agent does not work to construct a maximally stable distribution in $D\geq 6$. While this implies that $C_{min}=1$ for each agent is not always a maximally stable distribution, it does not imply that there is no such distribution in $D\geq 6$. We leave the construction of such a distribution to future work.. This corresponds to

[TABLE]

As earlier, integration of Eq. (53) with the above $\bar{C}_{1}$ results in an expression for $K_{c}(q)$ which is again minimized at $q=0$ (see Fig. ; For this maximally stable distribution $K_{c}(q)$ has been labelled as $K_{c}^{(max)}(q)$, shown in green). This gives $K_{c}=K_{c}(0)=2/(\pi h(0))=K_{c}^{(+)}$.

In addition to yielding an upper bound on $K_{c}$, 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 $F_{0}(\mathbf{W},\bm{\sigma})$ are initialized (e.g., when the uniform $\bm{\sigma}$ distribution is initialized); see Sec. IV.3. Note that it is not necessary for a maximally stable distribution to have $C_{min}=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 $C_{min}=1$ for each agent.

The largest possible value of the critical coupling constant, $K_{c}^{(+)}$, beyond which no stable incoherent equilibria exist, corresponds to the calculation of $K_{c}$ performed in Ref. Chandra, Girvan, and Ott, 2018 for $D\geq 4$, as shown via the arrows marked in Fig.1. Thus, for $D\geq 4$ we obtain

[TABLE]

In particular, for the choice of the distribution of rotations matrices Eq. (19), Ref. Chandra, Girvan, and Ott, 2018, presents an expression for $h(0)$, which we use to give values of $K_{c}^{(-)}$, $K_{c}^{(u)}$ and $K_{c}^{(+)}$ for the cases of even $D\leq 8$ in Table 1.

In order to demonstrate that any $K_{c}$ value between $K_{c}^{(-)}$ and $K_{c}^{(+)}$ can occur depending on the equilibrium, we consider a particular simple example: For every $\mathbf{W}_{i}$ in our randomly chosen $\mathbf{W}$-ensemble, we determine $\bm{\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^{(-)}$ (for the minimally stable case), with

[TABLE]

Using the expected value interpretation of $\bar{C}_{1}$, we thus obtain from Eqs. (58), (61) and (62)

[TABLE]

corresponding to

[TABLE]

Hence for $D=4$, by choosing values of $p^{(u,+,-)}$, we can construct a distribution to have any given value of $K_{c}$ between $K_{c}^{(-)}$ and $K_{c}^{(+)}$ (We expect a similar construction to exist for all even $D\geq 4$). Furthermore, for any given $K_{c}$ in the range Eq. (57), Eqs. (64) and (65) represent only two constraints on the three parameters, $p^{(u)}$, $p^{(-)}$ and $p^{(+)}$. Thus, for each value of $K_{c}$ in the range $K_{c}^{(-)}<K_{c}<K_{c}^{(+)}$ there are an infinity of possible choices for $p^{(u)}$, $p^{(+)}$ and $p^{(-)}$ (i.e., an infinite number of distribution functions) satisfying Eq. (65).

IV Macroscopic bursts and Instability-Mediated Resetting

In the previous section we considered $N\to\infty$ and showed that, for even $D\geq 4$, the stability of incoherent equilibria depends on their associated equilibrium distribution function, $F$. In particular, we observe a range of critical parameter values $K_{c}$ for instability onset from $K_{c}^{(-)}$ to $K_{c}^{(+)}=2K_{c}^{(-)}$, where $K_{c}$ depends on $F$. By definition, for any $K<K_{c}^{(-)}$ all incoherent states are stable, and for $K>K_{c}^{(+)}$ all incoherent states 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 state at a critical coupling value of $K_{c}<K_{c}^{(+)}$ with Fig. 1, which shows the attracting value of the magnitude of the order parameter, $|\bm{\rho}|$, characterizing the coherence of the agent population for $t\to\infty$? In particular, in Fig. 1, how is the time-asymptotic value for $|\bm{\rho}|$ maintained at zero for $K_{c}^{(-)}<K<K_{c}^{(+)}$ despite multiple incoherent equilibria losing their stabilities at critical coupling strengths $K_{c}<K$?

IV.1 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 $C_{min}=0$ for each agent (see Sec.III). This initial setup will be invariant to evolution with a coupling strength of $K<K_{c}^{(-)}$. We then consider a sudden increase in $K$ to a value satisfying $K_{c}^{(-)}<K<K_{c}^{(+)}$.

The dynamics observed following this change of $K$ is represented in Fig. 3(a), for a numerical simulation of $N=10^{6}$ agents in $D=4$ dimensions, initialized to the minimally stable distribution with $C_{min}=0$, and then numerically integrated according to Eq. (2) with a coupling strength of $K=1.4>K_{c}^{(-)}\approx 1.064$ (see Table 1). In Fig. 3(a) we plot two quantities — in the orange solid curve we present $|\bm{\rho}(t)|$, and in the blue dashed curve we show the average value of $C_{min}$ over all agents, $\langle C_{min}\rangle$. Note the rapid rise and fall of $|\bm{\rho}|$ which is accompanied by a change in value of $\langle C_{min}\rangle$. This rapid change in $\langle C_{min}\rangle$ indicates the evolution of the system away from the initialized incoherent distribution (constructed to have $\langle C_{min}\rangle\approx 0$) to a different different incoherent distribution with a larger value of $\langle C_{min}\rangle$.

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

For almost every perturbation $\delta F$, the distribution $F+\delta F$ will no longer lie in the manifold of incoherent states $\mathcal{M}$. Since the initially chosen incoherent state 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 Fe^{st}$, $\operatorname{Re}(s)>0$. This corresponds to increasing distance away from the manifold of incoherent states, $\mathcal{M}$, and hence appears as the sharp increase in $|\bm{\rho}|$ described earlier (orange curve in Fig. 3(a)). Note, however, that for $K\leq K_{c}^{(+)}$ the analysis in Ref.Chandra, Girvan, and Ott, 2018 shows that are no time-independent attractors with $|\bm{\rho}|>0$, and, further, our numerical experiments indicate that there are no $|\bm{\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 $\mathcal{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 in $\mathcal{M}$, moves away from $\mathcal{M}$, and is then attracted back towards $\mathcal{M}$, but to a different incoherent steady state (corresponding to some distribution $F_{1}$) 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 $|\bm{\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\in\mathcal{M}$ to the distribution $F_{1}\in\mathcal{M}$ with $F_{1}\neq F$ is not distinguishable solely by observation of the time-asymptotic values of $\bm{\rho}$, since both distributions correspond to incoherent steady states. However, a signature of this transition is displayed in the transient dynamics of the macroscopic observable $\bm{\rho}$ in the form of a rapid short-lived burst of $|\bm{\rho}|$ away from its steady state value near zero.

IV.2 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 $F_{1}$, denoted $K_{c}^{(1)}$, is necessarily greater than $K$. Hence, due to the evolution of the system from $F\in\mathcal{M}$ to $F_{1}\in\mathcal{M}$ the critical coupling strength of the system has been reset from $K_{c}<K$ to $K_{c}^{(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 $\mathcal{M}$ corresponding to $|\bm{\rho}|=0$ at the initial state and at the asymptotic final state) 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<K_{c}^{(+)}\approx 2.128$, corresponding to Fig. 3(b), and (ii) evolution with $K=2.0<K_{c}^{(+)}$, corresponding to Fig. 3(c). Note that in Fig. 3(b) $|\bm{\rho}|$ and $\langle C_{min}\rangle$ do not change significantly, whereas in Fig. 3(c) for $K=2.0$ we see a characteristic short burst of $|\bm{\rho}(t)|$, accompanied by a change in $\langle C_{min}\rangle$. Thus we infer that $1.6<K_{c}^{(1)}<2.0$, hence indicating this instability-mediated resetting of the critical coupling constant for instability. To more precisely pin down the value of $K_{c}^{(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 $|\bm{\rho}|$ attained during the evolution as a function of $K$. We interpret Fig. 3(d) as follows: For all values of $K<K_{c}^{(1)}$ there is no burst in $|\bm{\rho}|$ and hence the maximum value is near zero; for $K>K_{c}^{(1)}$ the burst in $|\bm{\rho}|$ results in a large value of this maximum, and this transition from zero indicates a value of $K_{c}^{(1)}\approx 1.75$.

We use a similar setup to verify Eq. (63). We consider three series of numerical simulations, corresponding to initial conditions of the minimally stable distribution (constructed with $\langle C_{min}\rangle=0$), the uniform $\bm{\sigma}$ distribution (setup as described in Sec. III), and the maximally stable distribution (constructed with $\langle C_{min}\rangle=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 $|\bm{\rho}|$ attained during the evolution $t\geq 0$. 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 $|\bm{\rho}(t)|$ as a function of $K$. As earlier, for $K$ below the corresponding $K_{c}$ this maximum value will be approximately zero, and for $K$ above $K_{c}$ the rapid macroscopic burst of $|\bm{\rho}|$ will be apparent with a larger maximum value of $|\bm{\rho}|$. Thus we expect the onset of such transient bursts for the three cases at the theoretically described values $K_{c}^{(-)}$, $K_{c}^{(u)}$ and $K_{c}^{(+)}$, respectively, according to Eq. (63). 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>K_{c}^{(+)}$ the maximum attained value corresponds to the stable state of $|\bm{\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 $K_{c}$, we examined the case of an abrupt increase in $K$ from a value of $K<K_{c}$ to a value $K>K_{c}$. The distribution $F$ remains invariant to evolution for $K<K_{c}$, 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 $F_{1}\in\mathcal{M}$.

While the system state is away from $\mathcal{M}$ and is being attracted towards $F_{1}$, if the value of $K$ is altered again to one greater than $K_{c}^{(1)}$, then the system will again be repelled away from $\mathcal{M}$. As the system is then attracted towards another distinct distribution $F_{2}$ (with a critical coupling strength of $K_{c}^{(2)}>K_{c}^{(1)}$), the coupling strength can be varied again to $K>K_{c}^{(2)}$, resulting in additional delay in the attraction towards $\mathcal{M}$. In this fashion, if we consider a slowly increasing coupling strength, then we can delay this attraction towards the manifold $\mathcal{M}$ for large amounts of time, resulting in $|\bm{\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 $|\bm{\rho}|$ as a function of $K$ and time. We observe $|\bm{\rho}|$ to show a small fluctuating increase at $K\approx K_{c}^{(-)}$, which is sustained until $K\approx K_{c}^{(+)}$, after which $|\bm{\rho}|$ approaches the steady state value of $|\bm{\rho}|>0$ shown in Fig. 5. If we consider successively slower rates of increase of $K$, the resulting plot of $|\bm{\rho}|$ as a function of $K$ displays smaller fluctuations from $|\bm{\rho}|=0$ sustained through $K_{c}^{(-)}<K<K_{c}^{(+)}$. In the limit of an infinitely slow rate of increase of $K$, we expect that $|\bm{\rho}|$ will remain at zero for all $K<K_{c}^{(+)}$, reproducing Fig. 1.

IV.3 Resetting due to finite size

So far we have restricted our analysis to the $N\to\infty$ limit, wherein several incoherent equilibria in the manifold $\mathcal{M}$ can simultaneously be stable to perturbations orthogonal to $\mathcal{M}$ and are neutrally stable to perturbations in $\mathcal{M}$. Hence, a small perturbation within $\mathcal{M}$ can move an incoherent equilibrium in $\mathcal{M}$ to another nearby incoherent equilibrium in $\mathcal{M}$, and many such small perturbations can cumulatively cause a large change from an initial incoherent state. As we have observed earlier, transient dynamics away from $\mathcal{M}$ appear to shift the critical coupling strength for loss of stability towards $K_{c}^{(+)}$. Thus, we suspect that perturbations away from $\mathcal{M}$ are biased towards maximally stable states. 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\to\infty$ as small, we can regard the system with $N$ finite but large as being akin to the $N\to\infty$ limit system with small added perturbations. Thus we might suspect finite, large $N$ to induce a slow secular evolution of the $N\to\infty$ incoherent equilibria towards a maximally stable state within $\mathcal{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 state. We demonstrate this effect in Fig. 6(a), where we plot $\langle C_{min}\rangle$ as a function of time for evolution of $N=10^{3}$ agents initialized to the minimally stable distribution with $\langle C_{min}\rangle=0$, evolved with $K=0.9<K_{c}^{(-)}$. Note that $\langle C_{min}\rangle$ undergoes slow growth and eventually asymptotes to a large value of $\langle C_{min}\rangle\approx 0.7$ at very long times. After this long time, the large value of $\langle C_{min}\rangle$ 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 $K_{c}^{(-)}<K<K_{c}^{(+)}$, we do not observe any characteristic burst in the value of $|\bm{\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 states restored as $N\to\infty$. This picture is confirmed numerically Fig. 6, where it can be clearly seen that for a larger value of $N=10^{4}$ initialized as earlier with $\langle C_{min}\rangle=0$ and evolved at $K=0.9<K_{c}^{(-)}$ takes about ten times longer time to reach an asymptotic state for $\langle C_{min}\rangle$ (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 state with $K<K_{c}$ (where $K_{c}$ is calculated in the $N\to\infty$ limit) and wait for sufficiently long time, then one could continuously increase $K$ without the incoherent state becoming unstable until $K$ reaches $K_{c}^{(+)}$.

V Conclusions

In this paper we look at a $D$-dimensional generalization of the Kuramoto modelChandra, Girvan, and Ott (2018). Unlike the case of the standard ($D=2$) Kuramoto model, we have shown that for even $D\geq 4$ there are an infinite number of time-independent distributions of agents (defining the manifold $\mathcal{M}$) that correspond to the completely incoherent state (i.e., having $|\bm{\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 $K_{c}$ corresponding to the given distribution (Sec. III). Further, for each value of $K_{c}$ within a range $K_{c}^{(-)}<K_{c}<K_{c}^{(+)}=2K_{c}^{(-)}$ there exists an infinite number of distributions that become unstable as $K$ is increased through $K_{c}$. In Sec. IV we show that these properties result in transitions within the $|\bm{\rho}|=0$ manifold $\mathcal{M}$ of steady states as $K$ is increased in the range $[K_{c}^{(-)},K_{c}^{(+)}]$, which leave their signatures as short-lived macroscopic bursts in the value of $|\bm{\rho}|$ (Fig. 3). These transitions imply a change in the microscopic state of the system, with the system state 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<K_{c}^{(+)}$ the only stable steady states are on $\mathcal{M}$, we demonstrate (Fig. 5) that considering a linearly increasing $K$ results in a small positive fluctuating value of $|\bm{\rho}(t)|$ (and hence indicating evolution not on $\mathcal{M}$) which can be sustained for long periods of time as $K$ is linearly increased through the range $K_{c}<K<K_{c}^{(+)}$ (where $K_{c}$ refers to the originally initialized distribution); also, these fluctuations in $|\bm{\rho}|$ become smaller as the rate of increase of $K$ with time becomes slower. Since there are a multitude of stable states on $\mathcal{M}$, with neutral stability to perturbations in $\mathcal{M}$, noise can cause slow evolution of states in $\mathcal{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.

VI Acknowledgement

We thank Michelle Girvan and Thomas M. Antonsen for useful discussion. This work was supported by ONR grant N000141512134 and by AFOSR grant FA9550-15-1-0171