A Coupled Oscillator Model for the Origin of Bimodality and Multimodality

TL;DR

This paper demonstrates that bimodal and multimodal distributions can naturally arise from repulsive coupling in oscillator models, challenging the assumption that natural distributions are typically unimodal.

Contribution

It introduces a rigorous analysis of how multimodality emerges from inhibitory coupling in variants of the Kuramoto model, expanding understanding of distribution shapes in coupled systems.

Findings

Multimodality can result from repulsive coupling dynamics.

The analysis applies broadly to various coupling functions.

Multimodal distributions emerge naturally in oscillator systems.

Abstract

Perhaps because of the elegance of the central limit theorem, it is often assumed that distributions in nature will approach singly-peaked, unimodal shapes reminiscent of the Gaussian normal distribution. However, many systems behave differently, with variables following apparently bimodal or multimodal distributions. Here we argue that multimodality may emerge naturally as a result of repulsive or inhibitory coupling dynamics, and we show rigorously how it emerges for a broad class of coupling functions in variants of the paradigmatic Kuramoto model.

Figures

Equations

$${\dot{\theta }}_i={\omega }_i+\frac{K}{N}\sum _{i=1}^{N}f({\theta }_j-{\theta }_i),i=1,\dots ,N.$$

$$f(0)$$

$$f^{\prime }(0)$$

$$f(u)$$

$$f^{\prime }(u)$$

$$f(\pi )$$

$${\dot{\theta }}_i={\omega }_0+\frac{K}{N}\sum _{j=1}^{N}f({\theta }_j-{\theta }_i),i=1,\dots ,N.$$

$${\dot{\theta }}_1$$

$$={\omega }_0+K(1-x)f({\theta }_2-{\theta }_1)$$

$${\dot{\theta }}_2$$

$$={\omega }_0-Kxf({\theta }_2-{\theta }_1).$$

$$\dot{\psi }=-Kf(\psi ).$$

$$\dot{ϵ}$$

$$=-Kxf(-{\psi }_0+ϵ)-K(1-x)f(ϵ).$$

$$\dot{ϵ}\approx -ϵK\left[x{f}^{\prime }({\psi }_0)+(1-x)f^{\prime }(0)\right].$$

$$x{f}^{\prime }({\psi }_0)+(1-x)f^{\prime }(0)<0.$$

$$(1-x)f^{\prime }({\psi }_0)+x{f}^{\prime }(0)<0.$$

$$f^{\prime }(0)<(1-x)[f^{\prime }(0)-f^{\prime }({\psi }_0)]<-f^{\prime }({\psi }_0).$$

$$\frac{f^{\prime }(0)}{f^{\prime }(0)-f^{\prime }({\psi }_0)}<x<\frac{-f^{\prime }({\psi }_0)}{f^{\prime }(0)-f^{\prime }({\psi }_0)}.$$

$$f(u;a)=\frac{1}{{\pi }^2{a}^2}u\left({\pi }^2-u^2\right)\left(a^2-u^2\right).$$

$$\frac{{\pi }^2}{3{\pi }^2-2a^2}<\frac{2{\pi }^2-2a^2}{3{\pi }^2-2a^2},$$

$$|a|<\pi /\sqrt{2}\equiv a_{\text{crit}}.$$

$$f^{\prime }(0)<$$

$$+(x+z)f^{\prime }({\psi }_1+{\psi }_2)]$$

$$f^{\prime }(0)<-2f^{\prime }(\frac{2\pi }{3}).$$

$$a<\frac{2}{3}\sqrt{\frac{14}{3}}\pi \equiv a_{\text{tricrit}}\approx 1.44\pi .$$

$$\dot{\psi }=Kxf(-\psi )-K(1-x)f(\psi ).$$

$$xf(-{\psi }^{\ast })=(1-x)f({\psi }^{\ast })$$

$${\dot{\theta }}_1$$

$${\dot{\theta }}_2$$

$${\dot{\theta }}_3$$

$${\dot{\psi }}_1$$

$${\dot{\psi }}_2$$

$$f({\psi }_2)$$

$$f({\psi }_2+{\psi }_1)$$

$$x{f}^{\prime }(0)+y{f}^{\prime }({\psi }_1)+z{f}^{\prime }({\psi }_1+{\psi }_2)$$

