Transient phases in the Vicsek model of flocking

The Vicsek model of flocking is studied by computer simulation. We confined our studies here to the morphologies and the lifetimes of transient phases. In our simulation, we have identified three distinct transient phases, namely, vortex phase, colliding phase and multi-domain phase. The mapping of Vicsek model to XY model in the $v \to 0$ limit prompted us to explore the possibility of finding any vortex kind of phases in the Vicsek model. We have obtained rotating vortex phase and measured the lifetime of this vortex phase. We have also measured the lifetimes of other two transient phases, i.e., colliding phase and multi-domain phase. We have measured the integrated lifetime ($\tau_t$) of all these transient phases and studied this as function of density ($\rho$) and noise ($\eta$). In the low noise regime, we proposed here a scaling law $\tau_t N^{-a} = F(\rho N^{-b})$ where $F(x)$ is a scaling function like $F(x) \sim x^{-s}$. By the method of data collapse, we have estimated the exponents as $a=-0.110\pm0.010$, $b=0.950\pm0.010$ and $s=1.027\pm0.008$. The integrated lifetime $\tau$ (defined in the text differently) was observed to decrease as the noise approaches the critical noise from below. This behaviour is quite unusual and contrary to the critical slowing down observed in the case of equilibrium phase transitions. We have provided a possible explanation from the time evolution of the distribution of the directions of velocities.


Introduction:
The non-equilibrium phase transition under driven noise, shown by the natural systems possessing collective behaviour, is being studied [1,2,3] with much interest. In nature, we find the collective motion of large groups of motile agents (e.g., collective motion of animal groups [4], flocking of birds, schooling of fish [5], motion of bacteria and microorganisms [6] etc.) which move over a distance much larger than the size of the motile agents. Despite the difference in microscopic interactions of different systems, such collective motions can be described by the common macroscopic behaviours of simple physical models. In 1995, T.
Vicsek and coworkers proposed [7] such a model, which consists of point particles moving with constant absolute velocity but imperfectly aligned in the direction with their neighbours. In the limit of low imperfection in alignment, a phase transition is observed, where all agents move in a particular direction coherently. Needless to say, for higher values of noise, such kind of unidirectional coherent motion is not found. This flocking transition gives rise to symmetry broken ordered phase with long range order even in two dimensions (which was forbidden at equilibrium by Mermin and Wagner theorem [8]). This is the key reason of interest taken by the modern researchers to study [9,10] the Vicsek model of flocking. A considerable amount of investigations were made in recent years which prompted intense interests of active research in Physics [11] even today. We briefly mention a few of those here.
A variant of original Vicsek model was studied [12] with binary interactions among the flocks and a topological distance dependent transition is observed. They [13] have also observed coherent and interesting cyclic states (circular motion of the agents) with topological distances. Very recently, the original Vicsek model is modified [14] by introducing the time delay in the interaction and angular restriction (vision cone). This modification leads to a state where the agents spontaneously condense into a dense drop. The effects of inertia (microorganisms having finite size and mass) of the model flocks in a turbulent environment are also studied [15] recently.
A considerable amount of work has been done to study the nature of the phases and the transitions (continuous or discontinuous) in Vicsek model of flocking. It is reported [16] that the phase transition in the Vicsek model is best understood as a liquid-gas like transition (which shows metastability, hysteresis etc.) which is in contrast to the bulk phase separation observed [17] in active Ising model. From the hydrodynamic equations, three different kinds of solutions are obtained [18], namely, periodic orbits (microphase separation), homoclinic orbits (solitonic objects) and heteroclinic trajectories (phase separation). Although, the scaling was proposed by Vicsek similar to equilibrium critical phenomena, later, it was shown [19] that this transition is indeed first order.
Since the Vicsek model can provide a special non-equilibrium phase transition which breaks the symmetry in ordered phase even in two dimensions, it would be interesting to know whether any vortex kind of phase can be observed in the limit of v → 0 [7]. This may be a common question, what kind of transient phases can be observed, in general. In this paper, we have systematically studied the temporal evolution of the transient phases in the Vicsek model and their lifetimes as function of density and noise. The paper is organised as follows: In section-2 the Vicsek model is described, in section-3 the methodology and the numerical results are reported and the paper ends with a summary in section-4.

