Oscillators that sync and swarm

TL;DR

This paper introduces 'swarmalators', systems where oscillators synchronize their phases while moving through space, revealing five possible collective states with applications in biological and physical systems.

Contribution

It proposes a coupled phase and spatial dynamics model called swarmalators, and predicts five distinct long-term collective states.

Findings

Five collective states predicted by the model

Applicable to biological systems like sperm and frogs

Potential observations in colloidal suspensions

Abstract

Synchronization occurs in many natural and technological systems, from cardiac pacemaker cells to coupled lasers. In the synchronized state, the individual cells or lasers coordinate the timing of their oscillations, but they do not move through space. A complementary form of self-organization occurs among swarming insects, flocking birds, or schooling fish; now the individuals move through space, but without conspicuously altering their internal states. Here we explore systems in which both synchronization and swarming occur together. Specifically, we consider oscillators whose phase dynamics and spatial dynamics are coupled. We call them swarmalators, to highlight their dual character. A case study of a generalized Kuramoto model predicts five collective states as possible long-term modes of organization. These states may be observable in groups of sperm, Japanese tree frogs, colloidal suspensions of magnetic particles, and other biological and physical systems in which self-assembly and synchronization interact.

Figures

Equations

$$\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{]},$$

$$\dot{{\theta }_i}={\omega }_i+\frac{K}{N}\sum _{j=1}^N{H}_{\mathrm{att}}({\theta }_j-{\theta }_i)G({\mathbf{x}}_j-{\mathbf{x}}_i)$$

$$\displaystyle\dot{\mathbf{x}}_{i}=\mathbf{v}_{i}+\frac{1}{N}\sum_{j\neq i}^{N}\Bigg{[}\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\mathbf{x}_{j}-\mathbf{x}_{i}|}\Big{(}A+J\cos(\theta_{j}-\theta_{i})\Big{)}$$

$$\displaystyle\hskip 142.26378pt-B\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\mathbf{x}_{j}-\mathbf{x}_{i}|^{2}}\Bigg{]}$$

$$\dot{{\theta }_i}={\omega }_i+\frac{K}{N}\sum _{j\ne i}^N\frac{\sin ({\theta }_j-{\theta }_i)}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}.$$

$$\displaystyle\int_{0}^{R}\Big{[}(s-r)\mathcal{K}\left(\frac{4rs}{(r+s)^{2}}\right)+(r+s)\mathcal{E}\left(\frac{4rs}{(r+s)^{2}}\right)$$

$$\displaystyle+\frac{\pi^{2}}{2J}(r-s)\Big{]}\frac{2Js}{r}g_{1}(s)\;ds=0$$

$$\displaystyle g_{1}(r)=\frac{2(1+J)}{\pi}\int_{0}^{R}\mathcal{K}\Bigg{(}\frac{4sr}{(r+s)^{2}}\Bigg{)}\frac{g_{1}(s)}{s+r}s\;ds,$$

$$R_{\mathrm{sync}}=(1+J{)}^{-1/2}.$$

$$R_{\mathrm{async}}=1$$

$$g_2(r)=1-\frac{{\Gamma }_J}{r},R_1\le r\le R_2$$

$$R_1={\Delta }_J\frac{-\sqrt{3}J-3\sqrt{12-5J}\sqrt{J+4}+12\sqrt{3}}{12J},$$

$$R_2=\frac{{\Delta }_J}{2\sqrt{3}}$$

$$W_{\pm }=S_{\pm }{e}^{i{\Psi }_{\pm }}=\frac{1}{N}\sum _{j=1}^N{e}^{i({\varphi }_j\pm {\theta }_j)},$$

$$\rho (\mathbf{x},\theta ,t)$$

$${\dot{b}}_1(\mathbf{x},t)$$

$$\frac{(J+K)}{2}{\rho }_0(\mathbf{x})\int \frac{1}{|\tilde{\mathbf{x}}-\mathbf{x}|}{b}_1(\tilde{\mathbf{x}},t)d\tilde{\mathbf{x}}.$$

$${\lambda }_m{c}_m(r)$$

$${\lambda }_0{c}_0(r)$$

$$H_0(r,s)$$

$$+\frac{J(r+s)g^{\prime }(r)}{4{\pi }^{2}r}\mathcal{E}\left(\frac{4rs}{(r+s{)}^2}\right)$$

$$K_c\approx -1.2J.$$

$$\displaystyle\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{j\neq i}^{N}\Bigg{[}\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\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{]}$$

$$+{\xi }_i^{\mathbf{x}}(t),$$

$$\dot{{\theta }_i}={\omega }_i+\frac{K}{N}\sum _{j\ne i}^N\frac{\sin ({\theta }_j-{\theta }_i)}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}+{\eta }_i(t),$$

$$\displaystyle\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{j\neq i}^{N}\Bigg{[}\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\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}|^{3}}\Bigg{]},$$

$$\dot{{\theta }_i}=\frac{K}{N}\sum _{j\ne i}^N\frac{\sin ({\theta }_j-{\theta }_i)}{|{\mathbf{x}}_j-{\mathbf{x}}_i|},$$

$$R_{\mathrm{sync}}$$

$$R_{\mathrm{async}}$$

$$\displaystyle\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{j\neq i}^{N}\Bigg{[}\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\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{]}$$

$$+{\xi }_i^{\mathbf{x}}(t)+v_0\hat{n},$$

$$\dot{{\theta }_i}={\omega }_i+\frac{K}{N}\sum _{j\ne i}^N\frac{\sin ({\theta }_j-{\theta }_i)}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}+{\eta }_i(t),$$

$$\dot{{\beta }_i}=-{\beta }_i+\frac{1}{|{\Lambda }_i|}\sum _{j\in {\Lambda }_i}{\beta }_j+{\zeta }_i(t),$$

$$\rho (r,\varphi ,\theta ,t)=\frac{1}{4{\pi }^2}{g}_1(r),0\le r\le R.$$

$$\underline{v}$$

$$\nabla .\underline{v}$$

$$\displaystyle\mathbf{v}_{\mathbf{x}}(\mathbf{x},\theta,t)=\int\Bigg{(}\Big{(}\mathbf{\tilde{x}}-\mathbf{x}\Big{)}\Big{(}1+J\cos(\tilde{\theta}-\theta)\Big{)}-$$

$$\displaystyle\hskip 113.81102pt\frac{\mathbf{\tilde{x}}-\mathbf{x}}{|\mathbf{\tilde{x}}-\mathbf{x}|^{2}}\Bigg{)}\rho(\mathbf{\tilde{x}},\;\tilde{\theta},t)\;d\mathbf{\tilde{x}}\,d\tilde{\theta},$$

