Collection of polar self-propelled particles with a modified alignment interaction

We study the disorder-to-order transition in a collection of polar self-propelled particles interacting through a distance dependent alignment interaction. Strength of the interaction, $a^{d}$ ($0<a<1$) decays with metric distance $d$ between particle pair, and the interaction is short range. At $a = 1.0$, our model reduces to the famous Vicsek model. For all ${\it a}>0$, the system shows a transition from a disordered to an ordered state as a function of noise strength. We calculate the critical noise strength, $\eta_c(a)$ for different $a$ and compare it with the mean-field result. Nature of the disorder-to-order transition continuously changes from discontinuous to continuous with decreasing $a$. We numerically estimate tri-critical point $a_{TCP}$ at which the nature of transition changes from discontinuous to continuous. The density phase separation is large for ${\it a}$ close to unity, and it decays with decreasing $a$. We also write the coarse-grained hydrodynamic equations of motion for general ${\it a}$, and find that the homogeneous ordered state is unstable to small perturbation as ${\it a}$ approaches to $1$. The instability in the homogeneous ordered state is consistent with the large density phase separation for ${\it a}$ close to unity.


I. INTRODUCTION
Flocking [1][2][3][4], the collective and coherent motion of large number of organisms, is one of the most familiar and ubiquitous biological phenomena. In the last one decade, there have been an increasing interest in the rich behaviors of these systems that are different from their equilibrium counterparts [5][6][7]. One of the key features of these flocks is that there is a transition from a disordered state to a long ranged ordered state in two-dimensions with the variation of system parameters (e.g., density, noise strength) [8][9][10]. The study of the phase transition in these systems is an active area of research, even after many years since the introduction of the celebrated model by Vicsek et. al. [8]. Many studies have been performed with different variants of metric distance model [11] and topological distance model [12][13][14]. In the Novel work of Vicsek, it is observed that the disordered to ordered state transition is continuous [8], but later other studies [9,10] confirmed that the transition is discontinuous. Some studies on the topological distance model claim the transition to be discontinuous [14], whereas other studies [12,13] find it continuous. Therefore, the nature of the transition of polar flock is still a matter of debate.
In our present work we ask the question, whether the nature of transition in polar flock can be tuned by tuning certain system parameters. And how do the characteristics of system change for the two types of transitions (discontinuous / continuous) ? To answer this, we introduce a distance dependent parameter a such that the strength of interaction decays with distance. For a = 1, the interaction is same as that in the Vicsek model. For all non-zero distance dependent parameter (a > 0), the system is in a disordered state at small density and high noise strength, and in an ordered state at high density and low noise strength. We calculate the critical noise strength η c (a) for different a and compare it with the mean-field result. The nature of the disorder to order transition continuously changes from discontinuous to continuous with decreasing a. We estimate the tri-critical point in the noise strength η and a plane, where the nature of the transition changes from discontinuous to continuous. We also calculate the density phase separation in the system. The density phase separation order parameter is large for a close to unity, and it monotonically decays with decreasing a . Linear stability analysis of the homogeneous ordered state shows an instability as a approaches to 1, which is consistent with large density phase separation for a 1.
This article is organised as follows. In section II, we introduce the microscopic rule based model for distance dependent interaction. The results of numerical simulation are given in section III. In section IV, we write the coarse-grained hydrodynamic equation of motion, calculate the mean field estimate of critical η c (a), and discuss the results of linear stability analysis. Finally in section V, we discuss our results and future prospect of our study. Appendix A is at the end, that contains the detailed calculation of the linear stability analysis.

II. MODEL
We study a collection of polar self-propelled particles on a two-dimensional substrate. The particles interact through a short range alignment interaction, which decays with the metric distance.
Each particle is defined by its position r i (t) and orientation θ i (t) or unit direction vector n i (t) = [cos θ i (t), sin θ i (t)]. Dynamics of the particles are given by two update equations. One for the position and other for the orientation. Self-propulsion is introduced as a motion towards its orientation with a fixed step size(v 0 in unit time). Hence, the position update equation of the particles and the orientation update equation with a distance dependent short range alignment interaction where the sum is over all the particles within the interaction radius (R 0 ) of the i th particle, i.e., is the number of particles within the interaction radius of the i th particle at time t, and d is the metric distance between a pair of particles (i, j). W i (t) is the normalisation factor. The strength of the noise η is varied between zero to 1, and ζ i (t) is a random unit vector. Note that this model reduces to the celebrated Vicsek model for a = 1.0.

