# Multistability in coupled oscillator systems with higher-order interactions
and community structure

Per Sebastian Skardal<EMAIL_ADDRESS>Department of
Mathematics, Trinity College, Hartford, CT 06106, USA  Sabina Adhikari
Department of Applied Mathematics, University of Colorado at Boulder, Boulder,
Colorado 80309, USA  Juan G. Restrepo  Department of Applied Mathematics,
University of Colorado at Boulder, Boulder, Colorado 80309, USA

###### Abstract

We study synchronization dynamics in populations of coupled phase oscillators
with higher-order interactions and community structure. We find that the
combination of these two properties gives rise to a number of states
unsupported by either higher-order interactions or community structure alone,
including synchronized states with communities organized into clusters in-
phase, anti-phase, and a novel skew-phase, as well as an incoherent-
synchronized state. Moreover, the system displays a strong multistability,
with many of these states stable at the same time. We demonstrate our findings
by deriving the low dimensional dynamics of the system and examining the
system’s bifurcations using a stability analysis and perturbation theory.

Complex Networks, Synchronization

###### pacs:

05.45.Xt, 89.75.Hc

> Spontaneous entrainment and pattern formation in large ensembles of coupled
> oscillator units play a role in a wide range of applications in mathematics,
> physics, engineering, and biology Strogatz2003 ; Pikovsky2003 . Recent
> insights from both neuroscience and physics Petri2014Interface ;
> Ashwin2016PhysD  point to the presence of higher-order interactions in
> populations of coupled dynamical systems. While initial studies have begun
> to uncover the effects that higher-order interactions have on macroscopic
> system dynamics, the role that both small- and large-scale structural
> properties of networks play in shaping collective dynamics in the presence
> of higher-order interactions remains unexplored and poorly understood. Here
> we study the dynamics of systems of coupled phase oscillators with community
> structure and higher-order interactions. Using dimensionality reduction
> methods to derive the macroscopic system dynamics, we uncover and
> analytically characterize a range of new, coexisting states that are not
> supported by higher-order interactions or community structure alone.

##  I Introduction

Synchronization of coupled oscillators is one of the most prominent examples
of emergent collective behavior, with examples including circadian rhythms
Lu2016Chaos , brain rhythms KitzbichlerPLOS2009 , cardiac excitation Glass1988
, pedestrian synchronization in suspension bridges Strogatz2005Nature ,
Josephson junctions Wiesenfeld1996PRL , and power grids Rohden2012PRL . The
Kuramoto model Kuramoto1984  of coupled phase oscillators has been
particularly useful for shedding light on synchronization phenomena due to its
analytic tractability and versatility for modeling many different systems with
various physical properties. While the Kuramoto model describes the pairwise
interactions of oscillators, it does not take into account the possibility of
interactions involving multiple nodes simultaneously, i.e., higher-order
interactions. In recent years, the role of higher-order interactions in
shaping the collective response of coupled dynamical systems has received much
interest Battiston2020PhysRep ; Majhi2022Interface . Higher-order interactions
have been shown to qualitatively change the dynamics of epidemic and opinion
dynamics Iacopini2019NatComms ; Matamalas2019 , the behavior of ecological
systems GrilliNature2017 , and other systems that display emergent collective
dynamics Carletti2020JPhys . For coupled oscillator systems, higher-order
interactions are motivated by applications in neuroscience Petri2014Interface
; Giusti2016JCN ; Reimann2017 ; Sizemore2018JCN  and physics Ashwin2016PhysD ;
Leon2019PRE . Higher-order interactions can cause dramatic changes to the
dynamics of coupled oscillators, including bistability, explosive transitions,
and extensive multistability Skardal2019PRL ; Skardal2020CommPhys ; Xu2020PRR
; Skardal2020JPhys ; Xu2021PRR ; Skardal2021PRR ; Anwar2022Chaos ;
Skardal2022Chaos ; Leon2022PRE ; Kachhvah2022NJP ; Kundu2022PRE . The Kuramoto
model has also been generalized to include phases defined not only on the
nodes of a network, but on each edge defining a higher-order interaction
Millan2020PRL ; Ghorbanchian2021Communications ; DeVille2021Chaos  (a
hyperedge when the interactions are encoded in a hypergraph, or a simplex when
they are encoded in a simplicial complex).

For the network Kuramoto model with pairwise interactions, the role of network
community structure is well studied, as it gives rise to a variety of
phenomena including hierarchical synchronization Skardal2012PRE , chimera
states Abram2008PRL , and nonmonotonic synchronization transitions
Restrepo2019PRR . Despite the prevalence of community structure in a wide
range of real-world networks Girvan2002PNAS ; Newman2006PNAS , the effect of
community structure in networks with higher-order interactions is not yet
studied. In this paper we study a generalization of the Kuramoto model that
includes higher-order interactions and community structure. We apply the Ott-
Antonsen ansatz Ott2008Chaos ; Ott2009Chaos  to obtain a low-dimensional
system of equations. Focusing on the case of two communities, we find that,
depending on the strength of the community structure, the system admits
multiple coexisting stable synchronization states. In addition to both
communities being incoherent and both communities being synchronized and in-
phase synchronized with one another, i.e., entrained with a macroscopic phase
difference $\phi=0$, a handful of more complicated states exist. In one such
state one community is synchronized while the other is (nearly) incoherent.
Another state is characterized by both communities being synchronized, but
each community is anti-phase synchronized with one another, i.e., forming two
clusters with a phase difference $\phi=\pi$. Finally, when inter-community
coupling is negative a surprising state emerges where both communities are
synchronized, but are neither in-phase nor anti-phase, rather the two
communities form two clusters with a phase difference $0<\phi<\pi/2$. These
states might find relevance in systems of oscillators where higher-order
interactions are present, such as brain rhythms.

The paper is organized as follows. In Sec. II we present our model and
relevant notation. In Sec. III we derive a low-dimensional description of the
dynamics by using the Ott-Antonsen ansatz. In Sec. IV we study in detail the
case of two communities by analyzing particular cases and validate our
theoretical analysis with numerical results. In Sec. V we present our
conclusions.

##  II Governing equations

We consider the dynamics of a large ensemble of $N$ phase oscillators
organized into $C$ communities. In general we assume that each community
$\sigma$ consists of $N_{\sigma}$ oscillators, with $\sigma=1,\dots,C$
indexing each community, so that $\sum_{\sigma=1}^{C}N_{\sigma}=N$. For
simplicity, we consider both dyadic and triadic (i.e., pair-wise and triplet)
interactions, although our results could be extended to interactions involving
more nodes. Following Ref. Skardal2020CommPhys , the dynamics of each
oscillator is given by

