Synchronization in leader-follower switching dynamics

The features of animal population dynamics, for instance, flocking and migration, are often synchronized for survival under large-scale climate change or perceived threats. These coherent phenomena have been explained using synchronization models. However, such models do not take into account asynchronous and adaptive updating of an individual's status at each time. Here, we modify the Kuramoto model slightly by classifying oscillators as leaders or followers, according to their angular velocity at each time, where individuals interact asymmetrically according to their leader/follower status. As the angular velocities of the oscillators are updated, the leader and follower status may also be reassigned. Owing to this adaptive dynamics, oscillators may cooperate by taking turns acting as a leader or follower. This may result in intriguing patterns of synchronization transitions, including hybrid phase transitions, and produce the leader-follower switching pattern observed in bird migration patterns.


INTRODUCTION
An intricate leader-follower dynamics associated with the rankings of the agents is often observed in diverse areas, including sports, such as team pursuit in speed skating [1] and track cycling, and biological systems, such as bird flocks [2,4,5], insect swarms [6,7], or fish schools [8].In team pursuits in speed skating, a team of three skaters races with another team with the goal of overtaking the other team, which starts on the opposite side of the rink.Team pursuit is technically demanding; skaters in a team need to follow each other closely and synchronously in a line to minimize the drag.The leader faces the wind, whereas the followers encounter less wind.To reduce the energy expended, the leader leaves the front of the team and rejoins the team at the rear.This switching behavior helps minimize the total energy cost of the team.Otherwise, the leader will be exhausted quickly.Bird flocks also exhibit similar switching dynamics during seasonal long-distance migrations.It was recently found that juvenile Northern Bald Ibis cooperate by taking turns as leaders and precisely matching their flying times in the trailing and leading positions.This strategy of shuffling positions enables the flock to withstand the physical exertion during migration flights.
Attempts have been made to understand the movement of groups of animals such as fish schools, bird flocks, and insect swarms in terms of statistical mechanics using maximum entropy methods [9], agent-based simulations [10], and analytically tractable models [11], including synchronization models.The Vicsek model [2] provides useful insights into collective swarming behavior in nature, for example, in bird flocks [4,8].Agents can swarm by adapting the average motion of nearby agents, which may be interpreted as a synchronization phenomenon.Moreover, the presence of an actual leader agent can generally enhance the cohesive movements of the Vicsek agents [15,16].The Vicsek model shares some features with the Kuramoto model (KM) in the aspect that, as the interaction strength increases, a phase transition occurs from the incoherent state to the coherent state, i.e., a collective synchronous state.By contrast, whereas the Vicsek model can exhibit rich spatiotemporal patterns in two [2,7] or three dimensions [3], the KM is more suitable for complex dynamics on a ring.Nevertheless, the two models play complementary roles in the study of collective behavior in natural systems.However, to the best of our knowledge, these synchronization models do not take into account asynchronous updating of an individual's status at each time.

MODEL
Inspired by the leader-follower switching dynamics described before, here we modify the original KM with leader-follower switching dynamics.In the model, all oscillators are classified into two groups, the leader group L and follower group F according to their angular velocity rankings at each time step.Each of followers interact with all the other oscillators, whereas each of leaders interact with only followers.According to this rule, circular motions of the leaders are suppressed, in the sense that the leaders, unlike the followers, do not have the advantage of being dragged by others.Thus, this model is called the restricted KM (r-KM).We calculate the angular velocity of each oscillator at each instance as follows: where θ i (t) and θi (t) denote the angle and angular velocity of oscillator i, respectively.The dot represents the time derivative.The intrinsic frequencies {ω i } of each oscillators follow the Lorentzian distribution g(ω) = γ/[π(γ 2 + ω 2 )] with γ = 0.5.K is the interaction strength.N is the system size.hN (0 ≤ h ≤ 1) is the population of followers.We repeat the update dynamics in (1) until the system reaches a steady state.As a result of this updating, the leaders drag the followers, and the followers are dragged by the leaders.The detailed rule is presented in Fig. 1.
It should be noted that switching between leaders and followers occurs.If a follower's angular velocity exceeds that of the leader, the follower and leader switch groups because the number of elements in each group is fixed.Thus, even though the leader and follower groups are separate, their microscopic compositions are continuously updated by switching oscillators according to their angular velocity rankings.We find that a synchronization transition (ST) from an incoherent state to a coherent state occurs at a transition point K c .Multiple transitions and hybrid STs occur depending on the control parameter h.The oscillators with intermediate values of the intrinsic frequency exhibit rapid switching behavior between the two groups.By contrast, those with extreme values of the intrinsic frequency rarely switch unless the coupling constant is sufficiently large beyond the transition point K c .The combined switching intervals span a broad range of values and have a power-law distribution.
The synchronization transition is characterized by complex order parameters defined as  of the oscillators, followers, and leaders, respectively.In addition, In the coherent phase (R = 0), either follower group or both leader and follower groups are synchronized.It is notable that the two groups do not necessarily entrain at a same mean angular velocity in a coherent state.The leader-follower entrainment may follow at a delayed transition point, where the collective phases rotate at a common mean angular velocity Ω = Ω F = Ω L , but with different phase separations given by ψ, ψ L , and ψ F .

