The kinetic regime of the Vicsek model

# The kinetic regime of the Vicsek model

## Abstract

We consider the dynamics of the system of self propelling particles modeled via the Vicsek algorithm in continuum time limit. It is shown that the alignment process for the velocities can be subdivided into two regimes: “fast“ kinetic and “slow“ hydrodynamic ones. In fast kinetic regime the alignment of the particle velocity to the local neighborhood takes place with characteristic relaxation time. So that the bigger regions arise with the velocity alignment. These regions align their velocities thus giving rise to hydrodynamic regime of the dynamics. We propose the mean-field like approach in which we take into account the correlations between density and velocity. The comparison of the theoretical predictions with the numerical simulations is given. The relation between Vicsek model in the zero velocity limit and the Kuramoto model is stated. The mean-field approach accounting for the dynamic change of the neighborhood is proposed. The nature of the discontinuity of the dependence of the order parameter in case of vectorial noise revealed in Gregorie and Chaite, Phys. Rev. Lett., 92, 025702 (2004) is discussed and the explanation of it is proposed.

self-propelling particles, Vicsek model, self organization
###### :
\layoutstyle

8x11single

address=Department for Theoretical Physics, Odessa National University, Dvoryanskaya 2, 65026 Odessa, Ukraine

address=Department for Theoretical Physics, Odessa National University, Dvoryanskaya 2, 65026 Odessa, Ukraine

## Introduction

The equilibrium statistical mechanics and thermodynamics of Hamiltonian systems are well developed areas of Statistical Physics. There are also a lot of remarkable results for open systems, which considered far from equilibrium [Haken(1983), Prigogine and Nicolis(1977)]. In general one can consider the systems with some constraints which bounds the coordinates of the particles since the general formalism remains unchanged (e.g. Liouville equation etc.). It is well known that due to instability of trajectories the Hamiltonian systems reach the thermodynamic equilibrium which can be characterized several macroscopic parameters of state despite huge number of microscopic degrees of freedom. In addition total momentum and angular momentum are conserved. Though for every specific configuration the formation of stationary local vortical structures [Erdmann et al.(2005)Erdmann, Ebeling, and Mikhailov, D’Orsogna et al.(2006)D’Orsogna, Chuang, Bertozzi, and Chayes, Chen and Leung(2006)] may occur due to conservation of the angular momentum. Obviously the inclusion of potential forces has little grounds for the systems of intellectual particles (individuals in flocks, crowds etc.), which differ very much in this respect from the molecular system for which either all forces have mainly potential character or dissipative.

In [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet] the minimal model, the so called Vicsek model (VM) of such type was introduced. The dynamic rule for the alignment of the particle’ velocities constructed in such a way that at high density the kinetic energy of disordered motion is transformed into the one of ordered motion so that the total kinetic energy is conserved. Then the system reaches the final state with nonzeroth total momentum even in the low (one- and two-) dimensional cases. In such a case the appearance of the ordered state is predetermined by the dynamic rule. Note that the VM does not take into account the potential interparticle forces.

In this paper we consider the kinetic regime for the VM when the particle aligns along the velocity direction of its neighborhood and give the estimation for the critical noise amplitude of order-disorder transition. The structure of paper is as follows. In Section 1 we derive the continuum time analog of the VM and show that angular velocity of a particle consists of two terms which describe alignment. One of the terms describes the fast kinetic relaxation to the local direction the other one describes the hydrodynamic regime of alignment between the domains where local alignment is settled down. In Section 2 the process of the relaxation of one-particle velocity to the local value of the neighborhood is considered. The dependence of the rate of the relaxation on the density is obtained by the numerical experiment and the corresponding Fokker-Plank equation is derived. In Section 3 the influence of the dynamic nature of the environment is discussed. The results obtained are summarized in the conclusion.

## 1 Vicsek model in continuum time limit and the two regimes of the dynamic

The Vicsek model of the dynamic of self-propelling particles [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet] can be represented by the relation:

 vi(n+1)×ui(n)=0,∀i,n, (1)

Here

 ui=∑jH(rij)vj∣∣∣∑jH(rij)vj∣∣∣ (2)

is the unit vector corresponding to the averaged velocity of the neighborhood and is the averaging kernel with the characteristic averaging scale . The absolute value of the velocity of each particle is assumed to be constant, i.e.

 ∣vi(n+1)∣=∣vi(n)∣=vi. (3)

