# Order, chaos, and dimensionality transition in a system of swarmalators

Joao U. F. Lizárraga  Instituto de Física Gleb Wataghin, Universidade Estadual
de Campinas, Unicamp 13083-970, Campinas, São Paulo, Brazil  Kevin P. O’Keeffe
Senseable City Lab, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA  Marcus A. M. de Aguiar<EMAIL_ADDRESS>Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, Unicamp
13083-970, Campinas, São Paulo, Brazil

###### Abstract

Similar to sperm, where individuals self-organize in space while also striving
for coherence in their tail swinging, several natural and engineered systems
exhibit the emergence of swarming and synchronization. The arising and
interplay of these phenomena have been captured by collectives of hypothetical
particles named swarmalators, each possessing a position and a phase whose
dynamics are affected reciprocally and also by the space-phase states of their
neighbors. In this work, we introduce a solvable model of swarmalators able to
move in two-dimensional spaces. We show that several static and active
collective states can emerge and derive necessary conditions for each to show
up as the model parameters are varied. These conditions elucidate, in some
cases, the displaying of multistability among states. Notably, in the active
regime, individuals behave chaotically, maintaining spatial correlation under
certain conditions, and breaking it under others on what we interpret as a
dimensionality transition.

††preprint: APS/123-QED

##  I Introduction

Swarming and synchronization are emergent phenomena observed in various
natural systems. Swarming, more noticeable due to its occurrence in physical
space, is evidenced in groups of fishes, birds, insects, bacteria, among
others [1, 2, 3, 4, 5, 6, 7, 8, 9]. Synchronization, sometimes less apparent
due to its nature, is observed in both living and non-living systems,
including groups of humans, frogs, heart cells, neurons, among others [10, 11,
12, 13, 14, 15, 16, 17]. From a theoretical standpoint, both phenomena have
been extensively studied using distinct frameworks based on well-known models
named after Kuramoto [18], Stuart-Landau [19], Couzin [20], and Vicsek [21].
In these models, single particles are represented by their positions or by
periodic internal degrees of freedom dubbed phases. Then, individuals’ spatial
self-organization gives rise to swarming and, in case of phase coherence, to
synchronization.

Despite the success in studying swarming and synchronization independently, a
significant step up has taken place recently, influencing both fields. A
singular type of particles, named swarmalators, have been introduced in [22]
such that, collectively, these can synchronize and swarm as their spatial and
phase dynamics interplay. Several systems where individuals’ phases influence
their positions and vice versa have been spotted in nature. These include
magnetic quincke rollers [23], sperm [24, 25, 26], starfish embryos [27], tree
frogs [28], nematodes [29, 30], and more [31, 32, 33, 34, 35, 36].

The swarmalators model, introduced in [22], is defined by the equations

$\displaystyle\begin{split}\dot{\vec{p}}_{i}&=\vec{\varepsilon}_{i}+\frac{1}{N}\sum_{j=1}^{N}\left[I_{\alpha}(\vec{p}_{ji})F(\theta_{ji})-I_{\rho}(\vec{p}_{ji})\right],\\\
\dot{\theta}_{i}&=\omega_{i}+\frac{1}{N}\sum_{j=1}^{N}H_{\alpha}(\theta_{ji})G(\vec{p}_{ji}),\end{split}$
(1)

where $\vec{p}_{i}$ and $\theta_{i}$ represent the position and the phase of
the $i$-th swarmalator, respectively. Each individual is affected by intrinsic
spatial and angular velocities $(\vec{\varepsilon}_{i},\omega_{i})$, and by
the coupling with other individuals, determined by attraction
$(I_{\alpha},H_{\alpha})$, repulsion $(I_{\rho})$, and influence $(F,G)$
functions. We use $\vec{p}_{ji}$ and $\theta_{ji}$ as compact representations
of ($\vec{p}_{j}-\vec{p}_{i}$) and ($\vec{\theta}_{j}-\vec{\theta}_{i}$),
respectively. It was shown, in the same work, that computing a two-dimensional
(2D) instance of the model would lead to the emergence of several static and
active states, each with unique features. Afterwards, various modifications of
this 2D model were introduced, analyzing effects produced by an external
stimulus [37], chirality of the particles [38], different interactions [39,
40, 41], among others [42, 43, 44, 45, 46, 47]. One drawback shared by most of
these studies, however, is that despite their numeric findings, analytic
results are very limited given the complexity of the model.

The starting point of the work we present here is the one-dimensional (1D)
swarmalators model, also denominated “ring model”, introduced in [48].
Following the same terminology as in Eqs. (1), the 1D model is defined by

$\displaystyle\begin{split}\dot{x}_{i}&=\frac{1}{N}\sum_{j=1}^{N}I_{\alpha}(x_{ji})F(\theta_{ji}),\\\
\dot{\theta}_{i}&=\frac{1}{N}\sum_{j=1}^{N}H_{\alpha}(\theta_{ji})G(x_{ji}),\end{split}$
(2)

where attraction and influence functions are chosen to be sines and cosines,
respectively, weighted by scalar coupling constants, and
$(x_{i},\theta_{i})\in(\mathbb{S}^{1},\mathbb{S}^{1})$. Individuals are
assumed to be identical, so intrinsic velocities are not considered. Thus, the
dynamics of the system are defined by a couple of Kuramoto-like equations. In
these, the independent synchronization of positions and phases, promoted by
the sines, is strengthened by the interplay induced by the cosines. Given the
simplicity in handling periodic functions, this version of the 1D swarmalators
model is analytically tractable. Taking advantage of this, several solvable
variations have been presented afterwards, with a focus on analyzing the
effects generated by nonidentical swarmalators [49], distributed couplings
[50], and other factors [51, 52, 53]. One drawback of the 1D swarmalators
model in its ‘bare’ form [Eq. (2)] is that its structure prevents the
emergence of active states, where $(\dot{x}_{i},\dot{\theta}_{i}\neq 0)$, in
contrast with the 2D instance of Eq. (1), where active states emerge even for
identical particles. To observe such states in the 1D framework, the solutions
found so far involved adding frustration [52], external stimuli [53],
considering nonidentical individuals [49], or by splitting the population and
mixing coupling signs [51].

Here, we introduce a 2D swarmalators model that shares features with the ones
described in Eqs. (1) and (2), leveraging their advantages and overcoming
their drawbacks: the particles are able to move in a $2\pi$-periodic 2D space,
and the model is tractable analytically. We show that the model presents
several static and active states, some of them similar to ones found in
previous studies and some of them new. The active states, moreover, are part
of a chaotic regime which, under certain conditions, generates a
dimensionality transition.

We outline the general form of our model in Sec. II, along with the specific
equations used in this work. Next, in Sec. III, we present the static and
active states that emerge through numerical computations of the model. Upon
the intuition gained from the numerical results, in Sec. IV, we show the main
results concerning the stability of the static states. In Sec. V, we present
the conditions where the system shows chaotic behaviors and the dimensionality
transition. In Sec. VI, we summarize static and active states constructing
diagrams showing their stable regions. We conclude this work by discussing our
findings and future research suggestions in Sec. VII.

##  II The model

We aim to exploit the solvability of the 1D model introduced in [48], in a 2D
setup. However, before we describe our generalization, we find it useful to
comment on the limitations that arise when trying to recover the dynamics in
Eqs. (1) from Eqs. (2). First, and most evident, the 1D model lacks an
explicit repulsive component. This drawback hinders the emergence of states
that rely on the different scales of attraction and repulsion forces. For
instance, the 1D model cannot generate active states, where the particles keep
moving after a transient. Even if an additional position coordinate $y_{i}$ is
considered, defined to be symmetric with $x_{i}$, the Kuramoto-like feature,
as defined in [48], would lead to the synchronization in each axis, collapsing
the system to static analogs of the 1D states. Moreover, the attraction and
repulsion functions in Eqs. (1) depend on the spatial distances between the
particles $\vec{p}_{ji}$. For example, the 2D instance presented in [22] use a
power law in $p_{ji}$ as repulsion function, generating collective states
where particles distribute radially. Extensions of the ring model to 2D must,
therefore, include explicit dependence of repulsive interactions on particle’s
positions.

For our 2D setup, we define the position of the $i$-th swarmalator as
$\vec{p}_{i}=(x_{i},y_{i})$ but use the distances in each axis, $x_{ij}$ and
$y_{ij}$, independently in the dynamical equations, as using the modulus
$p_{ij}$ would undermine the reduction of complexity that we are looking for.
The coupling between $x$ and $y$ is done through the phases $\theta_{i}$,
included in the functions $F_{\alpha}$ and $F_{\rho}$ that affect attraction
$I_{\alpha}$ and repulsion $I_{\rho}$, respectively. The phase dynamics, on
the other hand, are defined by their mutual interaction $H$ and by the
interplay of both spatial coordinates $G$. In its general form the model is
described by

$\displaystyle\begin{split}\dot{x}_{i}&=u_{i}+\frac{1}{N}\sum_{j=1}^{N}\left[I_{\alpha}(x_{ji})F_{\alpha}(\theta_{ji})-I_{\rho}(x_{ji})F_{\rho}(\theta_{ji})\right],\\\
\dot{y}_{i}&=\upsilon_{i}+\frac{1}{N}\sum_{j=1}^{N}\left[I_{\alpha}(y_{ji})F_{\alpha}(\theta_{ji})-I_{\rho}(y_{ji})F_{\rho}(\theta_{ji})\right],\\\
\dot{\theta}_{i}&=\omega_{i}+\frac{1}{N}\sum_{j=1}^{N}H_{\alpha}(\theta_{ji})G(x_{ji},y_{ji}).\end{split}$

For the specific model we study here, we consider that
$(x_{i},y_{i},\theta_{i})\in(\mathbb{S}^{1},\mathbb{S}^{1},\mathbb{S}^{1})$,
and individuals are assumed to be identical, so
$(u_{i},\upsilon_{i},\omega_{i})=0$. Following [48], we choose combinations of
sines and cosines for all functions, weighted by coupling constants $J^{\pm}$
and $K^{\prime}$. This makes the problem amenable to analytical treatment, but
introduces repulsion in a weak sense:

$\displaystyle\begin{split}\dot{x}_{i}&=\frac{1}{N}\sum_{j=1}^{N}\left\\{J^{-}\sin(x_{ji})\cos(\theta_{ji})-J^{+}\left[1-\cos(x_{ji})\right]\sin(\theta_{ji})\right\\},\\\
\dot{y}_{i}&=\frac{1}{N}\sum_{j=1}^{N}\left\\{J^{-}\sin(y_{ji})\cos(\theta_{ji})-J^{+}\left[1-\cos(y_{ji})\right]\sin(\theta_{ji})\right\\},\\\
\dot{\theta}_{i}&=\frac{K^{\prime}}{N}\sum_{j=1}^{N}\sin(\theta_{ji})\left[2+\cos(2x_{ji})+\cos(2y_{ji})\right].\end{split}$
(3)

