# Collective behaviour of swarmalators on a 1D ring

Kevin O’Keeffe  Senseable City Lab, Massachusetts Institute of Technology,
Cambridge, MA 02139  Steven Ceron  Sibley School of Mechanical and Aerospace
Engineering, Cornell University, Ithaca, NY 14853, USA  Kirstin Petersen
Department of Electrical and Computer Engineering, Cornell University, Ithaca,
NY 14853, USA

###### Abstract

We study the collective behavior of swarmalators, generalizations of phase
oscillators that both sync and swarm, confined to move on a 1D ring. This
simple model captures some of the essence of movement in 2D or 3D but has the
benefit of being solvable: most of the collective states and their
bifurcations can be specified exactly. The model also captures the behavior of
real-world swarmalators which swarm in quasi-1D rings such as bordertaxic
sperm and vinegar eels.

##  I Introduction

Synchronization and swarming are universal phenomena Winfree (2001); Kuramoto
(2003); Pikovsky et al. (2003); Toner and Tu (1995); Vicsek et al. (1995) that
are in a sense spatiotemporal opposites. Synchronizing units self-organize in
time, but not in space; laser arrays fire simultaneously Jiang and McCall
(1993); Kozyreff et al. (2000), heart cells trigger all at once Peskin (1975),
but neither system has a notion of spontaneous group movement. Swarming units
flip the picture: they self-organize in space, not time. Birds fly in flocks
Bialek et al. (2012), fish swim in schools Katz et al. (2011), but neither
coordinates the timing of an internal state or rhythm.

The units of some systems appear to self-organize in both space and time. In
Biology, sperm Yang et al. (2008); Riedel et al. (2005), vinegar eels Quillen
et al. (2021a, b), and other microswimmers Taylor (1951); Tamm et al. (1975);
Huang et al. (2021) synchronize the wriggling of their tails which is
speculated to hydrodynamically couple to their motion. In Chemistry, magnetic
Janus particles Yan et al. (2012, 2015); Hwang et al. (2020), dieletric Quinke
rollers Zhang et al. (2020); Bricard et al. (2015); Zhang et al. (2021), and
other active entities Manna et al. (2021) lock their rotations enabling a kind
of sync-dependent self-assembly. In Engineering, land-based robots and aerial
drones can be programmed to swarm based on their synchronizable internal
clocks Barciś et al. (2019); Barciś and Bettstetter (2020); Monaco et al.
(2020). Sync and swarming are also suspected to interact in spatial cognition
Monaco et al. (2019), embryology Uriu and Morelli (2017, 2014); Tsiairis and
Aulehla (2016) and the physics of magnetic domain walls Hrabec et al. (2018).

The theoretical study of systems which both sync and swarm is just beginning.
Takana pioneered the initiative by formulating a model of chemotactic
oscillators Tanaka (2007) and found diverse phenomena Iwasa and Tanaka (2017,
2010). Leibchen and Levis generalized the Vicsek model and found new types of
long range synchrony Levis et al. (2019). O’Keeffe et al introduced a model of
‘swarmalators’ 111Short for ‘swarming oscillators’. whose collective states
have been realized in Nature and technology Barciś et al. (2019); Barciś and
Bettstetter (2020); Zhang et al. (2020) and is being actively extended. The
inclusion of noise Hong (2018), local coupling Lee et al. (2021); Jiménez-
Morales (2020), periodic forcing Lizarraga and de Aguiar (2020), mixed sign
interactions McLennan-Smith et al. (2020), and finite $N$ effects O’Keeffe et
al. (2018) have been studied. The potential of swarmalators in bio-inspired
computing has been explored O’Keeffe and Bettstetter (2019), as has the well-
posedness of $N\rightarrow\infty$ solutions of the swarmalator model Ha et al.
(2021, 2019).

Here we study swarmalators confined to move on a 1D ring. Our aim is two-fold.
First, to model swarmalators which swarm purely in 1D. Frogs, nematodes, and
other organisms are often bordertaxic, seeking out the ring-like edges of
their confining geometry Bau et al. (2015); Yuan et al. (2015); Nosrati et al.
(2016); Ketzetzi et al. (2021); Creppy et al. (2016); Aihara (2009); Aihara et
al. (2014). Janus particles, when acting as microrobots for precision medicine
Le et al. (2019); Agrawal and Agrawal (2019); Yi et al. (2016), will need to
navigate pseudo 1D grooves and channels Xiao et al. (2018). Toy models for
these systems will be useful for applied research.

Figure 1: Order parameters of the 2D swaramalator model
$S_{\pm}e^{i\phi_{\pm}}:=(N)^{-1}\sum_{j}e^{i(\phi_{j}\pm\theta_{j})}$, where
$\phi,\theta$ are the spatial angle and phase of swarmalators. $S_{+}$
bifurcates from 0 at $K_{1}$ as the static async state (Fig 2(b))
destabilizes, $S_{-}$ from $0$ at $K_{2}$ as the active phase wave (Fig 2(f))
destabilizes. The nature of these bifurcations, as well as analytic
expressions for both $K_{1}$ and $K_{2}$, are unknown.

Our second aim is theoretical, namely to investigate the original 2D
swarmalator model whose physics is not understood 222Note the Figure displays
the order parameters of the instance of the model; see Appendix B. Figure 1
shows its order parameters
$S_{\pm}:=(N)^{-1}\sum_{j}e^{i(\phi_{j}\pm\theta_{j})}$ – where $\phi,\theta$
denote the spatial angle and phase – dependence on the coupling strength $K$.
At a critical $K_{1}$, $S_{+}$ jumps from zero as the system transitions from
a static async state (Fig 2(b)) to an active phase wave state (Fig 2(e)) in
which swarmalators run in a space-phase vortex. As expected of order
parameters, $S_{+}$ grows as $K$ is increased. But then it begins to decline
at a second value $K_{2}$ as the swarmalator vortex bifurcates into a broken
band of mini-vortices (Fig 2(e)). The cause of this non-monotinicity, as well
as analytic values for $K_{1},K_{2}$, are unknown. Like the old puzzles to
understand the transitions of the Kuramoto model Kuramoto (2003); Strogatz
(2000); Strogatz and Mirollo (1991); Mirollo and Strogatz (2007); Crawford
(1994), the bifurcations of the swarmalator model ‘cry out for a theoretical
explanation’ Strogatz (2000).

Our hope is that a retreat to a 1D ring will give some first clues on how to
provide such an explanation.

Figure 2: Collective states of the 2D swarmalator model introduced in O’Keeffe
et al. (2017) where swarmalators are represented as colored dots where the
color refers to the swarmalators phase. As described in Appendix B, the states
displayed are from a slightly different intstance of the model to that
presented in O’Keeffe et al. (2017) which produces the same qualitative
behavior, but better emphasises the non-monotonic behavior of the order
parameters $S{\pm}$ shown in Figure 1. In all panels a Euler method was used
with timestep $dt=0.1$ for $T=1000$ units for $N=1000$ swarmalators. (a)
Static sync: $(J,K,\sigma)=(1,1,10)$ (b) Static async $(J,K,\sigma)=(1,1,10)$
(c) Static phase wave $(J,K,\sigma)=(1,1,10)$ (d) Splintered phase wave
$(J,K,\sigma)=(1,1,10)$ (e) Active phase wave. In the three static states
(a)-(c) swarmalators do not move in space or phase. In the splintered phase
wave, each colored chunk is a vortex: the swarmalators librate in both space
and phase. In the active phase wave, the librations are excited into
rotations; the swarmalators split into counter-rotating groups as indicated by
the black arrows.

##  II Model

We study a pair of modified Kuramoto models,

$\displaystyle\dot{x_{i}}$
$\displaystyle=\nu_{i}+\frac{J}{N}\sum_{j}^{N}\sin(x_{j}-x_{i})\cos(\theta_{j}-\theta_{i})$
(1) $\displaystyle\dot{\theta_{i}}$
$\displaystyle=\omega_{i}+\frac{K}{N}\sum_{j}^{N}\sin(\theta_{j}-\theta_{i})\cos(x_{j}-x_{i})$
(2)

where $(x_{i},\theta_{i})\in(S^{1},S^{1})$ are the position and phase of the
$i$-th swarmalator for $i=1,\dots,N$ and ($\nu_{i},\omega_{i}$), $(J,K)$ are
the associated natural frequencies and couplings. We consider identical
swarmalators $(\omega_{i},\nu_{i})=(\omega,\nu)$ and by a change of frame set
$\omega=\nu=0$ WLOG.

Eq. (2) models position-dependent synchronization. The familiar Kuramoto sine
term minimizes swarmalators’ pairwise phase difference (so they sync) while
the new cosine term strengthens the coupling between nearby swarmalators
$K_{ij}=K\cos(x_{j}-x_{i})$ (so the sync is position dependent). Eq. (1) is
Eq. (2)’s mirror-image: it models phase-dependent swarming. Now the sine term
minimizes swarmalators’ pairwise distances (so they swarm / aggregate) and the
cosine term strengthens the coupling between similarly phased swarmalators
$J_{ij}=J\cos(\theta_{j}-\theta_{i})$. You can also think of Eqs.(1) and (2)
as modelling synchronization on the unit torus.

Converting the trig functions to complex exponentials and rearranging makes
the model even simpler

$\displaystyle\dot{x_{i}}$
$\displaystyle=\frac{J}{2}\Big{(}S_{+}\sin(\Phi_{+}-(x_{i}+\theta_{i}))+S_{-}\sin(\Phi_{-}-(x_{i}-\theta_{i}))\Big{)}$
(3) $\displaystyle\dot{\theta_{i}}$
$\displaystyle=\frac{K}{2}\Big{(}S_{+}\sin(\Phi_{+}-(x_{i}+\theta_{i}))-S_{-}\sin(\Phi_{-}-(x_{i}-\theta_{i}))\Big{)}$
(4)

where

$W_{\pm}=S_{\pm}e^{i\Phi_{\pm}}=\frac{1}{N}\sum_{j}e^{i(x_{j}\pm\theta_{i})}.$
(5)

The terms $x_{i}\pm\theta_{i}$ occur naturally so we define

$\displaystyle\xi_{i}=x_{i}+\theta_{i}$ (6)
$\displaystyle\eta_{i}=x_{i}-\theta_{i}$ (7)

And find

$\displaystyle\dot{\xi_{i}}$
$\displaystyle=J_{+}S_{+}\sin(\Phi_{+}-\xi)+J_{-}S_{-}\sin(\Phi_{-}-\eta)$ (8)
$\displaystyle\dot{\eta_{i}}$
$\displaystyle=J_{+}S_{+}\sin(\Phi_{+}-\xi)+J_{+}S_{-}\sin(\Phi_{-}-\eta)$ (9)