The noise is not included. For the VM, is proportional to a Heaviside step function. One can also consider other models for the averaging kernel. One can say that the dynamics of individual particle is subjected to reduction of the difference between the direction of its velocity and that of the average velocity of the surrounding given by Eq. (2).

Note that at every step the direction of the velocity coincides exactly with the direction of . Another words, the vector remains unchanged during the velocities updating. This is specific for discrete formulation but this is not the case in the continuous time limit since both vectors are rotating during infinitesimal time interval . In such a limit the angular velocity consists of two parts. The first one is the angular velocity of the unit vector for the average velocity of the nearest neighbors. The second one is the relative angular velocity. When the time is continuous taking into account constraint Eq. (1) the equation of motion of such a particle can be written as:

 ddtvi=ωvi×vi. (4)

Here is the “angular velocity“ of -th particle.

This angular velocity depends on the velocities of neighboring particles. The self-propelling force and the frictional force are assumed to balance each other. The hydrodynamic model which is based on the equations of motion (1) and it continual analog (4) was proposed in [?, ?].

Note that in the discrete CVA at every step the direction of the velocity coincides exactly with the direction of . Another words, the vector remains unchanged during the velocities updating. This is specific for discrete formulation but this is not the case in the continuous time limit since both vectors are rotating during infinitesimal time interval . In such a limit the angular velocity consists of two parts. The first one is the angular velocity of the unit vector for the average velocity of the nearest neighbors. The second one is the relative angular velocity.

 ωvi=ωui+ωvu, (5)

where

 ωui= ui×˙ui (6) ωvu= Avi×ui. (7)

The quantity is inverse to characteristic length scale. The latter is the radius of interaction and is the parameter of the model. Indeed, in the limit each particle has the same neighborhood, provided that , i.e. does not depend on . Therefore, in such a limit they all has the same angular velocity, which is given by the first term of Eq. (5). From Eq. (5) and Eqs. (6), (7) for 2D case we obtain:

 ˙ωvi=−A(vi⋅ui)(ωvi−ωui)−(˙ωui+ω2ui). (8)

## 2 One-particle relaxation

Here we show that the equation (5) in 2D is closely related to well known Kuramoto model (KM) for the phase synchronization. Indeed, let the angle characterize the direction of the velocity of -th particle, then and the equation (5) takes the form:

 ˙θi=˙ψi+Asin(ψi−θi) (9)

Here denotes the angle which determines the direction vector . It is one of the variant of the short-range version of the KM [?] (see also [Acebron et al.(2005)Acebron, Bonilla, Vicente, Ritort, and Spigleri] and reference therein) of the form:

 ˙θi=ωi+K∑⟨i,j⟩sin(θj−θi), (10)

where the brackets stand for nearest-neighbor oscillators. Thus we can state that in the zero velocity limit the Vicsek model with continuum time belongs to the short-range Kuramoto model class [Acebron et al.(2005)Acebron, Bonilla, Vicente, Ritort, and Spigleri]. This allow to conclude that for low velocity the ordering in the Vicsek model is governed by the same mechanism as the synchronization in the KM. Since the synchronization transition is of continuum type one can expect that the continuum character of the transition take place for low enough velocity in Vicsek model too. This conclusion is in correspondence with the results of [Nagy et al.(2007)Nagy, Daruka, and Viscek].

According to its definition vector changes slower than the velocity of a particle . From Eq. (8) it follows that the second term in Eq. (5) governs the kinetic of the alignment process while the first term is of hydrodynamic character since it determines the behavior on scales greater than . Therefore Eq. (8) shows that in continuous time limit the CVA system has the stable state where the particles align along some direction, which is characterized by the director . Equations (4), (5) can serve as the continuum time analog for the CVA. Additional confirmation of that is the behavior of these terms with respect to reverting the time . The first term changes its sign and therefore produces the reversible contribution to the equation of motion, while the second term does not change the sign thus representing irreversible part of the CVA, which governs irreversible one-particle kinetics of the alignment. To study Eq. (8) analytically we use the approximation which takes into account that the variable is the “collective“ one, thus there is the time interval which we call “kinetic“ regime where it changes much slower so that and its derivative can be omitted. In addition we assume that the value of does not depend on time which reflects that the number of “interacting“ neighbors remains constant. In such an approximation Eq. (8) reduces to more simple form:

 ˙ωvi=−A(vi⋅ui.)ωvi. (11)

