Swarmalators with thermal noise

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.

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 ∈ R d and phase θ ∈ 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 [28] or the Kuramoto model for synchronization studies [29], which captures the essential ingredient while remaining simple enough to solve.This paper focuses on a simple one-dimensional (1D) swarmalator model [30,31] which could fill this 'modeling gap' in studies of swarmalation, as Verberk [32] 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 [31], 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 [31].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 [33], 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 [34][35][36][37].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 where x i and θ i denote the position and phase of the ith swarmalator, respectively, and are both periodic with a period of 2π.One may add an identical natural frequency ω for phase in Eq. ( 2), but it can be eliminated by taking a simple transformation of θ i → θ i + ωt.Likewise, the identical directed velocity v can be also set to zero without loss of generality.ξ x i and ξ θ i represent thermal noise with zero mean, and characterized by where D x and D θ denote the "temperature" of the heat bath for x and θ variables, respectively.The two noises are independent of each other, i.e. ξ x i ξ θ j = 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 [30] and for random nonidentical ones [31].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 , θ i ) ∈ T 1 .In this sense, it can be referred to as a generalized version of the thermal Kuramoto model [35] 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 [33], wherein spins possess two degree of freedoms, x and θ, that are coupled.For the special case of J = K and D x = D θ , 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 where an EQ (magnetic) phase transition is expected as J varies.For J = 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 = D θ , 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 [38] where θ 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 θ, (c)mixed between (b) and (d), and (d) two well-defined isotropic clusters separated by a distance of π in both directions.
One may consider the Kuramoto synchronization (or XY magnetic) order parameters [29] for each variable x and θ as where R x and R θ measure the magnitude of synchronization (real and non-negative) in x and θ, respectively, and Θ x and Θ θ are their mean angles (real).• • • 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 "π" transform of x j → x j + π and θ j → θ j + π for any j, as noted in Refs.[30] and [31].It implies that the dynamics of any (x, θ) state is equivalent to that of the (x+π, θ +π) 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 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 [30] and [31], as where S ± is real and non-negative by definition and Φ ± is real.This particular form of the correlation measure can be inferred from the linearity of correlations between x and θ with slopes of ±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 = 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 which imply that the new variables X and Y are both periodic with a period of 4π.Subsequently, the dynamic equations in Eqs. ( 1) and ( 2) can be simplified to with J ± = (J ± K)/2, and S ± and Φ ± are defined in Eqs. ( 8) and ( 9).The new noises, ξ X i (t) and ξ Y i (t), are defined as that are not independent of each other in general.Simple calculations lead to ξ where The correlation between the two noises D XY disappears when the two temperatures are equal (D x = D θ ), It is worth noting that the dynamic equations in Eqs.(11) and (12) are invariant under the transformation of X i → X i + 2π or Y i → Y i + 2π for any i, implying that the variables X and Y can be treated as periodic with a period of 2π, instead of 4π, if the initial distributions possess the same periodicity.
When J = K (J − = 0) with D x = D θ (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 = K) or a thermal gradient caused by multiple reservoirs (D x = D θ ).