where $J_{\pm}=(J\pm K)/2$ 333Note we could set $J=1$ WLOG by rescaling time
but we decided against this so as to make the facilate a clean comparison to
the 2D model for which both $(J,K)$ appear.

We see ring swarmalators obey a sum of two Kuramoto models where the
traditional order parameter $Re^{i\Phi}:=(N)^{-1}\sum_{j}e^{i\theta_{j}}$ has
been replaced by a pair of new order parameters
$W_{\pm}=S_{\pm}e^{i\Phi_{\pm}}=(N)^{-1}\sum_{j}e^{i(x_{j}\pm\theta_{j})}$.
The new $W_{\pm}$ measure the system’s total amount of ‘space-phase order‘.

What kind of order is this? The limiting cases are trivial: static sync,
$(x_{i},\theta_{i})=(x^{*},\theta^{*})$, which produces maximal order
$S_{\pm}=1$ (which follows from substitution into Eq.(5)), and static async,
in which positions $x_{i}$ are fully uncorrelated with phases $\theta_{i}$,
which produces minimal order $S_{\pm}=0$. Between these two extremes, however,
lies something more interesting: Perfect correlation between position and
phase $x_{i}=\pm\theta_{i}+c$ for constant $c$ which yields either
$(S_{+},S_{-})=(0,1)$ or $(S_{+},S_{-})=(0,1)$ 444If the correlation is
$x_{i}+\theta_{i}=c$, then $S_{+}=1$. If $x_{i}-=\theta_{i}+c$, then
$S_{-}=1$.. What does $x_{i}\pm\theta_{i}$ mean physically? Picture
swarmalators on a ring as a group of fireflies flying around a circular track,
flashing periodically. Let $\theta_{i}=0$ be the beginning of their phase
cycle (when they flash). Then $x_{i}\pm\theta_{i}+c$ means each firefly
flashes at the same point on its circular lap $x_{i}=c$. Plotting swarmalators
as colored dots in space (where color represents phase) illustrates this
behavior most evocatively. Then $x_{i}=\pm\theta_{i}+c$ corresponds to a
colored splay state (Fig 2(c), Fig 3(c)). Since one of $S_{\pm}$ are maximal
in this rainbow-like state, we call $S_{\pm}$ the rainbow order parameters.

As we will show, $S_{\pm}$ are natural order parameters for the ring model
(insofar as they can distinguish between each of its collective states).

Figure 3: Collective states of ring model. Top row, swarmalators are drawn in
as colored points on a unit circle in a dummy $(\tilde{x},\tilde{y})$ plane to
allow a comparison with the 2D swarmalator model in Figure 2. $x$ corresponds
to the angular position on this unit circle (drawn black in (a) and (b) for
greater clarity. The point sizes in these panels, and the corresponding panels
(f),(g) are also larger, again to make things clearer) while the color
corresponds to the phase. Bottom row, swarmalators are drawn as points in the
$(x,\theta)$ plane. All results were found by integrating Eqs. (1), (2) using
an RK4 solver with timestep $dt=0.1$ for $T=500$ time units for $N=500$
swarmalators. Initial positions and phases were drawn from $[-\pi,\pi]$ in all
panels except panel (a), which were drawn from $[0,\pi]$ (since this choice of
initial conditions realized the static sync state). (a),(f) Static sync
$(J,K)=(1,1)$ (b),(g) Static $\pi$-state $(J,K)=(1,1)$ (c),(h) Static phase
wave $(J,K)=(1,-0.5)$ (d),(i) Active async $(J,K)=(1,-1.05)$. Here the
swarmalators jiggle about in $(x,\theta)$ as indicated by the double ended
arrows with no global space-phase order as indicated by the scatter plot in
(i). The amount of motion / jiggling depends on the population size $N$ as
discussed in the main text and illustrated in Figure 7. (e), (j) Static async
$(J,K)=(1,-2)$.

Motivation for model. Before showing our analysis, we quickly show how the
ring model connects to the original 2D model O’Keeffe et al. (2017). The 2D
model is

$\displaystyle\dot{\mathbf{x}}_{i}=\mathbf{v}_{i}+\frac{1}{N}\sum_{j=1}^{N}\Big{[}\mathbf{I}_{\mathrm{att}}(\mathbf{x}_{j}-\mathbf{x}_{i})F(\theta_{j}-\theta_{i})-\mathbf{I}_{\mathrm{rep}}(\mathbf{x}_{j}-\mathbf{x}_{i})\Big{]},$
(10)
$\displaystyle\dot{\theta_{i}}=\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}H_{\mathrm{\mathrm{att}}}(\theta_{j}-\theta_{i})G(\mathbf{x}_{j}-\mathbf{x}_{i})$
(11)

We pick out the rotational component of the swarming motion – since that part
is analogous to movement on a ring – by converting Eqs.(10), (11) to polar
coordinates. For certain choices of
$\mathbf{I}_{\mathrm{att}}(\mathbf{x}),F(\theta)$ etc this yields

$\displaystyle\dot{r_{i}}$
$\displaystyle=\tilde{\nu}(r_{i})+\frac{J}{2}\Bigg{[}\tilde{S}_{+}\cos\Big{(}\Phi_{+}-\xi_{i}\Big{)}+\tilde{S}_{-}\cos\Big{(}\Phi_{-}-\eta_{i}\Big{)}\Bigg{]}$
(12) $\displaystyle\dot{\xi_{i}}$
$\displaystyle=\tilde{\omega}(r_{i})+\tilde{J}_{+}(r_{i})\tilde{S}_{+}\sin\Big{(}\Psi_{+}-\xi_{i}\Big{)}$
$\displaystyle\qquad\qquad+\tilde{J}_{-}(r_{i})\tilde{S}_{-}\sin\Big{(}\Psi_{-}-\eta_{i}\Big{)}$
(13) $\displaystyle\dot{\eta_{i}}$
$\displaystyle=\tilde{\omega}(r_{i})+\tilde{J}_{-}(r_{i})\tilde{S}_{+}\sin\Big{(}\Psi_{+}-\xi_{i}\Big{)}$
$\displaystyle\qquad\qquad+\tilde{J}_{+}(r_{i})\tilde{S}_{-}\sin\Big{(}\Psi_{-}-\eta_{i}\Big{)}$
(14)

where we have switched to
$(\xi_{i},\eta_{i})=(\phi_{i}+\theta_{i},\phi_{i}-\theta_{i})$ coordinates and
$\phi_{i}$ is the spatial angle of the $i$-th swarmalator (which is analogous
to $x_{i}$ in the ring model). We put the derivation and definitions of the
various new quantities in the Appendix because they are cumbersome and
uninformative to display here.

Eqs. (13), (14) reveal the ring model hidden in the 2D model’s core. The
$(\dot{\xi}_{i},\dot{\eta}_{i})$ equations have the same form as the
$(\dot{\xi_{i}},\dot{\eta_{i}})$ equations in the 1D model (Eqs. (8), (9));
both are a summed pair of Kuramoto models. The only difference is that in the
2D model, $\tilde{\omega}_{i}(r_{i})$ and $\tilde{J}_{\pm}(r)$ depend on
$r_{i}$. So the ring model is like the 2D model with the radial dynamics
turned off. This is why we think studying the ring model will yield hints on
how to study the 2D model (one of the paper’s aims).

##  III Results

Figure 4: Rainbow order parameters $S_{\pm}$ for the ring swarmalator model
with $\nu_{i}=\omega_{i}=0$ and $J=1$. We have assumed WLOG that $S_{+}>S_{i}$
(i.e. if $S_{+}<S_{-}$ in simulations we swapped
$(S_{+},S_{-})\rightarrow(S_{-},S_{+})$). Data were collected by integrating
the model (1), (2) using an RK4 method with $(dt,T)=(0.1,500)$. The first
$90\%$ of data were discarded as transients and the mean of the remaining
$10\%$ were taken. We chose just $N=10$ swarmalators to illustrate the active
async state as clearly as possible; the fluctuations in $S_{\pm}$ that
characterize the state decay to $0$ for larger $N$ (see main text).

###  III.1 Numerics

Simulations show the system settles into five collective states depicted in
Figure 3. We provide the Mathematica code used for the simulations at
555https://github.com/Khev/swarmalators/blob/master/1D/on-
ring/regular/sandbox.nb.

  1. 1.

Static sync (Fig 3(a),(f)). Swarmalators fully synchronize their positions
$x_{i}=x^{*}$ and phases $\theta_{i}=\theta^{*}$ resulting in maximal space-
phase order $S_{\pm}=1$.

  2. 2.

Static $\pi$ state (Fig 3(b),(g)). One group of swarmalators synchronize at
$(x^{*},\theta^{*})$ and the remaining fraction synchronize $\pi$ units away
$(x^{*}+\pi,\theta^{*}+\pi)$. $S_{\pm}=1$ here also.

  3. 3.

Static phase wave (Fig 3(c),(h)). Swarmalators form a static splay state with
$x_{i}=2\pi i/N$ and $\theta_{i}=\pm x_{i}+c$ where the offset $c$ is
arbitrary and stems from the rotational symmetry in the model 666Linear
transformations $x_{i}\rightarrow x_{i}+C$ and $\theta_{i}\rightarrow C$ do
not change the dynamics because only differences $x_{j}-x_{i}$ and
$\theta_{j}-\theta_{i}$ appear in the model.. In $(\xi_{i},\eta_{i})$
coordinates, either $\xi_{i}$ is splayed $\xi_{i}=2\pi i/N$ and $\eta_{i}$ is
locked $\eta_{i}=c$ or vice versa. The order parameters are either
$(S_{+},S_{-})=(1,0)$ (where phase gradient of the rainbow is clockwise) or
$(S_{+},S_{-})=(1,0)$ (where the phase gradient of the rainbow is counter-
clockwise).

  4. 4.

Active async (Fig 3(d),(i)). Swarmalators form a dynamic steady state, moving
in clean limit cycles for small $N$, but in erratic, jiggling patterns for
large $N$. The motion cools and ultimately freezes as $N\rightarrow\infty$.
There is little space-phase order as indicated by the low, time averaged
values of $S_{\pm}$ (Fig 1) so we call this state ’active async’.

  5. 5.

Static async ((Fig 3(e),(j))). Swarmalators form a static, asynchronous
crystal with $S_{\pm}=0$.

Figure 4 shows the curve $S_{\pm}(K)$ can distinguish between all but the
static sync and static $\pi$ states.

Now we analyze the stability of the states. The static sync, $\pi$, and static
async states are analyzed using standard techniques, but are useful as warm
ups to the analysis of the much harder static phase wave state (which is the
main analytic contribution of the paper). The active async state, being non-
stationary, is analyzed mostly numerically.

###  III.2 Analysis

Static sync. We calculate the stability of the state by linearizing around the
fixed point in $(\xi,\eta)$ space. We seek the eigenvalues $\lambda$ of the
Jacobian $M$