$$v_{\theta }(\mathbf{x},\theta ,t)=\int \frac{\sin (\tilde{\theta }-\theta )}{|\tilde{\mathbf{x}}-\mathbf{x}|}\rho (\tilde{\mathbf{x}},\tilde{\theta },t)d\tilde{\mathbf{x}}d\tilde{\theta }.$$

$$\displaystyle\nabla_{\mathbf{x}}.\mathbf{v}_{\mathbf{x}}=\int-2\Big{(}1+J\cos(\tilde{\theta}-\theta)\Big{)}\rho(\mathbf{\tilde{x}},\;\tilde{\theta},t)d\mathbf{\tilde{x}}\,d\tilde{\theta}$$

Full text

This year marks the fiftieth anniversary of a breakthrough in the study of synchronization. In 1967, Winfree proposed a coupled oscillator model for the circadian rhythms that underlie daily cycles of activity in virtually all plants and animals winfree1967biological . He discovered that above a critical coupling strength, synchronization breaks out spontaneously, in a manner reminiscent of a phase transition. Then Kuramoto simplified Winfree’s model and solved it exactly kuramoto1975self , leading to an explosion of interest in the dynamics of coupled oscillators strogatz_book ; pikovsky2003synchronization ; acebron2005kuramoto . Kuramoto’s model in turn has been generalized to other large systems of biological oscillators, such as chorusing frogs aihara2008mathematical , firing neurons montbrio2015macroscopic ; pazo2014low ; o2016dynamics ; luke2014macroscopic ; laing2014derivation , and even human concert audiences clapping in unison crowd . The analyses often borrow techniques from statistical physics, such as mean-field approximations, renormalization group analyses daido1988lower ; ostborn2009renormalization , and finite-size scaling hong2007entrainment ; hong2015finite . There has also been traffic in the other direction, from biology back to physics. For example, insights from biological synchronization have shed light on neutrino oscillations neutrinos , phase locking in Josephson junction arrays wiesenfeld1996synchronization , the dynamics of power grids motter2013spontaneous ; dorfler2012synchronization , and the unexpected wobbling of London’s Millennium Bridge on opening day strogatz2005theoretical .

A similarly fruitful interplay between physics and biology has occurred in the study of the coordinated movement of groups of animals. Fish schools, bird flocks, and insect swarms couzin2007collective ; buhl2006disorder ; sumpter2010collective ; herbert2016understanding ; Ballerini29012008 have been illuminated by maximum entropy methods Bialek27032012 , agent-based simulations reynolds1987flocks , and analytically tractable models based on self-propelled particles vicsek1995novel , and continuum limits bernoff2013nonlocal ; topaz2004swarming ; topaz2006nonlocal ; kolokolnikov2011stability .

Studies of swarming and synchronization have much in common. Both involve large, self-organizing groups of individuals interacting according to simple rules. Both lie at the intersection of nonlinear dynamics and statistical physics. Nevertheless the two fields have, by and large, remained disconnected. Studies of swarms focus on how animals move, while neglecting the dynamics of their internal states. Studies of synchronization do the opposite: they focus on oscillators’ internal dynamics, not on their motion. In the past decade, however, a few studies of “mobile oscillators,” motivated by applications in robotics and developmental biology, have brought the two fields into contact mobile1 ; belykh2004blinking ; stilwell2006sufficient ; frasca2008synchronization ; fujiwara2011synchronization . Even so, the assumption has been that the oscillators’ locations affect their phase dynamics, but not conversely. Their motion has been modeled as a random walk or as externally determined, without feedback from the oscillators’ phases.

We suspect that somewhere in nature and technology there must be mobile oscillators whose phases affect how they move. For instance, many species of frogs, crickets, and katydids call periodically, and synchronize in vast choruses walker1969acoustic ; greenfield1994synchronous ; aihara2008mathematical ; aihara2014spatio . The natural question is whether they tend to hop toward or away from others depending on the relative phases of their calling rhythms, and if so, what spatiotemporal patterns are produced.

A clue comes from the physics of magnetic colloids yan2012linking ; snezhko2011magnetic ; martin2013driving and microfluidic mixtures of active spinners nguyen2014emergent ; van2016spatiotemporal , both of which show rich collective behavior. In these systems, the particles or spinners attract or repel one another, depending on their orientations. Given that orientation is formally analogous to the phase of an oscillation (both being circular variables), a similarly rich phenomenology is expected for mobile oscillators whose phases affect their motion. We call these hypothetical systems ‘swarmalators’ because they generalize swarms and oscillators.

One possible instance of a swarmalator system is a population of myxobacteria, modeled in 2001 by Igoshin and colleagues igoshin2001pattern . The movements of these bacteria in space are thought to be influenced by an internal, biochemical degree of freedom, which appears to vary cyclically. Igoshin et al. igoshin2001pattern modeled it as a phase oscillator. Experimental evidence suggests that the evolution of this phase is influenced by the spatial density of neighboring cells; thus there appears to be a bidirectional coupling between spatial and phase dynamics, as required of swarmalators.

Tanaka and colleagues also made an early contribution to the modeling of swarmalators tanaka2007general ; iwasa2010hierarchical . They analyzed a broad class of models in the hope of finding phenomena which were not system-specific. They considered chemotactic oscillators, whose movements in space are mediated by the diffusion of a background chemical. The oscillators’ consumption of this chemical depends on their internal states, thereby completing the bidirectional space-phase coupling. Tanaka et al. tanaka2007general ; iwasa2010hierarchical began with a general model with these ingredients, from which they derived a simpler model by means of center manifold and phase-reduction methods.

Here we take a bottom-up approach. We propose a simple model of a swarmalator system which lets us study some of its collective states analytically. We hope our work will draw attention to this class of problems, and stimulate the discovery and characterization of natural and technological systems of swarmalators.

Results

The model. We consider swarmalators free to move in the plane. The governing equations are

[TABLE]

