Vectorial active matter on the lattice: polar condensates and nematic filaments

We introduce a novel lattice-gas cellular automaton (LGCA) for compressible vectorial active matter with polar and nematic velocity alignment. Interactions are, by construction, zero-range. For polar alignment, we show the system undergoes a phase transition that promotes aggregation with strong resemblance to the classic zero-range process. We find that above a critical point, the states of a macroscopic fraction of the particles in the system coalesce into the same state, sharing the same position and momentum (polar condensate). For nematic alignment, the system also exhibits condensation, but there exist fundamental differences: a macroscopic fraction of the particles in the system collapses into a filament, where particles possess only two possible momenta. Furthermore, we derive hydrodynamic equations for the active LGCA model to understand the phase transitions and condensation that undergoes the system. We also show that generically the discrete lattice symmetries—e.g. of a square or hexagonal lattice—affect drastically the emergent large-scale properties of on-lattice active systems. The study puts in evidence that aligning active matter on the lattice displays new behavior, including phase transitions to states that share similarities to condensation models.


Introduction
Self-organized patterns of self-propelled entities are found at every scale in biology: from in vitro cytoskeletal extracts [1,2,3], bacterial systems [4,5,6,7,8], and insect swarms [9] to schools of fish [10], herds of mammals [11], and human crowds [12].Collective phenomena are also observed in artificial active systems such as synthetic swimmers [13,14], phoretic colloids [15,16], and active rollers [17,18,19].Active systems in the absence of torques -often referred as scalar active matter -can phase separate [20,21,22].Other complex, self-organized, collective motion patterns exhibit either polar or nematic (orientational) order.Given the vectorial nature of orientational order, active systems sharing such a feature are often referred to as "vectorial active matter" [23].It is widely believed that self-organized pattern formation in vectorial active matter requires particles to interact via some type of velocity alignment Figure 1: LGCA dynamics.Particles reside in velocity channels within lattice sites.During the reorientation step (R), particles within each lattice site are assigned new velocity channels according to the transition probability T ϑ (s).Subsequently, during the migration step (M), particles are deterministically moved to neighboring lattice sites.mechanism [24,25], though alternative mechanisms exist [26], such as combined attraction-repulsion and active forces [27].The best known example of an active system with polar velocity alignment is the Vicsek model [28], where complex, active polar flows spontaneously emerge.A different class of active systems is obtained by changing the velocity alignment in the Vicsek model from polar to nematic [29,30].Self-propelled (point-mass) particles with nematic alignment do not generate polar flows, but self-organize into high-density, nematic bands [30].The theoretical understanding of vectorial active matter has been strongly based on the study of off-lattice Vicsek-like models by means of numerically costly, particle-based simulations [28,30,31,32,33,34] and by the derivation of hydrodynamic equations for these systems [35,36,37,38,39,40,41]; for comprehensive reviews, see [24,25].There is, however, a renewed interest in understanding the behavior of active matter on lattices given the suitability of efficient numerical schemes, and in particular due to the large repertoire of available analytical tools to study lattice models that are expected to help to improve our theoretical understanding of active matter physics [42,43,44,45,46,47].
Here, we introduce an active lattice-gas cellular automaton (LGCA) for compressible (i.e.without volume exclusion) vectorial active matter, including both polar and nematic velocity alignment, where interactions are, by construction, zero-range.We show that the introduced LGCA with polar alignment undergoes a non-equilibrium phase transition that is fundamentally different from the non-local, incompressible active LGCA, where the onset of order is characterized by the presence of travelling polar bands [48] as occurs in the off-lattice Vicsek model [49].In contrast, in the present model (with polar alignment), the emergence of order occurs via a process by which a macroscopic fraction of the particles coalesce into the same state, and thus sharing the same position and momentum, i.e. a polar condensate.For nematic aligment, we find that the system displays a condensation, but with a fundamental difference: a macroscopic fraction of the particles in the system collapses into a nematic filament.Note that recently, condensation behavior was found in scalar active matter in the presence of an external potential by a mechanism known edge-diffusivity [50,51].In scalar active matter -i.e. in self-propelled disks -it has been observed that phase separation into dense and dilute phases can result from volume exclusion interactions among the particles, a phenomenon known as motility-induced phase separation (MIPS); see [52] for a review).Our study shows that in vectorial active matter, condensation behavior is based on velocity alignment and does not require an external potential or external field, even in the absence of volume exclusion interactions.Thus, the emergence of polar and nematic order in the active LGCA is fundamentally different from the onset of order in other models [28,35,31,43,47,44,45,46,42].Similarly, phase separation in the active LGCA, that takes the form of a condensation, is fundamentally different from MIPS [52].At high interaction strengths and low effective temperatures, these models show a great variety of collective behaviors such as homogeneous, polar-ordered states, formation of spatiallyextended clusters, emergence of long-range order, filament-linked clusters, phase separation into extended jammed clusters, the formation of moving bands, and rippling; whereas our active LGCA models consistently show only two phases: a disordered, homogeneous state, and an aligned, highly condensed state in either single lattice sites, or in filaments of a single lattice site in width, depending on the nature of the alignment interaction.Furthermore, the hydrodynamic equations for the introduced LGCA allows us to show that generically the discrete lattice symmetries -at least for square and hexagonal lattices -have a strong impact on the emergent large-scale properties of on-lattice active systems.In summary, our study unveils the existence of condensed phases in vectorial active matter on the lattice, without attractive or excluded volume interactions.