In this form, $J^{-}$ is the weight of the Kurmamoto-like attractive
interaction, that is enhanced by phase synchronization. $J^{+}$, on the other
hand, is the weight of the weak repulsion term, whose sign depends on the
difference of phases. It is weak because it disappears when particles are on
top of each other, allowing for full synchronization in $x$, $y$ and $\theta$.
When phases are synchronized, $\theta_{ij}=0$, the equations for $x_{i}$ and
$y_{i}$ reduce to independent Kuramoto dynamics. Similarly, when
$x_{ij}=y_{ij}=0$, the equation for $\theta_{i}$ follows the Kuramoto model.
Repulsion and attraction terms in each direction depend only on the distances
in that direction. This is the key feature that allows for simplifications in
the analytical treatment while still supporting active states.

As a final step, for an easier understanding of the system, we rewrite
$J^{\pm}~{}=~{}(J_{A}\pm J_{R})/2$ and $K^{\prime}=K/2$. These new parameters
do not alter the structure shown in Eqs. (3), but facilitate manipulating
attraction and repulsion scales, affected respectively by $J_{A}$ and $J_{R}$.
Thus, the expanded model is given by

$\displaystyle\begin{split}\dot{x}_{i}&=\frac{1}{N}\sum_{j=1}^{N}\frac{J_{A}}{2}\left\\{\sin(x_{ji})\cos(\theta_{ji})-\left[1-\cos(x_{ji})\right]\sin(\theta_{ji})\right\\}-\frac{J_{R}}{2}\left\\{\sin(x_{ji})\cos(\theta_{ji})+\left[1-\cos(x_{ji})\right]\sin(\theta_{ji})\right\\},\\\
\dot{y}_{i}&=\frac{1}{N}\sum_{j=1}^{N}\frac{J_{A}}{2}\left\\{\sin(y_{ji})\cos(\theta_{ji})-\left[1-\cos(y_{ji})\right]\sin(\theta_{ji})\right\\}-\frac{J_{R}}{2}\left\\{\sin(y_{ji})\cos(\theta_{ji})+\left[1-\cos(y_{ji})\right]\sin(\theta_{ji})\right\\},\\\
\dot{\theta}_{i}&=\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_{ji})\left\\{1+\frac{1}{2}\left[\cos(2x_{ji})+\cos(2y_{ji})\right]\right\\}.\end{split}$
(4)

Notice that, in this form, the terms proportional to $J_{R}$ in the equation
for $\dot{x}_{i}$ include a truly repulsive interaction, $-\sin(x_{ji})$, and
half of the weakly repulsive term, $-[1-\cos(x_{ji})]$ (and similarly for
$\dot{y}_{i}$). The terms in $J_{A}$, however, include the attractive part
$\sin(x_{ji})$ but also the other half of the weakly repulsive term. By
considering this expanded 2D model, we can derive the original ring model
under intuitive convenient conditions (see Appendix A).

##  III Numerical simulations

Simulations performed for the model described in Eqs. (4) show the emergence
of states that gather features of several states observed in instances of Eqs.
(1) and (2). Additionally, we observe active states that, despite the
similarity with states presented in previous works, show interesting
properties. All these states are presented in Figs. 1, 2, and 3, where the
computations were performed for a population of $N=500$ swarmalators and
different values of $J_{A}$, $J_{R}$, and $K$. The particles start positioned
uniformly in a $2\pi$-length cube in the $(x,y,\theta)$ space, and their
states evolve along $10^{5}$ time-steps. From Fig. 1 to Fig. 3, each shows the
particles distributed in the $(x,y)$ space colored according to their phases
(top rows), and their respective distributions in the $(x,y,\theta)$ space
(bottom rows).

In Figs. 1(a) and 1(d), we present the point synchronous state, which shows
the phase synchronization of particles and their convergence to a fixed point
in the $(x,y)$ space. Once these reach a $(x_{0},y_{0},\theta_{0})$ state,
they remain steady. A similar state emerges on the 1D model [48], where
particles converge to a fixed position on the ring. In our model, the
emergence of this state requires phase synchronization $(K>0)$, and also
$(J_{A},J_{R}>0)$. Moreover, attraction must be stronger than repulsion,
$J_{A}>J_{R}$, so that particles, which start distributed across the $(x,y)$
space, can collapse to a fixed position. For $J_{R}>J_{A}$ the particles
converge to the distributed synchronous state shown in Figs. 1(b) and 1(e). In
this state, particles synchronize (as $K>0$) but they remain distributed in
the $(x,y)$ space given that now repulsion is stronger than attraction. This
state shows a first approach to spatial states generated by instances of Eqs.
(1), where particles tend to distribute radially [22, 39, 37]. Finally, in
Figs. 1(c) and 1(f), we show the distributed asynchronous state. In this, the
individuals keep their initial phases frozen $(K=0)$, and the repulsion
between them is strong enough to keep the particles distributed uniformly in
the $(x,y)$ space. Similar states can also be observed emerging in the 1D
model [48, 50], where particles remain asynchronized while distributed along
the ring, and in 2D and 3D instances of Eqs. (1), where the spatial
distribution is radial [22, 39].

Figure 1: Snapshots of the [(a),(d)] point synchronous, [(b),(e)] distributed
synchronous, and [(c),(f)] distributed asynchronous states after a transient.
The panels show the spatial distribution of the particles (top row) and the
respective correlation scatter plots (bottom row). The parameters
$(K,J_{A},J_{R})$ are set as $(1,1,0.5)$ for the state in [(a), (d)],
$(1,0.5,1)$ for the state in [(b), (e)], and $(0,-0.5,0.5)$ for the state in
[(c),(f)]. See movies S1, S2, and S3 in the Supplemental Material.

Panels in Fig. 2 show three types of static phase wave state. First, Figs.
2(a) and 2(d) show that individuals’ positions $(x_{i},y_{i})$ are positively
correlated and that phases $\theta_{i}$ are also positively correlated to
these positions. Second, Figs. 2(b) and 2(e) show that individuals’ positions
$(x_{i},y_{i})$ are positively correlated but that phases $\theta_{i}$ are
negatively correlated to $x_{i}$ and therefore $y_{i}$. Third, Figs. 2(c) and
2(f) show that individuals’ positions $(x_{i},y_{i})$ are negatively
correlated and that phases $\theta_{i}$ are positively correlated to $x_{i}$
only. Each of these states shows similarities with states found in previous
work. The static phase wave state in 1D shows the emergence of correlations
between particles’ positions and phases [48], and in 2D, the correlation
emerges between the polar angle that describes the radial position of each
particle and its respective phase [22]. A remarkable feature in our model,
furthermore, is that independent of the type of correlation, it features the
concept of ‘like attracts like’ introduced in [22]: despite individuals’
phases being frozen $(K=0)$, they end up grouping in space with similar phases
individuals. For simplicity, in the following sections, we will refer to each
type of static phase wave state as summarized in Table 1. An additional
observation is that, Figs. 2(c) and 2(f) show a small curvature in the
correlation, that reminds us of the buckled phase wave state introduced in
[48].

| Type of correlation
---|---
State | Positive | Negative
Static phase wave I | $x_{i},y_{i},\theta_{i}$ |
Static phase wave II | $x_{i},y_{i}$ | $x_{i},\theta_{i}$
Static phase wave III | $x_{i},\theta_{i}$ | $x_{i},y_{i}$
Table 1: Classification of the static phase wave states according to the type
of correlation between individuals’ positions and phases. Figure 2: Snapshots
of the static phase wave states after a transient. The panels show the spatial
distribution of the particles (top row) and the respective correlation scatter
plots (bottom row).[(a), (d)], [(b), (e)], and [(e), (f)] pairs correspond to
different types of correlation between $x_{i},y_{i}$, and $\theta_{i}$ states.
The parameters $(K,J_{A},J_{R})$ are set as $(0,0,-0.8)$ for the state in
[(a), (d)], $(0,0.8,0)$ for the state in [(b), (e)], and $(0,1,-1)$ for the
state in [(c),(f)]. See movies S4, S5, and S6 in the Supplemental Material.

In Fig. 3, we show instances of the most noteworthy states of this work: the
active ones. We call these ‘instances’, since all of them are part of a
chaotic regime that we explore in more detail in the following sections. The
first one, presented in Figs. 3(a) and 3(d), shows that after a transient,
there is an apparent correlation between $x_{i}$, $y_{i}$, and $\theta_{i}$,
similar to the ones presented for the static phase wave states in Fig 2.
However, individuals move, and despite the positive correlation between
$x_{i}$ and $y_{i}$, shown in Fig. 3(a), individuals’ positions and phases
generate unsteady irregular annular shapes, as noted in Fig. 3(d). These
swirls remind us of the ones formed by non-conformist individuals in [51], so
we name our state after it. We refer to the second instance of active state,
shown in Figures 3(b) and 3(e), as butterfly state. In this state, particles
move while positions $x_{i}$ and $y_{i}$ keep linearly correlated.
Furthermore, the population of particles continuously alternates between
splitting into two clusters and merging back to a single one. Similar to the
swirling state, individuals’ positions and phases are not correlated.
Furthermore, the shape generated by these, as shown in Fig. 3(e), reminds us
of the well-known Lorenz attractor projected in 2D. Finally, the bouncing
state, depicted in Figs. 3(c) and 3(f), shows similarities with the butterfly
state in terms of the correlation between individuals’ positions [Fig. 3(c)]
and the shapes generated between their positions and phases [Fig. 3(f)].
However, the population of particles goes back and forth, from a single
cluster to several ones, in an erratic way. Despite this behavior, the linear
correlation between $x_{i}$ and $y_{i}$ is preserved.

Figure 3: Snapshots of the [(a),(d)] swirling, [(b),(e)] butterfly, and
[(c),(f)] bouncing states after $10^{5}$ time-steps. The panels show the
spatial distribution of the particles (top row) and the respective correlation
scatter plots (bottom row). The arrows represent the motion of the particles.
The parameters $(K,J_{A},J_{R})$ are set as $(-1,0,-0.5)$ for the state in
[(a), (d)], $(-1,-0.5,-0.5)$ for the state in [(b), (e)], and $(-1,-1,-0.5)$
for the state in [(c),(f)]. See movies S7, S8, and S9 in the Supplemental
Material.

##  IV Static states

We perform perturbation analyses to study the stability of the static
equilibrium states reached by our system. The computations we carry out are
based on the analyses presented in [52] for an instance of the 1D model. Then,
since the procedures are known, only major considerations and results are
presented in the main text and complex calculations are left for the Appendix.

###  IV.1 Point synchronous (PS) state

This state is characterized by the convergence of the particles to a point
state $x_{i}=x$, $y_{i}=y$, and $\theta_{i}=\theta$, which proves to be a
solution of Eqs. (3). We add small perturbations $\delta x_{i}$, $\delta
y_{i}$, and $\delta\theta_{i}$ to each particle in their equilibrium state,
and find expressions for their time evolution. The perturbation dynamics and
their analysis as a linear system are shown in detail in Appendix B.

The eigenvalues that determine the stability of this state are

$\displaystyle\begin{split}\lambda^{PS}_{0}&=0,\\\ \lambda^{PS}_{1}&=-2K,\\\
\lambda^{PS}_{2}&=\frac{J_{R}-J_{A}}{2},\end{split}$ (5)