$\displaystyle\dot{\theta}_{i}^{\sigma}$
$\displaystyle=\omega_{i}^{\sigma}+\sum_{\sigma^{\prime}=1}^{C}\frac{K_{1}^{\sigma\sigma^{\prime}}}{N_{\sigma^{\prime}}}\sum_{j=1}^{N_{\sigma^{\prime}}}\sin(\theta_{j}^{\sigma^{\prime}}-\theta_{i}^{\sigma})$
$\displaystyle+\sum_{\sigma^{\prime}=1}^{C}\sum_{\sigma^{\prime\prime}=1}^{C}\frac{K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}}{N_{\sigma^{\prime}}N_{\sigma^{\prime\prime}}}\sum_{j=1}^{N_{\sigma^{\prime}}}\sum_{l=1}^{N_{\sigma^{\prime\prime}}}\sin(2\theta_{j}^{\sigma^{\prime}}-\theta_{l}^{\sigma^{\prime\prime}}-\theta_{i}^{\sigma}),$
(1)

where $\theta_{i}^{\sigma}$ and $\omega_{i}^{\sigma}$ denote the phase and
natural frequency of oscillator $i$ in community $\sigma$,
$K_{1}^{\sigma\sigma^{\prime}}$ is the dyadic coupling strength between a pair
of oscillators in communities $\sigma$ and $\sigma^{\prime}$, respectively,
and $K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}$ is the tetradic
coupling strength between a triplet of oscillators in communities $\sigma$,
$\sigma^{\prime}$, and $\sigma^{\prime\prime}$, respectively. We emphasize
that this all-to-all topology with community-wise uniformity in coupling
eliminates the complications of both intra- and inter-community structures
that may blur the fundamental effects of community structure. We will assume
that the natural frequencies of oscillators in each community $\sigma$ are
drawn from the distribution $g^{\sigma}(\omega)$, which we discuss later.

Next, we introduce the generalized community-wise order parameters

$\displaystyle
z_{q}^{\sigma}=r_{q}^{\sigma}e^{i\psi_{q}^{\sigma}}=\sum_{j=1}^{N_{\sigma}}e^{qi\theta_{j}^{\sigma}}.$
(2)

Note that for the case of $q=1$ we recover the classical Kuramoto order
parameter $z^{\sigma}=r^{\sigma}e^{i\psi^{\sigma}}$ for each community, whose
magnitude $r^{\sigma}$ quantifies the degree of synchronization in community
$\sigma$, while higher-order modes, e.g., $q=2$, measure cluster
synchronization Skardal2011PRE . Using the generalized order parameters in Eq.
(2) for $q=1$ and $2$ we may rewrite Eq. (1) as

$\displaystyle\dot{\theta}_{i}^{\sigma}=\omega_{i}^{\sigma}+\frac{1}{2i}\left(He^{-i\theta_{i}^{\sigma}}-H^{*}e^{i\theta_{i}^{\sigma}}\right),$
(3)

where

$\displaystyle
H=\left(\sum_{\sigma^{\prime}=1}^{C}K_{1}^{\sigma\sigma^{\prime}}z^{\sigma^{\prime}}\right)+\left(\sum_{\sigma^{\prime}=1}^{C}\sum_{\sigma^{\prime\prime}=1}^{C}K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}z_{2}^{\sigma^{\prime}}z^{\sigma^{\prime\prime}*}\right),$
(4)

where $*$ denotes the complex conjugate and we let
$z^{\sigma}=z_{1}^{\sigma}$.

##  III Reduced Dynamics

Next we derive the reduced macroscopic dynamics of the system. Taking the
thermodynamic limit $N\to\infty$ in such a way that the relative sizes of the
communities for each $\sigma$, $N_{\sigma}/N$, remain constant, we introduce
$C$ density functions $f^{\sigma}(\theta,\omega,t)$ that describe, for
community $\sigma$ at time $t$, the fraction of oscillators with phase and
frequency in the infinitesimal intervals $[\theta,\theta+d\theta)$ and
$[\omega,\omega+d\omega)$, respectively. First, due to conservation of
oscillators each $f^{\sigma}$ must satisfy the continuity equation

$\displaystyle\frac{\partial f^{\sigma}}{\partial
t}+\frac{\partial}{\partial\theta}\left(f^{\sigma}\dot{\theta}\right)=0.$ (5)

Second, since the frequency distribution $g^{\sigma}(\omega)$ is fixed the
Fourier series for $f^{\sigma}$ must take the form

$\displaystyle
f^{\sigma}(\theta,\omega,t)=\frac{g^{\sigma}(\omega)}{2\pi}\left(1+\sum_{n=1}^{\infty}\hat{f}_{n}^{\sigma}(\omega,t)e^{in\theta}+\text{c.c.}\right),$
(6)

where c.c. indicates the complex conjugate of the preceding term. Following
the work of Ott and Antonsen, we propose the ansatz where the Fourier
coefficients decay geometrically, i.e., there is some function
$\alpha^{\sigma}(\omega,t)$ for which
$\hat{f}_{n}^{\sigma}(\omega,t)=[\alpha^{\sigma}(\omega,t)]^{n}$. Remarkably,
using this ansatz in Eq. (6) and inserting it into the continuity equation,
Eq. (5), all modes collapse onto a single differential equation for
$\alpha^{\sigma}$, namely

$\displaystyle\frac{\partial}{\partial
t}\alpha^{\sigma}=-i\omega\alpha^{\sigma}+\frac{1}{2}\left(H^{*}-H\alpha^{\sigma*2}\right).$
(7)

To connect the dynamics of $\alpha^{\sigma}$ with those of the order parameter
$z^{\sigma}$, we now note that in the thermodynamic limit we may write

$\displaystyle
z_{q}^{\sigma}=\int_{-\infty}^{\infty}\int_{0}^{2\pi}f^{\sigma}(\theta,\omega,t)e^{qi\theta(t)}d\theta
d\omega,$ (8)

and after integrating over $\theta$ we have that

$\displaystyle
z_{q}^{\sigma}=\int_{-\infty}^{\infty}[\alpha^{\sigma*}(\omega,t)]^{q}g^{\sigma}(\omega)d\omega,$
(9)

Finally, by assuming that the natural frequency distributions are Lorentzian,
i.e.,
$g^{\sigma}(\omega,t)=\Delta^{\sigma}/\\{\pi[(\Delta^{\sigma})^{2}+(\omega-\omega_{0}^{\sigma})^{2}]\\}$
we may evaluate Eq. (9) using the simple pole at
$\omega=\omega_{0}^{\sigma}-i\Delta^{\sigma}$ and the Cauchy residue theorem ,
yielding

