Consequence of anisotropy on flocking: the discretized Vicsek model

We numerically study a discretized Vicsek model (DVM) with particles orienting in q possible orientations in two dimensions. The study investigates the significance of anisotropic orientation and microscopic interaction on macroscopic behavior. The DVM is an off-lattice flocking model like the active clock model (ACM; Chatterjee et al 2022 Europhys. Lett. 138 41001) but the dynamical rules of particle alignment and movement are inspired by the prototypical Vicsek model (VM). The DVM shows qualitatively similar properties as the ACM for intermediate noise strength where a transition from macrophase to microphase separation of the coexistence region is observed as q is increased. But for small q and noise strength, the liquid phase appearing in the ACM at low temperatures is replaced in the DVM by a configuration of multiple clusters with different polarizations, which does not exhibit any long-range order. We find that the dynamical rules have a profound influence on the overarching features of the flocking phase. We further identify the metastability of the ordered liquid phase subjected to a perturbation.


I. INTRODUCTION
Active matter systems represent a fascinating class of materials composed of self-propelled entities that convert energy into mechanical motion, leading to complex and often out-of-equilibrium behaviors [1][2][3][4].The emerging phenomena in active matter systems have gained significant attention in recent years due to their potential applications in various physical, biological, and engineering systems [2,5,6].Active matter exhibits dynamic behaviors such as collective motion [2], pattern formation [7,8], and even the ability to exhibit controlled transport [9].These systems encompass a wide range of physical, chemical, and biological entities, from swimming bacteria [10], mammalian herds [11], fish schools [12,13], and sterling flocks to amoeba and bacteria colonies [14], to the cooperative behavior of cytoskeletal filaments and molecular motors in living cells [9,15,16] or in vitro environments to synthetic colloidal particles equipped with motors [17,18].To understand and unravel the fundamental principles governing active matter systems, new models [8] have been developed in the last two decades.
The Vicsek model (VM), introduced by Vicsek and collaborators in 1995 [19], provides a fundamental framework for studying the collective behavior of particles under aligning interactions.In this model, particles adjust their velocities to align with the average velocities of neighboring particles, leading to the emergence of co-herent motion and ordered patterns.VM has played a crucial role in advancing our understanding of flock dynamics [20][21][22][23].At low particle density or high noise, the particles move in random directions, and no long-range order is observed.The transition from the gas phase at high noise and low density to the polar ordered Toner-Tu phase at low noise and high density, displaying longrange order (LRO) by a coherent motion of all particles, is first order [24].But, in contrast to conventional firstorder phase transition scenarios, the coexistence phase of the VM can manifest as either multiple bands of particles moving collectively, a phenomenon known as microphase separation [23,25], or a polar-ordered cross sea phase [26], primarily driven by giant number fluctuations (GNF) [23].
Nevertheless, it is important to note that the VM posits a continuous range of possible directions for the motion of particles.However, when considering a scenario in which particles are constrained to discrete, equidistant angular orientations within a twodimensional plane, such as in the active clock model (ACM) [27,28], the VM-inspired dynamical principles governing particle alignment and movement remain uncharted territory.In a recent study on the ACM [27], it was revealed that in large systems, any values of discrete orientations result in significant changes in phenomenology when compared to the VM.These changes include the loss of long-range correlations, the pinning of global order, and the transformation of coexistence bands into a single moving domain.Additionally, another study on the ACM [28] with different dynamical rules shows that for a small number of directions, the ACM mirrors the active Potts model (APM) [29,30], exhibiting macrophase separation of the coexistence region and reorientation transition of the ordered band from transverse to longi-tudinal motion as bias velocity is increased.Conversely, with more directions, the ACM transitions towards the VM, displaying microphase separation and transversely moving bands without the reorientation transition.Remarkably, the transition in the q → ∞ limit of ACM [28], known as the active XY model, shares the same universality class as the VM.Motivated by these findings, in this paper, we undertake an extensive computational investigation that examines in detail a q-state discretized version of the Vicsek model (DVM) where the rules of particle alignment and movement follow the Vicsek protocol.
We ask several intriguing questions that persist within the context of the DVM, e.g., (a) How does discretizing the directions of particles in the VM, affect the overall diverse collective dynamics and steady-state phases?(b) What is the impact of q and system size on the coexistence region (micro-or macro-phase separation)?(c) How do the behavior of the density fluctuations, the direction of the system's global order, and the behavior of correlation functions in the liquid phase correspond with the self-organized patterns in the phase-coexistence region?(d) What is the nature of the DVM liquid phase as a function of q?To answer these questions, we study the DVM in an off-lattice domain, focusing on the three key factors: the anisotropy and degeneracy parameter q, noise level η, and system size.
The paper is organized as follows: after introducing the model in Sec.II, we present our numerical results in Sec.III.Finally, we summarize and discuss the implications of our findings in Sec.IV.