$$y{f}^{\prime }(0)+x{f}^{\prime }({\psi }_1)+z{f}^{\prime }({\psi }_2)$$

Full text

††preprint: AIP/123-QED††preprint: AIP/123-QED

In this paper we employ oscillators as a test system for understanding how bimodality—the splitting of oscillators into two rather than one cluster—may emerge as a result of coupling between interacting units. We present numerical and analytical results showing that repulsive coupling can lead to bimodality (or multimodality) for a wide range of detailed interaction dynamics.

I Introduction

Synchronization is a widespread phenomenon observed in biological Saigusa et al. (2008); Myung et al. (2015); Taylor et al. (2009), chemical Toiya et al. (2010); Vaidyanathan (2015); Li et al. (2003), physical Pantaleone (2002); Yoshida et al. (2000); Henk and Alejandro (2003); Ulrichs, Mann, and Parlitz (2009), and social settings Repp (2005); Kirschner and Tomasello (2009); De Paula (2009); Pluchino, Latora, and Rapisarda (2006). A paradigmatic mathematical model that can explain synchronization in many contexts is the Kuramoto model Kuramoto (2003, 1975); Strogatz (2000); Acebrón et al. (2005); Moreno and Pacheco (2004). Much work has been done on understanding the complex and surprising dynamics of the Kuramoto model and its variants, but the vast majority of that research focuses on the case of attractive coupling; here we are interested in the case where the coupling is repulsive.

Repulsive (or inhibitory) coupling is of physical interest as it arises frequently in the context of neuronal networks (e.g., see refs. Myung et al., 2015; Wang, Chen, and Perc, 2011), chemical interactions (e.g., refs. Toiya et al., 2010; Epstein, 1991; Hohmann, Kraus, and Schneider, 1998), and many other systems (see refs. Koseska et al., 2007; Ullner et al., 2007, 2008; Marvel, Strogatz, and Kleinberg, 2009; Laje and Mindlin, 2002; Balázsi et al., 2001). Some coupled oscillator models have examined repulsive coupling: Giver et al. developed a local variant of the Kuramoto model with repulsive coupling based on the interaction between water micro-droplets with reactants of the Belousov-Zhabotinsky reaction Giver, Jabeen, and Chakraborty (2011). Hong and Strogatz developed two variants of the Kuramoto model that involved mixes of positive and negative coupling Hong and Strogatz (2011, 2012).

The relationship between network structure and repulsive coupling has also been analyzed, with LevnajićLevnajić (2011, 2012) showing that, given the network coupling structure, many different phase configurations can arise. Recently, it has been shown that synchronization can arise in both repulsive and attractive coupling scenarios subject to common noise Pimenova et al. (2016); Nagai and Kori (2010); Gil, Kuramoto, and Mikhailov (2010); Gong et al. (2019). Gong et al.Gong et al. (2019), inspired by the work of Gil et al.Gil, Kuramoto, and Mikhailov (2010), studied instances where common noise can lead to clustering in the phase distribution of oscillators for repulsive coupling.

Nakamura et al.Nakamura, Tominaga, and Munakata (1994) investigated the effect of time-delayed nearest-neighbor coupling in the Kuramoto model and found that it could lead to the development of clustered states for both attractive and repulsive coupling. Mishra et al.Mishra et al. (2015) demonstrated that “chimeralike” states could arise with globally coupled Liénard systems incorporating both attractive and repulsive mean-field feedback. Yeldesbay et al.Yeldesbay, Pikovsky, and Rosenblum (2014) established that chimeralike states can arise in the Kuramoto-Sakaguchi model. They also considered a model with oscillators that could be synchronous (attractive coupling) or asynchronous (repulsive coupling) depending on their natural frequencies. They found that in this case a chimera state arises.

Golomb et al.Golomb et al. (1992) showed that clustering is possible in a coupled oscillator model with repulsive coupling that is suited for strong interactions between the limit-cycle oscillators. They further provided theory for when a frequency locked stationary phase distribution and when a nonperiodic attractor can arise.

Tsimring et al.Tsimring et al. (2005) showed that heterogeneous globally coupled oscillators obeying the standard Kuramoto model can cluster with all configurations having a zero order parameter, but this clustering breaks down as the number of oscillators increases. They also showed that, with local coupling, clustering can occur for nonidentical oscillators given sufficiently large coupling strength.