for $i=1,\ldots,N$, where $N$ is the population size, $\mathbf{x}_{i}=(x_{i},y_{i})$ is the position of the $i$-th swarmalator, and $\theta_{i},\omega_{i}$ and $\mathbf{v}_{i}$ are its phase, natural frequency, and self-propulsion velocity. The functions $\mathbf{I}_{\mathrm{att}}$ and $\mathbf{I}_{\mathrm{rep}}$ represent the spatial attraction and repulsion between swarmalators, while the phase interaction is captured by $H_{\mathrm{att}}$. The function $F$ in equation (1) measures the influence of phase similarity on spatial attraction, while $G$ in equation (2) measures the influence of spatial proximity on the phase attraction.

Consider the following instance of this model:

[TABLE]

For simplicity, we chose power laws for $\mathbf{I}_{\mathrm{att}},\mathbf{I}_{\mathrm{rep}}$ and $G$ along with analytically convenient exponents. The sine function in $H_{\mathrm{att}}$ was similarly motivated, in the spirit of the Kuramoto model kuramoto1975self . We first consider identical swarmalators so that $\omega_{i}=\omega$ and $\mathbf{v}_{i}=\mathbf{v}$. Further, we assume propulsion with constant magnitude and direction $\mathbf{v}=v_{0}\hat{n}$ where $\hat{n}$ is a constant vector (we relax these simplifications later). Then by a choice of reference frame we can set $\omega=v_{0}=0$ without loss of generality. Finally, by rescaling time and space we set $A=B=1$. This leaves us with a system with two parameters $(J,K)$.

The parameter $K$ is the phase coupling strength. For $K>0$, the phase coupling between swarmalators tends to minimize their phase difference, while for $K<0$, this phase difference is maximized. The parameter $J$ measures the extent to which phase similarity enhances spatial attraction. For $J>0$, “like attracts like”: swarmalators prefer to be near other swarmalators with the same phase. When $J<0$, we have the opposite scenario: swarmalators are preferentially attracted in space to those with opposite phase. And when $J=0$, swaramalators are phase-agnostic, their spatial attraction being independent of their phase. To keep $\mathbf{I}_{\mathrm{att}}>0$, we constrain $J$ to satisfy $-1\leq J\leq 1$.

Before stating our results, we pause to discuss our model’s features. As mentioned above, the model’s purpose is to study the interplay between synchronization and swarming. But what precisely do we mean by swarming? While, to our knowledge, there is no unanimous classification, elements of a swarming system typically (i) attract and repel each other, leading to aggregation, and (ii) align their orientations so as to move in the same direction. Succinctly then, a swarming system models aggregation and/or alignment.

Our model accounts for aggregation, but not for alignment: the spatial dynamics (1) model phase-dependent aggregation, while the phase dynamics (2) model position-dependent synchronization. There are no alignment terms. Indeed, the particles of our system do not have an orientation so there is nothing to align! We chose to neglect an orientation state variable, and thus alignment, for two reasons. The first was simply because we believe there are swarmalator systems in which orientation does not play a role, such as the Japanese tree frogs aihara2008mathematical ; aihara2014spatio or chemotactic oscillators tanaka2007general ; iwasa2010hierarchical . The second was that modeling orientable swarmalators adds an additional layer of complexity; it gives each swarmalator an orientation $\beta$, increasing the number of state variables per swarmalator from three (a two-dimensional position $(x,y)$ and an internal phase $\theta$) to four.

In the interest of minimalism we wished to avoid this complication for now. Hence as it stands our model applies only to swarmalators without an orientation. However we later show that our results are robust to the inclusion of simple alignment dynamics, indicating their potential to hold for systems of orientable swarmalators as well.

Numerics. We performed numerical experiments to probe the behavior of our system. Unless otherwise stated, the simulations were run using python’s ODE solver ‘odeint’. We initially positioned the swarmalators in a box of length $2$ and drew their phases from $[-\pi,\pi]$, both uniformly at random. We found the system settles into five states (Supplementary Movies 1-5). In three of these states, the swarmalators are ultimately static in space and phase. In the remaining two, the swarmalators move. However in all states, the density of swarmalators $\rho(\mathbf{x},\theta,t)$ is time-independent, where $\rho(\mathbf{x},\theta,t)\,d\mathbf{x}\,d\theta$ gives the fraction of swarmalators with positions between $\mathbf{x}$ and $\mathbf{x}+d\mathbf{x}$, and phases between $\theta$ and $\theta+d\theta$ at time $t$. In Fig. 1 we show where these states occur in the $(J,K)$ parameter plane. We next discuss these five states.

1. Static synchrony. The first state is shown in Fig. 2(a). The swarmalators form a circularly symmetric, crystal-like distribution in space, and are fully synchronized in phase, as indicated by all of them having the same color in Fig. 2(a). Since the swarmalators are ultimately stationary in $\mathbf{x}$, and they all end up at the same phase $\theta$, we call this the ‘static sync’ state. It occurs for $K>0$ and for all $J$, as seen in Fig. 1.

In the continuum limit, this state is described by $\rho(r,\phi,\theta,t)=\frac{1}{2\pi}g_{1}(r)\delta(\theta-\theta_{0})$, where $\phi$ is the spatial angle $\phi=\tan^{-1}(y/x)$, and the final phase $\theta_{0}$ is determined from the initial conditions. In Supplementary Note 1 we uset a technique used by Kololnikov et al fetecau2011swarm when studying swarms to derive the following pair of integral equations for $g_{1}$:

[TABLE]

where $\mathcal{K},\mathcal{E}$ are the complete elliptic integral of the first and second kinds, and $R$ is the radius of the disk in the $(x,y)$ plane which must be determined. We were unable to solve these equations for $g_{1}(r)$ and $R$, so instead solve them numerically, and show the results in Supplementary Note 1. Analytic progress can however be made if a linear attraction kernel $\mathbf{I}_{\mathrm{att}}(\mathbf{x})=\mathbf{x}$, is used instead of the unit vector kernel we are currently considering. Then, as shown in Kolokolnikov et al. fetecau2011swarm , the radial density becomes $g_{1}(r)=1$, i.e swarmalators are uniformly distributed. In this special case we can also calculate $R$ analytically,

[TABLE]

We show a full derivation in the Methods section. In dimensionful units, this reads $R=\sqrt{B/(A+J)}$. Thus the radius is determined by the ratio of the strengths of the attractive to the repulsive forces $\mathbf{I}_{\mathrm{att}},\mathbf{I}_{\mathrm{rep}}$ (in the static sync state, the effective attraction force is $A+J\cos(\theta_{j}-\theta_{i})=A+J$, since all swarmalators have the same phase). Figure 3(a) shows the prediction (7) agrees with simulation results.