$\displaystyle
z^{\sigma}(t)=\alpha^{\sigma*}(\omega_{0}^{\sigma}-i\Delta^{\sigma},t),$ (10)

and

$\displaystyle z_{2}^{\sigma}(t)=[z^{\sigma}(t)]^{2}$ (11)

Thus, by evaluating Eq. (7) at $\omega=\omega_{0}^{\sigma}-i\Delta^{\sigma}$
and taking a complex conjugate, we obtain the following closed system for the
$C$ order parameters:

$\displaystyle\dot{z}^{\sigma}=i\omega_{0}^{\sigma}z^{\sigma}-\Delta^{\sigma}z^{\sigma}+$
$\displaystyle\frac{1}{2}\left\\{\left(\sum_{\sigma^{\prime}=1}^{C}K_{1}^{\sigma\sigma^{\prime}}z^{\sigma^{\prime}}\right)+\left(\sum_{\sigma^{\prime}=1}^{C}\sum_{\sigma^{\prime\prime}=1}^{C}K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}(z^{\sigma^{\prime}})^{2}z^{\sigma^{\prime\prime}*}\right)\right.$
$\displaystyle\left.-\left[\left(\sum_{\sigma^{\prime}=1}^{C}K_{1}^{\sigma\sigma^{\prime}}z^{\sigma^{\prime}*}\right)+\left(\sum_{\sigma^{\prime}=1}^{C}\sum_{\sigma^{\prime\prime}=1}^{C}K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}(z^{\sigma^{\prime}*})^{2}z^{\sigma^{\prime\prime}}\right)\right](z^{\sigma})^{2}\right\\}.$
(12)

Lastly, it is convenient to treat the macroscopic dynamics described by Eq.
(12) in polar coordinates, which are given by

$\displaystyle\dot{r}^{\sigma}$
$\displaystyle=-\Delta^{\sigma}r^{\sigma}+\frac{1-(r^{\sigma})^{2}}{2}\left[\sum_{\sigma^{\prime}=1}^{C}K_{1}^{\sigma\sigma^{\prime}}r^{\sigma^{\prime}}\cos(\psi^{\sigma^{\prime}}-\psi^{\sigma})+\sum_{\sigma^{\prime}=1}^{C}\sum_{\sigma^{\prime\prime}=1}^{C}K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}(r^{\sigma^{\prime}})^{2}r^{\sigma^{\prime}}\cos(2\psi^{\sigma^{\prime}}-\psi^{\sigma^{\prime\prime}}-\psi^{\sigma})\right],$
(13) $\displaystyle\dot{\psi}^{\sigma}$
$\displaystyle=\omega_{0}^{\sigma}+\frac{1+(r^{\sigma})^{2}}{2r^{\sigma}}\left[\sum_{\sigma^{\prime}=1}^{C}K_{1}^{\sigma\sigma^{\prime}}r^{\sigma^{\prime}}\sin(\psi^{\sigma^{\prime}}-\psi^{\sigma})+\sum_{\sigma^{\prime}=1}^{C}\sum_{\sigma^{\prime\prime}=1}^{C}K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}(r^{\sigma^{\prime}})^{2}r^{\sigma^{\prime}}\sin(2\psi^{\sigma^{\prime}}-\psi^{\sigma^{\prime\prime}}-\psi^{\sigma})\right].$
(14)

Before proceeding to a more in-depth analysis of the reduced equations given
in Eqs. (13) and (14), we discuss more precisely the nature of coupling
strengths in light of community structure and the effect on higher-order
interactions. The implication of community structure on dyadic interactions is
relatively straightforward, namely, coupling between pairs of oscillators
within the same community must be stronger than coupling between oscillators
in different communities, i.e., for $\sigma\neq\sigma^{\prime}$ we have that

$\displaystyle K_{1}^{\sigma\sigma}\geq K_{1}^{\sigma\sigma^{\prime}},$ (15)

where equality is only obtained in the limit where community structure
vanishes. On the other hand, the implication of community structure in the
context of triadic coupling is more complicated. Here we will assume that the
effects of community structure on triadic coupling act in such a way that the
triadic coupling strengths generically decrease as the oscillators belong to
different communities, i.e., for $\sigma\neq\sigma^{\prime}$,
$\sigma\neq\sigma^{\prime\prime}$, and
$\sigma^{\prime}\neq\sigma^{\prime\prime}$ we have that

$\displaystyle K_{2}^{\sigma\sigma\sigma}\geq
K_{2}^{\sigma\sigma\sigma^{\prime}}=K_{2}^{\sigma\sigma^{\prime}\sigma}\geq
K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime}}\geq
K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}.$ (16)

Figure 1: Community structure and higher-order coupling. In the case of two
communities, the coupling strength associated to an oscillator in community
$\sigma$ via a hyperedge of size $m$ is given by $\alpha^{j}K_{m-1}$, where
$j$ is the number of oscillators in the hyperedge not belonging to community
$\sigma$.

Note that the left- and right-most terms represent the triadic coupling
strength between oscillators that are, respectively, all in the same community
and all in different communities. On the other hand, the three middle terms
represent triadic coupling between a triplet of oscillators that collectively
belong to two communities. We take these to be all equal except for
$K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime}}$, which we assume is small since
the affected oscillator belongs to a different community than the other two.
To capture this hierarchy of coupling in a simple way we introduce a parameter
$\alpha\in[-1,1]$ that scales coupling strengths, namely for
$\sigma\neq\sigma^{\prime}$, $\sigma\neq\sigma^{\prime\prime}$, and
$\sigma^{\prime}\neq\sigma^{\prime\prime}$, dyadic coupling strengths are
given by

$\displaystyle
K_{1}^{\sigma\sigma}=K_{1},~{}~{}~{}\text{and}~{}~{}~{}K_{1}^{\sigma\sigma^{\prime}}=\alpha
K_{1},$ (17)

and triadic coupling strengths are given by

$\displaystyle
K_{2}^{\sigma\sigma\sigma}=K_{2},~{}~{}~{}K_{2}^{\sigma\sigma\sigma^{\prime}}=K_{2}^{\sigma\sigma^{\prime}\sigma}=\alpha
K_{2},$ $\displaystyle
K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime}}=\alpha^{2}K_{2},~{}~{}~{}\text{and}~{}~{}~{}K_{2}^{\sigma\sigma^{\prime}\sigma^{\prime\prime}}=\alpha^{3}K_{2}.$
(18)