In scalar form for the angle of the alignment between the vectors and taking into account that and

 vi⋅ui=cosα,

after integrating Eq. (11) we obtain:

 ˙α=−Asinα. (12)

where is some parameter which determines the alignment rate and obviously depends on the density and the average velocity of the neighbors. Here we put the following initial condition , which is in accordance with that in simulations. The solution of Eq. (12) is:

 tanα2=tanα02exp(−At), (13)

Thus the one-particle alignment process has relaxation character.

To compare this result with the simulation data we performed the simulation with the small initial disalignment of the directions of the particles in dense system . The results of simulation are represented on Fig. 1 and demonstrate the existence of such an interval, where the dependence given by Eq. (11) takes place.

The kinetic regime of the system subjected to the stochastic increment of the direction [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet, Czirök et al.(1996)Czirök, Ben-Jacob, Cohen, and Viscek] can be described by stochastic modification of Eq. (12):

 ˙α=−Asinα+L(t), (14)

where is the standard white noise term. Then Eq. (14) is equivalent to Fokker-Plank equation for the density distribution function [Van Kampen(2001)]:

 ∂fv∂t=∂∂α(Asinαfv)+D∂2fv∂α2, (15)

where is the diffusion coefficient and we use the approximation:

 A=λu (16)

for the alignment rate for the dimension reasons, where is the local density. Such density dependence of the relaxation rate is supported by the simulations (see Fig. 2).

The case corresponds to the deterministic case, with the solution:

 f(0)v(α,t)=f(eλuttanα2)sinα, (17)

which demonstrates the alignment process in accordance with Eqs. (12), (13). The stationary solution of Eq. (15) is:

 f(st)v(α)=C(D)eλuDcosα. (18)

Note that the distribution (18) was used in [Csahok and Viscek(1995)] for the lattice model as the analog of the Boltzman distribution. Thus the consideration presented above can serve as the ground for such representation.

## 3 The account of dynamical environment

In the zero velocity limit the Vicsek model can be considered in terms of the lattice model as the systems of interacting spins with the interaction favoring the alignment and the emergence of the long range order. Yet there is the major difference between the CVA and the equilibrium models. This is the coupling between the density and the velocity fields. Due to such coupling in the static () case the the equilibrium systems does not order for densities below some threshold value close to 1 (which corresponds to the percolation threshold of randomly distributed spheres) while in the SSP ordering is found for all velocities [Nagy et al.(2007)Nagy, Daruka, and Viscek]. The instant change of the environment in the neighborhood of the particle can be considered as the noise factor correlated with the local value of the order parameter - the average velocity of the neighbors:

 ui=1Ni∑⟨i,j⟩vi, (19)

where stands for the nearest neighbors of -th particle and is the number of such neighbors. Thus the vector is just the weighted sum of the random vectors. The number of summands is also random and describes the dynamic coupling between the density and the velocity of the neighbors. Let one-particle distribution function is . Since the order parameter is the collective variable as has been said above it changes more slowly than the one-particle function. Thus one can consider as the parameter for distribution function . Using the standard procedure analogous to the mean-field approach it is possible to get the self-consistent equation for the order parameter. Indeed, assuming that all the correlation between the particles are incorporated into we can find the relation between the characteristic functions for and :

 Missing or unrecognized delimiter for \right (20)

Here is the probability of neighbors to appear within the interaction range and

 Gv=⟨eikv⟩ (21)

is the characteristic function of the velocity distribution of a particle . It depends on the average value of the local order parameter (19). In order to get the specific results for the order parameter one need to specify the form of the density distribution for the number of neighbors and the distribution function for the velocity. The latter depends on the type of noise introduced into the motion. The original variant of the CVA [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet, Czirök et al.(1999)Czirök, Viscek, and T.Viscek] used so called scalar noise when the direction of a particle is updated with the random increment of the angle. In recent work [Grégoire and Chaté(2004)] another type of noise, so called vectorial noise, was proposed. We will consider these type of noise below.

To simplify the consideration we assume that the probability distribution does not depend on the distribution for the velocity. The simplest choice for the distribution of the number of particles in nearest surrounding is the Poisson distribution:

 pn=λnn!e−λ, (22)

where is mean number (density) of particles in the nearest surrounding. In order to get the analytic result in such a case it is expedient to use the non normalized order parameter instead of (19):

 ~ui=∑⟨i,j⟩vj. (23)

Then we get simple expression for the characteristic function :

 Extra open brace or missing close brace (24)