2. Static asynchrony. Swarmalators can also form a static async state, illustrated in Fig. 2(b). At any given location $\mathbf{x}$, all phases $\theta$ can occur, and hence all colors are present everywhere in Fig. 2(b). This is seen more clearly in a scatter plot of the swarmalators in the $(\phi,\theta)$ plane, depicted in Fig. 4(a). Notice that the swarmalators are distributed uniformly, meaning that every phase occurs everywhere. This completely asynchronous state occurs in the quadrant $J<0$, $K<0$, and also for $J>0$ as long as $J$ lies in the wedge $J<|K_{c}|$ shown in the phase diagram in Fig. 1. As for the static sync state, we were able to calculate the radius of the circular distribution when a linear attraction kernel $\mathbf{I}_{\mathrm{att}}(\mathbf{x})$ was used. In the Methods section we show this radius is given by

[TABLE]

which agrees with simulation as shown in Fig. 3(a).

3. Static phase wave. The final stationary state occurs for the special case $K=0$ and $J>0$. This means the swarmalators’ phases are frozen at their initial values. How, then, does the population evolve? Since $J>0$, ‘like attracts like’: swarmalators want to settle near others with similar phase. The result is an annular structure where the spatial angle $\phi$ of each swarmalator is perfectly correlated with its phase $\theta$, as seen in Fig. 2(c) and 4(b). Since the phases run through a full cycle as the swarmalators arrange themselves around the ring, we call this state the ‘static phase wave’.

In density space, this static phase wave state is described by $\rho(r,\phi,\theta)=g_{2}(r)\delta(\phi\pm\theta+C_{1})$ where the $\pm$ and the constant $C_{1}$, are determined by the initial conditions. In the Methods section we again consider the linear attraction kernel, and find that $g_{2}(r)$ can be obtained analytically,

[TABLE]

with \Gamma_{J}=2J(R_{2}^{3}-R_{1}^{3})\Big{(}3J(R_{2}^{2}-R_{1}^{2})+12\Big{)}^{-1}. This in turn lets us calculate the inner and outer radii $R_{1},R_{2}$ of the annulus:

[TABLE]

with $\Delta_{J}=\sqrt{\frac{3J-\sqrt{36-15J}\sqrt{J+4}-12}{J-2}}$. Figure 3(b) shows agreement between these predictions and simulation.

4. Splintered phase wave. Moving from $K=0$ into the $K<0$ half-plane, we encounter the first non-stationary state, shown in Fig. 5(a) and Fig. 4(c). As can be seen, the static phase wave splinters into disconnected clusters of distinct phases. Accordingly we call this state the ‘splintered phase wave’. It is unclear what determines the number of clusters. Fewer are found when smaller length scales for the interaction functions $\mathbf{I}_{\mathrm{att}},\mathbf{I}_{\mathrm{rep}},G$ are used. However the parameters $J,K$ also play a role, although how precisely has not yet been determined. Within each cluster, the swarmalators “quiver,” executing small amplitude oscillations in both position and phase about their mean values.

5. Active phase wave. As $K$ is further decreased, these oscillations increase in amplitude until the swarmalators start to execute regular cycles in both spatial angle and phase. This motion is best illustrated in Fig. 4(d), in which shear flow about the $\phi_{i}=\theta_{i}\pm C$ axis is evident. This type of flow follows from a conserved quantity in the model: $\langle\dot{\phi}\rangle=\langle\dot{\theta}\rangle=0$, which can be seen by averaging equations (3) and (4) over the population. There are also oscillations in the radial position, where each swarmalator travels from the inner rim to the outer rim and back, in one orbit around the annulus.

This new, and final, state is similar to the double milling states found in biological swarms carrillo2009double , where populations split into counter-rotating subgroups. It is also similar to the vortex arrays formed by groups of sperm riedel2005self , where the angular position $\phi$ of each sperm is correlated with the phase $\theta$ associated with the rhythmic beating of its tail.

At the density level, the state is like a blurred version of the static phase wave, insofar as the spatial angle and phase of a given swarmalator are roughly correlated, as evident in Fig. 4(d). However unlike the static phase wave, the swarmalators are non-stationary. To highlight this difference, we name this state the ‘active phase wave’.

Order parameters. Having described the five states of our system, we next discuss how to distinguish them. We define the following order parameter,

[TABLE]

where $\phi_{i}:=\tan^{-1}(y_{i}/x_{i})$. As shown in Fig. 6, the magnitude $S_{\pm}$ varies from $1$ to [math] as we decrease $K$ from 0, passing through all the states in the upper left quadrant of the $(J,K)$ plane. (Note that all states except for static sync occur in this part of parameter space, so we hereafter confine our attention to just this region.)

To see why $S_{\pm}$ varies in this manner, recall that in the static phase wave, the spatial angle and phase of each swarmalator are perfectly correlated, $\phi_{i}=\pm\theta_{i}+C_{1}$ (recall that the $\pm$ and $C_{1}$ are determined by the initial conditions. This means either $S_{+}$ or $S_{-}$ is non-zero). Therefore $S_{\pm}=1$ at $K=0$, where the static phase wave state is realized. Moving into the $K<0$ plane we encounter the splintered phase wave. Here the correlation between $\phi_{i}$ and $\theta_{i}$ is not perfect, and so $S_{\pm}<1$. As $K$ is decreased the decay of this correlation is non-monotonic, which induces a dip in $S_{\pm}$ as seen in Fig. 6. Once the active phase wave is reached however this non-monotonicity disappears. As a result $S_{\pm}$ declines uniformly until it finally drops to zero when the static async state is reached, in which $\phi_{i}$ and $\theta_{i}$ are fully uncorrelated.

To sum up, $S_{\pm}$ is zero in the static async state, bifurcates from zero at a critical coupling strength $K_{c}$, is non-zero in the non-stationary splintered and active phase wave states, and is one in the static phase wave state.

Notice however that since $S_{\pm}$ is non-zero for both the splintered and active phase wave, it cannot distinguish between these states. To do this, we use another order parameter $\gamma$. We define this to be the fraction of swarmalators that have executed at least one full cycle in phase and position, after transients have been discarded. Then $\gamma$ is zero for the splintered phase wave, and non-zero for the active phase wave. Using $\gamma$ in concert with $S_{\pm}$ then allows us to discern all the macroscopic states of our system as illustrated in Fig. 6.