revealing that the PS state emerges only when $K>0$ and $J_{A}\geq J_{R}$, as
expected from the remarks in the previous section: phase synchronization is
driven by $K>0$ and particles’ clustering result from $J_{A}>J_{R}$. Moreover,
from Eqs. (3), we can see that phase synchrony strengthen the attractive
effects in spatial dynamics, reason why particles converge to a single point
in space instead of cluttering while keeping distance from each other.

###  IV.2 Static phase wave (SPW) states

The main feature shared by these states is the linear correlation among
$x_{i}$, $y_{i}$, and $\theta_{i}$, as classified in Tab. 1. It is easy to
check that equilibrium states based on these correlations satisfy the
dynamical equations.

In the static phase wave I, the equilibrium state is defined by $x_{i}=2\pi
i/N+x_{0}$, $y_{i}=2\pi i/N+y_{0}$, and $\theta_{i}=2\pi i/N+\theta_{0}$. The
perturbation analysis performed for this equilibrium is described in Appendix
B.2. The eigenvalues that determine the stability of this state are

$\displaystyle\begin{split}\lambda^{(I)}_{0}&=0,\\\
\lambda^{(I)}_{1}&=0.5J_{R},\\\ \lambda^{(I)}_{2}&=0.5J_{R}+0.25J_{A},\\\
\lambda^{(I)}_{3}&=0.25J_{R}+0.125K\\\ &\qquad\pm
0.5(0.25J_{R}^{2}+0.0625K^{2}-1.25J_{R}K)^{0.5},\\\
\lambda^{(I)}_{4}&=0.25J_{R}+0.375K\\\ &\qquad\pm
0.5(0.25J_{R}^{2}+0.562K^{2}+0.75J_{R}K\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad+0.5J_{A}K)^{0.5}.\\\ \end{split}$
(6)

The simplest case to frame, using these eigenvalues, is when $(K=0)$, which
implies that individual phases will remain distributed randomly as in their
initial state. Then, the conditions for stability are $J_{R}\leq 0$ and
$J_{R}\leq-0.5J_{A}$. The negative repulsion and the positive attraction terms
indicate a primarily attractive nature of the particles. However, because
phases are different, particles cluster with others of similar phases.
Additionally, considering $(K>0)$ would lead to phase synchronization and,
consequently, disrupt the clustering behavior among different phases. This
behavior is also aligned with our expectation for this state based on the
nature of our model. The case where $(K<0)$ is more complex given that the
phase asynchrony is now weighted, incorporating an active ingredient to the
system. We explore this active state in the following section.

The static phase wave II is given by $x_{i}=2\pi i/N+x_{0}$, $y_{i}=2\pi
i/N+y_{0}$, and $\theta_{i}=-2\pi i/N+\theta_{0}$. The details of the
perturbation analysis for this equilibrium state are shown in Appendix B.3.
The eigenvalues that determine its stability are

$\displaystyle\begin{split}\lambda^{(II)}_{0}&=0,\\\
\lambda^{(II)}_{1}&=-0.5J_{A},\\\ \lambda^{(II)}_{2}&=-0.5J_{A}-0.25J_{R},\\\
\lambda^{(II)}_{3}&=-0.25J_{A}+0.125K\\\ &\qquad\pm
0.5(0.25J_{A}^{2}+0.0625K^{2}+1.25J_{A}K)^{0.5},\\\
\lambda^{(II)}_{4}&=-0.25J_{A}+0.375K\\\ &\qquad\pm
0.5(0.25J_{A}^{2}+0.562K^{2}-0.75J_{A}K\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad-0.5J_{R}K)^{0.5}.\\\ \end{split}$
(7)

The conditions for stability are similar to those found in the previous case
(static phase wave I): for $K=0$, we need $J_{A}\geq 0$ and $-J_{A}\geq
0.5J_{R}$. Both conditions suggest the attractive nature of the particles,
which will lead to the clustering of particles with similar phases. When
$K>0$, diversity among phases is disrupted, whereas $K<0$ drives the emergence
of active states. The main difference between this state and the previous one
lies in the regions of the $(J_{A}-J_{R})$ plane where each state appears, as
we will demonstrate in a following section.

Finally, the equilibrium states in the static phase wave III is defined as
$x_{i}=2\pi i/N+x_{0}$, $y_{i}=-2\pi i/N+y_{0}$, and $\theta_{i}=2\pi
i/N+\theta_{0}$. Despite the similarity of these expressions to the equilibria
in the two previous cases, the analysis is more complex. The negative
correlation between $x_{i}$ and $y_{i}$ prevent simplifications in the
perturbation analysis. Consequently, finding the eigenvalues that characterize
the state’s stability becomes more intricate. For details of this analysis, we
refer the reader to the Appendix B.4. Although we do not present the
corresponding eigenvalues in the main text, in a following section we will
show the stability regions of this state in the $(J_{A}-J_{R})$ plane.

###  IV.3 Distributed states

To study the stability of these states we assume a continuum of particles,
instead of a discrete set like in the previous cases. We define the fraction
of particles lying between $x$, $y$, $\theta$ and $x+\mathrm{d}x$,
$y+\mathrm{d}y$, $\theta+\mathrm{d}\theta$ at time $t$, by the density
$\rho(x,y,\theta,t)\mathrm{d}x\mathrm{d}y\mathrm{d}\theta$. Moreover, $\rho$
will satisfy the normalization condition

$\int\rho(x,y,\theta)\mathrm{d}x\mathrm{d}y\mathrm{d}\theta=1$ (8)

where all integrals run from 0 to $2\pi$.

To perform the stability analysis we introduce order parameters of the form

$S_{\sigma}e^{\mathrm{i}\mkern
1.0mu\phi_{\sigma}}=\frac{1}{N}\sum_{j=1}^{N}e^{\mathrm{i}\mkern
1.0mu\sigma_{j}}.$ (9)

In this expression, $\sigma$ can represent any individual variable $x$, $y$,
$\theta$, or a linear combination of these. So, for instance, if we consider
$\sigma\mathrel{\mathop{\mathchar 58\relax}}\theta$, the order parameter will
measure phase coherence, giving $S_{\theta}=1$ for a fully synchronized state.
Another instance, a bit more interesting is the consideration of
$\sigma\mathrel{\mathop{\mathchar 58\relax}}(\theta+x)$, which measures the
correlation between individuals’ positions $x_{j}$ and phases $\theta_{j}$. In
this case, getting $S_{(\theta+x)}=1$ is an indicator of a negative
correlation between $x_{i}$ and $\theta_{i}$ as in the static phase wave state
II.

In the continuum limit, considering Eq. (8), we can rewrite the order
parameters as

$S_{\sigma}e^{\mathrm{i}\mkern
1.0mu\phi_{\sigma}}=\int_{0}^{2\pi}e^{\mathrm{i}\mkern
1.0mu\sigma}\rho(x,y,\theta)\mathrm{d}x\mathrm{d}y\mathrm{d}\theta,$ (10)

so that Eq. (3) turns into the mean-field equations

$\displaystyle\begin{split}v_{x}&=\frac{J_{A}}{2}S_{(\theta+x)}\sin\left(\phi_{(\theta+x)}-\theta-x\right)+\frac{J_{R}}{2}S_{(\theta-x)}\sin\left(\phi_{(\theta-x)}-\theta+x\right)-\frac{J_{A}+J_{R}}{2}S_{\theta}\sin\left(\phi_{\theta}-\theta\right),\\\
v_{y}&=\frac{J_{A}}{2}S_{(\theta+y)}\sin\left(\phi_{(\theta+y)}-\theta-y\right)+\frac{J_{R}}{2}S_{(\theta-y)}\sin\left(\phi_{(\theta-y)}-\theta+y\right)-\frac{J_{A}+J_{R}}{2}S_{\theta}\sin\left(\phi_{\theta}-\theta\right),\\\
v_{\theta}&=KS_{\theta}\sin\left(\phi_{\theta}-\theta\right)+\frac{K}{4}\Big{[}S_{(\theta+2x)}\sin\left(\phi_{(\theta+2x)}-\theta-2x\right)+S_{(\theta-2x)}\sin\left(\phi_{(\theta-2x)}-\theta+2x\right)\\\
&\qquad\qquad\qquad\qquad\qquad\qquad+S_{(\theta+2y)}\sin\left(\phi_{(\theta+2y)}-\theta-2y\right)+S_{(\theta-2y)}\sin\left(\phi_{(\theta-2y)}-\theta+2y\right)\Big{]}.\end{split}$
(11)

A fully incoherent state, as depicted in the distributed asynchronous case,
would be portrayed by the convergence of all the order parameters, in Eqs.
(11), to zero. However, we can also describe the distributed synchronous state
by considering that incoherence happens in space only, letting the system
reach phase synchronization ($S_{\theta}=1$). For both cases, we define the
respective incoherent densities

$\displaystyle\begin{split}\rho_{0}^{A}&=\frac{1}{8\pi^{3}},\\\
\rho_{0}^{S}&=\frac{1}{4\pi^{2}},\end{split}$ (12)

which under conditions of the order parameters described above, prove to be
equilibrium states of the continuity equation

$\frac{\partial\rho}{\partial t}=-\vec{\nabla}(\rho\vec{v}),$ (13)

where $\vec{v}=(v_{x},v_{y},v_{\theta})$ is the velocity field defined in Eq.
(11). We perform perturbation analyses for both incoherent states (see
Appendix C). The eigenvalues that determine the stability of the distributed
asynchronous state are

$\displaystyle\begin{split}\lambda_{1}^{A}&=\frac{J_{A}}{16\pi^{2}},\\\
\lambda_{2}^{A}&=-\frac{J_{R}}{16\pi^{2}},\\\
\lambda_{3}^{A}&=\frac{K}{8\pi^{2}},\\\
\lambda_{4}^{A}&=\frac{K}{32\pi^{2}},\end{split}$ (14)

and, for the distributed synchronous state,

$\lambda^{S}=\frac{J_{A}-J_{R}}{8\pi^{2}},$ (15)

which additionally requires $(K>0)$ to reach phase coherence.

##  V Active states

In the previous section, we discussed all the conditions that allow the
emergence of the static states. Despite these analytical results, we found
only very few hints on the conditions driving the emergence of active states.
In order to get more clues on their behavior, we start our analysis by
considering distributed states with finite number of particles. Then, we will
demonstrate that all the active states shown in Sec. III are actually part of
the same chaotic regime.

###  V.1 Active asynchronous state

The analysis of distributed states performed in Sec. IV.3 relied on the
assumption that the number of particles is infinite, so that we could take the
continuum limit. This allowed us to consider that all order parameters in Eq.
(11) would converge asymptotically to zero. However, when dealing with a
finite number of particles, this last assumption is not quite true. For $K=0$,
in particular (frozen phases), phase incoherence is not perfect and
$S_{\theta}>0$. This drawback promotes the appearance of small velocities in
$\dot{x}_{i}$ and $\dot{y}_{i}$, proportional to $(J_{A}+J_{R})/2$, which
ultimately will drive the emergence of an active state. For an instance of
this state, see movie S10 in the Supplemental Material. If $J_{A}=-J_{R}$,
i.e., attraction and repulsion have the same intensity, the state still
converges to a static configuration.