II. MODEL
We consider N self-propelled particles within a twodimensional off-lattice domain of size L x × L y (L x > L y for rectangular domain and L x = L y = L for square domain) with periodic boundary conditions.Akin to the two-dimensional VM, each point particle i is endowed with an off-lattice position vector r i = (x i , y i ) and moves with a constant speed v 0 in individual directions given by a unit orientation vector σ i = (cos θ i , sin θ i ) with an orientation angle θ i ∈ (0, 2π) where and } denote discrete orientations of the particles.q denotes the ground state degeneracy where each particle can only have discrete orientations allowed by the q value and therefore, the continuous U (1) symmetry of the VM is replaced by the discrete Z q symmetry.
At each discrete time step ∆t = 1, a particle i with velocity v 0 moves a fixed distance v 0 ∆t and interacts with N i neighboring particles within a circle of unit radius around it.The position evolves in the following way: while the new orientation is determined by a projection of the updated orientation proposed by the Vicsek rule onto one of the q allowed directions: where P is the projection and θt i is the orientation angle of a spin-weighted sum of orientation vectors of neighboring particles.ξ t i ∈ [−π, π] is a scalar noise uniformly distributed and uncorrelated for all sites and times.Such noise is often called white since it has a flat Fourier spectrum.η is the noise amplitude.
We define the projection onto the allowed directions probabilistically by where θ 1 and θ 2 are the two allowed directions which are closest to θ, such that θ 1 = 2πn/q < θ and θ 2 = 2π(n + 1)/q > θ for some n.The probability p ∈ [0, 1] is given by p = q 2π (θ − θ 1 ), i.e., minimal (0) for θ close to θ 1 and maximal (1) for θ close to θ 2 (cf.Fig. 1).Note that then for all particles going into the direction, say, θ 1 , the probability to turn stochastically into another direction, say θ 2 , is of the order ∼ q 2 η, i.e. small for small q and noise η.
Notice that the dynamical rules governing DVM are different from the q-state ACM [28].The q-state active clock model (ACM) [28] is a natural discretization of the VM in 2d where particles move in q equidistant angular directions with an alignment interaction inspired by the ferromagnetic equilibrium clock model.In the ACM, the hopping rate of a particle in state θ along any direction ϕ is D(1 − ε) for ϕ ∈ [0, 2π] and Dε for ϕ = θ.Here ε is the self-propulsion "velocity" which indicates asymmetric diffusion and D is the diffusion constant.On the contrary, in the DVM, the hopping probability of a particle with orientation angle θ along another discrete direction ϕ is zero as the particle always follows its orientation [see Eq. ( 2)].Hence, while the likelihood of hopping in non-preferred directions remains nonzero within the ACM (and also depends on the self-propulsion velocity of the particle), such movement is impossible within the DVM.It should also be noted that for ε = 0, the ACM [28] transforms to a (diffusive) Brownian clock model whereas the v 0 = 0 limit makes the DVM purely passive and the q-state DVM reduces to the equilibrium q-state clock model.In light of the above discussion, it is then evident that the transverse fluctuations say in q = 4 ACM, are stronger than the q = 4 DVM.Transverse fluctuations in the ACM mainly originate from the nonzero hopping probability of a particle along its non-preferred directions where thermal fluctuation, through the inverse temperature β, also plays an important role as it controls the flipping dynamics.However, in the large β limit of the ACM, similar to the small η limit of the DVM, the probability that all particles moving in a particular direction will flip their orientation to another direction is also very small.But, unlike DVM, an ACM particle can move along a direction different than its orientation angle.This crucial difference in the hopping dynamics, as we will see, plays an essential role in the steady-state pattern formation of the DVM at low η and q.
DVM control parameters are the average particle density ρ 0 = N/L x L y , noise strength η, particle velocity v 0 = 0.5 (unless mentioned otherwise), and the measure of anisotropy q.According to Eq. (1) a large q signifies weak anisotropy while a small q signifies strong anisotropy.
We performed numerical simulations of the stochastic process with parallel updates.The initial configuration is prepared homogeneously by assigning random orientations and positions to the particles as defined in Eq. ( 1) and Eq. ( 2), respectively.After the initialization, we let the system evolve under various control parameters for t eq to reach the steady state followed by measurements of various quantities until the maximum simulation time t max .t eq and t max are functions of system size, η, and q.In our simulation, the maximum system size considered is 1024 2 and we have observed that for this square domain, t eq = 10 5 is sufficient for the system to reach the steady state irrespective of η and q.So we take t eq = 10 5 as the steady state time and perform steady state analysis up to a maximum time t max = 10 6 .the typical non-equilibrium steady-state configurations of the DVM in Fig. 2 for q = 9 and density ρ 0 = 2.