Stability analysis. To calculate the critical coupling strength $K_{c}$ at which the static async state loses stability, we consider perturbations $\eta$ in density space defined by

[TABLE]

where $\rho_{0}(\mathbf{x},\theta,t)=(4\pi^{2})^{-1}g_{1}(r)$ is density in the static async state. In Supplementary Note 2, we substitute this ansatz into the continuity equation, expand $\eta$ in a Fourier series, $\eta(\mathbf{x},\theta,t)=\sum_{n=0}b_{n}(\mathbf{x},t)e^{in\theta}+c.c.$, and derive evolution equations for the harmonics $b_{n}$. We show the critical mode is $b_{1}(\mathbf{x},t)$ (higher modes are zero, and the zeroth mode is stable) which obeys

[TABLE]

We next expand in an additional Fourier series: $b_{1}(r,\phi,t)=\sum_{m=0}^{\infty}f_{m}(r,t)e^{im\phi}+c.c.$. Substituting this ansatz into (14) leads to a evolution equation for each mode $f_{m}(r,t)$. We then set $f_{m}(r,t)=e^{\lambda_{m}t}c_{m}(r)$ and derive the following eigenvalue equation:

[TABLE]

where $R$ is the radius of the support of the density in the static async state. We focus first on the zeroth mode $f_{0}$ for which we can compute $H_{0}(r,s)$ analytically:

[TABLE]

where $\mathcal{K},\mathcal{E}$ are the complete elliptic integral of the first and second kinds. We were unable to solve (16) for $\lambda_{0}$ analytically. Instead, we found it numerically by approximating the integral using gaussian quadrature. This reduces (16) to the form $\lambda^{\prime}c_{i}=M_{ij}c_{j}$ where $M_{ij}=H_{0}(r_{i},r_{j})w_{j}$, $w_{j}$ are gaussian quadrature weights, $r_{i}=i*(R/N^{\prime})$ and $i=1\dots N^{\prime}$. The eigenvalues $\lambda_{0}^{\prime}(N^{\prime})$ of $M_{ij}$, which depend on the number of grid points $N^{\prime}$ used in the quadrature, then approximate $\lambda_{0}$.

The eigenvalues $\lambda_{0}^{\prime}(N^{\prime})$ have unexpected properties. The real part of the most unstable eigenvalue, denoted $\lambda^{*}_{0}\prime(N^{\prime})$, is positive for all $J,K$. This tells us that $f_{0}$ is always unstable, which in turn tells us that the static async state is always unstable! In Fig. 7 we plot $\lambda^{*}_{0}\prime(N^{\prime})$ versus $K$ for $J=0.5$ and $N^{\prime}=200$ grid points. As can be seen it is small but positive for sufficiently negative $K$. Note however that there is a transition-like point $K^{*}_{0}\approx-0.5$ beyond which $\lambda^{*}_{0}\prime(N^{\prime})$ increases sharply. Figure 7 also shows $\lambda^{*}_{m}\prime(N^{\prime})$ for $m=1,2,3,4$, which have the same behavior as $\lambda^{*}_{0}\prime(N^{\prime})$: they are small but positive for $K<K_{m}^{*}$, and grow sharply for $K>K_{m}^{*}$.

Small but positive eigenvalues for $K<K^{*}_{m}$ were a surprise. We were expecting them to be negative, since simulations show the static async state is stable. We were thus suspicious of these results, and doubted the accuracy of the approximation $\lambda^{*}_{0}\prime(N^{\prime})$ to the true $\lambda_{0}^{*}$. We therefore repeated the calculation for different values of $N^{\prime}$ up to $N^{\prime}=1600$ in Supplementary Note 2. Contrary to our expectations, we found that while the $\lambda^{*}_{m}\prime(N^{\prime})$ got smaller, they consistently remained positive for $K<K^{*}_{m}$.

We also crudely investigated the $N^{\prime}\rightarrow\infty$ limit by (i) fitting our data to curves of the form $a+b(N^{\prime})^{c}$ and (ii) using Richardson extrapolation. Due to the small magnitudes of the $\lambda^{*}_{m}\prime(N^{\prime})$ however, the results were rather unconvincing. Typical values for the best fit parameter $a$, which represents the limiting behavior of $\lambda_{m}^{*}\prime$, were $a\sim 10^{-6}$. The confidence interval for this parameter also contained positive and negative values. On top of that the approximations from methods (i) and (ii) gave inconsistent results. Hence we were unable to reliably determine the sign of $\lambda_{m}^{*}\prime$ when $K<K^{*}_{m}$ and $N^{\prime}\rightarrow\infty$, which preventing us from accurately ascertaining the stability of the static async state. We restate however that the fact that $\lambda_{m}^{*}\prime>0$ for the large but finite value of $N^{\prime}$ we used is significant evidence that the unanticipated instability of the static async state is genuine.

While a rigorous determination of the sign of $\lambda_{m}^{*}\prime$ when $K<K^{*}_{m}$ remains elusive, our analysis certifiably shows its magnitude is very small. Hence, whatever the stability or instability of the $m$-th mode $f_{m}$ turns out to be, it must be weak. In turn, then, the static async state has weak stability properties for $K<K_{c}$, where $K_{c}=\min_{m}K^{*}_{m}$ (i.e., at the point the most unstable $f_{m}$ loses stability). How can we find this $K_{c}$? In Fig. 7 we see the $f_{1}$ becomes unstable first. There are of course an infinite number of modes, but as can be seen, $\lambda^{*}_{m}$ appears to decrease with increasing $m$. Thus we assume $\min_{m}K^{*}_{m}=1$. In Supplementary Note 2 we approximate $K^{*}_{1}=\arg\!\max\frac{d^{2}\lambda_{1}^{*}\prime}{dK^{2}}$, calculate it for different $J$, and find the following linear relation:

[TABLE]

Summarizing our main result: in the continuum limit $N\rightarrow\infty$, the static async state is unstable for $K>K_{c}$, and either weakly stable, neutrally stable, or weakly unstable for $K<K_{c}$. Further, numerical evidence suggests that the third option, weakly unstable, is the most likely. While this result is perhaps unsatisfying from a technical perspective, in practice it has utility. For example as shown in Fig. 1, the approximation (18) for $K_{c}$ agrees reasonably well with finite $N$ simulations.