This coupling structure is illustrated schematically for the case of two
communities in Fig. 1. Note that in such a case with only two communities,
which will be our focus for the remainder of this paper, the weakest form of
triadic coupling plays no role.

Figure 2: Stable states. Schematic illustration of the five possible stable
states in the two community case: (a) in-phase synchronized, (b) anti-phase
synchronized, (c) incoherent-synchronized, (d) incoherent, and (e) skew-phase
synchronized.

##  IV Two communities

We now turn to present a detailed analysis of the two community case, for
which the macroscopic dynamics are described by the amplitudes, $r_{1}$ and
$r_{2}$, and phases, $\psi_{1}$ and $\psi_{2}$, of the respective communities.
In fact, the dynamics can be further reduced to three dimensions by
introducing the phase difference $\phi=\psi_{1}-\psi_{2}$, which yields the
following system of equations

$\displaystyle\dot{r}_{1}$
$\displaystyle=-\Delta_{1}r_{1}+\frac{1-r_{1}^{2}}{2}\left[K_{1}r_{1}+\alpha
K_{1}r_{2}\cos(\phi)+K_{2}r_{1}^{3}+\alpha
K_{2}r_{1}^{2}r_{2}\cos(\phi)+\alpha
K_{2}r_{2}^{2}r_{1}\cos(2\phi)+\alpha^{2}K_{2}r_{2}^{3}\cos(\phi)\right],$
(19) $\displaystyle\dot{r}_{2}$
$\displaystyle=-\Delta_{2}r_{2}+\frac{1-r_{2}^{2}}{2}\left[\alpha
K_{1}r_{1}\cos(\phi)+K_{1}r_{2}+\alpha^{2}K_{2}r_{1}^{3}\cos(\phi)+\alpha
K_{2}r_{1}^{2}r_{2}\cos(2\phi)+\alpha
K_{2}r_{2}^{2}r_{1}\cos(\phi)+K_{2}r_{2}^{3}\right],$ (20)
$\displaystyle\dot{\phi}$
$\displaystyle=\delta\omega-\frac{\alpha}{2}\left[K_{1}\left(\frac{r_{2}^{2}(1+r_{1}^{2})+r_{1}^{2}(1+r_{2}^{2})}{r_{1}r_{2}}\right)\sin(\phi)+K_{2}\left(\frac{r_{1}^{2}r_{2}^{2}(1+r_{1}^{2})+r_{1}^{2}r_{2}^{2}(1+r_{2}^{2})}{r_{1}r_{2}}\right)\sin(-\phi)\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+K_{2}\left(\frac{r_{2}^{3}r_{1}(1+r_{1}^{2})+r_{1}^{3}r_{2}(1+r_{2}^{2})}{r_{1}r_{2}}\right)\sin(2\phi)+\alpha
K_{2}\left(\frac{r_{2}^{4}(1+r_{1}^{2})+r_{1}^{4}(1+r_{2}^{2})}{r_{1}r_{2}}\right)\sin(\phi)\right],$
(21)

where $\delta\omega=\omega_{0}^{1}-\omega_{0}^{2}$. Since our focus is on the
effect of community structure, we assume that the aggregate local dynamics of
each community is the same, so that $\Delta_{1}=\Delta_{2}=\Delta$ and
$\omega_{0}^{1}=\omega_{0}^{2}$, yielding $\delta\omega$=0. In this case we
find five different qualitatively-different steady-states, depending on the
choice of parameters and initial conditions:

  1. (i)

An in-phase synchronized state with $r_{1}=r_{2}>0$ and $\phi=0$.

  2. (ii)

An anti-phase synchronized state with $r_{1}=r_{2}>0$ and $\phi=\pi$.

  3. (iii)

An incoherent-synchronized state with $r_{1}>0$, $r_{2}\approx 0$ and
$\phi=0$.

  4. (iv)

An incoherent state with $r_{1}=r_{2}=0$.

  5. (v)

A skew-phase synchronized state with $r_{1}=r_{2}>0$ and $0<\phi<\pi/2$.

These five states are illustrated in Fig. 2(a)–(e), respectively, with red and
blue circles denoting the locations of the order parameters
$r_{1}e^{\psi_{1}}$ and $r_{2}e^{i\psi_{2}}$ in each case, and purple circles
denoting overlapping order parameters. As we will see below, the four states
(a)–(d) are all supported by the system dynamics when interaction between
communities is cooperative, i.e., $\alpha>0$. On the other hand, we find that
when interaction between communities is contrarian, i.e., $\alpha<0$, the
incoherent and skew-phase synchronized states are supported.

Before moving to our analytical results, we illustrate the dynamics and the
multistability between these different states with some numerical simulations.
Finding the richest dynamics occurring for a combination of small dyadic
coupling and large triadic coupling, in Fig. 3(a) and (b) we plot the
trajectories of the dynamics of the two order parameters in the complex unit
disc for (a) four and (b) two different initial conditions using $\Delta=1$,
$K_{1}=0.1$, $K_{2}=10$, and (a) $\alpha=0.1$ and (b) $-0.1$. Initial and
final states are plotted as open and filled circles, respectively, and
trajectories that end up at the different states are labelled: (i) in-phase
synchronized, (ii) anti-phase synchronized, (iii) incoherent-synchronized,
(iv) incoherent, and (v) skew-phase synchronized. These trajectories are
plotted in solid blue, dashed red, dot-dashed green, dotted purple, and solid
orange, respectively. Note that for these parameters all four states (i)–(iv)
are stable for $\alpha=0.1$ and both states (iv) and (v) are stable for
$\alpha=-0.1$. Grey lines in Fig. 3(a) correspond to direct simulations of the
full system (1) with a total of $N=10^{5}$ oscillators, i.e.,
$N_{\sigma}=5\times 10^{4}$ in each community, which show good agreement with
the reduced dynamics in Eqs. (13)-(14).

Figure 3: Multistability in systems with higher-order interactions and
community structure. Illustration of the possible states using trajectories of
order parameters on the complex unit disc for (a) positive and (b) negative
$\alpha$. Parameters are given by $\Delta=1$, $K_{1}=0.1$, $K_{2}=10$, and
$\alpha=0.1$ or $-0.1$. The labels (i), (ii), (iii), (iv), and (v) correspond
to the in-phase synchronized state, the anti-phase synchronized state, the
incoherent-synchronized state, the incoherent state, and the skew-phase
synchronized state. Gray curves describe the dynamics of direct simulations
with $N=10^{5}$ total oscillators.