Collective
The system exhibits phases of polar ordered liquid (a), liquid-gas coexistence (b-c), and disordered gas (d) as noise strength η is increased from 0.1 to 0.5.The phasecoexistence region is characterized by a low-noise crosssea phase (b) and a high-noise band phase (c).The cross-sea phase has particle density much higher at the crossing points than anywhere else and has recently been reported as the fourth phase of the VM [26].We notice that such patterns assemble spontaneously without an external drive at certain parameter values.Conversely, the band phase is a collection of high-density bands moving parallelly along a specific direction at a constant speed.Polar flocks [the homogeneous ordered liquid phase shown in Fig. 2(a)] can be observed in a large class of active matter systems and have been considered robust to fluctuations.But recent studies have argued that liquid polar flocks are metastable to the presence of a small obstacle [31] or to the nucleation of an oppositephase droplet [32] and triggers counter-propagating dense clusters leading to the reversal of the liquid flow.In light of these observations, in the subsequent part of this paper, we will investigate the stability of the DVM liquid phase.Fig. 2 suggests that the system exists in distinct phases, which we will characterize next.The non-equilibrium steady-state behavior of the DVM is illustrated by representative late-stage snapshots as a As a function of η and q, we observe six distinct self-organized patterns: cluster (η = 0.1, q = 4), macrophase (η = 0.3, q = 4), microphase (η = 0.4, q = 6 → 16), cross-sea (η = 0.3, q = 8 → 10), ordered liquid (η = 0.2, q = 6 → 16), and disordered gas (η = 0.5, q = 4 → 16).
function of noise strength η and anisotropy parameter q (see Fig. 3).We observe six distinct self-organized patterns in the (q, η) plane which completely describe the DVM.At low noise and q, we observe a locally ordered cluster phase.Although each cluster is highly dense and polar, the system as a whole does not possess any net polarization (see Appendix A).This cluster phase, which has not been observed earlier in any flocking models, appears only for small q and η when the system is strongly discretized and the impact of fluctuation is insignificant.The appearance of a cluster phase in the q = 4 DVM for small (q, η) and high density instead of a polar ordered phase similar to the 4-state APM [29,30] and ACM [27,28] can be attributed to the absence of transverse fluctuations through hopping.In DVM, particle movement is solely controlled by the orientation θ similar to the VM [see Eq. ( 2)].Therefore, q = 2 DVM only manifests one-dimensional movement (along the xaxis) of high-density clusters of self-propelled particles having orientations θ = 0 and θ = π.We observe that these clusters never coalesce due to the lack of transverse fluctuations and thus never form a band or polar liquid phase.Similar observations are made when the constant transverse diffusion is switched off in the twodimensional active Ising model (AIM) [33], although the one-dimensional AIM [34] exhibits a flocking state consisting of a single dense ordered aggregate.Unsurprisingly, analogous to q = 2 and q = 4, for small η, we observe a cluster phase also for the q = 3 DVM.We do not observe any band formation similar to q = 4 even when noise is increased.For larger noise strengths, cluster size reduces, and the system exhibits a disordered gas phase.For fixed noise, an increase in density only increases the cluster size without changing the system morphology.This signifies that diffusion along the nonpreferred hopping directions plays a crucial role in forming large liquid domains.For APM and ACM at large β and small q and DVM with small η and q, the probability of transverse flipping is very small since fluctuations are weak.However, the nonzero finite hopping rates along the unbiased directions facilitate the formation of large liquid domains in the APM and ACM; whereas the absence of unbiased hopping gives rise to a cluster phase in the DVM.Thus, in principle, the interplay of q and η determines the fate of the cluster phase.If q is small but fluctuation is large, the rigid cluster phase can relax and form a large ordered domain (see the snapshot for q = 6 and η = 0.2).On the contrary, if η is small but q is large, the weak anisotropy helps the cluster phase merge into a large ordered domain (see the snapshot for q = 10 and η = 0.1).This also explains why the cluster size increases with q for a fixed η (see the snapshots for η = 0.1).
Beyond the cluster phase, for η = 0.2, we observe the emergence of the polar ordered liquid phase (q ⩾ 6).For intermediate noise strength (η = 0.3 − 0.4), a transition is observed from macrophase separation (q = 4, a single liquid band moving through the gaseous background) to microphase separation (multiple bands moving parallelly or in a cross-sea pattern through the gas phase) in the coexistence region where the number of bands increases with q.This is a consequence of having more particle orientations allowed through q.The cross-sea phase, where interactions become more intense due to the characteristics of the band structure, appears between the polar liquid phase and the parallel band state [35] for a fixed q.This phase is not simply a superposition of waves of inclined bands, but an independent self-organized complex pattern with an inherently selected crossing angle.It is worth noting that a single cross-sea pattern typically involves the crossing of at least two bands, where the crossing angle is approximately ∼ π/4 [26].We have observed that this crossing angle remains consistent regardless of q.For a fixed η = 0.3, as q is increased from q = 7 to q = 16, we observe that the self-organized patterns change from bands (q = 7) to cross-sea pattern (q = 8, 9, 10) followed by a polar liquid phase (q = 16).
For the same values of the control parameters, the VM (q → ∞) exhibits features similar to the q = 16 DVM where the phase point on the (η, ρ 0 ) diagram almost lies on the liquid binodal [23].This happens because the anisotropy becomes weaker with q and, for q ⩾ 16, the characteristic of the system becomes similar to the VM.The VM does exhibit a cross-sea phase but at a different parameter regime [26].For a large noise, e.g., η = 0.5, we observe a disordered gas phase irrespective of q.
Since the origin of different phases in the DVM depends on the spatial anisotropy, in Appendix B, we show that the cluster phase is ubiquitous for small (q, η) and present a cluster size analysis of the q = 4 DVM for varying noise.For finite system size, we found that the steadystate DVM exhibits bistability of the cross-sea phase and the band phase a within a range of the noise amplitude which results in a hysteresis (See Appendix C).In Appendix C, we further discussed the stability of the crosssea phase.See also Appendix D and Appendix E for more discussions on the DVM phases with spatial anisotropy (rectangular domain).
Based on the above discussion, we can quantify the DVM phase diagram.We observe four different phases in Fig. 3 [26]: the ordered phase, the cross-sea phase, the band phase, and the disordered gas phase.It is, however, challenging to define the phase boundaries by visual inspection.Previously, phase-separated density profiles were used to identify binodals that delimit the gas and liquid phases from the co-existence region [23].Yet, this technique can not differentiate between the distinct selforganized patterns we observe in the coexistence region (i.e.macrophase separation, microphase separation, and cross-sea) of the DVM.One might also want to distinguish different structures in terms of the global polar order parameter or the global magnetization defined as: The polar order parameter can be used to study the transition between the disorder gas phase and the band states but can not differentiate between the band and the crosssea states [26].The Binder cumulant constructed from this order parameter shows only the transition between the disordered phase and the band state [26].The polar order parameter m is also insensitive around the cross-sea state.Thus, for precise quantification of different phases and their boundaries, we compute the structural order parameter C 2 [26,35,36] as follows: where G 2 (r 1 , r 2 ) = P 2 (r 1 , r 2 )−P 1 (r 1 )P 2 (r 2 ) for one-and two-particle probability density functions P 1 and P 2 , Θ is the Heaviside step function, and R is the distance around an arbitrary fixed point in space.For macroscopically isotropic systems, G 2 can be expressed in terms of the pair correlation function g(r) which physically signifies the probability of finding a particle at a distance r relative to that of a given reference particle and provides a statistical description of the local packing and particle density of the system.It has been shown [26,36] that C 2 performs better than m in capturing the structural change and the Binder cumulant of C 2 is also more efficient in distinguishing different features than the Binder cumulant of the order parameter m.In practice, to compute C 2 , we take all pairs of particles, draw circles of radius R around them, and calculate the overlap area of the two circles.The overlap area is given by , where d is the distance between the centers of the circles.
The dependence of the time-averaged structural order parameter ⟨C 2 ⟩ with respect to η for different q is shown in Fig. 4. The value of ⟨C 2 ⟩ is the lowest for disorder gas and becomes maximum when particles are clustered.For the band phase, ⟨C 2 ⟩ values for macrophase separation are larger than the microphase separation and cross-sea phase.Between microphase separation and the cross-sea phase, ⟨C 2 ⟩ cross > ⟨C 2 ⟩ micro .Although a change in ⟨C 2 ⟩ is not very significant in these two phases, it is still a better candidate for distinguishing the cross-sea from the microphase separation than the traditional polar order parameter.In Fig. 5, we plot the η − q phase diagram by computing ⟨C 2 ⟩ for the six different phases where in the coexistence region, ⟨C 2 ⟩ macro > ⟨C 2 ⟩ cross > ⟨C 2 ⟩ micro .This phase diagram complements Fig. 3, which has been constructed using the density field and depicts the nature of phase separation in the DVM as a function of q.The phase diagram of Fig. 5 has been constructed on a large square domain of dimension 1024 2 .For this system size, one can expect that the finite size effect is relatively small and Fig. 5 will remain qualitatively the same in the thermodynamic limit.As discussed in Appendix C, the DVM might show the bistability of two different coexistence states at the steady state for a particular noise amplitude.However, we believe the system will evolve to a definite steady-state phase at the thermodynamic limit Now, the nature of the coexistence region for any finite, large q has been a subject of discussion in the context of q-state ACM [27,28].It was argued in Ref. [27] that spatial anisotropy plays a crucial role in determining the macro/microphase separation of the coexistence region in the ACM and one should observe a macrophase separation of the coexistence region for a finite q beyond a characteristic length scale which diverges for large q.ACM with a different set of dynamical rules than Ref. [27] was studied in Ref. [28] where the flocking transition in the ACM was argued to be equivalent to the VM at large q.The q-state DVM is governed by a completely differ- ent set of microscopic rules than both the models of the ACM [27,28] and we will therefore investigate next the impact of microscopic rules on the DVM steady-state as a function of q.We plan to do this through the analysis of number fluctuations, the pinning-unpinning property of the system's global order, and the structure factor manifesting the correlation of polarization.
Number fluctuations.-It is known that fluctuations play an essential role in selecting the phaseseparated patterns in flocking models.In the AIM [33], the density fluctuation ∆n 2 = ⟨n 2 ⟩ − ⟨n⟩ 2 (where n is the number of particles in a box of size < L) in the liquid phase was found to be normal (∆n 2 = n) and the AIM coexistence region shows a macrophase separation with a single liquid domain moving on a gaseous background.In the VM [23], on the contrary, giant density fluctuations (∆n 2 = n 1.6 ) are observed which break large liquid domains and prevent the bands from coarsening further.This results in a microphase-separated coexistence region.Giant density fluctuations of the homogeneous ordered phase have also been evaluated in other flocking models [24,[37][38][39][40][41][42].In the q-state ACM [27,28], which can be thought of as a bridge between the discrete AIM and the continuous VM, density fluctuations show a transition from normal to giant fluctuation as q increases and can be explained as a transition from AIM physics to VM physics.As the DVM also exhibits a transition from macrophase to microphase (and cross-sea) separation of the coexistence region as q increases, it is thus useful to investigate the density fluctuation in the DVM.
We show the number fluctuations ∆n 2 versus average particle number ⟨n⟩ in Fig. 6, computed in the ordered liquid phase of the DVM for various q.n(ℓ) is the number of particles in boxes of different sizes ℓ included in a 300 2 domain (with ℓ ⩽ 150), with ⟨n⟩ = ρ 0 ℓ 2 .As shown in Ta- ble I, the number fluctuation behaves like ⟨n⟩ ξ with the fluctuation exponent ξ increasing with q, from ξ ≃ 1.22 for q = 4 to ξ ≃ 1.64 for large q.This transition from normal fluctuations for small q to giant fluctuations for larger q was also observed in the ACM [27,28] although ξ for q = 4 and 5 are moderately larger in the DVM than the ACM.The increase in density fluctuations with q can be attributed to the fact that, for large q, particles have more rotational degrees of freedom due to the weak anisotropy and therefore more directional freedom to propel.The existence of giant number fluctuations (GNF) and its connection with microphase separation in the VM was hypothesized in Ref. [23].It was argued that GNF (ξ ≃ 1.6) breaks bulk liquid domains and produces a smectic-like microphase separation in the coexistence regime whereas the system stabilizes in the bulk phase when the density fluctuations are normal (ξ ≃ 1) [33].Using the same logic for ACM [28], it was argued that GNF is responsible for microphase separation in the coexistence regime for q ⩾ 8, although a direct relation between the existence of GNF in the ordered phase and the microphase separation in the coexistence phase is still ambiguous.Nonetheless, if we compare ξ in Table I and the snapshots in Fig. 3, we observe a correspondence between the fluctuation exponent and the pattern formation in the coexistence region of DVM.For q < 8, although the difference in ξ(q) is small, the exponents change continuously with varying q, similar to the active clock model [28].
One should also consider the finite size effect on the fluctuation exponents as discussed in Ref. [28].In Fig. 6, the data can be fitted to two different power-law regimes (the extracted exponents depend on the interval along the x-axis to which the fits are restricted): (a) ξ tabulated in Table I in the interval [10 3 , 5 × 10 4 ] and (b) a smaller ξ in the interval [5 × 10 4 , 10 5 ].Around the second interval, the plot shows a "saturation" because ξ must decrease with increasing ⟨n⟩ due to the finite-size cut-off at ⟨n⟩ = N = ρ 0 L 2 , where ∆n 2 vanishes.
Pinned property of the order parameter.-Theglobal order parameter defined in Eq. ( 6) quantifies the overall ordering of the particles in the system.In the ACM [27], at finite size, the direction of the global polar order Φ ≡ arg⟨m⟩ exhibits distinct behaviors in the liquid phase depending on the value of q: it displays AIMlike properties (pinned along an angle) for small q and VM-like behavior (unpinned over time) for large q.However, for a particular q, an unpinning to pinning transition is observed as the system size is increased which indicates the sensitivity of the ACM [27]  on system size and a possible transition from micro to macro phase separation of the coexistence region beyond a length scale.We therefore are interested in analyzing this property for the DVM.
In Fig. 7(a), we show the time evolution of Φ(t) in the DVM liquid phase for varying q.Similarly to the observation made in the ACM [27], Φ(t) begins wandering slowly with q and becomes an unpinned variable of t for large q.In other words, for weak anisotropy, the global ordering does not remain constrained to a specific orientation.While, for small q, Φ(t) remains pinned and tends to exhibit a stable global ordering.Microscopically, this refers to a picture in which, at large q, a proportional number of degrees of freedom allows the particles to choose between adjacent directions facilitated by fluctuation, while it is not the obvious choice for particles in small q that require a significantly larger jump to switch directions.Translating this to the global polar order parameter and comparing Fig. 7(a) with Fig. 6 we propose that GNF corresponds to the unpinned behavior of Φ(t).Likewise, the unpinning nature of the direction of the global polar order is a characteristic of microphase separation.
In addition to q, the finite system size also affects the evolution of Φ(t) (as was also shown in Ref. [27]), which is shown in Fig. 7(b) for q = 9.Similar to the ACM [27], we observe a transition from unpinned behavior to pinned behavior in Φ(t) as the system size increases.In larger systems with polar order, a particle interacts with more particles in the neighborhood and correlates over longer distances.Higher connectivity promotes stronger alignment and cooperative motion among the particles.As a result, the direction of the global order becomes more pronounced and persistent in larger systems, leading to the pinned state.Fig. 7(b) further indicates that if Φ(t) is pinned for L = 300, it must also be pinned for L = 1024.This is inconsistent with the microphase separation and the cross-sea patterns observed in the coexistence region of the q = 9 DVM (Fig. 3).Evidently, the correlation proposed earlier between the pinned property of the system's ordered liquid phase and the system morphology observed in the coexistence region (macro/micro/crosssea) is not conclusive in the DVM (see Appendix F for more details).Structure factor.-Tounderstand the system's behavior for different q, investigation of the magnetization correlation function which analyzes how the particle's orientations are correlated as a function of q is a good measure.In the ACM [27], it was observed that the structure factor (Fourier transform of the orientation correlation converges to finite values at large wavelengths beyond a crossover length scale.This length scale is a function of q which signifies the crucial impact of the system size on anisotropy. To explore the correlation among particle polarizations, we consider the transverse magnetization structure factor S ⊥ (k) = ⟨m ⊥ (k)m ⊥ (−k)⟩ [27] against wavelength k and plot it in Fig. 8.The structure factor has been calculated in the liquid phase on a 300 2 domain for (a) various q values (the same behavior is observed for the structure factor of the density field) and (b) for a fixed q = 9 but for various system sizes.The results presented in Fig. 8(a) show that for small q, the structure factor S(k) converges to finite values as the wave vector k → 0. This convergence indicates an AIM-like behavior or a macrophase separation of the coexistence region [27].However, one can notice that this convergence is achieved only beyond certain length scales and these length scales are functions of q.The structure factor captures the correlations between particle orientations and therefore the magnitude of ordering at different length scales.A saturation of S(k) for small q when k approaches zero thus signifies a strongly correlated liquid domain whereas for large q, due to weak anisotropy or more allowed orientations for ordering, particles inside the liquid domain are not as strongly correlated as for small q.Fig. 8(b) shows S(k) for several system sizes and manifests that with larger system sizes (L ≥ 300), S(k) tends to converge to a finite value when k → 0. Our earlier argument that stronger interactions between particles (with larger L) promote robust ordering also applies here.Comparing Fig. 8 with Fig. 7, we conclude that the pinning behavior (unpinning behavior) of Φ(t) and the saturation of S(k) (algebraic scaling of S(k)) compliment each other and convey the same physics.It was argued in Ref. [27] considering the pinning properties of the order parameter and behavior of the structure factor in the liquid phase that for large q ACM, VM behavior (microphase separation) will be observed only up to large finite sizes.But the asymptotic large length scale behavior will be AIM-like (macrophase separation) where the length scale diverges with q as exp(q 2 ).Ref. [27] also showed that in the phase coexistence region of the ACM, a microphase to macrophase transition occurs when the linear system size increases along the transverse direction at fixed q, which we do not observe (in the DVM, multiple bands do not merge to a single band when L y is increased for a fixed L x ).We observe a similar behavior of Φ(t) and S(k) as Ref. [27] but can not draw a conclusive correspondence between the large length-scale liquid phase behavior of Φ(t) and S(k) to the phase-coexistence behavior of the DVM as shown in Fig. 3.
In Fig. 3, the snapshots on a large square domain (without any spatial anisotropy) show the existence of a microphase separation and cross-sea patterns (for which one needs at least two bands) of the phase-coexistence region for large q values.The number fluctuation plotted in Fig. 6 corroborates this observation by exhibiting GNF for those large q values.The large length scale asymptotic behavior (for large q) of the direction of global order Φ(t) (Fig. 7) and the structure factor S(k) (Fig. 8) in the ordered liquid phase respectively shows a pinned behav- ior and saturation for k ∥ = 0 which as argued in Ref. [27] signifies an AIM phenomenology.However, our numerical investigation of the DVM does not show a cross-over from micro-to macrophase separation for higher q values as observed in the ACM [27] although matches with the observation of another model of ACM [28] with different dynamical rules.Therefore, we argue that the impact of dynamical rules governing flipping and hopping in a flocking model has a significant influence over the system dynamics.Order parameter distribution.-InVicsek-like models, where particles are active and move with a constant velocity, the ordered state exhibits a true longrange order (LRO) in two dimensions because of a spontaneous symmetry breaking due to the out-of-equilibrium nature of the models.In Fig. 9, we show the time-and ensemble-averaged distribution of the order parameter m = (m x , m y ) for increasing system sizes, where m x = 1 N N i=1 cos θ i and m y = 1 N N i=1 sin θ i .In Fig. 9(a) and Fig. 9(b), ring-like distributions (unpinned orientations) are the characteristic of the quasi long-range ordered (QLRO) phase.However, this occurs due to the finitesize effect and is similar to the impact of finite-size on Φ(t) and S(k).We recover the LRO for larger system sizes (L = 200 and L = 300) where the distributions display nine distinct isolated spots (pinned orientations) that correspond to the 9-fold degeneracy of the ordered liquid phase, each spot having equal probability.One can expect that the finite size effect will be much weaker for L = 1024 and the spread of the distribution in the LRO phase around the allowed ordering angles will also be more precise.
The DVM for large q exhibits a QLRO phase when v 0 = 0 (see Appendix G).We argue that DVM with immobile particles reduces to the two-dimensional q-state clock model (with a quenched bond-disorder as only particles within a fixed distance interact) which approaches the XY model for large q with vanishing LRO regime [43].For v 0 > 0 and a fixed L, as flocking directions increase with q, we again observe ring-like distributions for large q (see Appendix G) but beyond a length scale which is proportional to q, the order parameter distributions for large q show q isolated spots characteristic of the LRO phase.This is similar to the unpinned to the pinned transition of Φ(t) and convergence of S(k) to a finite value at k → 0 for large q values as the system size is increased.In the VM (q → ∞), even for a large L, the order parameter distribution shows a ring-like structure because of the continuous symmetry.
Stability of the ordered liquid phase.-Asdiscussed in the context of Fig. 2, here we present a brief analysis of the stability of the DVM ordered liquid phase by inserting a small high-density counter-propagating liquid droplet.Initially, the average direction of the particles in the polar liquid phase is aligned along the direction Φ liq = π, and the average particle orientation of the liquid droplet is Φ d = 0.The radius of the droplet is r d and density is ρ d 0 and it is inserted in an ordered phase of density ρ 0 (ρ d 0 ≫ ρ 0 ).We take several q values and calculate the probability (P rev ) that the droplet grows against the main order and reverses the ordered phase as a function of η, r d , and ρ d 0 .The simulation protocol follows Ref. [32], where the initial configuration is prepared by particles with θ = π and we let the system evolve up to a certain time t to reach the steady state.The system retains the average global polarization in the direction θ = π.Then, a circular region of radius r d centered at (L/2, L/2) is selected and an additional ∆N = (ρ d 0 − ρ 0 )πr 2 d number of particles are added to make a high-density circular droplet.Finally, the orientation of all the particles within the droplet is changed to θ = 0.
In Fig. 10, we study the fate of polar flocks in q = 6 and 16-state DVM by introducing a small high-density counter-propagating droplet against the initial polar ordered liquid phase of the main flow and observe the subsequent time evolution.One should note here that the perturbation through the droplet is very small i.e. the ratio of droplet diameter to the linear length of the simulation box is ∼ 10 −2 (r d = 5, L = 1024).We observe, similar to Ref. [31], that the droplet grows with time leaving behind a dilute region (t = 10 2 ) and adds more and more particles as it moves along forming a principal dense, curved band (followed by several other curved bands) that invades the whole system ballistically (t = 10 3 and t = 10 4 ).In the final stage, this principal curved band connects itself over the system boundaries and widens until a steady state liquid phase of a different Φ liq emerges (at t = 10 5 ), signifying the metastability of the DVM liquid phase.The time-evolution is similar for both small and large q, which signifies that both discrete and continuous-symmetry flocks are metastable.
The growth pattern of the DVM droplet is similar to the VM [31] but distinct from the AIM [32].In the DVM, after its introduction, the droplet front interacts with the liquid particles outside and creates a curved band of particles having several different orientations (impact of q).In this process, the droplet seizes to exist and it is this high-density curved band that destroys the initial flow.
In AIM [32], the droplet grows along the direction transverse to the propagation (due to the constant transverse diffusion) creating a comet-like trail of disordered particles that can not be observed in the DVM.
Fig. 11 quantifies P rev , the probability of reversing the main flow upon the introduction of a given droplet, for several control parameters.For each set of control parameters, we have taken 20 independent realizations to calculate P rev .Akin to Ref. [31], we observe that the noise strength η has a strong influence on the reversal dynamics and P rev increases from 0 to 1 as η is increased [Fig.11(a)].This is because for small η, the ordered phase is very stable, and thus, the counter-propagating dense bands find it difficult to reverse its flow.For large η, fluctuations are stronger, and therefore, the probability of reversal increases.Fig. 11(a) also exhibits that the transitional value of η (P rev = 1  2 ) above which a droplet triggers a reversal is a decreasing function of r d although there is a critical η (η ∼ 0. r d combined with large ρ d 0 are shown to facilitate the reversal of the initial liquid phase.The droplet-induced reversal of the liquid flow is found independent of q where P rev is found to behave similarly for each q under certain values of η [Fig.11(c)].
In addition, a study of the metastability of the DVM liquid phase on rectangular geometry (L x = 800, L y = 100), produces a similar outcome.One can also consider droplet movement in other directions than opposite to the global flow, e.g.transversely propagating droplets with Φ d = π/2 for q = 8 or Φ d = 2π/3 for q = 6 and perform a similar analysis to check whether the liquid phase is susceptible to droplets irrespective of their propagation direction.

IV. DISCUSSION
Our study motivated by the active clock model [27,28], considers a true q-state discrete version of the Vicsek model where q defines the strength of orientation anisotropy.At small q, the system is highly anisotropic, which, however, vanishes in the limit q → ∞ when we recover the Vicsek model.The DVM shows qualitatively similar features as the ACM [28] for intermediate noise strength η where a transition from macrophase to microphase separation is observed in the coexistence region as q increases.But for small q and η, the liquid phase appearing in the ACM at low temperatures is replaced in the DVM by a cluster phase.The cluster phase consists of multiple clusters with different polarization (see Fig. 12) which does not exhibit a long-range order.The clusters grow and merge with increasing q leading to a homogeneous ordered phase at large q.For small q, a long-range ordered phase can be achieved by increasing the noise strength.At low noise and small q, the flipping probability is very small, and in addition, transverse fluctuations through hopping are also absent.The combined influence of these factors results in clusters failing to grow continuously at small η and q, preventing the system from reaching a homogeneous liquid state.Therefore, the polar ordered phase which is ubiquitous in discretized flocking models (such as the AIM, APM, and the small q limit of the ACM) for small noise has been replaced by a cluster phase in the DVM and is a consequence of the strong anisotropy through q.The DVM recovers the ordered phase for small noise at the large q limit, which signifies the DVM to be in the same class as the VM (q → ∞ DVM).
The self-organized patterns in the coexistence region of the discretized VM indicate a transition from AIMlike patterns to VM-like patterns as anisotropy becomes weaker.This observation is corroborated by the giant density fluctuations for large q.However, the large length scale behavior of the direction of global order Φ(t), the structure factor S(k), and the order parameter distribution in the liquid phase do not correspond with the phasecoexistence patterns of the large q DVM.The DVM without any spatial anisotropy and at large length scales shows a transition from a macrophase-separated coexistence region to a coexistence region having microphase or cross-sea pattern with increasing q.This is similar to the observation made in the q-state active clock model [28] which also approaches the VM at the large q limit.
We also find that the DVM liquid phase is susceptible to perturbation applied through a counter-propagating droplet.The liquid phase reorients and propagates along the direction of the droplet.The reversal dynamics is significantly impacted by the noise strength η [31] but remains independent of q.The stability of the high-density flocking ordered phase at low noise is still an open problem and will be addressed in a subsequent study [44].
As a final remark, we add that the rotational flexibility of the particles and microscopic details of the dynamical rules can significantly impact the macroscopic properties of the ordered phase.It would be interesting to compare the model predictions of the DVM with suitable experiments where anisotropy in the particle orientation may be controlled by a finite number of motility directions.steady state remains in an ordered liquid phase (similar to the steady state behavior of the 4-state APM or ACM at low temperature) signifying that the fluctuation is weak to alter the broken symmetry phase into a cluster phase.However, with an increase in the noise (η = 0.2), the steady-state cluster phase appears again.This also suggests that for DVM with weak anisotropy and fluctuations, the cluster phase is the non-equilibrium steadystate and the well-known small η or large β ordered liquid phase can only be achieved by taking a strongly polarized ordered initial condition.It is worth noting that the number of clusters in (c) is higher than in (b) due to more relaxation via the noise.A further increase of the noise will lead to a single large ordered domain as shown in the phase diagram of Fig. 18.
Appendix B: Cluster size analysis for q = 4 The cluster size analysis for q = 4 is shown in Fig. 13 for a rectangular domain of size 800 × 100.We use a box-counting method to find a cluster which is described below: If L x and L y are the linear sizes of a rectangular domain, L x × L y is the total number of boxes we consider.For a box i, n i is the number of particles in that box and c i is the cluster label of the box.If the box does not contain any particles c i = 0 i.e. it is not part of any cluster.Then we use the Depth-first search (DFS) algorithm to find the connected boxes that are not void of particles and label them as a single cluster.For each cluster label c j , we calculate the size of the cluster S ci as following: where δ ci,cj is the Kronecker delta function that equals 1 if c i = c j , and 0 otherwise.Then we calculate the cluster size probability distribution denoted as P (S): Fig. 13 illustrates that for small noise, the probability of larger cluster formation is high.This is because reduced fluctuation in the system facilitates the formation of high-density clusters.However, with noise (η = 0.3), the cluster phase vanishes and the q = 4 DVM exhibits a macrophase separation, resulting in a decrease in the probability of obtaining large clusters.At large noise (η ⩾ 0.4), the system becomes disordered, leading to a lesser probability of formation of large clusters.We would also like to mention that for a fixed noise, P (S) versus S for various orientations (θ = 0, π/2, π, 3π/2) shows almost identical distributions signifying no preference in orientation in the cluster formation.In this section, we analyze the stability of the different DVM steady-state phases.Figure 14 (a) illustrates the stability of the cross-sea phase where we ask if one starts with a band state in the cross-sea regime, will the system evolve to a cross-sea pattern?The figure shows that the system initially configured in a highdensity single band (visualized as a vertical stripe of aligned particles), evolves into a cross-sea phase after a long time (t = 131072).A similar observation is made with microphase-separated initial bands in Fig. 14 (b), suggesting that the steady-state phases of the DVM are independent of the initial configurations.However, as observed in Ref. [35] for the VM, in the DVM also, phase boundaries between the two neighboring phases, the cross-sea phase and the microphase, are not always distinct.For a finite system size, bistability between these two neighboring phases can be observed, as shown in Fig. 15.At low noise (η = 0.1), the system exhibits a homogeneous high-density ordered phase, where particles move mainly in a uniform direction.With a small increase of the noise (η = 0.2), the cross-sea phase starts to appear, and at high noise amplitude (η = 0.45), the system transitions to a disordered phase, where particles show random orientations.Interestingly, in the intermediate noise amplitudes (∼ η = 0.28 − 0.35), the system exhibits bistability where both cross-sea and microphase features can be found depending upon the initial conditions.In the VM, bistability was observed between the ordered phase and cross-sea phase and also between the microphase and disordered phase [35].This bistability arises due to the fluctuation at relatively higher noise amplitudes and finite system size.We therefore expect the system to evolve to a specific stable steady-state phase for a particular noise amplitude at the thermodynamic limit (L → ∞).The bistability between the cross-sea phase and the microphase as a function of noise amplitude is further demonstrated by a hysteresis loop in Fig. 16.The hysteresis loop is obtained by plotting ϕ cs as a function of noise η keeping the density (ρ 0 ) and particle velocity (v 0 ) constant.Where ϕ cs is either 1 or 0, corresponding to cross-sea and microphase, respectively.In the low noise limit, we evolve an initially disordered configuration into a steady state coexistence phase by slowly varying the noise amplitude by (∆η = 0.02).The cross-sea phase demonstrated by ϕ cs = 1 starts ∼ η = 0.2 and is retained up to a large noise (red curve).At ∼ η = 0.35, the crosssea switches to the microphase and remains steady for up to ∼ η = 0.4.While reducing the noise from this high value (blue curve), the microphase bands (ϕ cs = 0) sustain for a wide range of η stretching significantly below the cross-sea to microphase transition point (∼ η = 0.35).Finally, the microphase switches back to the cross-sea phase around (∼ η = 0.28).As shown in Fig. 16, both the transitions (at small and large η) between these two phases appear to be discontinuous.
Appendix D: Impact of spatial anisotropy: steady-states for rectangular domain In Fig. 17, we investigate how spatial anisotropy influences the non-equilibrium steady-state behavior of the DVM if we switch from a square domain to a rectangular domain by analyzing late-stage representative snapshots as a function of noise strength (η) and discretization parameter (q).When noise is low and q is small, we observe the emergence of a locally ordered cluster phase similar to our finding for the square domain.However, when q is small and fluctuations are pronounced, those cluster phases relax and transform into a larger, organized domain (see the snapshot for q = 4 and η = 0.3) akin to the Fig. 3. Conversely, when η is small and q is large, weak anisotropy facilitates the merger of the cluster phase into a larger, well-organized domain (see the snapshot for q ⩾ 8 and η = 0.1) which is observed at q = 10 in the absence of spatial anisotropy (see Fig. 3).An increase in the fluctuation for small q might exhibit multiple bands but those are connected by the periodic boundary conditions and should be considered a single band (see the snapshot for q = 5 and η = 0.2, 0.3).However, no cross-sea phase is observed with the rectangular geometry, which is seen in Ref [26] probably due to larger particle velocity (v 0 = 1).An increase in anisotropy q also increases the no. of bands in the coexistence region at intermediate noise (η = 0.4) because with q, the density fluctuation increases along with the magnetization fluctuation which prompts the breaking of large domains.
Appendix E: (η − ρ0) phase diagram of the DVM In Fig. 18, we plot the noise-density (η − ρ 0 ) phase diagrams computed on a rectangular domain of 800 × 100 for various q.The clustering phase is very prominent for small q values (q = 4, 5) and exists for high densities.As q is increased, the conventional ordered liquid phase appears at high densities and for q ⩾ 8, the cluster phase disappears and we notice the emergence of the typical η − ρ 0 diagram observed for Vicsek-like systems [Fig.18(d)].These Vicsek-like phases are characterized by a more global alignment of the particles, leading to coherent motion and the absence of distinct clusters or bands.In this regime, the behavior of the system is predominantly governed by the alignment interactions between the particles rather than the specific value of q.It shows clear boundaries between different phases based on varying values of q and η.At small q and intermediate densities, a macrophase separation similar to AIM is observed.However, as q exceeds a threshold (e.g., q ⩾ 8), the system transitions to a Vicsek-like phase separation in the coexistence region.As a function of η and q, we observe five distinct self-organized patterns: cluster (η = 0.1, q = 4), macrophase separation (η = 0.3, q = 4), microphase separation (η = 0.4, q = 7 → 16), ordered liquid (η = 0.2, q = 7 → 16), and disordered gas (η = 0.5, q = 4 → 16).contrast, for v 0 > 0 (see Fig. 21), the order parameter distribution exhibits LRO through isolated points of phase ordering as activity facilitates the broken symmetry phase.For q > 7, a comprehensible LRO phase is observed at a large length scale limit as shown in Fig. 9.
To understand better the nature of ordering of the DVM liquid phase, we show the order parameter ⟨m⟩ against increasing system size L for several q in Fig. 22.The data presented are averaged over time and several initial configurations.We note that, ⟨m⟩ remains inde-pendent of the system size L (actually, m scales with L as m ∼ L −λ , decays much slower than a power law) for all q for v 0 = 0.5 [Fig.22(a)] signifying LRO.As a result, the liquid phase of the constant-speed DVM is LRO and the direction of the order parameter exhibits a pinned behavior.For v 0 = 0, shown in Fig. 22 is expected to algebraically decay to zero for L → ∞ and this effect is more pronounced for larger q because, for large q, the v 0 = 0 DVM approaches the two-dimensional XY model where the low-temperature phase is QLRO.

FIG. 1 .
FIG. 1. (Color online) Schematic of the DVM for q = 4. (a) Four allowed orientations for the particle are 0, π/2, π, and 3π/2 with θ t i = 0.The circular neighborhood (of radius R) of particle i contains neighboring particles to calculate θ t+∆t i .(b) The new orientation proposed by the Vicsek rule (blue dotted arrow) will be projected either along θ1 = 0 or θ2 = π/2 probabilistically.
22 1.28 1.41 1.48 1.67 1.64 1.64 1.64 TABLE I. (Color online) Number fluctuation exponents ξ for several values of q, reported from Fig. 6.The typical error on the fluctuation exponents is 0.03.

FIG. 9 .
FIG. 9. (Color online) Order parameter distributions of the q = 9 DVM for the liquid phase.Parameters: η = 0.3, v0 = 0.5 and ρ0 = 6.Ring-like distributions for smaller system sizes (L = 50 and L = 100) in (a) and (b) change to distinct isolated spots for larger system sizes (L = 200 and L = 300) in (c) and (d), which correspond to the ordered phase's q-fold degeneracy.
FIG. 10. (Color online) Time evolution snapshots of the orientation field showing the reversal of the initial ordered liquid phase in (a) q = 6 and (b) q = 16-state DVM following the introduction of a tiny counter-propagating high-density droplet (motion of direction is denoted by black arrows, t = 1) of radius r d = 5 and density ρ d 0 = 10ρ0.Colorbar represents particle orientation field.Parameters: L = 1024, η = 0.4, v0 = 1 and ρ0 = 2.

=
Number of clusters with size S Total number of clusters .(B2) Appendix C: Stability of the DVM steady-state phases

FIG. 18
FIG. 18. (Color online) Phase diagrams on η − ρ0 plane computed on a rectangular domain of 800 × 100 for various q.A macrophase separation similar to AIM in (a) and (b) transforms to a Vicsek-like phase separation (microphase) in the coexistence region in (c) and (d) at intermediate densities.v0 = 0.5.

FIG. 19
FIG. 19.(Color online) Time evolution of the direction of global order Φ in the coexistence region for (a) L = 512 and various q and (b) q = 9 and several values of L. (c-d) Snapshots showing multiple bands moving parallelly constituting a coexistence region for q = 9 and q = 16, respectively.System size L = 512.Parameters: η = 0.3, ρ0 = 1.5, v0 = 0.5.

FIG. 20
FIG. 20. (Color online) Zero-activity limit of the order parameter distributions of the DVM with varying q in the liquid phase.Parameters: L = 100, η = 0.3, and ρ0 = 6.Distinct isolated spots in (a), (b), and (c) indicate LRO while ringlike distributions in (d), (e), and (f) are characteristic of the QLRO phase.