Genericity. Our analysis so far has been for the instance (3), (4), of the model defined by (1), (2). This begs the question: are the phenomena we found generic to the model? Or specific to this instance of the model? To answer this question, we ran simulations for different choices of the functions $\mathbf{I}_{\mathrm{rep}},\mathbf{I}_{\mathrm{att}}$ and $G$; see Supplementary Note 4.

In all but one case, we found the same phenomena. The exception is when a linear attraction kernel $\mathbf{I}_{\mathrm{att}}(\mathbf{x})=\mathbf{x}$ is used. Here we found new states, which we call ‘non-stationary phase waves’. They are similar to the active phase wave, except now the phase $\Psi_{\pm}$ of the order parameter $W_{\pm}$ begins to rotate, reminiscent of the traveling wave states found in the Kuramoto model with distributed coupling strengths hong2011kuramoto ; hong2016phase . We further discuss this and other properties in Supplementary Note 4.

Noise and disordered natural frequencies. The swarmalators previously considered were identical and noiseless. We now relax these idealizations. Then the governing equations are

[TABLE]

where $\omega_{i}$ are random variables drawn from a Lorentzian g(\omega)=(\sigma/\pi)\Big{[}(\omega-\mu)^{2}+\sigma^{2}\Big{]}^{-1}. By a change of frame we set $\mu=0$, leaving just $\sigma$, which quantifies the strength of the disorder. We choose white noise variables $\eta_{i}(t)$ and $\xi^{\mathbf{x}}_{i}(t)$ with zero mean and strengths $D_{x},D_{y},D_{\theta}$ characterized by $\langle\xi^{x}_{i}(t)\xi^{x}_{j}(t^{\prime})\rangle=2D_{x}\delta_{ij}\delta(t-t^{\prime})$, etc.

Simulations show that when just phase noise $D_{\theta}$ is turned on, noisy versions of all the states are realized. The splintered phase wave however degenerates into the active phase wave for all but the smallest noise $D_{\theta}\gtrsim 10^{-3}$. In the remaining states, the spatial densities remain compact supported with the same radii, except now the swarmalators have noisy phase motion (this induces some spatial movement, which disappears when $N\rightarrow\infty$ as we show in Supplementary Note 3). Hence the following states, where we have swapped the descriptor ‘static’ with ‘noisy’, are robustly realized when $D_{\theta}>0$: noisy phase wave, active phase wave, noisy async.

Frequency disorder $\sigma>0$ has a more serious effect. Since $g(\omega)$ is symmetric about zero, there are equal numbers of swarmalators with oppositely signed natural frequencies. This turns the static/noisy phase wave into the active phase wave, in the sense that counter-rotating groups develop. This is not seen in the async state. Here, there are noisy spatial movements which vanish as $N\rightarrow\infty$, as in the noisy async state. In contrast however, the swarmalators execute noisy, but full, phase cycles. To highlight this distinction, we rename the state the active async state. The states realized are then the active phase wave and the active async.

Finally spatial noise $D_{x},D_{y}>0$ simply blurs the spatial densities of the states. No other phenomena are induced. Hence when $D_{\theta},\sigma,D_{x},D_{y}>0$, we again get the active phase wave and active async states.

In Fig. 8 we plot the order parameter $S(K)$ for different amounts of noise and frequency disorder. As for the original model, $S$ simply declines to zero as $K$ is decreased, with the noise and disordered frequencies changing just the shape of the curves and the value of $K_{c}$. Note the disappearance of the dip in $S$ for small $K$, which indicates the absence of the splintered phase wave state. Note also we do not plot the second order parameter $\gamma$ which discerns the splintered phase wave since this state does not robustly exist when $\sigma,D_{\theta}\neq 0$.

Swarmalators in 3D. So far we have considered swarmalators moving in two dimensions. While there are physical systems where this approximation is valid, such as certain active colloids pohl2014dynamic or sperm, which are often attracted to two dimensional surfaces maude1963non , this restriction was mostly for mathematical convenience. Here we explore the more physically realistic case of motion in three spatial dimensions (in Supplementary Note 6 we also explore motion in one dimension). For simplicity we consider the case of identical swarmalators with no noise, although we relax these idealizations in Supplementary Note 5. Our system is then

[TABLE]

where $\mathbf{x}_{i}=(x_{i},y_{i},z_{i})$. These are the same as equation  (3) and (4), except the exponent of the hard shell repulsion is now $3$ (we choose this because it yields simple formulas for the radii of certain states).

Simulations show that analogues of the states found in 2D are realized. We show these as scatter plots in the $(x,y,z)$ plane in Fig. 9. We also provide movies of the evolution to these states in Supplementary Movies 9-12. The static sync and async states become spheres (note we do not plot the static sync state due to space limitations) as seen in panel (a). As in the 2D case, we can calculate their radii when a linear attraction kernel is used,

[TABLE]

which agree with simulation as shown in Supplementary Note 5.

In panel (b) we show the static phase wave becomes a sphere with a cylindrical hole through its center. The orientation of this cylinder is determined by the initial conditions. The phase and azimuthal angle $\phi=\tan^{-1}(y/x)$ are correlated in the same way for each value of the polar angle $\alpha=\cos^{-1}(z/\sqrt{x^{2}+y^{2}+z^{2}})$ (when the azimuthal and polar angles are measured relative to the axis of the cylindrical hole). We show this more clearly in a scatter plot in the $(\theta,\phi)$ plane in Supplementary Note 5 .

As in the 2D model, this correlation between $\phi$ and $\theta$ persists for the splintered phase waves and active phase wave states as can be seen in panels (c) and (d) of Fig. 9. The motion of the swarmalators in these states are as before: in the splintered phase wave they ‘quiver’, executing small oscillations in space and phase, while in the active phase wave they execute full rotations (note the spatial component of these rotations are in the azimuthal direction $\hat{\phi}$ only, not in the polar direction $\hat{\alpha}$). In Supplementary Note 5 we show how the order parameters $S_{\pm},\gamma$ can also be used to differentiate these 3D states.

Alignment and self-propulsion. Up to now we have considered the trivial case of swarmalators that propel themselves with constant magnitude and direction, in a manner uninfluenced by their neighbors. This allowed us to set this term to zero via a change of reference. In many real systems, however, such behavior is unrealistic: individuals often adjust the direction of their motion to align with that of their neighbors. Vicsek studied this alignment effect in a seminal work vicsek1995novel .

