Synchronization of moving oscillators in three dimensional space
Abstract
We investigate macroscopic behavior of a dynamical network consisting of a timeevolving wiring of interactions among a group of random walkers. We assume that each walker (agent) has an oscillator and show that depending upon the nature of interaction, synchronization arises where each of the individual oscillators are allowed to move in such a random walk manner in a finite region of three dimensional space. Here the vision range of each oscillator decides the number of oscillators with which it interacts. The live interaction between the oscillators is of intermediate type ( i.e., not local as well as not global) and may or may not be bidirectional. We analytically derive density dependent threshold of coupling strength for synchronization using linear stability analysis and numerically verify the obtained analytical results. Additionally, we explore the concept of basin stability, a nonlinear measure based on volumes of basin of attractions, to investigate how stable the synchronous state is under large perturbations. The synchronization phenomenon is analyzed taking limit cycle and chaotic oscillators for wide ranges of parameters like interaction strength between the walkers, speed of movement and vision range .
pacs:
05.45.Xt, 87.75.FbRecent researches attest to the fact that synchronization in timevarying complex network is very essential due to its possible applications in various fields. Particularly, in fields like mobile communications, robotics etc. synchrony among moving and dynamic agents is certainly necessary. But only a few works on this subject have been reported so far. For some situations in coordinated motions of vehicles, flashing fireflies and swarming birds etc. interaction with any other agent lying in a specified zone may not always be possible. In this article, we try to resolve this issue by forming those zones in a slightly different manner and examine global synchrony in a network of mobile agents moving in a finite arena of three dimensional space, where an oscillator is assigned for each agent. In this context, both limit cycle and chaotic systems have been dealt with. The interaction between the oscillators depends on their proximity and the orientation of the movements as well. We elaborately study the influences of all the network parameters such as vision range, speed of motion and coupling strength both locally (using linear stability analysis) and nonlocally (using basin stability analysis).
I Introduction
In the last few decades, complex network theory has revealed that the type of interaction between the nodes is very important in a network as far as the emergence of collective behavior such as synchronization is concerned. Synchronization refers to a process wherein two (or many) systems adjust a given property of their motion to a common behavior due to interaction between each other or to a forcing syper (); sycha (). Recently, collective behavior of timevarying complex network has been an active research area due to its potential application in power transmission system pts (), consensus problem cons (), persontoperson communication mcom (); mcom2 (); mcom3 (); mcom4 (), to name only few. In timevarying networks, connection topology of the network changes over time and emerging behaviors can be realized in such timevarying complex networks due to that varying interaction topology. Timevarying networks can be based on different types of interaction platforms, e.g., nodes of the network can be static in space with timevarying edges (coupling) between them and on the other hand, nodes can be moving in the space with the formation of the edges depending on the movement of the nodes. Timevarying network of the first type arises in functional brain networks discussed in fbn () where the outcomes confirmed that the brain connectivity patterns develop with time and frequency while preserving a smallworld structure. Also in power transmission system pts (), where a random reconnection of links between nodes has been taken into account, in consensus problem cons () in which robustness to changes in multiagent network topology due to node or link failures is analyzed and this structure has also been investigated in mobile communication mcom (). These makes the analysis of the timevarying networks of the first type so important. Emergence of synchronization on these type of timevarying dynamical networks has been studied theoretically belykh (); sinha (); sychen1 (); sychen2 () and experimentally danas ().
On the other hand, there are so many examples of networks of movable physical agents that can interact with each other depending strictly on their proximity. For example, studies on group of animals, particularly where transition from disorder to order in mobile swarming desert locusts is reported in grani (). Coordinated motions of robots robot () where interactions within a group of mobile robots are investigated and proposed as a possible strategy for controlling the robots and also coordinated motions in vehicles vehicle () can be considered as such networks of movable physical entities in which regulated dynamics appear. Notably, there are also many examples where synchronization in moving oscillators play an essential role, e.g., in chemotaxis discussed in chemo () where the motion of each element is driven by the local gradient of the chemical density and the elements produce and consume this chemical in amounts that depend on its internal state. In wireless sensor networks sensor (), sensor nodes need to coordinate their operations and collaborate to achieve a complex sensing task, and also in mobile ad hoc networks adhoc ().
Therefore, systematic study on emerging behaviors in networks of moving oscillators deserve special attention. In this context, some early works include studying synchronization in networks of movable Kuramoto phase oscillators movkur () and also in mobile chaotic systems movch () where topology changes due to agent’s motion. Spectral analysis of synchronization in such networks has also been done spect () and it has been shown that the spectral pattern differs between synchronization mechanisms. A dynamical network of moving oscillator is investigated to model the spread of infectious disease infect (). Such a network is also considered to show the effectiveness of nonlinear coupling over linear one in the context of synchrony of the network of moving oscillators punit (). Emergence of synchrony in a network of moving integrateandfire oscillators interacting with immediate neighbors upon firing time has also been analyzed if1 (); if2 (). Synchronization properties in ensemble of random walkers carrying internal clock in a periodic space under normal diffusive motion together with superdiffusion has been studied comprehensively rev1 (). Influence of selfpropulsion on the emergence of synchronization in dynamical networks of mobile oscillators is investigated rev2 () recently. Pattern formation in chemically interacting active rotors with selfpropulsion rev3 () needs to be addressed as well.
Other notable works on synchronizability of a moving agent network includes wang (); pori (); chaos2016 (). B. Kim et al.kim () investigated synchronous behavior of movingagent network in three dimensional (3D) space based on the assumption of restricted interactions, where a certain number of fixed zones were preitemized and at any time, only those agents within these zones were supposed to interact with each other. This scheme may be related to the process of coordination in robotic networks.
But, most of these works discussed above have been done by considering the movement of the agents solely in a finite region of two dimensional plane. The movement of the agents in (3D) space could be more general platform of moving which is practically possible in many ways. In particular, the scenario of moving oscillators in 3D space has a direct application in biological systems such as flashing fireflies which are moving in 3D space and interchange the light signal to search the potential mates firefl1 (); firefl2 (). More importantly, in the cases of coordinated motions in vehicles, flashing fireflies etc., for a particular agent, interaction with other agents lying inside a specified region like sphere or circle (in case of 3D or 2D spaces respectively) centered at that agent (proposed in most of the previous works), may not be possible. Rather, the connection may be formed depending on the direction of the movement and the agent may interact with only those, the agent can see (observe) and not with any others, staying on its back. For all these situations our algorithm (discussed in the next section) will be significant and crucial.
So, in this paper, we put the scheme of mobility in a general way and propose a finite (cubical) region of three dimensional space and study networks of limit cycle as well as chaotic oscillators. We discuss one fundamental emerging behavior in complex network, namely, synchronization. We numerically investigate the synchronization process for a class of moving oscillators including chaotic systems which are moving in 3D space. Most of the earlier works movkur (); movch (); spect (); infect (); punit (); if1 (); if2 (); wang () reported the synchronization behaviors in moving oscillators which are moving in two dimensional plane but little attention has been paid on synchronization in timevarying network where each oscillator is moving in 3D space. Here, in our scheme, the vision range decides the number of oscillators with which a particular oscillator interacts at an instant of time. We rigorously study the variation of order parameter for synchronization by changing three physical parameters, namely speed of movement of the agents , vision range and interaction strength . We also derive the density dependent threshold of critical coupling strength for synchronization in the network using linear stability analysis. To explore the stability of the synchronous state against nonsmall perturbations, we use basin stability measure, which is a nonlocal and nonlinear measure of stability depending on the basin of attraction’s volume.
The remaining part of this paper is organized as follows. In Sec. II, we discuss how the oscillators move in three dimensional space. In Sec. III, we present the mathematical form of the network of moving oscillators where network topology changes due to oscillator’s motion. In Sec. IV, numerical results are shown for the moving oscillators network using limit cycle (LandauStuart) and chaotic (Rössler) oscillators and we show how the level of synchrony depends upon the parameters like vision range, speed of movement and interaction strength. Linear stability analysis of the synchronized state is discussed in Sec. V. In Sec. VI, we examine the stability of the synchronized state using basin stability approach. Finally, we summarize our results in Sec. VII.
Ii Random walk algorithm
In this section, we demonstrate the scheme of random walk movement of oscillators in a finite region of three dimensional space. Initially, in physical space, number of oscillators are randomly distributed in node mesh with a pseudo random number ( ) corresponding to the th oscillator having position coordinate , . The algorithm makes all the oscillators to move either to right or left, to up or down, to front or back depending on the values of the pseudo random number associated with the oscillators. Two parameters are associated with the random walk in physical space, namely, the speed of movement of the oscillators and the vision range .
The random walk movement of the oscillators in space is made by a simple change in the position coordinates after each and every time iteration in the following way:

If the pseudorandom number corresponding to the th oscillator lies in or , then the movements along the +ve axis (the right) or ve axis (the left) direction are given as :
(1) corresponding to or respectively.
Here corresponds to the characteristic time for the movement of the agents. Similarly, the movements along axes (to front or back) and axes (to up or down) corresponding to the cases of pseudorandom number respectively lying in or and or is defined. 
Updating time iteration to its next value, the described changes will be applied again for all the position coordinates accordingly and so on. For this, of course, varies with iteration.

Now according to the proposed algorithm without any suitable boundary conditions, the oscillators could be out of the physical space, hence we impose periodic boundary conditions. Doing this, we make sure that all the oscillators remain in the physical space for all the course of time.
As already stated earlier, no oscillator interacts with all the other oscillators, rather they interact with only those that belong to a specified region (subregion) of that particular oscillator. That means oscillators will interact with each other only when they are in a specified close proximity for a particular instant of time. Here in order to make it possible, inside the physical space in which oscillators are allowed to move, a much smaller region (subregion) of cubical structure is assigned to each and every oscillator. We call the subregion for an oscillator as vision size.
Suppose at a particular instant of time , an oscillator is moving to its right (say) which is the case of the corresponding pseudorandom number to be in , then a cube (the vision size) of volume is created to its right. Then the oscillator will interact only with those oscillators which lie in the created cube. But it does not interact with any oscillator to its left at that instant of time. Similarly, if an oscillator is moving up, which is the case of the corresponding pseudo random number to be in , then a cube of volume is created in upward direction. Then the oscillator will interact with the oscillators which belong to its vision size but it will not interact with any oscillator in its downward direction. Actually, the vision size (the cube) for a particular oscillator will be created in a direction along which the oscillator is moving. Creation of this type of vision size in the direction of motions is highly realistic in the cases of flashing fireflies, coordinated motions in vehicles etc. The same algorithm is applied to all the remaining cases. A graphical view of movement of oscillators and creation of vision size at a particular instant is shown in Fig. 1.
Iii Mathematical model of moving oscillators
Under the above assumptions, we consider number of mobile agents or walkers which are moving in 3D space. Each agent is a dynamical system whose state at time is described by , where is a dimensional vector of dynamical variables of the th oscillator and is a vector field. The interaction between the random walkers changes over time and the mathematical form of network of moving oscillators is
(2) 
where is the control parameter accounting for the strength of the interaction between the random walkers. Here is a timevarying zerorow sum Laplacian matrix of order modeling the network connections at time . In particular, if th oscillator lies in the vision size of th oscillator and otherwise zero. Also, where is the number of oscillators which lie in the vision size of th oscillator at time . The matrix is the coupling matrix indicating which variables of one oscillator are coupled with which variables of the others. For example, coupling with only the first components, the form of is given by
Finally, the network of moving oscillators has three control parameters: the interaction strength between the random walkers , the speed of movement of the agents and the vision range of the oscillators. We fixed the parameters and mainly, throughout this paper. We integrated the network equation (2) using fifthorder RungeKuttaFehlberg algorithm with a time step size of . The pseudorandom number is generated in FORTRAN 90. Initially, at time , all the oscillators are randomly distributed in the physical space using pseudorandom number In the following, we will explore the dynamics of the moving limit cycle and chaotic oscillators using the timevarying network (2), our main emphasis will be to identify the parameter range for three parameters and for synchronization.
Iv Results
We first demonstrate how perfect synchrony emerges due to random walk movement in three dimensional space using coupled network (2). To do this, we take LandauStuart oscillator, when uncoupled, exhibits a stable limit cycle near supercritical Hopf bifurcation and has an unstable focus at origin. In the last two decades, LandauStuart oscillator has been used in the context of many emerging behaviors. So it will be interesting to explore our scheme using coupled LandauStuart oscillators. We consider moving LandauStuart oscillators whose topologies change with the oscillators’ motion in the form (2) with
(3) 
where , the coupling matrix is , and is the intrinsic frequency of the individual LandauStuart oscillator. Only the components of LandauStuart oscillators are connected with each other via diffusive coupling.
In order to quantify the level of synchrony, we calculate the order parameter as,
(4) 
where , , the phase is calculated in plane as , and represents average over time. For complete synchronization, attains its value near unity while when the oscillators are desynchronized.
In Fig. 2, dependence of the order parameter initc () on different system parameters, namely interaction strength , speed of movement and vision range are shown. Figure 2(a) shows how order parameter varies with respect to interaction strength where we fix other two parameters at and . We observe that starting with a small value, starts to increase with increasing and finally gets saturated after reaching , which is the maximum value it can achieve that indicates complete synchrony between the moving agents. Here the order parameter starts saturating after a certain value which means is sufficient to get perfect synchrony. Next we fix the interaction strength and speed of movement and vary the vision range of the agents. The variation of the order parameter with respect to vision range of the agents are shown in Fig. 2(b). As the vision range increases, we note that starting again with a small value, gets saturated after reaching , which indicates complete synchrony among the moving oscillators. Thus as increases, the volume of the subregion and hence the possible number of oscillators in the subregion also increases, which makes the synchronization faster. In this case, the order parameter starts saturating after a certain value which means would be enough to get absolute synchrony. Finally, Fig. 2(c) shows change in the order parameter in terms of speed of movement for and . As the speed increases, a type of long distance jump of the nodes occurs and we see that starting again with the small value, starts developing with increasing and saturation arises after reaching indicating complete synchronization between the moving oscillators. In this case the order parameter for .
In order to understand the different synchronization regimes, we need to compare the timeseries of moving agents for different physical system parameter values that lead the system to its perfect synchronization state. For this, we take three different values of from Fig. 2(a) where the other two parameters are fixed at and Figures 3(a) and 3(b) show the time series and spatiotemporal plot respectively for components of moving LandauStuart oscillators for . From these figures, it is noted that all the moving oscillators are in incoherent state (as expected, since the value of is much lower than unity). Next we consider the higher value of , the time series and spatiotemporal plot of the oscillators are shown in Figs. 3(c) and 3(d) which clearly show that all the oscillators are not in complete synchronization state. But a coherence in the time series profiles can be seen for higher value of interaction strength . At this point, the value of is almost unity and complete synchrony occurs. Finally, Figs. 3(e) and 3(f) respectively illustrate the time series and spatiotemporal plot for complete synchrony state.
Next we show that perfect synchrony appears due to random walk movement in space not only for oscillators of limit cycle type but also for oscillators that are chaotic in nature by taking moving Rössler oscillators. We consider moving Rössler oscillators in the form (2) where has the form:
(5) 
with . Each oscillator is in chaotic state for , . Only the components of the Rössler oscillators are connected with each other, the coupling matrix is in the form .
Now we analyze how the other physical parameters of the system, namely interaction strength , speed of movement and vision range influence synchronization similar to the previous case. In Fig. 4, variation of the order parameter for synchronization of moving Rössler oscillators by varying for different are depicted. Figure 4(a) shows how order parameter varies with respect to interaction strength with fixed values of and for . Here the order parameter starts saturating with for a certain value . Figure 4(b) shows how order parameter varies with respect to for with , fixed as above. Similarly as before, the order parameter starts developing with increasing and finally gets saturated for (where ) which means would be enough to get absolute synchrony. Figure 5 characterizes the dependence of the order parameter on and with fixed. Figure 5(a) shows variation of order parameter with respect to vision range where the values of and are fixed. From this figure we see that for small values of , all the moving oscillators are in desynchronized state as there are almost no links between the oscillators. As we increase this parameter , we can see that synchronization of the moving oscillators appear for . Finally Fig. 5(b) shows change of order parameter in terms of speed of movement where and .
As already discussed, a saturation in the value of the order parameter can be seen no matter which physical system parameter, namely , or is varied. So it is important to look at the time series profile to judge upto what extent it changes depending on the physical system parameters before and after gets saturated. To see how much the time series profile changes for before and after its saturation, we plot time series and spatiotemporal evolution for components of moving Rössler oscillators for three different values of interaction strength at fixed values of and with in Fig. 6. At lower value of , incoherent states occur among the oscillators as in Fig. 6(a) and corresponding spacetime plot is shown in Fig. 6(b). Figures 6(c) and 6(d) illustrate the timeseries and spacetime plot for , that show some of the oscillators are near to inphase synchrony. By increasing the interaction strength to complete synchrony occurs as the value of order parameter is there. Figures 6(e) and 6(f) describe the time series and spatiotemporal plot respectively for , which show that all the moving oscillators achieve complete synchrony.
V Linear Stability Analysis
In this section, we analytically calculate the critical coupling strength for stable synchrony of moving oscillators. A time average Laplacian matrix is computed for the timevarying network using the procedure proposed by Stilwell et al. stilwell (). According to stilwell (), if the system of coupled oscillators given by
(6) 
with fixed interaction topology (i.e., with fixed Laplacian matrix ) admits a stable synchronization manifold and if there is a constant such that , then there exists such that for all fixed the system defined by
(7) 
(i.e., coupled according to timevarying Laplacian ) also carries a stable synchronization manifold. Thus, if the time average of the Laplacian matrix , given by reflects synchrony of the system, then the time varying network will also support synchronization (provided that the switchings between the network configurations are fast).
Here we note that in this work, the time scale of the agents’ motion is considered much shorter than that of the oscillator dynamics. To find the time average Laplacian matrix , we start with the simplest case of for which four network configurations are possible at each time t. These network realizations are based on the following:

There is no interaction between the two nodes.

There is a unidirectional interaction from the node to the node .

There is a unidirectional interaction towards node from node .

There exists bidirectional interaction between the nodes.
Then the timeaverage matrix is given by
(8) 
where is the probability that there is no interaction between the nodes and , , are the probabilities corresponding to the cases of unidirectional interaction towards node , towards node and bidirectional interaction among the nodes respectively. , , and are the corresponding Laplacian matrices:
and
As (say), becomes . Now, being the probability of activation of a link between the nodes, , where ( is the volume of the considered physical space), under static node assumption and for moving node consideration (small fluctuation may arise depending on ). Henceforth we will assume and so .
In fact, it can be shown that the similar result holds for any number of nodes and with and is the global (alltoall) Laplacian matrix with bidirectional interactions. In fact, the above calculations can be repeated for any (Confirmation using can be found in Appendix A) and finally we can conclude that
(9) 
where
As reported in msf1 (); msf2 (); msf3 (); msf4 (); msf5 (), a block diagonalized variational equation of the form represents the dynamics of the system around the synchronization manifold whenever the dynamics of the network is described by the equations
(10) 
where is the coupling function, is the th eigen value of , and , are the Jacobian matrices corresponding to and around the state of synchrony. In terms of complex eigenvalues (), the master stability equation becomes . Then the maximum Lyapunov exponent of this equation serves the master stability function (MSF) as function of and . Thus the stbility of the synchronized state can be investigated from the eigenvalues of and looking at the negative sign of at . Now considering type MSF (for which the synchronization manifold is stable in the interval of type ), as in our case, eigenvalues of are , , there is a critical value of , say that satisfies
(11) 
which implies
(12) 
with (for coupled Rössler systems).
Clearly, according to the considered values of and , for , and for , , which are also the values of that we found for synchrony in Figs. 4(a) and 4(b) respectively. In fact, decreases with increasing according to the relation (12). For complete understanding, the critical stability curve above which the network achieves stable synchrony is plotted in the parameter space in Fig. 7. Thus, the critical value for achieving synchronization decreases for increasing values of (i.e., the density of moving oscillators), which is quite expected as, if the number of nodes increases in the physical space of fixed volume, the possibility of interaction between nodes increases as well. Furthermore, pursuing the movement algorithm proposed in the article, the nodes move in the space at every iteration step, so that their neighbors (with which they interact) also get changed very frequently. And because of their movement in the whole space, effectively the nodes gain opportunity of interacting with so many (may be even all) other nodes during the entire course of time. Here these facts are inspected to be playing a deciding role in realizing global synchrony at a lower coupling strength as the density increases. Hence for fixed interaction strength , there will be critical value of population density for which synchrony in the network emerges. This densitydependent threshold for achieving synchronization is related to that in quorumsensing transition in indirectly coupled systems quo1 (); quo2 (); quo3 (). This mechanism plays a decisive role in bacterial infection, biofilm formation and bioluminescence quo4 (). However, in our algorithm the nodes are directly interacting with other nodes lying in the vision range.
Vi Basin stability analysis
The sections discussed above are made up of linearizationbased study on synchronous state in network of moving oscillators. In other words, the results revealed that the synchronization state is stable for small perturbations near synchronous state and in some cases small perturbations could be infinitesimal. Such stability condition for synchronization using linear stability analysis is necessary but not sufficient jost (). Recently, Menck et al. bs () proposed a nonlocal stability concept to quantify how stable a state is against nonsmall perturbations in terms of basin stability (BS) which reflects the fact that both linear stability and BS approaches should be dealt with in order to check out the stability of the synchronized state. So here, we are concerned with the basin stability analysis of the synchronized and desynchronized state of the moving oscillators’ network. BS of the synchronization state gives information about the relative volume of basin of attraction that approach to the synchronous state. In contrast to linear stability, basin stability is a nonlinear as well as a nonlocal measure that determines how stable a state is while facing nonsmall perturbations. In this case, basin’s volume quantifies the state’s stability.
To calculate the basin stability for Rössler system, we use the following procedure: We randomly choose sufficiently large number of initial conditions within the volume . For each fixed value of network parameter (e.g. coupling strength ), we integrate the network equation (2) together with (5) for number of random initial conditions from the prescribed volume and identify number of initial conditions for which tol (). Then the basin stability at that network parameter is
(13) 
The value of changes in the range from 0 to 1. means that the network exhibits stable synchronization state for large number of initial conditions.
We start with fixed values of and while vary interaction strength for moving Rössler oscillators’ network. As can be noticed in region I=, where is the critical coupling strength obtained analytically in Eq.(12), BS of the synchronized state (black colored curve) remains zero whereas the BS of the desynchronized state (in color magenta) remains unity up to (Fig. 8(a)). In this region, the moving oscillators are in desynchronized state whatever be the amount of perturbation from the synchronization manifold. With an increase in the value of beyond BS of synchrony starts developing in the region II= where kd () is the threshold coupling strength for the unit BS. This indicates that although linear stability predicts as the critical value of after which the synchrony of the network is stable, this is no more the case when nonsmall perturbations (random initial conditions) are taken into account. BS of synchronization (desynchronization) becomes unity (zero) in the region III=, which means synchrony of the network is stable here for even numerous choice of initial conditions from the basin volume and hence confirms stability in a global sense.
Next keeping and fixed, a similar transition from BS zero (unity) to BS unity (zero) for synchrony (desynchronization) is observed for varying vision range in Fig. 8(b). Here again, after having zero values in a certain range I, BS of synchrony (in color red) gets nonzero, nonunit values in region II and finally gets unit value in the range III. BS of desynchronization (in magenta color) takes the complementary values in the respective regions. Finally, Fig. 8(c) characterizes the variation of BS of synchronization (desynchronization) with respect to speed of movement figured in blue (magenta) colored curve where and are fixed. Here also a similar type of increment (decrement) in the value of BS of synchrony (desynchronization) is observed as before. Similarly the variations of BS for LandauStuart oscillators by changing the physical system parameters are discussed in Appendix B.
Vii Conclusion
In summary, we have presented a network model of interaction between moving oscillators where each oscillator is moving in three dimensional space. Depending upon the network parameters, one of the most intriguing emergent phenomenon in coupled oscillators system, namely synchronization for networks of spatially moving limit cycle and chaotic oscillators in three dimensional space is observed. We have investigated the parameter regions for synchronization state by varying the three physical systems’ parameters, namely the interaction strength , speed of motion of the moving oscillators and vision range . Density dependent threshold of coupling strength for synchrony using linear stability analysis has also been determined. Apart from the numerical investigation and linear stability analysis, we probe the notion of basin stability, a nonlinear measure based on volumes of basin of attractions, to investigate how stable the synchronous state is under arbitrary perturbations.
The present study is important in the context of timevarying complex networks where the interactions between elements (species, individuals etc.) that occur in real systems (biological, social etc.) changes with respect to time due to spatial movement of the elements. Since we have considered the collective behaviors of moving oscillators in threedimensional space, obtained results has a direct application in biological systems such as flashing fireflies firefl1 (); firefl2 () which are moving in 3D space.
APPENDIX A:
CONFIRMATION OF RESULTING FORMULA (9) WITH
To confirm the concluding formula given in (9), here we consider . In this case, the possible network configurations are based on the following cases:

There is no interaction between the nodes.

Two nodes are interacting with each other but the third one is not interacting with any of these two nodes.

A single node is interacting with the other two nodes but there is no connection between these two nodes.

There exists connection between all the nodes.
Then the timeaverage matrix is given by
(14) 
where is the probability that there is no connection between any node. is the probability that node and node are connected with each other but node is isolated. and are similarly defined. is the probability of the case when node is interacting with both node and node but node and node are not connected. and have similar meanings. Finally, stands for the probability that all the three nodes are interacting with each other. , , , , , , and are the corresponding Laplacian matrices.
Combining all the types of interactions (unidirectional and bidirectional), can be written as
(15) 
where the probabilities (Laplacian matrices) () or () corresponds to the cases when there is a unidirectional connection from node to node or from node to node and () is the probability (Laplacian matrix) that node and node are connected bidirectionally.
Since, (say) and
and
So, from Eq.(15) we obtain Similarly,
and
Again since (say),
, where
In a similar way, taking all the types of interactions into account becomes
(16) 
where the probabilities (Laplacian matrices) () or () are for the cases when there is a unidirectional connection from node to both node and or from both node and to node . () or () corresponds to the cases when there is a unidirectional connection from node to node and from node to node or from node to node and from node to node and finally () is the probability (Laplacian matrix) that node connects to both node and bidirectionally. But, node and are not connected in these cases.
Since, (say) and
and
So, Eq.(16) gives with
Proceeding in the same way, we have
with
and
with
Again since (say),
(say),
(say),
and (say) so
In this way, taking all possible type of interactions among all the three nodes, we will have , where is the sum of all the probabilities corresponding to different type of interactions (considering every possible combination of unidirectional and bidirectional interactions). Thus from Eq.(14), becomes,
(say), where is the probability of activation of a link between three nodes and
Thus where and is the zerorow sum Laplacian matrix of order three for global bidirectional interactions.
APPENDIX B:
BASIN STABILITY FOR MOVING LANDAUSTUART OSCILLATORS
Figure 9(a) shows how the BS of the synchronized (desynchronized) state changes with interaction strength where and are fixed for moving LandauStuart oscillators’ network. In region I of Fig. 9(a), BS of the synchronized state (black colored curve) remains zero whereas the BS of the desynchronized state (in color magenta) retains maximum possible value (the unity) upto . In region II, BS of synchrony starts increasing as passes this value, but here this BS is nonunit as well. BS of synchronization (desynchronization) turns to be unity (zero) in the region III. Next we keep and fixed, for which a similar transition from BS zero (unity) to BS unity (zero) for synchrony (desynchronization) is observed for increasing vision range as in Fig. 9(b). Here red and magenta colors are used to figure curves representing BS of synchronization and desynchronization respectively. Finally, Fig. 9(c) depicts the change in BS of synchronization (desynchronization) with respect to speed of movement where and are kept fixed.
Acknowledgments: Authors thank the anonymous referees for valuable suggestions to improve the manuscript in the present form. D.G. was supported by SERBDST (Department of Science and Technology), Government of India (Project no. EMR/2016/001039).
References
 (1) I. Blekman, Synchronization in Science and Technology, Nauka, Moscow, 1981 (in Russian); ASME Press, New York, 1988 (in English)
 (2) S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, C. S. Zhou, Phys. Rep. 366, 1101 (2002).
 (3) M. L. Sachtjen, B. A. Carreras, and V. E. Lynch, Phys. Rev. E 61, 4877 (2000).
 (4) R. OlfatiSaber, J. A. Fax, and R. M. Murray, Proc. IEEE 95, 215 (2007).
 (5) J.P. Onnela, J. Saramaki, J. Hyvonen, G. Szabo, D. Lazer, K. Kaski, J. Kertesz, and A. L. Barabasi, Proc. Natl. Acad. Sci. USA 104, 7332 (2007).
 (6) Y. Wu, C. Zhou, J. Xiao, J. Kurths, and H. J. Schellnhuber, Proc. Natl. Acad. Sci. U. S. A. 107, 18803 (2010).
 (7) J. L. Iribarren and E. Moro, Phys. Rev. Lett. 103, 038702 (2009).
 (8) R. Guimerà, L. Danon, A. DíazGuilera, F. Giralt, and A. Arenas, Phys. Rev. E 68, 065103 (2003).
 (9) M. Valencia, J. Martinerie, S. Dupont, and M. Chavez, Phys. Rev. E 77, 050905(R) (2008).
 (10) I. V. Belykh, V. N. Belykh, and M. Hasler, Physica D 195, 188 (2004).
 (11) V. Kohar, P. Ji, A. Choudhary, S. Sinha and J. Kurths, Phys. Rev. E 90, 022812 (2014).
 (12) J. Lü, G. Chen, IEEE Trans. Autom. Control 50, 841 (2005).
 (13) M. Chen, Phys. Rev. E 76, 016104 (2007).
 (14) S. K. Bhowmick, R. E. Amritkar, and S. K. Dana, Chaos 22, 023105 (2012).
 (15) J. Buhl, D. J. T. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson, Science 312, 1402 (2006).
 (16) A. Buscarino, L. Fortuna, M. Frasca, and A. Rizzo, Chaos 16, 015116 (2006).
 (17) H. G. Tanner, A. Jadbabaie, and G. J. Pappas, in Proceedings, 42nd IEEE Conference on Decision and Control, pp. 20162021 (2003).
 (18) D. Tanaka, Phys. Rev. Lett. 99, 134103 (2007).
 (19) F. Sivrikaya, and B. Yener, IEEE Network 18, 45 (2004).
 (20) K. Römer, in Proceedings of the 2nd ACM International Symposium on Mobile ad hoc Networking & Computing (ACM, NY, USA, 2001), pp. 173182.
 (21) N. Fujiwara, J. Kurths, A. DÂ´Ä±azGuilera, Phys. Rev. E 83, 025101(R) (2011).
 (22) M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna and S. Boccaletti, Phys. Rev. Lett. 100, 044102 (2008).
 (23) N. Fujiwara, J. Kurths, A. DíazGuilera, AIP Conf. Proc. 1389, 1015 (2011).
 (24) M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna and S. Boccaletti, Phys. Rev. E 74, 036110 (2006).
 (25) L. Janagal, P. Parmananda, Phys. Rev. E 86, 056213 (2012).
 (26) L. Prignano, O. Sagarra, A. DíazGuilera, Phys. Rev. Lett. 110, 114101 (2013).
 (27) L. Prignano, O. Sagarra, P. M. Gleiser, A. DíazGuilera, Int. J. Bifurcat. Chaos 22, 1250179 (2012).
 (28) R. Großmann, F. Peruani, and M. Bär, Phys. Rev. E 93, 040102(R) (2016).
 (29) D. Levis, I. Pagonabarraga, and A. DíazGuilera, Phys. Rev. X 7, 011028 (2017).
 (30) B. Liebchen, M. E. Cates and D. Marenduzzo, Soft Matter 12, 7259 (2016).
 (31) L. Wang, H. Shi, and YX. Sun, Phys. Rev. E 82, 046222 (2010).
 (32) M. Porifiri, D. J. Stilwell, E. M. Bollt and J. D. Skufca, Physica D 224, 102 (2006).
 (33) N. Fujiwara, J. Kurths, and A. DíazGuilera, Chaos 26, 094824 (2016).
 (34) B. Kim, Y. Do, and YC. Lai, Phys. Rev. E 88, 042818 (2013).
 (35) R. E. Mirollo and S. H. Strogatz, SIAM (Soc. Ind. Appl. Math.) J. Appl. Math. 50, 1645 (1990).
 (36) G. M. RamírezAvila, J.L. Deneubourg, J.L. Guisset, N. Wessel and J. Kurths, Euro. Phys. Lett. 94, 60007 (2011).
 (37) While calculating all the order parameters (throughout the paper, in Figs. 2 , 4 and 5), we have taken the initial conditions very near to the synchronization manifold so that the order parameter reaches 1 (the unity) exactly when i.e., at the earliest possible.
 (38) D. J. Stilwell, E. M. Bollt, and D. G. Roberson, SIAM J. Appl. Dyn. Syst. 5, 140 (2006).
 (39) L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
 (40) M. Barahona and L. M. Pecora, Phys. Rev. Lett. 89, 054101 (2002).
 (41) L. Huang, Q. Chen, YC. Lai, and L. M. Pecora, Phys. Rev. E. 80, 036204 (2009).
 (42) J. Sun, E. M. Bollt and T. Nishikawa, Euro. Phys. Lett. 85, 60011 (2009).
 (43) A. E. Hramov, A. E. Khramova, A. A. Koronovskii, and S. Boccaletti, Int. J. Bifurcat. Chaos 18, 845 (2008).
 (44) A. Camili and B. L. Bassler, Science 311, 1113 (2006).
 (45) A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang, K. Showalter, Science 323, 614 (2009).
 (46) S. Majhi, M. Perc, and D. Ghosh, Sci. Rep. 6, 39033 (2016).
 (47) C. D. Nadell, J. Xavier, S. A. Levin, and K. R. Foster, PLoS Comput. Biol. 6, e14 (2008).
 (48) J. Jost and M. P. Joy, Phys. Rev. E 65, 016201 (2001).
 (49) P. J. Menck, J. Heitzig, N. Marwan and J. Kurths, Nature Phys. 9, 89 (2013).
 (50) In this context, satisfying has been taken into consideration that corresponds to synchronizability of the network.
 (51) The value of may vary depending on the choice of .