Moving forward, we aim to understand these states and their dynamics using
analytical methods. To this end, we note that the collection of five possible
states lies on two manifolds within the full phase space. The first of these
we call the $r$-manifold, and is characterized by $r_{1}=r_{2}=r$, i.e., both
communities are characterized by exactly the same degree of synchronization.
The second manifold, which we call the $\phi$-manifold, is characterized by
$\phi=0$, i.e., the two communities share the same mean phase. Note that
states (i), (ii), (iv), and (v) lie on the $r$-manifold, while states (i),
(iii), and (iv) lie on the $\phi$-manifold. We will now focus on these two
manifolds to give analytical insight into the system dynamics.

Figure 4: Stable and unstable branches: in-phase and anti-phase synchronized
states. Fixed points $r_{a}^{+}$, $r_{b}^{+}$ (solid lines) and $r_{a}^{-}$,
$r_{b}^{-}$ (dashed lines) as a function of $K_{2}$ for $K_{1}=0$, $\Delta=1$,
and $\alpha=0.1$, $0.45$, and $0.7$.

###  IV.1 The $r$-manifold

Beginning with the case $r_{1}=r_{2}=r$, we also make the simplifying
assumption that $K_{1}=0$, i.e., coupling is solely triadic. As we will see,
this suffices to recover all five steady-states illustrated above. On this
manifold the dynamics then reduce to

$\displaystyle\dot{r}$ $\displaystyle=-\Delta r+\frac{K_{2}r^{3}(1-r^{2})}{2}$
$\displaystyle\times\left[1+\alpha\cos(\phi)+\alpha\cos(2\phi)\alpha^{2}\cos(\phi)\right],$
(22) $\displaystyle\dot{\phi}$
$\displaystyle=-K_{2}r^{2}\alpha(1+r^{2})\left[-\sin(\phi)+\sin(2\phi)+\alpha\sin(\phi)\right],$
(23)

for which, in addition to the incoherent state $r=0$, we find fixed points
$(r_{a}^{-},0)$, $(r_{a}^{+},0)$, $(r_{b}^{-},\pi)$, $(r_{b}^{+},\pi)$,
$(r_{b}^{-},\pm\phi^{\text{skew}})$, and $(r_{b}^{+},\pm\phi^{\text{skew}})$,
where

$\displaystyle
r_{a}^{\pm}=\sqrt{\frac{1}{2}\pm\frac{\sqrt{(\alpha+1)^{2}K_{2}-8}}{2(\alpha+1)\sqrt{K_{2}}}},$
(24) $\displaystyle
r_{b}^{\pm}=\sqrt{\frac{1}{2}\pm\frac{\sqrt{\left(\alpha^{2}-1\right)K_{2}+8}}{2\sqrt{\alpha^{2}-1}\sqrt{K_{2}}}},$
(25)
$\displaystyle\phi^{\text{skew}}=\pm\arccos\left(\frac{1-\alpha}{2}\right).$
(26)

Note first that when $\alpha=0$ (i.e., the communities are isolated from one
another), the branches $r_{a}^{\pm}$ and $r_{b}^{\pm}$ are equivalent. Next,
when $0<\alpha<1$ the branches $r_{a}^{\pm}$ and $r_{b}^{\pm}$ appear at
saddle-node bifurcations, respectively, at

$\displaystyle
K_{2}^{a}=\frac{8}{(1+\alpha)^{2}}~{}~{}~{}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}~{}~{}~{}K_{2}^{b}=\frac{8}{1-\alpha^{2}}.$
(27)

After the saddle-node bifurcation occurs, the branches $r_{a}^{+}$ and
$r_{b}^{+}$ are stable fixed points that correspond to the in-phase
synchronized and anti-phase synchronized states, respectively, with
$r_{a}^{-}$ and $r_{b}^{-}$ being their unstable counterparts. (By
linearization, see Appendix A, it can be shown that the $r_{a}^{+}$ and
$r_{b}^{+}$ are stable and $r_{a}^{-}$ and $r_{b}^{-}$ are unstable.) These
branches are illustrated in Fig. 4, depicting $r\cos(\phi)$ as a function of
$K_{2}$ for $\alpha=0.1$, $0.45$, and $0.7$ (blue, red, and green,
respectively). Note that the two branches occur at similar values of $K_{2}$
for small $\alpha$, i.e., strong community structure, and as $\alpha$
increases, i.e., community structure becomes weaker, the branches occur at
$K_{2}$ values that become more different, as larger $\alpha$ promotes in-
phase synchronization but pushes back the anti-phase synchronized branch.

Figure 5: Stability diagram for the $r$-manifold. Diagram illustrating the
stability of different states in the $r$-manifold for $K_{1}=0$ and
$\Delta=1$.

Lastly, for $-1<\alpha<0$ the system displays different synchronized behavior.
Rather than supporting the in-phase and anti-phase synchronized states, these
become unstable and the skew-phase synchronized state becomes stable, taking
the same local degrees of synchronization as the anti-phase synchronized
state, i.e., $r=r_{b}^{\pm}$ (again, with $r_{b}^{+}$ being the stable fixed
point), but interestingly, the two organize into two clusters that lie at a
phase difference $\phi^{\text{skew}}\in(0,\pi/2)$. The system transitions
directly from the incoherent state to a stable state where both communities
are synchronized with the same value of the order parameter $r_{b}^{+}$, but
have a phase difference $0<\phi^{\text{off}}<\pi$ [see Eq. (26)].

The stability diagram for the $r$-manifold for the case of $K_{1}=0$ is
summarized in Fig. 5, which shows the stable steady states as a function of
the parameters $K_{2}$ and $\alpha$. For positive values of $\alpha$, as
$K_{2}$ is increased, the incoherent state is first the only stable state,
after which the in-phase synchronized state is born at $K_{2}^{a}$ in a
saddle-node bifurcation, and subsequently the anti-phase synchronized state is
born at $K_{2}^{b}$, again in a saddle-node bifurcation. On the other hand,
for negative $\alpha$, as $K_{2}$ increases, the incoherent state is, again,
first the only stable state, after which only one saddle-node bifurcation
occurs at $K_{2}^{b}$, giving rise to the skew-phase synchronized state.

###  IV.2 The $\phi$-manifold

Next, we proceed to analyze the $\phi$-manifold, specifically with the aim of
gaining some understanding of the incoherent-synchronized state, which is the
only observed steady-state not on the $r$-manifold analyzed above. To this
end, focusing again on the case of $K_{1}=0$, we assume that $\phi=0$,
yielding the reduced set of equations

