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 Q=2 Potts (Ising) model. The critical exponents linked to the distribution function of global observables vary when α changes.


Introduction
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.

The Model
In the ordinary VM, the agents (labeled by i ∈ [1, N], where N = ρL 2 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 (x i (t)) is given by where θ i (t) is the direction of motion of the ith particle at time t, v i (t) = v i (cos θ i (t), sin θ i (t)), v i ∆t is the distance that the ith particle traverses in the time interval [t, t + ∆t], {R i } shows the set of particles at a distance less than or equal to R from the ith particle, and . ζ i is a uniform random number in the interval [− 1 2 , 1 2 ], and η is the strength of the disorder.
In the ordinary VM, where the particles have a constant velocity v VM the coherence order parameter defined as φ η (t) 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 η c , above (below) which φ(η) ≡ φ η (t) t and the systems are in the disordered (ordered) phase. t 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 Θ(x) is a step function (Θ(x) = 1 for x ≥ 0, 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]. Θ(x) is considered in our model to prevent unphysical rare events, such as flights above the threshold speed v max imposed by physical conditions (l max ≡ v max ∆t), 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 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 1 < α < 2, the average length and its variance diverge as the limits of → 0 or l max → ∞ 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 α 2 > 0.5 [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 p(φ) 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 selfsimilar 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 α.

Simulation Results
We simulated the system for L R = 32, 64, 128 and 256, α ∈ [0.8 − 1.95], and ρ = 2. Additionally, R and ∆t 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 δη = 3 × 10 −6 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 χ 2 value for all other graphs by measuring the cumulative distance of the points from this base graph. The distribution of χ 2 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.  To quantify the order of the transition, we show the time series of φ and the corresponding probability distribution p(φ) in the Figure 2a for η > η c , η = η c and η < η c for α = 1.5 (left, with a single peak) and α = 1.95 (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 α = 1.5). 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 G η 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 η c in the thermodynamic limit. The method relies on identifying the intersection point (η c (L, L )) of two graphs corresponding to two successive Ls and extrapolating the obtained point to L, L → ∞. The resulting transition point η c 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 χ φ ≡ L 2 φ 2 − φ 2 satisfy the finite size scaling hypothesis [51]: where ≡ η c −η η c , β and ν and γ are some exponents, and F φ and F χ are universal functions with the asymptotic behaviors lim These exponents are related via a hyperscaling relation [51] νd = γ + 2β, 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 p(φ) at η c for α = 1.5, respectively. While the peak goes to the left as L increases, L β/ν φ is L-independent exactly at η = η c as is shown in the lower inset of Figure 2c. The values of η c 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 χ 2 of the fitting increases as α increases (the worse fitting is for α = 1.95 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. χ max (the amount of χ φ at its peak) also scales with L by the exponent γ/ν. Figure 3b displays the α dependence of χ φ and χ max , 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 densitydriven 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 (ρ = 0.3). In this case, similar properties, such as high density limit ρ = 2, 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 α = 1.5 and η = 3. The critical density in this case was found to be ρ c = 1.16 ± 0.01, and the exponents of the transition were β ρ = 0.16 ± 0.02 and ν ρ = 1.32 ± 0.14. 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 η c 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 α.

Mean Field Arguments
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 n out (n in ), and define the density of the population within (or outside of) the swarmed cluster as ρ in (ρ out ). These mean field quantities are obtained through a straightforward argument illustrated in Figure 5. n in 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 n in as follows (noting that on average the fraction 1 4 of particles go to a required direction): where p acu (l > r) ≡ ∞ r p Levy (l)dl is the Levy accumulated probability density. To calculate n out , we allow the swarmed cluster to move by a distance ofl α during a single time step. The preferred direction θ can be calculated by the previous step given in Figure 5b. To calculate n out , we note that the particles in the horizontal distance y (on the green bar) with a flight lengthl α − (r + y ) < l <l α + (r − y ) 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 p acu (r,l α , y ) ≡ l α +(r−y ) max{l α −(r+y ), } p Levy dl, and also (noting that on average the fraction 1 2 of particles go to a required direction), we finally find the relations where x ≡ r/l max and˜ ≡ /l max . To obtain the mean field relations for ρ in and ρ out , we use the following argument: In the ordered phase, there are typically multiple particles inside a disk with area πR 2 , 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 l α R. This leads to the relation ρ in ≈ l α πR ρ out . 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 r 2 ρ in + (l 2 max − r 2 )ρ out = ρl 2 max , gives the following relation: Noting that n ≡ n in − n out 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 r * ) is the fixed point of the dynamical behavior of n, i.e., n(r * ) = 0. The average size of the swarmed clusters, indicated by r * , 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 r * and η c , we first note that r * → 0 ⇒ η c → 0 and r * → ∞ ⇒ η c → ∞. Therefore, we expect that η c 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 r * → ∞ in our MF arguments. For the percolation probability, we should seek the conditions that lead to r * → ∞, or equivalently r * /L → 1. 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 η c and r * , which is monotonic increasing function. If one fits this relation using a power-law function, the exponent would be ≈ 0.07, although the range of the quantities are too small to deduce a power-law form (it commonly should be more than one decade).

Geometrical Observables
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.

Density Fractal Analysis
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 blackand-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ρ ≡ N −1 pixels ∑ i ρ i over that sample, where N pixels is the total number of pixels in the system. The filling fraction of the ith box is µ i ≡ N i (δ) N pixels , where N i (δ) is the number of white pixels in the ith box [54,55]. Note that ∑ N box i=1 N i = N pixels , where N box is the total number of boxes. If we attribute a local mass for ith box as m i (δ) ≡ 1 − δ N i (δ),0 and a total mass to the cluster as M(δ) ≡ ∑ i m i (δ), where δ m,n 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 µ i , and is defined as where q is a moment. For scale-invariant systems, Z q scales with δ in a power-law form with the exponent γ q , but the exponent may not be a unique number in all scales: The generalized q-dimension is then defined as so that D f = lim q→0 D q . Note that if one interprets µ i as a probability associated with a small segment (δ) of the system, then D q plays the role of a normalized q-Renyi entropy (Re q (δ)) in the thermodynamic limit δ → 0 so that Therefore, the mass fractal dimension of samples is related to q = 0 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 log δ and not δ d (d = 2 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 q → 1 the fact that relates it to D 1 so that SH(δ)| δ→0 = −D 1 log δ.
Finally the correlation dimension is defined as where where R k is the position of the kth white pixel (not box), and Θ is a step function. It is shown that [54] We generated 10 3 configurations (density snapshots) at the transition points to investigate the anomalous dimensions. Figure 7a shows log Z q (δ) in terms of log(δ) for α = 1.5, and the inset shows γ q 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.

Contour Line Analysis
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] d f = 1 + κ 2 (these extended objects are fractal paths, or sometimes loops). As another example for the Potts models, it is shown that ν(κ) = 2 3(2−κ) , γ(κ) = 4+3κ 2 6κ(2−κ) and β(κ) = 3κ−2 12κ , whereκ ≡ 1/κ is the dual of κ, which is given by the equation √ Q = −2 cos(π/κ) such that, for example,κ = 4 3 for Q = 2 [62].
To investigate this, we start with the fractal dimension of loops (d loop f ) 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 { r 1 , r 2 , . . . , r l } as is then defined by the relation log(l) = d f log(r) + cont., 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 sm. The scale-invariant properties of the distribution function of sm, l, and r are also evident (excluding finite size effects), with P(x) ∝ x −τ x . Here, τ x 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 d loop f is observed to remain constant at approximately 1.40 for all values of α, indicating its robustness. This exponent is associated with κ = 0.8 (in the case of conformal invariance). However, the scaling exponents τ x (x ≡ sm, l, r) 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 Q = 2 Potts (Ising) model, the diffusivity parameters, which are κ G = 3 4 and κ FK ≡κ G = 1/κ G = 4 3 , respectively [62]. It also was previously shown that τ Ising r ≈ 3.4, and τ Ising l = d D f + 1 ≈ 2.5, and also d Geometrical [64]. The exponents that we found for our model for small α values are in agreement with the exponents explored above. By setting κ = 0.8 in our model, we obtain d f (κ) = 1.4, ν(κ) ≈ 0.9, γ(κ) ≈ 1.55 and β(κ) ≈ 0.12. 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 Q = 2 Potts model but is not a perfect fit for this model.

Concluding Remarks
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 p(φ) are in close proximity and difficult to distinguish as shown in the right inset of Figure 2a. Since the values of η c 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 α ≥ 2, 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 Q = 2 Pottes (Ising) model, while the β exponent does not match.

Data Availability Statement:
The research data supporting this publication are provided within this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Crossover to Continuous Transition
In this section, we elaborate on the crossover between first-and second-order transition regimes. In Figure A1, the alterations in the Binder cumulant and the probability density function (PDF) are illustrated as α changes. For the largest α value, the Binder cumulant drops sharply at the transition point to a minimal value as a function of L, signifying the first-order transition. As α decreases, the depth of this drop also decreases, indicating a crossover to a continuous transition. The same trend is apparent in the PDF (the right figure). When α is the largest, the PDF exhibits a well-defined bi-modal distribution, but as α decreases, this bi-modality becomes less distinct. In Figure A2 we represent the results for the dilute phase, i.e., small density regime. We see that in this case, we observe the same results as ρ = 2, i.e., the order-disorder transition which takes place at smaller η values. Figure A2. φ in terms of η (main) and the corresponding data collapse analysis (inset) in the low density limit ρ = 0.3 and α = 1.5. The data collapse analysis shows that η c = 2.098 ± 0.01, β = 0.38 ± 0.1, and ν = 0.9 ± 0.1

Appendix B. Mean Field Theory
In this appendix, we describe some details of the mean field theory. As stated in the text, n in (r) = πρ out 2 l max +r r+ r P(l > r − r)dr where Based on this, we calculate the number of particles that leave this area. For determining n out , we note that the required accumulated probability is given this time by Equation (8) is obtained by inserting these equations into Equations (6) and (7). For the analysis of the average r, we need to determine the difference between n in and n out . Let n be such a quantity n ≡ n in − n out (A4) which counts the rate of change of the average number of active particles inside the swarmed cluster. Then, one can determine the dynamical behaviors of the model in terms of r by studying n.
To be more precise, the average r (which we call r * ) is the fixed point of the dynamical behavior of n, i.e., n(r * ) = 0. In Figure A3, we plot ρ in and ρ out for α = 0.8, α = 1.01, α = 1.25, α = 1.5, α = 1.75, α = 1.95. The graphs for n in (blue line) and n out (orange line) shown are represented in terms of r and α in Figure A4. r * is the point at which these two graphs meet, i.e., the fixed point of the dynamics. The behavior of r * is of central importance in this analysis since the average size of the swarmed clusters shows which phase is stable in which α. n is shown in Figure A5 and the corresponding r * . The resulting graph is Figure 6. Figure A5. n (defined in Equation (A4)) in terms of r for l max = 100 for various α values. r * has been defined as the point n(r * ) = 0 with bold circles with the same color as the main graph (r * decreases as α increases).

Appendix C. Similarities with Q-State Potts Model
This section provides a description of the Q-state Potts model as a framework for second-order transitions. The model exhibits the tricritical point at Q trc = 4, above (under) which the ordered-disordered transition is of the first (second) order. Its relation to the other models is well known; the examples are the O(n) model [62], Coulomb gas [62,65,66], eight vertex model [67], Schramm-Loewner evolution [62], XY model [68][69][70], and conformal field theories [66]. Given that our fractal dimension D loop f ≈ 1.4 for α < α Tcr is consistent with the fractal dimension of the external perimeter of the geometrical clusters of the Ising (D Ising f = 11 8 ), we explore its properties in this section, i.e., the Fortuin-Kasteleyn (FK) clusters. To this end, suppose that the distribution of clusters with m sites is given by P(m) so that P(m) ∝ m −τ exp (−ϑ m) where τ is its exponent and is the decay rate which diverges as T → T c (the critical order-disorder transition temperature), and σ is its corresponding exponent. Then the standard exponents of the transition are given in terms of τ and σ. Importantly, [62] is the "heat capacity" exponent, which is zero for the Ising (Q = 2) model. We do not have this exponent in our model. The two other exponents that we calculated in our model are β (the order parameter exponent) and γ (the fluctuations exponent). The relation between these two exponents and τ and σ is Additionally, the exponent ν (noise renormalization exponent) is given as follows: where d is the dimension of space (here is two). This is the first hyper-scaling relation: The other two exponents (not calculated in our project) are the Fisher exponent and massfractal dimension, which are, respectively, We also have the following hyperscaling relation between β, γ and ν: All of these exponents tell us that 2 < τ < 3. Table A1 shows the exponents for the case Q = 2 and our model. In this table, we see that the best similarity is obtained for Q = 2, which is the critical Ising model. The only difference is concerning the β exponent, which is almost three times bigger than the one for the Ising model. The corresponding diffusivity parameter κ in the SLE theory is given in terms of the basic exponents σ and τ as follows: It is more convenient to represent these exponents in terms of β and γ whereκ ≡ 1/κ as described in the main text.