NUMERICAL RESULTS
Here, the overall numerical simulation results are presented.The time-averaged coherence R is an order parameter that determines the ST from the incoherent state to the coherent state in Fig. 2(a), (d), and (g).It is observed that the STs may be categorized into three cases depending on the magnitude of the control parameter h: (i) h < h , (ii) h < h ≤ h u , and (iii) h > h u .We demonstrate that for cases (i) and (iii), two STs occur consecutively, whereas for case (ii), a hybrid ST occurs.h ≈ 0.60 and h u ≈ 0.72 are estimated.We first consider case (i), specifically h 0.5.In this case, all of the followers initially have negative natural frequencies ω i , whereas the majority of the leaders have positive ω i .When the coupling constant K is small, switching rarely occurs, and the population arrangement is maintained.Thus, the group angular velocity of the followers, Ω F , is negative, whereas that of the leaders, Ω L , is positive; however, Ω L < |Ω F |, as shown in Fig. 2(b).The two groups rotate in opposite directions.It should be noted that hN followers interact with all the other oscillators, whereas N − hN leaders interact only with the followers.Owing to this rule, the synchronization of followers is reinforced but it is weakly disturbed by the leader.Thus, R F is larger than R L when K is small.In particular, when h is small, the population of followers is smaller than that of leaders, and the disturbance effect is stronger than when h is large.Thus, the transition point K c1 is larger for smaller h.The order parameter begins to increase gradually from zero at K c1 .In the subcritical regime K < K c1 , the average angular velocities |Ω F | and Ω L decrease slowly as K increases to K c1 .K > K c1 is a traveling-wave phase.After K c1 , another transition point exists at K c2 .Within the interval [K c1 , K c2 ], Ω F , and Ω L are entrained, and they then converge rapidly to a stagnant value of the average angular velocity Ω, as shown in Fig. 2(b).At K c2 , a cusp exists in the order parameter R as the system is entrained [Fig.2(a)].As h is increased in (i) but remains below 0.5, the population of followers becomes more similar to that of leaders.Thus, K c1 decreases.However, K c2 increases because |Ω F | becomes more similar to Ω L .The cusp gradually fades as h increases.When h is close to 0.5, the follower and leader populations are balanced, and |Ω F | ≈ Ω L in the subcritical regime.
(ii) When h exceeds 0.5, the population balance breaks down, and |Ω F | becomes smaller than Ω L in the subcritical regime, in contrast to case (i) [Fig.2(e)].Thus, as h is increased, the disturbance strength becomes larger; i.e., the synchronization is more strongly suppressed, and thus, a larger transition point K c is needed.As we observed in the restricted Erdős-Rényi (r-ER) model, this r-KM exhibits a hybrid ST in interval (ii).The order parameter R behaves as R − R 0 ∼ (K − K c ) β , where R 0 is the jump in the order parameter, and K c is the transition point [20,22,27,28].The exponent β is approximately β ≈ 0.75, independent of h, which is in contrary to the behavior in the r-ER model for hybrid percolation transitions induced by cluster merging dynamics [18].However, K c depends on h: K c increases as h is increased.The phase and frequency STs occur at the same transition point, K c .Ω F and Ω L suddenly drop to zero at K c [Fig.2(d)].These behaviors are associated with the pattern of sudden angular-velocity locking, as shown in the inset.At K > K c , the leader-follower angular velocity may vary depending on the parameter values.It is a traveling-wave phase but has small Ω.
(iii) For h > h u , the order parameter R exhibits multiple transition behaviors [Fig.2(e)].In this regime, because h is sufficiently large, the ST of the entire system is governed mainly by the formation of a synchronized cluster of followers.A continuous ST of the follower group occurs at K c1 ; however, a ST of the leader group follows at a different transition point, K c2 , where K c2 > K c1 .This is because the leaders do not interact with each other directly, but gather indirectly with the help of a sufficiently large mass of synchronized followers.Although Ω F drops to zero at K c1 , the leader population runs at a different angular velocity, Ω L = Ω F , until K c3 , where is a chimera phase, and [K c2 , K c3 ] is a standing-wave phase.Finally, the leader group is entrained to the follower group; Ω L suddenly drops to zero at K c3 , and the entire system is frequency-synchronized.The coherent, non-travelingwave phase lies beyond this point.The order parameter R shows a corresponding increase by a small amount, exhibiting a discontinuous transition.A similar consecutive transition followed by a small jump transition has been observed in the percolation transition in multilayer networks when the interaction between layers is weak [36].However, the transition point K c1 of regime (iii) weakly varies on h, which may be attributed to the finite-size effect of the transition point, because this transition point is determined mainly within the follower group and its population depends on h.
Next, we investigate microscopically the leaderfollower switching dynamics.The angular velocity and rank of each oscillator change continuously over time, as shown in Fig. 3(a) and (c), respectively.This figure shows the angular velocity of a given oscillator i in a system of size N = 1000 versus time in the steady state.Note that for the oscillators with intermediate ω ranks, although the changes in the angular velocity, θi , are small in (a), and the fluctuations are also small in (b), the changes in their rankings are quite large in (c), and the fluctuations are also large in (d), particularly at both boundaries of the middle group.Moreover, in Fig. 4(b) and (e), the oscillators with intermediate ω ranks in the shaded region switch their status consecutively with a probability close to unity, spending half of the total time as lead- ers and the other half as followers.Their population increases as K is increased in the supercritical regime.This result seems to resemble the pattern of bird migration, where all the birds act as a leader at roughly the same rate [5,35].By contrast, oscillators on the low or high side of the ω rankings rarely switch their leader-follower status unless the coupling constant K is large enough that the middle class takes over a wide range of population.Between the unchanging high or low class and the switching middle class, there are thin intervals of oscillators that switch frequently over the follower-leader separation with probability q.Let us denote the switch probability distribution as P (q) and assume that it follows a power law, P (q) ∼ q −x .If we assume that the oscillator motions are uncorrelated, the interevent time (IET) distribution is calculated as , where P (z|q) is the conditional probability, and B(2 − x, z) is the beta function.The case x = 0 corresponds to the uniform probability distribution, where the exponent of P (z) is two.A nonuniform distribution or the effects of the correlations may change the precise value of the exponent.Consequently, the IETs between two consec- distribution of leader-follower switches for the r-KM, (middle column) each oscillator's frequency of swings, and (right column) participation rate as a leader for a subcritical coupling K = 0.5 and a supercritical coupling K = 1.5, in regimes (ii) h = 0.6 and (iii) h = 0.8.Note that there are an intermediate swinging class and populations on the sides, which rarely switch.The IET distribution shows crossover behavior between two power-law behaviors.Data in the short IET region are contributed by the thin intervals between the lower and middle and between the middle and upper groups.Data in the long IET region are contributed by the lower group, which also switch, but with small probability.The thin intervals, which build the power laws in the IET distribution, correspond to the kinked regions in the swing and leadership curves.For the subcritical K, all the swing curves share a universal soliton-like pattern, and the IET exponents are also independent of h.The height of the soliton increases rapidly near the transition point.In the supercritical regime, the swing curve is flat-topped [34] and may be accompanied by additional soliton-like features on the sides.The flat-top width tends to grow with increasing K.In regimes (i) and (ii), the swinging population may grow indefinitely until it is nearly equal to the system size, whereas in regime (iii), the maximum width is bounded.In the supercritical regime, the swing curve pattern depends on the detailed parameter values, and supercritical IET exponent values are not universal.
utive switching events are heterogeneous, and their distribution has a power-law tail, P (z) ∼ z −α [Fig.4(a) and (d)].This heavy tail in the IET distribution arises from the dichotomous dynamics and global sorting process.The IET distribution exhibits crossover behavior between short and long IETs.Short IETs are due to the oscillators in the thin intervals between the lower and middle groups and between the upper and middle groups.Long IETs are due to the oscillators in the lower group which switch with small probability.Oscillators in the upper group, on the other hand, do not switch.We find empirically that the exponent of the IET distribution depends on the shape of the frequency distribution g(ω), and also on the model parameters K and h.