$\displaystyle\dot{r}_{1}$ $\displaystyle=-\Delta r_{1}+\frac{1-r_{1}^{2}}{2}$
$\displaystyle\times(K_{2}r_{1}^{3}+\alpha K_{2}r_{1}^{2}r_{2}+\alpha
K_{2}r_{2}^{2}r_{1}+\alpha^{2}K_{2}r_{2}^{3}),$ (28)
$\displaystyle\dot{r}_{2}$ $\displaystyle=-\Delta r_{2}+\frac{1-r_{2}^{2}}{2}$
$\displaystyle\times(\alpha^{2}K_{2}r_{1}^{3}+\alpha
K_{2}r_{1}^{2}r_{2}+\alpha K_{2}r_{2}^{2}r_{1}+K_{2}r_{2}^{3}).$ (29)

While an exact analytical expression for the incoherent-synchronized state is
difficult to obtain, an approximation is tractable using a perturbation
analysis. First, we note that, provided that $K_{2}\geq 8$, in the limit
$\alpha\to 0$ (i.e., when community structure is so strong that the two
communities are isolated from one another), both $r_{1,2}=0$ and
$\sqrt{\frac{1}{2}+\frac{\sqrt{K_{2}-8\Delta}}{2\sqrt{K_{2}}}}$ are stable.
Moreover, the incoherent-synchronized state is only stable for sufficiently
strong community structure, thereby making $\alpha$ a suitable perturbation
parameter. We then consider a perturbation of the steady-state $r_{1}>0$,
$r_{2}=0$ by the parameter $\alpha$, i.e., searching for solutions with
expansions

$\displaystyle
r_{1}=r_{1}^{(0)}+r_{1}^{(1)}\alpha+r_{1}^{(2)}\alpha^{2}+r_{1}^{(3)}\alpha^{3}+\mathcal{O}(\alpha^{4}),$
(30) $\displaystyle
r_{2}=r_{2}^{(0)}+r_{2}^{(1)}\alpha+r_{2}^{(2)}\alpha^{2}+r_{2}^{(3)}\alpha^{3}+\mathcal{O}(\alpha^{4}).$
(31)

Inserting Eqs. (30) and (31) into Eqs. (28) and (29) and setting
$\dot{r}_{1},\dot{r}_{2}=0$ yields

$\displaystyle
2\Delta\left(r_{1}^{(0)}+r_{1}^{(1)}\alpha+r_{1}^{(2)}\alpha^{2}+r_{1}^{(3)}\alpha^{3}\right)$
$\displaystyle=\left[K_{2}r_{1}^{(0)3}(1-r_{1}^{(0)2})\right]$
$\displaystyle+\left[r_{1}^{(1)}(3K_{2}r_{1}^{(0)2}-5K_{2}r_{1}^{(0)4})\right]\alpha$
$\displaystyle+\left[K_{2}r_{1}^{(0)}\left((3-10r_{1}^{(0)2})r_{1}^{(1)2}+r_{1}^{(0)}(3r_{1}^{(2)}-5r_{1}^{(0)2}r_{1}^{(2)}+r_{2}^{(1)}-r_{1}^{(0)2}r_{2}^{(1)})\right)\right]\alpha^{2}$
$\displaystyle+\left[K_{2}((1-10r_{1}^{(0)2})r_{1}^{(1)3}+2r_{1}^{(0)}r_{1}^{(1)}((3-10r_{1}^{(0)2})r_{1}^{(2)}+r_{2}^{(1)}-2r_{1}^{(0)2}r_{2}^{(1)})\right.$
$\displaystyle~{}~{}~{}~{}+\left.r_{1}^{(0)}(r_{2}^{(1)2}+r_{1}^{(0)}((3-5r_{1}^{(0)2})r_{1}^{(3)}+r_{2}^{(2)}-r_{1}^{(0)}(r_{2}^{(1)2}+r_{1}^{(0)}r_{2}^{(2)}))))\right]\alpha^{3},$
(32) $\displaystyle
2\Delta\left(r_{2}^{(1)}\alpha+r_{2}^{(2)}\alpha^{2}+r_{2}^{(3)}\alpha^{3}\right)$
$\displaystyle=\left[K_{2}r_{1}^{(0)2}(r_{1}^{(0)}+r_{2}^{(1)})\right]\alpha^{2}$
$\displaystyle+\left[K_{2}(r_{2}^{(1)3}+r_{1}^{(0)}r_{2}^{(1)}(2r_{1}^{(1)}+r_{2}^{(1)})+r_{1}^{(0)2}(3r_{1}^{(1)}+r_{2}^{(2)}))\right]\alpha^{3},$
(33)

where we’ve used that $r_{2}^{(0)}=0$. Collecting terms at different orders of
$\alpha$, we solve to obtain

$\displaystyle r_{1}^{(0)}$
$\displaystyle=\sqrt{\frac{1}{2}+\frac{\sqrt{K_{2}-8\Delta}}{2\sqrt{K_{2}}}},$
(34) $\displaystyle r_{1}^{(1)}$ $\displaystyle=0,$ (35) $\displaystyle
r_{1}^{(2)}$
$\displaystyle=-\frac{K_{2}r_{1}^{(0)2}(1-r_{1}^{(0)2})r_{2}^{(1)}}{2\Delta+K_{2}r_{1}^{(0)2}(3-5r_{1}^{(0)2})}=0,$
(36) $\displaystyle r_{1}^{(3)}$
$\displaystyle=-\frac{K_{2}r_{1}^{(0)}(1-r_{1}^{(0)2})(r_{2}^{(1)2}+r_{1}^{(0)}r_{2}^{(2)})}{2\Delta+K_{2}r_{1}^{(0)2}(3-5r_{1}^{(0)2})}=\frac{\sqrt{K_{2}}}{2\sqrt{2}\sqrt{K_{2}-8\Delta}}\sqrt{1+\frac{\sqrt{K_{2}-8\Delta}}{\sqrt{K_{2}}}}.$
(37) $\displaystyle r_{2}^{(0)}$ $\displaystyle=0,$ (38) $\displaystyle
r_{2}^{(1)}$ $\displaystyle=0,$ (39) $\displaystyle r_{2}^{(2)}$
$\displaystyle=\frac{K_{2}r_{1}^{(0)2}(r_{1}^{(0)}+r_{2}^{(1)})}{2\Delta}=\frac{K_{2}}{4\sqrt{2}\Delta}\left(1+\frac{\sqrt{K_{2}-8\Delta}}{\sqrt{K_{2}}}\right)^{3/2},$
(40) $\displaystyle r_{2}^{(3)}$
$\displaystyle=\frac{K_{2}(r_{2}^{(1)2}(r_{1}^{(0)}+r_{2}^{(1)})+r_{1}^{(0)2}r_{2}^{(2)})}{2\Delta}=\frac{K_{2}^{2}}{16\sqrt{2}\Delta^{2}}\left(1+\frac{\sqrt{K_{2}-8\Delta}}{\sqrt{K_{2}}}\right)^{5/2}.$
(41)

