Oscillatory behavior in a system of swarmalators with a short-range repulsive interaction

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is an Accepted Manuscript of an article published in Physical Review E on June 2020, available at: https://doi.org/10.1103/PhysRevE.101.062202 .” Oscillatory Behavior in a System of Swarmalators with a Short1 Range Repulsive Interaction2 F. Jim´enez-Morales∗ 3 Dept. de F´ısica de la Materia Condensada. Universidad de Sevilla4 P.O. Box 1065, 41080-Sevilla, Spain5 (Dated: February 1, 2020)6 Abstract7 One of the computational models proposed to study emerging phenomena in spatially distributed8 systems is the swarmalator model [Kevin P.O’Keeffe, Hyunsuk Hong and Steven H. Strogatz,9 Nat. Commun. 8, 1504 (2017) ] where only long-range interactions are considered. But in living10 systems many of the collective behaviors observed arise from the combined effect of long and11 short range interactions. In this work an extension to the swarmalator model is presented which12 includes a Gaussian short-range repulsive term of the type 1 σe−|x|2/σ and present the results of13 numerical simulations. Using several order parameters we can distinguish between static and14 dynamic aggregations and between synchronous and asynchronous states. We have found six long-15 term collective states some of them not previously reported being the most remarkable one an16 oscillatory collective behavior. The results obtained show a multiplicity of complex behaviors that17 extend the applicability of the swarmalator model.18 PACS: 05.45-a, 89.75-k19 ∗[email protected] 1 I. INTRODUCTION20 Nature is full of systems composed of a multitude of identical and spatially extended21 elements that show many and diverse collective behaviors. Some examples frequently cited22 in the literature are the collective migration of bacterial colonies [2], the complex structures23 developed by social insects [19] or flocks of birds that move together [7]. In a broad sense, all24 collective behaviors are characterized by the existence of a spatial self-organization (swarm-25 ing) and/or a temporal coordination of the individual units (synchronization) [21].26 27 One of the most interesting computational models proposed to study swarming and28 synchronization phenomena in spatially extended systems is the concept of “swarmalator”29 introduced in [16] and later expanding in [8, 14, 15]. An experimental application of the30 model has been achieved in [1] where it was implemented in a group of mobile robots and31 in [13] where it has been proposed for the control of decentralized systems such as mobile32 multi-agent robotic platforms.33 34 Swarmalators are small mobile units in the plane or in space which have an internal35 variable –the phase– that affect how the units move. Early contribution to those kind36 of models were proposed in [9, 10, 18], and in [12] where the same concept was used to37 explain synchronized rotations in aggregates of magnetotactic bacteria. Here we focus on38 the swarmalator model [16] which stands out for its flexibility and simplicity, their dynamical39 equations are the following:40 ˙ xi=1 N N X j6=i (xj−xi) [Iatt(1 + Jcos(θj−θi)) −Irep] (1) 41 ˙ θi=K N N X j6=i sin(θj−θi)Hatt (2) Nis the number of units with spatial coordinates xiand θiare the corresponding phase.42 The function Iatt represents the spatial attraction and Irep the repulsion. The influence of43 the phase interaction is considered of the Kuramoto type. The system can be studied in44 terms of only two parameters: Kwhich is the phase coupling strength and Jwhich measures45 the extent to which phase similarities enhances spatial attraction. In the original model the46 attractive and the repulsive terms were selected as power laws:47 2 Iatt =Hatt =1 |xj−xi|, Irep =1 |xj−xi|2(3) Depending on the values of the pair (J, K) the long term state reached by the swar-48 malators system can be classified according to two main characteristics: mobility and phase49 coherence of the units. In this way the reported final states were: static-synchronous,50 static-asynchronous, static phase wave, a splintered state in which the units make small51 oscillations in each cluster and a dynamic state in which the swarmalators are divided into52 counter-rotating groups [16].53 54 Potential long-range laws like that of equations 3 are frequent in the study of physical55 systems. But in biological systems, although there are long-range interactions there are56 also short-range interactions. In fact it is these short-range interactions that make the most57 biological sense since most of the functions between cells or animals occur by direct contact58 between them [3, 20]. Therefore, in order to understand many of the collective behaviors59 observed in living beings, both types of interaction must be taken into account. This is the60 case, for example, in the adhesion of bacterial colonies to a surface caused by a combination61 of long and short range forces [6]; in the flocks of birds [4] where the local interaction be-62 tween individuals decays exponentially with the distance between them; or in insect swarms63 of flying midges [17] that shows a long-range attractive force and a short-range repulsion.64 65 In this paper we study a swarmalator system that includes both types of interaction, one66 attractive long-range, as in the original model, and another short-range repulsive force. The67 question to be solved is whether the system of units is grouped or dispersed as a result of68 their phase difference and in each case the different spatio-temporal patterns that arise.69 3 II. THE MODEL70 The dynamical equations of the model are:71 ˙ xi=1 N−1 N X j6=i (xj−xi)"(1 + Jcos(θj−θi)) |xj−xi|−1 σe−|xj−xi|2/σ#(4) ˙ θi=K N−1 N X j6=i sin(θj−θi) |xj−xi|(5) The attractive force term is the same as in the original model but the repulsive term is72 selected as a Gaussian function Irep =1 σe−|xj−xi|2/σ where σis a measure of the range of the73 local interaction. Also for simplicity we take J= 1 and the study of the possible collective74 behaviors is made in terms of the pair of parameters σand K.75 76 Considering the set of swarmalators as if it were a system of particles we are going to77 describe its mobility by means of the modulus of the linear momentum Pand the kinetic78 energy Tdefined as follows:79 P=1 N N X i=1 ˙ xi , T =1 2N N X i=1 |˙ xi|2(6) The most direct way to detect a collective movement is by calculating the velocity of80 the center of mass of the system or the lineal momentum P. Usually this quantity is used81 normalized but no here because we want to compare its value as we change the parameters82 σand K. Non-zero values of Pindicate a net displacement of the system. In the original83 model the long-range attractive and repulsive terms make Pnull. But in our case due to84 the short range term it is possible to have net displacements of the system. Also and in85 order to distinguish between a dynamic and a static situation, it is necessary to know the86 value of T, for example when the assembly is rotating around a fixed center of mass the87 value of Pis null but not T.88 89 To measure the coherence of the phases it is customary through the module Rof the90 complex order parameter [11] defined as:91 Z=ReiΦ=1 N N X j eiθj(7) 4 where Φ is the average angle. Also for the swarmalator model it has been used another order92 parameter Swhich measures the correlation between the phases and the spatial distribution.93 The definition of S[16] is given by:94 W±=S±eiΨ±=1 N N X j=1 ei(φj±θj)(8) where φi=tan−1(yi/xi) is the spatial angle. The parameter Sis the maximum value95 between S+and S−.96 III. SIMULATIONS97 We take a system of N= 500 units in two dimensions. The initial state is selected98 randomly with the phases in the range [−π, +π], and the spatial coordinates x and y in the99 range [−2,2], which are the same range of values used in [8, 14–16]. Kvalues are in the100 range [−1,1] and the values of the σparameter range from 0.1 to 1.2. The set of differential101 equations 4, 5 is solved by the standard Runge-Kutta method with dt = 0.1 and we study102 the the long term behavior reached for the system of swarmalators for each pair of values103 (K, σ) after 1000 time steps.104 105 Figure 1 summarizes our results. Each point is an average over 5 different initial situa-106 tions. The different types of behaviors are classified considering whether the configuration107 is static or dynamic; and if the phase of the elements is synchronized or not. In this way the108 final state can be: static and synchronized (SS), dynamic and synchronized (DS) or dynamic109 asynchronous (DA). Other cases are observed when the phases are distributed among the110 elements as a wave, in that case we can have a static wave (SW) or a dynamic wave state111 (DW). The last type of behavior found is when the system is formed by several disconnected112 zones that can be also dynamic or static, we call this splintered state (Sp).113 114 The results of Figure 1 show the state the set arrives at the final step of the simulation.115 A complementary information is obtained measuring the temporal average of the different116 parameters R, S, P and T. For doing that we let a transient time of 500 time steps and the117 average is obtained accumulating the value in the remaining τ= 500 time steps. So for118 example the temporal average of Ris given by R=1 τPtR(t) where R(t) is the instanta-119 5 0 0.5 1 1.5 2 -1 -0.5 0 0.5 1 σ K SS SW DW DS DA Sp FIG. 1. The behaviors shown by the model in terms of the Kparameter and the radius σ. Each point is an average over five different initial conditions. The number of swarmalators is N= 500, dt = 0.1 and the final state is obtained after 1000 time steps. We have found six kinds of behaviors: (SS)static synchronized, (SW) static wave, (DW) dynamic wave, (DS) dynamic synchronized, (DA) dynamic asynchronous and (Sp) splintered state. neous value. Similar expressions are for the average values of the other parameters. Figure120 2a shows the contour map of the average parameter R. It is observed that regardless of121 the radius of the interaction σthe coherent states are found for positive values of coupling122 K > 0 where R > 0 being the degree of coherence greater the closer is Rto the value 1.123 For K < 0 and also regardless of σthe parameter R= 0, so that in that region the set does124 not reach a synchronized state. However, that does not mean that there is an absence of125 symmetry. Figure 2b shows a contour map of the time average of the Sparameter that126 gives us a measure of the existence of correlations between the spatial angle and the phase127 of each element. In the region K > 0, where the coherent state are found, the values of128 Sare null and there is no correlation between the spatial variables and the phase. But129 now it is in the zone of K < 0 where S > 0 indicating the existence of states that al-130 though not all of its elements have the same phase there are strong correlations between the131 phase and the spatial positions. It is in this region where the states called waves are located.132 133 The Rand Sparameters are useful for discriminating between synchronized, asyn-134 chronous and wave states. To study whether the final states are dynamic or static we use135 6 a) b) 0.1 -1 -0.5 0 0.5 1 K 0.5 1 1.5 2 σ 0 0.5 1 0.1 -1 -0.5 0 0.5 1 K 0.5 1 1.5 2 σ 0 0.5 1 FIG. 2. Contour plot of the time averages of the Kuramoto order parameter Rand the parameter S versus the interaction radius σand K. The values, in arbitrary units, are an average over τ= 500 time steps after a transient time of another 500 time steps. a) The plot of Rshows that coherent states are only obtained for positive values of Kwhere Ris higher than 0. For negative values of Kcoherent states are not found as R= 0. b) The plot of S. Correlations between the phase and the spatial angle is found only for K < 0 whereas for K > 0 there is not correlation between the phase and S= 0. the temporal averages of the modulus of the linear momentum Pand the kinetic energy136 T. Figure 3 shows a contour map of P, and Tin arbitrary units as a function of Kand137 σ. The zone colored in black corresponds to states where the set remains static while the138 one colored in red and yellow indicates those parameters for which the set is moving. For139 high values of Pthe whole system moves with the same speed and the center of mass of140 the system defined by rCM =1 NPixiis not fixed. This happens for example when K= 0141 regardless of σ. Other swarmalator model in which there are collective displacements has142 been found by Hong [8] where a short-range term is introduced into the equation of evolution143 of the phases in addition to a thermal noise. As stated before, in our case these collective144 behaviors are due to the short-range of the repulsive term in the spatial variable and arise145 in a more direct and simpler way.146 147 Apart from the dynamic states in which a net displacement occurs, there are other148 dynamic states in which the center of mass is fixed while the individual swarmalators are149 rotating around it and also there are other cases in which the units are performing small150 movements in several clusters. In Figure 3 that region is located in the lower left margin of151 7 a b 0.1 -1 -0.5 0 0.5 1 K 0.5 1 1.5 2 σ 0 0.5 1 0.1 -1 -0.5 0 0.5 1 K 0.5 1 1.5 2 σ 0 0.5 1 FIG. 3. Contour plots of the time average of (a:)Pand (b:)Tversus σand K. The number of units is N= 500 and dt = 0.1. The values, in arbitrary units, are an average over 500 time steps after a transient time of another 500 time steps. Dynamic states are mainly located for values of K= 0 and small values of the interaction radius σ. the plot in the range −1≤K≤ −0.6 and 0.1≤σ≤0.6.152 153 Now we are going to show the different spatial and angular diagrams that correspond to154 the states established in the Figure 1.155 156 The static synchronized state (SS) is obtained only for positive values of the coupling157 parameter K > 0 and consists of various concentric circles as it is shown in Figure 4-a. The158 colors corresponds to the phase angle θwhich is the same for all the cells so that the plot159 of the phase θversus the spatial angle φis a straight horizontal line ( Figure 4-b ). It is160 observed that as the σparameter increases the most external radius grows. In [15, 16] the161 static synchronous state was a circle whose surface is completely occupied by the elements162 all of them with the same phase but without the existence of empty spaces. For the range of163 parameters studied here that dense and circular distribution has not been found but instead164 this unreported and new configuration.165 166 Figure 5 shows the scatter plot in the x-y plane with the velocities of each element and167 the plot in the (φ, θ) space corresponding to the dynamic synchronized state (DS). This168 state is also found only for K > 0. The DS state consists of groups of swarmalators with169 a circular geometry with all the elements in phase, at least most of them, that can move170 8 a) b) 0 π 2π 0π2π θ φ 0 π 2π 0π2π θ φ FIG. 9. Scatter plots in (φ, θ) space corresponding to Figure 8-a,f. For the static wave state the relationship between the phase θand the spatial angle φis a straight line. Similar figures are obtained for the other cases of Figure 8. the X-Y space and the scatter plot in the φ, θ plane. In this case the aggregate consists of two222 ellipses with a wave-like distribution of the phases. Both ellipses move synchronously so that223 the major and the minor axis, denoted by Dxand Dy, lengthen and contract harmoniously.224 The time series of the mass center of the system, Rx(t) = 1 NPixi(t), Ry(t) = 1 NPiyi(t), is225 shown in Figure 11c. After a long transient time Rxand Ryoscillate both with the same low226 frequency f0≈0.003 s−1, and also they have a high oscillation frequency f1≈0.062 s−1. The227 trajectory followed by Rxand Ryis shown in Figure 11d here the color palette indicates the228 evolution time after transient time has been discarded. Initially the center of mass rotates229 forming small spirals until the final state is established with a larger circle.230 The oscillation of the set is also shown by measuring the evolution of the size of axes of231 the ellipse Dxand Dywhich are shown in Figure 12a. The power spectrum of Dxshows the232 high oscillation frequency f1(Dyalso shows the same frequency ). The oscillatory behavior233 is also observed in the components of the velocity of the center of mass Vx(t) = 1 NPi˙xi(t)234 and Vy(t) = 1 NPi˙yi(t) plotted in Figure 12c. The components of the velocity have also the235 same oscillation frequency f1.236 237 If the swarmalator model were useful to simulate biological systems it is interesting to238 point out that a similar movement is obtained in dense bacterial suspensions [5] where the239 the global oscillation is achieved despite the chaotic movement of each individual bacterium.240 Figure 13 shows the position and the velocity of a single swarmalator selected randomly.241 The plot of the component Y versus X shows various elliptical rotation around a central242 point. The components of the velocity shows a star like picture as can be seen in Figure 13.243 15 0π2π a) -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Y X -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Y X b) -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Y X -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Y X c) -0.2 0 0.2 -0.2 0 0.2 Y X -0.2 0 0.2 -0.2 0 0.2 Y X d) -0.2 0 0.2 -0.2 0 0.2 Y X -0.2 0 0.2 -0.2 0 0.2 Y X FIG. 10. Dynamic wave state for the pair of parameters (K, σ): a) (−0.90,0.2), b) (−0.90,0.4), c)(−0.60,0.8), d)(−0.6,1.2). The velocity of each particle has been plot with an arrow with the same color as the phase. To better visualize the velocities a scale has been used for a) and b), and a different one for c) and d). In cases a) and b) all the elements move with almost the same velocity. In cases c) and d) the inner strip is rotating with a velocity higher than the velocities of the outer elements. 16 a) b) -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Y X 0 π 2π 0 π 2π 0π2π θ φ c) d) 0.04 0.06 0.08 0.1 0.12 0.14 500 1000 1500 2000 2500 Rx, Ry t 0.095 0.1 0.105 0.11 0.115 0.12 0.05 0.06 0.07 t Ry Rx 0*100 1*103 2*103 3*103 FIG. 11. Swarmalator with a collective oscillation for the parameters K=−0.1, σ = 0.4. a) The plot in the plane X, Y .b) The scatter plot in (φ, θ). c) Time series of the coordinates of the center of mass Rx(t) and Ry(t). d) The center of mass trajectory. The color of the points represent the time. For other particles the movement is quite the same. The collective oscillatory motion seen244 has its origin in the coordination of the oscillations of the whole set. In our case the245 individual units do not have a chaotic behavior but instead is a harmonic function of the246 type: x=Acos(ω1·t) + Bsin(θ)cos(ω2·t), y=Asin(ω1·t) + Bsin(θ)sin(ω2·t) with247 ω2>> ω1.248 IV. CONCLUSION249 Swarming and synchronization is an ubiquitous natural phenomena that can be simulated250 by simple mathematical models. In this paper we studied a two-dimensional swarmalator251 model including a repulsive short-range interaction. Depending on the coupling parameter252 Kand the radius σof the interaction it is found six long-term behaviors most of them not253 17 a) b) 0.6 0.65 0.7 0.75 0.8 0.6 0.65 0.7 0.75 0.8 Dy Dx 0 0.1 0.2 0.3 0 0.1 0.2 Power f (Hz) c) d) -0.001 0 0.001 0.002 0.003 2000 2250 2500 Vx, Vy t 0 0.5 1.0 1.5 0 0.1 0.2 x10-4 Power f (Hz) FIG. 12. Swarmalator with a collective oscillation for the parameters K=−0.1, σ = 0.4. a) The plot of the major axis of the ellipse Dyversus the minor axis Dx.b) The power spectrum corresponding to Dx.c) Time series of the components Vxand Vyof the center of mass velocity. d) Power spectrum of Vx. previously reported. They are classified according according to the synchrony of the phases254 of the elements and their mobility in: a) static synchronous b) static wave c) dynamic255 wave d) dynamic synchronous e) dynamic asynchronous and f) splintered state. Especially256 interesting are those states that show a collective motion (P > 0) in which the whole set257 moves in a given direction. And among them those that show a collective oscillation for258 its potential use to simulate and understand many biological processes in which swarming,259 synchronization and rhythmicity play a fundamental role. A pertinent future goal would be260 to investigate if these states are realized in biological systems.261 262 18 a) b) -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 500 1000 1500 2000 x t -0.05 -0.025 0 0.025 0.05 500 1000 1500 2000 vx t c) d) -0.2 0 0.2 0.4 -0.2 0 0.2 0.4 y x -0.04 -0.02 0 0.02 0.04 -0.04-0.02 0 0.02 0.04 t vy vx 0*100 1*103 2*103 3*103 FIG. 13. Oscillatory behavior of a single unit for the parameters K=−0.1, σ = 0.4. Transient time has been discarded. a) Time series of the X coordinate. b) Time series of the x component of the velocity. c) The trajectory of in the plane x, y; and d) the plot of the components vyversus vx. The color refers to the time. ACKNOWLEDGMENT263 Many thanks to Kevin P. O’Keeffee for the encouragement and for reading the prelim-264 inary manuscript and making numerous suggestions that have greatly improved the final265 result.266 267 This work was supported by project TIN2017-80842-P cofinanced by Ministerio de268 Econom´ıa, Industria y Competitividad (MINECO) of Spain, through the Agencia Estatal269 de Investigaci´on (AEI), and by Fondo Europeo Regional (FEDER) of the EU.270 271 Notes: The simulations were done in an Intel I7 processor with ubuntu 18.04. Programs272 were written in C using the gcc 7.4.0 compiler. The contour plots of Figure 2 and Figure 3273 19 were obtained processing the data of Figure 1 with the gnuplot command dgrid3d 50,50274 gauss .2.275 REFERENCES276 [1] Agata Barci´s, Michal Barci´s, and Christian Bettstetter. Robots that sync and swarm: A proof277 of concept in ros-2. arXiv:1903.06440v4, 2019.278 [2] Avraham Be’er and Gil Ariel. A statistical physics view of swarming bacteria. Mov. Ecol.,279 7:9, 2019.280 [3] Abdelghani Bellouquid and Marcello Delitala. Mathematical Modeling of Complex Biological281 Systems. A Kinetic Theory Approach. Birk¨auser, Boston, 2001.282 [4] Andrea Cavagna, Lorenzo Del Castello, Supravat Dey, Irene Giandina, and Stefania Melillo.283 Short-range interactions versus long-range correlations in bird flocks. Physical Review E,284 92:012705, 2015.285 [5] Chong Chen, Song Liu, Xia-Qing Shi, Hugues Chat´e, and Yilin Wu. Weak synchronization286 and large-scale collective oscillation in dense bacterial suspensions. Nature, 242:210–227, 2017.287 [6] Yun Chen, Henk J. Busscher, Henry C. van der Mei, and Willem Norde. Statistical analysis of288 long- and short-range forces involved in bacterial adhesion to substratum surfaces as measured289 using atomic force microscopy. Applied and Environmental Microbiology, 77(15):5065–5070,290 2011.291 [7] John T. Emlen. Flocking behavior in birds. The Auk, 69-2:160–170, 1952.292 [8] Hyunsuk Hong. Active phase wave in the system of swarmalators with attractive phase293 coupling. Chaos, 28:103112, 2018.294 [9] Masamoto Iwasa, Kazuki Ida, and Dan Tanaka. Juggling motion in a system of motile coupled295 oscillators. Phys. Rev. E, 83:036210, 2011.296 [10] Masamoto Iwasa and Dan Tanaka. Dimensionality of clusters in a swarm oscillator model.297 Phys. Rev. E, 81:066214, 2010.298 [11] Y. Kuramoto. Chemical Oscillations, Waves and Turbulence. Springer, Berlin, 1984.299 [12] M.Belovs, R.Livanovics, and A. Cebers. Synchronized rotation in swarms of magnetotactic300 20 bacteria. Phys. Rev. E, 96:042408, 2017.301 [13] Joseph D. Monaco, Grace M. Hwang, Kevin M. Shultz, and Kechen Zhang. Cognitive swarm-302 ing: An approach from the theoretical neuroscience of hippocampal function. Proc. of SPIE,303 10982:109822D, 2019.304 [14] Kevin O’Keeffe and Christian Bettstetter. A review of swarmalators and their potential in305 bio-inspired computing. Proc. of SPIE, 10982:109822, 2019.306 [15] Kevin P.O’Keeffe, Joep H.M. Evers, and Theodore Kolokolnikov. Ring states in swarmalator307 systems. Phys. Rev. E, 98:022203, 2018.308 [16] Kevin P.O’Keeffe, Hyunsuk Hong, and Steven H. Strogatz. Oscillators that sync and swarm.309 Nat. Commun., 8:1504, 2017.310 [17] James G. Pukett, Douglas H. Kelley, and Nicholas T. Ouellete. Searching for effecive forces311 in laboratory insect swarms. Scientific Reports, 4:4776, 2014.312 [18] Dan Tanaka. General chemotactic model of oscillators. Phys. Rev. Lett., 99:134103, 2007.313 [19] Guy Theraulaz, Jacques Gautrais, and Jean-Louis Deneubourg. The formation of spatial314 patterns in social insects: from simple behaviours to complex structures. Philosophical Trans-315 actions. Series A, 361,1807:1263–82, 2003.316 [20] Chad M. Topaz and Andrea L. Bertozzi. Swarming patterns in a two dimensional kinematic317 model for biological groups. SIAM J. Appl. Math, 65:152–174, 2004.318 [21] Tams Vicsek and Anna Zafeiris. Collective motion. Phys.Rep, 517:71–140, 2012.319 21