Swarmalators with thermal noise

TL;DR

This paper studies a one-dimensional model of swarmalators influenced by thermal noise, revealing four collective states and providing a detailed phase diagram that aligns well with simulations, offering insights into thermal self-synchronizing and self-assembling systems.

Contribution

It introduces a tractable XY-type model of swarmalators with thermal noise, identifying multiple collective states and deriving an accurate phase diagram.

Findings

Identification of four distinct collective states.

Discovery of a mixed state blending two behaviors.

High-accuracy phase diagram matching simulations.

Abstract

We investigate a population of swarmalators, a mobile version of phase oscillators that both sync in time and swarm through space. We focus on a XY-type model of identical swarmalators running on a one-dimensional ring and subject to thermal noise. We uncover four distinct collective states, some of which capture the behavior of real-world swarmalators such as vinegar eels and sperm. Among these, the most intriguing is the `mixed state', which blends two of the other states. We present a comprehensive phase diagram from the Fourier mode analysis with a high accuracy, which is in excellent agreement with numerical simulation results. Our model serves as a tractable toy model for thermal systems that both self-synchronize and self-assemble interdependently.

Equations

$$\dot{x_i}$$

$$\dot{{\theta }_i}$$

$$⟨{\xi }_i^x(t){\xi }_j^x(t^{\prime })⟩$$

$$⟨{\xi }_i^{\theta }(t){\xi }_j^{\theta }(t^{\prime })⟩$$

$$H=-\frac{J}{2N}\sum _{i,j=1}^N\cos (x_j-x_i)\cos ({\theta }_j-{\theta }_i),$$

$$Z_x$$

$$Z_{\theta }$$

$$W_{+}$$

$$W_{-}$$

$$X_i\equiv x_i+{\theta }_i,\text{and}{Y}_i\equiv x_i-{\theta }_i,$$

$${\dot{X}}_i$$

$${\dot{Y}}_i$$

$${\xi }_i^X={\xi }_i^x+{\xi }_i^{\theta },\text{and}{\xi }_i^Y={\xi }_i^x-{\xi }_i^{\theta },$$

$$⟨{\xi }_i^X(t){\xi }_j^X(t^{\prime })⟩$$

$$⟨{\xi }_i^Y(t){\xi }_j^Y(t^{\prime })⟩$$

$$⟨{\xi }_i^X(t){\xi }_j^Y(t^{\prime })⟩$$

$$D_X=D_Y=D_x+D_{\theta }\equiv D,$$

$$D_{XY}=D_x-D_{\theta }\equiv d.$$

$$\frac{\partial \rho }{\partial t}=-\frac{\partial }{\partial X}{\mathcal{J}}_X-\frac{\partial }{\partial Y}{\mathcal{J}}_Y$$

$${\mathcal{J}}_X$$

$${\mathcal{J}}_Y$$

$$F_{+}\equiv S_{+}\sin ({\Phi }_{+}-X),F_{-}\equiv S_{-}\sin ({\Phi }_{-}-Y).$$

$$\rho (X,Y,t)=\frac{1}{(2\pi {)}^2}\sum _{n,m=-\infty }^{\infty }{\alpha }_{n,m}(t)e^{in(X-{\Phi }_{+})}{e}^{im(Y-{\Phi }_{-})},$$

$$S_{+}$$

$$S_{-}$$

$${\dot{\alpha }}_{n,m}-i(n{\dot{\Phi }}_{+}+m{\dot{\Phi }}_{-}){\alpha }_{n,m}$$

$$=-[D(n^2+m^2)+2dnm]{\alpha }_{n,m}$$

$$+\frac{1}{2}(n{J}_{+}+m{J}_{-})S_{+}({\alpha }_{n-1,m}-{\alpha }_{n+1,m})$$

$$+\frac{1}{2}(n{J}_{-}+m{J}_{+})S_{-}({\alpha }_{n,m-1}-{\alpha }_{n,m+1}).$$

$$0=-Dn{\alpha }_{n,0}+\frac{1}{2}{J}_{+}{S}_{+}({\alpha }_{n-1,0}-{\alpha }_{n+1,0}),$$

$${\alpha }_{n,o}=\frac{I_{|n|}\left(\frac{J_{+}{S}_{+}}{D}\right)}{I_0\left(\frac{J_{+}{S}_{+}}{D}\right)},$$

$$S_{+}\approx \sqrt{2(J_{+}/2D-1)},$$

$${\alpha }_{n,m}=\frac{I_{|n|}\left(\frac{JS}{D}\right)I_{|m|}\left(\frac{JS}{D}\right)}{[I_0\left(\frac{JS}{D}\right){]}^2}.$$

$${\dot{c}}_{1,0}=-D\left(1-\frac{J_{+}}{2D}\right)c_{1,0},$$

$${\dot{c}}_{0,1}=-D\left(1-\frac{J_{+}}{2D}\right)c_{0,1},$$

$${\dot{c}}_{n,m}=-\left[D(n^2+m^2)+2dnm\right]c_{n,m}\text{for other}(n,m).$$

$${\dot{\alpha }}_{n,m}=-[u(n^2+m^2)+2znm]{\alpha }_{n,m}$$

$$+(n+vm){\alpha }_{1,0}({\alpha }_{n-1,m}-{\alpha }_{n+1,m})$$

Full text

I Introduction

Swarmalators, short for “swarming” oscillators, are mobile phase oscillators that can both synchronize in time and swarm through space ref:swarmalators . They were recently introduced to model the numerous systems in which synchronization and swarming co-exist and interact. Examples are sperm ref:sperms , vinegar eels ref:vinegareels , magnetotactic bacteria, starfish embroys ref:livingcrystals , Japanese Tree frogs ref:Japanesefrogs , Janus particles ref:Janusparticles , Quincke rollers ref:Quinckerollers , and robotic swarms ref:robots .

Research on swarmalators began to gain traction about 15 years ago when Iwasa and Tanaka introduced a universal model of chemotactic oscillators which produced diverse states ref:chemo ; ref:diverse . O’Keeffe, Hong, and Strogatz later introduced a generalized Kuramoto model which produced five collective states that are commonly observed in nature ref:swarmalators . This swarmalator model is being used as a spring-board to further study the collective behavior of swarmalators. The effects of pinning ref:pinning , forcing ref:forcing , and various types of coupling (delayed ref:delay , finite-range ref:finiterange stochastic ref:stochastic , mixed sign ref:coupling ; ref:gourab ; ref:adv ; ref:pre-coupling ) as well as other phenomena/scenarios ref:ring ; ref:mfl ; ref:topological ; ref:jimenz ; ref:ha ; ref:steven have been explored. For reviews of swarmalators, see ref:review1 ; ref:review2 .

Theoretical results on swarmalators are scarce due to the complexity of the system, which is characterized by nonlinear couplings of numerous degrees of freedom. For a system of $N$ particles, each with position $x\in\mathbb{R}^{d}$ and phase $\theta\in S^{1}$, the analysis of the nonlinear coupled ordinary differential equations (ODEs) is challenging. What is missing in the swarmalator field is the “right” or natural toy model such as the Ising model for equilibrium phase transition study ref:Ising or the Kuramoto model for synchronization studies ref:Kuramoto , which captures the essential ingredient while remaining simple enough to solve.

This paper focuses on a simple one-dimensional ($1$D) swarmalator model ref:1D ; ref:1Dwithfrequency which could fill this ‘modeling gap’ in studies of swarmalation, as Verberk ref:verberck has dubbed the interplay between swarming and synchronizing oscillations. The model confines the swarmalators to run on a 1D ring, which makes it analytically tractable yet it still produces real-world behavior; its collective states capture the behavior of real-world 1D swarmalators, such as bordertaxic vinegar eels and sperm, and also capture the rotational analogue of 2D & 3D swarmalators, such as magnetic Janus particles and dielectric Quinke rollers.

In a recent study ref:1Dwithfrequency , the 1D swarmalator model with “quenched” disorder was explored, wherein the natural frequencies of each swarmalator were chosen randomly and then fixed for all times (nonidentical oscillators). The authors utilized a “toroidal” Ott-Antonsen(OA) ansatz on a special submanifold of state space to identify four distinct long-lived states ref:1Dwithfrequency . Meanwhile, systems in nature are typically influenced by their surrounding environment, corresponding to a heat bath, which implies that natural systems are stochastic and affected by thermal noise. In this context, it is pertinent to investigate the effects of “temporal/thermal” disorder on the collective properties of the system.

To this end, we introduce thermal noise into the 1D swarmalator model and explore its collective behavior. Here, we restrict our analysis to the case of identical swarmalator only, leaving the generalization to nonidentical ones with thermal noise for future study, as this would significantly complicate the analysis. We note that our case with identical oscillators can be also interpreted as a coupled XY model for a magnetic spin system ref:XY , which provides various interesting nonequilibrium steady states.

We analyse the model both numerically and analytically. It is well known that the analytic OA ansatz fails in the presence of thermal noise ref:OA ; ref:KS ; ref:HJHP ; ref:Pikovsky . Consequently, the Fourier mode analysis does not close with a finite number of terms, which leads to an infinite hierarchy of coupled equations. Nevertheless, as the contributions from higher modes become increasingly negligible to the synchronization order parameters, the coupled equations can be suitably truncated and solved numerically to evaluate the order parameters highly accurately. Our results compare very well with numerical simulations by integrating the equations of motion directly. We also analytically derive the instability condition for the completely disordered phase and employ a perturbation theory in order to reach out to a partially ordered phase. Our analysis is consistent with numerical results.

II Swarmalators on a 1D ring

We consider a population of $N$ coupled identical swarmalators on a 1D ring, of which the dynamics are governed by

[TABLE]

where $x_{i}$ and $\theta_{i}$ denote the position and phase of the $i$th swarmalator, respectively, and are both periodic with a period of $2\pi$. One may add an identical natural frequency $\omega$ for phase in Eq. (2), but it can be eliminated by taking a simple transformation of $\theta_{i}\rightarrow\theta_{i}+\omega t$. Likewise, the identical directed velocity $v$ can be also set to zero without loss of generality. $\xi_{i}^{x}$ and $\xi_{i}^{\theta}$ represent thermal noise with zero mean, and characterized by

[TABLE]

where $D_{x}$ and $D_{\theta}$ denote the “temperature” of the heat bath for $x$ and $\theta$ variables, respectively. The two noises are independent of each other, i.e. $\langle\xi_{i}^{x}\xi_{j}^{\theta}\rangle=0$. $J$ and $K$ are the coupling parameters.

This 1D swarmalator model has been previously studied in the absence of thermal noises for identical swarmalators ref:1D and for random nonidentical ones ref:1Dwithfrequency . The novelty of this paper lies in the presence of thermal noises without quenched disorder for nonidentical swamalators.

Equation (2) models position dependent synchronization, wherein the sine terms induce the swarmalators to reduce their phase difference for $K>0$, while the cosine term implies that the synchronization tendency is stronger for oscillators which are closer in space. Equation (1) is the mirror-image of Eq. (2), which models phase-dependent swarming or aggregation for $J>0$. This can be thought as synchronization occurring on a torus, as both position and phase are circular variables $(x_{i},\theta_{i})\in\mathbb{T}^{1}$. In this sense, it can be referred to as a generalized version of the thermal Kuramoto model ref:KS for two coupled populations of identical oscillators.

This model can be also seen as a nonequilibrium generalization of the widely-recognized XY model for a magnetic spin system ref:XY , wherein spins possess two degree of freedoms, $x$ and $\theta$, that are coupled. For the special case of $J=K$ and $D_{x}=D_{\theta}$, the dynamics described by Eqs. (1) and (2) is equivalent to that of an equilibrium (EQ) system in contact with a thermal reservoir of temperature $D_{x}$ with the Hamiltonian

[TABLE]

where an EQ (magnetic) phase transition is expected as $J$ varies. For $J\neq K$, the Hamiltonian structure is lost and the system can not relax to an equilibrium state. This is a typical nonequilibrium (NEQ) situation in the presence of external non-conservative forces, where various NEQ (magnetic) phases can emerge. For $D_{x}\neq D_{\theta}$, this model may describe a heat engine under the thermal gradient.

Yet another way to interpret the model is as a 1D Vicsek-type model ref:Vicsek where $\theta$ corresponds to an orientation, as opposed to an internal phase variable.

III Order parameters

In order to gain insight into the order parameters characterizing collective behaviors, numerical simulations are carried out by integrating Eqs. (1) and (2) for various parameter values. Four distinct steady states were identified, and their typical states are visualized in Fig. 1; (a) uniformly distributed, (b) linearly correlated between $x$ and $\theta$, (c)mixed between (b) and (d), and (d) two well-defined isotropic clusters separated by a distance of $\pi$ in both directions.

One may consider the Kuramoto synchronization (or XY magnetic) order parameters ref:Kuramoto for each variable $x$ and $\theta$ as

[TABLE]

where $R_{x}$ and $R_{\theta}$ measure the magnitude of synchronization (real and non-negative) in $x$ and $\theta$, respectively, and $\Theta_{x}$ and $\Theta_{\theta}$ are their mean angles (real). $\langle\cdots\rangle$ denotes the thermal (noise) average. However, they are not effective order parameters, as they vanish in all four phases. This is due to the special coupling structure in Eqs. (1) and (2), which are invariant under the “$\pi$” transform of $x_{j}\rightarrow x_{j}+\pi$ and $\theta_{j}\rightarrow\theta_{j}+\pi$ for any $j$, as noted in Refs. ref:1D and ref:1Dwithfrequency . It implies that the dynamics of any $(x,\theta)$ state is equivalent to that of the $(x+\pi,\theta+\pi)$ state. In particular, the number of swamalators in each cluster in (c) and (d) appears to be equal, as is expected for large $N$ in accordance with combinatorial theory. Thus, the contributions from each cluster cancel out, leading to $Z_{x}=Z_{\theta}=0$ in the $N\rightarrow\infty$ limit. Likewise, the same conclusion can be drawn for (b).

Thus, it is necessary to explore alternative order parameters that can be used to quantify the correlation between two variables, which were proposed previously in ref:1D and ref:1Dwithfrequency , as

[TABLE]

where $S_{\pm}$ is real and non-negative by definition and $\Phi_{\pm}$ is real. This particular form of the correlation measure can be inferred from the linearity of correlations between $x$ and $\theta$ with slopes of $\pm 1$, as seen in Fig. 1. With these order parameters, the four distinct steady states can be clearly distinguished as follows: (a) Disordered (incoherent) state with $(S_{+},S_{-})=(0,0)$; (b) Phase wave state with $(S,0)$ or $(0,S)$, where the value of $S$ depends on coupling parameters. The $(S,0)$ and $(0,S)$ states are equally probable, analogous to the magnetically ordered up and down states in the Ising model; (c) Mixed state with $(S_{1},S_{2})$ or $(S_{2},S_{1})$, where $S_{1}\neq S_{2}$ and both are finite. Finally, (d) Synchronized (ordered) state with $(S,S)$, where both correlations in Eqs. (8) and (9) are equal in magnitude.

IV Coordinate transformation

It is convenient to make coordinate transformations as

[TABLE]

which imply that the new variables $X$ and $Y$ are both periodic with a period of $4\pi$. Subsequently, the dynamic equations in Eqs. (1) and (2) can be simplified to

[TABLE]

with $J_{\pm}=(J\pm K)/2$, and $S_{\pm}$ and $\Phi_{\pm}$ are defined in Eqs. (8) and (9).

The new noises, $\xi^{X}_{i}(t)$ and $\xi^{Y}_{i}(t)$, are defined as

[TABLE]

that are not independent of each other in general. Simple calculations lead to $\langle\xi_{i}^{X}\rangle=\langle\xi_{i}^{Y}\rangle=0$ and

[TABLE]

where

[TABLE]

The correlation between the two noises $D_{XY}$ disappears when the two temperatures are equal ($D_{x}=D_{\theta}$),

It is worth noting that the dynamic equations in Eqs. (11) and (12) are invariant under the transformation of $X_{i}\rightarrow X_{i}+2\pi$ or $Y_{i}\rightarrow Y_{i}+2\pi$ for any $i$, implying that the variables $X$ and $Y$ can be treated as periodic with a period of $2\pi$, instead of $4\pi$, if the initial distributions possess the same periodicity.

When $J=K$ ($J_{-}=0$) with $D_{x}=D_{\theta}$ ($d=0$), the dynamic equations decouple completely, leading to the EQ XY model with two decoupled degree of freedoms for $X$ and $Y$, which can be solved exactly. By varying $J$, we find that the EQ magnetic transition occurs at $J_{c}=2D_{X}=4D_{x}$ from the disordered state $(0,0)$ (a), directly to the ferromagnetic ordered state $(S,S)$ (d) without an intermediate phase wave or mixed state. Intermediate states (b) and (c) may emerge only when the system is driven out of EQ, either by an external non-conservative force ($J\neq K$) or a thermal gradient caused by multiple reservoirs ($D_{x}\neq D_{\theta}$).

V Fourier mode analysis

In the $N\rightarrow\infty$ limit, the system state can be described by a continuous probability density function (PDF) $\rho(X,Y,t)$ at time $t$ in terms of new variables $X$ and $Y$. As discussed in the preceding section, these variables can be treated as periodic with a period of $2\pi$, thus their range is both set to $(0,2\pi)$ with the normalization condition as $\int_{0}^{2\pi}\int_{0}^{2\pi}dXdY\rho(X,Y,t)=1$.

The Fokker-Planck equation ref:Risken , derived from the dynamic equations of (11) and (12) with Eqs. (17) and (18), is given by

[TABLE]

where the probability currents ${\cal J}_{X}$ and ${\cal J}_{Y}$ are

[TABLE]

with

[TABLE]

The Fourier series of the PDF $\rho(X,Y,t)$ in $X$ and $Y$ with period $2\pi$ is given by

[TABLE]

where the extra factors $\Phi_{+}$ and $\Phi_{-}$ are inserted for convenience. As the PDF should be real, $\alpha^{*}_{n,m}=\alpha_{-n,-m}$ holds for all $(n,m)$ and the normalization condition yields $\alpha_{0,0}=1$.

The order parameters in Eqs. (8) and (9) can be expressed in terms of Fourier coefficients as

[TABLE]

It should be noted that, as $S_{\pm}$ is real, $\alpha_{1,0}$ and $\alpha_{0,1}$ must also be real. Substitution of Eq. (23) into the Fokker-Planck equation of Eq. (19) yields time-dependent mode-coupled equations as

[TABLE]

where $S_{+}=\alpha_{1,0}$ and $S_{-}=\alpha_{0,1}$. These equations form an infinite hierarchy of coupled equations, which can not be reduced to a finite number of equations unlike the Kuramoto model without thermal noise, where the OA ansatz holds.

VI steady state solutions

The steady-state solutions (fixed points) $\rho^{s}(X,Y)$ can be obtained by setting $\dot{\alpha}_{n,m}=\dot{\Phi}_{\pm}=0$ in Eq. (26). There exists the trivial solution $\alpha_{n,m}=0$ for any pair of $(n,m)$ except for $\alpha_{0,0}(=1)$, resulting in $S_{+}=S_{-}=0$ and $\rho^{s}_{\textrm{(a)}}=1/(2\pi)^{2}$ (disordered state (a)).

The solutions of the phase wave state (b) can be obtained exactly by the choice of $S_{+}=0$ or $S_{-}=0$. With $S_{-}=0$, we find the steady-state equations for $\alpha_{n,0}$ as

[TABLE]

which are decoupled from all other coefficients $\alpha_{n,m}$ for $m\neq 0$. With the use of the recurrence relation of the modified Bessel function, it is straightforward to find the solution of Eq. (27) as

[TABLE]

where $I_{n}$ is the modified Bessel function. The order parameter value $S_{+}$ can be determined self-consistently by the above equation for $n=1$ with $\alpha_{1,0}=S_{+}$. Note that $S_{+}$ depends only on $J_{+}/D$, independent of $J_{-}/D$. For small $S_{+}$, we obtain

[TABLE]

implying that the phase wave state is physically meaningful only for $J_{+}/D\geq 2$. All other coefficients are set to zero, which satisfy the steady-state condition, resulting in $\rho^{s}_{\textrm{(b)}}=\frac{1}{(2\pi)^{2}}e^{\frac{J_{+}S_{+}}{D}\cos(X-\Phi_{+})}/I_{0}(\frac{J_{+}S_{+}}{D})$. There exists also the symmetric solution with $S_{+}=0$.

The steady-state solutions for the mixed state (c) and the ordered state (d) cannot be derived analytically. Instead, we truncate higher-order Fourier modes in the dynamic equation (26) by keeping the modes of $(n,m)$ up to order $\ell(=|n|+|m|)$. Then, steady-state solutions are calculated by numerical iterations with $\dot{\Phi}_{\pm}=0$. The solutions well converge with increasing $\ell$ and saturate around $\ell=10$. This fast convergence is due to the exponential decay of $\alpha_{n,m}$ with $\ell$ in the steady state. These iteration results are shown by the black solid lines ($S_{+}$) and the purple ones ($S_{-}$) in Fig. 2 (a)-(d) for various values of $J$ and $K$ with $D_{x}=D_{\theta}=0.15$. Both the mixed state and the ordered state are well established and the exact solution of the phase wave state matches with numerical results perfectly well. We note that all coefficients $\alpha_{n,m}$ turn out to be real in all steady states, even if one starts from a complex initial state.

In the EQ situation ($J_{-}=0$ and $d=0$), the exact solution is possible due to decoupling of the dynamic equations, yielding $S_{+}=S_{-}\equiv S$, which has the same value of $S_{+}$ in the phase wave state. The steady-state PDF is given by a simple product of those of the phase wave states for $X$ and for $Y$, and the steady-state solutions are given by

[TABLE]

Finally, we remark that the vanishing Kuramoto order parameters $R_{x}$ and $R_{\theta}$ in Eqs. (6) and (7) are guaranteed in our Fourier mode analysis, as $R_{x}$ and $R_{\theta}$ correspond to half-integer Fourier modes such as $\alpha_{\frac{1}{2},\pm\frac{1}{2}}$, which do not exist as a Fourier mode with periodicity of $2\pi$.

VII linear stability analysis of the disordered state

The disordered (incoherent) state is characterized by $\alpha_{n,m}=0$ for any pair of $(n,m)\neq(0,0)$. We take a small perturbation to the disordered state such that $\alpha_{n,m}=\epsilon c_{n,m}$ for $(n,m)\neq(0,0)$ with a small parameter $\epsilon$. From the dynamic equation (26) with $S_{+}=\epsilon c_{1,0}$ and $S_{-}=\epsilon c_{0,1}$, we find, up to the linear order in $\epsilon$,

[TABLE]

Thus, the stability condition is simply given by $J_{+}/D<2$ as $D(n^{2}+m^{2})+2dnm=D_{x}(n+m)^{2}+D_{\theta}(n-m)^{2}>0$. As the other steady states are defined only for $J_{+}/D>2$, the disordered state should be globally stable for $J_{+}/D<2$, which agrees perfectly well with numerical data obtained by iterations of the dynamic equation (26).

VIII Numerical simulations

We perform numerical simulations by integrating Eqs. (11) and (12) for various values of $J$ and $K$ with temperatures $D_{x}=D_{\theta}=0.15$ for $N=10^{4}$. We used the Euler method with step size $\delta t=0.01$ for $M_{t}=10^{5}$ time steps, where the initial $M_{t}/2$ steps of each run were discarded as a transient, after which $S_{+}$ and $S_{-}$ were measured and averaged over remaining time for 10 independent samples.

Figure 2 shows simulation data (symbols) for $S_{\pm}$ versus $K$ when $J=0.2$, 0.4, 3.0 and when $J=K$. All data are in full agreement with exact solutions and iteration results in Sec. VI. Note that the transition from the disordered state to other steady states should occur at $J_{+}=2D$, derived as the instability threshold in the preceding section. With $D=0.3$, this transition point is located at $J+K=4D=1.2$, which are fully consistent with data. Small finite-size effects in simulation data can be seen near the transition from the disordered state.

IX Phase diagram

Figure 3 presents the phase diagram in the $(e^{-J_{+}},e^{-J_{-}})$ plane for convenience, when the temperatures are equal ($D_{x}=D_{\theta}$, i.e. $d=0$). In this case, the dynamic equations (1) and (2) are symmetric under the interchange of $x\leftrightarrow\theta$ and $K\leftrightarrow J$. Thus, the phase diagram for $e^{-J_{-}}>1$ can be easily deduced from that for $e^{-J_{-}}<1$. In order to locate the phase boundary more precisely, we use the exact solutions and the iteration results with terms up to the order $\ell=|n|+|m|=20$.

The phase diagram reveals several noteworthy features. First, the mixed phase $(S_{1},S_{2})$ colored in red clearly exists and separates the phase wave phase $(S,0)$ and the ordered phase $(S,S)$ completely, though the mixed state region itself is quite small in size. In Fig. 4, the order parameter values are plotted along a vertical line in the phase diagram at $e^{-J_{+}}=0.3$ in the range of $0.40\leq e^{-J_{-}}\leq 0.46$. The value of $S$ remains constant in the phase wave phase $(S,0)$, as predicted by the exact solution in Sec. VI. This value is also observed at the EQ point ($e^{-J_{-}}=1$) in the ordered phase. As $e^{-J_{-}}$ is increased, the system goes through the mixed and the ordered phase in succession before arriving at the EQ point.

Second, the disordered phase $(0,0)$ appears for $J_{+}/D<2$ for all values of $J_{-}$, consistent with the linear stability analysis result presented in Sec. VII. Third, all four phases converge at a single multicritical point where $J_{+}/D=2$ and $J_{-}/D=1$. The value of $J_{-}$ at the multicritical point will be derived exactly through a perturbation theory near the disordered state in the following section.

Finally, we point out that direct transitions from the disordered phase to the ordered phase are possible not only along the EQ line but also at $J_{+}=2D$ (blue line in Fig. 3) in nonequilibrium situations. This feature represents a major difference from the system with quenched disorder only ref:1Dwithfrequency , where the direct transition is possible only at the symmetric point $J=K$ ($e^{-J_{-}}=1$). It would be of interest to study the existence and location of a multicritical point when both quenched and thermal disorder are present simultaneously, which is left for future study.

We also studied the case with two different temperatures ($D_{x}\neq D_{\theta}$). In this case, the dynamic equations are symmetric with $D_{x}\leftrightarrow D_{\theta}$, in addition to $x\leftrightarrow\theta$ and $K\leftrightarrow J$. The qualitative features of the phase diagram are similar to those depicted in Fig. 3. We find the multicritical point, the $J_{-}$ value of which is given by $J_{-}=D+d$, while the $J_{+}$ value remains unchanged at $2D$ (see Sec. X).

X perturbation theory near the disordered state

For convenience, we rewrite the mode-coupled equations (26) with $\dot{\Phi}_{\pm}=0$ as

[TABLE]

with new parameters

[TABLE]

where time $t$ is scaled by a factor of $2/J_{+}$ for simplicity. We consider $u,v\geq 0$ only, because only the trivial disordered phase is possible for $u<0$ and there is a reflection symmetry under $v\leftrightarrow-v$. Also the positivity of diffusion constants ($D_{x}>0,D_{\theta}>0$) demands $|z|<u$.

As the ordering begins to emerge for $u\lesssim 1$, we introduce a small positive parameter $\delta$ such that $u=1-\frac{\delta}{2}$. From the exact relation of Eq. (28) for the phase wave state, we obtain the steady-state values for the Fourier coefficients $\alpha_{n,0}$ in lower orders of $\delta$ as

[TABLE]

where we utilized the expansion of $I_{n}(a)$ for small $a$ and $S_{+}^{2}\approx\delta(1-\delta/3+\delta^{2}/72)+O(\delta^{3})$. Along the EQ line, it is also straightforward to obtain from Eq. (30)

[TABLE]

which demonstrates the rapid exponential decay of the steady-state Fourier coefficients with $\ell=|n|+|m|$. These analytic results motivate us to assume that $\alpha_{n,m}$ in lower orders for general steady states can be expressed as

[TABLE]

with $a_{0,0}=1$ and $b_{0,0}=0$ for normalization and $a_{n,m}=a_{-n,-m}$ and $b_{n,m}=b_{-n,-m}$.

X.1 steady states

Insert this leading order into Eq. (32) with ${\dot{\alpha}}_{n,m}=0$ for steady states, we can then determine $a_{n,m}$ order by order in $\delta$. Up to the order of $\delta^{3/2}$, we obtain relations between $a_{n,m}$ by

[TABLE]

There are also relations for $b_{1,0}$ and $b_{0,1}$, which are shown in the Appendix A.

By combining the first four equations, four distinct solutions can be obtained, each of which represents a fixed point such as $(a_{1,0},a_{0,1})=$ $(0,0)$, $(1,0)$, $(0,1)$, $(\sqrt{a},\sqrt{a})$ with $a=[1+4v(v-z)/(1-z^{2})]^{-1}$. The $(0,0)$ fixed point is associated with the disordered phase (a), while the $(1,0)$ and $(0,1)$ points denote two equally probable states of the phase wave phase (b). The $(\sqrt{a},\sqrt{a})$ point corresponds to the symmetric ordered phase (d). All other coefficients $a_{n,m}$ up to $\ell=3$ can be calculated by inserting the fixed point solutions into the above equations. Up to this order, the mixed phase (c) does not appear.

We also report the fixed point values of $(b_{1,0},b_{0,1})$, which represent the $O(\delta^{3/2})$ term of $(\alpha_{1,0},\alpha_{0,1})$ as $(0,0)$, $(-1/6,0)$, $(0,-1/6)$, and $(b,b)$ for four fixed points, respectively. The value of $b_{n,m}$ is discussed in the Appendix A. All other terms are higher-order than $O(\delta^{3/2})$, thus can be set to zero.

X.2 linear stability analysis of the phase wave state

The steady state solution of the phase wave state is quite simple as $a_{1,0}=1,a_{2,0}=1/2,a_{3,0}=1/6,b_{1,0}=-1/6$ with all other coefficients $a_{n,m}=0$, up to the order of $\delta^{3/2}$. In order to check its stability, we take a small perturbation to this state by adding $\epsilon c_{n,m}$ to the steady state values of $\alpha_{n,m}$ except for $(0,0)$. Then, it is straightforward to find up to the linear order in $\epsilon$ from Eq. (32) as

[TABLE]

in the lower order of $\delta$.

Note that most of $c_{n,m}$ vanish in the long time limit for small $\delta$, as $n^{2}+m^{2}+2znm>0$. In contrast, $c_{0,1}$ may diverge depending on details of Eq. (47), thus $c_{1,1}$ and $c_{1,-1}$ may also diverge as their dynamics involve $c_{0,1}$, as in Eqs. (50) and (51). Consequently, we can drop vanishing terms as $c_{3,0}$, $c_{1,-2}$, $c_{1,2}$, $c_{2,1}$, and $c_{2,-1}$ in Eqs. (46)-(51), leading to the closed equations for $c_{1,0}$, $c_{0,1}$, $c_{2,0}$, $c_{0,2}$, $c_{1,1}$, and $c_{1,-1}$. These linear equations can be analyzed by the standard eigenvalue analysis, yielding five negative eigenvalues and one eigenvalue $\lambda=[(2v-z)^{2}-1]\delta$. Hence, the stability condition is given by

[TABLE]

implying $|J_{-}-d|<\frac{1}{2}|J_{+}|$ in terms of original parameters, which predicts the location of the multicritical point precisely as $J_{-}=D+d$ with $J_{+}=2D$ (see Fig. 3 and Sec. IX).

XI Examples of real word swarmalators

In this section, we list some examples of real world swarmalators which exhibit some of the collective states we have found (see ref:review1 for an exhaustive list of real world swarmalators). The async(disordered state), sync(ordered state), and phase wave states have all been observed in these systems. The mixed state has not been explicitly observed, but due to its close resemblance to the phase wave, it may in fact have been realized but not observed.

(i) Vinegar eels are a type of nematode (a family of worm microswimmers) commonly found in beer mats and tree slime ref:vinegareels . When suspensions of these eels are prepared on glass disks they swarm around the 1D edges and synchronize the beating of their tails, thereby satisfying the defintion of a 1D swaramlators. Under certain conditions they form metachronal waves, which are equivalent to the phase waves states here found but winding number $k>1$ ref:vinegareels .

(ii) Sperm is another type of microswimmer which can be considered a 1D swarmalator. When sperm from ram semen is confined to 1D circular geometeries, they show a transition from an isotropic state (equivalent to our async state) to an uniformly rotating vortex (equivalent to our phase wave state) ref:sperms . Sperms are free to move in 2D and also form synchronous clusters.

(iii) C. elegans are another type of microswimmers which also swarm and sync the gait of their tails. When confined to 1D channels they form synchronous clusters, which is analogous to the sync state ref:yuan2014 .

(iv) Magnetic domain walls are characterized by a center of mass $x$ and a magnetic dipole vector with orientation $\theta$. When subject to external forcing, both $x$ and $\theta$ undergo periodic motion as so the walls can be considered swarmalators ref:magneticdomainwalls . In an experiment with $N=2$ walls, different varieties of sync phenomena can be observed.

(v) Bristle bots are homemade ‘automata’ made from a toothbrushes and powered by a simple cell-phone motor ref:bristle . When placed on a circular drum, they self-organize around the ring-like edge. The incoherent state is in  Fig. 9 a, phase wave state (Fig. 9 b), and the sync state (Fig. 9 c) are observed (By Fig. 9 we mean Fig. 9 in Ref. ref:bristle ).

XII Discussions

The joint action of swarming and synchronization defines a new type of emergence about which little is known. The paper is part of effort to explore this unchartered terrain by studying a simple 1D model of swarmalators subject to thermal noise. We found four collective states, encapsulated their stabilties in a phase diagram, illustrated their transitions with order parameter curves $S_{\pm}(J,K)$, and discussed realizations in nature. These states were previously reported in a study of the 1D swarmlator model with quenched disorder ref:1Dwithfrequency , but the bifurcation structure for the thermal noise is different: As shown in the phase diagram in Sec. IX, through the blue solid line, the phase transition directly occurs from the $(0,0)$ state to the $(S,S)$ state without going through the $(S,0)$ or $(S_{1},S_{2})$, which is very different from the 1D model with quenched disorder ref:1Dwithfrequency . In other words, in the absense of thermal noise, there is only one “tetracritical” point there, which means that the direct phase transition from the $(0,0)$ to the $(S,S)$ state occurs only at that point. However, when the thermal noise comes in the system, such region for the direct phase transition (same as that for the EQ line) is found to enlarge, not only at the multicrossing point, which is induced by the thermal noise.

We here considered identical swarwalators with same natural frequency, and showed that the four states can be induced by means of only the thermal noise and nonequilibrium interactions. In other words, the nonequilibrium features such as the four states stem from the thermal noise and the interaction between the units in the system. We here claim that the quenched disorder is not the inevitable component to induce the four states. Specifically, we have unequivocally detected the presence of the mixed state in the system. We have conducted a stability analysis of the phase wave state using the perturbation theory, and determined the multicritical point where the four states converge. Our findings demonstrate a great level of agreement with numerical results, indicating the originality of our work.

The main theoretical takeaway of our work that allowing oscillators to swarm greatly enriches their collective behavior. The Kuramoto model of regular oscillators with thermal noise, for example, has a single transition: from incoherence to synchrony ref:footnote . The 1D swarmalator model on the other hand has four collective states and four distinct transitions, as depicted on Fig. 2. An interesting implication of this novel bifurcation structure is that synchrony does not increase monontonically with the phase coupling $K$ (as happens in the Kuramoto model). Holding $J$ constant and turning up $K$ can either leave you stuck in the phase wave (Fig 2(a)) or take you out of the sync state into the mixed state (Fig 2(b)) – in both cases the overall phase synchrony decays.

An tantalizing goal for future work would be to find the mixed state in a realworld system of swarmalators. While to our knowledge it has not directly reported, we suspect the mixed state might be lurking within the recently reported metachronal waves of vinegar eels ref:vinegareels . As we mentioned, metachronal waves are mimicked by our phase wave; since the phase wave is visually similar to the mixed state, the phase wave might have been misidentified. If experimenters could extract the $(x_{i},\theta_{i})$ of each eels, then the $S_{\pm}$ could be measured and in principle the mixed state could be detected. The search for the mixed state could also include other biological microwswimmers, magentic domain walls, or bristlebots. There are also theoretical avenues to explore in the future. One could study the effects of delayed coupling, external forcing, or mixed sign interactions.

XIII Acknowledgements

This research was supported by NRF Grants No. 2021R1A2B5B01001951 (H.H) and No. 2017R1D1A1B06035497 (H.P.), and individual KIAS Grants No. PG064901 (J.S.L.), and No. QP013601 (H.P.) at Korea Institute for Advanced Study.

Appendix A relations for higher-order coefficients

Relations for higher-order coefficients are shown as

[TABLE]

For the fixed point $(a_{1,0},a_{0,1})=(1,0)$, we obtain $a_{3,0}=1/6$ and $a_{0,3}=0$ from Eq. (41) and all other coefficients are zero. With these values, it is easy to derive that $b_{1,0}=-1/6$ and all other $b_{n,m}=0$. Note that this value is consistent with Eq. (34).

In the ordered phase with $(a_{1,0},a_{0,1})=(\sqrt{a},\sqrt{a})$, the above equations for $b_{n,m}$ possess a symmetry under the exchange of indices ($n\leftrightarrow m$) and also linearity, implying that $b_{1,0}=b_{0,1}$, $b_{2,0}=b_{0,2}$, and $b_{1,1}(v,z)=b_{1,-1}(-v,-z)$. This observation excludes the possibility of the mixed phase ($S_{+}\neq S_{-}$), at least up to this order. It is straightforward to get explicit expressions for $b_{n,m}$ for the symmetric ordered phase, but they are too complicated to be displayed here.