Model definition
We study an hexagonal LGCA, where each node possesses six velocity channels, which can be occupied by a number s ϑ of cells moving in the direction specified by c ϑ .In contrast to standard passive [53] and active [48,54] LGCA that consider that each velocity channel of a node can either be unoccupied or occupied by at most one cell, i.e. s ϑ ∈ {0, 1} -similar to the Pauli exclusion principle -here we assume particles without volume exclusion by letting any velocity channel c ϑ of a node contain an unrestricted number of particles s ϑ ∈ N: i.e. many particle can share the same state (position and momentum).At each discrete time step, as illustrated in Fig. 1, two processes take place: (i) Particles modify their velocity by transitioning to a channel ϑ, according to the transition probability T ϑ : with H ϑ (r, k) defined as: where s(r,k) = [s 1 (r,k), s 2 (r,k), • • • , s 6 (r,k)] is the configuration of node r at time k, which is given by the number of particles s α (r, k) occupying channel α (with α ∈ [1, 2, • • • , 6]), the strength of the velocity alignment is denoted by β ∈ R + , and Z is the node normalization constant.Finally, m = 1 corresponds to polar and m = 2 to nematic velocity alignment.These transition probabilities are the volume exclusion-free analogues of those derived previously [55], which were obtained as discretized, steady state distributions of flocking off-lattice Langevin equation models with pairwise velocity alignment.Note that Eq. 1 is a Maxwell-Boltzmann distribution, where the parameter β is inversely proportional to the effective temperature.Indeed, as seen in [55], β is proportional to the interaction strength, and inversely proportional to the noise amplitude of off-lattice, stochastic differential equation models.Thus, in the limit β → 0 (corresponding to infinitely high temperatures or noise amplitudes), Eq. 1 reduces to a discrete, uniform distribution, with particles exhibiting no preference in their reorientation.On the other hand, in the limit β → ∞ (corresponding to an absolute zero temperature or no noise), T ϑ = 1 for a single orientation ϑ and T ϑ = 0 for all other orientations, resulting in completely deterministic behavior, for m = 1, and two opposite orientations with T ϑ = 1 2 , and T ϑ = 0 for all other orientations, resulting in a one-dimensional, symmetric random walk along a single axis, for m = 2.We stress that Eqs. ( 1) and ( 2) define a purely local velocity alignment process: the reorientation dynamics does not consider the configuration of neighboring nodes.Note that the number of particles at node r at time k is given by n(r, k) = α s α (r,k).This implies that Eq. ( 2) for m = 1 can be expressed as (ii) Particles move spatially according to the following rule: particles in node r at channel ϑ are transported into the same channel into the neighboring node r + c ϑ ; see Fig. 1.It should be noted that lattice-gas cellular automata are mesoscopic, rather than microscopic, models.Therefore, time and space are coarse-grained.At this level of resolution, the reorientation processes that occur continuously in microscopic models are assumed to have reached a steady state at every discrete time step [55], and thus can only move at a fixed speed of one lattice site per time step.