III. NUMERICAL STUDY
We numerically simulate the microscopic model introduced by Eqs.1 and 2 for different distance dependent parameter a. For a = 1, the particle interacts with the same strength with all the particles inside its interaction radius (Vicsek's model [8]). As we decrease a, interaction strength decays with distance. a is varied from 1.0 to small value 0.001. For a = 0.0 the particles are non-interacting. Speed of the particles is fixed to v 0 = 0.5. We start with random orientation and homogeneously distributed particles on a 2−dimensional substrate of size L × L with periodic boundary conditions. For all the simulations, we keep mean density ρ 0 = N L 2 = 2.0. Number of particles were varied from N = 1000 to 10000. We start from a random state and each particle is updated using Eqs. 1 and 2. One simulation step is counted after sequential update of all the particles. All the measurements are performed after 10 5 simulation steps, and a total of 10 6 steps are used in simulations.

A. Disorder-to-order transition
First we study the disorder-to-order transition in the system for different a. Ordering in the system is characterised by the global velocity, In the ordered state, i.e., when large number of particles are oriented in the same direction, then V is close to 1, and it is close to zero for a random disordered state. In Fig. 1 (a-d)  V (t) shows fluctuations, but there is no switching behaviour. We further calculate probability distribution P (V ) of the global velocity for the same set of a and η values as used for the time series plots. As shown in Fig. 2(a), P (V ) is bimodal for a = 1.0, i.e., there are two distinct peaks for P (V ). Two finite values of V corresponds to two states of the system. Two peaks come closer with decreasing a, and for small a(= 0.01), P (V ) shows only one broad peak in Fig. 2(d). The bimodal distribution of the V confirms that the transition is discontinuous for a 1.
To further characterise the nature of the transition, we calculate the fourth order cumulant or the Binder cumulant, i.e., U (η) vs. η plot is shown in Fig. 3. It shows strong discontinuity from U = 1/3 (for disordered state) to U = 2/3 (for ordered state) as we approach critical η c (a) for a = 1 in Fig. 3 (a), and discontinuity decreases with a. It smoothly goes from a disordered state (U = 1/3) to an ordered state (U = 2/3) for a = 0.01 in Fig. 3 (d). For a 0.4, U vs. η plot shows strong discontinuity at large N , but for a 0.4 it becomes continuous.
Therefore, The nature of the transition continuously changes from discontinuous to continuous on decreasing a. The critical noise strength η c (a) also decreases with decreasing a. We plot η c (a) vs. a in the of Fig. 4. The solid line indicates the nature of the disorder-to-order transition is discontinuous, and the dashed line indicates the continuous transition. The value of a at which the above transition changes from discontinuous to continuous one, we call it as tri-critical-point (TCP) a T CP . For a > a T CP the transition is discontinuous, and for a < a T CP it is continuous. TCP shows a small dependence on N for any fixed v 0 and ρ 0 . We define the TCP for any system size as the point where the Binder cumulant U starts to show discontinuous variation. In the upper inset of Fig. 4, we plot 1 − a T CP vs. 1/N , and extrapolate the TCP for N → ∞ or 1/N → zero. As 1/N approaches to zero, 1 − a T CP ≈ 0.61. Hence the a T CP is ≈ 0.39. Hence, the extrapolated value of a T CP matches well with the a T CP in phase diagram, which is marked as blue square in Fig. 4. In the lower inset of Fig. 4, we plot the critical η c (a) vs. a on semi-log scale and compare the results with the mean field result in Eq. 13. Mean field approximation is good when density distribution is homogeneous. In such limit, density at each point is close to the mean density of the system.
As shown in Fig. 5 density distribution becomes more and more inhomogeneous as we increase a.
Hence, for the small a values numerical estimate of η c (a) should be more close to MF. We show in lower inset of Fig. 4 the numerical η c (a) matches very well with MF for small a < 0.1.

B. Density phase separation
The density distribution of particles also changes as we vary a. Density fluctuation plays an important role in determining the nature of the transition in polar flock [10,15,16,[20][21][22][23]. In Fig.   5 we show the real space snapshot of particle density for different a(= 1, 0.5, 0.4 and 0.01) close to critical noise strength η c (a). Clusters are small and homogeneously distributed for small a, but as a approaches to 1 we find large, dense and anisotropic clusters. We quantify the density distribution by calculating the density phase separation order parameter in Fourier space defined as, where k = 2π(m,n) L is a two dimensional wave vector and m, n = 0, 1, 2 ...., L − 1 . The reference frame is chosen so that the orthogonal axes (1, 0) and (0, 1) are along the boundary of the substrate, and (1, 1) represents diagonal direction. We calculate the first non-zero value of Q(k) in all three directions Q(1, 0), Q(0, 1) and Q(1, 1). The average density phase separation order parameter < Q > is (Q(1, 0) + Q(0, 1) + Q(1, 1))/3.
We also characterize the density phase separation using the standard deviation in particle number ∆φ in a unit size sub-cell. It is defined as where φ j is the number of particles in the j th sub-cell. To calculate ∆φ we first divide the whole system into N c (= L 2 ) unit sized sub-cells, then calculate the number of particles in each sub-cell, and from there we calculate the standard deviation in particle distribution. Q(t) and ∆φ(t) are calculated at different times in the steady state, and then average over a large time to obtain < Q > and < ∆φ > respectively. Plots of < Q > and < ∆φ > vs. a on log-log scale are shown on Fig. 6 (a) and (b) respectively. For a 1 both < Q > and ∆φ are large; however , as we decrease a, they decay monotonically. For a close to unity both < Q > and < ∆φ > show fast decay (exponential), and for smaller a they decay algebraically with a. In the insets of Fig. (a) and (b), we show the exponential decay of the density phase separation order parameter < Q > (∼ e 0.46a ), and the standard deviation in particle distribution < ∆φ > (∼ e 0.33a ) for a ≈ 1. We find that for a ≈ 1, the density phase separation is high, and the nature of the disorder-to-order transition is also first order. Hence, the change in the nature of both the disorder-to-order transition and the density phase separation shows variation on decreasing a.

IV. HYDRODYNAMIC EQUATIONS OF MOTION
We estimate the η c (a) and also study the linear stability of homogeneous ordered state with varying a. The coarse-grained hydrodynamic variables are coarse-grained density ρ(r, t) and velocity V (r, t) and they are defined as, We can write the coupled hydrodynamic equations of motion for density and velocity as obtained in Toner and Tu [18] and for velocity For our distance dependent model we have introduced an additional general a dependence to alignment parameter α(ρ, η, a) in the velocity equation 10. In [18] α is treated as a constant. But in general α is a function of microscopic parameters (e.g. density, noise strength etc.) when derived from microscopic model. For a = 1, our model reduces to the Vicsek's model, and α = α 0 (ρ − ρ c ).
ρ c in general depends on system parameters (viz: noise strength, speed etc.) On increasing density large noise is required to break the order or ρ c increases with η. Using mean-field-like argument it can be shown that ρ c η 2 v 2 0 [10] or α = α 0 (ρ − 4η 2 ). α shows linear dependence on ρ for a = 1, when all the particles within the coarse-grained radius interact with same strength. In general for a < 1, strength of interaction decays with distance. Again using the mean-field limit when density inside the coarse-grained radius is homogeneous, following form of α is obtained Hence α changes sign at critical η c .
The homogeneous solution for the disordered state is V 0 = 0 (for η > η c ), and for the ordered state In Fig. 4 (lower inset) we plot the function η c (a) vs. a as given in Eq.13 on semi-log scale and its comparison to numerically estimated η c (a). We find that the data matches very well with numerical result for small a limit. Deviation from the MF expression increases with increasing a when the density distribution becomes more inhomogeneous Fig. 5.
Now we study the linear stability analysis of Eqs. 9 and 10 about the homogeneous ordered state for general a. Detail steps of linear stability analysis are given in the appendix A. We find that for large a homogeneous ordered state is unstable with respect to small perturbation. The condition for the instability is obtained in Eq. A13. where ] . Hence, using the expression for α from Eq. 11 we get condition for instability of the hydrodynamic mode, We plot F (q, a) = α 0 (a ln(a)+1−a) (ln(a)) 2 Fig. 7, and find that the instability of the hydrodynamic mode increases with a. Unstable homogeneous state for a ≈ 1 is consistent with the large density phase separation obtained in numerical simulation. System shows first order disorder-to-order transition for large a. As we decrease a the nature of the transition changes continuously, and also the density phase separation decays. We construct the phase diagram in the noise strength and the distance dependent parameter (η, a) plane. The nature of the disorder-to-order transition is first order for a 1, and it changes to continuous type with decreasing a, and at a tri-critical point the nature of the transition changes from discontinuous to continuous. Earlier studies of [15,16] find that the disorder-to-order transition in polar flock can be mapped to the liquid-gas transition. In our study, we find that the density plays an important role and the large density inhomogeneity leads to the discontinuous transition in these systems. The effect of density is characterized by the phase separation order parameter < Q > and the standard deviation in number of particles in unit sized sub-cells < ∆φ > for different a. We find that the density phase separation is large for a 1, and it decays with decreasing a. Hence, the discontinuous disorder-to-order transition and the large density phase separation are common for a approaching to unity.
Our study concludes that the nature of the disorder-to-order transition in collection of polar flock is not always necessarily first order, and it strongly depends on the interaction amongst the parti-cles. The study of [11] shows that the transition from random to collective motion changes from continuous to discontinuous with decreasing restriction angle. The critical noise amplitude also decreases monotonically on decreasing the restriction angle. In our model we propose a parameter a, which can also tune the nature of such transition. Our model would be useful to study the disorder-to-order transition in biological and granular systems, where interaction between close-by neighbours is stronger than the interaction of particles with other neighbours.
Now we introduce fluctuations in hydrodynamic equation for density and if we consider only linear terms then Eq.9 will reduce to, We consider the velocity fluctuation only in the direction of orientational ordering. So δV y and q y is zero in our analysis. Now density Eq. A2 we can write as, Similarly we introduce fluctuations in velocity Eq. 10 and we are writing velocity fluctuation equation for ordering direction. We also introduce functional density dependency in α(ρ). We have done Taylor series expansion of α(ρ) in Eq.10 at ρ = ρ 0 , and consider upto first order derivative term of α(ρ). Now velocity equation will reduces to, where α 1 = ∂α ∂ρ | ρ 0 also λ is combination of three λ s(λ = λ 1 + λ 2 + 2λ 3 ) terms. Now considering no fluctuation along perpendicular direction of velocity field, equation along ordering direction(x-direction) reduces to, Now we are introducing Fourier component, ∆Y (q, S) = dr exp(iq.r) exp(St)dt in above two fluctuation equations A3, A5 . Then we are writing the coefficient matrix for the coupled equations.
Here we are writing q x = q.
Earlier study [17,19] finds horizontal fluctuation or fluctuation in the direction of ordering is important when system is close to transition. Here important thing is that unlike isotropic problem d > 2 there is no transverse mode, we always have just two longitudinal Gold-stone modes associated with δρ and V x . We get solution for hydrodynamic modes in symmetry broken state, where the sound speeds, and the damping ε ± in the Eq. A7 are O(q 2 ) and given by, So real part of the modes are − ± . Now we know the instability conditions are 1) If Re[S ± ] > 0 we will get homogeneous polarized state, which is unstable. 2) If Re[S ± ] < 0 we will get homogeneous polarized state, which is stable to small perturbation. We know the expression for ± , Close to transition point α 0. So we can write, We have checked Re[S − ] = − − < 0 always holds, so this mode is always stable. Re[S + ] = − + > 0 for and then this mode becomes unstable.