DISCUSSION
For the ordinary KM (h → 1), the phase transition type is determined by the shape of the frequency distribution g(ω).A unimodal or bimodal g(ω) yields a continuous or discontinuous transition, respectively [32].A hybrid ST appears in some special marginal cases, such as flat g(ω) classes [27,28,30] and a degree-frequencycorrelated model on a scale-free network with the degree exponent λ = 3 [31].In these cases, the maximum competition arises among the oscillators and suppresses the formation of a synchronization seed cluster.Consequently, a macroscopic number of oscillators are suddenly entrained at a transition point.In contrast to the previous models, the r-KM can produce a hybrid ST if the parameter h is tuned to an appropriate value, even for a unimodal g(ω).The reason is the asymmetric interactions between the leader and follower groups and the tug-of-war critical switching dynamics.Moreover, the r-KM exhibits various other ST patterns, e.g., multiple consecutive transitions.The reason is that the leader and follower groups may be regarded as double layers or communities, where multiple transitions occur naturally when the interaction between the two layers is asymmetric [36,37].By contrast, traveling waves are absent in the ordinary KM with symmetric g(ω).In the r-KM, however, traveling waves are possible owing to the asymmetric interactions between leader and follower groups, which also appear in the KM with asymmetrically competing interactions [28].It is notable that the r-KM takes into account the temporal adaptations of individual oscillators by allowing oscillators to switch their leader/follower status.Hence, the model produces novel self-organized phenomena, where we find a rapidly switching middle class that is half followers and half leaders, as shown in Fig. 4.This type of temporal organization equalizes the distribution of the duration times in the follower and leader positions; this behavior is similar to that observed in the flight patterns of migrating flocks of birds [5].Finally, we remark that Eq. ( 1) is invariant under (θ, ω, g) ↔ (−θ, −ω, 1 − g).Hence, an inverse model of the r-KM, with a limitation applied to the followers instead of the leaders, yields an identical result except that g → 1 − g.In addition, the model has rotational symmetry, and the order parameter curve R(K) is not affected by the translation of the frequency distribution g(ω) → g(ω − Ω 0 ), where the distribution is centered at a nonzero mean, Ω 0 .
In conclusion, we proposed an oscillator model (the r-KM) that includes leader-follower switching dynamics in addition to synchronization.This rank-based status reassignment, which mimics adaptive behavior in complex systems, is unprecedented in synchronization models.This new paradigmatic model exhibits rich dynamical behavior that encompasses a hybrid ST, a traveling-wave phase, standing waves, chimera states, and a powerlaw IET distribution.These diverse patterns of STs arise from the rich behavior of the model, which includes asymmetric interactions, a multilayer leader-follower perspective, and temporal shuffling.Here, we studied the case of Lorentzian distribution g(ω), but the results can be extended to other distributions as well.Methods For numerical integrations we applied the Heun's method with discrete time step size of ∆t = 0.01.Natural frequencies ω i were regularly sampled from the Lorentzian frequency distribution g(ω) = γ/[π(γ 2 + ω 2 )] so that G(ω i ) are equally spaced, where G(ω) = ω g(ω )dω .The width of the distribution was fixed to γ = 0.5, which gives a continuous transition at K c = 1 for the ordinary KM (which corresponds to h = 1 of the r-KM).The traveling wave order parameters Ω, Ω L , and Ω F were obtained by measuring the time averaged rate of change in the complex phases ψ, ψ L , and ψ F .Those collective phases were sampled every two consecutive time steps to get rid of the irrelevant high frequency noises.Order parameters were time averaged for a sufficiently long time after reaching a steady state.