Emergence of polar condensates
For polar alignment (m = 1), we observe, in strong contrast to other polar active systems whose polar ordered phase consists of either a traveling polar band or an spatially homogeneous polar phase [28,35,31], the emergence of polar, moving condensates; see Fig. 2a.In the following, we explain the mechanism that leads to the emergence of such condensates.We make use of a mean-field argument and define the probability of finding a particle in channel ϑ of a certain node at time k as q ϑ (k) = s ϑ (k)/n(k), with n(k) = α∈Θ s α (k) being the number of particles in the node.The node configuration, after reorientation, takes the form: Let us assume a steady state condition: q ϑ (k + 1) = q ϑ (k) for all ϑ.Under this assumption, the transcendental Eq. ( 3) can be easily solved; see Figs. 2b and c.Note that it is possible to recast Eq.( 3) as q ϑ = exp [βnP • c ϑ ] /Z in order to obtain an equation for the polar order parameter of the node: This transcendental equation describes a phase transition.The disorder state, given by ||P|| = 0, and corresponding to q ϑ = 1/6 for all ϑ, is a solution of the above transcendental equation.This solution is stable for a number of particle n below the critical n * value.This is evidenced by the behavior of ||P|| in Fig. 2b.Above n * , q ϑ = 1/6 is not stable and particles accumulate preferentially in one velocity channel and polar order in the node emerges.It is very important to realize that the probability of finding all particles in a node occupying the same velocity channel, for n > n * grows exponentially with the number of particles n in the node.This implies that once all particles in a node occupy the same velocity channel, the probability that single particles transition to other velocity channels decays exponentially with n, as does the probability of not finding the n particles in the same node.This mechanism leads to a condensation of particles.To understand this, assume, without loss of generality, that all particles are initially located in channel 0; then q 0 = e βn /Z = (1+e βn(γ1−1) +• • •+e βn(γ5−1) ) −1 , with γ j = cos(−π j/3) so γ j −1 < 0. The probability that particles change channel is 1 − q 0 = 5 j=1 e βn(γj −1) /[1 + 5 j=1 e βn(γj−1) ], which decays exponentially fast with n, while the probability of finding all particles again in channel 0 after one time step behaves as (1 + e βn(γ1−1) + • • • + e βn(γ5−1) ) −n .Consequently, the probability for a particle in a condensate of mass n ≫ n * to leave the condensate decays exponentially with n; Fig. 2c.This feature is shared by the standard zero-range process (ZRP) [56].It should be mentioned that, in this case, randomness decreases exponentially with particle number n due to the pairwise, additive definition of H ϑ (r, k) (Eq.2), which increases linearly as more aligned particles are present.Conversely, other population models may exhibit similar fluctuations originating from multiplicative noise terms [57], which we do not consider in this model.Particle condensation can be evidenced by looking at the distribution of particles per node P (n, t); Fig. 2d (cf.[58]).The second peak that shifts over time is indicative of an aggregation dynamics.Note that here, it is the alignment mechanism that holds particles together.The emergent condensate moves ballistically in contrast to the processes studied in [56,58].This means that the condensation dynamics is not solely mediated by particle exchange with the bulk, but results from the motion of the condensate: as the condensate moves cross the system collects clusters in its way and grows in size, Figs.2d and e.This is further supported by the fact that the mean condensation time increases with system size (see Fig. S8 of the Supplementary Material), since clusters need to travel longer distances to accumulate particles in distant lattice sites.Note that, as argued in [59], in the absence of velocity correlations among clusters (and only in this scenario), polar order becomes intimately connected with the number and distribution of moving clusters.As expected, in Fig. 2f we observe that as particles coalesce and the mean cluster size grows with time, the global order parameter P g -defined as P g = || r α∈Θ q α (r)c α || also does it.Furthermore, note that the mean cluster size scales with time with an exponent larger than 1/2; cf.[56].Moreover, Fig. 2f shows that the condensate size N c , defined as N c = max {s ϑ (r) : ϑ ∈ Θ, r ∈ L} , and P g exhibit almost identical dependency with β ρ, with ρ = n (r) = N/L 2 , where N is the number of particles in a system of linear size L.An additional observable obtained from simulations was the local polar order, defined as P ℓ = r || α∈Θ q α (r)c α ||/N o , where N o is the number of nodes r ∈ L with at least one particle in any channel.Thus, P ℓ is the average of the polar order within individual nodes, across all occupied nodes in the lattice.It is evident that, if particles are homogeneously distributed among all channels in every non-empty node, then P ℓ = 0. Conversely, if for every node, all particles are contained within a single channel, then P ℓ = 1.Fig. 2g shows a plot of the global order P g and the local polar order P ℓ , as a function of the model parameters after several time steps.After the transition, the local polar order P ℓ saturates up to its maximum value, while the global polar order P g exhibits non-monotonic behavior, reaching a maximum exactly at the transition point.This means that, after the transition, most occupied nodes accumulate particles within a single node.At the transition point, few condensates form and therefore the condensate with the most particles increases the overall global polar order.After the transition, several condensates with their own orientation form at once.The difference among the condensates' orientation drives the goal polar order down.Thus, disconnected condensates are present for all parameter values after the transition.

