Bifurcation analysis and structural stability of simplicial oscillator populations
aa r X i v : . [ n li n . AO ] M a y Bifurcation analysis and structural stability of simplicial oscillator populations
Can Xu, ∗ Xuebin Wang, and Per Sebastian Skardal † Institute of Systems Science and College of Information Science and Engineering, Huaqiao University, Xiamen 361021, China Department of Mathematics, Trinity College, Hartford, Connecticut 06106, USA
We present an analytical description for the collective dynamics of oscillator ensembles withhigher-order coupling encoded by simplicial structure, which serves as an illustrative and insightfulparadigm for brain function and information storage. The novel dynamics of the system, includ-ing abrupt desynchronization and multistability, are rigorously characterized and the critical pointsthat correspond to a continuum of first-order phase transitions are found to satisfy universal scalingproperties. More importantly, the underlying bifurcation mechanism giving rise to multiple clus-ters with arbitrary ensemble size is characterized using a rigorous spectral analysis of the stablecluster states. As a consequence of SO group symmetry, we show that the continuum of abruptdesynchronization transitions result from the instability of a collective mode under the nontrivialantisymmetric manifold in the high dimensional phase space. I. INTRODUCTION
Spontaneous synchronization in populations of inter-acting units is a ubiquitous phenomena in complex sys-tems [1, 2]. Describing such collective behaviors in termsof coupled phase oscillators and studying phase transi-tions between incoherence and synchrony have provenuseful in of physics, biology and social systems [3–5]. Ex-ploring the routes towards synchrony and uncovering theunderlying mechanism in various levels have attractedincreasing theoretical and experimental interests [6–8].Most existing literatures focus on pairwise interactionbetween oscillators that typically contain the first har-monic based on a phase reduction [9, 10]. However, inmany cases, one needs to go beyond such setups allow-ing for non-pairwise coupling that incorporates high or-der structures, i.e., simplexes [11–15], which have gainedmore attention in recent years. Recent advances demon-strate that the simplicial interactions, sometimes calledhypernetworks, play essential roles in many systems rang-ing from signal transmission in neural networks to struc-tural function correlation in brain dynamics [16, 17].As shown in [18–20], notable features of higher-orderinteractions in coupled oscillator ensembles include theemergence of extensive multistability and a continuumof abrupt desynchronization transitions (ADTs), whichlead to potential applications in memory and informa-tion storage. Despite these findings, a rigorous analysisof the bifurcations and stability properties that charac-terize both the macroscopic and microscopic dynamics islacking, thereby leaving our understanding of these phys-ical phenomena incomplete.In this paper, we provide an analysis of synchronizationtransitions induced by simplicial interactions and nonlin-ear higher order coupling in oscillator ensembles. Usinga self-consistent approach we perform a bifurcation anal-ysis that uncovers a general phenomenon that gives rise ∗ [email protected] † [email protected] to ADTs and extensive multistability. Furthermore, weestablish scaling relations that describe the critical pointsfor desynchronization that are universal in the sense thatthey do not depend on the functional form of the naturalfrequency distribution. Most importantly, we performa rigorous stability analysis for finitely- and infinitely-many partially synchronized states with arbitrary sys-tem size using a spectral analysis of the stable clusterstates. We reveal that a continuum of ADTs originatefrom the instability of a collective mode under the non-trivial antisymmetrical manifold in the high dimensionalphase space.The remainder of this paper is organized as follows. InSec. II we present a bifurcation analysis for the globally-coupled system. In Sec. III we analyze the stabilityof cluster states and characterize the properties of theADTs. In Sec. IV we briefly present an additional ex-ample of a system with frequency-weighted coupling. Fi-nally, in Sec. V we conclude with a discussion of ourresults. II. BIFURCATION ANALYSIS
We consider oscillator ensembles with three-way non-linear coupling whose evolution is given by˙ θ i = ω i + 1 N N X m =1 N X n =1 K imn sin( θ m + θ n − θ i ) , (1)where θ i is the phase of oscillator i with i = 1 , . . . , N , N is the system size, and ω i is the natural frequencyof oscillator i , assumed to be drawn from a distribution g ( ω ), which is assumed to be even throughout this paper.Unlike the classical Kuramoto-like models, the couplingterm in Eq. (1) is not pairwise, but involves triplets. Fi-nally, K imn is the coupling strength among triplet of os-cillators i , m , and n . While we will consider more generalcases below, we first restrict our attention to the simplestsetup K mni = K ≥ Z k = R k e i Θ k = N P ∞ m =1 e ikθ m for k = 1 and 2. In the case where a singlesynchronized cluster emerges the classical Kuramoto or-der parameter Z is sufficient to describe the macroscopicsynchronization dynamics. However, when two clustersemerge the Daido order parameter Z measures the de-gree of (cluster) synchronization [22], while Z measurethe degree of asymmetry between the two clusters. R k and Θ k correspond to amplitudes and arguments, respec-tively, of the order parameters. By exponentiating thesine term and distributing the summation, these defini-tions allow us to rewrite Eq. (1) as˙ θ i = ω i + KR sin(2Θ − θ i ) , (2)where KR acts as an effective force on each oscillatoraiming to entrain it via the mean field Θ . The pres-ence of a higher order harmonic ( ∼ sin 2 θ i ) and nonlin-ear coupling ( ∼ KR ), gives rise to a series of nontrivialdynamical features.Proceeding with our analysis, it is convenient to in-troduce the parameter q = KR . Given the SO groupsymmetry (i.e., rotation and reflection) [23] of Eq. (2),we enter the appropriate rotating frame, i.e., consider thetransformation θ θ + Ω t (where Ω is the mean naturalfrequency), to obtain a non-rotating solution ( ˙Θ k = 0)and further set Θ k = 0 by appropriately shifting initialconditions. In this case, the whole population can be di-vided into those phase locked oscillators ( | ω i | < q ) andthat of drifting ones ( | ω i | > q ). For the phase locked case,the oscillators are entrained by the mean field, yieldingfixed points ˙ θ i = 0 satisfyingsin 2 θ i = ω i q and cos 2 θ i = s − ω i q . (3)In fact the second-order harmonic in yields two stablefixed points (as well as two unstable fixed points), whichmanifest in the formation of two clusters with phase dif-ference π [24, 25]. Trigonometric identities allow us towritesin θ i = sin 2 θ i θ i and cos θ i = ± r θ i , (4)where the + or − in the cosine term determine whether θ i falls in the cluster near θ = 0 or π , respectively. Herewe replace the ± in the cosine term with the parameter p i , which takes the value 1 with probability η ( ω i ) and is − − η ( ω i ). Thus, η ( ω i ) is an indi-cator function satisfying 0 ≤ η ( ω i ) ≤ ω i thatbecome entrained to the θ = 0 cluster.Turning briefly to the drifting oscillator population,each oscillator rotates non-uniformly on the unit circle with period T i = 2 π/ p ω i − q . Moreover, the contri-bution of drifiting oscillators to the order parameter Z k is given by N P | ω i | >q T i R π e ikθi | ˙ θ i | dθ i . In fact, the sym-metry of the natural frequency distribution causes thiscontribution to vanish, implying that drifting oscillatorsdo not contribute to either order parameter Z k whether N is finite or not. The contribution of the locked oscil-lators to the order parameters can be simplified by firstnoting that Z k = N P ∞ m =1 e ikθ m = N P ∞ m =1 cos( kθ m ) + i sin( kθ m ). Similar to the drifting oscillators, the evenand odd symmetries of the natural frequency distributionand sine term, respectively, cause the imaginary part tovanish, yielding R k = N P ∞ m =1 cos( kθ m ). Writing theseout explicitly, we have that R = 1 N ∞ X i =1 cos( θ i ) = 1 N ∞ X i =1 p i p (1 + cos 2 θ i ) / , (5) R = 1 N ∞ X i =1 cos(2 θ i ) = 1 N ∞ X i =1 p − ( ω i /q ) , (6)from which we see that R is determined by the param-eter q = KR that controls the range of phase locking,whereas R is further restricted by the indicator func-tion η that describes the distribution of p i ’s, reflectingthe asymmetry of phase locked oscillators between twoclusters, and therefore is of major importance. Convert-ing sum in Eq. (5) to a sum over frequencies (i.e., treatingthe population as a distribution) and replacing the cosineterm, we find that the parameterized equation for deter-mining R is given by1 √ K = 1 N X | ω i | q c ), respectively. η = 1 (red), 0 .
95 (green), 0 . .
85 (cyan), 0 . .
75 (purple). Dashed line denotes the theoretical dependencebetween R c and K c . For each value of η , the initial phases are set to 0 and π with probability η and 1 − η , respectively,then K decreases to 0 and restores to initial value with ∆ K = 0 .
01 adiabatically. In the simulation, N = 10 and g ( ω ) = , ω ∈ ( − , and the corresponding critical order parameters are R c ( η ) = (2 η − √ q c f ( q c ) , (10) R c ( η ) = 1 N X | ω i |
0) state and the incoherent (i.e., R k = 0). More-over, the incoherent state remains stable for all K in the N → ∞ limit. These two criteria imply that the systemundergoes, for a given value of η, an ADT characteriz-ing a discontinuous jump from a synchronized state tothe incoherent state, but no complementary transitionwhere the incoherent state spontaneously synchronizes.Moreover, a different ADT will occur at each differentvalue of η , so in fact a continuum of ADTs exist. Thisphenomenon differs essentially from the traditional first-order phase transition characterized by a finite hysteresisloop [26, 27] where (i) typically a single desynchroniza-tion transition occurs (often via a saddle-node bifurca-tion) and (ii) a complementary synchronization transi-tion typically follows at a larger value coupling strength.The self-consistent analysis and scaling formulas pre-sented above capture the macroscopic dynamics exhib-ited by Eq. (1) with arbitrary g ( ω ) and N . Note thatboth K c and R c exhibit a monotonic dependence on η (decreasing and increasing, respectively), whereas R c remarkably remains constant. For instance, for thecase of a standard normal frequency distribution, i.e., g ( ω ) = √ π e − ω [19], the ADTs occur at q c = 1 . f ( q c ) = 0 . g ( ω ) = π ( ω +1) [20], we havecritical points q c = 1 . f ( q c ) = 0 . f ( q ) is smooth. However, forthe non-smooth and finitely-supported uniform distribu-tion g ( ω ) = for ω ∈ ( − , q c = 1 since f ( q ) = √ q if q < , √ q q +( q − (1+ q − q √ q − q if q ≥ . (12)Moreover, f ′ ( q c ) does not exist since f ′ ( q + c ) = f ′ ( q − c ),then the corresponding critical points are K c = η − with R c = 2 √
23 (2 η −
1) = 1 √ K c , and R c = π . (13)In Fig. (1) we illustrate our results for the uniform dis-tribution, plotting the order parameters R and R inpanels (a) and (b), respectively, using η = 1 (red), 0 . . .
85 (cyan), 0 . .
75 (pur-ple). Analytical predictions are plotted as solid curves,while results from simulations with N = 10 oscillatorsare plotted with circles, both as the coupling strengthis decreased from a synchronized cluster state and in-creased from the incoherent state. the analytical pre-diction matches simulations. Critical values R c and R c are plotted as dotted curves, indicating the continuumof ADTs. Analytical predictions match the simulationresults very well. III. STABILITY
To better understand the system dynamics, we con-sider the structural stability question of what clusterstates emerge from an arbitrary configuration. A lin-ear stability analysis demonstrates that the incoherentstate can never lose its stability for any K . This followsfrom the fact that the continuous spectrum of the inco-herent state, i.e., R k = 0, is purely imaginary ( iω ) andthe nonlinear mean field has no contribution to discretespectrum [28]. Thus due to the purely nonlinear inter-action term in Eq. (2) it is impossible for a spontaneousphase transition from the incoherent state. We thereforeaim to investigate the stability of partially synchronizedstates, first focusing on finite N case, which then turnsout to generalize to the thermodynamic limit N → ∞ directly.The SO group symmetry ensures that the driftingoscillators generate a purely imaginary spectrum in thelinear stability analysis, thereby having no contributionsto the mean field [29, 30]. In other words, they effectivelydecouple from the locked oscillators in Eq. (2), so we canonly consider N l oscillators locked to the mean field forsome q excluding the drifting ones. We note, however,that in some cases, such as uniformly-distributed frequen-cies, all oscillators become locked in a synchronized state,yielding N l = N [31]. To study stability of such partiallysynchronized states, we consider a small perturbation x i , | x i | ≪
1, on each phase locked oscillator θ i in Eq. (2)and, after neglecting high order terms, the evolution forperturbation satisfies ˙ x = Jx where x = ( x , . . . , x N l )and J is the Jacobian matrix with entries J ij = ∂ ˙ θ i ∂θ j . Infact, J is of the form J = 2 KR N M − q D , (14)where the matrix M is given by M ij = cos( θ j − θ i )and a diagonal matrix D is given by D ij = cos(2 θ i ) δ ij .The stability properties for a given cluster state is thendetermined by the eigenvalues of J . First, note that P N l j =1 J ij = 0, which gives a trivial eigenvalue λ = 0(with constant eigenvector v ∝ = (1 , . . . , N l − J , which can be expressedas B ( λ ) = det( λ I − J ). Next, by defining E = λ I KR + R D , (15)one may write 2 KR E ( I − E − M /N ) . (16)Thus, the characteristic polynomial takes the form B ( λ ) = 2 KR det( E )det( I − N E − M ) . (17)Finally, assuming that E must is in fact invertible, wehave that det ( E ) = 0, and the key task for calculatingthe eigenvalue spectrum depends on finding the roots ofthe last term alone.Next we ease notation by defining four vectors c ( k ) with c ( k ) i = cos kθ i and s ( k ) with s ( k ) i = sin kθ i that satisfy theorthogonality property c ( k ) · s ( k ′ ) = 0 ( k ′ = 1 , k = 2 corresponds to purely deterministicvectors, while k = 1 is random vectors due to clusteringand multi-branches. That is, the entries of c (2) and s (2) do not depend on whether a given oscillator falls in the0 or π cluster, but the entries of c (1) and s (1) do. More-over, we note that the rank of M is only 2 since, for all x ∈ R N l , we have Mx = ( c (1) · x ) c (2) +( s (1) · x ) s (2) . Intro-ducing the orthonormal basis a = c (1) k c (1) k and a = s (1) k s (1) k ,the matrix E − M /N restricted to the subspace that isspanned by { a i } turns out to be a 2 × Q , i.e.,( E − M /N ) a , a = Q ( λ ). The entries of this matrix aredefined as Q ij = a i · N E − Ma j , whose diagonal entriesare Q = 1 N N l X i =1 [2 η ( ω i ) −
1] 2 KR | c (1) i | c (2) i λ + 2 qc (2) i , (18) Q = 1 N N l X i =1 [2 η ( ω i ) − KR ( s (2) i ) / | c (1) i | λ + 2 qc (2) i , (19)and the off-diagonal Q ij ∝ N P N l m =1 c ( i ) m s ( j ) m E mm with i = j ,while the other N l − E − M /N . Orthogonality implies that Q and Q arezero, so the characteristic polynomial simplifies to B ( λ ) = 2 KR N l Y i =1 E ii [1 − Q ( λ )][1 − Q ( λ )] . (20)In fact, the rational functions Q jj ( λ ) ( j = 1 ,
2) arestrictly decreasing away from their N l poles λ i = − qc (2) i (lim λ → λ ± i Q jj ( λ ) = ±∞ ), so Q jj ( λ ) = 1 must have oneroute between two consecutive poles. The necessary con-dition for stable locked state requires all c (2) i >
0. Inaddition, we find that Q (0) = 1 corresponds to a triv-ial eigenvalue λ = 0. The remaining one eigenvalue isuniquely determined by Q ( λ ) = 1 for λ > λ i . Since Q ( λ ) is decreasing, the root is negative if and only if Q (0) <
1, which leads to the structural stability con-dition for the configuration of population1 N N l X i =1 [2 η ( ω i ) − KR ( s (2) i ) q | c (1) i | c (2) i < , (21)or equivalently, F ′ ( q ) <
0. Therefore, we conclude thatthe eigen-spectrum of J is made up of three parts, a heav-ily populated part consisting of N l − λ = 0 correspond-ing to rotation invariance, and a lone eigenvalue outsidethe poles. In the thermodynamic limit N → ∞ , thefirst part merges into continuous spectrum with λ ( ω ) = − p q − ω , while the last part turns out to be a non-trivial discrete spectrum determined by the continuumlimit Q ( λ ) = 1 ( λ = λ ( ω )).It is also instructive to characterize the eigenvectors of J . Considering a frequency-dependent vector x ∈ R N l with x i = χ ( ω i ), the space R N l can be split into two sub-spaces V even and V odd . If χ is an even (or odd) function,then x ∈ V even (or V odd ). This definition allows x tobe random vector, such as c ( k ) ∈ V even , s ( k ) ∈ V odd , so V even ⊥ V odd and V = V even ⊕ V odd . The eigenspace of R c1 K decreasing K decreasing K increasing K (a) R K increasing K (b) R (c) c c R K K c +- - c q q (d) c q q FIG. 2: Bifurcation diagram of the order parameters R (a), R (b) with K for the frequency weighted coupling, g ( ω ) = π ( ω +1) .Circles are numerical simulations and solid lines are theoretical results (stable branches, q > q c ). η = 1 (red), 0 .
95 (green),0 . .
85 (cyan), 0 . .
75 (purple). (c) Analytical relation between R c and K c (black solid line) with g ( ω ) = √ π e − ω , blue circles represent numerical simulations with different η , which is in agreement with the dashed line (a)representing universal phase transition. (d) Sketch map of the stable configuration, the blue arrows (red) stand for a/an stable(unstable) perturbed direction. Numerical strategy is the same as Fig. (1). J can be described via the basis of V even and V odd . For x e ∈ V even , we have Jx e = KR N c (2) ( c (1) · x e ) − q Dx e .Setting c (1) · x e = 1, the eigenvector equation Jx e = λ x e implies that x ei = c (2) i / [ N E ii ( λ )] in which λ is a rootof Q ( λ ) = 1. Similarly, if x o ∈ V odd , imposing s (1) · x o = 1, we obtain x oi = s (2) i / [ N E ii ( λ )] corre-sponding to Q ( λ ) = 1. It is worth mentioning thatRe[ ∇ ( Z k ) · x e ] = 0 and Im[ ∇ ( Z k ) · x o ] = 0.The structural stability of the clusters can then be in-terpreted as follows. Recall that, after entering a rotatingframe and shifting initial conditions we have the meanphases Θ k = 0. This corresponds to two clusters: onecentered at θ = 0 and the other at θ = π . Preserving thisparticular symmetry, perturbation may then contract orspread oscillators in each cluster, effectively pulling os-cillators along the unit circle in the real or imaginarydirections, respectively. Our analysis above reveals thatthe eigenvalues determined by Q ( λ ) = 1 correspondto even eigenvectors, i.e., symmetric perturbation, thatinduce purely imaginary displacement of the centroid ofconfiguration in the linear approximation, thereby induc-ing a small spread in each cluster. On the other hand,the eigenvalues for Q ( λ ) = 1 correspond to odd eigen-vectors, i.e., antisymmetric perturbation, that leads topurely real displacement of the centroid, thereby con-tracting each cluster. As the structure parameter q isdecreasing, the nontrivial solitary eigenvalue is shiftedtowards positive real axis and collides with the origin at q = q c . Thus, the continuum of ADTs originate from theinstability of a collective mode under the nontrivial an-tisymmetric perturbation in the high dimensional phase space associated with this solitary eigenvalue. IV. ANOTHER EXAMPLE:FREQUENCY-WEIGHTED COUPLING
Lastly, we show that our theory generalizes by con-sidering the coupling heterogeneity of the form K mni = K | ω i | , i.e., establishing the correlation between oscilla-tors frequency and coupling strength as in Refs. [32,33] that generates explosive synchronization and quan-tized time-dependent clustering. We show that such afrequency-weighted coupling scheme captures the essen-tial properties of simplicial dynamics and can achieve cer-tain cluster arrangement. In this case, the SO groupsymmetry still holds, and the mean-field equation can bewritten as ˙ θ i = ω i − q | ω i | sin(2 θ i ) . (22)If q <
1, no oscillators are entrained by the mean-fieldindicating the asynchronous state. Interestingly, once q ≥
1, all oscillators become phase-locked ( N l = N ) with c (2) fi = p − q − and s (2) fi = ω i / ( q | ω i | ). The paramet-ric function F ( q ) is a smooth function for a constant η ,i.e., F ( q ) = (2 η − q [1 + (1 − q − ) ] / (2 q ) , (23)and the fold bifurcation characterizing a continuum ofuniversal ADTs occurs at q c = 2 / √ K c = 83 √ η − , R c = √ η − , and R c = 12 . (24)The results are presented in Fig. 2 (a)–(c) (details simi-lar to those for Fig. 1), displaying remarkable agreementbetween theory and simulation.The stability analysis for the multi-clusters follows thesimilar way, where the Jacobian matrix is J f = W ( 2 KR N M − q D ) (25)with the frequency matrix W being W ij = | ω i | δ ij . Theresulting characteristic polynomial is B f ( λ ) = 2 KR det( E f )[1 − Q f ( λ )][1 − Q f ( λ )] , (26)where E f = W − λ I KR + R D , (27) Q f = 1 N N X i =1 c (1) fi c (2) fi E fii ( λ ) , and (28) Q f = 1 N N X i =1 s (1) fi s (2) fi E fii ( λ ) . (29)The eigenvalue spectrum together with the associatedeigenvectors has the same appearance as uniform cou-pling. In particular, the phase locked states exhibit fixedcluster arrangement regardless of g ( ω ), namely, the os-cillators are recruited into two symmetric groups ± θ ∗ ( η = 1) or four groups ± θ ∗ , π ± θ ∗ ( η < Q f (0) < θ ∗ · tan 2 θ ∗ < θ ∗ < π . Decreasing q makes the symmetric clusters move away from each other(Fig. (2)d) and the configuration becomes unstable once q < q c . The ADTs occur owing to a saddle node bifurca-tion where all oscillators become unlocked correspondingto the resting state. V. DISCUSSION
Summarizing, we have presented rigorous analyticaldescriptions for the collective dynamics in a populationof globally coupled phase oscillators with higher-orderinteractions via simplicial structure. Extensive multista-bility of clusters and ADTs emerge directly from thesehigher order interactions and the nonlinear mean field.These results were obtained using a self-consistency anal-ysis and rigorous bifurcation theory and reveal scalingproperties of the transitions in the form of dependenceof the critical points on the indicator and structural con-stants. The rigorous characterization of the eigenvaluespectrum presented above demonstrates that the stabilityproperties of infinite many partially synchronized statesare controlled by a nontrivial solitary eigenvalue, whilethe occurrences of ADTs correspond to the instability ofa collective mode.
ACKNOWLEDGEMENTS
C. X. and X. W. acknowledge support from the Na-tional Natural Science Foundation of China (Grants Nos.11905068 and 11847013). [1] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchroniza-tion: a Universal Concept in Nonlinear Sciences (Cam-bridge University Press, Cambridge, England 2001)[2] S.H. Strogatz,
Sync: The Emerging Science of Sponta-neous Order (Hypernion, New York, 2003).[3] M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Self-organized synchronization in decentralized power grids,Phys. Rev. Lett. , 064101 (2012).[4] E. Montbri´o, and D. Paz´o, Kuramoto Model forExcitation-Inhibition-Based Oscillations, Phys. Rev.Lett. , 244101 (2018).[5] B. Ottino-L¨offler, and S. H. Strogatz, Volcano Transitionin a Solvable Model of Frustrated Oscillators, Phys. Rev.Lett. , 264102 (2018).[6] I.Z. Kiss, Y. Zhai, and J.L. Hudson, Emerging Coherencein a Population of Chemical Oscillators, Science ,1676 (2002).[7] J.A. Acebr´on, L.L. Bonilla, C.J. P´erez Vicente, F. Ritort,and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena. Rev. Mod. Phys. , 137(2005).[8] P.S. Skardal and A. Arenas, Control of coupled oscillatornetworks with application to microgrid technologies, Sci.Adv. 1, e1500339 (2015).[9] Y. Kuramoto, in International Symposium on Mathemat-ical Problems in Theoretical Physics , edited by H. Araki,Lecture Notes in Physics No. 30 (Springer, New York,1975), p. 420.[10] F.A. Rodrigues, T.K.D. M. Peron, P. Ji, and J. Kurths,The Kuramoto model in complex networks, Phys. Rep. , 1 (2016).[11] V. Salnikov, D. Cassese, and R. Lambiotte, Simpli-cial complexes and complex systems, Eur. J. Phys. ,014001 (2019).[12] C.C. Gong, and A. Pikovsky, Low-dimensional dynamicsfor higher order harmonic globally coupled phase oscilla-tor ensemble, arXiv preprint arXiv:1909.07718 (2019)[13] C. Bick, P. Ashwin, and A. Rodrigues, Chaos in generi- cally coupled phase oscillator networks with nonpairwiseinteractions, Chaos , 094814 (2016).[14] P. Ashwin and Ana Rodrigues, Hopf normal form with S N symmetry and reduction to systems of nonlinearlycoupled phase oscillators, Physica D , 14 (2016).[15] I. Le´on and D. Paz´o, Phase reduction beyond the firstorder: The case of the mean-field complex Ginzburg-Landau equation, Phys. Rev. E , 012211 (2019).[16] M.W. Reimann, M. Nolte, M. Scolamiero, K. Turner, R.Perin, G. Chindemi, P. Dlotko, R. Levi, K. Hess, andH.Markram,Cliques of neurons bound in to cavities pro-vide a missing link between structure and function, Front.Comput. Neurosci. , 48 (2017).[17] J. C. W. Billings, M. Hu, G. Lerda, et al, Simplex2Vecembeddings for community detection in simplicial com-plexes. arXiv preprint arXiv:1906.09068 (2019).[18] T. Tanaka and T. Aoyagi, Multistable Attractors in aNetwork of Phase Oscillators with Three-Body Interac-tions, Phys. Rev. Lett. , 224101 (2011).[19] M. Komarov, and A. Pikovsky, Finite-size-induced tran-sitions to synchrony in oscillator ensembles with nonlin-ear global coupling, Phys. Rev. E , 020901(R) (2015)[20] P.S. Skardal, and A. Arenas, Abrupt Desynchronizationand Extensive Multistability in Globally Coupled Oscil-lator Simplexes, Phys. Rev. Lett. , 248301 (2019)[21] J.-W. Wang, L.-L. Rong, Q.-H. Deng, and J.-Y. Zhang,Evolving hypernetwork model, Eur. Phys. J. B , 493(2010).[22] H. Daido, Multibranch Entrainment and Scaling in LargePopulations of Coupled Oscillators, Phys. Rev. Lett. ,1406 (1996).[23] J.D. Crawford, Amplitude expansions for instabilitiesin populations of globally-coupled oscillators, J. Statist. Phys. , 1047 (1994).[24] P.S. Skardal, E. Ott, and J.G. Restrepo, Cluster syn-chrony in systems of coupled phase oscillators withhigher-order coupling, Phys. Rev. E , 036208 (2011).[25] M. Komarov and A. Pikovsky, Multiplicity of SingularSynchronous States in the Kuramoto Model of CoupledOscillators, Phys. Rev. Lett. , 204101 (2013).[26] W.S. Lee, E. Ott, and T.M. Antonsen, Large coupledoscillator systems with heterogeneous interaction delays,Phys. Rev. Lett. , 044101 (2009).[27] J. G˜omez-Garde˜nes, S. G´omez, A. Arenas, and Y.Moreno, Explosive synchronization transitions in scale-free networks, Phys. Rev. Lett. , 128701 (2011).[28] S.H. Strogatz and R.E. Mirollo, Stability of incoherencein a population of coupled oscillators, J. Stat. Phys. ,613 (1991).[29] O.E. Omel’chenko,and M. Wolfrum, Nonuniversal tran-sitions to synchrony in the SakaguchiCKuramoto model,Phys. Rev. Lett. , 164101 (2012).[30] O.E. Omel’chenko,and M. Wolfrum, Bifurcations in theSakaguchiCKuramoto model, Physica D , 74-85(2013).[31] R.E. Mirollo, and S.H. Strogatz, The spectrum of thelocked state for the Kuramoto model of coupled oscilla-tors, Physica D , 249-266 (2005).[32] X. Zhang, X. Hu, J. Kurths, and Z. Liu, Explosive syn-chronization in a general complex network, Phys. Rev.E, 88, 010802 (2013).[33] H. Bi, X. Hu, S. Boccaletti, X. Wang, Y. Zou, Z. Liu,and S. Guan, Coexistence of quantized, time dependent,clustersing lobally coupled oscillators, Phys. Rev. Lett.117