Figure 1 .
Figure1.Schematic illustration of the r-KM with three leaders and four followers.All oscillators are classified into leaders and followers with fraction h and 1 − h, respectively, according to their angular velocity rankings at each time.Each of followers interact with all the other oscillators, whereas leaders interact with only followers.As a result, the leaders drag the followers, and the followers are dragged by the leaders.Hence, a leader-follower status switching may occur if a follower overtakes a leader, as denoted by a pair of arrows.
, andZ L (t) = R L e iΩ L t+ψ L = 1 N − hN i∈L e iθi ,where R, R F , and R L are the magnitudes of phase coherences of all

Figure 2 .
Figure 2. Three cases of ST in r-KM.Plots of the coherence order parameter R (left column) and traveling wave order parameter Ω (right column) versus K. (a),(b) When the number of followers is small (h = 0.3), two consecutive transitions occur at Kc1 and Kc2, from which R is nonzero, and ΩF and ΩL converge to zero, respectively.(c),(d) When the number of followers or leaders has an intermediate value (h = 0.6), the hybrid ST occurs with critical exponent β ≈ 0.66 at Kc, from which R is nonzero, and Ω is zero.(e),(f) When the number of followers is large (h = 0.8), consecutive transitions occur at Kc1, Kc2, and Kc3, from which RF increases, RL increases, and Ω is zero, respectively.In the insets of (b), (d), and (f), the thin gray curves represent the time-averaged angular velocity θi of each oscillator i in steady state versus K.The curve is drawn for every 20 oscillators for visualization.Thick black curve represents Ω(K) for all of the oscillators.Note that at Kc in (d), sudden angular-velocity locking of a macroscopic number of oscillators occurs.By contrast, in (b) and (f), a giant angular velocity cluster grows continuously by merging with other oscillators one by one.For all the panels, g(ω) = γ/[π(ω 2 + γ 2 )], with γ = 0.5, N = 10 3 ; the data for a single configuration are used.