Closest to the work we present here, OkudaOkuda (1993) looked at the effect that an arbitrary coupling function may have on oscillators and developed theory as to when an $n$-cluster state, with all clusters being the same size, can arise. He found that harmonics in the coupling function are necessary for clusters to arise.

The central limit theorem Billingsley (1995) may influence us to expect that distributions in nature should tend to a singly-peaked, unimodal shape akin to the Gaussian normal distribution. Yet, bimodality and multimodality can be observed in biological Chaudhary, Saeedi, and Costello (2016); Seeman et al. (1984a); Lewis and Rinzel (2003), social Wu et al. (2010a); DiMaggio, Evans, and Bryson (1996); Flay (1978), and chemical Bonifacio et al. (2015); Reiter and Epstein (1987); Ren-Chang, Jian-Chun, and Qiu-Jun (2011); Sun and Nesbitt (1977) contexts and beyondDekel and Birnboim (2006); Zhang, Mapes, and Soden (2003); Kuhl and Meltzoff (1982) (see Fig. 1 for selected examples). In this paper we demonstrate that multimodality may arise as a result of repulsive or inhibitory coupling dynamics and we give an in-depth explanation of how it can arise for a range of coupling functions.

II Model with antisymmetric repulsive coupling

We begin by considering a system of $N$ phase oscillators characterized by natural frequencies $\omega_{i},~{}i=1\ldots N$. The oscillators are globally coupled with coupling strength $K$ through an interaction function $f$ that depends only on the phase difference between each pair of oscillators:

[TABLE]

Here $K>0$ represents attractive coupling and $K<0$ represents repulsive coupling.

We consider interaction functions $f(u)$, $u\in(-\pi,\pi]$, that satisfy the following conditions:

[TABLE]

These conditions impose: (2a) no coupling effects between oscillators in sync; (2b) locally attractive (repulsive) coupling near sync state for $K>0$ ($K<0$); (2c) odd interaction function; (2d) no discontinuities in $f^{\prime}(u)$; (2e) $2\pi$-periodic interaction function on $(-\pi,\pi]$ domain. We point out that conditions $(\ref{Oddfun})$ and $(\ref{PeriodicFun})$ lead to $f(\pi)=\lim_{u\to-\pi^{+}}f(u)=0$.

II.1 Identical Oscillators

We assume that oscillators frequencies are drawn from a known frequency distribution $g(\omega)$. For simplicity we first consider the case of identical oscillators, i.e., we set the distribution to be $g(\omega)=\delta(\omega-\omega_{0})$, so the system becomes

[TABLE]

II.2 Bimodal equilibria

We assume that the number of oscillators is large, $N\gg 1$, and we look for bimodal equilibria by making the ansatz of an oscillator phase distribution $h(\theta)=x\delta(\theta-\theta_{1})+(1-x)\delta(\theta-\theta_{2})$, where $0<x<1$ describes the fraction in cluster 1. Note that this constitutes an explicit restriction to a bimodal manifold within the broader space of all possible oscillator phase distributions. Then system (3) can be reduced to two coupled ordinary differential equations (ODEs):

[TABLE]

We define a new phase-difference variable $\psi=\theta_{2}-\theta_{1}$ and write its dynamical system by subtracting Eq. (4) from Eq. (5):

[TABLE]

We observe that the fixed points of the system for $\psi$ are fully determined by the zeros of $f(\psi)$. From the assumptions above $f(\psi)$ must have zeros at $\psi=0$ and $\psi=\pi$. Furthermore, if conditions (2a–2e) hold and $f(\psi)$ has no other zeros (as in the case of the red dashed curve from Fig. 2), then it is implied that $f^{\prime}(\pi)\leq 0$. Hence, within the bimodal manifold, the fixed point at $\psi=\pi$ should be stable with $\psi=0$ being unstable. $\psi=\pi$ corresponds to a bimodal equilibrium with two clusters of oscillators separated by 180∘ of phase.

