Morteza Nattagh Najafi,Rafe Md. Abu Zayed,Seyed Amin Nabavizadeh
Swarming Transition in Super-Diffusive Self-Propelled Particles
super-diffusive Vicsek model,Levy flights,second-order phase transition,critical exponents,05.,05.20.-y,05.10.Ln,05.45.Df
A super-diffusive Vicsek model is introduced in this paper that incorporates Levy flights with exponent . The inclusion of this feature leads to an increase in the fluctuations of the order parameter, ultimately resulting in the disorder phase becoming more dominant as increases. The study finds that for values close to two, the order–disorder transition is of the first order, while for small enough values of , it shows degrees of similarities with the second-order phase transitions. The article formulates a mean field theory based on the growth of the swarmed clusters that accounts for the decrease in the transition point as increases. The simulation results show that the order parameter exponent , correlation length exponent , and susceptibility exponent remain constant when is altered, satisfying a hyperscaling relation. The same happens for the mass fractal dimension, information dimension, and correlation dimension when is far from two. The study reveals that the fractal dimension of the external perimeter of connected self-similar clusters conforms to the fractal dimension of Fortuin–Kasteleyn clusters of the two-dimensional Potts (Ising) model. The critical exponents linked to the distribution function of global observables vary when changes.
Swarming Transition in Super-Diffusive Self-Propelled Particles

A super-diffusive Vicsek model is introduced in this paper that incorporates Levy flights with exponent . The inclusion of this feature leads to an increase in the fluctuations of the order parameter, ultimately resulting in the disorder phase becoming more dominant as increases. The study finds that for values close to two, the order–disorder transition is of the first order, while for small enough values of , it shows degrees of similarities with the second-order phase transitions. The article formulates a mean field theory based on the growth of the swarmed clusters that accounts for the decrease in the transition point as increases. The simulation results show that the order parameter exponent , correlation length exponent , and susceptibility exponent remain constant when is altered, satisfying a hyperscaling relation. The same happens for the mass fractal dimension, information dimension, and correlation dimension when is far from two. The study reveals that the fractal dimension of the external perimeter of connected self-similar clusters conforms to the fractal dimension of Fortuin–Kasteleyn clusters of the two-dimensional Potts (Ising) model. The critical exponents linked to the distribution function of global observables vary when changes.