Figure 3 .
Figure 3. Trajectories of (a) angular velocity and (c) ranking as a function of time for given oscillators i, where the index i is ordered by increasing natural frequency ωi.(b) Time average and (d) standard deviation of the quantities in (a) and (c) versus oscillator index i.In (a) and (c), a finite fraction of oscillators in the intermediate rankings have angular velocities with small and fluctuating amplitudes.However, the rankings of those oscillators fluctuate greatly.However, the oscillators with small and large indices i show large fluctuations in angular velocity but small fluctuations in rank.Simulations are performed for N = 10 3 .The control parameters are set to h = 0.7 and K = 1, which are near the transition point for the hybrid ST.In (b) and (d), the index of i is taken in 10oscillator steps.

Figure 4 .
Figure 4. Plots of (left column) the interevent time (IET) distribution of leader-follower switches for the r-KM, (middle column) each oscillator's frequency of swings, and (right column) participation rate as a leader for a subcritical coupling K = 0.5 and a supercritical coupling K = 1.5, in regimes (ii) h = 0.6 and (iii) h = 0.8.Note that there are an intermediate swinging class and populations on the sides, which rarely switch.The IET distribution shows crossover behavior between two power-law behaviors.Data in the short IET region are contributed by the thin intervals between the lower and middle and between the middle and upper groups.Data in the long IET region are contributed by the lower group, which also switch, but with small probability.The thin intervals, which build the power laws in the IET distribution, correspond to the kinked regions in the swing and leadership curves.For the subcritical K, all the swing curves share a universal soliton-like pattern, and the IET exponents are also independent of h.The height of the soliton increases rapidly near the transition point.In the supercritical regime, the swing curve is flat-topped[34] and may be accompanied by additional soliton-like features on the sides.The flat-top width tends to grow with increasing K.In regimes (i) and (ii), the swinging population may grow indefinitely until it is nearly equal to the system size, whereas in regime (iii), the maximum width is bounded.In the supercritical regime, the swing curve pattern depends on the detailed parameter values, and supercritical IET exponent values are not universal.