Emergence of nematic filaments
The strict local nature imposed by zero-range interactions seems to prevent nodes to exchange orientational information in order to establish extended ordered structures.This is true for polar alignment: the polar order of two neighboring nodes is not coupled.Counterintuitively, for nematic alignment, and in contrast to polar alignment, nodes can exchange orientational information and particles self-organize into stable, elongated, nematic filaments that lead to global nematic order; Fig. 3a.We stress that these nematic filaments should not be confused with, and are different from the nematic bands observed in off-lattice models [30]: i) the width of these filaments is always of one node and does not vary with the distance to the critical point as occurs in nematic bands [30], and ii) as consequence of this, the transversal instability observed for nematic bands, precluding global nematic order, is not present [60,39].
How can an extended nematic filament emerge from zero-range interactions?Picture an initial configuration where a large number n of particles are placed on a single channel of a lattice site.As in the polar case, q ϑ = 1/6 is not stable for n > n * and particles accumulate in two opposite velocity channels.In the next time step, particles move to neighboring sites associated to these two velocity channels.If the density in each new lattice site is high, the same process occurs in the new sites, and thus information on the initial orientational order, propagates forward and backward, as shown in the sketch of Fig. 3c.If we consider that this process starts taking place in an area surrounded by a disordered population of particles, then as the nematic filament expands enlarging its scattering cross-section, these particles feed the growing filament until it eventually percolates.The described mechanism leads to global nematic order as evidenced by the global nematic order parameter |Q| = r α∈Θ q α (r) e 2iα .Note that the emergent filament contains a macroscopic fraction of the particles in the system, as evidenced by the behavior of the filament density N f with β, which is similar to the one of N c for polar alignment; Fig. 3d.Furthermore, we observe that for fixed ρ, the N grows as L 2 , and thus N f as L, which implies that for L → ∞, this quantity diverges: the nematic filament acts as a condensate with a macroscopic particle fraction located at it moving in two opposite directions.The local nematic order, Q ℓ = r α∈Θ q α (r) e 2iα /N o was also observed in simulations.Here as well, we observed a monotonic increase of the local nematic order and a non-monotonic behavior of the global nematic order; Fig. 3b.As before, this indicates the formation of few filaments with most particles at the transition point, driving polar order up, and several filaments with different orientations far above the transition point, mantaining a high local nematic order, but inhibiting nematic order at a global level.