$M=\left[\begin{array}[]{cc}Z_{\xi}&Z_{\eta}\\\ N_{\xi}&N_{\eta}\\\
\end{array}\right]$ (15)

where

$\displaystyle(Z_{\xi})_{ij}$
$\displaystyle=\frac{\partial\dot{\xi}_{i}}{\partial\xi_{j}}$ (16)
$\displaystyle(Z_{\eta})_{ij}$
$\displaystyle=\frac{\partial\dot{\xi}_{i}}{\partial\eta_{j}}$ (17)
$\displaystyle(N_{\xi})_{ij}$
$\displaystyle=\frac{\partial\dot{\eta}_{i}}{\partial\xi_{j}}$ (18)
$\displaystyle(N_{\eta})_{ij}$
$\displaystyle=\frac{\partial\dot{\eta}_{i}}{\partial\eta_{j}}$ (19)

Evaluating the derivatives in the above using Eqs. (8), (9) results in a clean
block structure:

$M=\left[\begin{array}[]{cc}J_{+}A(\xi)&J_{-}A(\eta)\\\
J_{-}A(\xi)&J_{+}A(\eta)\\\ \end{array}\right]$ (20)

where $A(y)_{i,i}=-\frac{1}{n}\sum_{j=1}^{n}\cos(y_{j}-y_{i})$ and
$A(y)_{i,j}=\frac{1}{n}\cos(y_{j}-y_{i})$ for a dummy variable $y$. Evaluated
at the fixed points of the static sync state, $\xi_{i}=c_{1}$ and
$\eta_{i}=c_{2}$ for constants $c_{1},c_{2}$, this becomes even simpler,

$M_{SS}=\left[\begin{array}[]{cc}J_{+}A_{0}&J_{-}A_{0}\\\
J_{-}A_{0}&J_{+}A_{0}\\\ \end{array}\right]$ (21)

where

$A_{0}:=\begin{bmatrix}-\frac{N-1}{N}&\frac{1}{N}&\dots&\frac{1}{N}\\\
\frac{1}{N}&-\frac{N-1}{N}&\dots&\frac{1}{N}\\\ {4}\\\
\frac{1}{N}&\frac{1}{N}&\dots&-\frac{N-1}{N}\end{bmatrix}$ (22)

Notice that Jacobian of the entire system $M$ has $\dim(N)=2N$ since there are
two state variables $(x,\theta)$ for each of the $N$ swarmalators, but that
$\dim(A_{0})=N$ since it is a subblock of $M$. Now, the eigenvalues
$\hat{\lambda}$ of $A_{0}$ are well known: there is $1$ with value
$\hat{\lambda}=0$ stemming from the rotational symmetry of the model and $N-1$
$\hat{\lambda}=-1$. We use these to find the desired eigenvalues $\lambda$ of
$M_{SS}$ using the following identity for symmetric block matrices

$\det E:=\det\left[\begin{array}[]{cc}C&D\\\ D&C\\\
\end{array}\right]=\det(C+D)\det(C-D)$ (23)

This implies the eigenvalues of $E$ are the union of the eigenvalues of $C+D$
and $C-D$. Applying this identity to $M_{SS}$ (which has the required
symmetric structure) yields

$\displaystyle\lambda_{0}$ $\displaystyle=0$ (24) $\displaystyle\lambda_{1}$
$\displaystyle=-J$ (25) $\displaystyle\lambda_{2}$ $\displaystyle=-K$ (26)

with multiplicities $2,-1+N,-1+N$ (which sum to the required $2N$). This tells
us the static sync is stable for $J>0$ and $K>0$ consitent with simulations.

Static $\pi$-state. The fixed points here are
$(x_{i},\theta_{i})=(c_{1},c_{2})$ for $i=1,2,\dots,N/2$ and
$(x_{i},\theta_{i})=(c_{1}+\pi,c_{2}+\pi)$ for $i=N/2,\dots,N$. Conveniently,
the shift in $\pi$ for exactly half of swarmalators does not change the form
of the Jacobian, $M_{\pi}=M_{SS}$, so the stability is the same as before.
This means the static sync and $\pi$ states are bistable for all $J>0,K>0$.
Appendix A discusses the basins of attraction for each state.

Static phase wave. We calculate the stability of the static phase wave using
the same strategy as before: linearize around the fixed points and exploit the
block structure of the Jacobian $M$. This time, however, the calculations are
tougher.

The fixed points of the static phase wave take two forms: Either $\xi$ is
splayed $\xi_{i}=2\pi(i-1)/N+c_{1}$ and $\eta_{i}$ is synchronized
$\eta_{i}=c_{2}$ (clockwise rainbow) or $\xi$ is synchronized $\xi_{i}=c_{1}$
and $\eta$ splayed $\eta_{i}=2\pi(i-1)/N+c_{2}$ (counter clockwise rainbow).
Here $i=1,\dots,N$ and the constants $c_{1},c_{2}$ are offsets. WLOG we
analyze the fixed point with $\xi$ splayed and $\eta$ sync’d. The Jacobian is

$M_{SPW}=\left[\begin{array}[]{cc}J_{+}A_{1}&J_{-}A_{0}\\\
J_{-}A_{1}&J_{+}A_{0}\\\ \end{array}\right]$ (27)

where $A_{0}$ is as before Eq.(22) but $A_{1}$ is new:

$\displaystyle(A_{1})_{ii}$ $\displaystyle=-\frac{1}{N}\sum_{j\neq
i}\cos\frac{2\pi}{N}(j-i):=-\frac{1}{N}\sum_{j\neq i}c_{ij}$ (28)
$\displaystyle(A_{1})_{ij}$
$\displaystyle=\frac{1}{N}\cos\frac{2\pi}{N}(j-i):=\frac{1}{N}c_{ij}$ (29)

The following notation will be useful

$\displaystyle c_{ij}$ $\displaystyle:=\cos 2\pi(i-j)/N$ (30) $\displaystyle
s_{ij}$ $\displaystyle:=\sin 2\pi(i-j)/N$ (31) $\displaystyle\beta_{ij}$
$\displaystyle:=c_{ij}+\mathcal{I}s_{ij}$ (32) $\displaystyle c_{i}$
$\displaystyle:=\cos 2\pi i/N$ (33) $\displaystyle s_{i}$ $\displaystyle:=\sin
2\pi i/N$ (34)

where $\mathcal{I}=\sqrt{-1}$ is the imaginary unit 777We choose the non-
standard notation $\mathcal{I}$ since we have already used $i$ and $j$ as
indices., $\beta_{0,0}:=\beta=e^{2\pi\mathcal{I}N}$ is the primitive root of
unity, and $\beta_{k}=\beta^{k}$ is the $k$-th root of unity. The diagonal
element $A_{ii}$ may be simplified. Recalling the sum of roots of unities are
zero (you can think of $\beta_{k}$ as a vector pointing to the beginning of
the $k$-th segment of size $1/N$ of the unit circle; then summing all the
vectors around the unit circle results in zero)

$\sum_{k=0}^{N-1}\beta_{k}=\sum_{k=0}^{N-1}c_{k}+\mathcal{I}s_{k}=0$ (35)

which implies

$\displaystyle\sum_{k=0}^{N-1}c_{k}=0$ (36)
$\displaystyle\sum_{k=0}^{N-1}s_{k}=0$ (37)

(In other words, the discrete sum of $\cos(2\pi k/N)$ and $\sin(2\pi kN)$
around the unit circle is zero, which can be seen by symmetry). Applying these
identities to Eq. (28) for $(A_{1})_{ii}$,

$\displaystyle(A_{1})_{ii}$ $\displaystyle=\frac{1}{N}\sum_{j\neq
i}\cos\frac{2\pi}{N}(j-i)=-\frac{1}{N}\sum_{k=1}^{N-1}c_{k}$ (38)
$\displaystyle=-\frac{1}{N}(-c_{0}+\frac{1}{N}\sum_{k=0}^{N-1}c_{k})=-\frac{1}{N}(-c_{0}+0)$
(39) $\displaystyle=\frac{c_{0}}{N}$ (40)

Notice in the first sum over $i$ the $j-th$ term is excluded which means the
second sum over $k$ begins at $1$. Since $c_{0}=c_{0,0}=1$, we get

$\displaystyle(A_{1})_{ij}=\frac{1}{N}c_{i,j}$ (41)

With $A_{0}$ (Eq. (22)) and $A_{1}$ (Eq. (41)) in hand, we can begin finding
the eigenvalues of $M_{SPW}$ Eq. (27) by using another identity for block
matrices. If the consituent matrices $A,C$ of a general block matrix $P$,

$P=\left[\begin{array}[]{cc}A&B\\\ C&D\\\ \end{array}\right]$ (42)

commute then

$\det(P)=\det(AD-BC)$ (43)

Luckily, the constituent matrices $A_{0}$ and $A_{1}$ of $M_{SPW}$ do commute.
So we apply the above identity to the characteristic equation for $\lambda$

$\det(M_{SPW}-\lambda I)=\det\left[\begin{array}[]{cc}J_{+}A_{1}-\lambda
I&J_{-}A_{0}\\\ J_{-}A_{1}&J_{+}A_{0}-\lambda I\\\ \end{array}\right]=0$ (44)

And find

$\displaystyle\det(M_{SPW}-\lambda I)$
$\displaystyle=\det\Big{[}(J_{+}A_{1}-\lambda I)(J_{+}A_{0}-\lambda I)$ (45)
$\displaystyle\hskip 28.45274pt-J_{-}^{2}A_{1}A_{0}\Big{]}$
$\displaystyle:=\det(G)$ (46)

where to compactify the RHS we have defined

$G:=\begin{bmatrix}g_{0}&g_{1}&\dots&g_{N-1}\\\ g_{N-1}&g_{0}&\dots&g_{N-2}\\\
{4}\\\ g_{1}&g_{2}&\dots&g_{0}\end{bmatrix}$ (47)

where

$\displaystyle
g_{0}=\lambda^{2}+\frac{J_{+}}{N}\lambda(N-1-c_{0})+\frac{J_{-}^{2}-J_{+}^{2}}{N}c_{0}$
(48) $\displaystyle
g_{(k>0)}=-\frac{J_{+}}{N}\lambda(1+c_{k})+\frac{J_{-}^{2}-J_{+}^{2}}{N}c_{k}$
(49)

We want the determinent of $G$ which the product of its eigenvalues
$\hat{\lambda}_{j}$ (not to be confused $M_{SPW}$’s eigenvalues $\lambda_{i}$
which we are trying to find),

$\displaystyle\det(G)$ $\displaystyle=\prod_{j=0}^{N-1}\hat{\lambda}_{j}=0$
(50)