###  V.2 Chaos in 3N dimensions

In Sec. IV we studied the stability of static states by computing the
eigenvalues of the linearized dynamics in scenarios where phases were frozen,
$K=0$, or driven towards synchronization by $K>0$. Negative phase couplings
($K<0$), however, not only makes computations more complex but also drives the
emergence of active states (as shown in Figs. 3).

An interesting feature of these states is that particles are sensitive to
small changes in their initial conditions. We highlight this property in Figs.
4(a), 4(b), and 4(c), where trajectories in the $x$-axis are shown for the
same particle in scenarios where $K=-1$. In each of these, we compare the
trajectories followed by a single particle when the system starts on the
$N$-dimensional initial states $(\vec{x}_{0},\vec{y}_{0},\vec{\theta}_{0})$
and $(\vec{x}_{0}+\delta^{(1)}_{x},\vec{y}_{0},\vec{\theta}_{0})$, where
$\delta^{(1)}_{x}=10^{-4}$ represents that only the particle of interest’s
state is perturbed (for instance $x_{1(0)}+\delta^{(1)}_{x}$). It is clear
then that, when considering $(J_{A},J_{R})=(0.5,-0.5)$ [Fig. 4(a)], the
perturbation does not affect the trajectory of the particle considerably. In
fact, from Eqs. (6), we can infer that these parameters drive the emergence of
the static phase wave I. However, once $J_{A}$ becomes negative
[$(J_{A},J_{R})=(-0.5,-0.5)$ for Fig. 4(b) and $(J_{A},J_{R})=(-0.5,-3)$ for
Fig. 4(c)], the perturbed trajectories suggest the existence of chaos.

In order to explore the presence of chaotic states we compute the maximum
Lyapunov exponent as a function of $J_{A}$ for $K=-1$ and $J_{R}=-0.5$. We
take two $3N$-dimensional trajectories of the system, one with initial
conditions $(\vec{x}_{0},\vec{y}_{0},\vec{\theta}_{0})$, and the other one
with $(\vec{x}_{0}+\delta^{(1)}_{x},\vec{y}_{0},\vec{\theta}_{0})$. We remark
that, although $\delta^{(1)}_{x}$ affects the initial condition of a single
particle in the $x$ direction, it perturbs the entire $3N$-dimensional state
through the couplings. We then calculate the evolution of the distances
$d_{(t)}=d_{0}e^{\lambda t}$ between the $3N$-dimensional trajectories and
infer the Lyapunov exponents $\lambda$ from the transitioning slope of
$\log[d_{(t)}/d_{0}]$. These are shown in Fig. 4(d) and, as expected, there is
a threshold, at about $(J_{A}=0.2)$, where the system jumps from regular to
chaotic (positive Lyapunov exponent). Additionally, from Eqs. (6), we can see
that $J_{A}\in[0,0.2]$, is a region where the static phase wave I state is
also stable, meaning that it is a region of bistability between the static
phase wave I and the chaotic regime.

State trajectories in the chaotic regime form strange attractors whose
structures differ according to the parameters $(K,J_{A},J_{R})$. Fig. 5 shows
a bidimensional projection of the chaotic trajectory followed by a single
particle in the $x-\theta$ torus in four scenarios. We can see that, depending
on the scale difference between attracting $(J_{A})$ and repulsing $(J_{R})$
terms, the trajectories patterns grow from a small disc [Fig. 5(a)], to
butterfly-like [Fig. 5(b)], to a disordered single-loop [Fig. 5(c)], and end-
up with a scribble-like structure that covers the whole torus [Fig. 5(d)].

Another interesting phenomenon happens when the scale difference between
$J_{A}$ and $J_{R}$ surpass a threshold, leading to a transition of
dimensionality. We can see in the insets of Fig. 4(d) that, despite being
inside the chaotic regime, when $(J_{A}>-2)$, the linear correlation between
$x_{i}$ and $y_{i}$ is held. However, once $J_{A}$ gets below this threshold,
this correlation breaks. This behavior is strikingly unexpected given that
expressions for $\dot{x}_{i}$ and $\dot{y}_{i}$ are symmetric, and even the
variations in $(J_{A},J_{R})$ are the same. This dimensionality jump partially
explains the scribble-like structure in Fig. 5(d), given that more dimensions
might be required to show a clearer attractor structure.

Figure 4: Temporal behavior of single particle trajectories along the $x$-axis
with non-perturbed (lavender) and perturbed (red) initial conditions (top
row), and maximum Lyapunov exponents for different values of $J_{A}$ (bottom
row). Parameters $(K,J_{R})$ are set as $(-1,-0.5)$ for all the figures.
Trajectories are generated for (a) $J_{A}=0.5$, (b) $J_{A}=-0.5$, and (c)
$J_{A}=-3$. The square in the top corner of (a) shows a magnification of the
circled trajectories, and the ones in (d) show correlations in $x-y$
coordinates depending on the values of $J_{A}$ (See movies S11 and S12 in the
Supplemental Material). The green dashed lines represent the change from
static phase wave to chaotic states (left) and the breaking of $x-y$
coordinates correlation (right).

###  V.3 Low-dimensional chaos

The chaotic behavior of the system can also be identified through the
correlation functions defined in Sec IV.3. Specifically, we use the metrics

$\displaystyle S_{(\theta\pm x)}^{max}$
$\displaystyle=\max\left[S_{(\theta+x)},S_{(\theta-x)}\right],$ $\displaystyle
S_{(x\pm y)}^{max}$ $\displaystyle=\max\left[S_{(x+y)},S_{(x-y)}\right],$

so we can get a measure of the linear correlation between the respective
variables, at any specific time-step, independent on their sign. As shown in
the bottom row of Fig. 5, perturbing the initial conditions, generates order
parameters that evolve differently. Some of these [Figs. 5(e), (f), (g)]
fluctuate around a fixed value, which indicates that we can use them to
discern between states. In Fig. 5(h), however, the fluctuations have a much
larger amplitude, hindering any attempt to classify the state based on its
value. Additionally, the dimensionality jump is also spotted by $S^{max}_{x\pm
y}$, given that it converges to $1$ when $J_{A}$ is over the threshold [Figs.
5(e),(f), and (g)], indicating that $x_{i}$ and $y_{i}$ are fully correlated,
but it fluctuates once $J_{A}$ falls below it [Fig. 5(h)].

Figure 5: Top row: trajectories (blue) followed by single particles in the
$x-\theta$ torus. The red circles are positioned at the start of each
trajectory. Bottom row: evolution of $S^{max}_{(\theta\pm x)}$ and
$S^{max}_{(x\pm y)}$ for the original initial conditions (purple and gold,
respectively) and the perturbed initial conditions (orange and green,
respectively). Parameters $(K,J_{A},J_{R})$ are set as [(a), (e)]
$(-1,0,-0.5)$, [(b), (f)] $(-1,-0.5,-0.5)$, [(c), (g)] $(-1,-1,-0.5)$, and
[(d), (h)] $(-1,-3,-0.5)$. The small rectangle in (a) shows a magnification of
the circled lines. See movies S13, S14, S15, and S16 in the Supplemental
Material (for trajectories of two particles see movies S17, S18, S19, and
S20).

##  VI Attraction-Repulsion phase diagram

In Fig. 6, we present a summary of the regions where each state emerges in the
$J_{A}-J_{R}$ plane. These are based on the eigenvalues obtained in Sec. IV
and the Lyapunov exponents calculated in Sec. V.

Note that the active states analyzed in previous sections assumed a positive
correlation between $x_{i}$, $y_{i}$, and $\theta_{i}$, corresponding to an
extension of static phase wave I for $(K<0)$. We refer to this chaotic regime
as Chaos I. Similarly, active effects remain consistent when considering the
correlation of static phase wave II, leading us to name that region Chaos II.
We can also spot several regions of bistability and multistability. Not only
between static states, but also between static and active ones. Remarkably,
the regions shown in Fig. 6, consider that initial conditions are distributed
randomly in $x$, $y$, and $\theta$. Starting in full synchronization
$(\theta_{i}=\theta)$ and considering frozen phases $(K=0)$, would lead to the
emergence of a synchronous state depending on $(J_{A},J_{R})$. More exotic
behaviors, that we do not study here, can be obtained when considering non-
random initial conditions.

Figure 6: Regions where states emerge on the $J_{A}-J_{R}$ plane for (a)
$K<0$, (b) $K=0$, and (c) $K>0$. Arrows and gradients in (a) represent state
variations in the chaotic regime, transitioning from swirling to bouncing
states. Colors labeling different states are an indicator of multistability.

##  VII Discussion

We have studied a swarmalators system where particles move in the 2D plane
with periodic boundary conditions. The model, that generalizes the 1D instance
introduced in [48], considers attraction and repulsion terms, which ultimately
lead the emergence of distributed and active states, some of them similar to
the ones introduced in [22] and [48]. All interaction functions involve only
sines and cosines, making stability analysis relatively simple, but leading to
a weak form of repulsion. In contrast with previous 2D models, we obtained
analytical conditions for the emergence of all the static states we presented
and also analyzed the nature of the active ones.

Two states presented in this work, the distributed synchronous and
asynchronous states, closely resemble their counterparts described in [22],
the only difference being that, in our findings, particles are not radially
distributed. For these, we determined analytically the conditions that drive
their emergence considering only control parameters $(K,J_{A},J_{R})$. Along
this line, we have also introduced an active asynchronous state, that emerges
under specific conditions for a finite number of particles. Furthermore, we
derived the conditions that guarantee the emergence of a point synchronous
state, similar to the one introduced in [48] for the 1D model.

We also described static phase wave states that share features with their 1D
[48, 50] and 2D counterparts [22]. The main difference with the former is
that, in our model, linear correlations show up in three different ways, as
classified in Tab. 1. Moreover, as pointed before for the distributed
synchronous and asynchronous states, our model leads to square, instead of
radial, symmetry in particle distribution. Despite these differences, the
phase wave states share the ‘like attracts like’ feature, characteristic of
their definition. Our analytical findings allowed us to state all the
conditions in $(K,J_{A},J_{R})$ that drive the emergence of each type of
correlation. Strikingly, we found that for $(K\leq 0)$ there is a region of
multistability in the $J_{A}-J_{R}$ plane, where the three types of phase wave
state could emerge.

Despite the static states and their counterparts found in previous works, the
primary distinction between our model and these lies in the active states,
which have remained elusive in 1D simplifications of the original swarmalators
model. Previous studies where active states emerged in 1D were based on the
inclusion of additional parameters, such as frustration [37], external forcing
and pinning [53], or mixed coupling signs [51]. Here we demonstrated that
active states can emerge solely by adjusting the model parameters
$(K,J_{A},J_{R})$ when repulsion is included. This is coherent with the
original idea of swarmalator systems where unsteady states emerge as a result
of scale differences between attractive and repulsive effects. Except for the
radial symmetry, the active states we found are similar to the splintered and
active phase wave states from [22]. However, we showed that the active states
in our model are part of a spectrum of chaotic states. Even more interesting,
we showed that under certain conditions in the chaotic regime, the system
suffers a dimensionality transition, where the spatial correlation breaks and
particles behave chaotically in $x$, $y$, and $\theta$ axes.