The standard relation between the moments of the distribution and the derivatives of the characteristic function at leads to:

 ⟨~u⟩=λ⟨v⟩⇔⟨u⟩=⟨v⟩. (25)

Note that such a simple relation between the order parameter and the particle velocity is due to specific form of the Poisson distribution. It allows to justify the expression for the relaxation rate Eq. (16) in low density limit. The velocity distribution function depends on the specific way of introducing the source of noise. Below we consider two types of noise which are widely used in simulations.

### 3.1 Scalar noise

Here we consider the case of so called scalar noise. This is the most obvious way to introduce the stochastic source into the dynamic of a particle. At every step the random increment of the angle is added. Using the results obtained above for the velocity distribution function, we choose it in accordance with Eq. (18)

 fv=CeλuDcosα,C=12πI0(λuD). (26)
 u=I1(λuD)I0(λuD) (27)

It’s obvious, that it has trivial solution , which losses its stability depending on the average density and the diffusion coefficient (i.e. the intensity of scalar noise). Expanding Eq. (27) near the trivial solution we obtain :

 0=(λ2D−1)u+λ16D3u3+o(u3). (28)

From here the critical density value is as following:

 λc=2D (29)

The comparison of the solution of Eq. (27) and the results of numerical simulation obtained in [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet, Czirök and T.Viscek(2000)] is on Fig. 3.

From Eq. (20) it clear that behavior of the order parameter near the critical value is determined by the analytical properties of the characteristic function . Near the critical threshold the dependence of the order parameter has typical for the mean-field approximation square root dependence:

 u∝(λc−λ)1/2. (30)

The differentiability properties of the characteristic function Eq. (20) and therefore the applicability of the expansion Eq. (28) depends on the character of the distribution of the neighbors. If the fluctuations of the density near the threshold value are big this can lead to slow convergence of . In such a case one can expect non smooth dependence of on the parameters of the distribution function , in particular on the average value of the local order parameter .

In general one need to construct the kinetic equation for the distribution function. Some attempts to derive such equation have been made in a way similar to classic Boltzmann approach [Bertin et al.(2006)Bertin, Droz, and Grégoire] though only binary collisions were taken into account. This approximation is valid only for the system of low density. The applicability of these result to the systems of Vcsek’s type is problematic because of the multiparticle character of the “collision“ process.

### 3.2 Vector noise

The vectorial noise was introduced in [Grégoire and Chaté(2004)] as another realistic model of the noisy environment. In such a case the random vector is added, either to the local order parameter [Grégoire and Chaté(2004)] or [Aldana et al.(2007)Aldana, Dossetti, Huepe, Kenkre, and Larralde, Pimentel et al.(2008)Pimentel, Aldana, Huepe, and Larralde]. Then the corresponding direction for the velocity of the particle is determined:

 θi(n+1)=Arg(ui(n)+ξi) (31)

In addition, the amplitude of the noise can be chosen so that [Grégoire and Chaté(2004)]. The results obtained in [Grégoire and Chaté(2004)] revealed the difference between the VM with the scalar noise and raised the intensive discussion (see [Nagy et al.(2007)Nagy, Daruka, and Viscek, Aldana et al.(2007)Aldana, Dossetti, Huepe, Kenkre, and Larralde, Pimentel et al.(2008)Pimentel, Aldana, Huepe, and Larralde, Chaté et al.(2007)Chaté, Ginelli, and Grégoire]).

Here we derive the one-particle velocity distribution function for the case of vectorial noise and show that it has essentially nonlinear character which leads to the apparent discontinuity in the dependence of the order parameter on the noise intensity.

We assume that the distribution of the direction of the vector is uniform and independent on the the distribution of the number of neighbors, which is characterized by the corresponding probabilities .

From simple geometrical consideration of vector noise algorithm it is easy to get the relation:

 cosα=1 ⎷1+sin2φ(ξu+cosφ)2,ξ

If then the distribution function for the direction is as following:

 f(α)=12π⎛⎜ ⎜ ⎜ ⎜⎝1+cosα√(ξu)2−sin2α⎞⎟ ⎟ ⎟ ⎟⎠ (33)

Thus the self-consistent equation for the order parameter is as following:

 u=⟨v⟩=F(ξu) (34)

where

 Unknown environment '%'