Casting an eye back to Eq. (47), we see $G$ is a circulent matrix so its
eigenvalues $\hat{\lambda}_{j}$ are known exactly! –

$\displaystyle\hat{\lambda}_{j}$
$\displaystyle=\sum_{k=0}^{N-1}g_{k}\beta_{j*k}$ (51)
$\displaystyle\hat{\lambda}_{j}$ $\displaystyle=\lambda^{2}+\lambda
J_{+}-\frac{J_{+}}{N}\lambda\sum_{k=0}^{N-1}(1+c_{k})\beta_{j*k}$
$\displaystyle\hskip
56.9055pt+\frac{J_{-}^{2}-J_{+}^{2}}{N}\sum_{k=0}^{N-1}c_{k}\beta_{j*k}$ (52)

where $\beta_{k}$ is the primitive root of unity as before. Note we mean the
product $j*k$ in $\beta_{j*k}=\cos 2\pi j*k/N+\mathcal{I}\sin 2\pi j*k/N$.
Note also that while $G$ has $N$ eigenvalues $\hat{\lambda}$ (running from
$j=0,\dots N-1$), $M_{SPW}$ still has the required $2N$ eigenvalues $\lambda$
since each $\hat{\lambda_{j}}$ is quadratic in $\lambda$.

One last push remains: we simplify the summands in Eq. (52) using some basic
trig identities, find each $\hat{\lambda}_{j}$ for $j=0,1,\dots N-1$, and then
set $\hat{\lambda}_{j}=0$ (since $\det(G)=0$ so each term in the product must
be zero) which yields a quadratic equation for our target $\lambda_{j}$.

We begin with $j=0$ since it is distinguished from the other values $j$ takes.
The first summand on the RHS of Eq.(52) becomes
$\sum_{k=0}^{N-1}(1+c_{k})\beta_{0}=\sum_{k=0}^{N-1}(1+c_{k})(1)=N+0=N$. The
second summand becomes
$\sum_{k=0}^{N-1}c_{k}\beta_{j*k}=\sum_{k=0}^{N-1}c_{k}(1)=0$ (which follows
from from Eq. (36)). Plugging these in yields

$\displaystyle\hat{\lambda}_{0}$ $\displaystyle=\lambda^{2}+\lambda
J_{+}-\frac{J_{+}}{N}\lambda(N)=\lambda^{2}$ (53)

Setting $\hat{\lambda}_{0}=0$ gives our first eigenvalue

$\displaystyle\lambda_{0}=0$ (54)

with multiplicity $2$.

Next we analyze $\hat{\lambda}_{j>0}$ all at once. First note

$\displaystyle\sum_{k=0}^{N-1}\beta_{j*k}$
$\displaystyle=\sum_{k=0}^{N-1}\cos\frac{2\pi jk}{N}+\mathcal{I}\sin\frac{2\pi
jk}{N}$ (55) $\displaystyle=0+0\times\mathcal{I}=0$ (56)

for all $j>0$ (when $j=0$ we get the simple sum $N-1$ which is why we
considered the case $j=0$ separately). This implies the second summand in Eq.
(52) simplifies into the third summand
$\sum_{k=0}^{N-1}(1+c_{k})\beta_{j*k}=\sum_{k=0}^{N-1}c_{k}\beta_{j*k}$ when
$j>0$ which in turn becomes

$\displaystyle\sum_{k=0}^{N-1}c_{k}\beta_{j*k}=\sum_{k=0}^{N-1}c_{k}c_{j*k}+\mathcal{I}c_{k}s_{j*k}$
(57)

Using the standard formulas,

$\displaystyle\cos(a)\cos(b)$
$\displaystyle=\frac{1}{2}\Big{(}\cos(a+b)+\cos(a-b)\Big{)}$ (58)
$\displaystyle\cos(a)\sin(b)$
$\displaystyle=\frac{1}{2}\Big{(}\sin(a+b)-\sin(a-b)\Big{)}$ (59)

we see

$\displaystyle\sum_{k=0}^{N-1}c_{k}c_{j*k}$
$\displaystyle=\frac{1}{2}\sum_{k=0}^{N-1}\Big{(}\cos\frac{2\pi(1+j)k}{N}+\cos\frac{2\pi(1-j)k}{N}\Big{)}$
(60) $\displaystyle\sum_{k=0}^{N-1}c_{k}s_{j*k}$
$\displaystyle=\frac{1}{2}\sum_{k=0}^{N-1}\Big{(}\sin\frac{2\pi(1+j)k}{N}+\sin\frac{2\pi(1-j)k}{N}\Big{)}$
(61)