If additional roots of $f(\psi)$ exist between [math] and $\pi$, these will also correspond to bimodal fixed points with alternating stability (again restricted to the bimodal manifold). We focus on the cases where there are no other fixed points or there is exactly one other fixed point $\psi_{0}$ in $(0,\pi)$; other cases are similarly tractable. Figure 2 illustrates the typical general shapes of the interaction functions that we consider.

II.2.1 Stability of bimodal equilibrium

To investigate the broader stability of solutions to perturbations outside the bimodal manifold, we consider the perturbation of a single oscillator by a small amount $\epsilon$. Because $N\gg 1$, we approximate the dynamics of the two clusters as unaffected by this perturbation. We examine the evolution of distance between the perturbed oscillator and the group from which it was perturbed, $\epsilon(t)$, to evaluate whether the system returns to its initial state.

For convenience, we move into a rotating frame by redefining $\theta_{i}\rightarrow\theta_{i}+\omega_{0}t$, which is equivalent to setting $\omega_{0}=0$. Without loss of generality we choose oscillator index $N$ from the $\theta_{2}$ cluster for the perturbation and assume $\theta_{1}=0$, and thus $\theta_{2}=\psi_{0}\leq\pi$ (assuming for now that our interaction function has only one or zero fixed points in $(0,\pi)$). Then $\theta_{N}=\theta_{2}-\epsilon=\psi_{0}-\epsilon$, and

[TABLE]

We expand the functions $f$ in a Taylor series to linear order:

[TABLE]

Assuming that $K<0$ (repulsive coupling, our case of interest in this manuscript), this implies stability if and only if

[TABLE]

A nearly identical calculation starting with the perturbation of a single oscillator from the $\theta_{1}$ (zero phase) cluster leads to a similar equation,

[TABLE]

Since Eqs. (7) and (8) must be simultaneously satisfied for stability of the full bimodal distribution, the following inequality must hold:

[TABLE]

Interestingly, this implies that the slope of the interaction function $f(\psi)$ must be steeper at $\psi=\psi_{0}$ compared to $\psi=0$ if the bimodal state is to be stable. We can also compute explicit bounds on the proportion of the oscillators in each group by isolating fraction $x$ in inequality $(\ref{PreBiModCond})$ :

[TABLE]

III Concrete example

As a concrete example, we consider a simple class of interaction functions

[TABLE]

These functions have roots on $(-\pi,\pi]$ at [math], $\pi$, and $\pm a$, and satisfy all the conditions set forth earlier in section II. As long as $0<|a|<\pi$ there are three roots in $0\leq u\leq\pi$, and one can check that $f^{\prime}(0)=1$ for all choices of $a$ (see Fig. 3 for example plots). For inequality (10) to be satisfiable, we require

[TABLE]

which reduces to

[TABLE]

We note that symmetry of the roots allows us to consider positive $a$ without loss of generality. Figure 4 shows the results of numerical experiments where we test this predicted stability threshold. In each panel, Eq. (3) is implemented with the interaction function from Eq. (11). We initialize $xN$ oscillators at $\theta_{1}=0$ and $(1-x)N$ at $\theta_{2}=a$, then add a small random perturbation $\xi_{i}$ to each oscillator’s initial phase, where $\xi_{i}$ is drawn from the normal distribution $\mathcal{N}(0,\delta^{2})$, with $\delta=0.1$ used in Fig. 4. We numerically integrate the system using a 4th/5th order Runge-Kutta scheme and consider evidence for stability if it approaches the unperturbed state, i.e. $\psi=\theta_{2}-\theta_{1}\rightarrow a$ with $x_{\text{final}}=x_{\text{initial}}$. We note that in these experiments we set coupling strength $K=-1000$.

In panels (a) and (b), we use $N=100$ oscillators, $\omega_{0}=0$, and set $a=\pi/2$, consistent with the stability threshold from Eq. (12), $a<a_{\textrm{crit}}=\pi/\sqrt{2}$. The stable band of fractionation according to inequality (10) is then $2/5<x<3/5$. In panel (a), we set $x_{\text{initial}}=0.55$, below the band’s upper bound; in panel (b), we set $x_{\text{initial}}=0.65$, above the band’s upper bound. As expected, the bimodal equilibrium appears stable in panel (a), but unstable in panel (b), where eleven oscillators move between clusters to establish a different equilibrium within the stable fractionation band ($2/5<x_{\text{final}}=0.54<3/5$).