From a theoretical perspective, our model incorporates intriguing
characteristics that may pique curiosity as a dynamical system in its own
right. Given the system’s dimensionality and the somewhat simplistic method we
employed to calculate the maximum Lyapunov exponent, it would be worthwhile to
explore the entire spectrum of Lyapunov exponents. We suspect that the chaotic
regime might exhibit hyperchaotic behaviors under certain conditions, where
more than one Lyapunov exponent is positive. It may be also of interest
understanding the bifurcations triggered by varying $K$. In Fig 6, we have
shown three specific cases that allowed us to discern between the states when
$(K<0)$, $(K=0)$, and $(K>0)$. However, it may be worthwhile to study in
detail the intermediate states of the system as $K$ transitions from being
positive-valued to negative-valued and vice versa. This investigation should
provide valuable insights into two open questions: why certain states, like
the swirling one, have a very long transient before exhibiting their active
characteristics, and what mechanisms underlie the emergence of chaos.
Regarding the latter question, we also recommend a rigorous study of the
routes to chaos and dimensionality breaking when $(K<0)$. Preliminary
numerical studies on this subject have given us hints of emergent
intermittency and crises when examining the evolution of the order parameters
for different $(J_{A},J_{R})$. Finally, it would also be worthwhile to explore
whether a threshold on the number of particles $N$ exists, and if so, at what
point chaotic behaviors are triggered.

Our model merges two lines of research on the same topic: 1D systems based on
the ring model and 2D modifications of the original model. Therefore, its
adaptation to previous work on these two lines can occur smoothly. For
instance, future work could take our model as a base to consider chirality
[38], different phase interactions [41, 39, 40], finite cutoff ranges [54],
clusters with different attractive and repulsive strengths [47], or short
range repulsive interactions [55]. Then, given the structure of our model,
finding analytical solutions on each of these modifications would be less
challenging than doing it for instances of the 2D model introduced in [22].

Given that the interactions in our model are based on Kuramoto-like terms,
modifications performed in the Kuramoto model could also be adapted to ours
[56]. A particularly interesting take would be to consider frustration, given
that the addition of lag parameters are responsible for the emergence of a
turbulent-like state in the ring model [52], and also drive the formation of
interesting spatial patterns in swarming only systems [57, 58]. Furthermore,
the inclusion of external stimulus discussed in previous work, implies
affecting the dynamics in a periodic nature [53] and mainly in the phase
dynamics [37, 59]. One suggestion would be to consider external stimuli in
space, mimicking a shear, as has been done for bacterial suspensions [60]. It
would also be of interest considering finite cutoff interaction distances
between individuals [54] so that different concentrations of individuals can
be managed, emulating the study on vinegar eels in [61]. If a collective of
swarmalators works as a medium, it would also be worth studying how it is
affected by external particles and vice versa, similar to the study in [33]
for nematodes. Additionally, we could explore changing the mobility space of
the swarmalators, for instance, by imposing discontinuities in the $(x,y)$
space, or also considering self-propulsion of the particles in the direction
of their phases, as presented in [62].

States emerging in our model could also fit collective swarming-only behaviors
if we consider that their three spatial degrees of freedom interplay. A clear
example of this approach is depicted by the behavior of swarming mosquitoes
[63, 64], which shows similarities with the active asynchronous state in a
spherical description. For this setup, we consider switching the model
variables such that $\theta_{i}$ determines the radial distances, and
$(x_{i},y_{i})$, the azimuthal and polar angles, respectively. The resulting
state shows clustered individuals in the middle of a cloud of loosely behaved
ones, coherent with a description given in [63] (see movie S21 in the
Supplemental Material). We also expect the chaotic states, shown in Fig. 5, to
be useful in the study of populations that exhibit this type of behaviors (see
movie S22 in the Supplemental Material). Finally, our model could be
considered as a basis for engineered systems, similar to the ones presented in
[65, 66, 67].

###### Acknowledgements.

This work was partly supported by FAPESP Grant No. 2021/14335-0 and CNPq Grant
No. 301082/2019-7 (M.A.M.A.) and FAPESP Grant No. 2021/04251-4 (J.U.F.L.).

##  Appendix A Derivation of the ring model

We can collapse the 2D system to a single dimension in space by considering
$\dot{x}_{i}=\dot{y}_{i}$. Hence, we get rid of the second expression in Eqs.
(4). Additionally, if we assume that there is no scale difference between
attraction and repulsion, and that individuals are purely attracting (i.e.
$J=J_{A}=-J_{R}$, for $J_{A}>0$), the equations governing the system’s
dynamics become

$\displaystyle\begin{split}\dot{x}_{i}&=\frac{J}{N}\sum_{j=1}^{N}\sin(x_{ji})\cos(\theta_{ji}),\\\
\dot{\theta}_{i}&=\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_{ji})\left[\cos(x_{ji})\right]^{2}.\end{split}$

Notice that, without loss of generality, we have neglected the scaling factor
in the phase dynamics ($2K\rightarrow K$). This 1D system fits with the
structure described in Eqs. (2), and, the only difference with the 1D model
introduced in [48] is the squared cosine in the phase dynamics. This variation
ultimately improves the driving towards synchronous or asynchronous states
(depending on the value of $K$), as the influence of the individuals’
positions will always be null or positive.

##  Appendix B Perturbation analysis of coherent states

For each state described below, we will be able to rearrange the respective
perturbation dynamics as

$\delta\dot{\vec{r}}_{i}~{}=~{}\mathbf{R}\delta\vec{r}_{i},$ (16)

where $\delta\vec{r}_{i}=\left(\delta x_{i},\delta
y_{i},\delta\theta_{i}\right)$ and

$\mathbf{R}=\begin{pmatrix}\mathbf{R_{1}}&\mathbf{R_{2}}&\mathbf{R_{3}}\\\
\mathbf{R_{4}}&\mathbf{R_{5}}&\mathbf{R_{6}}\\\
\mathbf{R_{7}}&\mathbf{R_{8}}&\mathbf{R_{9}}\end{pmatrix}$ (17)

is a $3N\times 3N$ block matrix where each $N\times N$ block is circulant.
Thus, the stability of Eq. (16) can be analyzed by finding the eigenvalues of
$\mathbf{R}$.

We follow the usual procedure to solve an eigenvalue problem, which starts by
defining
$\mathbf{R_{*}}~{}=~{}\mathbf{R}~{}-~{}\vec{\lambda}\mathbf{I}_{3N\times 3N}$
and then concludes by finding the eigenvalues $\vec{\lambda}$ from the
determinant of $\mathbf{R_{*}}$. Thus, considering that $\mathbf{R_{*}}$ is
also a block matrix, we can rewrite its determinant as

$\det\left(\mathbf{R}_{*}\right)=\det\left(\mathbf{M}\right),$

where

$\displaystyle\begin{split}\mathbf{M}&=\mathbf{R_{1}^{*}}(\mathbf{R_{5}^{*}}\mathbf{R_{9}^{*}}-\mathbf{R_{6}}\mathbf{R_{8}})-\mathbf{R_{2}}\left(\mathbf{R_{4}}\mathbf{R_{9}^{*}}-\mathbf{R_{6}}\mathbf{R_{7}}\right)\\\
&\qquad+\mathbf{R_{3}}\left(\mathbf{R_{4}}\mathbf{R_{8}}-\mathbf{R_{5}^{*}}\mathbf{R_{7}}\right),\end{split}$

since blocks commute, and

$\displaystyle\mathbf{R^{*}_{1}}$
$\displaystyle=\mathbf{R_{1}}-\vec{\lambda}_{1}\mathbf{I}_{N},$
$\displaystyle\mathbf{R^{*}_{5}}$
$\displaystyle=\mathbf{R_{5}}-\vec{\lambda}_{5}\mathbf{I}_{N},$
$\displaystyle\mathbf{R^{*}_{9}}$
$\displaystyle=\mathbf{R_{9}}-\vec{\lambda}_{9}\mathbf{I}_{N}.$

Given that blocks composing $\mathbf{R}$ are circulant, the $N~{}\times~{}N$
matrix $\mathbf{M}$ will also be circulant. Then, we can use the general
solution to find the determinant of circulant matrices as

$\det\left(\mathbf{M}\right)=\prod_{k=0}^{N-1}\sum_{r=0}^{N-1}M_{r+1}\zeta^{kr},$
(18)

from where we are able to obtain the eigenvalues $\vec{\lambda}$. The term
$\zeta=\exp\left\\{2\pi\mathrm{i}\mkern 1.0mu/N\right\\}$ is a primitive
$N$-root of unity. Notice that the structure of $\mathbf{M}$ follows that of
$\mathbf{R}$ in Eq. (17) (i.e. elements’ sub-indexes represent the same
position in the matrix).

###  B.1 Point Synchronous

We add small individual perturbation to the equilibrium states as

$\displaystyle x_{i}$ $\displaystyle=x+\delta x_{i}$ $\displaystyle y_{i}$
$\displaystyle=x+\delta y_{i}$ $\displaystyle\theta_{i}$
$\displaystyle=x+\delta\theta_{i}.$

By plugging these into Eq. (4), we find the perturbation dynamics, governed by

$\displaystyle\begin{split}\delta\dot{x}_{i}&=\frac{J_{A}-J_{R}}{2N}\sum_{j=1}^{N}(\delta
x_{j}-\delta x_{i}),\\\
\delta\dot{y}_{i}&=\frac{J_{A}-J_{R}}{2N}\sum_{j=1}^{N}(\delta y_{j}-\delta
y_{i}),\\\
\delta\dot{\theta}_{i}&=\frac{2K}{N}\sum_{j=1}^{N}(\delta\theta_{j}-\delta\theta_{i}),\end{split}$

which can be arranged as the linear system shown in Eq. (16). In this
particular case, we have that
$\left\\{\mathbf{R_{2}},\mathbf{R_{3}},\mathbf{R_{4}},\mathbf{R_{6}},\mathbf{R_{7}},\mathbf{R_{8}}\right\\}=\mathbf{0}$
and $\mathbf{R_{1}}=\mathbf{R_{5}}$, which simplifies the eigenvalue problem
considerably. Finally, using Eq. (18), we obtain the eigenvalues shown in Eqs.
(5).

###  B.2 Static phase wave I

In this case, the perturbed equilibrium states are defined as

$\displaystyle x_{i}$ $\displaystyle=\frac{2\pi\mathrm{i}\mkern
1.0mu}{N}+x_{0}+\delta x_{i},$ $\displaystyle y_{i}$
$\displaystyle=\frac{2\pi\mathrm{i}\mkern 1.0mu}{N}+y_{0}+\delta y_{i},$
$\displaystyle\theta_{i}$ $\displaystyle=\frac{2\pi\mathrm{i}\mkern
1.0mu}{N}+\theta_{0}+\delta\theta_{i},$

