Binary Mixtures of Locally Coupled Mobile Oscillators
BBinary Mixtures of Locally Coupled Mobile Oscillators
Gon¸calo Paulo
1, 2 and Mykola Tasinkevych
1, 2, ∗ Departamento de F´ısica, Faculdade de Ciˆencias,Universidade de Lisboa, 1749-016 Lisboa, Portugal Centro de F´ısica Te´orica e Computacional, Faculdade de Ciˆencias,Universidade de Lisboa, 1749-016 Lisboa, Portugal
We study synchronization dynamics in binary mixtures of locally coupled Kuramoto oscillatorswhich perform Brownian motion in a two-dimensional box. We introduce two models, where inmodel I there are two type of oscillators, say A and B , and any two similar oscillators tend tosynchronize their phases, while any two dissimilar ones tend to be out of phase. In model J , incontrast, the oscillators in subpopulation A behave as in model I , while the oscillators in subpop-ulation B tend to be out of phase with all the others. In the real space all the oscillators in bothmodels interact via a soft-core repulsive potential. Both subpopulations of model I and subpopu-lation A of model J , by their own, exhibit a phase coherent attractor in a certain region of modelparameters. The approach to the attractor, after an initial transient regime, is exponential withsome characteristic synchronization time scale τ . Numerical analysis reveals that the attractorsof the two subpopulations survive within model I , regardless of the composition of the mixture φ and the strength of the cross-population negative coupling constant H , and that τ sensitivelydepends on φ , H and the packing fraction. In particular, the ability of the oscillators to move andexchange neighbours can significantly decrease τ . In contrast, model J predicts suppression of thesynchronized state in subpopulation A and emergence of the coherent attractor in the “contrarians”subpopulation B for strong and weak cross-population coupling, respectively. I. INTRODUCTION
Synchronization phenomena has been studied for along time and is ubiquitous in numerous natural and arti-ficial systems, for example, in sociology, when people claptheir hands together [1], in ecology with the synchroniza-tion of frog choirs [2] or firefly blinking [3], in neurosciencewith the periodic spiking of autonomous pacemaker neu-rons [4, 5] and neural synchronization [6], in physiologywith biological rhythms [7], in active matter with the bac-teria self-propulsion [8], in nanotechnology with the syn-chronization of mechanical nanodevices [9], in condensedmatter physics with spin Hall nano-oscillators [10] etc.Synchronization can be thought about as the ”adjust-ment of rhythms of oscillating objects due to their weakinteraction” [11]. Oscillating objects which have beenmost studies are so-called self-sustained oscillators, i.e.,oscillators which can sustain their natural rhythm due tosome internal energy source and that are normally stableto small perturbations, returning to their original rhythmwhen left by themselves. These objects most frequentlyare modeled as limit-cycle oscillators whose behavior isgoverned by autonomous non-linear differential equations[12]. Depending on the topology of the interactions net-work of the oscillating objects as well as their physicalnature, various types of synchronization are possible suchas complete [13] and generalized [14] synchronization. In ∗ [email protected] this study we deal with a weaker form of synchroniza-tion, so-called phase synchronization [15, 16] which canbe realized in a system of coupled phase oscillators (theiramplitude is irrelevant in this case) and described by thecelebrated Kuramoto model [15].A basic understanding of the general conditions for theexistence and stability of the synchronized state of cou-pled oscillators is gained by applying the Master Stabil-ity Function formalism, which was developed by Pecora,Carroll and Barahona [17–19]. This method addressesa linear stability of the synchronized state of a genericsystem of linearly and symmetrically coupled dynamicalobjects, and allows to formulate a general criterion forthe synchronizability of the network links between thedynamical objects [forming the network nodes] indepen-dently of their specific physical nature [20]. For a givennature of the objects, the criterion sets a constraint onthe eigenvalues of the Laplacian matrix of the couplingnetwork in order to exhibit synchronized behavior [18].Besides considering oscillators with fixed interactionnetworks, researchers investigated synchronization intime-varying network topologies [21, 22] meant to ac-count for the effects of motility on coherent behavior ofautonomous objects, for example in swarming. It wasfound that networks with instantaneously disconnectedtopologies do allow for the synchronization of oscilla-tors, provided the topology varies in time in a certainway [22]. A dynamic network topology realized throughthe motion of locally coupled mobile agents was consid-ered in refs. [23–27]. A main common conclusion of these a r X i v : . [ c ond - m a t . s o f t ] F e b studies is that the mobility of the agents enhances syn-chronization. However, much richer behavior can be ob-tained when the interactions between the agents in thereal space are included, e.g., the synchronization timeexhibit a non-monotonic behavior as a function of themobility of agents with excluded volume interaction [27].In numerous instances however synchronization is un-desired [28, 29] as it can result in system malfunctioning[28, 29]. Examples include certain dysfunctions of thenervous system, such as epilepsy, Alzheimer and Parkin-son diseases, which are thought to be caused by excessiveneural synchronization [5, 6, 30–32], or dysfunctions in acomputer network due to an unintended synchronizationof the network routers [33]. Researchers addressed theissue of suppression of an undesired synchronization ofphase oscillators in the context of the Kuramoto modelwith non-identical interactions [28, 34–38]. Thus, by in-troducing to the system a given number of “contrarian”oscillators, i.e., the oscillators which tend to dephase withall the other oscillators in the system, it is possible tocompletely suppress the global synchronization [34], or tochange the nature of the synchronization transition froma first-order to second-order type [39]. Furthermore, bycoupling the phase and spatial dynamics of oscillators itis possible to achieve very rich spatiotemporal behaviorsuch as particle segregation according to their phases andphase waves, as reported recently in [40].Here we consider systems of phase oscillators whichperform Brownian motion in a two dimensional box withperiodic boundary conditions. Regarding the oscillatorphase dynamics, we subdivide the system into two sub-populations, which we name subpopulation A and B , re-spectively, and describe the dynamics of the oscillatorphases in terms of a generalized Kuramoto model withlocal coupling of the oscillators [related to a finite rangeof interactions in the real space] and with a three-valuedcoupling constant, where each of the values describes thetwo intra- and one cross-population phase interactions.The main goal of this study is to understand how the syn-chronization within a target subpopulation, say of type A , can be influenced by introducing “controlling agents”of type B . In particular, we are interested in identifyingunder which conditions the synchronization in the targetsubpopulation can be significantly suppressed or elimi-nated completely.The rest of the paper is organized as follows: in thenext section we introduce a model for Brownian dynam-ics of the oscillators in the real space and two differentmodels for the phase dynamics. We also define synchro-nization order parameters and discuss some analyticalresults of the continuum mean-field approximation forthe standard Kuramoto model. In section III we discussthe simulation results emphasizing qualitative differencesbetween the two models. In the last section we summa-rize our results and discuss the potential extension of the present work. II. MODELS: DYNAMICS OF THEOSCILLATORS IN THE REAL AND PHASESPACESA. Brownian motion of oscillators with excludedvolume
We consider N oscillators of mass m in a two-dimensional (2D) box of size L × L and assume peri-odic boundary conditions in both spacial directions. Thetime-dependent position of the center of mass of oscillator i is denoted by r i ( t ) ≡ ( x i ( t ) , y i ( t )) T , where i = 1 , .., N and t stands for time. Without loss of generality we as-sume that the oscillators are dispersed in some fluid envi-ronment and that their mutual repulsion is described byinteraction potential U ( r , ..., r N ). The time evolutionof r i ( t ) is described by a system of coupled Langevinequations: m d r i dt = − γ d r i dt − ∇ i U ( r , ..., r N ) + (cid:112) γk B T ξ i ( t ) , (1)where the first term on the r.h.s. describes the frictionforce (due to, e.g., interaction with the surrounding am-bient fluid) with the friction coefficient γ , k B denotesthe Boltzmann constant and T is the absolute temper-ature of the ambient fluid. γ is related to the diffu-sion coefficient D by the fluctuation-dissipation theorem D = k B T /γ . The second term gives the force on oscil-lator i originating from the repulsive interactions withthe other oscillators characterized by the interaction po-tential U , ∇ i stands for the gradient operator actingupon r i . Thermal fluctuations in surroundings are ac-counted for by a random force ξ i ( t ), which is a zero-mean delta-correlated Gaussian process, i.e., (cid:104) ξ i ( t ) (cid:105) = 0, (cid:104) ξ i ( t ) ξ Tj ( t ) (cid:105) = δ ij δ ( t − t ) with the identity matrix, δ the Dirac delta function, and (cid:104) ... (cid:105) denoting an ensem-ble average. We approximate U as a pairwise sum ofLennard-Jones truncated and shifted potential U ( r , ..., r N ) = (cid:15) N (cid:80) i
2, and σ is the effective diameter of theoscillators.We have integrated Eqs. (1) numerically by employingLarge-scale Atomic/Molecular Massively Parallel Simula-tor (LAMMPS) package [41] at constant volume V = L and temperature T = 1 . (cid:15)/k B . The molecular dynamicssimulations were run with a time step ∆ t = 0 . (cid:112) mσ /(cid:15) .The temperature was kept constant by using Langevinthermostat as described in [42] with the Langevin fric-tion coefficient γ = 1 . (cid:112) m(cid:15)/σ modeling a backgroundimplicit solvent. These values of T and γ set the diffusionconstant D = 1 . (cid:112) σ (cid:15)/m . B. Generalized Kuramoto model
1. The standard Kuramoto model and the synchronizationorder parameter
Yoshiki Kuramoto introduced his celebrated model in1975 [15] in response to the earlier work on synchroniza-tion of biological rhythms by Winfree in 1967 [43]. TheKuramoto model can be solved exactly within certainlimits, exhibiting a synchronization phase transition insystems of phase oscillators. This model is widely usedto study synchronization phenomena in both dynamicaland static networks of interacting oscillators [44, 45]. Fora system of N all-to-all interacting phase oscillators, thedynamics of the phase θ k ( t ) of oscillator k is postulatedin the following form dθ k ( t ) dt = ω k + K N N (cid:88) l =1 sin( θ l ( t ) − θ k ( t )) , (3)where ω k is the natural frequency of the k th oscillatorwhen it is not interacting with any other oscillator and K is a coupling constant. The synchronization oftenis analyzed in terms of a complex order parameter z ( t )defined as z ( t ) = r ( t ) e i Ψ( t ) = 1 N N (cid:88) k =1 e iθ k ( t ) , (4)with r ( t ) quantifying the global oscillator coherence, andΨ( t ) being the average phase; here i = √− r ( t ) = 1,with r ( t ) = 0 in the opposite limit of totally dephasedoscillators.In the limit N → ∞ and for ω k = 0 , k = 1 , .., N Ottand Antonsen derived the following mean-filed equationfor r ( t ) [46] dr ( t ) dt = K (cid:0) r ( t ) − r ( t ) (cid:1) , (5)which for a given initial condition r (0) has the solution r ( t ) = 1 √ − K ( t − t ∗ ) , (6) with t ∗ = ln (cid:16) r (0) − (cid:17) K . (7)Eq. (6) is the square root of the logistic function andtheir shapes are quit similar. After some initial transientregime, K − defines a characteristic synchronization timescale on which r ( t ) reaches its saturation value of 1. In-deed, for t (cid:29) K − r ( t ) ∼ −
12 e − K ( t − t ∗ ) . (8) t ∗ marks the location of the inflection point of r ( t ). Be-low we show, that in certain regimes the synchronizationcurves in our systems can be very accurately describedby Eq. (6). (a) t . . . . . r A Sim.
F it (b)(c) t . . . . . r A Sim.Fit t r A (d) Figure 1. Panels (a) and (c) show r A ( t ) for systems at ζ = 1and ζ = 100, respectively, other parameters are N = 1800, L = 50 σ , φ = 0 . K = 0 . (cid:112) (cid:15)/mσ . Different curvesrepresent different system realizations, i.e., different initial θ σk ( t = 0). Large fluctuations of r A ( t ) for intermediate timesin (c) are due to the annihilation dynamics of topologicaldefects in the θ field, see Figs. 2(f) and (g). The red curvesin (a) and in the inset of (c) represent fits to the mean-fieldmodel r A ( t ) = 1 / √ e − κ ( t − t ) . In (c) the fitting was carriedout only within the time range that corresponds to r A ≥ . κt (cid:28)
1, behaviorof 1 − r A ( t ) with exponential decay 1 − r A ( t ) ∼ e − t/τ . Wefind that κ − (cid:54) = τ . Time is given in units of (cid:112) mσ /(cid:15) .
2. Models for binary mixtures of phase oscillators
We subdivide the system of N phase oscillators in N A oscillators of type A , whose phases we denote by θ A k , andin N B oscillators of type B , whose phases we denote by θ B k , and N = N A + N B . We consider that A− oscillatorsform the target subpopulation whose synchronization be-havior will be manipulated through the addition of thecontrol oscillators of type B . The governing equations ofour model are as follows dθ σk dt = ω σk + 1 N R (cid:88) σ (cid:48) = A , B K σσ (cid:48) N σ (cid:48) (cid:88) l =1 ,r kl
H < K drivesynchronization within a given subpopulation, while thenegative H forces distinct oscillators to proceed out ofphase.In model J we assume that each B− oscillator has neg-ative coupling to all the other N oscillators in the systemand set (cid:20) K AA K AB K BA K BB (cid:21) = (cid:20) K HH H (cid:21) . (11)We define composition φ of the binary mixture as φ = N B / ( N A + N B ) which controls the average number ofrepulsive pairwise contacts in the system. We also in-troduce packing fraction η = N πσ / L , and parame-ter ζ ≡ | H | /K quantifying the interaction asymmetryof mixture. Below we report the effects of these threeparameters on the synchronization of subpopulation A which is quantified by the order parameter r A ( t ) = Re N σ N A (cid:88) k =1 e iθ k ( t ) . (12)The order parameter for subpopulation B can be definedin a similar way. We integrate equations (9) numerically by using theEuler method with the same time step as is used tointegrate Eqs. (1), i.e., ∆ t = 0 . (cid:112) mσ /(cid:15) . Specifi-cally, we have implemented this integration in a ficti-tious LAMMPS pair style function, which is then su-perimposed with the particle pair potential in the realspace. The Kuramoto dynamics, as given by Eqs. (9), isswitched on only after the system of soft disks, whose dy-namics in the real space is governed by Eqs. (1), reachesthermal equilibrium. Initial phases θ k ( t = 0) , k = 1 , .., N are set to random values drawn from uniform distributionon [0 , π ). III. RESULTSA. Model I We introduce the time scale of local synchronization t θ = min( K − , H − ), and the second timescale t r asso-ciated with the rate of change of the network topologydue to the motion of oscillators t r = R /D . t r char-acterizes an average time that an oscillator requires todiffuse over a distance of the interaction range R . Equiv-alently t r can be thought of as an average duration ofthe interaction between the phases of two oscillators be-fore they diffuse apart beyond the interaction range. Inour case t r = 9 (cid:112) mσ /(cid:15) . The limit t θ (cid:29) t r correspondsto so-called fast switching [FS] regime in the language ofdynamic networks, when the network links change muchfaster as compared to the local phase reorientation. Pre-vious studies of the synchronization of identical mobileoscillators [locally coupled] revealed that in the FS regimeone can carry out time averaging of the original time-dependent adjacency matrix, which renders a mappingof the phase dynamics of mobile oscillators onto thatof an effective all-to-all static Kuramoto model [23, 25].Therefore, in the FS limit the model can be treated atthe mean-field level [27].Model I exhibits similar mean-field like behavior inthe FS regime. Indeed, numerical results obtained forvalues of coupling constants 0 < K, | H | < .
05 revealthat the time evolution of the order parameter r A ( t ) hasa functional form which agrees semi-quantitatively withthe mean-field solution of Ott and Antonsen in Eq. (6),i.e., r A ( t ) ∼ / √ e − κ ( t − t ) as shown in figure 1(a).The parameter κ is related [but not equal] to the inverseof the synchronization time τ , which characterizes thelate time exponential behavior of 1 − r A ( t ) [see Fig. 1(b)]; t defines the location of the inflection point. For the spe-cific case depicted in Fig. 1(a) we find by fitting the nu-merical data to the mean-field model that κ ≈ . × − , t ≈ . × . Figures 2(b)-(d) show snapshots of thesystem evolution in the FS regime. The network con- (a) t . . . . . r A (b) (c) (d)(e) t . . . . . r A (f) (g) (h) sin (2 ) Figure 2. (a) r A ( t ), (b)-(d) the instantaneous oscillator configurations at times matching those marked by vertical lines in(a). The represented system corresponds to fast switching regime, with N = 1800, L = 50 σ , φ = 0 . K = 0 . (cid:112) (cid:15)/mσ and ζ = 1. Panels (e)-(h), the same as above but at ζ = 100. Both A− and B− oscillators are shown and colored accordingto sin ( θ σk ). We choose this color scheme in order to emphasize the presence of the topological defects. Both systems reachglobally synchronized configurations, r ≈ (a) .
00 0 .
25 0 .
50 0 .
75 1 . φ . . . . κ ζ = 1 ζ = 5 ζ = 10 (b) .
00 0 .
25 0 .
50 0 .
75 1 . φ t ζ = 1 ζ = 5 ζ = 10 Figure 3. (a) κ , (b) t as functions of the composition φ of thebinary mixture for several values of ζ . κ and t are plotted inunits of (cid:112) (cid:15)/mσ and (cid:112) mσ /(cid:15) , respectively. In all the cases, K = 0 . (cid:112) (cid:15)/mσ , N = 1800 and L = 50 σ . The results havebeen obtained after averaging over 100 independent runs. nectivity changes rapidly due to the fast diffusion of theoscillators, and the phases “experience” a self-averagedeffective all-to-all topology. All phases approach the com- (a) . . . . . ζ . . . . κ φ = 0 . φ = 0 . φ = 0 . (b) . . . . ζ t φ = 0 . φ = 0 . φ = 0 . Figure 4. (a) κ and (b) t as functions of ζ for severalvalues of the mixture composition φ . κ is plotted in units of (cid:112) (cid:15)/mσ , and t in units of (cid:112) mσ /(cid:15) . In all the case K =0 . (cid:112) (cid:15)/mσ , N = 1800 and L = 50 σ . The results areobtained after averaging over 100 independent runs. plete synchronization with similar rates. No significantfluctuations in the local ordering are visible and the orderparameter approximately follows the mean-field dynam- (a) (b) (c) (d)(e) (f) (g) (h)(i) (j) (k) (l) sin (2 ) Figure 5. Typical instantaneous oscillator configurations obtained in the moderate switching regime at K = 0 . (cid:112) (cid:15)/mσ , whichcorresponds to t r /t θ = 0 .
9. (a)-(d) correspond to ζ = 20 and φ = 0 .
05; (e)-(h) to ζ = 20 and φ = 0 .
5; and (i)-(l) to ζ = 0 . φ = 0 .
5. Panels (a), (e), (i) depict the late time behavior of 1 − r A ( t ) ∼ e − t/τ . The snapshots in each of the rows are taken attimes which match those marked by vertical blue lines on panels (a), (e) and (i), respectively. For all the cases N = 7200 and L = 100 σ , which corresponds to the packing fraction η ≈ .
57. Time is given in units of (cid:112) mσ /(cid:15) . ics of Eq.(6). Similar behavior was found in ref. [25] andwas named “global synchronization”.Different behavior is observed for 0 . < K, | H | whichis depicted in Figs. 1(c) and (d). Now t r /t θ > .
45 suchthat the local phase dynamics competes with the dynam-ics of the network reconfiguration. In the limit t r /t θ (cid:29) D XY model after quenching from the infinite tem-perature to the values below the Kosterlitz-Thouless crit-ical point [27, 47]. The system snapshots in Figs. 2(f)-(h) demonstrate the spatially non-homogeneous struc-tures with locally synchronized regions which grow withtime similarly to the coarsening dynamics of XY model.We also observe several topological defects, which ap-pear when growing regions with different average phasesmeet. Three different time regimes can be clearly iden-tify in r A ( t ) in Fig. 1(c). Approximately linear regimefor t ≤
200 where all the curves approximately collapseon top of each other. This is the regime when the lo-cally coherent regions appear in different locations andgrow in an overall disordered background. When theseregions start touching, the orientational topological de-fects emerge and initiate the second dynamical regime for200 ≤ t ≤ r A ( t ) curves. At later times, t ≥ r A ( t ) can be very well approximated by themean-field profile in Eq. (6), see red dashed curves in theinset of Fig. 1(c). (a) .
00 0 .
25 0 .
50 0 .
75 1 . φ τ ζ =100 ζ =0.1 ζ =1 (b) − ζ τ φ α ( φ ) φ =0.50 φ =0.15 φ =0.05 Figure 6. (a) Synchronization time τ as a function of φ forseveral values of ζ . (b) τ as a function of ζ for several values of φ . Solid dotted lines represent fits to a power-law τ ∝ ζ α ( φ ) ,only data points for ζ >
10 are fitted. The resulting exponent α is shown as a function of φ in the inset. In all the cases N = 1800, L = 50 σ , K = 0 . (cid:112) (cid:15)/mσ . The results areobtained after averaging over 100 independent runs. As we discussed above, in the FS regime the system ex-hibits a mean-field behavior, and it is expected that forpure A− system, φ = 0, κ ∝ K and t ∝ K − [25]. Theaddition of the oscillators of type B changes this behaviorsuch that both κ and t become non-trivial functions of K , ζ , φ and η . We start by plotting in Fig. 3 κ and t as functions of φ at several values of ζ and at constant K and η . All the systems represented belong to the FSregime. Surprisingly, κ and t do not depend of the num-ber of B− oscillators introduced [within the numerical un-certainty] for ζ = 1. In this case both subpopulationssynchronize at the same rate, but with the average phasesthat are diametrically opposed to each other. By increas-ing ζ the phase repulsion between locally synchronized A− and B− domains increases. This evokes a kind of apositive feedback mechanism when a given locally coher-ent cluster acquires new members mostly because thoseare being repelled by a neighbouring coherent cluster ofanother type. This effect is more pronounced the largerthe two oppositely synchronized clusters are, such thateventually the global synchronization is driven by therepulsion, rapidly [for large ζ ] driving the two subpopu- lations into diametrically opposed states, decreasing theeffective synchronisation time τ ∼ κ − of A− oscillators,see blue triangles and red squares in Fig. 3(a). The lo-cation of the inflection point t also decreases with in-creasing ζ . The results show that there is an optimumcomposition of the mixture φ min ≈ . κ − and t attain their minimum values. Based on the symmetryarguments, one would expect rather symmetric shapesfor both τ ( φ ) and t ( φ ) curves, with the respective ex-trema located at φ min ≈ .
5, similarly to the case of theslow switching regime discussed below [see Fig. 6]. Werelate this slight asymmetry in the present case to thefinite size effects. (a) . . . . η . . . . τ \ N N = 1800 N = 7200 (b) . . . . η τ ζ =0.1 ζ =1 ζ =10 ζ =100 Figure 7. (a) Synchronization time τ as a function of thepacking fraction η for two values of N . Rescaling the verticalaxis by N − results in an approximate data collapse. Thepacking fraction is varied by changing the system size L ; K =0 . (cid:112) (cid:15)/mσ , ζ = 1, φ = 0 .
5. (b) τ as a function of η forseveral values of ζ . K = 0 . (cid:112) (cid:15)/mσ , φ = 0 . N = 1800. τ is given in units of (cid:112) mσ /(cid:15) . The results are obtained afteraveraging over 100 independent runs. Figure 4 summarizes the dependence of κ and t onthe interactions asymmetry ζ for several values of φ . Weemphasize that the presence of just a 5% of controllingagents of type B with the interaction asymmetry ζ = 10reduces κ − of the control population A by a factor of 2,see green circles in Fig. 4(a).In the moderate switching regime, when K or | H | islarger than 0 . (cid:112) (cid:15)/mσ we compute the synchronizationtime τ by fitting the late time behaviour of 1 − r A ( t ) tothe exponential time decay ∝ e − t/τ , as is demonstratedin Fig. 1(d) and in the first column of Fig.5. In Fig. 5we compare the synchronization dynamics of a A− richsystem at φ = 0 . , ζ = 20 [1st row], and 50 /
50 mix-tures of A− and B− type oscillators at ζ = 20 [2nd row]and ζ = 0 . ζ = 20] or increase [ ζ = 0 .
1] the synchronizationtime τ of the target A− subpopulation. Similar to thecase depicted in Fig. 2, here we also observe the coars-ening dynamics driven by vortex-antivortex annihilationreminiscent of the 2 D XY model. The defects are morepronounced in these simulations because we have used alarger number of oscillators. (a) t . . . . . r A H = -0.001 H = -0.001 H = -0.1 H = -0.001 H = -0.1 H = -1 H = -0.001 H = -0.1 H = -1 H = -10 (b) t . . . . . r B H = -0.001 H = -0.001 H = -0.1 H = -0.001 H = -0.1 H = -1 H = -0.001 H = -0.1 H = -1 H = -10 Figure 8. (a) r A , (b) r B as functions of time for several valuesof H . N = 1800, L = 100 σ , K = 0 . (cid:112) (cid:15)/mσ , and φ = 0 . A− oscillators exhibit complete synchronization with r ∗A = 1,and B− oscillators we find r ∗B > . | H | < . (cid:112) (cid:15)/mσ . r ∗A decreases systematically with | H | for | H | > . (cid:112) (cid:15)/mσ .Time is shown in units of (cid:112) mσ /(cid:15) . We find that in moderate switching regime the syn-chronization time τ , in contrast to the case of the FSregime shown in Fig. 3(a), is a symmetric function of φ with the minimum [maximum] at φ = 0 . ζ > ζ < φ with 1 − φ does not affect the model be-havior. We emphasize, however, that the mechanismsdriving the synchronization of A− oscillators are differ-ent for small and large φ . Indeed, for φ (cid:28) A− subpopulation is driven by the in-trapopulation attractions, while for φ ≈ B− subpopulation whichpromotes the synchronization of A− oscillators. For ex-ample, for the system of Fig. 6(a) the average distancebetween A− oscillators is ≈ . σ at φ = 0 .
95 which isalmost twice the interaction range, therefore the contri-bution of the attractive intrapopulation interactions isinsignificant.We observe that τ does not depend [up to numericaluncertainty] on the composition of the mixture at ζ = 1.Our numerical results also suggest that for ζ (cid:38)
10 thesynchronization time as a function of ζ has a power-lawform with an exponent α that depend on φ and with theminimum at φ = 0 .
5, as shown in Fig. 6(b).In the slow switching regime, when the network is ef-fectively static, one can argue using the dynamic scalinghypothesis [48] that the synchronization time τ ∼ N [27].Surprisingly, we find that this scaling also holds approx-imately in the moderate switching regime. In Fig. 7(a)we plot τ /N as a function of the packing fraction for twodifferent N , showing an approximate data collapse ontoa single master curve. Figure 7(b) also shows that there (a) .
05 0 .
30 0 .
50 0 .
80 0 . φ − . − . − − H . . . . r ∗ A (b) .
05 0 .
30 0 .
50 0 .
80 0 . φ − . − . − − H . . . . r ∗ B Figure 9. Heat map of the late time values of the order param-eter for (a) A− oscillators and (b) B− oscillators. The mapsshow the existence of configurational transitions from fullysynchronized r ∗ σ = 1 to disordered r ∗ σ (cid:28) φ or H . All simulations were carried out at K = 0 . (cid:112) (cid:15)/mσ , N = 1800, L = 100. is a ζ − dependent optimum value η min ( ζ ) of the packingfraction at which τ is minimal. At the lowest value ofthe packing fraction considered, η = 0 .
06, the averagedistance between the oscillators is ≈ . τ . Contrary, at higher values of η the mobility of the oscillators decreases due to the stericjamming of the system, resulting in a slower mixing andlarger τ . Indeed, we find numerically that for this systemthe diffusion coefficient varies over an order of magnitudefor 0 . ≤ η ≤ .
6. Surprisingly, the effect of the stericjamming can be overcome by increasing ζ , see red squaresin Fig. 7(a) B. Model J While model I is symmetric with respect to swappingthe type of the oscillators and replacing φ by 1 − φ , inmodel J this symmetry is lifted. B− oscillators are con-trarians to all N particles in the system. Similar models,but on static random networks, were considered by sev-eral groups [28, 34–38] with a common conclusion thatit is possible to fully suppress the global synchroniza-tion, provided the number of contrarians or their nega-tive coupling strength exceeds certain thresholds. Ourresults are in agreement with this general conclusion. (a) .
25 0 .
50 0 .
75 1 φ . . . . . r ∗ A H = − . H = − H = − (b) .
25 0 .
50 0 .
75 1 φ . . . . . r ∗ B H = − . H = − H = − Figure 10. (a) r ∗A , (b) r ∗B as functions of φ for several valuesof H . N = 1800, L = 100 σ , H = − (cid:112) (cid:15)/mσ and K =0 . (cid:112) (cid:15)/mσ . Figure 8 shows the time evolution of the order param-eters r A ( t ) and r B ( t ) for both subpopulations at φ = 0 . H . Surprisingly, for low values of | H | we observe anemergence of a partially synchronized states for contrar-ians with the late time values of the order parameter r B > .
8, see H = − .
001 curve in figure 8(b). Increasingthe repulsive strength | H | reduces the asymptotic value r ∗ σ ≡ r σ ( t → ∞ ) for both subpopulations.Next we discuss the dependence of r ∗ A and r ∗ B upon H , K and φ . On Fig. 9 we show heat maps of the orderparameters r ∗ A and r ∗ B in the ( φ, H ) plane. The mapsreveal the presence of continuous configurational transi-tions governed by either H or φ , at which r ∗ σ , σ = A, B decreases from 1 to 0 with decreasing H at a certain fixed φ , or with increasing φ at a certain fixed H .In Fig. 10 we plot r ∗A and r ∗B for several cuts throughthe heat map at several fixed values of H . Curiously, r ∗A ( φ ) [ r ∗B ( φ )] resembles the behavior of the spontaneous[field-driven] magnetization as a function of temperaturein the Ising model, with φ playing the role of tempera-ture. IV. CONCLUSIONS
We have carried out extensive numerical simulationsof the synchronization dynamics of two models of locallycoupled mobile oscillators of two types. The first modeldeals with symmetric binary mixtures, where alike oscil-lators tend to synchronize, while unlike ones tend to beout of phase. The second model is asymmetric in this re-spect as it contains a given fraction of contrarians whichtend do be out of phase with all the other oscillators inthe system.We have focused in model I on the characteristic syn-chronization time τ describing the asymptotic exponen- tial approach of the order parameter in Eq. (12) to unity,see also Figs. 1(b) and (d). We have found that τ de-creases with increasing the composition φ of the mix-ture and for | H | > K , i.e., when the repulsive interac-tions are more intense as the attractive ones [Fig. 3(a)and Fig. 6(a)]. τ is also very sensitive to the ratio ζ = | H | /K , and exhibits a power low decay for largeenough ζ [Fig. 6(b)]. The synchronization dynamics pro-ceed approximately as predicted by the mean-field the-ory, Eq. (6), in the fast switching regime which in ourcase is realized for K , and | H | < . (cid:112) (cid:15)/mσ [Fig. 1(a),and Fig. 2(a)]. For larger values of the coupling con-stants, the dynamics of the phase synchronization canbe mapped onto the coarsening relaxation dynamics ofthe 2 D XY model [Figs. 2(e)-(h) and Figs. 5], which inthe intermediate times is dominated by the motion andannihilation of topological defect with opposite windingnumbers. Controlling the overall number density of theoscillators provides additional means to enhance the syn-chronization, highlighting the role of the mobility of theoscillators [Fig. 7]. The general conclusion for model I isthe following: if a pure subpopulation exhibits a coherentattractor it will also be present in the mixtures.In contrast, model J allows for a complete suppressionof the coherent state, which can be achieved by increas-ing either ζ or φ , as shown in the configuration diagramsin Fig. 9. For K >
0, the late time values of both or-der parameters, r ∗A , and r ∗B vary continuously between 1and 0 with the increasing of the mixture composition φ [Fig. 10].Here we have used a simplifying assumption that allthe intrinsic frequencies ω σk = 0, see Eq. (9). It is ex-pected that the results obtained here will also hold forthe case of non-zero, but uniform frequencies ω σk = ω .This however will not be the case if the frequencies arechosen at random from some distributions with densities g σ ( ω ), σ = A , B . Another possible extension of this workis to study the effects of coupling between the phasesand the coordinates of the oscillators, e.g., by requiringthat the oscillators move in the directions determined bytheir phase variables. Then the goal is to discover somenovel swarming-like behaviors in the binary mixtures ofself-propelled oscillators with the heterogeneous aligninginteractions. One can also search for the conditions whenbinary systems demix into A− rich and B− rich subpop-ulations, when for example the means of their respective g σ ( ω ) have opposite signs. ACKNOWLEDGMENTS
We acknowledge financial support from the Por-tuguese Foundation for Science and Technology (FCT)under Contracts no. PTDC/FIS-MAC/28146/2017(LISBOA-01-0145-FEDER-028146), UIDB/00618/2020,UIDP/00618/2020, and IF/00322/2015.0 [1] Z. N´eda, E. Ravasz, T. Vicsek, Y. Brechet, and A. L.Barab´asi, Physics of the rhythmic applause, Phys. Rev.E , 6987 (2000).[2] I. Aihara, H. Kitahata, K. Yoshikawa, and K. Aihara,Mathematical modeling of frogs’ calling behavior and itspossible application to artificial life and robotics, Artifi-cial Life and Robotics , 29 (2008).[3] G. M. Ram´ırez- ´Avila, J. Kurths, S. Depick`ere, and J.-L. Deneubourg, Modeling fireflies synchronization, in AMathematical Modeling Approach from Nonlinear Dy-namics to Complex Systems (Springer International Pub-lishing, 2018).[4] D. Plenz and S. T. Kital, A basal ganglia pacemakerformed by the subthalamic nucleus and external globuspallidus, Nature , 677 (1999).[5] D. J. Surmeier, J. N. Mercer, and C. S. Chan, Au-tonomous pacemakers in the basal ganglia: who needsexcitatory synapses anyway?, Curr. Opin. Neurobiol. ,312 (2005).[6] R. G. Erra, J. L. P. Velazquez, and M. Rosenblum, Neu-ral synchronization from the perspective of non-linear dy-namics, Front. Comput. Neurosci. (2017).[7] L. Glass, Synchronization and rhythmic processes inphysiology, Nature , 277 (2001).[8] O. A. Igoshin, A. Mogilner, R. D. Welch, D. Kaiser,and G. Oster, Pattern formation and traveling waves inmyxobacteria: Theory and modeling, Proc. Natl. Acad.Sci. U. S. A. , 14913 (2001).[9] D. Antonio, D. A. Czaplewski, J. R. Guest, D. L´opez,S. I. Arroyo, and D. H. Zanette, Nonlinearity-inducedsynchronization enhancement in micromechanical oscil-lators, Phys. Rev. Lett. (2015).[10] A. A. Awad, P. D¨urrenfeld, A. Houshang, M. Dvornik,E. Iacocca, R. K. Dumas, and J. ˚Akerman, Long-rangemutual synchronization of spin hall nano-oscillators, Na-ture Physics , 292 (2016).[11] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchroniza-tion: A Universal Concept in Nonlinear Sciences , Cam-bridge Nonlinear Science Series (Cambridge UniversityPress, 2001).[12] A. Barrat, M. Barth´elemy, and A. Vespignani,
DynamicalProcesses on Complex Networks (Cambridge UniversityPress, 2008).[13] L. M. Pecora and T. L. Carroll, Synchronization inchaotic systems, Phys. Rev. Lett. , 821 (1990).[14] N. F. Rulkov, M. M. Sushchik, L. S. Tsimring, andH. D. I. Abarbanel, Generalized synchronization of chaosin directionally coupled chaotic systems, Phys. Rev. E , 980 (1995).[15] Y. Kuramoto, Lecture notes in physics, internationalsymposium on mathematical problems in theoreticalphysics., Springer-Verlag , 420 (1975).[16] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, Phasesynchronization of chaotic oscillators, Phys. Rev. Lett. , 1804 (1996).[17] L. M. Pecora and T. L. Carroll, Master stability functionsfor synchronized coupled systems, Phys. Rev. Lett. ,2109 (1998).[18] M. Barahona and L. M. Pecora, Synchronization in small-world systems, Phys. Rev. Lett. , 054101 (2002). [19] L. M. Pecora and T. L. Carroll, Master stability func-tion for globally synchronized systems, in Encyclopediaof Computational Neuroscience , edited by D. Jaeger andR. Jung (Springer New York, 2013) pp. 1–13.[20] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, andC. Zhou, Synchronization in complex networks, PhysicsReports , 93 (2008).[21] J. D. Skufca and E. M. Bollt, Communication andsynchronization in disconnected networks with dynamictopology: Moving neighborhood networks, MBE , 347(2004).[22] D. J. Stilwell, E. M. Bollt, and D. G. Roberson, Suf-ficient conditions for fast switching synchronization intime-varying network topologies, SIAM , 140 (2006).[23] M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna, andS. Boccaletti, Synchronization of moving chaotic agents,Phys. Rev. Lett. , 044102 (2008).[24] F. Peruani, E. M. Nicola, and L. G. Morelli, Mobilityinduces global synchronization of oscillators in periodicextended systems, New Journal Phys. , 093029 (2010).[25] N. Fujiwara, J. Kurths, and A. D´ıaz-Guilera, Synchro-nization in networks of mobile oscillators, Phys. Rev. E (2011).[26] K. Uriu and L. G. Morelli, Collective cell movement pro-motes synchronization of coupled genetic oscillators, Bio-phys. J.l , 514 (2014).[27] D. Levis, I. Pagonabarraga, and A. D´ıaz-Guilera, Syn-chronization in dynamical networks of locally coupledself-propelled oscillators, Phys. Rev. X , 011028 (2017).[28] V. H. P. Louzada, N. A. M. Ara´ujo, J. S. Andrade, andH. J. Herrmann, How to suppress undesired synchroniza-tion, Sci. Rep , 658 (2012).[29] E. L. Lameu, F. S. Borges, R. R. Borges, K. C. Iarosz,I. L. Caldas, A. M. Batista, R. L. Viana, and J. Kurths,Suppression of phase synchronisation in network basedon cat’s brain, Chaos , 043107 (2016).[30] P. J. Uhlhaas and W. Singer, Neural synchrony inbrain disorders: Relevance for cognitive dysfunctions andpathophysiology, Neuron , 155 (2006).[31] M. L. Kringelbach, N. Jenkinson, S. L. Owen, and T. Z.Aziz, Translational principles of deep brain stimulation,Nat. Rev. Neurosci. , 623 (2007).[32] G. Deuschl, C. Schade-Brittinger, P. Krack, J. Volkmann,H. Sch¨afer, K. B¨otzel, C. Daniels, A. Deutschl¨ander,U. Dillmann, W. Eisner, D. Gruber, W. Hamel, J. Her-zog, R. Hilker, S. Klebe, M. Kloß, J. Koy, M. Krause,A. Kupsch, D. Lorenz, S. Lorenzl, H. M. Mehdorn, J. R.Moringlane, W. Oertel, M. O. Pinsker, H. Reichmann,A. Reuß, G.-H. Schneider, A. Schnitzler, U. Steude,V. Sturm, L. Timmermann, V. Tronnier, T. Trottenberg,L. Wojtecki, E. Wolf, W. Poewe, and J. Voges, A ran-domized trial of deep-brain stimulation for parkinson ' sdisease, N. Engl. J. Med. , 896 (2006).[33] S. Floyd and V. Jacobson, The synchronization of pe-riodic routing messages, IEEE/ACM Trans. Netw. ,122–136 (1994).[34] D. H. Zanette, Synchronization and frustration in oscil-lator networks with attractive and repulsive interactions,EPL , 190 (2005).[35] H. Hong and S. H. Strogatz, Kuramoto model of coupledoscillators with positive and negative coupling parame- ters: An example of conformist and contrarian oscilla-tors, Phys. Rev. Lett. , 054102 (2011).[36] H. Hong and S. H. Strogatz, Conformists and contrariansin a kuramoto model with identical natural frequencies,Phys. Rev. E , 046202 (2011).[37] M. Mirchev, L. Basnarkov, F. Corinto, and L. Kocarev,Cooperative phenomena in networks of oscillators withnon-identical interactions and dynamics, IEEE Trans.Circuits Syst. I, Reg. Papers , 811 (2014).[38] I. Ratas and K. Pyragas, Eliminating synchronization inbistable networks, Nonlinear Dyn. , 1137 (2016).[39] X. Zhang, S. Guan, Y. Zou, X. Chen, and Z. Liu, Sup-pressing explosive synchronization by contrarians, EPL , 28005 (2016).[40] K. P. O’Keeffe, H. Hong, and S. H. Strogatz, Oscillatorsthat sync and swarm, Nature Communications (2017).[41] S. Plimpton, Fast parallel algorithms for short-rangemolecular dynamics, J. Comp. Phys. , 1 (1995).[42] T. Schneider and E. Stoll, Molecular-dynamics study ofa three-dimensional one-component model for distortive phase transitions, Phys.l Rev. B , 1302 (1978).[43] A. T. Winfree, Biological rhythms and the behavior ofpopulations of coupled oscillators, J. Theor. Biol. , 15(1967).[44] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths,The kuramoto model in complex networks, Physics Re-ports , 1 (2016).[45] R. Schmidt, K. J. R. LaFleur, M. A. de Reus, L. H.van den Berg, and M. P. van den Heuvel, Kuramotomodel simulation of neural hubs and dynamic synchronyin the human cerebral connectome, BMC Neurosci. (2015).[46] E. Ott and T. M. Antonsen, Low dimensional behavior oflarge systems of globally coupled oscillators, Chaos ,037113 (2008).[47] J. M. Kosterlitz, The critical properties of the two-dimensional xy model, J. Phys. C , 1046 (1974).[48] P. C. Hohenberg and B. I. Halperin, Theory of dynamiccritical phenomena, Rev. Mod. Phys.49