where we have inverted the $c_{j}$ notation for clarity. Both the cosine terms
in Eq. (60) and sine terms in Eq. (61) are zero (as per Eq.(56)) except when
the arguments are $0$. This occurs for $j=1,-1$ (where $j=-1$ is interpreted
modulo $N$ and equivalent to $j=N-1$. Since $\cos 0=1$ and $\sin=0$, this
yields

$\displaystyle\sum_{k=0}^{N-1}c_{k}c_{jk}$
$\displaystyle=\frac{N}{2}(\delta_{j,1}+\delta_{j,N-1})$ (62)
$\displaystyle\sum_{k=0}^{N-1}c_{k}s_{jk}$ $\displaystyle=0$ (63)

where $\delta_{i,j}=1$ is the Kronecker delta. Applying the above to Eq. (57)
gives

$\displaystyle\sum_{k=0}^{N-1}c_{k}\beta_{j*k}=\frac{N}{2}(\delta_{j,1}+\delta_{j,N-1})$
(64)

Plugging this into Eq. (52) for $\hat{\lambda}_{j}$

$\displaystyle\hat{\lambda}_{j}=\lambda^{2}+\lambda
J_{+}+\Big{(}-\frac{J_{+}}{N}\lambda+\frac{J_{-}^{2}-J_{+}^{2}}{N}\Big{)}(\delta_{j,1}+\delta_{j,N-1})=0$
(65)

which holds for $j=1,\dots,N-1$ since, recall, we have analyzed $j=0$
separately and in that case found $\lambda=0$ with multiplicity $2$ (Eq.
(54)). The desired $\lambda$ are the roots of the above equations. Both their
values and multiplicities depend on $N$. For $N\geq 4$, a general pattern
holds, but $N=2,3$ are special cases.

$N=2$ is special because both of the Kronecker delta functions trigger at the
same time: $\delta_{j,N-1}=\delta_{j,2-1}=\delta_{j,1}=1$. Eq. (65) then
becomes

$\displaystyle\hat{\lambda}_{j}=\lambda^{2}+\lambda
J_{+}+2\Big{(}-\frac{J_{+}}{N}\lambda+\frac{J_{-}^{2}-J_{+}^{2}}{N}\Big{)}=0$
(66) $\displaystyle\Rightarrow\lambda=\pm\sqrt{J_{+}{-}^{2}-J_{-}^{2}}$ (67)

Table 1 row $1$ reports these eigenvalues along with the two zero eigenvalues
$\lambda=0$ which gives the required $2N=4$ eigenvalues total.

$N=3$ is special because precisely one of the two kronecker $\delta_{i,j}$
functions is on for all $j$; there is no $j$ for which both
$\delta_{j,1},\delta_{j,N-1}$ are zero simultaneously. For
$j=1,\delta_{j,1}=1$ and for $j=N-1,\delta_{j,N-1}=\delta_{j,2}=1$ (recall $j$
from from $1$ to $N-1$; so when $N=2,j=1,2$ only). So we get the following for
both $j=1,2$

$\displaystyle\hat{\lambda}_{j}=\lambda^{2}+\lambda
J_{+}+\Big{(}-\frac{J_{+}}{N}\lambda+\frac{J_{-}^{2}-J_{+}^{2}}{N}\Big{)}=0$
(68)
$\displaystyle\Rightarrow\lambda_{\pm}=\frac{1}{4}\left(-J_{+}\pm\sqrt{9J_{+}^{2}-8J_{-}^{2}}\right)$
(69)

and so the pair $\lambda_{\pm}$ has multiplicity two as shown Table 1 row 2
(i.e. the full set is $\lambda_{+},\lambda_{+},\lambda_{-},\lambda_{-}$).

Finally, for $N\geq 4$, there are now intermediary values of $j=2,\dots N-2$
for which $\delta_{j,1}=\delta_{j,N-1}=0$ simultaneously in which case

$\displaystyle\hat{\lambda}_{j}=\lambda(\lambda+J_{+})=0$ (70)
$\displaystyle\Rightarrow\lambda=0,-J_{+}$ (71)

Table 1 row $3$ summarizes the $\lambda$’s and their multiplicities in this
general case which finally completes our calculation. We have confirmed the
expressions for $\lambda$ are correct by using Mathematica (notebook provided
at 888https://github.com/Khev/swarmalators/blob/master/1D/on-
ring/regular/stability-static-phase-wave.nb) to compute the eigenvalues of
$M_{SPW}$ for $N=2,\dots 6$ (for larger $N$, Mathematica starts to struggle).

$N$ | (Multiplicity, eigenvalue)
---|---
$2$ | $\lambda=(2,0),(1,\pm\sqrt{J_{+}^{2}-J_{-}^{2}})$
$3$ | $\lambda=(2,0),(2,\left(-J_{+}\pm\sqrt{9J_{+}^{2}-8J_{-}^{2}}\right)/4)$
$4$ | $\lambda=(3,0),(1,-J_{+}),(2,\left(-J_{+}\pm\sqrt{9J_{+}^{2}-8J_{-}^{2}}\right)/4)$
$N\geq 4$ | $\lambda=(0,N-1),(N-3,-J_{+}),(2,\left(-J_{+}\pm\sqrt{9J_{+}^{2}-8J_{-}^{2}}\right)/4)$
Table 1: Spectrum of static phase wave for different population sizes. The
cases $N=2,3$ are distinguished. For $N\geq 4$ the pattern is fixed. Notice a
multiplicity $1$ for eigenvalues that are complex conjugates $\lambda_{\pm}$
means one instance of the pair; a multiplicity of $2$ means two instances of
the pair: $\lambda_{+},\lambda_{+},\lambda_{-},\lambda_{-}$.

Now we use the expressions for $\lambda$ to predict the bifurcations of the
static phase wave. For convenience we write out the expressions below (which
are valid for $N>2$. For $N=2$, $\lambda=\sqrt{J_{-}^{2}-J_{+}^{2}}>0$ for
$K<0$ meaning the static phase wave is unstable in this case. We will analyze
the $N=2$ case fully in the next section).

$\displaystyle\lambda_{0}$ $\displaystyle=0$ (72) $\displaystyle\lambda_{1}$
$\displaystyle=-J_{+}$ (73) $\displaystyle\lambda_{2}$
$\displaystyle=\frac{1}{4}\left(-J_{+}\pm\sqrt{9J_{+}^{2}-8J_{-}^{2}}\right)$
(74)

$\lambda_{2}$ exists for all $N\geq 3$ and triggers a Hopf bifurcation at
$K_{c}=-J$ which comes from solving $Re(\lambda_{2}^{\pm})=0$ (the determinant
$9J_{+}^{2}-8J_{-}^{2}$ is negative here so $\lambda_{2}^{\pm}$ are complex
conjugates). It also triggers a saddle node bifurcation at $K=0$ which comes
from solving $\lambda_{2}^{+}=0$ ($\lambda_{2}^{-}$ is never $0$).
$\lambda_{1}$ exists for all $N\geq 4$ and triggers a saddle node bifurcation
also at $K_{c}=-J$. These results imply the static phase wave is stable when

$-J<K<0$ (75)

for all population sizes $N>2$. To recap, at the left boundary $K=-J$, it
destabilizes via a Hopf bifurcation for $N=3$ and a simultaneous Hopf and
saddle node for $N\geq 4$. At the right boundary $K=0$, it destabilizes via a
saddle node for all $N\geq 3$.

Figure 5: Dynamics in active async state. Top: Position $x(t)$ for a typical
swarmalator for $N=2,5,50$ swarmalators showing oscillatory behavior. For
small $N$, the oscillations are smooth and fast. For larger $N$, they are
irregular and slow. The dynamics of the swarmalator‘s phase $\theta_{i}$
(unplotted) are similar. Bottom: time series of order parameter $S_{+}$. The
behavior of $S_{-}$ is qualitatively the same. Simulation parameters were
$(dt,T)=(0.1,100)$ and $(J,K)=(1,-1.05)$.

Analysis of active async. The destabilization of the static phase wave via a
Hopf bifurcation for $N\geq 2$ implies a limit cycle is born at $K<K_{c}$.
Here swarmalators oscillate about their mean position in space and phase (Fig
3(d)) with little overall order $S_{\pm}\approx 0$ (Fig 4) for most $N$.
Figure 5 depicts the oscillations in $x(t)$ for a typical swarmalator for
different population sizes. For small $N$, the oscillations are smooth and
fast. For larger $N$, the oscillations are irregular and slow. (The
oscillations in phase, unplotted, behave the same way). The bottom panel shows
the dynamics of the rainbow order parameters $S_{\pm}$ are similar.

Figure 6: Phase space for $N=2$ swarmalators in $(\Delta x,\Delta\theta)$
coordinates (Eqs. (76), (77)) for $(J,K)=(1,-1)$. Saddle point are shown as
open red circles, nonlinear centers as open green circles. The family of
periodic orbits, which correspond to the active async state, are given by Eq.
(80) for $C=-6,\dots,6$.

In the simple case $N=2$ the active state can be studied using the standard
transformation to mean and difference coordinates $(\langle y\rangle,\Delta
y)=(y_{1}+y_{2},y_{2}-y_{1})$ where $y$ is a dummy variable. The mean
positions and phases are invariant $\dot{\langle
x\rangle}=\dot{\langle\theta\rangle}=0$ because of the pairwise oddness of
Eqs. (1), (2) which implies $x_{1}=-x_{2}$ and $\theta_{1}=-\theta_{2}$. The
differences, however, evolve according to

$\displaystyle\dot{\Delta x}$ $\displaystyle=J\sin\Delta x\cos\Delta\theta$
(76) $\displaystyle\dot{\Delta\theta}$
$\displaystyle=K\sin\Delta\theta\cos\Delta x$ (77)

These ODEs have 8 fixed points. For $K>0$, the two fixed points $(\Delta
x,\Delta\theta)=(0,\pi),(0,\pi)$ are stable nodes (which corresponds to static
synchrony) while the rest are saddles. For $K<0$, four nonlinear centers exist
while the other fixed points are again saddles. The family of periodic
solution surrounding the centers correspond to the active async state and can
be found explicitly by dividing Eq. (76) by Eq. (77) and integrating,

$\displaystyle\frac{d{\Delta x}}{d{\Delta\theta}}=\frac{J\sin\Delta
x\cos\Delta\theta}{K\sin\Delta\theta\cos\Delta x}$ (78) $\displaystyle
K\int\frac{\cos x}{\sin x}dx=J\int\frac{\cos\theta}{\sin\theta}d\theta$ (79)
$\displaystyle\sin x^{K}=C\sin\theta^{J}$ (80)

For some constant $C$ determined by initial conditions. Figure 6 plots the
contours of Eq. (80), which correspond to the limit cycle of the active async
state, along with the fixed point for $K=-1$. These are consistent with the
oscillations in $x(t)$ in Figure 5. They also tell us that the active async is
stable for all $K<0$ when $N=2$, which confirms our result from the last
section that the static phase wave is unstable when $N=2$.

Figure 7: Mean velocity $\langle v\rangle=N^{-1}\sum_{i=1}^{N}v_{i}$. where
$v_{i}=(v_{x}^{2}+v_{\theta}^{2})^{1/2}$ is the velocity of the $i$-th
swarmalator, in the $(N,K)$ plane for $J=1$ and $\nu_{i}=\omega_{i}=0$.
Simulation parameters were $(dt,T)=(0.1,5000)$. The first $90\%$ of data were
discarded as transients and the mean of the remaining $10\%$ are plotted.

For $N>2$, analysis of the active async state becomes unwieldy so we
numerically examine the state by plotting the mean velocity or activity
$\langle v\rangle=N^{-1}\sum_{i=1}^{N}v_{i}$ where
$v_{i}=(v_{x}^{2}+v_{\theta}^{2})^{1/2}$ is the velocity of the $i$-th
swarmalator. Since all of the other collective states are stationary, $\langle
v\rangle$ is a natural order parameter for the active async state. Figure 7,
in which $J=1$, shows the $\langle v\rangle$ persists for large $N$ when
$K\approx K_{c}=-1$ (recall, we know the state must exist for some $K$ for all
finite $N$ via the Hopf bifurcation). For small coupling strengths
$K\approx-2$, however, the motion dies out for moderately large systems
$N\gtrapprox 5$ and the final static async state is born.

Figure 8: Pseudo-bifurcation diagram for model (1), (2) in $(K,N)$ space with
$J=1$ (pseudo because the population size $N$, a discrete quantitiy, is
plotted on a continous axis; see main text). The $N\rightarrow\infty$ limit is
represented by plotting $...$ on the $y$-axis. The dotted line schematically
represents the bifurcation curve of the active async state and thus merges
with the line $K=J=-1$ as $N\rightarrow\infty$.

Analysis of static async. Here the swarmalators sit at fixed points
$(x_{i}^{*},\theta_{i}^{*})$ scattered uniformly in space and phase (Fig 3(e))
implying no global order $S_{\pm}=0$ (Fig 4). This state is hard to analyze
for finite $N$ because for most population sizes $N\geq 3$, multiple
configurations of fixed points $(x_{i}^{*},\theta_{i}^{*})$ exist for a given
set of parameters $(J,K)$. Simply enumerating this family of fixed points is a
hard problem – it has not been done in the regular Kuramoto model – never mind
analyzing their stability. In the continuum limit $N\rightarrow\infty$,
however, the stability may be analyzed since the state has a simple
representation: $\rho(x,\theta,t)=1/(4\pi^{2})$ where
$\rho(x,\theta,t)dxd\theta$ gives the fraction of swarmalators with positions
between $x+dx$ and phases between $\theta+d\theta$ at time $t$. The density
obeys the continuity equation

$\dot{\rho}+\nabla.(v\rho)=0$ (81)

And the velocity $v=(v_{x},v_{\theta})$ is given by the $N\rightarrow\infty$
limit of Eqs (1), (2)

$\displaystyle v_{x}$
$\displaystyle=J\int\sin(x^{\prime}-x)\cos(\theta^{\prime}-\theta)\rho(x^{\prime},\theta^{\prime},t)dx^{\prime}d\theta^{\prime}$
(82) $\displaystyle v_{\theta}$
$\displaystyle=K\int\sin(\theta^{\prime}-\theta)\cos(x^{\prime}-x)\rho(x^{\prime},\theta^{\prime},t)dx^{\prime}d\theta^{\prime}$
(83)

The stability is analyzed by plugging a general perturbation

$\rho=\rho_{0}+\epsilon\eta=(4\pi^{2})^{-1}+\epsilon\eta(x,\theta,t)$ (84)

into the continuity equation and computing the spectrum. Since such analyses
are standard, we put the details in Appendix C. The eventual result is

$K<K_{c}=-J$ (85)

Recall this is only valid as $N\rightarrow\infty$. Interestingly, this result
also proves the active async state disappears for $N\rightarrow\infty$ via a
‘squeeze’ argument: since the static phase wave also loses stability at
$K_{c}=-J$, there is no room left for the active async state to exist as
$N\rightarrow\infty$.

Figure 9: Bifurcation diagram in $(J,K)$ space as $N\rightarrow\infty$.

This completes our analysis. Figure 8 summarizes our findings in a pseudo
bifurcation diagram in $(K,N)$ space for $J=1$(pseudo because $N$ is a
discrete quantity and so should plot by plotted on a continuous axis, but we
wanted to show the active async state dependence on $N$). The boundary of the
active async state for $N\geq 3$ is found numerically by finding the first
time $\langle v\rangle=0$ which indicates the swarmalators are no longer
active. Figure 9 reports the bifurcation diagram in $(J,K)$ space when
$N\rightarrow\infty$.

##  IV Match to real world swarmalators

Swarmalators are defined as entities with a two-way interaction between
swarming and synchronization O’Keeffe et al. (2017). Below we list examples
which (i) appear to meet this definition 999We say ‘appear’ to meet these
definition because what is means to prove a bidirectional space-phase coupling
in an experimental system is somewhat ambiguous. Most experimental studies of
swarmalators either infer a space-phase coupling exists based on the
observations (like the microswimmers Quillen et al. (2021a); Yang et al.
(2008); Yuan et al. (2014) we describe; Yuan et al. (2014) do a particularly
comprehensive study), but do not specify if the coupling bidirectional. Other
studies demonstrate a two-way coupling in the sense that they encoding it in a
model that produces behavior similar to observed data (like magnetic domain
walls Hrabec et al. (2018) or myxobacteria Igoshin et al. (2001)). In other
words, we’re trying to be conservative in any claims we make., (ii) swarm in
ring-like geometries, and (iii) display collective behavior similar to the
ring model.

Sperm are a classic microswimmer which swarm in solution and sync their tail
gaits Yang et al. (2008). Sperm collected from ram semen and confined to 1D
rings bifurcate from an isotropic state (analogues to the static async state)
to a vortex state in which sperm rotate either clockwise or anticlockwise
Creppy et al. (2016) which implies their positions and orientations are
splayed just like the static phase wave (recall the static phase wave is
really a state of uniform rotation and only static in the frame co-moving with
the natural frequencies $\omega,\nu$) Creppy et al. (2016). Moreover, the
transition has associated transient decay of rotation velocity (see Fig 4a in
Creppy et al. (2016)) consistent with a Hopf bifurcation, just like the ring
model (74). Note unlike other studies of syncing sperm Yang et al. (2008),
here the phase variable is the sperm’s orientation, not its tail rhythm.

Vinegar eels are a type of nematode found in beer mats and the slime from tree
wounds Quillen et al. (2021a, b). They can be considered swarmalators because
they sync the wriggling of their heads, swarm in solution, and it seems likely
based on their behavior said sync and swarming interact Quillen et al. (2021a,
b) (neighbouring eels sync more easily than distant eels, so sync interacts
with swarming, and sync’d eels presumable affect each local hydrodynamic
environment and thereby affect each other’s movements, so swarming interacts
with sync). When confined to 2D disks, they seek out the 1D ring boundary
forming metachronal waves in which the phase of their gait and their spatial
positions around the ring are splayed similar to the static phase wave Quillen
et al. (2021a, b) (although note the winding number for the metachronal waves
is $k>1$; a full rotation in physical space $x$ produces $k>1$ rotations in
phase $\theta$).

C. elegans are another type of microswimmer which also swarm and sync the gait
of their tails. When confined to 1D channels they form synchronous clusters
analogous to the static sync state Yuan et al. (2014). Though not strictly
consistent with the ring model, a channel being a 1D line as opposed to a 1D
ring, we mention them here because it is natural to expect sync clusters would
persist in a 1D ring too.

##  V Discussion

Aim one of the paper was to model real-world swarmalators swarming in 1D. This
was a success. The model captured the behavior of vinegar eels, sperm, and C.
elegans. Even so, the model failed to capture the phenomenology of other 1D
swarmalators such as the the two-cluster states of bordertaxic Japanese tree
frogs Aihara et al. (2014) or the cluster dynamics of synthetic microswimmers
Guzmán-Lastra et al. (2016). Crafting a model that mimics these systems is a
challenge for future research.

Aim two was to provide a stepping stone – i.e. the ring model – to an analytic
understanding of the 2D swarmalator model’s phenomenology, in particular its
bifurcations. To this end, we studied the destabilization of the static async
state, deriving a 1D version of $K_{1}$. But an analogue of the active phase
wave state did not appear in the ring model 101010We note the splintered phase
wave and active phase waves are realized in a 1D ring model presented in the
Supplementary Information of O’Keeffe et al. (2017). But that model does not
have the clean format of the current model Eqs.(8),(9) which, being so simple,
seem the natural model to study 1D swarmalator phenomena. Moreover, the
analogues of the splintered phase wave and active phase wave states there are
unsteady, the order parameters varying in time. So in that sense that are not
simpler warm-up versions of the states in 2D., so we could not study the
second bifurcation at $K_{2}$ beyond which $S_{+}$ declines (Figure 1). On the
upside, a 1D analogue of the static phase wave was observed, and were able to
specify its stability (in the 2D model this state was only observed for $K=0$
and its stability was not calculated). Moreover, this stability calculation
certified the existence of the new active async state (we say certify because
we originally thought the state was just a long transient or perhaps a
numerical artefact; our results prove it exists for all finite $N$) in which
the swarmalators execute noisy, Brownian-like motion. What’s interesting here
is that this motion occurs for identical, noise-free swarmalators. This
suggest any jittery behavior observed in real-world swarmalator systems may
arise purely from the cross-talk between units’ tendency to sync and swarm,
and not from thermal agitation or other forms of noise.

Taken together, our analytic findings take us one step closer to understanding
the bifurcations of the swarmalator model which for now remain unexplained.

We would love to see future work analyze the ring model using OA theory Ott
and Antonsen (2008); Pikovsky and Rosenblum (2015); Engelbrecht and Mirollo
(2020); Chen et al. (2017). A breakthrough from Ott and Antonsen, the theory
states that the density of the infinite-$N$ Kuramoto model $\rho(\theta,t)$
has an invariant manifold of poisson kernels which allows dynamics for the
classic sync order parameter $Z:=N^{-1}\sum_{j}e^{i\theta_{j}}$ to be derived
explicitly. This is a big win. It reduces an $N>>1$-dimensional nonlinear
system to a simple 2D system (one ODE for the complex quantity $Z$) – a
drastic simplification which effectively solves the Kuramoto model. Given the
ring model in $(\xi,\eta)$ coordinates (Eqs.(8),(9)) resembles the Kuramoto
model so closely, we suspect it may be solved with OA theory’s magic: If
regular oscillators are defined on the unit circle $\theta_{i}\in S_{1}$ and
have an invariant manifold of poisson kernels, could ring swarmalators,
defined on the unit torus $(x,\theta)\in(S^{1},S^{1})$, have an invariant
manifold of some ‘toroidal’ poisson kernel? If so, explicit dynamics for the
rainbow order parameters $W_{\pm}$ may be derivable and in that sense the ring
model solved exactly.

## References

  * Winfree (2001) A. T. Winfree, _The geometry of biological time_ , vol. 12 (Springer Science & Business Media, 2001).
  * Kuramoto (2003) Y. Kuramoto, _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * Pikovsky et al. (2003) A. Pikovsky, J. Kurths, M. Rosenblum, and J. Kurths, _Synchronization: a universal concept in nonlinear sciences_ , 12 (Cambridge university press, 2003).
  * Toner and Tu (1995) J. Toner and Y. Tu, Physical review letters 75, 4326 (1995).
  * Vicsek et al. (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Physical review letters 75, 1226 (1995).
  * Jiang and McCall (1993) Z. Jiang and M. McCall, JOSA B 10, 155 (1993).
  * Kozyreff et al. (2000) G. Kozyreff, A. Vladimirov, and P. Mandel, Physical Review Letters 85, 3809 (2000).
  * Peskin (1975) C. S. Peskin, _Mathematical Aspects of Heart Physiology_ (Courant Institute of Mathematical Sciences, New York, 1975), pp. 268–278.
  * 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).
  * Katz et al. (2011) Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, Proceedings of the National Academy of Sciences 108, 18720 (2011).
  * Yang et al. (2008) Y. Yang, J. Elgeti, and G. Gompper, Physical review E 78, 061903 (2008).
  * Riedel et al. (2005) I. H. Riedel, K. Kruse, and J. Howard, Science 309, 300 (2005).
  * Quillen et al. (2021a) A. Quillen, A. Peshkov, E. Wright, and S. McGaffigan, arXiv preprint arXiv:2101.06809 (2021a).
  * Quillen et al. (2021b) A. Quillen, A. Peshkov, E. Wright, and S. McGaffigan, arXiv preprint arXiv:2104.10316 (2021b).
  * Taylor (1951) G. I. Taylor, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 209, 447 (1951).
  * Tamm et al. (1975) S. L. Tamm, T. Sonneborn, and R. V. Dippell, The Journal of cell biology 64, 98 (1975).
  * Huang et al. (2021) M. Huang, W. Hu, S. Yang, Q.-X. Liu, and H. Zhang, Proceedings of the National Academy of Sciences 118 (2021).
  * Yan et al. (2012) J. Yan, M. Bloom, S. C. Bae, E. Luijten, and S. Granick, Nature 491, 578 (2012).
  * Yan et al. (2015) J. Yan, S. C. Bae, and S. Granick, Soft Matter 11, 147 (2015).
  * Hwang et al. (2020) S. Hwang, T. D. Nguyen, S. Bhaskar, J. Yoon, M. Klaiber, K. J. Lee, S. C. Glotzer, and J. Lahann, Advanced Functional Materials 30, 1907865 (2020).
  * Zhang et al. (2020) B. Zhang, A. Sokolov, and A. Snezhko, Nature communications 11, 1 (2020).
  * Bricard et al. (2015) A. Bricard, J.-B. Caussin, D. Das, C. Savoie, V. Chikkadi, K. Shitara, O. Chepizhko, F. Peruani, D. Saintillan, and D. Bartolo, Nature communications 6, 1 (2015).
  * Zhang et al. (2021) B. Zhang, H. Karani, P. M. Vlahovska, and A. Snezhko, Soft Matter (2021).
  * Manna et al. (2021) R. K. Manna, O. E. Shklyaev, and A. C. Balazs, Proceedings of the National Academy of Sciences 118 (2021).
  * Barciś et al. (2019) A. Barciś, M. Barciś, and C. Bettstetter, in _2019 International Symposium on Multi-Robot and Multi-Agent Systems (MRS)_ (IEEE, 2019), pp. 98–104.
  * Barciś and Bettstetter (2020) A. Barciś and C. Bettstetter, IEEE Access 8, 218752 (2020).
  * Monaco et al. (2020) J. D. Monaco, G. M. Hwang, K. M. Schultz, and K. Zhang, Biological cybernetics 114, 269 (2020).
  * Monaco et al. (2019) J. D. Monaco, G. M. Hwang, K. M. Schultz, and K. Zhang, in _Micro-and Nanotechnology Sensors, Systems, and Applications XI_ (International Society for Optics and Photonics, 2019), vol. 10982, p. 109822D.
  * Uriu and Morelli (2017) K. Uriu and L. G. Morelli, Development, growth & differentiation 59, 351 (2017).
  * Uriu and Morelli (2014) K. Uriu and L. G. Morelli, Biophysical journal 107, 514 (2014).
  * Tsiairis and Aulehla (2016) C. D. Tsiairis and A. Aulehla, Cell 164, 656 (2016).
  * Hrabec et al. (2018) A. Hrabec, V. Křižáková, S. Pizzini, J. Sampaio, A. Thiaville, S. Rohart, and J. Vogel, Physical review letters 120, 227204 (2018).
  * Tanaka (2007) D. Tanaka, Physical review letters 99, 134103 (2007).
  * Iwasa and Tanaka (2017) M. Iwasa and D. Tanaka, Physics Letters A 381, 3054 (2017).
  * Iwasa and Tanaka (2010) M. Iwasa and D. Tanaka, Physical Review E 81, 066214 (2010).
  * Levis et al. (2019) D. Levis, I. Pagonabarraga, and B. Liebchen, Physical Review Research 1, 023026 (2019).
  * Hong (2018) H. Hong, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 103112 (2018).
  * Lee et al. (2021) H. K. Lee, K. Yeo, and H. Hong, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 033134 (2021).
  * Jiménez-Morales (2020) F. Jiménez-Morales, Physical Review E 101, 062202 (2020).
  * Lizarraga and de Aguiar (2020) J. U. Lizarraga and M. A. de Aguiar, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 053112 (2020).
  * McLennan-Smith et al. (2020) T. A. McLennan-Smith, D. O. Roberts, and H. S. Sidhu, Physical Review E 102, 032607 (2020).
  * O’Keeffe et al. (2018) K. P. O’Keeffe, J. H. Evers, and T. Kolokolnikov, Physical Review E 98, 022203 (2018).
  * O’Keeffe and Bettstetter (2019) K. O’Keeffe and C. Bettstetter, in _Micro-and Nanotechnology Sensors, Systems, and Applications XI_ (International Society for Optics and Photonics, 2019), vol. 10982, p. 109822E.
  * Ha et al. (2021) S.-Y. Ha, J. Jung, J. Kim, J. Park, and X. Zhang, Kinetic & Related Models (2021).
  * Ha et al. (2019) S.-Y. Ha, J. Jung, J. Kim, J. Park, and X. Zhang, Mathematical Models and Methods in Applied Sciences 29, 2225 (2019).
  * Bau et al. (2015) H. H. Bau, D. Raizen, and J. Yuan, in _Worm_ (Taylor & Francis, 2015), vol. 4, p. e1118606.
  * Yuan et al. (2015) J. Yuan, D. M. Raizen, and H. H. Bau, Journal of The Royal Society Interface 12, 20150227 (2015).
  * Nosrati et al. (2016) R. Nosrati, P. J. Graham, Q. Liu, and D. Sinton, Scientific reports 6, 1 (2016).
  * Ketzetzi et al. (2021) S. Ketzetzi, M. Rinaldin, P. Dröge, J. de Graaf, and D. J. Kraft, arXiv preprint arXiv:2103.07335 (2021).
  * 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).
  * Aihara (2009) I. Aihara, Physical Review E 80, 011918 (2009).
  * Aihara et al. (2014) I. Aihara, T. Mizumoto, T. Otsuka, H. Awano, K. Nagira, H. G. Okuno, and K. Aihara, Scientific reports 4, 1 (2014).
  * Le et al. (2019) T. C. Le, J. Zhai, W.-H. Chiu, P. A. Tran, and N. Tran, International journal of nanomedicine 14, 6749 (2019).
  * Agrawal and Agrawal (2019) G. Agrawal and R. Agrawal, ACS Applied Nano Materials 2, 1738 (2019).
  * Yi et al. (2016) Y. Yi, L. Sanchez, Y. Gao, and Y. Yu, Analyst 141, 3526 (2016).
  * Xiao et al. (2018) Z. Xiao, M. Wei, and W. Wang, ACS applied materials & interfaces 11, 6667 (2018).
  * Strogatz (2000) S. H. Strogatz, Physica D: Nonlinear Phenomena 143, 1 (2000).
  * Strogatz and Mirollo (1991) S. H. Strogatz and R. E. Mirollo, Journal of Statistical Physics 63, 613 (1991).
  * Mirollo and Strogatz (2007) R. Mirollo and S. H. Strogatz, Journal of Nonlinear Science 17, 309 (2007).
  * Crawford (1994) J. D. Crawford, Journal of Statistical Physics 74, 1047–1084 (1994).
  * O’Keeffe et al. (2017) K. P. O’Keeffe, H. Hong, and S. H. Strogatz, Nature communications 8, 1 (2017).
  * Yuan et al. (2014) J. Yuan, D. M. Raizen, and H. H. Bau, Proceedings of the National Academy of Sciences 111, 6865 (2014).
  * Guzmán-Lastra et al. (2016) F. Guzmán-Lastra, A. Kaiser, and H. Löwen, Nature communications 7, 1 (2016).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 037113 (2008).
  * Pikovsky and Rosenblum (2015) A. Pikovsky and M. Rosenblum, Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 097616 (2015).
  * Engelbrecht and Mirollo (2020) J. R. Engelbrecht and R. Mirollo, Physical Review Research 2, 023057 (2020).
  * Chen et al. (2017) B. Chen, J. R. Engelbrecht, and R. Mirollo, Journal of Physics A: Mathematical and Theoretical 50, 355101 (2017).
  * Igoshin et al. (2001) O. A. Igoshin, A. Mogilner, R. D. Welch, D. Kaiser, and G. Oster, Proceedings of the National Academy of Sciences 98, 14913 (2001).