The Vicsek model:
The original Vicsek model of flocking is described as follows: identical point-wise particles move continuously on an off-lattice plane represented by a two dimensional square shaped cell of linear size L with periodic boundary conditions. At t=0, N number of particles (i) were randomly distributed in the cell with the position of i th particle denoted by {x i , y i }, (ii) had the same absolute velocity v and (iii) had randomly distributed directions with the direction of i th particle denoted by θ i . The particles move in discrete time steps, i.e., the time interval (∆t) between two updating of directions and positions is considered to be unity (∆t = 1). At each time step, the velocities {vcosθ i , vsinθ i } of the particles are determined simultaneously (by parallel updating of all the particles) and the position of the i th particle is updated according to the relations and the angle is obtained from the expression where θ i (t) r denotes the average direction of the velocities of the particles (including the i th particle) surrounding the given particle within a circle of radius r. The average direction is given by the angle tan −1 [ sin(θ i (t)) r / cos(θ i (t)) r ] and ∆θ is a random number chosen with a uniform probability from the interval [−η/2, η/2] responsible for the noise in the velocity direction. This noise plays the role of indeterminacy in choosing the direction of motion. The noise tends to prevent the ordering in the system. ∆θ is delta-correlated scalar white noise ranging maximally in the interval [-π, π]. The equation (3) shows that each particle tends to self-organize in the same average direction of motion (first term) while the behaviour being randomly perturbed (second term). The flock synchronously evolves through an iteration of these rules and that is studied in detail by determining the absolute value of the average normalised velocity of the entire system of particles. The velocity v a is approximately 0 if the direction of motion of individual particles is distributed randomly (with uniform probability), while for the case of coherent motion (with ordered direction of velocities), v a ∼ = 1. Here, the average normalised velocity is considered as an order parameter of the system. In the steady state, the system undergoes a kinetic phase transition from no transport to finite net transport as the noise (η) and density (ρ = N/L 2 ) of the system are tuned individually. The critical value of noise corresponding to the phase transition is observed in this study by the singular nature (sharply peaked for finite system) of the variance of order parameter in the critical region. This variance (or, standard deviation) has played important role throughout this study in determining several properties of various phases which is to be discussed in the subsequent sections.

Simulations and results:
Inspired by the anticipated analogy of Vicsek model in the v → 0 limit to the wellknown XY-model, we initially speculated possible existence of vortices in the low velocity limit corresponding to analogous Kosterlitz-Thouless transition [20]. However, even while working with v = 0.03 which is used by Vicsek et. al. to study the kinetic phase transition, we have observed not only metastable vortices, but also two other novel transient phases with interesting dynamics. We have estimated their properties and lifetimes individually and have studied the dependence of the integrated lifetime of all these transient phases on density (ρ) and noise (η) of the system. We have explained some novel phenomena obtained in this study with the help of standard statistical tools and techniques.
Throughout all our simulations, the absolute velocity v and the radius of influence r are kept fixed at values v = 0.03 and r = 1.0 respectively. Except for the cases reported in section 3.1, all the characterisations of the transient phases are done using 1000 ensembles of the single realisation.

Morphology and characterisation of different transient phases at low noise:
As reported by Vicsek et. al. and reproduced in this study (Fig. 6a), the said system, when subjected to low noise (η = 0.1) well below the critical noise η c (Fig. 6b), exhibits spontaneous ordering. It starts from completely random configuration at t = 0 (Fig. 1a) and eventually achieves a global dynamically ordered phase exhibiting long range order in the steady state (at t = 1000 in Fig. 1b). The time taken to achieve collective behaviour is an intrigue play of noise, velocity and density. However, during the transient time well before achieving the steady state, the system exhibits different types of phenomena depending upon the density of the system. In the Fig. 1c, for a system with N = 300, L = 9 (ρ = 3.70), we have obtained existence of clear vortex which unbinds after a finite lifetime τ v (vortex region is marked by box in the figure). When both the number of particles and the system size are increased to N = 4000 and L = 33, keeping the density more or less same (ρ = 3.67), the number of vortices is also found to increase (Fig. 1d). We have obtained the second type of transient phase when the density was further lowered to ρ = 2.48 by keeping N fixed and increasing the cell size to L = 11. Fig. 1e depicts that apart from the vortex forming at earlier times, a certain phase exists for a lifetime τ c where two groups of oppositely moving particles collide and change their directions (collision region is marked by box in the figure). When density is further decreased to ρ = 1.78 (N = 300, L = 13), the third type of transient phase is observed where domains of particles having different directions are formed for a finite lifetime τ d at different regions in the box and possessing different strengths (number of particles in a single domain) (Fig. 1f).