V. FOURIER MODE ANALYSIS
In the N → ∞ limit, the system state can be described by a continuous probability density function (PDF) ρ(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π, thus their range is both set to (0, 2π) with the normalization condition as The Fokker-Planck equation [39], derived from the dynamic equations of ( 11) and (12) with Eqs. ( 17) and (18), is given by where the probability currents J X and J Y are with The Fourier series of the PDF ρ(X, Y, t) in X and Y with period 2π is given by where the extra factors Φ + and Φ − are inserted for convenience.As the PDF should be real, α * n,m = α −n,−m holds for all (n, m) and the normalization condition yields α 0,0 = 1.
The order parameters in Eqs. ( 8) and ( 9) can be expressed in terms of Fourier coefficients as It should be noted that, as S ± is real, α 1,0 and α 0,1 must also be real.Substitution of Eq. ( 23) into the Fokker-Planck equation of Eq. ( 19) yields time-dependent modecoupled equations as where S + = α 1,0 and S − = α 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.
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 α n,0 as which are decoupled from all other coefficients α n,m for m = 0.With the use of the recurrence relation of the modified Bessel function, it is straightforward to find the solution of Eq. ( 27) as 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 α 1,0 = S + .Note that S + depends only on J + /D, independent of J − /D.For small S + , we obtain implying that the phase wave state is physically meaningful only for J + /D ≥ 2. All other coefficients are set to zero, which satisfy the steady-state condition, resulting in ρ s 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 (= |n|+|m|).Then, steady-state solutions are calculated by numerical iterations with Φ± = 0.The solutions well converge with increasing and saturate around = 10.This fast convergence is due to the exponential decay of α n,m with 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 θ = 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 α 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 − ≡ 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 Finally, we remark that the vanishing Kuramoto order parameters R x and R θ in Eqs. ( 6) and ( 7) are guaranteed in our Fourier mode analysis, as R x and R θ correspond to half-integer Fourier modes such as α 1 2 ,± 1 2 , which do not exist as a Fourier mode with periodicity of 2π.

VII. LINEAR STABILITY ANALYSIS OF THE DISORDERED STATE
The disordered (incoherent) state is characterized by α n,m = 0 for any pair of (n, m) = (0, 0).We take a small perturbation to the disordered state such that α n,m = c n,m for (n, m) = (0, 0) with a small parameter .From the dynamic equation ( 26) with S + = c 1,0 and S − = c 0,1 , we find, up to the linear order in , ċn,m = − D(n 2 + m 2 ) + 2dnm c n,m for other (n, m).
Thus, the stability condition is simply given by 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 θ = 0.15 for N = 10 4 .We used the Euler method with step size δ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 ± 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.tures are equal (D x = D θ , i.e. d = 0).In this case, the dynamic equations ( 1) and ( 2) are symmetric under the interchange of x ↔ θ and K ↔ 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 = |n|+|m| = 20.

IX. PHASE DIAGRAM
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 ≤ e −J− ≤ 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 [31], 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 = D θ ).In this case, the dynamic equations are symmetric with D x ↔ D θ , in addition to x ↔ θ and K ↔ 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 Φ± = 0 as with new parameters where time t is scaled by a factor of 2/J + for simplicity.We consider u, v ≥ 0 only, because only the trivial disordered phase is possible for u < 0 and there is a reflection symmetry under v ↔ −v.Also the positivity of diffusion constants (D x > 0, D θ > 0) demands |z| < u.
As the ordering begins to emerge for u 1, we introduce a small positive parameter δ such that u = 1 − δ 2 .From the exact relation of Eq. ( 28) for the phase wave state, we obtain the steady-state values for the Fourier coefficients α n,0 in lower orders of δ as where we utilized the expansion of I n (a) for small a and ).Along the EQ line, it is also straightforward to obtain from Eq. ( 30) which demonstrates the rapid exponential decay of the steady-state Fourier coefficients with = |n|+|m|.These analytic results motivate us to assume that α n,m in lower orders for general steady states can be expressed as 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 .

A. steady states
Insert this leading order into Eq.( 32) with αn,m = 0 for steady states, we can then determine a n,m order by order in δ.Up to the order of δ 3/2 , we obtain relations between a n,m by 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), ( √ a, √ 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 ( √ a, √ a) point corresponds to the symmetric ordered phase (d).All other coefficients a n,m up to = 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.

B. 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 δ 3/2 .In order to check its stability, we take a small perturbation to this state by adding c n,m to the steady state values of α n,m except for (0, 0).Then, it is straightforward to find up to the linear order in from Eq. ( 32) as in the lower order of δ.

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 [26] 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 [3].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 [3].
(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) [2].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 [40].
(iv) Magnetic domain walls are characterized by a center of mass x and a magnetic dipole vector with orientation θ.When subject to external forcing, both x and θ undergo periodic motion as so the walls can be considered swarmalators [41].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 [42].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. [42]).

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 ± (J, K), and discussed realizations in nature.These states were previously reported in a study of the 1D swarmlator model with quenched disorder [31], 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 [31].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 [43].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 ( 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 [3].As we mentioned, metachronal waves are mimicked by our phase wave; since the phase wave is visually simi-lar to the mixed state, the phase wave might have been misidentified.If experimenters could extract the (x i , θ i ) of each eels, then the S ± 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.

Figure 3 FIG. 2 .
Figure3presents the phase diagram in the (e −J+ , e −J− ) plane for convenience, when the tempera-

FIG. 4 .
FIG. 4. (Color online) Plots of S± versus e −J − for e −J + = 0.3.Symbols represent the data obtained from iteration results with = 20, and the lines are guide to the eyes.
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.