##  Appendix A Basin of attractions for static sync and $\pi$-state

Are the static sync and $\pi$ states the only stable collective states when
$K>0$? Answering this question conclusively is difficult since it would mean
integrating the equations of motion for every initial condition in the phase
space $T^{N}=S^{N}\times S^{N}$. Figure 10 explores a subset of $T^{N}$:
initial position $x_{i}$ spaced uniformly on $(0,a\pi)$ and initial phases
$\theta_{i}$ spaced uniformly on $(0,b\pi)$ where $0\leq a,b\leq 2$ and then
perturbed slightly with noise of order $10^{-3}$ to both $x_{i},\theta_{i}$.
When $(a,b)=(2,2)$ these initial conditions correspond to swarmalators spread
out evenly over the unit circle in both space and phase. When $a,b<2$, the
swarmalators are spread out over subsets of the unit circle.

We use the Daido order parameters
$R_{1}e^{i\phi_{1}}:=N^{-1}\sum_{j}e^{i\theta_{j}}$ and
$R_{2}e^{i\phi_{2}}=N^{-1}\sum_{j}e^{2i\theta_{j}}$ to distinguish between the
static and $\pi$ states. The conditions for each state are

  1. 1.

Condition for static sync: $(R_{1},R_{2})=(1,1)$

  2. 2.