Rotating Vortex phase:
In the particular ensemble described in Fig. 1c, the vortex is located roughly in the region marked by the black box. The vortex forms immediately after the system starts updating according to the described rule and lives for around 40 discrete time steps. The vortex melts completely after that and a curvature in the velocity portrait morphology exists for sometime before the entire system attains completely directed order (A video of the simulation can be found in YouTube with the link https://youtu.be/730R0weRacM). The contribution due to vortex is measured by separately studying the designated region. We have calculated the time evolution of the sum of the angles of velocity directions of all the particles participating in the vortex and the average of the same (normalised by proper number of particles falling inside the box) and have compared these to the quantities obtained in the complementary cases of particles excluding the vortex region and the total intact system. Fig. 2a shows that the sum of the angles is almost constant for the entire lifetime of the vortex for the 'only-vortex' case whereas the complementary 'vortex-excluding' part and the total system both show an identical bump in the curve till the vortex exists. This is due to the fact that, when the vortex rolls out with time, the angles of the particles oppositely directed within the circular vortex cancel out each other and make the sum constant over time.
The vortex lifetime τ v is measured by the width of the bump (τ v ≈ 40). The entire time evolution scenario is shown in the inset plot where the 'only-vortex' and 'vortex-excluding' curves show mirror symmetric behaviour to each other, expectedly, whereas the sum achieves a constant value for the total system, since all the particles become uni-directional. The vortex contribution is further verified by calculating the time evolution of normalised average of the angles of directions of the particles (Fig. 2b). The 'only-vortex' contribution and the complementary 'vortex-excluding' contributions show complementary fall and rise in the curves respectively during the transient phase before all of them achieve the same average value denoting collective motion. This deep fall due to the 'only-vortex' contribution reflects that the vortices play non-trivial roles in several important studies. Moreover, the lifetime of the entire transient phase (τ t ) can be found from Fig. 2b where the normalised average of the angles of directions of the entire system of particles attains steady value. Further studies on the integrated lifetime of the transient phase will again be addressed in Section 3.2.
Another method of finding the vortex lifetime is to measure the temporal evolution of the total angular momentum M (Fig. 2c) and the average angular momentum M (Fig. 2d) [21]. Both the Fig.  2c and Fig. 2d show that the angular momentum for the 'only-vortex' case increases and reaches maximum value at the time up to when the vortex lives and then it starts decreasing as soon as the vortex starts melting. From this, the vortex lifetime can be approximately found to be τ v ≈ 40 which is consistent with the previous result. However, for the total system, the angular momentum does increase up to t ≈ 40, but due to the presence of still-existing curvatures in the velocity portrait of the entire system, it falls at a much slower rate than that for the 'only-vortex' case and attains small values when the system eventually achieves steady state ordered phase. This can be explained by the argument that the average angular momentum of the particles tends to achieve higher value when the flock is rotating in a vortex, whereas its values will be less when the particles exhibit unidirectional motion.

Colliding phase:
Although decreasing the density lowers the probability of obtaining vortices, it opens door to new kind of transient phases. One such interesting phenomenon observed in the lower density (ρ = 2.48) regime for the low noise case is collision of particles (A video of the simulation can be found in YouTube with the link https://youtu.be/uW8TXOZnRs8) as described in the Fig. 1e (colliding particles are marked by black box). This collision phase for this particular ensemble was identified by the normalised probability distribution of the angles of the directions of the particles which shows two distinct peaks corresponding to the two oppositely moving particle groups (Fig. 3a). And these two peaks are approximately separated by π. The difference in the heights of these two peaks is justified by the strengths of the two groups. We have measured the lifetime of this collision phase by calculating n(t) denotes the number of colliding particles at the time instant t. Fig. 3b exhibits the duration of the collision time (τ c ≈ 75) by the deep well in the curve (zoomed-in plot in inset). During the collision, the value of R(t) falls (roughly at t = 135) due to two opposite contributions and then increases again (roughly at t = 210, i.e., till the collision persists) and reaches steady state value equal to unity as the colliding particles adopt the mutual uni-directional flow.

Multi-domain phase:
On decreasing the density further to ρ = 1.78, the third type of transient phase is obtained which is the existence of several domains (Fig. 1f)

Dependence of integrated lifetime on density:
In order to find the dependence of the lifetime of these transient phases on controlling parameters of the system, we steer our attention to the integrated lifetime of all these transient phases instead of treating them individually. It is quite difficult to estimate the differentiated lifetimes of all these transient phases. In this and the consecutive section, all the obtained results and prescribed relationships are tested on configurations possessing 100, 200, 400, 800 particles (N ) and 1000 ensembles of each N .
The integrated lifetime of the transient phases in the low noise regime (η = 0.1) is calculated from the time evolution of the standard deviation of the distributed angles of directions of the particles σ 2 θ = θ 2 − θ 2 . This σ 2 θ of all the particles goes to zero (values 0.004 in our simulations) when the system achieves steady ordered state; corresponding time is the integrated lifetime (τ t ) of the total transient phase (Fig. 4a, the plot represents a single realisation of a system having N = 400, L = 10 (ρ = 4.0)). In the figure, there is a knee-like bending in the curve which is in close correspondence with the vortices, but prescribing any definite relationship between the two is beyond the scope of this current study.
However, this knee-like bending is found to appear stochastically when single realisations are considered. Upon averaging over many ensembles, it washes out this knee-like bending, due to continuous symmetry and stochastic nature. Hence, in order to incorporate ensembleaverage picture, the distribution of τ t (for any fixed set of values of density and noise) is found over 1000 different initial configurations. Similarly, different distributions are found for different values of densities (Fig. 4b). Hence, the average values of τ t are calculated from these normalised distributions, for different densities, in the low noise regime. Fig. 4c shows the dependence of τ t on density ρ for N = 100, 200, 400 and 800 which reveals a power law scaling. In the log scale, τ t exhibits fair straight line, i.e., linear dependence on ρ reflecting upon the scale invariance nature of the integrated lifetime (Fig.  4d). By the method of data collapse with assumed scaling relation, τ t N −a = F (ρN −b ) and best-fitting the curves with (τ t N −a ) = q(ρN −b ) −s (Fig. 5a), the scale invariant power law exponent is obtained to be s = 1.027 ± 0.008 and the most precise values of the data collapse exponents are obtained from the error minimisation calculations of the best linear fit of τ t = m(ρ −s N (a+bs) )+c (Fig. 5b). The best fitting error is defined as Q = i (y i −mx i −c) 2 , where y i and x i represent τ t and ρ −s N (a+bs) respectively (Fig. 5b). The surface plot of the error Q(a, b) (taking a and b as axes) is shown in Fig. 5c. The global minimum value of Q(= 11581.219) is obtained for the values of the collapse exponents a = −0.110(±0.010) and b = 0.950(±0.010) (marked by a cross in Fig. 5d).
The same methodology of finding the dependence of transient lifetime on density does not work for high noise regime since the second moment (σ 2 θ ) achieves saturation at values much larger than 0.004 (which is set as the cut-off value to represent the ordered phase). Hence, in order to scan the entire range of η [0,2π], another methodology is introduced which we have further used rigorously for studying the dependence of the lifetime of transient phases on noise when the density is kept fixed (in section 3.3). However, the same is also verified for the already-mentioned low-noise, density-dependent studies of the transient lifetime as well. Fig. 5e shows the time evolution of the order parameter v a for different densities (noise fixed at η = 0.1), best-fitting of which by the relation v a = c + d{1 − exp(−t/τ t )} reveals same power law dependence of lifetime on density reported in this section.

Dependence of integrated lifetime on noise:
For studying the dependence of integrated lifetime of the transient phases on noise, density is kept fixed at ρ = 4.0 and the time evolution of the order parameter v a is observed for various values of noise ( Fig. 6c shows a representative system with N = 400 and L = 10). The saturation of the order parameter denotes the non-equilibrium steady state condition. Note that the steady state values of the order parameter are in accordance with the phase transition curve shown in Fig. 6a. The critical value of noise here is obtained from Fig.  6b in which the fluctuation in the order parameter (σ va ) attains maximum value (which is believed to diverge eventually in the thermodynamic limit) at η c (L = 10) = 3.0 for this representative case. These curves in Fig. 6c are best fitted according to the relation v a = c + d{1 − exp(−t/τ )} where c and d are fitting parameters and τ gives the integrated lifetime of the transient phases, obtained from the time evolution of the order parameter. These best fitted τ 's, when plotted against noise η, present a remarkable result. It shows that the integrated lifetime of the transient phases decreases with higher values of noise as the noise approaches the critical value (Fig. 6d), which is contrary to the usual case of critical slowing down, observed in the equilibrium kind of phase transitions. In the figure, only the representative N = 400 case is shown while this reverse critical slowing down phenomenon is similarly obtained for other values of N as well mentioned in the previous section.
We have tried to provide a possible explanation of this contrast behaviour (the time scale decreases as one approaches the critical noise) of reverse critical slowing down, with the help of the time evolution of the profile of the normalised distribution of angles of directions of the particles (Fig. 7a). When the system starts from an initial random configuration (Fig.  7b), the probability distribution of angles is uniform (denoted by bold green line in Fig.  7a), irrespective of low or high noise. But in the low noise case (η = 0.1), the distribution is sharply peaked since all the particles become uni-directional (Fig. 7c) upon reaching the steady state (the blue peak in the Fig. 7a denotes low noise steady state case for a single ensemble), whereas when the system is subjected to higher noise (η = 2.0) the angular distribution exhibits smeared peak (red line in the Fig. 7a) reflecting that the morphology of velocity portrait is quasi-ordered (Fig. 7d), it is neither entirely random nor entirely ordered, yet has a finite non-zero value of the order parameter (Fig. 6a) in the steady state. So, it is quite natural to expect that the time taken to achieve a distribution having a sharp peak (with very small width) from an initial uniform distribution (with maximum width 2π) will naturally be larger than the time taken to achieve a distribution having a comparably smeared peak (with a comparatively larger width). This can be a justifiable argument for the decrease of the integrated lifetime with the increase of noise. A video of the simulation can be found in YouTube with the link https://youtu.be/uGOhvfsyc6M.

Summary:
In this article, we have studied the dynamical evolution of the transient phases of Vicsek model of flocking, by computer simulation. Although the Vicsek model is widely studied to explore the nature of non-equilibrium ordered steady state and the transitions, the behaviours of the transient phases are somehow overlooked. We have thoroughly investigated the dynamical evolutions of the morphologies of various transient phases designed by the direction of velocity vector of each agent. T. Vicsek, in his originally proposed model of flocking (in 1995), mentioned that this model may map to XY model in the v → 0 limit. This inspired us to search for any possible vortex phase (observed in two dimensional XY ferromagnet leading to Kosterlitz-Thouless transition [20]). In the low velocity limit, we have indeed observed the vortex phase. This vortex persists (in the low noise limit) with coherent rolling over a considerable scale of time. Moreover, for the different values of low density, we have also found the colliding phase, where two groups of oppositely moving agents collide 8 before getting ordered. We have also proposed a method of characterising this colliding phase and the method of estimating the lifetime of the colliding phase. A third kind of transient phase, namely, the multi-domain phase, was also observed. In this phase, the multiple domains of agents (marked by a variety of net directions of velocity vectors) were found to form in different regions of the entire cell of study. We have measured the number of domains, strength of each domain and the lifetime of the multi-domain phase. The integrated lifetime of all these transient phases is measured statistically and found to follow a scaling relation with the density, characterised by power law variation with some exponents. We have also studied the integrated lifetime as function of the noise. The integrated lifetime τ was observed to decrease as the noise approaches the critical noise from below. This is contrary to the critical slowing down, usually observed in the case of equilibrium critical phenomena. We have tried to provide a possible explanation of this unusual behaviour, from the time evolution of the distribution of the directions of velocities as follows: the temporal evolution of the Vicsek dynamics starts from an initial random (uniformly distributed) configuration of the directions of motion. This distribution is uniform and has the maximum width (2π). In the case of low noise, the final non-equilibrium steady state is an almost uniquely directed motion of all agents having a sharply peaked distribution of the directions of motion. On the other hand, for high noise (where the value of the order parameter is quite low; near the critical noise), the steady state distribution of the angles has a wide spread. So, it is quite natural to expect that, to achieve a very sharp distribution, the system will take more time. As a result, the transient phases will be longer lasting for low noise than that for high noise. A live demonstration can be found in YouTube with the link https://youtu.be/uGOhvfsyc6M, for better understanding. Extensive numerical studies reported [19] the discontinuous nature of this flocking transition which also supports this finding, since one might not expect the critical slowing down in the case of discontinuous transitions.    Figure 7: (a) The bold green line denotes the normalised uniform distribution of angles of velocity-directions of particles for N=400 at t=0. The blue sharp unimodal peak represents the ordered steady state at t=300 for a particular realisation of the sample at low noise, the corresponding angle distribution at t=300 at higher noise is represented by the red line which exhibits a smeared peak denoting quasi-randomness; (b) Velocity morphology of the particles at t = 0 (completely random); (c) Velocity morphology of the particles at t = 300 when the system is subjected to low noise (η = 0.1) (completely ordered); (d) Velocity morphology of the particles at t = 300 when the system is subjected to high noise (η = 2.0) (quasi-random).