Hydrodynamic equations
The onset of the (macroscopic) symmetry breaking, while the system remains near the disordered phase, can be described in an approximate manner by calculating the expected occupation numbers of lattice velocity channels.However, note that well-developed condensates and filaments cannot be described by smooth PDE solutions.As a first step, let us define f ϑ (r, k) := s ϑ (r, k) .By performing an innode mean-field approximation and employing the binomial statistics of velocity channel occupation numbers, one obtains the lattice-Boltzmann equation (LBE): is the density, equal to the expected value of the total number of particles in each lattice site, and m is the mean-field alignment probability.We perform a discrete Fourier transform on ϑ to arrive to hydrodynamic fields that depend on position and time.For polar alignment (m = 1), the relevant hydrodynamic fields reduce to the local density ρ := ϑ∈Θ f ϑ and local polar order (or momentum) M := ϑ∈Θ f ϑ (cos ϑ, sin ϑ).The temporal evolution of these fields, in dimensionless form, leading to unit active speed, is given by: where For detailed derivation, see [61].The term u can be approximated, assuming β ≪ 1, by a Taylor series.If we approximate it up to third order in β, we find u = β 2 M 1 − β 2 8 M 2 , and by replacing this expression into Eq.(4b), we obtain an equation that is invariant under arbitrary rotations of M as occurs in Toner-Tu theory [35].However, if the expansions is performed up to 5th order, we uncover that steady solution of the form M = M (cos θ, sin θ), with M a scalar constant and θ an angle, cannot exist for arbitrary values of θ.The only possible solutions correspond to θ = aπ 3 , with a ∈ [0, 1, • • • , 5], i.e. to one of the possible 6 directions of the lattice; for details see [61].
Differences with the standard off-lattice Vicsek model can be shown even using the third order expansion in β mentioned above.There are, as expected, two spatially homogeneous solutions, corresponding to M = 0 (disorder) and M > 0 (order).Performing a linear stability analysis around the disordered homogeneous steady state, and assuming perturbations along a single direction, we find that the system is stable against uniform perturbations for an initial average density β ρ < 2, showing that the product β ρ is the effective bifurcation parameter.The second steady state corresponds to spatially homogeneous polar order state with . The ordered homogeneous state only exists when β ρ > 2, where the disordered state becomes unstable.For non-uniform perturbations, the spatially homogeneous ordered state, that is never observed in the LGCA, is stable only for the parameter range 2.564 β ρ < 3. Otherwise, the real parts of the eigenvalues grow with increasing wavenumber, indicating that small wavelength perturbations dominate, as confirmed by direct integration of the hydrodynamic equations.This instability indicates a drastic density increase in localized areas of the space, which can be associated to the emergence of condensates at this level of description.Numerical evidence strongly suggests that spatially homogeneous ordered states, and ordered polar bands are not possible in the LGCA model, even when when forcing the system by using these configurations as initial conditions (see SI).
For nematic alignment (m = 2), the relevant hydrodynamic fields include nematic order, defined for simplicity as Q := ϑ∈Θ f ϑ (cos 2ϑ + i sin 2ϑ), and the hydrodynamic equations read: where d := ∂ x + i∂ y , * indicates complex conjugate, R (•) denotes real part, and fast relaxation of M was assumed.It is worth stressing, as we explain below, that according to Eq. ( 5) nematic order is constrained to exist, as occurs for M, on the orientations allowed by the lattice structure.Eq. ( 5) possesses spatially homogeneous disordered and nematically ordered solutions.The disordered state is stable when β ρ < 4, where the bifurcation depends only on the product β ρ, as before.The order in the spatially homogeneous steady state is given by Q = 2 ± 2 9 − 32 β ρ e i 2kπ 3 , k ∈ {0, 1, 2}.Two features are worth discussing.First, the ordered state exists when β ρ > 32  9 , i.e. it exists even when the disordered state is stable, indicating a subcritical bifurcation.Second, while for 32  9 < β ρ < 4 three pairs of solutions with different moduli, but same argument arg (Q) ∈ 2kπ 3 , k ∈ N exist, for β ρ > 4, six different solutions exist, each one with a different argument corresponding to one of the lattice directions ϑ ∈ Θ.For β ρ > 4, solutions with arg (Q) = 2kπ 3 are stable, while those with arg (Q) = (2k+1)π 3 are unstable.This is a consequence of the constraints imposed on Q by the lattice, such that nematic structures can only emerge along the three hexagonal lattice axes, as observed in the LGCA simulations.Again, simulations suggest that spatially homogeneous ordered states, as well as spatially extended bands are not stable, and consistently break into nematic filaments, or revert back to the homogeneous disordered state.In this case, however, ordered states are much less unstable than in the polar case, and need long times in order to break apart (see SI).