In panels (c) and (d), we again use $N=100$ oscillators and $\omega_{0}=0$, but here we examine the predicted stability threshold $a_{\text{crit}}=\pi/\sqrt{2}$ from Eq. (12). We expect the bimodal state with $\psi^{*}=a$ to be unstable for all positive $a>a_{\text{crit}}$ (but note that this state ceases to exist when $a>\pi$). We set $x_{\text{initial}}=1/2$ since this is within the fractionation stability band from inequality (10) for all $a<a_{\text{crit}}$. In panel (c), we set $a=a_{\text{crit}}-0.1$, just below the threshold for stability; in panel (d), we set $a=a_{\text{crit}}+0.1$, just barely in the unstable domain. As expected, the bimodal equilibrium again appears stable in panel (c), but it appears unstable in panel (d). Since no fractionation $x$ will lead to a stable bimodal equilibrium, the system must move to an entirely different state, and it appears to converge to a trimodal distribution of oscillator phases.

We are able to understand why the system converges to a trimodal state by performing a similar analysis for the stability of three-cluster, or trimodal, oscillator distributions. One can show that a necessary condition for stability is:

[TABLE]

where $\psi_{i}$ is the angle separating clusters at $\theta_{i}$ and $\theta_{i+1}$ ($\theta_{4}$ identified with $\theta_{1}$), and $x$, $y$, and $z$ are the fractionations of the three clusters at $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ respectively. With equal spacing between the clusters $\psi_{1}=\psi_{2}=2\pi-\psi_{1}-\psi_{2}=2\pi/3$, the necessary condition simplifies to

[TABLE]

For the example function shown in Eq. (11) this is

[TABLE]

This implies that a trimodal state remains stable for all $a<\pi$. It stably coexists with the bimodal state for $a<\pi/\sqrt{2}$, and may coexist with other multimodal states for $\pi/\sqrt{2}<a<\pi$. In general different multimodal states may stably coexist over various parameter ranges. More details of the analysis for trimodality can be found in the appendix.

IV Generalization to asymmetric interaction functions

We can relax assumption (2c) of an antisymmetric coupling function and still find stability boundaries for multimodal states. In place of Eq. $\eqref{eq:psidot}$ (which used oddness of the coupling function), we find instead

[TABLE]

Clearly $\psi^{*}=0$ and $\psi^{*}=\pi$ both remain fixed points. Other fixed points exist if

[TABLE]

has a solution on $-\pi<\psi^{*}\leq\pi$. Figure 5 shows an example of an asymmetric interaction function. Geometrically this condition can be understood as identifying intersections of $f(\psi)$ and its reflection $f(-\psi)$ when $x=1/2$ (or scaled versions when $x\neq 1/2$). Once multimodal fixed points are identified, stability analysis is analogous to that presented earlier.

V Generalization to non-identical oscillators

We argue that real-world bimodal or multimodal distributions may result from similar dynamics to those presented in this paper. Of course, heterogeneity is inevitable in most real-world systems, yet we have focused thus far on the case of identical oscillators. While we leave the more general analysis for future work, we have conducted numerical experiments that appear to show that the predicted behavior occurs even in the presence of oscillator heterogeneity.

Again using the same example interaction function from Eq. (11), we now draw frequencies, $\omega_{i}$, from a normal distribution $\mathcal{N}(0,\sigma^{2})$ and set the initial phases of the oscillators to $\theta_{i}=\xi_{i}$ (fraction $x$) or $\theta_{i}=a+\xi_{i}$ (fraction $1-x$), where $\xi_{i}$ is a small perturbation draw from the distribution $\mathcal{N}(0,\delta^{2})$. Figure 6 shows the results of perturbation experiments analogous to those presented in Fig. 4, with analogous results except that the final phase distributions have phases that cluster about the modes rather than all converging to them precisely (right panels show histograms of final states).