We here partially explore the effect of alignment on swarmalator systems. Accordingly we endow each swarmalator with an orientation $\beta$, which characterizes the direction of its self-propulsion. The inclusion of alignment makes our model complicated; there are now four state variables $(x,y,\theta,\beta)$ per swarmalator, which could interact with each other in potentially many ways. Furthermore, there are six parameters $(J,K,\sigma,D_{\theta},D_{x},D_{y})$, not to mention any additional parameters governing the evolution of $\beta$. An exhaustive study of orientable swarmalators is thus beyond the scope of the present work. Hence, we restrict ourselves to answering a simple question: are the states of our swarmalator system robust to the inclusion of simple alignment dynamics?

To this end, we study the simplest possible extension to our current model: we choose Vicsek type interactions between $\mathbf{x}$ and $\beta$, and leave $\beta$ and the phase $\theta$ uncoupled (although they are indirectly coupled through the position $\mathbf{x}$). Our system then reads

[TABLE]

where $\hat{n}=(\cos{\beta},\sin{\beta})$, $\Lambda_{i}$ is the set of swarmalators within a distance $\delta$ of the $i$-th swarmalator, and $|\Lambda_{i}|$ is the number of such neighbors. The $\zeta_{i}(t)$ is a white noise variable with zero mean and strength $D_{\beta}$ characterized by $\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=2D_{\beta}\delta_{ij}\delta(t-t^{\prime})$.

Simulations show that for certain parameter values aligned versions of all our states persist. We plot two of these in panels (a) and (b) of Fig. 10, where each swarmalator is depicted as a colored arrow, oriented according to $\beta$, and colored according to phase. As can be seen the swarmalators are aligned, with their space-phase distributions being the same as before. In contrast to the original model, however, the center of mass of each distribution now moves (in a direction determined by the initial conditions). In this sense the states are mobile. They are however equivalent to their static versions via a change of reference frame, $\mathbf{x}\rightarrow\mathbf{x}+\mathbf{v}_{0}t$. For larger $D_{\beta}$, the unaligned versions of the same states are realized, as illustrated in panels (c) and (d) of Fig. 10.

We have demonstrated that the phenomena of our system are insensitive to the inclusion of simple alignment dynamics. We restate however that we have not comprehensively explored the space defined by the other parameters $(J,K,\sigma,v_{0},D_{x},D_{y})$ given its large size. Thus it remains to be seen if new states will be found.

Discussion

We have examined the collective dynamics of swarmalators. These are mobile particles or agents with both phase and spatial degrees of freedom, which lets them sync and swarm. Furthermore, their phase and spatial dynamics are coupled. By studying simple models, we found this coupling leads to rich spatiotemporal patterns which we explored analytically and numerically. These patterns were robust to modifications to the model, namely motion in one, two, and three spatial dimensions, distributed natural frequencies, noisy interactions, and alignment dynamics. We thus believe they could be realized in nature or technology.

A pertinent future goal, then, is to investigate the behavior of real-world systems of swarmalators. As mentioned in the introduction, colloidal suspensions of magnetic particles yan2012linking ; snezhko2011magnetic ; martin2013driving or active spinners nguyen2014emergent ; van2016spatiotemporal are promising candidates. For example, structures equivalent to the static phase wave state have been experimentally realized by Snezhko and Aranson, when studying the behavior of ferromagnetic colloids at liquid-liquid interfaces snezhko2011magnetic (the particles comprising the colloids can be considered swarmalators if we interpret the angle subtended by their magnetic dipole vectors as their phase). As shown in Fig. 4 of snezhko2011magnetic , the colloids can form asters. These are structures composed of radial chains of magnetically ordered particles, which “decorate slopes of a self-induced circular standing wave” snezhko2011magnetic , analogous to the annular pattern of correlated phases and positions of the static phase wave shown in Fig. 2(c).

Perhaps colloidal equivalents of the splintered and active phase wave states could also be realized. Aside from being theoretically interesting, the ability to engineer these states could have practical application. For instance, Snezhko and Aranson also show that asters can be manipulated to capture and transport target particles. The non-stationary behavior of the splintered and active wave states might also have locomotive utility. Tentative evidence for this claim is provided by populations of cilia, whose collective metachronal waves, similar to the motion of swarmalators in the aforementioned states, are known to facilitate biological transport elgeti2013emergence ; okada1999abnormal ; wong1993nature .

Other plausible systems of real-world swarmalators are biological microswimmers, self-propelled micro-organisms capable of collective behavior lauga2009hydrodynamics . One such contender is populations of spermatoza, which exhibit rich swarming behavior such as trains immler2007hook ; moore2002exceptional and vortex arrays riedel2005self , the latter of which is reminiscent of the active phase wave state, as mentioned in the Results section. The phase variable for each sperm is associated with the rhythmic beating of the sperm’s tail, which can synchronize with that of a neighboring sperm taylor1951analysis ; fauci1995sperm . It has been theorized that this can induce spatial attraction yang2008cooperation , leading to clusters of synchronized sperm, consistent with experimentally observed behavior hayashi1996insemination .

There are also theoretical avenues to explore within our proposed model of swarmalators. For instance the curious stability properties of the static async state deserve further study. Another route would be to include more realism by including heterogeneity in the coupling parameters $K,J$, or by choosing more complicated interaction functions $\mathbf{I}_{\mathrm{att}},\mathbf{I}_{\mathrm{rep}},G,H_{\mathrm{att}}$. For example we chose $H_{\mathrm{att}}(\theta)=\sin(\theta)$ to mimic the Kuramoto model, but as we saw, it led to just the trivial static sync state when $K>0$. Perhaps choosing the more realistic Winfree model for the phase dynamics, which gives rise to richer collective behavior, would lead to more interesting swarmalator phenomena in this parameter regime.

Perhaps the most important direction for future work is to more fully explore the interplay among aggregation, alignment, and synchronization—or put another way, to explore the collective behavior of particles with a position $\mathbf{x}$, an orientation $\beta$, and an internal phase $\theta$. The primary goal of our work is to draw attention to this class of problems, which we believe define a wide landscape of emergent behavior. In this work, we have started to map out this landscape by studying a simple model that contains a subset of these three effects, namely aggregation and synchronization.