Discussion
The original LGCA was introduced with the intention of reproducing basic fluid mechanics properties.However, the initial study performed on a square lattice [62] displayed a series of anomalies due to the use a finite number of velocity channels.Later on, it was found that such anomalies vanish in hexagonal lattices [53], where the phenomenology of Navier-Stokes equations is recovered, despite the use of only 6 velocity channels.Similarly, it is tempting to speculate that an active LGCA with polar velocity alignment in an hexagonal lattice (and short-range interactions) has to reproduce the large-scale properties predicted by the Toner-Tu equations [35] for a (compressible) active polar fluid.After all, it is expected that most microscopic details are, at macroscopic scales, irrelevant except for the range and symmetry of the interactions.However, our study shows that spatially homogeneous ordered phases (cf.Toner-Tu phase [35]) or polar bands (cf.[31,46]) cannot exist in the LGCA model defined in this work.Moreover, in this LGCA, order emerges via a fundamentally different phase transition from the ones previously reported [28,35,31,47,46,48,30] in active matter.In LGCA, the emergence of global order is intrinsically related to a condensation transition -that occurs in the absence of an external potentialwhich shares similarities with the ZRP.For the polar LGCA, it is associated to the formation of a polar macroscopic condensate, while for the nematic LGCA to the collapse of a macrocopic number of particles in the systems into a nematic filament (that leads to global nematic order; cf.[30,39]).
After their introduction as fluid-dynamical models, LGCA have been used to model biological systems, specifically for modeling collective cell migration (see [63] for a review).The original LGCA formulation considers an exclusion principle, which limits the maximum number of cells which may be present at a single spatial point at a given time.On the one hand, our model could be considered as a 2D projection of a 3D system.For example, the accumulation of particles at certain lattice sites could correspond to the formation of an elongated structure by the piling up of agents, such as during the fruiting body and slug formation in D. discoideum [64].On the other hand, cells are quite compressible, and are able to accumulate in great numbers.Experimentally, it has been found that cells can greatly reduce their volume during tissue growth [65,66].Cells also regularly change their size to maintain epithelial homeostasis [67].For example, the photoreceptor epithelium shows a wide variety of receptor packing densities, corresponding to photoreceptor cells of varying sizes [68].Furthermore, invasive cancer cells can be ten times as compressible as their healthy counterparts [69], and cell compressibility has been shown to correlate with invasive potential [70].Moreover, it has been shown that densely packed cells confined to narrow environments can still move collectively, as long as they mantain sufficiently high adhesive contacts and velocity correlation levels [71].The present LGCA model extension, which does not consider an exclusion principle, represents a first step towards a better mathematical representation of these cellular aggregation and compression processes.
Finally, we strongly suspect that the condensation reported here is related to the recent findings that indicate that here exist fundamental differences between additive and non-additive interactions in vectorial active matter [72,38].In particular, we speculate that the emergent macroscopic behavior can be affected, not only by modifying the symmetry or range of the interactions (see [48] for an analogous model with extended range and an exclusion principle), as it is usually expected, but by eliminating the dependency on n(r, k) in Eq. ( 2), for instance by defining for m = 1, H(r, k) ϑ = P(r, k) • c ϑ .This will be the focus of future investigations.
AVB and LB: supervision, methodology, and writing (review and editing).AD: Funding acquisition, supervision, project administration, and writing (review and editing).FP: Conceptualization, formal analysis, investigation, methodology, project administration, supervision, and writing (review and editing).

Figure 2 :
Figure2: Dynamics of the polar, BA-LGCA.a. Simulation snapshot; node occupancy is color coded, with dark color indicating large particle number.The inset shows the state of a node.See[61] for movies.b.Global polar order P vs. average occupation per node n. c.Probability of a particle escaping a polarized condensate given its particle number.Inset: Probability of a particle staying in a polarized cluster.d Histogram of the number of particles in all lattice sites at two different time steps.e. Mean cluster size (blue) and polar order (red) as function of time; the dotted line corresponds to t 1/2 , the characteristic growth in classic ZRPs[56], while the dashed line corresponds to t 3/4 .Average over 200 simulations with β = 4.3 and ρ = 12.f.Polar order P g and condensate size N c as a function of β ρ for average ρ = 12.g.Global polar order P g and local polar order P ℓ as a function of β ρ for average ρ = 12. 150 simulations were averaged.A lattice size of 60 × 60 was used for all simulations.Parameter values were recorded after 1000 time steps.

Figure 3 :
Figure 3: Dynamics of the nematic, BA-LGCA.a. Snapshot of a simulation, with high node occupation numbers displayed in dark color; inset shows the state of a node.See [61] for movies.b.Global nematic order |Q| and local nematic order Q ℓ as a function of β ρ for average ρ = 12.c.Schematic depicts the filament formation process.Arrows indicate the direction of particle groups, with colors, from red to magenta, decreasing particle number.d.Filament density vs. β ρ averaged over 150 simulations for fixed ρ = 0.6.a Lattice of size 60 × 60 was used for all simulations.Parameter values were recorded after 1000 time steps.