and the perturbation dynamics are governed by

$\displaystyle\begin{split}\delta\dot{x}_{i}&=\frac{1}{2N}\sum_{j}\left[\delta
x_{j}F_{a}^{(I)}(i,j)+\delta\theta_{j}F_{b}^{(I)}(i,j)\right]\\\
&\qquad\qquad+\frac{J_{R}}{2}(\delta x_{i}-\delta\theta_{i}),\\\
\delta\dot{y}_{i}&=\frac{1}{2N}\sum_{j}\left[\delta
y_{j}F_{a}^{(I)}(i,j)+\delta\theta_{j}F_{b}^{(I)}(i,j)\right]\\\
&\qquad\qquad+\frac{J_{R}}{2}(\delta y_{i}-\delta\theta_{i}),\\\
\delta\dot{\theta}_{i}&=\frac{K}{2N}\sum_{j}\Big{[}\delta
x_{j}F_{c}^{(I)}(i,j)+\delta y_{j}F_{c}^{(I)}(i,j)\\\
&\qquad\qquad+\delta\theta_{j}F_{d}^{(I)}(i,j)\Big{]},\end{split}$ (19)

where

$\displaystyle F_{a}^{(I)}(i,j)$
$\displaystyle=J_{A}\cos\left(\frac{4\pi}{N}(j-i)\right)-J_{R},$
$\displaystyle F_{b}^{(I)}(i,j)$
$\displaystyle=J_{R}+J_{A}\cos\left(\frac{4\pi}{N}(j-i)\right)$
$\displaystyle\qquad\qquad-(J_{A}+J_{R})\cos\left(\frac{2\pi}{N}(j-i)\right),$
$\displaystyle F_{c}^{(I)}(i,j)$
$\displaystyle=\cos\left(\frac{6\pi}{N}(j-i)\right)-\cos\left(\frac{2\pi}{N}(j-i)\right),$
$\displaystyle F_{d}^{(I)}(i,j)$
$\displaystyle=\cos\left(\frac{6\pi}{N}(j-i)\right)+3\cos\left(\frac{2\pi}{N}(j-i)\right).$

The positive correlation between $x_{i}$ and $y_{i}$, in the equilibrium,
generates a clear symmetry in position dynamics. Then, as noticed, the
coefficients $F_{a}^{(I)}(i,j)$ and $F_{b}^{(I)}(i,j)$ repeat in the
expressions defining $\delta\dot{x}_{i}$ and $\delta\dot{y}_{i}$, and
$F_{c}^{(I)}(i,j)$ shows up twice in the definition of
$\delta\dot{\theta}_{i}$. These features allow for the simplification of the
linear system when structuring it as in Eq. (16). Thus, for this case we have
that $\mathbf{R_{1}}=\mathbf{R_{5}}$, $\mathbf{R_{3}}=\mathbf{R_{6}}$,
$\mathbf{R_{7}}=\mathbf{R_{8}}$, and
$\left\\{\mathbf{R_{2}},\mathbf{R_{4}}\right\\}=\mathbf{0}$. These
considerations allow us to split the eigenvalue problem into two problems of
lower dimensions such that

$\mathbf{M}^{(I)}=\mathbf{R_{1}^{*}}\left(\mathbf{R_{1}^{*}}\mathbf{R_{9}^{*}}-2\mathbf{R_{3}}\mathbf{R_{7}}\right).$

Then, using Eq. (18), we have that

$\displaystyle\det(\mathbf{R_{1}^{*}})$
$\displaystyle=\prod_{k=0}^{N-1}\Bigg{\\{}-\lambda+\frac{J_{R}}{2}+\frac{1}{2N}\sum_{j=1}^{N}\zeta^{(j-1)k}F_{a}^{(I)}(1,j)\Bigg{\\}},$
$\displaystyle\det\left(\mathbf{R_{1}^{*}}\mathbf{R_{9}^{*}}-2\mathbf{R_{3}}\mathbf{R_{7}}\right)$
$\displaystyle=\prod_{k=0}^{N-1}\Bigg{\\{}\lambda^{2}-\lambda\left[\frac{J_{R}}{2}+\frac{1}{2N}\sum_{j=1}^{N}\zeta^{(j-1)k}\left[F_{a}^{(I)}(1,j)+KF_{d}^{(I)}(1,j)\right]\right]$
$\displaystyle\qquad+\frac{K}{2N}\sum_{j=1}^{N}\zeta^{(j-1)k}\left[J_{R}F_{c}^{(I)}(1,j)+\frac{J_{R}}{2}F_{d}^{(I)}(1,j)-\frac{J_{A}+J_{R}}{2}\cos\left(\frac{2\pi}{N}(j-1)\right)\right]\Bigg{\\}},$

whose solution allows us to find the eigenvalues presented in Eqs. (6).

###  B.3 Static phase wave II

The only difference between this state and the previous one, in the
equilibrium, is that the correlation between $x_{i}$ and $\theta_{i}$ is now
negative. Then, after perturbing the equilibrium states individually, we have
that

$\displaystyle x_{i}$ $\displaystyle=\frac{2\pi\mathrm{i}\mkern
1.0mu}{N}+x_{0}+\delta x_{i},$ $\displaystyle y_{i}$
$\displaystyle=\frac{2\pi\mathrm{i}\mkern 1.0mu}{N}+y_{0}+\delta y_{i},$
$\displaystyle\theta_{i}$ $\displaystyle=-\frac{2\pi\mathrm{i}\mkern
1.0mu}{N}+\theta_{0}+\delta\theta_{i},$

and the perturbation dynamics are governed by

$\displaystyle\begin{split}\delta\dot{x}_{i}&=\frac{1}{2N}\sum_{j}\left[\delta
x_{j}F_{a}^{(II)}(i,j)+\delta\theta_{j}F_{b}^{(II)}(i,j)\right]\\\
&\qquad\qquad-\frac{J_{A}}{2}(\delta x_{i}+\delta\theta_{i}),\\\
\delta\dot{y}_{i}&=\frac{1}{2N}\sum_{j}\left[\delta
y_{j}F_{a}^{(II)}(i,j)+\delta\theta_{j}F_{b}^{(II)}(i,j)\right]\\\
&\qquad\qquad-\frac{J_{A}}{2}(\delta y_{i}+\delta\theta_{i}),\\\
\delta\dot{\theta}_{i}&=\frac{K}{2N}\sum_{j}\Big{[}\delta
x_{j}F_{c}^{(II)}(i,j)+\delta y_{j}F_{c}^{(II)}(i,j)\\\
&\qquad\qquad+\delta\theta_{j}F_{d}^{(II)}(i,j)\Big{]},\end{split}$

where

$\displaystyle F_{a}^{(II)}(i,j)$
$\displaystyle=J_{A}-J_{R}\cos\left(\frac{4\pi}{N}(j-i)\right),$
$\displaystyle F_{b}^{(II)}(i,j)$
$\displaystyle=J_{A}+J_{R}\cos\left(\frac{4\pi}{N}(j-i)\right)$
$\displaystyle\qquad\qquad-(J_{A}+J_{R})\cos\left(\frac{2\pi}{N}(j-i)\right),$
$\displaystyle F_{c}^{(II)}(i,j)$
$\displaystyle=\cos\left(\frac{2\pi}{N}(j-i)\right)-\cos\left(\frac{6\pi}{N}(j-i)\right),$
$\displaystyle F_{d}^{(II)}(i,j)$
$\displaystyle=3\cos\left(\frac{2\pi}{N}(j-i)\right)+\cos\left(\frac{6\pi}{N}(j-i)\right).$

Notice that these equations also show the recurrence of some coefficients in
the definition of $\delta\dot{x}_{i}$, $\delta\dot{y}_{i}$, and
$\delta\dot{\theta}_{i}$. This consideration leads to the same simplifications
described for the previous state when structuring the system as in Eqs. (16).
So we have that, $\mathbf{R_{1}}=\mathbf{R_{5}}$,
$\mathbf{R_{3}}=\mathbf{R_{6}}$, $\mathbf{R_{7}}=\mathbf{R_{8}}$, and
$\left\\{\mathbf{R_{2}},\mathbf{R_{4}}\right\\}=\mathbf{0}$.

$\mathbf{M}^{(II)}=\mathbf{R_{1}^{*}}\left(\mathbf{R_{1}^{*}}\mathbf{R_{9}^{*}}-2\mathbf{R_{3}}\mathbf{R_{7}}\right).$

$\displaystyle\det(\mathbf{R_{1}^{*}})$
$\displaystyle=\prod_{k=0}^{N-1}\Bigg{\\{}-\lambda-\frac{J_{A}}{2}+\frac{1}{2N}\sum_{j=1}^{N}\zeta^{(j-1)k}F_{a}^{(II)}(1,j)\Bigg{\\}},$
$\displaystyle\det\left(\mathbf{R^{*}_{1}R^{*}_{9}}-2\mathbf{R_{3}R_{7}}\right)$
$\displaystyle=\prod_{k=0}^{N-1}\Bigg{\\{}\lambda^{2}-\lambda\left[-\frac{J_{A}}{2}+\frac{1}{2N}\sum_{j=1}^{N}\zeta^{(j-1)k}\left(KF_{d}^{(II)}(1,j)+F_{a}^{(II)}(1,j)\right)\right]$
$\displaystyle\qquad+\frac{K}{2N}\sum_{j=1}^{N}\zeta^{(j-1)k}\left[J_{A}F_{c}^{(II)}(1,j)-\frac{J_{A}}{2}F_{d}^{(II)}(1,j)+\frac{(J_{A}+J_{R})}{2}\cos\left(\frac{2\pi}{N}(j-1)\right)\right]\Bigg{\\}},$

which lead to the eigenvalues presented in Eqs. (7).

###  B.4 Static phase wave III

For this case, the perturbed equilibrium states are defined as

$\displaystyle x_{i}$ $\displaystyle=\frac{2\pi\mathrm{i}\mkern
1.0mu}{N}+x_{0}+\delta x_{i},$ $\displaystyle y_{i}$
$\displaystyle=-\frac{2\pi\mathrm{i}\mkern 1.0mu}{N}+y_{0}+\delta y_{i},$
$\displaystyle\theta_{i}$ $\displaystyle=\frac{2\pi\mathrm{i}\mkern
1.0mu}{N}+\theta_{0}+\delta\theta_{i},$

and the perturbation dynamics are governed by

$\displaystyle\begin{split}\delta\dot{x}_{i}&=\frac{1}{2N}\sum_{j}\left[\delta
x_{j}F_{a}^{(III)}(i,j)+\delta\theta_{j}F_{b}^{(III)}(i,j)\right]\\\
&\qquad\qquad+\frac{J_{R}}{2}(\delta x_{i}-\delta\theta_{i}),\\\
\delta\dot{y}_{i}&=\frac{1}{2N}\sum_{j}\left[\delta
y_{j}F_{c}^{(III)}(i,j)+\delta\theta_{j}F_{d}^{(III)}(i,j)\right]\\\
&\qquad\qquad-\frac{J_{A}}{2}(\delta y_{i}+\delta\theta_{i}),\\\
\delta\dot{\theta}_{i}&=\frac{1}{2N}\sum_{j}\Big{[}\delta
x_{j}F_{e}^{(III)}(i,j)-\delta y_{j}F_{e}^{(III)}(i,j)\\\
&\qquad\qquad+\delta\theta_{j}F_{f}^{(III)}(i,j)\Big{]},\end{split}$