Figure 6: Incoherent-synchronized state: perturbation theory. Comparison of
the perturbation theory from the text (solid curves) to direct simulation of
the low dimensional dynamics for the incoherent-synchronized state using
$K_{1}=0$, $K_{2}=10$, and $\Delta=1$.

Note that, with $r_{1}^{(1)}$, $r_{1}^{(2)}$, and $r_{2}^{(1)}$ all zero, the
leading order correction to the order parameters $r_{1}$ and $r_{2}$ for the
synchronized and incoherent come at cubic and quadratic orders in $\alpha$,
respectively. In Fig. 6 we compare the perturbation theory obtained above
(solid curves) to direct simulations of the low dimensional dynamics (circles)
using $K_{2}=10$ and $\Delta=1$ for both $r_{1}$ and $r_{2}$ in panels (a) and
(b). We note strong agreements until near the critical value of $\alpha$ where
the incoherent-synchronized state disappears, denoted by the vertical dashed
lines.

##  V Discussion

We have studied the synchronization of coupled phase oscillators with
community structure and higher-order interactions. When there is no community
structure present, as in Refs. Skardal2020CommPhys , higher-order interactions
result in bistable states where both the incoherent and synchronized states
are stable. Therefore, one would naturally expect that for sufficiently strong
community structure there would be stable states where both communities are
synchronized, both are incoherent, and one is synchronized while the other is
incoherent. While these states are, indeed, observed for sufficiently strong
community structure, we have found that two additional states can appear. The
first is a state where both communities are synchronized but are clustered
opposite one another. This anti-phase synchronized state requires both strong
enough community structure and higher-order coupling strength (see Fig. 5).
Since community structure and higher-order interactions are observed in brain
oscillations, we hypothesize that this mechanism could be relevant for
producing anti-phase oscillations in brain rhythms. In addition, we found an
additional novel state where both communities are synchronized but oscillate
with a phase difference that can be tuned by the amount of community structure
[see Eq. (26)]. This skew-phase state is observed only when the sign of
coupling via triadic interactions alternates depending on the membership of
oscillators in the triad to different communities.

In order to focus on the more novel results, we studied particular cases of
the full model in Eq. (1). The most important simplifications were that we
considered only triadic interactions, two communities, and that we assumed
that the two communities had identical frequency distributions
($\Delta_{1}=\Delta_{2}$ and $\delta\omega=0$). Here we hypothesize on what
should be expected if these assumptions are removed, but leave a detailed
study for future work. The addition of a nonzero mean frequency difference
between communities ($\delta\omega>0$) should produce bifurcations to standing
wave solutions with two counterrotating groups of oscillators (as in, e.g.,
Martens2009PRE ). The hierarchical coupling structure in Eqs. (18) could be
generalized in a straightforward manner to cases with larger interactions and
more communities, and the subsequent analysis would follow in a similar way.
The presence of more communities should result in stable states where some
communities are synchronized while the rest are almost incoherent. How higher-
order interactions and community structure would affect the relative phases of
the synchronized communities or if additional stable states could be present
is not clear, however, and is left for future study. In summary, by analyzing
the simplest phase oscillator model that includes community structure and
higher-order interactions, we have found that multiple synchronized states can
coexist, including anti-phase and skew-phase synchronized states. We
anticipate that the interplay of community structure and higher-order
interactions in less simplified scenarios will result in complex oscillation
dynamics.

##  Appendix A Stability of the in-phase and anti-phase states

In this Appendix we provide detailed calculations for the linear stability of
the in-phase and anti-phase states discussed in Sec. IV. We start by
calculating the eigenvalues Jacobian associated to Eqs. (19)-(21) with
$K_{1}=0$ at the in-phase and anti-phase fixed points.

At $r_{1}=r_{2}=r_{a}^{+}$, $\phi=0$, the eigenvalues of the Jacobian are,
letting $k\equiv K_{2}/K_{2}^{a}$,

$\displaystyle\lambda_{1}$
$\displaystyle=\frac{2\alpha}{\alpha+1}\left(1-4k\right)-\alpha\sqrt{8K_{2}}\sqrt{k-1},$
(42) $\displaystyle\lambda_{2}$
$\displaystyle=\frac{4-2\alpha}{(\alpha+1)}-4k-(\alpha+1)\sqrt{2K_{2}}\sqrt{k-1},$
(43) $\displaystyle\lambda_{3}$
$\displaystyle=4\left(1-k\right)-\frac{1}{2}(\alpha+1)\sqrt{8K_{2}}\sqrt{k-1},$
(44)

which are negative for $\alpha>0$ and $K_{2}>K_{2}^{a}$ ($k>1$), confirming
the linear stability of the fixed point $r_{1}=r_{2}=r_{a}^{+}$, $\phi=0$ .

At $r_{1}=r_{2}=r_{a}^{-}$, $\phi=0$, the eigenvalues are

$\displaystyle\lambda_{1}$
$\displaystyle=\frac{2\alpha\left(-4k+4\sqrt{k-1}\sqrt{k}+1\right)}{\alpha+1},$
(45) $\displaystyle\lambda_{2}$
$\displaystyle=\frac{6}{\alpha+1}-4k+4\sqrt{k-1}\sqrt{k}-2,$ (46)
$\displaystyle\lambda_{3}$ $\displaystyle=4(1-k)+4\sqrt{k-1}\sqrt{k}.$ (47)

For $k>1$ the eigenvalue $\lambda_{3}$ is positive, and therefore the fixed
point $r_{1}=r_{2}=r_{a}^{-}$, $\phi=0$ is unstable. The linear stability and
instability of the fixed points $r_{1}=r_{2}=r_{b}^{+}$, $\phi=\pi$ and
$r_{1}=r_{2}=r_{a}^{-}$, $\phi=\pi$, respectively, can be checked similarly.

###### Acknowledgements.

PSS acknowledges NSF grant MCB-2126177.

