A kinetic model and scaling properties for non-equilibrium clustering of self-propelled particles

We demonstrate that the clustering statistics and the corresponding phase transition to non-equilibrium clustering found in many experiments and simulation studies with self-propelled particles (SPPs) with alignment can be obtained from a simple kinetic model. The key elements of this approach are the scaling of the cluster cross-section with the cluster mass -- characterized by an exponent $\alpha$ -- and the scaling of the cluster perimeter with the cluster mass -- described by an exponent $\beta$. The analysis of the kinetic approach reveals that the SPPs exhibit two phases: i) an individual phase, where the cluster size distribution (CSD) is dominated by an exponential tail that defines a characteristic cluster size, and ii) a collective phase characterized by the presence of non-monotonic CSD with a local maximum at large cluster sizes. At the transition between these two phases the CSD is well described by a power-law with a critical exponent $\gamma$, which is a function of $\alpha$ and $\beta$ only. The critical exponent is found to be in the range $0.8<\gamma<1.5$ in line with observations in experiments and simulations.

In many of these systems, a peculiar phase characterized by the existence of remarkably large moving clusters has been observed. We designate this phase of collective motion as nonequilibrium clustering. This phase often appears as an intermediate phase between a completely disordered phase, with homogeneous density, and a phase with global orientational order. An experimental case (myxobacteria) and a simulation example (self-propelled rods) are shown in figure 1. The occurrence of large clusters is strongly correlated with a crossover in the shape of the cluster size distributions (CSD) (see figures 1(b) and (d)). Figure 2 illustrates recent findings from simulations of self-propelled hard rods [24,25] and from simulations with a modified Vicsek model (VM) with nematic alignment [14][15][16]. The collective clustering phase is characterized by a non-monotonic CSD with a characteristic peak at large cluster sizes, see CSD for high densities in figure 1. At the onset of the collective clustering phase, the CSD follows a power law with a characteristic exponent. It has been observed in SPP experiments [5,6] and simulations [17,18] that the CSD can be power law distributed with the exponent typically varying between 0.85 and 1. 35. While the range of exponents indicates the absence of a universal scaling at the onset of the collective clustering phase, it is not yet clear what determines the actual value of the exponent.
Note that the collective clustering phase can be observed in the absence of global orientational order, as reported in myxobacterial experiments [6] and self-propelled rod simulations [17] (figure 1), as well as in the presence of global orientational order [11][12][13]26]. The connection between the collective clustering phase and the onset of orientational order, in a system with proven genuine long-range (orientational) order, is not known. Empirically, we observe that the onset of the collective clustering phase is typically found in between the disordered and the ordered phases, as illustrated in figure 2. We note that most investigations with where it has been found that the cluster size distribution (CSD) is a function of the cell density as shown in (b) (see [6]). (c) Snapshot of moving clusters in self-propelled rod simulations, where the CSD is known to be a function of the particle density, as illustrated in (d) (see [17]). The onset of the collective (clustering) phase (blue) is often found in between a disordered and an ordered (orange) phase. The control parameter depends on the model we look at: in (a), η denotes the fraction of covered area, while in (b), is angular noise. Panel (a) corresponds to simulations with self-propelled hard rods in two dimensions with aspect ratio κ = 13 [24,25]. Panel (b) shows a modified Vicsek model with nematic alignment with density ρ = 0.25 [14]. The noise value for the transition to global order was taken from a simulation with >10 6 self-propelled particles [15], while the noise value for the onset of clustering in a system with 16.384 self-propelled particles is based on the results of [16]. a focus on the cluster statistics (both in experiments and simulations) were carried out with small to intermediate numbers (around 10 2 -10 3 ) of self-propelled particles, whereas investigations on the onset of global order were conducted with rather large numbers (> 10 5 ). Hence, little is known about the onset of the collective clustering phase in systems (and simulations) with a large number of particles.
Here, we focus on the emergence of the collective clustering phase in the absence of global order, i.e. we assume that cluster velocities are uncorrelated. We propose a kinetic clustering theory based on the approach presented by us in [17], and extend our earlier treatment to coagulation and fragmentation kernels that depend, respectively, on the scaling of the cluster cross-section with cluster 'mass', i.e. the number of particles in the cluster, and the scaling of cluster perimeter with cluster mass. The scaling of the cluster cross-section and of the cluster perimeter are characterized by exponents α and β, respectively. Furthermore, we present a comprehensive analysis of this kinetic cluster model, including a finite size (FS) study of it. The FS study of the kinetic model reveals that: (i) the transition to the collective clustering phase-characterized by a non-monotonic CSD-is generic to all SPPs with either a polar or a nematic alignment mechanism, and (ii) there exists a critical asymptotic CSD exponent γ which is a function of α and β only. By studying the physically meaningful parameter space α-β, we found that γ always falls in the range 0.8 < γ < 1.5. Note that in the absence of cluster-cluster correlations, this transition marks the onset of collective motion as defined in [6]. Although the simplified clustering theory is, strictly speaking, only valid in the absence of global order, we show through simulations-using a large number of particles-that a comparable picture holds when global order is observed. This paper is organized as follows. In section 2 we present our simple kinetic clustering model under the assumption that there are no cluster-cluster correlations, and show that for any finite system there are two clustering phases. We perform a system size analysis of the kinetic model in section 3, where we report on the scaling properties of the obtained CSD, which we found to depend on the scaling of the cluster cross-section and cluster perimeter. In section 4 we study cluster formation in an SPP model that is known to exhibit cluster-cluster correlations, and contrast the obtained results with our cluster-cluster uncorrelated clustering theory. The implications and limitations of the proposed kinetic approach are discussed in section 5.

A kinetic model of clustering
We look for a description of the clustering process in terms of n 1 (t) , n 2 (t) , . . . , n N (t) , where n i (t) represents the average value at time t of clusters formed by i particles and N the number of particles in the system. To ease the notation, n i (t) is hereafter referred to as n i (t). Our strategy consists in deriving a coagulation Smoluchowski equation with fragmentation for these n i (t) objects. For an introduction to this kind of coagulation equation, see [27]. In SPP systems, the only conserved quantity is the number of particles N . Neither the overall orientational order nor the number of clusters is conserved. Thus, our generalized coagulation Smoluchowski equation should conserve the number of particles. We start by simplifying the clustering dynamics. Because we are dealing with self-propelled particles, we assume that clusters are the result of either an explicit or an effective alignment mechanism such that particles inside a cluster move coherently in the same direction. This implies that clusters move at speeds comparable with that of the individual particles, independently of the size of the cluster. Below we also discuss what can be expected if this condition is relaxed. In contrast, clusters of passive particles, driven by thermal fluctuations, are such that their mobility decreases as a function of their size. We stress that the assumption of a size-independent cluster speed evidences the non-equilibrium nature of our simple cluster theory. We further assume that (binary) collisions among clusters may result in cluster-cluster fusion, and neglect the possibility of cluster fragmentation induced by cluster-cluster collisions. This assumption is justified as long as the dynamics is overdamped, and is in accordance with what is observed in experiments with gliding bacteria [6] and simulations with self-propelled rods [17]. We simplify the cluster fragmentation dynamics by assuming an evaporation-like process by which clusters shrink in size by losing one by one those particles that are on the cluster boundary. In summary, we are assuming an irreversible clustering process in which a cluster of mass j can undergo the following 'reactions': Note that because the reaction C j+k → C j + C k does not occur, the process is, in this sense, irreversible. Now, we look for a description of the process in terms of the (average) number n j (t) of clusters with j particles (i.e. the number of C j ) at time t, whose time evolution takes the forṁ for j = 2, . . . , N − 1, where the dot denotes the time derivative, B j represents the rate at which a cluster of mass j loses particles, defined as and A j,k is the collision rate between clusters of mass j and k, defined by where L 2 is the area of the two-dimensional (2D) space where particles move and v c is the cluster speed, which we assume to be v c ∼ v 0 with v 0 the speed of individual particles. Note that equations (2) are such that mṅ m = 0 and thus the number of particles N = m n m (t) is conserved. In equation (3), d 2 /D is the typical time a particle located on the cluster boundary needs to detach from a cluster, with d the maximum distance two particles can be apart to still be considered as connected and D the diffusion coefficient with respect to the center of mass of the cluster. We stress that in equation (2) we have assumed that clusters are uncorrelated. Several statistical features of the model can be given explicitly if one assumes that the SPPs follow a dynamics where their translation is determined by a constant speed and the direction of motion is subject to alignment interactions with neighboring SPPs and an angular noise with amplitude η. An example of such dynamics is given, in two dimensions, by the following equations of motion: where j is the index of the particle, θ t j defines (in two dimensions) the particle moving direction, x t j denotes the position at time t, the sum is taken over all particles within a unit distance of j, which includes j, and η t j is a uniformly distributed random variable such that −η/2 < η t j < η/2, η t j = 0 and η t j η t k = (η 2 /12)δ t,t δ j,k . The dynamics given by equation (5) is frequently used in simulation studies, e.g. in [10-16, 28, 29] (for a review see [1]). More specifically, f (e iθ t k , e iθ t j ) = e iθ t k defines a polar alignment rule as in the VM [10], ]e iθ t k defines SPP with a nematic alignment rule as used [14,15]. Assuming that the SPP obey equation (5), the diffusion around the center of mass of a polar cluster of SPP can be approximated, as detailed in [16] where v 0 is the speed of individual particles and t sets the discrete time step for typical SPP dynamics such as the one given by equation (5). By making a Taylor expansion we find that D ∼ v 2 0 /12η 2 t. Note that for non-interacting random walkers subject to the same kind of noise, D ∝ v 2 0 /η 2 . In addition, we have to consider that the splitting rate B j has to be proportional to the number of particles on the cluster boundary. In equation (3), we have assumed that the perimeter l of a cluster scales with its mass j as l ∝ j β , with β constrained to 1/2 β 1. Let us provide a physical context to this assumption, always assuming that we are in a 2D space. Surface tension would tend to minimize cluster perimeters, and so clusters would be round and β = 1/2, as observed in liquid-vapor drops [30]. On the other hand, for random deposition of particles as in classical percolation, β ≈ 1. If clusters are 'chains' of particles, then also β = 1. In short, we expect 0.5 < β < 1, with β = 0.5 representing the most compact way we can arrange j particles to minimize the number of particles at the cluster interface, and β = 1 representing the opposite extreme, i.e. a situation where all particles are at the cluster boundary.
The collision rate A j,k is derived in analogy to the collision rate in kinetic gas theory [31], which means that we assume that clusters move ballistically in between collisions, which leads to a collision rate proportional to the sum of the scattering cross-section of the involved clusters, the relative cluster speed and the cluster density. Due to the latter assumption, A j,k ∝ L −2 . The relative cluster speed can be assumed to be not too different from that of the individual clusters, v c . In SPP systems, we often find that the individual cluster speed is close to that of individual particles, i.e. v c ∼ v 0 . If v c exhibits a dependence on cluster size of the form v c ∝ j −ω , the exponent ω can be absorbed into α. For simplicity, here we assume that v c ∼ v 0 . For instance, the speed of SPP clusters can be approximated by The scattering cross-section of a 2D cluster is the (average) projection of the cluster on a given axis. Thus, the scattering cross-section has to be always less than or equal to the cluster perimeter. If the scattering cross-section of a cluster scales as ∝ σ 0 j α , the previous observation implies that α β is the only meaningful physical scenario. Below, we will see that there are two qualitatively very different physical scenarios, α = β and α < β. Finally, q s represents the probability that a cluster-cluster collision resulted in a successful fusion of clusters. Here, we assume it is a constant but certainly it can depend on various (intensive) variables of the actual system. In particular, it depends on the alignment symmetry. A very rough assumption would be to assume that if q s = q 0 for the ferromagnetic alignment, q s = q 0 /2 for the nematic one. Equations (2) are scaled and transformed into dimensionless form by dividing equations (3) and (4) by D/d 2 , which leads to the following dimensionless parameter: Thus, equation (3) reduces to B j = j β , and equation (4) to A j,k = P( j α + k α ). The parameter P controls the relative weight of fragmentation with respect to coagulation. For the 2D SPP model type considered here, we obtain P ∝ (Lη) −2 . Note that while in a three-dimensional space we can expect A j,k and B j to follow a similar scaling as the one proposed above, one-dimensional flocking models are essentially different, exhibiting a homogeneous density phase and a condensed phase [32]. Finally, computation of the critical point P can be done, alternatively, by directly studying the stability of the individual phase and ignoring the actual shape of the CSD, as recently proposed in [33].
Except for a few special cases, for instance α = β = 0, equation (2) cannot be solved analytically. Nevertheless, direct numerical integration of equation (2), for instance through a fourth-order Runge-Kutta method, can be easily implemented. Note that the numerical implementation has to ensure particle number conservation at every integration step, which often requires a renormalization checkpoint in the integration loop. We integrate the dimensionless version of equations (2), with the initial condition n j (t = 0) = N δ j,1 , and find that the (weighted) CSD, defined as reaches a steady state, i.e. p(m, t → ∞) = p(m; P). This also implies that the number of clusters in the system, M(t) = N m −1 p m (t), also reaches a steady-state value as t → ∞.
In the literature the term CSD is frequently used to refer top(m, t) = n m (t)/M(t). The disadvantage of this definition is that its normalization constant, M(t), varies with time. Whilẽ p(m) refers to the probability of finding a cluster of size m, p(m) indicates the probability of a randomly selected particle to be in a cluster of size m. Given α and β, the CSD p(m) depends on N and the value of P. Figure 3 summarizes the clustering behavior with P for a given N . Note that there exists a critical P c (N ) that separates two different types of clustering behavior. For small values of P, the distribution p(m) is monotonically decreasing, dominated by an exponential tail. As P → P c , p(m) approaches a power law, with a system size cutoff. For P > P c , p(m) is non-monotonic and exhibits a peak at large cluster sizes. This dramatic change of behavior at P c unveils a phase transition, which is evidenced by figures 3(c) and (d) that show m = m p(m) and σ 2 = (m − m ) 2 p(m), respectively. Below, we will see that this transition is a genuine phase transition-in the thermodynamical sense-only for α = β.

Scaling properties
To study the asymptotic behavior of the CSD, equation (7), we need to solve the set of N coupled equations given by equation (2). These equations depend on the dimensionless parameter P, the  (2). We conduct a FS study of equation (2) to extrapolate the behavior of p(m) for N → ∞. Let us recall that the thermodynamical limit is given by N → ∞ as well as L → ∞, with the constraint that ρ = N /L 2 is a constant. In the previous version we have learned that P ∝ L −2 . In this section, we show that the critical point P c (N ; α, β)-defined in the previous section-depends on N . For every N , we numerically estimate P c (N ), using a bisection method, as the point at which the CSD p(m) is no longer monotonically decreasing. We find that at the critical point P c (N ), p(m) scales as where γ is a critical exponent that depends on α and β. This finding results from the evident collapse of the various system size curves at the critical point shown in figure 4(a). The figure indicates that for a given N , p(m; P c ) follows the scaling given by equation (8) up to a certain cluster size m * (N ) above which the FS of the system becomes evident. The exponent γ results from extrapolating the behavior for N → ∞ and it is numerically obtained by minimizing the error defined as W (γ ; N ) = N m m γ p(m; P c ). This means that to obtain γ as a function of α and β (figure 5), for every pair (α, β) a FS study, such as the one described above, has been realized. This implies first the estimation of P c for every system size N before computing the actual value of γ . Below the critical point, the following scaling is obeyed by p(m): where K is a constant. This is confirmed by the obtained collapse of the various system size curves where P was rescaled as indicated in equation (9), see figure 4(b). Below we explain that ξ is connected to the scaling of P c with N . For α = β, we find that K = 1, while for α = β, K depends on b, N and β. All this means that if, for instance, we double the system size and reduce P by 2 −ξ , we fall on the same distribution p(m) for α = β, while for α < β there is a multiplicative constant K . For P < P c , we find that the CSD is well fitted by p(m) ∼ m −γ exp(−m/m(P)). For P > P c , we observe that there is a shift with N of the peak that emerges at large cluster sizes ( figure 3(b)). The scaling given by equation (9) works up to a given size m m (N ) which is given by the minimum of p(m). We are interested in knowing the behavior of the peak with the system size N . To answer these questions, we study the scaling of Q(N ) = ∞ m m (N ) p(m). We find that for α = β, Q(N → ∞) → 1. This indicates the presence of a collective phase in which most of the particles are part of large clusters. Note that this does not mean the formation of a single giant cluster; the width of the peak does not shrink to zero. For α < β, Q does not converge to 1 as N → ∞. This is because for α < β, in the limit of N → ∞ the transition does not occur. Nevertheless, for any finite N , we always observe the transition to the collective phase. This observation can be understood by studying the system size behavior of P c (N ). We At the critical point P c , the CSD p(m; P c ) follows, asymptotically with N , a power law characterized by an exponent γ , equation (8). The critical exponent γ is a function of the exponents α and β, associated with the scaling of the coagulation and fragmentation kernels, A j,k and B j , respectively, with the corresponding cluster sizes. Panel (a) shows the behavior of γ for α = β, while panel (b) corresponds to β > α. The scaling of the critical point P c (N ) with the system size N follows a power law with exponent ξ as indicated by equation (10).
The value of ξ also depends on α and β. The inset in (b) shows the behavior of the exponent ξ as a function of β for α = 0.5. Note that while ξ = 1 for α = β, ξ < 1 for β > α. This result turns out to be extremely important to understand the behavior of the system in the thermodynamical limit, see the text.
find that with ξ = 1 for α = β and ξ < 1 for α < β. For a fixed value of α, ξ becomes a function of β as illustrated in the inset of figure 5(b) for α = 1/2. In the thermodynamical limit, with N → ∞ and L → ∞, while ρ = N /L 2 is constant, we can estimate the critical point. Combining equations (6) and (10) we obtain, for a given set of parameters D, v 0 , σ 0 and d, the critical density ρ c above which a peak at large cluster sizes emerges: From equation (11) it is clear that when ξ < 1, ρ c diverges with the system size N . Only for ξ = 1 is ρ c a finite quantity in the thermodynamical limit. Thus, for α = β there is a critical density ρ c above which the collective transition occurs. In contrast, for α < β, ρ c → ∞ and any finite density ρ 0 corresponds to the individual phase. Equation (11) predicts the critical noise intensity η c below which one expects the collective clustering phase to occur as η c ∝ ρ N ξ −1 , from which it follows that for ξ = 1, the critical noise scales as η c ∝ √ ρ. Interestingly, a similar scaling was reported for the onset of collective motion in simulations with the VM [10,26,29].

Clustering in the presence of global orientational order
The simple set of Smoluchowski-type equations (2) has been derived assuming that clusters move in an uncorrelated fashion, as observed for instance in bacterial experiments [6]. However, Figure 6. Simulation snapshots of the Vicsek model (VM) for various noise intensity values. Particles belonging to the same cluster are labeled with the same color. As can be seen in figure 7, for η = 1.0 the system is disordered, while for η < 1 there is a preferred direction of motion. Although the snapshots for η = 1.0 and η = 0.9 seem comparable, their CSD is remarkably different: an exponential distribution for the former and a power law for the latter. Traveling bands in the VM are composed of many correlated clusters (η = 0.7). Finally, for very low η values, most particles form part of very large clusters (η = 0.05).
several SPP models exhibit long-range orientational order and macroscopic structures such as bands, which indicates that there are strong correlations among cluster velocities. Thus, a priori we cannot expect equation (2) to describe quantitatively the cluster dynamics of such systems deep in phases with strong cluster-cluster correlations. Interestingly, it has already been observed that the VM in its (global) ordered phase displays CSDs that can be power law distributed, with γ roughly in the range 0.8 < γ < 1.3, according to [11,12,26]. Here, we carried out simulations in the VM, i.e. equations (5) with a choice of f = e iθ t k . Our simulation results confirm that indeed the VM exhibits power law distributed CSDs with exponents in this range, and show that this occurs close to the well-known disorder-order transition of the VM. Furthermore, we found that the cluster statistics in the VM, below and close to the critical point-i.e. when cluster-cluster correlations are weak or absent-resembles that obtained with the kinetic clustering model. As expected, far away from the critical point and well in the ordered phase, cluster-cluster correlations induce effects that cannot be accounted for by the kinetic clustering model. In short, the simulation data strongly suggest that close to the disorder-order critical point, the VM exhibits a transition to a collective clustering phase like the one described with the kinetic model.
The order-disorder transition in the VM is observed when the noise intensity η is decreased below a critical value η * . The orientational order is characterized by φ = | j exp(iθ t k )|/N , which defines the usual polarization order parameter using complex number notation [10]. In the disorder phase, the clusters cover the space homogeneously ( figure 6) and the CSD is dominated by an exponential tail, see figure 7. If η is decreased below η * , velocity-velocity correlations among particles become important, φ increases and clusters grow significantly in size. Below the critical point η * , the CSD is a power law with an exponent in the range [0.8, 1.3], which results also in a change in behavior of m /N , as shown in figure 7. The power law distributed CSD indicates that the system displays arbitrarily large clusters. Indirectly, this also means that particle-particle correlations are long-ranged, i.e. arbitrarily large like the cluster sizes. Note The average (normalized) cluster size m /N as a function of η is shown in (c), and the orientational order parameter φ versus η in (d). While the system is disordered, the CSD is dominated by an exponential tail. Below the critical point, the CSD is a power law with an exponent γ (η). Close to the critical point γ = 1.3 and when η is decreased γ approaches 1. For very low η values, the CSD decreases even more slowly than 1/m (see figure 6).
that this does not necessarily imply the existence of cluster-cluster correlations. Nevertheless, in the VM at low values of the noise intensity, cluster-cluster correlations become evident.
Arguably, due to these cluster-cluster correlations, we observe a clear deviation from what the theory predicts for η η * . Moreover, in the VM, γ is a function of η. We observe that close to η * , γ ≈ 1.3, and that when η is decreased, γ approaches 1. It is close to η * that we observe the emergence of a traveling band. Note that this band is not a connected component, but a cloud of highly correlated moving clusters. We find that the CSD is power law distributed in the band regime. For η η * the system is highly ordered and clusters move roughly in the same direction, which means that cluster collisions are less frequent. On the other hand, low values of η also imply that spreading of a cluster around its center of mass is also small. In summary, at low noise values the cluster dynamics slows down significantly in the VM, which results in a non-monotonic response of m /N with η ( figure 7). Interestingly, we observe that at extremely low η values the CSD decreases even more slowly than 1/m and m /N experiences a sharp increase.

Concluding remarks
We have argued that SPP systems exhibit for a finite system size two phases: an individual phase, where the CSD is dominated by an exponential tail, and a collective (clustering) phase characterized by a non-monotonic CSD. Assuming that the moving directions of clusters are uncorrelated, equation (2) justifies the existence of these two phases. At the transition point, the CSD is a power law characterized by an exponent γ that depends, according to the proposed kinetic clustering model, on the scaling with cluster size of the cluster crosssection-characterized by an exponent α-and the cluster perimeter-characterized by an exponent β. A systematic study of the α-β parameter space revealed that γ always falls in the range between 0.8 and 1.5, which is consistent with experimental observations [5,6] and simulations [17,18]. However, we learn that in the thermodynamical limit, the above-mentioned transition occurs only for α = β, while otherwise, the system remains in the individual phase. It is worth pointing out that for the special case α = β the kinetic model predicts a critical density value that is independent of the system size. Hence, the properties for the transition to nonequilibrium clustering found in the kinetic model for intermediate values of N hold true in the thermodynamic limit.
In addition to the study of the kinetic clustering theory given by equation (2), we characterized the cluster statistics of the VM [10]. Despite the fact that cluster-cluster correlations and global order emerge in the simulations of the VM that are not treated in equation (2), we found that a transition to a collective clustering phase is also present in the VM. Moreover, γ falls in the expected range if simulations are performed with ferromagnetic alignment, as shown here and in [11,12,26], as well as with nematic alignment [16]. Notably, we found that the CSD is power law distributed for a range of noise intensity values η close to the order-disorder critical point, with γ a function of η. In summary, the existence of two clustering phases is found as well in SPP systems in the presence or absence of global orientational order.
Our kinetic approach is complementary to the ongoing effort at deriving coarse-grained nonlinear field equations for the description of active matter [34][35][36][37][38][39][40][41][42] (for a review, see [2,43]) that are better suited to capture large-scale structures or the emergence of global orientational order in large systems. In contrast, the observed clustering phenomena are typically studied in the absence of long-range order, on a small scale, and for intermediate system sizes. It is also worth noting that the non-equilibrium clustering phase leads to (apparent) giant number fluctuations (GNF) [6], which often have similar properties as the ones predicted by Toner and Tu [34] in the phase with long-range orientational order. The kinetic clustering theory indicates that the CSD can be power law distributed, in the thermodynamical limit, in the absence of global orientational order. Since a power law distributed CSD necessarily leads to apparent GNF, the theory supports recent observations regarding the presence of GNF in the absence of orientational order, see [6,45]. Moreover, in section 4 we showed that in large-scale simulations, the CSD exhibits fat tails in the presence of global order. This certainly leads again to GNF. While the presence of fat tails in the CSD leads to apparent GNF, the contrary is not necessarily true. A solid connection between GNF and fat tails in the CSD remains to be established.
Finally, we expect a similar collective (clustering) transition and cluster dynamics to occur in all SPP systems where the particle speed is not strongly affected by the local density. The presence of a density-dependent speed, as in [44][45][46][47], can dramatically change the above-given clustering picture, since in these models, large clusters, when they are not strongly aligned, are prone to slowing down significantly. As result of this effect, the SPPs can form-when the alignment mechanism is either weak or inexistent-a single large cluster which coexists with a background gas of particles in the individual phase, as found in [44][45][46][47]. The cluster dynamics of these systems is likely to be related to an equilibrium-like phase separation, as suggested in [44,46].