In Fig. 6 panels (a) and (b), we use $N=1000$ oscillators and set $a=\pi/2<a_{\text{crit}}$ and $x_{\text{initial}}=1/2$. Even with perturbed initial phases and heterogeneous natural frequencies, the oscillators still remain in the bimodal state as predicted for $a<a_{\text{crit}}$. Specifically, panel (b) shows that the steady state distribution of oscillators has finite-width clustering about the fixed point positions predicted from the identical-oscillator case. In panels (c) and (d), since $a=\pi/\sqrt{2}+0.1>a_{\text{crit}}$, the bimodal state breaks down (consistent with the prediction of the identical-oscillator theory) and the system appears to converge to a trimodal equilibrium with three finite-width clusters.

VI Discussion

Coupled oscillators are an excellent testbed for models of synchronization or clustering. Even though real-world variables (e.g., sediment grain sizeSun et al. (2004), salmon body sizeThorpe (1977); Damsgrd et al. (2019), human communication frequencyWu et al. (2010b), dopamine receptor densitySeeman et al. (1984b), neutron star massSchwab, Podsiadlowski, and Rappaport (2010), galaxy colorBaldry et al. (2004), gamma ray burst durationMao, Narayan, and Piran (1994), tree heightEichhorn (2010), animal ornament sizeClifton, Braun, and Abrams (2016a)) may not be oscillatory or confined to a periodic domain, bimodality may emerge for qualitatively similar reasons. In our model, the coupling of one unit’s dynamical behavior to that of others is key to the phenomenon.

For clarity of presentation we have focused on a single example of interaction function (Eq. (11)), but evaluation of two other classes of interaction functions (triangle waves and antisymmetrized von-Mises kernels) also supports our analytical results—see supplementary material for details. In supplementary material, we also present further results regarding dependence of bimodal equilibria on coupling strength $K$, as well as some numerical evidence regarding sizes of basins of attraction; each of these topics merits further in-depth study. The analysis we present here focuses exclusively on the case of all-to-all coupling; we leave further investigation of the impact of network structure for future work.

For real-world scenarios where bimodality or multimodality is of interest, the interaction function may not be known exactly. Nevertheless, we expect that it will often be possible to assess whether the conditions expressed in Eqns. (2) hold in a particular case. It also seems plausible that functions describing real-world interactions between coupled systems will have no more than a handful of roots, making bimodality and trimodality likely outcomes when repulsive or inhibitory coupling is imposed.

One particularly important case occurs when the interaction function has only roots at zero and $\pi$, with the root at zero having larger or equal magnitude slope. That is the case in the standard Kuramoto model with sinusoidal coupling. In such a case we expect that the incoherent splay state will be stable. In general, the splay state should be stable when the tendency to cluster (due to long-distance interactions) cannot overcome the oscillators’ locally repulsive interactions.

VII Conclusions

We have shown that, when coupling is repulsive, multi-modality of the oscillator distribution can be a stable configuration for a wide range of interaction functions. We showed that bimodality can be expected under repulsive coupling when the slope of the interaction function at the origin is shallower than at the other root(s). We performed numerical experiments for both identical and nonidentical oscillators and observed results consistent with theory.

This demonstration that repulsive coupling can produce clustering under reasonable assumptions about the interaction dynamics is important as repulsive coupling is present in many natural systems. Hence, the theory we present in this paper provides an argument as to why one might expect multi-modality instead of unimodality or incoherence in systems known to have repulsive coupling.

Supplementary Material

See accompanying supplementary material for numerical experiments using a selection of additional interaction functions, for discussion of basins of attraction for different states, and for discussion of coupling strength dependency.

Appendix A Trimodal equilibria

We again consider a function $f(u)$ that satisfies conditions (2a)–(2e). We look for solutions with oscillators distributed according to $h(\theta)=x\delta(\theta-\theta_{1})+y\delta(\theta-\theta_{2})+(1-x-y)\delta(\theta-\theta_{3})$, where $x,y>0$, $x+y<1$ so that the oscillators will be in three clusters at $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ (we again assume that the natural frequencies are identical):

[TABLE]

Here, $z=1-x-y$. We define two variables $\psi_{1}=\theta_{2}-\theta_{1}$ and $\psi_{2}=\theta_{3}-\theta_{2}$, so that the system reduces to

[TABLE]

We set $\dot{\psi}_{i}=0$, $i=1,2$ and arrive at the following system of equations:

[TABLE]