where

$\displaystyle F_{a}^{(III)}(i,j)$
$\displaystyle=J_{A}\cos\left(\frac{4\pi}{N}(j-i)\right)-J_{R},$
$\displaystyle F_{b}^{(III)}(i,j)$
$\displaystyle=J_{R}+J_{A}\cos\left(\frac{4\pi}{N}(j-i)\right)-(J_{A}+J_{R})\cos\left(\frac{2\pi}{N}(j-i)\right),$
$\displaystyle F_{c}^{(III)}(i,j)$
$\displaystyle=J_{A}-J_{R}\cos\left(\frac{4\pi}{N}(j-i)\right),$
$\displaystyle F_{d}^{(III)}(i,j)$
$\displaystyle=J_{A}+J_{R}\cos\left(\frac{4\pi}{N}(j-i)\right)-(J_{A}+J_{R})\cos\left(\frac{2\pi}{N}(j-i)\right),$
$\displaystyle F_{e}^{(III)}(i,j)$
$\displaystyle=K\cos\left(\frac{6\pi}{N}(j-i)\right)-K\cos\left(\frac{2\pi}{N}(j-i)\right),$
$\displaystyle F_{f}^{(III)}(i,j)$
$\displaystyle=3K\cos\left(\frac{2\pi}{N}(j-i)\right)+K\cos\left(\frac{6\pi}{N}(j-i)\right).$

Despite the slight difference in the definition of equilibrium states compared
to the previous cases (types I and II), simplifications are not possible in
the perturbation equations of motion. As noticed, the only repeating
coefficient is $F_{e}^{(III)}(i,j)$ in the definition of
$\delta\dot{\theta}_{i}$. Hence, when structuring the perturbation dynamics as
in Eqs. (16), the only consideration we can make is that
$\mathbf{R_{7}}=-\mathbf{R_{8}}$. This drawback forces us to write the
determinant, using Eq. (18), as

$\displaystyle\det(\mathbf{M}^{(III)})$
$\displaystyle=\prod_{k=0}^{N-1}\Bigg{\\{}-\lambda^{3}+\lambda^{2}\left[\frac{J_{R}-J_{A}}{2}+\frac{1}{2N}S_{1}\right]+\lambda\left[\frac{J_{A}J_{R}}{4}+\frac{J_{A}J_{R}}{4N}T_{1}+\frac{J_{A}}{4N}S_{2}-\frac{J_{R}}{4N}S_{3}\right]$
$\displaystyle\qquad\qquad+\frac{K}{16N}\left(J_{A}+J_{R}\right)^{2}-\frac{J_{A}J_{R}}{8N}S_{4}\Bigg{\\}},$

where

$\displaystyle S_{1}$
$\displaystyle=\sum_{j=1}^{N}\zeta^{(j-1)k}\Big{(}F^{(III)}_{a}(1,j)+F^{(III)}_{c}(1,j)+KF^{(III)}_{f}(1,j)\Big{)},$
$\displaystyle S_{2}$
$\displaystyle=\sum_{j=1}^{N}\zeta^{(j-1)k}\Big{(}F_{a}^{(III)}(1,j)+KF_{e}^{(III)}(1,j)+KF_{f}^{(III)}(1,j)\Big{)},$
$\displaystyle S_{3}$
$\displaystyle=\sum_{j=1}^{N}\zeta^{(j-1)k}\Big{(}F_{c}^{(III)}(1,j)+KF_{e}^{(III)}(1,j)+KF_{f}^{(III)}(1,j)\Big{)},$
$\displaystyle S_{4}$
$\displaystyle=\sum_{j=1}^{N}\zeta^{(j-1)k}\left(2F_{e}^{(III)}(1,j)+F_{f}^{(III)}(1,j)\right),$
$\displaystyle T_{1}$
$\displaystyle=\sum_{j=1}^{N}\zeta^{(j-1)k}+0.25\left(\zeta^{(j-1)(k+2)}+\zeta^{(j-1)(k-2)}\right),$
$\displaystyle T_{2}$
$\displaystyle=0.5\sum_{j=1}^{k}\zeta^{(j-1)(k+1)}+\zeta^{(j-1)(k-1)}.$

To find the eigenvalues, we must address $N$ cubic equations of the form
$-\lambda^{3}+\lambda^{2}\beta+\lambda\xi+\eta$. We solve these using the
general formula for

$\beta=\begin{cases}0,&\text{for }k=0;\\\ 0.5(J_{R}-J_{A})+1.5K,&\text{for
}k=1,N-1;\\\ 0.25(J_{R}-J_{A}),&\text{for }k=2,N-2;\\\
0.5(J_{R}-J_{A})+0.25K,&\text{for }k=3,N-3;\\\
0.5(J_{R}-J_{A}),&\text{otherwise};\\\ \end{cases}$ (20)
$\xi=\begin{cases}0,&\text{for }k=0,\\\
0.25(J_{A}J_{R}+K(J_{A}-J_{R})),&\text{for }k=1,N-1;\\\
0.125(J_{A}^{2}+J_{R}^{2}+2.5J_{A}J_{R}),&\text{for }k=2,N-2;\\\
0.25(J_{A}J_{R}+K(J_{A}-J_{R})),&\text{for }k=3,N-3;\\\
0.25KJ_{A}J_{R},&\text{otherwise};\\\ \end{cases}$ (21)

and

$\eta=\begin{cases}0,&\text{for }k=0;\\\
0.125K[0.25(J_{A}+J_{R})^{2}&{}\textstyle\\\ \qquad-J_{A}J_{R}],&\text{for
}k=1,N-1;\\\ 0,&\text{for }k=2,N-2;\\\ -0.1875KJ_{A}J_{R},&\text{for
}k=3,N-3;\\\ 0,&\text{otherwise}.\\\ \end{cases}$

##  Appendix C Perturbation analysis of incoherent states

In general, we perturb the equilibrium state $\rho_{0}$ by a small quantity
$\delta\rho$, such that $\rho=\rho_{0}+\delta\rho$. Then, the temporal
evolution of the perturbation is governed by

$\frac{\partial}{\partial t}\delta\rho=-\vec{\nabla}(\delta\rho)\vec{v}.$ (22)

From Eq. (8), we know that

$\int_{0}^{2\pi}\delta\rho(x,y,\theta,t)\mathrm{d}x\mathrm{d}y\mathrm{d}\theta=0,$

and to first order in $\delta\rho(x,\theta,t)$, we have that

$S_{\sigma}^{1}e^{\mathrm{i}\mkern
1.0mu\phi_{\sigma}}=\int_{0}^{2\pi}e^{\mathrm{i}\mkern
1.0mu\sigma}\delta\rho(x,y,\theta,t)\mathrm{d}x\mathrm{d}y\mathrm{d}\theta.$
(23)

We expand $\delta\rho$ in Fourier series as

$\delta\rho=\sum_{m,n,l}f_{m,n,l}(t)e^{\mathrm{i}\mkern
1.0mu(mx+ny+l\theta)},$ (24)

and we will solve equations of the type $f(t)=\bar{f}e^{\lambda t}$ to get the
eigenvalues $\lambda$.

###  C.1 Distributed asynchronous

We consider the equilibrium $\rho_{0}^{A}$ defined in Eqs. (12). Using Eqs.
(22) and (23), we obtain the perturbation dynamics

$\displaystyle\begin{split}\frac{d}{dt}\delta\rho&=\frac{J_{A}}{16\pi^{3}}\left[S_{(\theta+x)}^{1}\cos\left(\phi_{(\theta+x)}-\theta-x\right)+S_{(\theta+y)}^{1}\cos\left(\phi_{(\theta+y)}-\theta-y\right)\right]\\\
&\qquad-\frac{J_{R}}{16\pi^{3}}\left[S_{(\theta-x)}^{1}\cos\left(\phi_{(\theta-x)}-\theta+x\right)+S_{(\theta-y)}^{1}\cos\left(\phi_{(\theta-y)}-\theta+y\right)\right]\\\
&\qquad\qquad+\frac{K}{32\pi^{3}}\bigg{[}S_{(\theta+2x)}^{1}\cos\left(\phi_{(\theta+2x)}-\theta-2x\right)+S_{(\theta-2x)}^{1}\cos\left(\phi_{(\theta-2x)}-\theta+2x\right)\\\
&\qquad\qquad\qquad+S_{(\theta+2y)}^{1}\cos\left(\phi_{(\theta+2y)}-\theta-2y\right)+S_{(\theta-2y)}^{1}\cos\left(\phi_{(\theta-2y)}-\theta+2y\right)+4S_{\theta}^{1}\cos\left(\phi_{\theta}-\theta\right)\bigg{]}.\end{split}$
(25)

We expand this equation following the form of Eq. (24) and we will see that
the only relevant terms correspond to

$\displaystyle\begin{split}\\{(m,n,l)\\}&=\\{(-1,0,1);(0,-1,-1);(1,0,-1);\\\
&\qquad\quad(0,1,-1);(0,0,-1);(2,0,-1);\\\
&\qquad\qquad(-2,0,-1);(0,2,-1);(0,-2,-1)\\}.\end{split}$ (26)

Then, the solutions give us the eigenvalues shown in Eqs. (14).

###  C.2 Distributed synchronous

Given that in this state phases are synchronized and only positions are
uniformly distributed, we uncouple the system’s dynamics. The incoherence
between phases and positions will lead to zeroing order parameters
$S_{(\theta\pm 2x)}$ and $S_{(\theta\pm 2y)}$. That allows us to rewrite

$v_{\theta}=KS_{\theta}\sin(\phi_{\theta}-\theta),$

which is exactly the mean-field description of the Kuramoto model without
considering natural frequencies [18]. Then, the condition to reach
synchronization is given by $(K>0)$.

In the equilibrium, the synchronization of phases will lead to
$\theta_{i}=\theta$. This condition allows for the equivalences

$\displaystyle\begin{split}S_{(\theta\pm x)}e^{\mathrm{i}\mkern
1.0mu\phi_{(\theta\pm x)}}&=S_{x}e^{\mathrm{i}\mkern
1.0mu(\theta\pm\phi_{x})},\\\ S_{(\theta\pm y)}e^{\mathrm{i}\mkern
1.0mu\phi_{(\theta\pm y)}}&=S_{y}e^{\mathrm{i}\mkern
1.0mu(\theta\pm\phi_{y})},\end{split}$

and therefore, the velocity field for this state is defined as

$\displaystyle\begin{split}v_{x}&=\frac{J_{A}-J_{R}}{2}S_{x}\sin\left(\phi_{x}-x\right),\\\
v_{y}&=\frac{J_{A}-J_{R}}{2}S_{y}\sin\left(\phi_{y}-y\right).\end{split}$ (27)