There has been a significant increase in attention towards emergent behaviors, particularly collective motions in active systems [1]. These systems exhibit large-scale cooperative phenomena and sometimes scale-invariant patterns, resembling the standard theory of statistical mechanics and critical phenomena. To draw a connection, the concept of “phase” is assigned to each collective mode of behavior of agents, and the term “phase transition” is used. Swarming is an example that is observed in various biological systems, e.g., fish schools with milling dynamics [2], bacterial systems [3,4,5,6], biofilm formation [7], and marching locusts [8]. Observations of self-organized critical behavior in systems such as midge swarms [9] prompt the question of what the significance of criticality is in these systems. Among the mathematical models [10], the Vicsek model (VM) [11] serves as a minimal prototypical example of active systems with self-propelled constituents that show swarming transition. The VM is a simple yet extensively studied model due to its ability to exhibit new features and diverse modes of collective behavior. For example, it was mapped onto complex networks [12] and XY spin chains [13,14,15]. Other studies have focused on the impact of angular noise [16] and hierarchical societies [17], with the aim of uncovering the key physical factors that contribute to emergent behaviors and its response to additional interactions/physical parameters; see [18]. Initial studies asserted that the VM exhibits critical behavior at the transition point [11,13,14,15], but later research showed that the transition is actually of the first order [19,20]. This discrepancy was attributed to the small finite size effects of the system. However, other studies suggest that the order of the transition may depend on how noise is added to the system [21] or could be an artifact of finite size and boundary conditions [22]. This leaves an important unanswered question: How do the long-range interactions [9], which lead to criticality in many systems [23,24,25], impact the collective motion when these individuals are swarming together? A recent study [26] found that the lack of fluctuations, or complete synchronization with the flock, in the VM and other classical models, such as the Cucker–Smale model, can reduce the adaptive response of a flock, which may be undesirable in certain scenarios. To increase the adaptive sensitivity to external threats, criticality and correspondingly large fluctuations are essential, as information corresponding to a local perturbation moves faster due to the scale-free correlations [27,28]. It is now well known that criticality is a vital ingredient for a flock state to survive in the presence of external threats [29,30]. We argue that the presence of scale-free stochasticity in, for instance, the flight distance of the constituents is a crucial factor for achieving criticality, which is absent in classical counterparts. To establish a closer connection with real-world scenarios, it is essential to incorporate scale-free stochasticity that generates anomalous diffusion, which is frequently observed in active systems. Examples include Levy flights of wandering albatrosses [31] and other animals [32], super-diffusive intracellular transport [33], entangled F-actin networks [34], microtubule-associated motors within a living eukaryotic cell [35], living yeast cells [36], mRNA molecules inside live E. coli cells [37,38], telomeres in the nucleus of mammalian cells [39], biomolecules in solution and living cells [40], the pathway of an Adeno-associated virus [41], and epithelial cell migration (with Levy flights as the asymptotics of the q-Weibull distributions) [42]. In most cases, super-diffusion is observed, which enhances the exploration process (e.g., foraging for animals) [31,32]. In this paper, we investigate the impact of introducing Levy flights into the system dynamics, resulting in super-diffusive stochasticity. This leads to the criticality of the system for sufficiently small step exponents () of Levy flights. The paper is organized as follows: In the next section we describe our model. Section 3 presents the simulation results, while Section 4 presents the mean field results. Section 5 investigates the geometrical and global features of the model and is divided into two subsections. The first subsection, Section 5.1, describes the mass fractal dimension and higher-order dimensions of the density field. In the second subsection, Section 5.2, we present our contour line analysis for the density field. We conclude the paper with a summary of our findings in the final section.

In the ordinary VM, the agents (labeled by , where represents the total number of active particles, is the density of the particles and L is linear size of the system) undergo correlated random walks in the system with an interaction range R such that any active particle inside a disk of radius R is fully seen by the central particle. The time evolution of the position of the ith particle at time t () is given by where is the direction of motion of the ith particle at time t, , is the distance that the ith particle traverses in the time interval , shows the set of particles at a distance less than or equal to R from the ith particle, and . is a uniform random number in the interval , and is the strength of the disorder. In the ordinary VM, where the particles have a constant velocity (), the coherence order parameter defined as is zero (non-zero) in the disordered (ordered) phase. This indicates that the particles move in a spatiotemporally coherent (incoherent) fashion [11], and provide information about the degree of orientation of the particles’ motion. There is a transition point , above (below) which and the systems are in the disordered (ordered) phase. is defined as the time average over a considerable time interval. In our model, we let the particles’ velocity obey the Levy flight distribution: where is a “step index” which generates correlations, and is a step function ( for , and zero otherwise). The occurrence of rare events (long-range) flights in the process introduces correlations by allowing the particles to access and persist in specific regions of space. Consequently, correlations in the particles’ positions over time emerge [43]. is considered in our model to prevent unphysical rare events, such as flights above the threshold speed imposed by physical conditions (), which serves as an IR cut off in the problem. In the analytical calculations, we require also a UV cutoff to ensure that the average flight is well-defined, where is a normalization constant. It is customary to consider this cutoff as the smallest scale in the problem, e.g., a lattice constant in the models on the lattice. We set it to R, which is the lattice constant in our setup. For , the average length and its variance diverge as the limits of or are approached, which results in the central limit theorem being invalid. In this situation, the distribution of the Levy random walkers is described by -stable distributions as reported in [44]. We refer to the particles in this regime as super-diffusive active particles because the diffusion exponent is [45]. For the ordinary VM, there is a coexistence of ordered and disordered phases at the transition point (which is of the first order). This leads to the existence of well-defined phase boundaries, known as phase coexistence [17], and results in a bimodal distribution function for the order parameter. The bimodal distribution can be either spatial or temporal, indicating that the two phases may either spatially coexist or dominate during different time intervals [17,46]. The existence of a bimodal distribution function may directly lead to a gap for at the transition point [47], around which the hysteresis effect is observed. In this state, the system stays during most of the observation time in the vicinity of one peak of or, in other words, a metastable branch is observed. A hysteresis loop results from the system’s resistance to entering the new peak [17]. The gap is defined as the difference between these two peak points exactly at the transition point. This is in accordance with the discontinuous (first order) phase transitions, where the phases coexist with well-defined boundaries. However, with the second-order phase transitions, the system shows self-similar patterns at the transition point, making it challenging to attribute a specific phase to a part of the system [47]. In our model, the scale-free stochasticity induces the evaporation of the two phases by promoting “tunneling” between the peaks. The amplitude of this effect is determined by the value of .