Condition for $\pi$-state: $(R_{1},R_{2})=(0,1)$

While $R_{2}=1$ for both the static sync and $\pi$ states, and thus can’t
distinguish between them, we include it to rule out the existence of some
other collective state for which $R_{2}\neq 0$.

Figure 10 shows that the static sync is stable when $(0\leq a\leq 1)\cap(0\leq
b\leq 1)$ and the $\pi$-state is stable in the remaining space. No other
collective states were observed.

Figure 10: Left panel: $R_{1}(a,b)$. Right panel: $R_{2}(a,b)$. Data were
collected by integrating Eqs.(1), (2) using an RK4 method with
$(dt,T)=(0.1,500)$ for $N=100$ swarmalators.

##  Appendix B Connection of ring model to 2D swarmalator model

Figure 11: Left panel: $S_{+}$ for the unit vector model which has the same
qualitative shape as that of the linear parabolic model presented in Figure 1.
Results for $T=500,1000$ are shown to establish convergence. Right panel: zoom
in of $S_{+}(K)$ showing the non-monotonic dip. Results were generated using
an RK4 solver with $dt=0.1$.

Here we show how the ring model is contained within the 2D swarmalator model
which is given by

$\displaystyle\dot{\mathbf{x}}_{i}=\mathbf{v}_{i}+\frac{1}{N}\sum_{j=1}^{N}\Big{[}\mathbf{I}_{\mathrm{att}}(\mathbf{x}_{j}-\mathbf{x}_{i})F(\theta_{j}-\theta_{i})-\mathbf{I}_{\mathrm{rep}}(\mathbf{x}_{j}-\mathbf{x}_{i})\Big{]},$
(86)
$\displaystyle\dot{\theta_{i}}=\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}H_{\mathrm{\mathrm{att}}}(\theta_{j}-\theta_{i})G_{\sigma}(\mathbf{x}_{j}-\mathbf{x}_{i})$
(87)

In O’Keeffe et al. (2017), the choices $I_{att}=x/|x|$, $I_{rep}=x/|x|^{2}$,
$F(\theta)=1+J\cos(\theta)$, $G(x)=1/|x|$, $H_{att}(\theta)=\sin(\theta)$ were
made. However, choosing linear spatial attraction $I_{att}(x)=x$, inverse
square spatial repulsion $I_{rep}(x)=x/|x|^{2}$ and truncated parabolic space-
phase coupling $G(x)=(1-|x|^{2}/\sigma^{2})H_{heaviside}({\sigma-|x|})$

$\displaystyle\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{j\neq
i}^{N}\Bigg{[}\mathbf{x}_{j}-\mathbf{x}_{i}\Big{(}1+J\cos(\theta_{j}-\theta_{i})\Big{)}-\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\mathbf{x}_{j}-\mathbf{x}_{i}|^{2}}\Bigg{]}$
(88) $\displaystyle\dot{\theta_{i}}=\frac{K}{N}\sum_{j\neq
i}^{N}\sin(\theta_{j}-\theta_{i})\Big{(}1-\frac{|\mathbf{x}_{j}-\mathbf{x}_{i}|^{2}}{\sigma^{2}}\Big{)}H_{heaviside}(\sigma-|\mathbf{x}_{j}-\mathbf{x}_{i}|)$
(89)

gives the same qualitative behavior but is nicer to work with analytically.

First we show it demonstrates the same behavior. Figure 2 shows its collective
states are the same as those of the original model, and Figure 1 shows its
order parameters $S_{\pm}(K)$ have the same shape: monotonic increase from
$(K_{1},K_{2})$ as $K\rightarrow 0^{-1}$, monontonic decrease from $(K_{2},0)$
as $K\rightarrow 0^{-1}$, and then a discontinuous jump to $S_{+}=1$ at $K=0$
($K<0$ is a singular perturbation). Note the plot of $S_{\pm}$(K) in Figure 6
of O’Keeffe et al. (2017), the monotonic decrease on $(K_{2},0)$ of $S_{+}$ is
slight and hard to see, so we plot in for finer $K$ and show a zoom in for
small $K$ in Figure LABEL:order-parameters-2D-unit-vector.

Now we show the ‘linear parabolic‘ model, so called because $I_{att}=x$ and
$G(x)$ is a parabolic, is cleaner analytically. In polar coordinates it takes
form

$\displaystyle\dot{r_{i}}$
$\displaystyle=H_{r}(r_{i},\phi_{i})-Jr_{i}R_{0}\cos\Big{(}\Psi_{0}-\theta_{i}\Big{)}+\frac{J}{2}\Bigg{[}\tilde{S}_{+}\cos\Big{(}\Phi_{+}-(\phi_{i}+\theta_{i})\Big{)}+\tilde{S}_{-}\cos\Big{(}\Phi_{-}-(\phi_{i}-\theta_{i})\Big{)}\Bigg{]}$
$\displaystyle\dot{\phi_{i}}$
$\displaystyle=H_{\phi}(r_{i},\phi_{i})+\frac{J}{2r_{i}}\Bigg{[}\tilde{S}_{+}\sin\Big{(}\Psi_{+}-(\phi_{i}+\theta_{i})\Big{)}+\tilde{S}_{-}\sin\Big{(}\Psi_{-}-(\phi_{i}-\theta_{i})\Big{)}\Bigg{]}$
$\displaystyle\dot{\theta_{i}}$
$\displaystyle=K\Big{(}1-\frac{r_{i}^{2}}{\sigma^{2}}\Big{)}R_{0}\sin(\Phi_{0}-\theta_{i})-\frac{K}{\sigma^{2}}R_{1}\sin(\Phi_{1}-\theta_{i})+\frac{Kr_{i}}{\sigma^{2}}\Bigg{[}\tilde{S}_{+}\sin\Big{(}\Psi_{+}-(\phi_{i}+\theta_{i})\Big{)}-\tilde{S}_{-}\sin\Big{(}\Psi_{-}-(\phi_{i}-\theta_{i})\Big{)}\Bigg{]}$