To set bounds on the fractionation of the clusters, we assume that there exists points $\psi_{1},\psi_{2}\in(-\pi,\pi)$ such that Eqns. (21) and (22) are satisfied. Additionally, we put our system of coupled oscillators into a rotating frame so that $\theta_{i}\rightarrow\theta_{i}+\omega_{0}t$. In the rotating frame, we set $\theta_{1}=0$, $\theta_{2}=\psi_{1}$, and $\theta_{3}=\psi_{1}+\psi_{2}-2\pi$. As before we perturb an oscillator from one of the three groups. We do this for all three groups and get a system of inequalities

[TABLE]

All these must be simultaneously satisfied for stability of a trimodal state. Adding, we find

[TABLE]

This states that the weighted sum of the slopes of the coupling function at $\psi=\psi_{i}$, where the weights are the proportions for the groups separated by $\psi_{i}$, is greater in magnitude than the slope at the origin. This condition reduces to

[TABLE]

if $\psi_{1}=\psi_{2}=2\pi-\psi_{1}-\psi_{2}=2\pi/3$. As an example, we return to the class of interaction functions that we introduced in Section II. We relax the assumption that $|a|<\pi$ and consider the case when $\psi_{1}=\psi_{2}=2\pi-\psi_{1}-\psi_{2}$. To satisfy inequality $(\ref{StabTriMod})$, this means that

[TABLE]

which reduces to

[TABLE]

Figure 7 shows the results of a numerical experiment where we test this threshold. In both panels we use $N=99$, $x=y=z=1/3$, $\psi_{1}=\psi_{2}=2\pi-\psi_{1}-\psi_{2}=2\pi/3$, and set $\omega_{0}=0$. We expect the trimodal state to be unstable for $a>a_{\text{tricrit}}$. In panel (a) we set $a=a_{\text{tricrit}}+0.1$ and perturb the oscillators by amount $\xi_{i}$, with values drawn from the distribution $\mathcal{N}(0,0.01)$ . We can see that this perturbation leads to the system leaving the trimodal state and going to a bimodal state with $180^{\circ}$ phase difference.

One might be interested in why the bimodal state is stable in panel (a). Since there are only zeros at $\psi=0$ and $\psi=\pi$, one may check the stability by evaluating the derivative of $f(\psi)$ at these points. One can show that if

[TABLE]

the $180^{\circ}$ antiphase state is stable. Thus, when $a=a_{\text{tricrit}}+0.1>\sqrt{2}\pi$ the trimodal state becomes unstable and perturbations lead to the stable bimodal state.

For the case, when $\sqrt{2}\pi<a<a_{\text{tricrit}}$, both the trimodal state and the bimodal state are stable configurations. Figure 8 shows the result of the numerical experiment where we place the parameter inside the previously stated interval and outside of the interval. In all panels we use $N=300$, and set $\omega_{0}=0$. As before, in all panels we perturb the oscillators by amount $\xi_{i}$ from the predicted fixed points, whose values are drawn from the distribution $\mathcal{N}(0,0.01)$. In panels (a) and (b) we set $a=1.43\pi\in(\sqrt{2}\pi,a_{\text{tricrit}})$. In these cases we expect both the bimodal state and the trimodal state to be stable for this value of $a$. In panels (a) and (b), we set the fractionation to be equal in all groups, and we set the spacing between groups to be equal. As expected, we see that the trimodal state and the bimodal state are stable under perturbation.

In panel (c) we set $a=1.43\pi-0.1<\sqrt{2}\pi<a_{\text{crit}}$. As expected, we see that the bimodal state is unstable and the system goes in to trimodal state. Given the proximity of the clusters to $\pm\pi$, we have added black dashed lines that at $\pm\pi$, so that one can see that the difference between the final state and $\pm\pi$. In panel (d), we set $a=1.43\pi+0.1>a_{\text{tricrit}}>\sqrt{2}\pi$. We also observe an expected result, as trimodality appears to be unstable and the system converges to a bimodal equilibrium, which is stable given that $a>\sqrt{2}\pi$.

In summary, we have a necessary condition for the stability of the trimodal equilibrium. Although, this condition is only necessary for stability, not sufficient, numerical experiments seems to point to it being an accurate threshold in examples we have considered. Also, theory and numerical experiments demonstrate that multistability of different multimodal equilibria is possible over parameter space. The theory for the stability of higher modes we leave for future work.