We simulated the system for and 256, , and . Additionally, R and were set to 1. The active particles were initially distributed randomly in the system with random uncorrelated movement orientations. For the systems that exhibit a hysteresis effect, we moved forward (increase ) and backward (decrease ), and for both cases, we changed by in each time step. To control statistical fluctuations, we generated three samples at each time step. We used the maximum likelihood estimation (MLE) method [48] to estimate the best values for the exponent and corresponding error bars. For the data collapse analysis, we selected the graph of the largest system size as a base. We calculated the value for all other graphs by measuring the cumulative distance of the points from this base graph. The distribution of values was then used to determine the best exponents and corresponding error bars. This technique was used to better understand the relationship between different variables in the system. The density pattern of the model in the transition point is strongly dependent on the value of as shown in Figure 1. Figure 1 investigates two different scenarios of a physical system with different values of . In the first case, , corresponding to Figure 1a, the system displays a self-similar pattern without any phase separation. The second case, where , corresponding to Figure 1b, exhibits a strip-like dense (ordered) phase in a background of a dilute (disordered) phase. The self-similar pattern for is reminiscent of the second-order transition, while the phase separation for is a fingerprint of the first-order transition. To quantify the order of the transition, we show the time series of and the corresponding probability distribution in the Figure 2a for , and for (left, with a single peak) and (right, with a bi-modal structure). The single peak structure was observed for all small enough values, suggesting that the transition in our model for these values does not follow the first-order phase transition paradigm. In contrast, the transition is discontinuous for values around two. In Figure 2b, the - plot for different values demonstrates that increasing shifts the graphs to the left, suggesting that the disordered phase is stabilized by higher values. The Binder cumulant method was used to extract the transition points. This method is especially an effective tool for identifying the type of transition, as it is a continuous quantity for second-order phase transitions. It is defined as This function shows a sudden deep minimum at the transition point, which is characteristic of first-order phase transitions (upper inset of Figure 2b for ). This function is inspected further in Figure A1a in Appendix A. There, it is shown that the depth of the valley decreases as decreases, which reveals the system crossover to regimes with properties similar to second-order phase transitions in sufficiently small values. The observation is consistent with the modal structure of the distribution function, Figure A1b (see Appendix A for more details). While is theoretically expected to be L-independent at the second-order phase transitions [49], statistical uncertainties prevent the points from precisely intersecting at a single point. The method proposed in [50] is utilized to estimate in the thermodynamic limit. The method relies on identifying the intersection point () of two graphs corresponding to two successive Ls and extrapolating the obtained point to . The resulting transition point is shown in the lower inset of Figure 2b, exhibiting a decreasing behavior in terms of , i.e., the disordered phase dominates as increases. Although the observations mentioned suggest similarities with second-order phase transition for small enough values, additional statistical evidence is required to validate this hypothesis. According to the standard theory of second-order (continuous) phase transitions, the order parameter and the order parameter fluctuation satisfy the finite size scaling hypothesis [51]: where , and and are some exponents, and and are universal functions with the asymptotic behaviors , , and . These exponents are related via a hyperscaling relation [51] where d is the Euclidean dimension of space, which is two here. In the main part of Figure 2c and its upper inset, we show the finite size dependence of in terms of and at for , respectively. While the peak goes to the left as L increases, is L-independent exactly at as is shown in the lower inset of Figure 2c. The values of obtained through this data collapse analysis agree remarkably well with those determined via the Binder cumulant method. The results of the data collapse analysis presented in Figure 2d confirm that our model obeys finite size scaling relations for sufficiently small values in accordance with second-order phase transitions. Note that the exponents and are robust against changes in , and also of the fitting increases as increases (the worse fitting is for in our set), showing that the scaling hypothesis works for sufficiently small values, and the finite-size scaling Equation (4) is not appropriate for large values. One may use the theory of finite-size scaling for first-order transitions [47] for those cases, which is outside the scope of the present study. is reported in Figure 3a, exhibiting a pronounced peak at the transition point. (the amount of at its peak) also scales with L by the exponent . Figure 3b displays the dependence of and , indicating that an increase in corresponds to an increase in fluctuations of , which is the reason why increasing favors the disordered phase. A data collapse according to Equation (4) for is presented in Figure 3c, giving for all values. The exponent is found to be in the interval [1.4–1.6] for all values, which agrees with the hyperscaling relation Equation (5). The behavior of the system in terms of density is crucial in understanding the emergent collective behaviors of self-propelled particles. We found evidence regarding the density-driven order–disorder transition for super-diffusive active particles by inspecting the properties of in terms of . We first considered the transition in terms of for low density regime (). In this case, similar properties, such as high density limit , were observed, shown in Appendix A (Figure A2). We also found also a density-driven order–disorder transition, which is depicted in Figure 4, showing this transition in terms of for and . The critical density in this case was found to be , and the exponents of the transition were and . Again, the transition showed similarities with the second-order phase transition paradigm (such as the finite size scaling, based on which we extracted the exponents using the data collapse analysis) but was not completely fitted to it. To conclude this section, we observed that our model undergoes an order–disorder phase transition for all values, and decreases as increases. Our findings reveal that the phase transition is of the first order for large values, similar to the ordinary Vicsek model. Conversely, when decreases sufficiently, the transitions align better with second-order phase transitions, or at least show similarities with second-order phase transitions. We observed a crossover from first-order to second-order transition by decreasing .

To comprehend the impact of on our model’s characteristics, especially to address why the disordered phase becomes increasingly stable with increasing , we develop a mean field theory. We consider a swarmed cluster with an average radius r, and study the dynamics of the particles that leave or enter the cluster as an impact of Levy flights. We quantify the rate of particle movement out of (or into) the system with (), and define the density of the population within (or outside of) the swarmed cluster as (). These mean field quantities are obtained through a straightforward argument illustrated in Figure 5. is not sensitive to the movement of the swarmed cluster since the Levy distribution is nearly invariant under a small boost. Therefore, we use the setup in Figure 5a for determining as follows (noting that on average the fraction of particles go to a required direction): where is the Levy accumulated probability density. To calculate , we allow the swarmed cluster to move by a distance of during a single time step. The preferred direction can be calculated by the previous step given in Figure 5b. To calculate , we note that the particles in the horizontal distance (on the green bar) with a flight length remain in the coherent swarmed cluster. Based on this, we calculate the number of particles that leave this area. Noting that the required accumulated probability is , and also (noting that on average the fraction of particles go to a required direction), we finally find the relations where and . To obtain the mean field relations for and , we use the following argument: In the ordered phase, there are typically multiple particles inside a disk with area , whereas in the disordered phase, particles are unlikely to encounter each other within a single Levy flight, resulting in no particles within an area of . This leads to the relation . It is worth noting that the relationship is not qualitatively dependent on the power used in the equation. This, along with conservation of the total number of particles , gives the following relation: Noting that is the rate of change of the average number of active particles inside the swarmed cluster, one can determine the dynamical behaviors of the model in terms of r. The average r (which we call ) is the fixed point of the dynamical behavior of n, i.e., . The average size of the swarmed clusters, indicated by , is used to determine the stability of each phase for a given value of . The relationship between the average size of the swarmed clusters and the parameter is presented in Figure 6. The figure shows that as increases, the disordered phase becomes more stable. To make connection between and , we first note that and . Therefore, we expect that is a decreasing function of , which agrees with the simulation results. As pointed out in [52] for a standard Vicsek model with large densities, the system may percolate and form a huge cluster. This happens when in our MF arguments. For the percolation probability, we should seek the conditions that lead to , or equivalently . In our model, the probability of percolation is greater for the smaller values. Given that in this limit, the fluctuations are much greater than for large values (because of the criticality of the system at this limit), it is consistent with the percolation theory which tells us that at the percolation threshold, that fluctuations are maximal. In the inset of Figure 6, we show the relation between and , which is monotonic increasing function. If one fits this relation using a power-law function, the exponent would be , although the range of the quantities are too small to deduce a power-law form (it commonly should be more than one decade).

In addition to local characteristics, systems undergoing continuous phase transitions display global geometric features, which have been the subject of numerous analytical and simulation studies. This characterization can reveal previously hidden aspects of the models that are not apparent with only local observables. Section 5.1 focuses on the fractal analysis of the density of active particles at the transition point. In Section 5.2, the critical loop ensemble (CLE) of the iso-density loops on the system is analyzed.

We employ the fractal analysis method for the density configurations [53] that we obtain at the transition points. The density configurations are first converted into black-and-white images, after which the white pixels are statistically analyzed; see Figure 1. The system is meshed using boxes of a specific size linear size , and the statistics of the filling fraction of each box is calculated. A pixel is considered white or occupied if the density of active particles at that site is greater than the spatial average of density over that sample, where is the total number of pixels in the system. The filling fraction of the ith box is , where is the number of white pixels in the ith box [54,55]. Note that , where is the total number of boxes. If we attribute a local mass for ith box as and a total mass to the cluster as , where is a Kronecker delta, then the box counting fractal dimension is obtained as In a multifractal system, this exponent depends on the scale that we are considering or changes from region to region. A unified standard theory called multifractal analysis was previously developed, which employs a generalized partition function that yields a spectrum of exponents, including the fractal, information, and correlation dimensions [53]. This q-generalized partition function is related to the qth moment of the fluctuations of , and is defined as where q is a moment. For scale-invariant systems, scales with in a power-law form with the exponent , but the exponent may not be a unique number in all scales:The generalized q-dimension is then defined as so that . Note that if one interprets as a probability associated with a small segment () of the system, then plays the role of a normalized q-Renyi entropy () in the thermodynamic limit so that Therefore, the mass fractal dimension of samples is related to Renyi entropy It is worth noting that the hypothesis of scale invariance, as described in Equation (12), has led to the fact that Renyi entropy is proportional to and not ( in our case) as expected for the ordinary systems. This serves as an important characteristic of the scale-invariant systems, for which the system is not extensive [56]. The information dimension associated with the Shannon entropy is obtained in the limit the fact that relates it to so that Finally the correlation dimension is defined as where where is the position of the kth white pixel (not box), and is a step function. It is shown that [54] We generated configurations (density snapshots) at the transition points to investigate the anomalous dimensions. Figure 7a shows in terms of for , and the inset shows in terms of q, which is well described by a linear function. The numerical values for the dimensions are reported in the inset of Figure 7b. For the small values, the exponent remains constant and stable across various values. Interestingly, the mass fractal dimensions are lower than one, which is generally possible for the fractals with fractional filling boxes. However, as approaches two, power-law fittings fail to fit, and the resulting exponents deviate from the others. This is not surprising since, at these points, the system does not display a fractal structure.

The critical (conformal) loop ensemble (CLE) theory enabled another type of classification of two-dimensional (2D) critical models based on their global geometrical properties [57,58]. When considering this approach, the focus is on random curves that can be transformed into dynamic stochastic paths, or exploration processes, within a connected domain in the plane. This idea was originally suggested by Loewner and is known as stochastic Loewner evolution (SLE) [59,60], which is now widely recognized as a means of characterizing the interfaces of two-dimensional statistical models using growth processes. The method can classify these interfaces into one-parameter classes, with the diffusivity parameter serving as the representative parameter, called SLEκ [59,60]. Many other exponents are related to , such as the fractal dimension of level lines or interfaces, which is [61] (these extended objects are fractal paths, or sometimes loops). As another example for the Potts models, it is shown that , and , where is the dual of , which is given by the equation such that, for example, for [62]. To investigate this, we start with the fractal dimension of loops () for a single connected cluster, which is defined as the scaling exponent between the loop length l and the loop gyration radius r. The latter is defined for a closed path as , where is the loop center-of-mass. is then defined by the relation , where 〈〉 denotes the ensemble average. We apply the Hoshen–Kopelman algorithm [63] to identify the connected components of the clusters. This algorithm involves coloring the entire cluster and assigning different colors to separate clusters while traversing the sample. We analyze these connected clusters by measuring their external boundaries with length l, as well as their gyration radius r and mass, which we denote as . The scale-invariant properties of the distribution function of , l, and r are also evident (excluding finite size effects), with . Here, represents the corresponding scaling exponent. We consider clusters that are associated with the white-and-black density pattern. The scaling hypothesis in the transition points is supported by the findings presented in Figure 8a,b. It should be noted that when the value of approaches two, the system loses its scale-invariance property and also becomes anisotropic, indicating that the critical or scaling exponents cannot be considered dependable in this range. Given the system and the associated exponents, the fractal dimension is observed to remain constant at approximately for all values of , indicating its robustness. This exponent is associated with (in the case of conformal invariance). However, the scaling exponents () exhibit changes as varies. This is the first time in this study that we see the set of exponents change as a function of , suggesting that a range of systems can be visited within this interval of . It would be intriguing to compare the results with the ones for the geometric and Fortuin–Kasteleyn (FK) clusters of the critical Potts (Ising) model, the diffusivity parameters, which are and , respectively [62]. It also was previously shown that , and , and also [64]. The exponents that we found for our model for small values are in agreement with the exponents explored above. By setting in our model, we obtain , , and . These values are consistent with our simulation results, except for the exponent, which shows a discrepancy. This suggests that the interfaces of our model do not exhibit conformal invariance. Hence, we conclude that our model is self-similar for sufficiently small values, and exhibits some similarity to the Potts model but is not a perfect fit for this model.

In this paper, we studied a super-diffusive variant of the Vicsek model by introducing scale-free (Levy) stochasticity to the flights taken by the active particles during each time step. As a result of this modification, the transition, which is of the first order for the conventional Vicsek model, shows similarities with second-order phase transitions for small values. Since we observed some features of first-order phase transitions, we denote this regime as weakly second-order phase transitions. In contrast, for values around two, our model displays first-order phase transitions, but like the Vicsek model, it also exhibits some characteristics of scale invariance. The latter led Vicsek et al. to incorrectly conclude that their model was of the second order and calculate “fictitious exponents” [11]. This occurs when two peaks of are in close proximity and difficult to distinguish as shown in the right inset of Figure 2a. Since the values of and the “fictitious exponents” for our model (for values close to two) are similar to those of the ordinary Vicsek model, we can infer that our model exhibits the same characteristics as the ordinary Vicsek model at the transition point for values close to two. Additionally, for , the -stable Levy systems are unstable towards a fixed point that is characterized by a Gaussian distribution function [44], indicating that perturbing the ordinary Vicsek model with Gaussian-distributed flights does not alter the fundamental properties of the Vicsek model, and is therefore an “irrelevant perturbation”, while for the small values, it is similar to second-order phase transitions and the corresponding perturbation is relevant. We observe a crossover from first-order phase transition (large values) to weakly second-order transitions (small values). We developed a mean field theory for our model, which successfully describes why the disordered phase becomes more stable as increases. The geometrical properties of the model at the transition points were also investigated. We found a series of anomalous dimensions, including the mass dimension, the information dimension and the correlation dimension. Our critical loop ensemble study shows that this system has similarities to the Pottes (Ising) model, while the exponent does not match.