In this case we use the equilibrium $\rho_{0}^{S}$ from Eqs. (12), and
considering (23), the perturbation dynamics are governed by

$\frac{d}{dt}\delta\rho=\frac{J_{A}-J_{R}}{8\pi^{2}}\left[S_{x}^{1}\cos\left(\phi_{x}-x\right)+S_{y}^{1}\cos\left(\phi_{y}-y\right)\right].$
(28)

Then, using Eq. (24) for this equation, we find the eigenvalue shown in Eq.
(15).

## References

  * Larsson [2012] M. Larsson, Current Zoology 58, 116 (2012).
  * Marras  _et al._ [2015] S. Marras, S. S. Killen, J. Lindström, D. J. McKenzie, J. F. Steffensen, and P. Domenici, Behavioral ecology and sociobiology 69, 219 (2015).
  * Emlen [1952] J. T. Emlen, The Auk 69, 160 (1952).
  * Bialek  _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Proceedings of the National Academy of Sciences 109, 4786 (2012).
  * Sullivan [1981] R. T. Sullivan, The Florida Entomologist 64, 44 (1981).
  * Chittka and Rossi [2022] L. Chittka and N. Rossi, Trends in Cognitive Sciences (2022).
  * Kearns [2010] D. B. Kearns, Nature Reviews Microbiology 8, 634 (2010).
  * Ariel  _et al._ [2015] G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M. Harshey, and A. Be’Er, Nature communications 6, 8396 (2015).
  * Be’er and Ariel [2019] A. Be’er and G. Ariel, Movement ecology 7, 1 (2019).
  * Néda  _et al._ [2000] Z. Néda, E. Ravasz, Y. Brechet, T. Vicsek, and A.-L. Barabási, Nature 403, 849 (2000).
  * Strogatz  _et al._ [2005] S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott, Nature 438, 43 (2005).
  * Yaniv  _et al._ [2014] Y. Yaniv, I. Ahmet, J. Liu, A. E. Lyashkov, T.-R. Guiriba, Y. Okamoto, B. D. Ziman, and E. G. Lakatta, Heart Rhythm 11, 1210 (2014).
  * Bychkov  _et al._ [2020] R. Bychkov, M. Juhaszova, K. Tsutsui, C. Coletta, M. D. Stern, V. A. Maltsev, and E. G. Lakatta, Clinical Electrophysiology 6, 907 (2020).
  * Aihara  _et al._ [2006] I. Aihara, H. Kitahata, K. Aihara, and K. Yoshikawa, METR 35, 1 (2006).
  * Aihara  _et al._ [2007] I. Aihara, S. Horai, H. Kitahata, K. Aihara, and K. Yoshikawa, IEICE transactions on fundamentals of electronics, communications and computer sciences 90, 2154 (2007).
  * Aihara  _et al._ [2008] I. Aihara, H. Kitahata, K. Yoshikawa, and K. Aihara, Artificial Life and Robotics 12, 29 (2008).
  * Uhlhaas and Singer [2006] P. J. Uhlhaas and W. Singer, neuron 52, 155 (2006).
  * Strogatz [2000] S. H. Strogatz, Physica D: Nonlinear Phenomena 143, 1 (2000).
  * Aoyagi [1995] T. Aoyagi, Physical review letters 74, 4075 (1995).
  * Couzin  _et al._ [2002] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Journal of theoretical biology 218, 1 (2002).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Physical review letters 75, 1226 (1995).
  * O’Keeffe  _et al._ [2017] K. P. O’Keeffe, H. Hong, and S. H. Strogatz, Nature communications 8, 1504 (2017).
  * Zhang  _et al._ [2020] B. Zhang, A. Sokolov, and A. Snezhko, Nature communications 11, 4401 (2020).
  * Creppy  _et al._ [2016] A. Creppy, F. Plouraboué, O. Praud, X. Druart, S. Cazin, H. Yu, and P. Degond, Journal of The Royal Society Interface 13, 20160575 (2016).
  * Maggi  _et al._ [2023] C. Maggi, F. Saglimbeni, V. C. Sosa, R. Di Leonardo, B. Nath, and A. Puglisi, PRX Life 1, 013003 (2023).
  * Creppy  _et al._ [2015] A. Creppy, O. Praud, X. Druart, P. L. Kohnke, and F. Plouraboué, Physical Review E 92, 032722 (2015).
  * Tan  _et al._ [2022] T. H. Tan, A. Mietke, J. Li, Y. Chen, H. Higinbotham, P. J. Foster, S. Gokhale, J. Dunkel, and N. Fakhri, Nature 607, 287 (2022).
  * Aihara  _et al._ [2014] I. Aihara, T. Mizumoto, T. Otsuka, H. Awano, K. Nagira, H. G. Okuno, and K. Aihara, Scientific reports 4, 3891 (2014).
  * Peshkov  _et al._ [2022] A. Peshkov, S. McGaffigan, and A. C. Quillen, Soft Matter 18, 1174 (2022).
  * Yuan  _et al._ [2014] J. Yuan, D. M. Raizen, and H. H. Bau, Proceedings of the National Academy of Sciences 111, 6865 (2014).
  * Ērglis  _et al._ [2007] K. Ērglis, Q. Wen, V. Ose, A. Zeltins, A. Sharipo, P. A. Janmey, and A. Cēbers, Biophysical journal 93, 1402 (2007).
  * Sarfati  _et al._ [2020] R. Sarfati, J. C. Hayes, É. Sarfati, and O. Peleg, Journal of The Royal Society Interface 17, 20200179 (2020).
  * Quillen  _et al._ [2022] A. Quillen, A. Peshkov, B. Chakrabarti, N. Skerrett, S. McGaffigan, and R. Zapiach, Physical Review E 106, 064401 (2022).
  * Barotta  _et al._ [2023] J.-W. Barotta, S. J. Thomson, L. F. Alventosa, M. Lewis, and D. M. Harris, Communications Physics 6, 87 (2023).
  * Manna  _et al._ [2021] R. K. Manna, O. E. Shklyaev, and A. C. Balazs, Proceedings of the National Academy of Sciences 118, e2022987118 (2021).
  * Ali  _et al._ [2023] N. Ali, S. Nand, A. Kiran, M. Mishra, and V. Mehandia, Physics of Fluids 35 (2023).
  * Lizarraga and de Aguiar [2020] J. U. Lizarraga and M. A. de Aguiar, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2020).
  * Ceron  _et al._ [2023] S. Ceron, K. O’Keeffe, and K. Petersen, Nature Communications 14, 940 (2023).
  * Hong [2018] H. Hong, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (2018).
  * McLennan-Smith  _et al._ [2020] T. A. McLennan-Smith, D. O. Roberts, and H. S. Sidhu, Physical Review E 102, 032607 (2020).
  * Sar  _et al._ [2022] G. K. Sar, S. N. Chowdhury, M. Perc, and D. Ghosh, New Journal of Physics 24, 043004 (2022).
  * Belovs  _et al._ [2017] M. Belovs, R. Livanovičs, and A. Cēbers, Physical Review E 96, 042408 (2017).
  * Japón  _et al._ [2022] P. Japón, F. Jiménez-Morales, and F. Casares, Cells & Development 169, 203726 (2022).
  * Adorjani  _et al._ [2023] B. Adorjani, A. Libal, C. Reichhardt, and C. Reichhardt, arXiv preprint arXiv:2309.10937 (2023).
  * Blum  _et al._ [2022] N. Blum, A. Li, K. O’Keeffe, and O. Kogan, arXiv preprint arXiv:2210.11417 (2022).
  * Kongni  _et al._ [2023] S. J. Kongni, V. Nguefoue, T. Njougouo, P. Louodop, F. F. Ferreira, R. Tchitnga, and H. A. Cerdeira, Physical Review E 108, 034303 (2023).
  * Ghosh  _et al._ [2023] S. Ghosh, G. K. Sar, S. Majhi, and D. Ghosh, Physical Review E 108, 034217 (2023).
  * O’Keeffe  _et al._ [2022] K. O’Keeffe, S. Ceron, and K. Petersen, Physical Review E 105, 014211 (2022).
  * Yoon  _et al._ [2022] S. Yoon, K. O’Keeffe, J. Mendes, and A. Goltsev, Physical Review Letters 129, 208002 (2022).
  * O’Keeffe and Hong [2022] K. O’Keeffe and H. Hong, Physical Review E 105, 064208 (2022).
  * Hao  _et al._ [2023] B. Hao, M. Zhong, and K. O’Keeffe, arXiv preprint arXiv:2309.02342 (2023).
  * Lizárraga and de Aguiar [2023] J. U. Lizárraga and M. A. de Aguiar, Physical Review E 108, 024212 (2023).
  * Sar  _et al._ [2023] G. K. Sar, D. Ghosh, and K. O’Keeffe, Physical Review E 107, 024215 (2023).
  * Lee  _et al._ [2021] H. K. Lee, K. Yeo, and H. Hong, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021).
  * Jiménez-Morales [2020] F. Jiménez-Morales, Physical Review E 101, 062202 (2020).
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, Reviews of modern physics 77, 137 (2005).
  * Kruk  _et al._ [2018] N. Kruk, Y. Maistrenko, and H. Koeppl, Physical Review E 98, 032219 (2018).
  * Kruk  _et al._ [2020] N. Kruk, J. A. Carrillo, and H. Koeppl, Physical Review E 102, 022604 (2020).
  * Childs and Strogatz [2008] L. M. Childs and S. H. Strogatz, Chaos: An Interdisciplinary Journal of Nonlinear Science 18 (2008).
  * López  _et al._ [2015] H. M. López, J. Gachelin, C. Douarche, H. Auradou, and E. Clément, Physical review letters 115, 028301 (2015).
  * Quillen  _et al._ [2021] A. Quillen, A. Peshkov, E. Wright, and S. McGaffigan, Physical Review E 104, 014412 (2021).
  * Levis  _et al._ [2019] D. Levis, I. Pagonabarraga, and B. Liebchen, Physical Review Research 1, 023026 (2019).
  * Cavagna  _et al._ [2023] A. Cavagna, I. Giardina, M. A. Gucciardino, G. Iacomelli, M. Lombardi, S. Melillo, G. Monacchia, L. Parisi, M. J. Peirce, and R. Spaccapelo, Scientific Reports 13, 8745 (2023).
  * Facchinelli  _et al._ [2015] L. Facchinelli, L. Valerio, R. S. Lees, C. F. Oliva, T. Persampieri, C. M. Collins, A. Crisanti, R. Spaccapelo, and M. Q. Benedict, Malaria Journal 14, 1 (2015).
  * Togashi [2019] Y. Togashi, The Journal of Physical Chemistry B 123, 1481 (2019).
  * Zhou  _et al._ [2020] C. Zhou, N. J. Suematsu, Y. Peng, Q. Wang, X. Chen, Y. Gao, and W. Wang, ACS nano 14, 5360 (2020).
  * Chen  _et al._ [2023] B. Chen, H. Tan, M. Ding, L. Liu, S. Wang, X. Peng, H. Tian, J. Jiang, J. Gao, W. Huang, _et al._ , ACS nano 17, 13826 (2023).