Others have considered the remaining subsets. For example, Leon and Liverpool have explored the interaction between alignment and synchronization leoni2014synchronization . They introduced a new class of soft active fluids whose units have an orientation and phase. They found this mixture can either enhance or inhibit the transition from disordered states to states with polar order. The latter states are roughly similar to the aligned static async states. They also found transitions from disordered states to states with phase order, which are analogous to unaligned static sync states. Yet counterparts of the static, splintered, and active phase waves were not reported.

The final combination, aggregation and alignment, is perhaps the most well studied, in both new models and old. For instance, Starnini et al. starnini2016emergence recently introduced a model of mobile particles capable of aggregating and aligning their opinions, and found the emergence of echo chambers. Even in the classic Vicsek model and its numerous extensions, new phenomena are still being found. For instance, Kruk et al. found that delayed alignment in the Vicsek model produces self-propelled chimeras kruk2015self ; perhaps delayed phase interactions could lead to similar states for swarmalators. Liebchen and Levis liebchen2016rotating considered units with an intrinsic rotation, and found ‘phase separated droplets’: clusters of rotation-synchronized particles surrounded by a sea of incoherent particles (multiple droplets are also possible). These droplets are similar to our static sync states, but they differ in the crucial respect that the entire population is synchronized in our static sync state. Here too, the counterparts of our static, splintered, and active phase waves were not seen.

Thus, to the best of our knowledge, no other models display states analogous to the splintered phase waves and active phase waves found in our swarmalator model. In that sense, those two states are unprecedented.

Methods

Properties of static sync and async state. We here use techniques used by Fetecau et al. fetecau2011swarm when studying swarm dynamics to study the static sync and static async states with a linear attraction kernel $\mathbf{I}_{\mathrm{att}}(\mathbf{x})=\mathbf{x}$ is used. We start with the async state whose density is

[TABLE]

We wish to solve for the radial density $g_{1}(r)$ and the radius $R$ of its support. In this state the swarmalators are at rest and their phases are unchanging, so $\underline{v}\equiv\underline{0}$, where $\underline{v}=(v_{x},v_{y},v_{\theta})$. As we will show, it is also useful to consider the divergence of the velocity, which must also be zero (from the continuity equation for the conservation of swarmalators, and by applying the assumptions that the density for the static async state is stationary and the velocity is zero). This gives us a pair of simultaneous equations,

[TABLE]

We begin with divergence term given by (30). In cartesian coordinates the velocity reads

[TABLE]

The divergence has a spatial and phase component: $\nabla.\underline{v}=\nabla_{\mathbf{x}}.\mathbf{v}_{\mathbf{x}}+\partial_{\theta}v_{\theta}$. The phase component $\partial_{\theta}v_{\theta}$ is trivially zero, since the swarmalators’ phases are uniformly distributed in phase in the static async state. We find the spatial component by applying $\nabla_{\mathbf{x}}$ to (31):

[TABLE]

Here we have used the identity (expressed most cleanly in cartesian coordinates)

[TABLE]

Simplifying this, and substituting $\partial_{\theta}v_{\theta}=0$ gives the full divergence

[TABLE]

By (30) we require this to be zero, which gives a self-consistent equation for $\rho$:

[TABLE]

Finally substituting the ansatz given by (28) into this and performing the integration over $\phi$ gives

[TABLE]

This tells us $\rho$ is constant inside a disc of radius $R$. The radius $R$ can be determined via self-consistency: $M=\int_{0}^{R}rg(r)dr=\int_{0}^{R}rMdr\Rightarrow R=1$. By normalizing $\rho$ as per (28) we find $M=1$ which means $g(r)=2$. Putting this all together gives

[TABLE]

We must now check if the solutions given by (39) and (40) imply $\underline{v}\equiv\underline{0}$ as required by (29). We do this in cartesian coordinates, in which

[TABLE]

where $\theta_{0}$ is the final, common phase of each swarmalator. Substituting this into equations (31) and (32) for $v_{\mathbf{x}},v_{\theta}$ and performing the integration gives

[TABLE]

where we have used the identity

[TABLE]

We see that $\mathbf{v}_{\mathbf{x}}=0$ at $\theta=\theta_{0}$ if $R=1$, as required. Hence we have shown that the solutions (39), (40) satisfy equations (29) and (30).

Carrying out the same analysis for the static sync state leads to

[TABLE]

where $\theta_{0}$ is the final common phase of the swarmalators in the static sync state.

Properties of static phase wave state. Here we calculate the density of swarmalators, and inner and outer radii $R_{1},R_{2}$ of the annulus, in the static phase wave state, when a linear attraction kernel is used.

The calculation is the same as for the static sync and async states: we assert

[TABLE]

The density of the static phase wave state is

[TABLE]

We first calculate the divergence, which in polar coordinates is given by

[TABLE]

The velocity $\underline{v}=(v_{r},v_{\phi},v_{\theta})$ is given by

[TABLE]

[TABLE]

[TABLE]

Taking the derivatives on these, plugging in Eq. (49) for $\rho$, and substituting the result into (50), gives

[TABLE]

Setting this to zero, we see $g(r)$ satisfies

[TABLE]

which means it can be determined self-consistently in terms of $R_{1}$ and $R_{2}$. The result is

[TABLE]

with \Gamma_{J}=2J(R_{2}^{3}-R_{1}^{3})\Big{(}3J(R_{2}^{2}-R_{1}^{2})+12\Big{)}^{-1}.

Next we use the result (56) in $\underline{v}=\underline{0}$ to compute the inner and outer radii $R_{1},R_{2}$. We first evaluate $v_{r}$ by substituting (56) into (51). Performing the integration we get

[TABLE]

with

[TABLE]

Since $\underline{v}\equiv\underline{0}$, the coefficients $C,D$ must be zero. This yields two equations for $R_{1},R_{2}$, with solutions

[TABLE]

with

[TABLE]

and small-$J$ expansion given by

[TABLE]

Acknowledgments

Research supported by United States NSF Grant Nos. DMS-1513179 and CCF-1522054 (S.H.S), and by South Korean NRF Grant No. NRF-2015R1D1A3A01016345 (H.H).

Author contributions

K.P.O and S.H.S conceived the research. K.P.O and H.H performed the numerics. K.P.O. performed the analysis and drafted the manuscript. All authors discussed the results, drew conclusions and edited the manuscript.

Additional information

Competing financial interests. The authors declare no competing financial interests.

Data availability. The data that support the findings of this study (simulation source code and figure raw data) are available from K.P.O upon request. Email: kevin.p.okeeffe@gmail.com