where

$\displaystyle H_{r}(r_{i},\phi_{i})$
$\displaystyle=\frac{1}{N}\sum_{j}\Big{(}r_{j}\cos(\phi_{j}-\phi_{i})-r_{i}\Big{)}(1-d_{ij}^{-2})$
(90) $\displaystyle H_{\phi}(r_{i},\phi_{i})$
$\displaystyle=\frac{1}{N}\sum_{j}\frac{r_{j}}{r_{i}}\sin(\phi_{j}-\phi_{i})(1-d_{ij}^{-2}),$
(91) $\displaystyle Z_{0}=R_{0}e^{i\Psi_{0}}$
$\displaystyle=\frac{1}{N}\sum_{j}e^{i\theta_{j}},$ (92)
$\displaystyle\hat{Z}_{0}=\hat{R}_{0}e^{i\hat{\Psi}_{0}}$
$\displaystyle=\frac{1}{N}\sum_{j\in N_{i}}e^{i\theta_{j}},$ (93)
$\displaystyle Z_{2}=R_{2}e^{i\Psi_{2}}$
$\displaystyle=\frac{1}{N}\sum_{j}r_{j}^{2}e^{i\theta_{j}},$ (94)
$\displaystyle\hat{Z}_{2}=\hat{R}_{2}e^{i\hat{\Psi}_{2}}$
$\displaystyle=\frac{1}{N}\sum_{j\in N_{i}}r_{j}^{2}e^{i\theta_{j}},$ (95)
$\displaystyle\tilde{W}_{\pm}=\tilde{S}_{\pm}e^{i\Psi_{\pm}}$
$\displaystyle=\frac{1}{N}\sum_{j}r_{j}e^{i(\phi_{j}\pm\theta_{j})}$ (96)
$\displaystyle\hat{W}_{\pm}=\hat{S}_{\pm}e^{i\hat{\Psi}_{\pm}}$
$\displaystyle=\frac{1}{N}\sum_{j\in N_{i}}r_{j}e^{i(\phi_{j}\pm\theta_{j})}$
(97)

where the $\hat{Z_{0}},\dots$ order parameters are summed over all the
neighbours $N_{i}$ of the $i$-th swarmalator: those within a distance
$\sigma$. Notice that rainbow order parameters $\tilde{W}$ here are weighted
by the radial distance $r_{j}$, which is not the case for the ring model
Eqs.(10), (11) (that’s why we put a tilde over the W). Assuming
$\sigma>max(d_{ij})$, we can set
$\hat{Z_{0}}=Z_{0},\hat{Z}_{1}=Z_{1},\hat{W_{\pm}}=W_{\pm}$. Then
$S_{\pm}\sin(\Phi_{\pm}-(\phi\pm\theta))$ etc of the ring model starting to
emerge. If we assume there is no global synchrony $Z_{0}=Z_{2}=0$, which
happens generically in the frustrated parameter regime $K<0,J>0$, and
transform to $\xi_{i}=\phi_{i}+\theta_{i}$ and $\eta_{i}=\phi_{i}-\theta_{i}$
coordinates the ring model is revealed (the terms in the square parentheses in
the latter two equations.)

$\displaystyle\dot{r_{i}}$
$\displaystyle=\tilde{\nu}(r_{i})+\frac{J}{2}\Bigg{[}\tilde{S}_{+}\cos\Big{(}\Phi_{+}-\xi_{i}\Big{)}+\tilde{S}_{-}\cos\Big{(}\Phi_{-}-\eta_{i}\Big{)}\Bigg{]}$
(98) $\displaystyle\dot{\xi_{i}}$
$\displaystyle=\tilde{\omega}(r_{i},\phi_{i})+\Bigg{[}J_{+}(r_{i})\tilde{S}_{+}\sin\Big{(}\Psi_{+}-\xi_{i}\Big{)}+J_{-}(r_{i})\tilde{S}_{-}\sin\Big{(}\Psi_{-}-\eta_{i}\Big{)}\Bigg{]}$
(99) $\displaystyle\dot{\eta_{i}}$
$\displaystyle=\tilde{\omega}(r_{i},\phi_{i})+\Bigg{[}J_{-}(r_{i})\tilde{S}_{+}\sin\Big{(}\Psi_{+}-\xi_{i}\Big{)}-J_{+}(r_{i})\tilde{S}_{-}\sin\Big{(}\Psi_{-}-\eta_{i}\Big{)}\Bigg{]}$
(100)

where

$\displaystyle\tilde{\nu}(r_{i},\phi_{i})$
$\displaystyle=H_{r}(r_{i},\phi_{i})$ (101)
$\displaystyle\tilde{\omega}(r_{i},\phi_{i})$
$\displaystyle=H_{\phi}(r_{i},\phi_{i})$ (102) $\displaystyle J_{\pm}(r_{i})$
$\displaystyle=\frac{J}{2r_{i}}\pm\frac{Kr_{i}}{\sigma^{2}}$ (103)

##  Appendix C Stability of static async state

The density obeys the continuity equation

$\dot{\rho}+\nabla.(v\rho)=0$ (104)

And the velocity $v=(v_{x},v_{\theta})$ is given by the $N\rightarrow\infty$
limit of Eqs (1), (2)

$\displaystyle v_{x}$
$\displaystyle=J\int\sin(x^{\prime}-x)\cos(\theta^{\prime}-\theta)\rho(x^{\prime},\theta^{\prime},t)dx^{\prime}d\theta^{\prime}$
(105) $\displaystyle v_{\theta}$
$\displaystyle=K\int\sin(\theta^{\prime}-\theta)\cos(x^{\prime}-x)\rho(x^{\prime},\theta^{\prime},t)dx^{\prime}d\theta^{\prime}$
(106)

Consider a perturbation around the static async state $\rho_{0}=1/(4\pi^{2})$

$\rho=\rho_{0}+\epsilon\eta=(4\pi^{2})^{-1}+\epsilon\eta(x,\theta,t)$ (107)

Normalization requires $\int\rho(x,\theta)=1$ which implies

$\int\eta(x,\theta,t)dxd\theta=0$ (108)

The density ansatz (107) decomposes the velocity

$v=v_{0}+\epsilon v_{1}=\epsilon v_{1}$ (109)

where $v_{0}$ is velocity in the static async state $v_{0}=0$. The perturbed
velocity $v_{1}$ is given by Eqs (112) with $\rho$ replaced by $\eta$.
Plugging Eqs (107), (109) into (104) yields

$\dot{\eta}+\rho_{0}(\nabla.v_{1})=0$ (110)

To tackle this, first write the $v_{1}$ in terms of the order parameters
$W_{\pm}$

$\displaystyle v_{x}^{1}$
$\displaystyle=\frac{J}{2}\text{Im}\Big{(}W_{+}^{1}e^{-i(x+\theta)}+W_{-}^{1}e^{-i(x-\theta)}\Big{)}$
(111) $\displaystyle v_{\theta}^{1}$
$\displaystyle=\frac{K}{2}\text{Im}\Big{(}W_{+}^{1}e^{-i(x+\theta)}-W_{-}^{1}e^{-i(x-\theta)}\Big{)}$
(112)

where the perturbed order parameters are

$W_{\pm}^{1}=\int
e^{i(x^{\prime}\pm\theta^{\prime})}\eta(x^{\prime},\theta^{\prime})dx^{\prime}d\theta^{\prime}$
(113)

The divergence term is
$\nabla.v_{1}.=\partial_{x}v_{x}+\partial_{\theta}v_{\theta}$ is easy to
compute since only exponential terms

$\nabla.v_{1}=-\frac{J+K}{2}\text{Re}\Big{(}W_{+}^{1}e^{-i(x+\theta)}+W_{-}^{1}e^{-i(x-\theta)}\Big{)}$
(114)

Plugging this into the evolution equation for $\eta$ Eq. (110) yields

$\dot{\eta}(x,\theta,t)=\frac{J+K}{8\pi^{2}}\text{Re}\Big{(}W_{+}^{1}e^{-i(x+\theta)}+W_{-}^{1}e^{-i(x-\theta)}\Big{)}$
(115)

Using $\text{Re}\;z=\frac{1}{2}(z+\bar{z})$ we get

$\displaystyle\dot{\eta}$
$\displaystyle=\frac{J+K}{16\pi^{2}}\Big{(}W_{+}^{1}e^{-i(x+\theta)}+\bar{W}_{+}^{1}e^{i(x+\theta)}+W_{-}^{1}e^{-i(x-\theta)}$
$\displaystyle+\bar{W}_{-}^{1}e^{i(x-\theta)}\Big{)}$ (116)

Only a few Fourier modes are distinguished so we expand $\eta(x,\theta,t)$ in
terms of complex exponentials

$\displaystyle\eta$
$\displaystyle=(2\pi)^{-1}\Big{(}\alpha_{0,0}(t)+\alpha_{0,1}(t)e^{i\theta}+\alpha_{1,0}(t)e^{ix}$
(117) $\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\bar{\alpha}_{n,m}(t)e^{i(nx+m\theta)}+\bar{\beta}_{n,m}(t)e^{i(nx-m\theta)}$
$\displaystyle+c.c.\Big{)}$ (118)

where $c.c$ denotes the complex conjugate. Unconventionally, we associate the
complex conjugate $\bar{\alpha},\bar{\beta}$ with the first harmonic. This is
so that the order parameters are expressed in terms of
$W_{+}^{1}=\alpha_{1,1}(t)$ and $W_{-}^{1}=\beta_{1,1}(t)$ (i.e. without the
bar overhead). Also, notice the normalization condition (108) implies
$\alpha_{0,0}(t)=0$.

Integrating Eq. (116) with respect to $\int(.)e^{i(nx\pm n\theta)}$ extracts
the evolution equations for each mode. The ones of interest are

$\displaystyle\dot{W}_{\pm}^{1}=\frac{J+K}{4\pi}W_{\pm}^{1}$ (119)

and $\dot{\alpha}_{n,m}=\dot{\beta}_{n,m}=0$ for all $(n,m)\neq(1,1)$. Setting
$W_{\pm}^{1}=w_{\pm}^{1}e^{\lambda_{\pm}t}$ yields

$\displaystyle\lambda_{\pm}=\frac{J+K}{4\pi}$ (120)

which implies the static async state is stable for

$\displaystyle K<K_{c}=-J$ (121)

consistent with simulations.