## References

  * [1] S. H. Strogatz, Sync: the Emerging Science of Spontaneous order (Hypernion, 2003).
  * [2] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, 2003).
  * [3] Petri, G. et al. Homological scaffolds of brain functional networks, J. R. Soc. Interface 11, 20140873 (2014).
  * [4] Ashwin, P. & Rodrigues, A. Hopf normal form with $S_{N}$ symmetry and reduction to systems of nonlinearly coupled phase oscillators, Physica D 325, 14 (2016).
  * [5] Z. Lu et al., Resynchronization of circadian oscillators and the east-west asymmetry of jet-lag, Chaos 26, 094811 (2016).
  * [6] M. G. Kitzbichler, M. L. Smith, S. R. Christensen, and Ed Bullmore, Broadband Criticality of Human Brain Network Synchronization, PLOS Computational Biology 5, e1000314 (2009).
  * [7] L. Glass and M. C. Mackey, From Clocks to Chaos: The Rhythms of Life (Princeton University Press, Princeton, 1988).
  * [8] S. Strogatz et al., Theoretical mechanics: crowd synchrony on the Millennium Bridge, Nature 438 43 (2005).
  * [9] K. Wiesenfeld, P. Colet, and S. H. Strogatz, Synchronization transitions in a disordered Josephson series array, Phys. Rev. Lett. 76, 404 (1996).
  * [10] M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Self-organized synchronization in decentralized power grids, Phys. Rev. Lett. 109, 064101 (2012).
  * [11] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, New York, 1984).
  * [12] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Networks beyond pairwise interactions: Structure and dynamics, Phys. Rep. 874, 1 (2020).
  * [13] S. Majhi, M. Perc, and D. Ghosh, Dynamics on higher-order networks: a review, J. R. Soc. Interface 19, 20220043 (2022).
  * [14] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Simplicial models of social contagion, Nat. Commun. 10, 2485 (2019).
  * [15] J. T. Matamalas, S. Gómez, and A. Arenas, Abrupt phase transition of epidemic spreading in simplicial complexes, Phys. Rev. Res. 2, 012049(R) (2019).
  * [16] J. Grilli, G. Barabás, M. Michalska-Smith, and Stefano Allesina, Higher-order interactions stabilize dynamics in competitive network models, Nature 548 210 (2017).
  * [17] T. Carletti, D. Fanelli, and S. Nicoletti, Dynamical systems on hypergraphs, J. Phys. Complex. 1, 035006 (2020).
  * [18] C. Giusti, R. Ghrist, and D. S. Bassett, Two’s company, three (or more) is a simplex, J. Comput. Neurosci. 41, 1 (2016).
  * [19] M. W. Reimann et al. Cliques of neurons bound into cavities provide a missing link between structure and function, Frontiers in Comp. Neuro. 11, 48 (2017).
  * [20] A. E. Sizemore et al. Cliques and cavities in the human connectome, J. Comput. Neurosci. 44, 115 (2018).
  * [21] I. Léon and D. Pazó, Phase reducation beyond the first-order: The case of the mean-field complex Ginzburg-Landau equation, Phys. Rev. E 100, 012211 (2019).
  * [22] P. S. Skardal and A. Arenas, Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes, Phys. Rev. Lett. 122, 248301 (2019).
  * [23] P. S. Skardal and A. Arenas, Higher order interactions in complex networks of phase oscillators promote abrupt synchronization switching, Commun. Phys. 3, 218 (2020).
  * [24] C. Xu, X. Wang, and P. S. Skardal, Bifurcation analysis and structural stability of simplicial oscillator populations, Phys. Rev. Res. 2, 023281 (2020).
  * [25] P. S. Skardal and A. Arenas, Memory selection and information switching in oscillator networks with higher-order interactions, J. Phys. Complex. 2, 015003 (2020).
  * [26] C. Xu and P. S. Skardal, Spectrum of extensive multiclusters in the Kuramoto model with higher-order interactions, Phys. Rev. Res. 3, 013013 (2021).
  * [27] P. S. Skardal, L. Arola-Fernández, D. Taylor, and A. Arenas, Higher-order interactions improve optimal collective dynamics on networks, Phys. Rev. Res. 3, 043193 (2021).
  * [28] M. S. Anwar and D. Ghosh, Intralayer and interlayer synchronization in multiplex network with higher-order interactions, Chaos 32, 033125 (2022).
  * [29] P. S. Skardal and C. Xu, Tiered synchronization in coupled oscillator populations with interaction delays and higher-order interactions, Chaos 32, 053120 (2022).
  * [30] I. León and D. Pazó, Enlarged Kuramoto model: Secondary instability and transition to collective chaos, Phys. Rev. E 105, L042201 (2022).
  * [31] A. D. Kachhvah and S. Jalan, Hebbian plasticity rules abrupt desynchronization in pure simplicial complexes, New J. Phys. 24, 052002 (2022).
  * [32] S. Kundu and D. Ghosh, Higher-order interactions promote chimera states, Phys. Rev. E 105, L042202 (2022).
  * [33] A. P. Millán, J. J. Torres, and G. Bianconi, Explosive higher-order Kuramoto dynamics on simplicial complexes, Phys. Rev. Lett. 124, 218301 (2020).
  * [34] R. Ghorbanchian, J. G. Restrepo, J. J. Torres, and G. Bianconi, Higher-order simplicial synchronization of coupled topological signals, Communications Physics 4 120 (2021).
  * [35] L. DeVille, Consensus on simplicial complexes: Results on stability and synchronization, Chaos 31, 023137 (2021).
  * [36] P. S. Skardal and J. G. Restrepo, Hierarchical synchrony of phase oscillators in modular networks, Phys. Rev. E 85, 016208 (2012).
  * [37] D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, Solvable model for chimera states of coupled oscillators, Phys. Rev. Lett. 101, 084103 (2008)
  * [38] J. G. Restrepo and P. S. Skardal, Competitive suppression of synchronization and nonmonotonic transitions in oscillator communities with distributed time delay, Phys. Rev. Res. 1, 033042 (2019).
  * [39] M. Girvan and M. E. J. Newman, Community structure in social and biological networks, Proc. Natl. Acad. Sci. U.S.A. 99, 7821 (2002).
  * [40] M. E. J. Newman, Modularity and community structure in networks, Proc. Natl. Acad. Sci. U.S.A. 103, 8577 (2006).
  * [41] E. Ott and T. M. Antonsen, Low dimensional behavior of large systems of globally coupled oscillators, Chaos 18, 037113 (2008).
  * [42] E. Ott and T. M. Antonsen, Long time evolution of phase oscillator systems, Chaos 19, 023117 (2009).
  * [43] P. S. Skardal, E. Ott, and J. G. Restrepo, Cluster synchrony in systems of coupled phase oscillators with higher-order coupling, Phys. Rev. E 84, 036208 (2011).
  * [44] E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, and T. M. Antonsen, Exact results for the Kuramoto model with a bimodal frequency distribution, Phys. Rev. E 79, 026204 (2009).