The solution of Eq. (34) is shown on Fig. 4. There is the interval of the noise intensity , where two nontrivial solutions for the order parameter exist with the hysteretic (or subcritical) jump. It is clear that the branch where represents the unstable state in analogy with the situation typical for the first order phase transitions. For the model considered and . The situation here is analogous with that for the Kuramoto model [Acebron et al.(2005)Acebron, Bonilla, Vicente, Ritort, and Spigleri]. In the latter case the type of bifurcation of the partially synchronized phase depends on the properties of the frequency distribution function of the oscillators, namely the sign of (see [Acebron et al.(2005)Acebron, Bonilla, Vicente, Ritort, and Spigleri]). Thus we state that the difference of the one-particle distribution function in case of scalar and vector noise is the source of the change the type of the bifurcation from supercritical for scalar noise to the subcritical for vector noise in the Vicsek model. This can explain the difference in the results of works [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet] and [Grégoire and Chaté(2004)].

## 4 Conclusion

In this paper we consider the continuum time limit for the Vicsek’s algorithm [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet] of the dynamics of the self propelling particles. It is shown that there is the time interval where the kinetic regime of the relaxation of the particle velocity to the local value of the average velocity of the neighbors takes place. The relaxation rate depends on the density linearly at least for not too big number of neighbors. The cases of vectorial and the scalar noises are considered. Within the proposed mean field model it is shown that the for the case of vectorial noise [Grégoire and Chaté(2004)] the subcritical bifurcation of the stationary solution takes place. This is in contrast with the case of scalar noise originally considered in [Viscek et al.(1995)Viscek, Czirök, Ben-Jacob, Cohen, and Shochet], where the supercritical transition occurs.

### References

1. H. Haken, Synergetics: An Introduction (Springer Verlag, Berlin, 1983).
2. I. Prigogine and G. Nicolis, Self-Organization in Non-Equilibrium Systems: From Dissipative Structures to Order Through Fluctuations (J. Wiley & Sons, New York, 1977).
3. U. Erdmann, W. Ebeling, and A. Mikhailov, Phys. Rev. E 71, 051904 (2005).
4. M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Phys. Rev. Lett. 96, 104302 (2006).
5. H.-Y. Chen and K.-t. Leung, Phys. Rev. E 73, 056107 (2006).
6. T. Viscek, A. Czirök, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
7. V. Kulinskii, V. Ratushnaya, A. Zvelindovsky, and D. Bedeaux, Europhys. Lett. 71, 207 (2005).
8. V. Ratushnaya, V. Kulinskii, A. Zvelindovsky, and D. Bedeaux, Physica A: Statistical Mechanics and its Applications 366, 107 (2006), ISSN 03784371, URL http://dx.doi.org/10.1016/j.physa.2005.11.002.
9. H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, Prog. Theor. Phys. 77, 1005 (1987).
10. J. A. Acebron, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigleri, Rev. Mod. Phys. 77, 137 (2005).
11. M. Nagy, I. Daruka, and T. Viscek, Physica A 377, 445 (2007).
12. A. Czirök, E. Ben-Jacob, I. Cohen, and T. Viscek, Phys. Rev. E 54, 1791 (1996).
13. N. G. Van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland Personal Library) (North Holland, 2001), ISBN 0444893490.
14. Z. Csahok and T. Viscek, Phys. Rev. E 52, 5297 (1995).
15. A. Czirök, M. Viscek, and T.Viscek, Physica A 264, 299 (1999).
16. G. Grégoire and H. Chaté, Phys. Rev. Lett. 92, 025702 (2004).
17. A. Czirök and T.Viscek, Physica A 281, 17 (2000).
18. E. Bertin, M. Droz, and G. Grégoire, Physical review. E, Statistical, nonlinear, and soft matter physics 74 (2006), ISSN 1539-3755, URL http://view.ncbi.nlm.nih.gov/pubmed/17025488.
19. M. Aldana, V. Dossetti, C. Huepe, V. M. Kenkre, and H. Larralde, Physical Review Letters 98, 095702 (pages 4) (2007), URL http://link.aps.org/abstract_prlv98/e095702.
20. J. A. Pimentel, M. Aldana, C. Huepe, and H. Larralde, Physical Review E (Statistical, Nonlinear, and Soft Matter Physics) 77 (2008), URL http://dx.doi.org/10.1103/PhysRevE.77.061138.
21. H. Chaté, F. Ginelli, and G. Grégoire, Physical Review Letters 99, 229601 (pages 1) (2007), URL http://link.aps.org/abstract_prlv99/e229601.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters