Fernando Peruani

Lutz Schimansky-Geier

Markus Bär

Cluster dynamics and cluster size distributions in systems of self-propelled particles

Systems of self-propelled particles (SPP) interacting by a velocity alignment mechanism in the presence of noise exhibit a rich clustering dynamics. It can be argued that clusters are responsible for the distribution of (local) information in these systems. Here, we investigate the statistical properties of single clusters in SPP systems, like the asymmetric spreading of clusters with respect to their moving direction. In addition, we formulate a Smoluchowski-type kinetic model to describe the evolution of the cluster size distribution (CSD). This model predicts the emergence of steady-state CSDs in SPP systems. We test our theoretical predictions in simulations of SPP with nematic interactions and find that our simple kinetic model reproduces qualitatively the transition to aggregation observed in simulations.

Examples of large-scale self-organized patterns in systems of self-propelled particles with short-range interactions are found at all scales, from group of animals [START_REF]Three Dimensional Animals Groups[END_REF][START_REF] Cavagna | Proc. Natl. Acad. Sci[END_REF][START_REF] Bhattacharya | [END_REF] and human crowds [4] down to insects [5,6], bacteria [START_REF] Zhang | Proc. Natl. Acad. Sci[END_REF], or actin filaments [START_REF] Schaller | [END_REF]. Such patterns are also found in non-living system like in driven granular media [9,10,11,12]. Despite the fact that the interaction mechanisms between individual elements are of a different nature, it is possible to determine some common requirement to achieve large-scale (spatial) self-organization. Particularly important for the emerging macroscopic patterns are the self-propulsion of the agents, and their velocity alignment mechanism. Such velocity alignment mechanism is often the result of purely physical forces which give rise to an effective alignment interaction as observed in systems of self-propelled rods interacting by volume exclusion interactions [10,13,14,15,16]. In other cases, the alignment mechanism emerges from the complex behavior of the moving entities like in birds [START_REF] Cavagna | Proc. Natl. Acad. Sci[END_REF]. Simple individual-based models like the Vicsek model [17] have helped to reveal the relevance of these two elements, self-propulsion and alignment, by reducing the problem to the competition between a local aligning interaction and some noise [18]. In two dimensions, self-propelled particles moving at constant speed with a ferromagnetic-like velocity alignment exhibit at low noise a phase characterized by true long-range polar order which translates into a net flux of particles [17,19]. This ordered phase exhibits several remarkable features like the spontaneous formation of elongated high density bands that move at roughly constant speed in the direction perpendicular to the long axis of the band, and anomalous density fluctuations for low noise levels [19,20]. When the alignment is replaced by a nematic velocity alignment, particles display a phase characterized by true long-range nematic order at low noise intensity [21]. Interestingly, spontaneous density segregation into bands is also observed for these particles, for both, ordered and disordered phase [21]. These bands are formed by particles moving in opposite directions along the long axis of the band. In the ordered phase, bands are wide and straight, while bands in the disordered phase are thinner and highly fluctuating, leading to a zero average global nematic order. For very low noise, bands dissappear and an ordered phase with anomalous density fluctuations is observed [21].

The two examples we discussed above, self-propelled particles with ferromagnetic and nematic velocity alignment, display large-scale high-density patterns and seggregation in absence of any attracting force [19,21]. Interestingly, this also occurs at densities below the percolation threshold. Clearly, a prerequisite for the emergence of ordered phases is the establishment of some sort of effective interaction/communication among the particles. Such interaction/communication should allow information on the local order to be distributed across the system in order to reach global order. On the other hand, this local information is transported by the particles themselves as they move. At this point it is important to notice that local (velocity) alignment leads to local polar order which results in the formation of clusters of particles that move in the same direction. In other words, information is often rather transported by these clusters than by isolated particles, as discussed in [21]. Altogether, the cluster dynamics of self-propelled particles plays a key role in the formation of globally ordered phases.

Here, we study first the properties of single clusters and second propose a simple kinetic theory that describes the steady state cluster size distribution of self-propelled particles. Our approach is a modification of the Smoluchowki kinetic equations that were developped to describe the aggregation of colloids [22]. Clustering effects and the emergence of steady state cluster size distributions in self-propelled particle systems were observed in the Vicsek model [23], selfpropelled rods [13], sperm-like swimmers [24], and swimming [START_REF] Zhang | Proc. Natl. Acad. Sci[END_REF] and gliding bacteria [25]. While transient clusters and aggregation were reported in colloidal self-propelled rods [26] and in a 1d-system of self-propelled particles [27]. Our theoretical approach aims at finding a common mechanism behind these scattered observations. We start by defining the equation of motion of generic self-propelled particles in section 1. Then, we focus on the statical properties of single clusters in section 2, to finally propose a simple set of equation for the evolution of the cluster size distribution, section 3. A summary of the obtained results is given in section 4.

Evolution equation of self-propelled particles

We consider point-like particles moving at constant speed in a two dimensional space, as proposed in [14], and assume an over-damped situation such that the state of particle i at time t is given by its position x i and its direction of motion θ i . The evolution of i-th particle is governed by the following equations:

ẋi = v 0 v(θ i ) (1) θi = -γ(v 0 ) ∂U ∂θ i (x i , θ i ) + ηi (t) (2)
where γ(v 0 ) ∝ v -1 0 is a relaxation constant, and U the interaction potential between particles, and hence ∂U ∂θi (x i , θ i ) defines the velocity alignment mechanism. The parameter v 0 represents the active speed of the particles, v(θ i ) is defined as v(θ i ) = (cos(θ i ), sin(θ i )), and ηi (t) is an additive white noise applied to the direction of motion.

In analogy to spin systems, a ferromagnetic velocity alignment mechanism is given by a potential defined as:

U F (x i , θ i ) = - |xi-xj|≤ǫ cos(θ i -θ j ) ( 3
)
where ǫ is the interaction radius of the particles. For the nematic alignment mechanism, the potential takes the form:

U LC (x i , θ i ) = - |xi-xj|≤ǫ cos 2 (θ i -θ j ) . (4)
One can add a coupling strength coefficient in front of the sum of Eqs. ( 3) and (4). We assume that the coupling strength is absorbed in γ(v 0 ) in Eq. ( 2). Notice that the potential given by Eq.

(3) exhibits one minimum, while Eq. ( 4) has two minima. The latter situation corresponds to particles moving in opposite directions.

In the limiting case of very fast angular relaxation we obtain from Eqs. ( 1) and ( 2) the updating rules:

x t+∆t i = x t i + v 0 v θ t i ∆t (5)
θ t+∆t i = arg    |x t i -x t j |≤ǫ f (v(θ t j ), v(θ t i ))    + η t i (6)
where arg (b) indicates the angle of a vector b in polar coordinates, and η t i is a delta-correlated white noise of strength η (η t i ǫ -η 2 , η 2 ). Given two vectors a and b, f (a, b) is defined as follows. For ferromagnetic alignment, f (a, b) = a and Eqs. ( 5) and ( 6) reduce to the Vicsek model [17]. For nematic alignment, f takes the form f (a, b) = sign(a.b) a and Eqs. ( 5) and ( 6) define a minimal self-propelled rod model [14,21].

Order parameters

Ordered phases can be characterized by the following order parameters. The ferromagnetici.e., polar -order parameter is defined by:

S F = 1 N N k=1 exp(i θ k ) , (7)
where N stands for the total number of particles in the system, and the direction θ k is represented as a phase in the complex plane. This definition is equivalent to S F = 1 N N k=0 v(θ k ) . S F takes the value 1 when all particle move in the same direction, while in the disordered phase, i.e., when particles move in any direction with equal probability, it vanishes.

On the other hand, the nematic ordered parameter takes the form:

S LC = 1 N N k=1 exp(i 2 θ k ) . (8)
Formally, S LC can be derived from the order parameter matrix Q of liquid crystals (LC) [START_REF] Doi | The Theory of Polymer Dynamic[END_REF], as the largest eigenvalue (which here we have normalized such that S LC ∈ [0, 1]). When the system is perfectly nematically ordered, i.e., when particles move in opposite directions along the same axis, S LC takes the value 1. A genuine nematic phase implies half of the particles moving in one direction, and the other half in the opposite direction. If this symmetry is broken, for instance, by having 3/4 of the particles moving in one direction, and 1/4 in the opposite, S LC will still be S LC = 1, but S F > 0. In summary, a perfectly polarly ordered phase is characterized by S F = S LC = 1, while for genuine nematically ordered phase, S F = 0 and S LC = 1.

Evolution of a single (isolated) cluster

This section is mainly devoted to the understanding of the evolution of a single, isolated cluster. We consider a situation in which initially all particles are located at the origin and move in direction + x with speed v 0 . Since initially each particle can see all the others, the problem can be described initially by a mean-field. According to Eq.( 6), and assuming ferromagnetic interactions, particles calculate the same common direction of motion. The additive noise η t i acts just as a perturbation around the global common direction of motion. This average vector can be thought as an external field that guides the particles. Thus, particles perform a directed random walk [START_REF] Peruani | [END_REF]. Since initially each particle can see all the others, all of them calculate the same average direction α 0 . The initial condition in our example is such that all particles move in direction + x at t = 0, and thus α 0 = 0. This assumption is only true for η < π, if particles interact by nematic alignment. Then, the angular dynamics of the ith-particle can be approximated by:

θ t+∆t i = α 0 + η t i , (9)
where again, for our example, α 0 = 0. This means that the probability of finding a randomly chosen particle pointing in direction θ at t + ∆t is

P (θ) = 1 η g(θ, α 0 , η) , (10)
where g(θ, α 0 , η) is defined to be 1 when α 0 -η/2 ≤ θ ≤ α 0 + η/2, and 0 otherwise, i.e.,

g(θ, α 0 , η) = H(α 0 + η/2 -θ).H(η/2 -α 0 + θ)
, where H(x) is a Heaviside function. In Eq. ( 9) we are assuming that the average direction α 0 does not change over time, which is not true for a finite cluster. Nevertheless, this simplification makes the problem analytically tractable, and, as we will see, the obtained results will be of interest. Under these assumptions, the position of the i-particle at time t n = n∆t -using the discrete time description, i.e., Eqs. ( 5) and ( 6) -is given by:

x i (t n ) = n k=0 cos(θ i (t k ))v 0 ∆t (11)
y i (t n ) = n k=0 sin(θ i (t k ))v 0 ∆t . (12)
From Eq. ( 11) and (12), and using Eq. (10), it is possible to derive x(t n ) , y(t n ) , x 2 (t n ) , and y 2 (t n ) , where ... denotes an average over all particles and realizations of the noise. For the calculation, it is important to assume that θ i (t k )θ j (t l ) = σδ k,l δ i,j , with σ the second moment of distribution given by Eq. (10). The calculation becomes straightforward after computing cos(θ i (t k )) cos(θ i (t k )) , cos(θ i (t k )) cos(θ i (t l )) , sin(θ i (t k )) sin(θ i (t k )) , and sin(θ i (t k )) sin(θ i (t l )) . With this at hand, we derive the diffusion coefficients of the spreading process. To gain intuition on the problem, we express the time evolution of the cloud of particles in terms of the (continuum time) particle density ρ(x, t) which obeys the following equation (assuming propagation of the particles along + x, and thus α 0 = 0):

∂ t ρ(x, t) = -V (η, v 0 )∂ x ρ(x, t) + ▽. (D(η, v 0 )▽ρ(x, t)) . (13)
In the convective term, V (η) is the mean projection of the instantaneous velocity of the particles on the + x semi-axis. V (η) is by definition V (η) = v 0 2π 0 dθP (θ) cos(θ), that takes the form:

V (η) = v 0 2 sin(η/2) η . ( 14
)
This is nothing else than lim tn→∞ x(t n ) /t n . Notice that in the limit of η → 0, V (η) → 1, and thus there is a deterministic transport of particles without any diffusion. In the limit of η → 2π, V (η) → 0, and particles experience diffusive motion and no convective flux. Now we focus on the diffusive term of Eq.( 13), where the diffusion matrix D takes the form: Symbols correspond to measurements of A(t) in simulations with N = 1000 self-propelled particles with nematic alignment, which initially were located at the origin and all pointing in direction + x. The dashed lines correspond to the approximation given by Eq.( 18)

D(η, v 0 ) = D x (η, v 0 ) 0 0 D y (η, v 0 ) . ( 15
)
with D x and D y defined as:

D x (η, v 0 ) = 1 2 - sin(η/2) η/2 2 + sin(η) 2η v 2 0 ∆t (16)
D y (η, v 0 ) = 1 2 - sin(η) 2η v 2 0 ∆t (17)
These terms are obtained from lim tn→∞ ( x 2 (t n ) -x(t n ) 2 )/t and ( y 2 (t n ) -y(t n ) 2 )/t n . Notice that for small noise η, D y > D x . This implies that the cluster spreads more in the direction orthogonal to the common moving direction. As mentioned above, particles with ferromagnetic alignment form macroscopic bands that move in the direction perpendicular to the longest axis of the band. Eqs.( 16) and (17) indicate that at the level of an individual (isolated) cluster a similar situation arises: the cluster moves in the direction perpendicular to its longest axis. Notice that ∆t in Eqs. ( 16) and ( 17) plays the role of the inverse of the turning rate, i.e., α -1 in Eq. ( 10) and Eq. ( 11) in [START_REF] Peruani | [END_REF]. The description given by Eq. ( 13) is valid while the cloud of particles remains one coherent giant cluster. For longer times, this picture fails and particle motion looses its coherence. The interaction among particles is such that particles move in the same direction as long as they can see each other. It cannot, however, prevent the particles from slowly moving apart due to the fluctuations the direction of motion.

In the following, we compute the time evolution of the spreading coefficient A of the swarm around its center of mass. The spreading area is defined by

A(t n ) = x(t n ) 2 -x(t n ) 2 .
This quantity, A(t n ), can be derived directly from the discrete (time) process by making use of Eqs. ( 11) and (12) as explained above. Alternatively, the continuum time equivalent A(t) can be obtained from Eq.(13). To obtain x(t) , both sides of Eq. ( 13) are multiplied by x and integrated over space to get a simple expression for ∂ t ( x(t) ) from which finally x(t) is obtained. One proceeds similarly to get x(t) 2 , but here both sides of Eq. ( 13) are multiplied by x 2 . The resulting expression is:

A(t n ) = 1 - sin(η/2) η/2 2 (v 2 0 ∆t)t n = D ef f (η) t n . (18)
A similar expression can be obtained by proposing a continuum-time process for the orientation dynamics, where ∆t is replaced by the inverse of the stochastic turning rate, see [START_REF] Peruani | [END_REF]. Notice . The solid black line refers to the prediction for S F given by Eq. ( 19), while the dashed line is the approximation for S LC given by Eq. ( 20). (c) shows the trajectory of the center of mass of the particle ensemble corresponding to the numerical experiment with η = 0.3 (solid line) and η = 1.5 (dot-dashed line).

that in the limit of η → 0, D ef f (η → 0) vanishes, while V (η → 0) = v 0 . Thus, in this limit, the motion becomes purely ballistic. On the other hand, the diffusive limit corresponds to η → 2π, where V (η → 2π) = 0 and D ef f (η → 2π) = v 2 0 ∆t. This corresponds to the diffusion coefficient of an ensemble of regular random walkers that make at each ∆t a step with equal probability in any direction. Fig. 1 shows a comparison between Eq.( 18) and simulations with N = 1000 self-propelled particles interacting by nematic alignment at various noise intensities η. As shown in Fig. 1, the agreement between simulations and the predictions of Eq. ( 13) is fairly good. Now, we turn our attention to the orientational order parameters exhibited by this cloud of moving particles. At first glance, it might seem that, since the bunch of particles moves coherently, S F and S LC have to be both equal to one. Figs. 2(a) and (b) show S F and S LC as function of time in simulation with N = 1000 self-propelled particles with nematic alignment at two different noise strength values. Fig. 2 indicates that S F and S LC inside the coherently moving swarm are still functions of η. We compute S F in terms of P (θ) by inserting Eq.( 10) into the average given by Eq. [START_REF] Zhang | Proc. Natl. Acad. Sci[END_REF], and obtain:

S F = 2 η sin(η/2) . (19)
As observed by Dossetti et al., this approximation corresponds in fact to the limit of N → ∞ -see Eq. (A8) and its derivation in [START_REF] Dossetti | [END_REF]. Similarly for S LC , by inserting Eq.( 10) into Eq.( 8) we obtain:

S LC = sin(η) η 2 . ( 20
)
Fig. 3 compares the predictions of Eqs. ( 19) and ( 20) and simulations performed with N = 1000 particles. The symbols correspond to temporal averages of simulated time series of S F and S LC as shown in Figs. 2(a) and (b). The agreement is remarkably good for the 10 4 integration steps that the simulations span. It is worth to notice that during this period the center of mass of the swarm travels a distance d much larger than the interaction radius ǫ, i.e., d ≫ ǫ (see Fig. 2(c)), while the swarm's expansion around its center of mass is comparatively very small. Despite the fact that the swarm moves in a quite coherent way, the value of the order parameters drops with increasing η (see Figs. 2(a) and (b)), and the center of mass performs a more pronounced meandering trajectory as shown in Fig. 2(c). Interestingly, the fluctuations of the center of mass around the y axis are much larger than the swarm spreading around its center. This is due to large fluctuations of the average total momentum vector that are not reflected in the fluctuations of S F . The total momentum vector fluctuates in direction and modulus, but S F fluctuations are only related to fluctuations of its modulus. The same applies to the orientation tensor Q and to its associated scalar liquid crystal order parameter S LC . 19) (solid curve) and (20) (dashed curve), respectively. Symbols correspond to S F (triangles) and S LC (squares) measured in simulations performed with N = 1000 self-propelled particles with nematic alignment initially located at the origin and all pointing in direction + x.

As said above, this simple description of the swarm is valid until the spreading around the center of mass is such that the density of the moving cluster falls below percolation. The critical density is obtained by assuming overlapping discs of radius r = ǫ/2 = 1/2, ρ p ∼ 1.45 [32]. Now, let us imagine that the swarm does not evolve in an infinite space, but in a box with periodic boundary conditions. Let us assume in addition that the size L 2 of the box is such that N/L 2 is much larger than ρ p . Can this system be described in terms of such simple equations as Eq. ( 19) and ( 20) for all times? If the initial condition of the system is completely random, does the system reach the same steady state? For a random initial condition, we cannot expect a steady state with a unique direction of motion. We assume that two main opposite direction of motion emerge, α 0 and α 0 + π, and that the angular dynamics of the i-th particle is simply given by:

θ t+∆t i = α 0 + η t i with probability p + (α 0 + π) + η t i with probability p - (21)
In consequence, P (θ) = (p + /η)g(θ, α 0 , η) + (p -/η)g(θ, α 0 + π, η). Using this expression to compute the order parameters S F and S LC , we find that S F = (2/η) sin(η/2)(p + -p -), while S LC is given by Eq. ( 20). This hypothesis is tested in Fig. 4. The simulations were performed with N = 2 12 particles at high density (ρ = 4) in a box with periodic boundary conditions. The initial condition was random and the simulations were carried out for 10 6 time steps. The solid curve corresponds to the approximation given by Eq. ( 20), where no fitting parameter is used. The dot-dashed curve is a fitting of the first 8 data points to the left of η c = 2 assuming S LC ∼ (η c -η) β , where β = 0.46±0.03. This assumption may turn to be wrong for large system sizes. The figure shows that Eq. ( 20) provides a good approximation of S LC for small values of η, but fails to describe orientational order at large values of the noise intensity, i.e., close to the transition point. In summary, Fig. 4 shows that at least for small system sizes the ordering dynamics can be described through the above approach for low noise intensity, i.e., far away from the the critical point. The obtained results cannot be used to extrapolate the behavior of the system to very large system sizes, though well in the ordered phase they may still hold. This would require a systematic analysis of finite-size effects in large-scale simulations, which is beyond the scope of this work. Moreover, we known from [21] that in large systems of particles obeying Eqs.( 5) and ( 6), fluctuations play a crucial role in the emerging macroscopic dynamics.

Cluster size distribution at low density

In this section we derive a simple theory to understand the emergence of steady state cluster size distributions in self-propelled particle systems. We illustrate the phenomenon of cluster formation in these systems by performing simulations of self-propelled particles with nematic alignment. Fig. 5 shows simulation snapshots for three values of the angular noise η. Notice that as η is decreased, clusters become significantly larger.

Formally, the cluster dynamics of these systems can be described by deriving a master equation for the evolution of the probability p(n(t)), where n(t) = n 1 (t), n 2 (t), ..., n N (t), with n 1 (t) being the number of isolated particles, n 2 (t) the number of two-particle clusters, n 3 (t) the number of three-particle clusters, etc. This kind of approach has also been used to understand equilibrium nucleation in gases, where the transition probabilities between states are function of the associated free energy change [START_REF] Schimansky-Geier | Self-organization by Nonlinear Irreversible Processes[END_REF][START_REF] Schweitzer | [END_REF]35]. Here, however, we will adopt an alternative more phenomenological strategy for our non-equilibrium problem. Instead of looking for the complete description of the clustering process in terms of p(n(t)), we will derive equations for the time evolution of the mean value of n 1 (t), n 2 (t), ..., n N (t). To ease the notation, we will refer to n 1 (t) , n 2 (t) , ..., n N (t) simply as n 1 (t), n 2 (t), etc.

The simulations were performed with N = 2 14 particles at low density, ρ = 0.25. We observe a transition from an apparently homogeneous density at large noise values to pronounced formation of clusters at low noise (see Fig. 5). Fig. 6 shows the (weighted) cluster size distribution p(m) = (m n m (t))/N for various values of the noise intensity η. We have considered that a cluster is an ensemble of connected particles. Two particles are connected whenever they are separated a distance less or equal to d -in Fig. 6, d = 2ǫ. The figure indicates that there is transition in the cluster size distribution from a monotonically decreasing distribution for large values of η to a distribution with a peak at large cluster sizes for small enough noise values. Interestingly, at the transition point the cluster size distribution follows a power-law.

We look for an explanation for the observed clustering phenomena. Through the study of the coherence of an initially perfectly oriented swarm, we have learned that clusters have a finite life-time before they break into parts. Now we incorporate the fact that a moving cluster collects particles whose relative direction of motion is such that |∆θ| < π/2, provided the noise intensity η is low enough. If the topology of the system is a finite torus, these two effects, spreading of particles due to fluctuations in the direction of motion and collection of particles due to random collisions of clusters, reach an equilibrium and the cluster size distribution (CSD) becomes a steady distribution. As mentioned above, we look for a description of the process in terms of the number n i (t) of clusters with i particles at time t. The time evolution equations for the n i (t) have the following form:

ṅ1 = 2B 2 n 2 + N k=3 B k n k - N -1 k=1 A k,1 n k n 1 ṅj = B j+1 n j+1 -B j n j - N -j k=1 A k,j n k n j + 1 2 j-1 k=1 A k,j-k n k n j-k for j = 2, ....., N -1 ṅN = -B N n N + 1 2 N -1 k=1 A k,N -k n k n N -k (22)
where the dot denotes the time derivative, B j represents the rate for a cluster of mass j to loose a particle, and is defined as 22) with N = 128 for various values of the dimensionless parameter P , where P = 2ǫd 2 v 0 aD ef f (η) . Notice the transition from a monotonically decreasing distribution for small values of P to a non-monotonic distribution with a peak at large cluster sizes for large values of P . and A j,k is the collision rate between clusters of mass j and k, defined by

B j = D ef f (η) d 2 j , (23)
A j,k = v 0 2ǫ a j + √ k , (24)
where a = L 2 is the area of the two-dimensional space where particles move. In Eq. ( 23), d denotes, as before, the maximum distance that two particles can be separated apart to be considered as connected. In this way, d 2 /D ef f (η) is the characteristic time a particle on a cluster boundary needs to detach from a cluster. The splitting rate B j is proportional to the inverse of this characteristic time multiplied by the number of particles on the boundary, which we estimate as √ j. On the other hand, the collision rate A j,k is derived in analogy to the collision rate in kinetic gas theory [START_REF] Reif | Fundamentals of Statistical and Thermal Physics[END_REF] between two disk-like particles A and B which is known to be proportional to the relative velocity of the particles and the sum of their diameters. We have approximated the diameter of a cluster of mass j by ǫ √ j. Direct numerical integration of these equations shows that Eq.( 22) produces qualitatively similar distributions as the one observed in individual-based simulations, see Fig. 7. The different curves correspond to various values of the dimensionless parameter P , defined as P = 2ǫd 2 v0 aD ef f (η) . For small values of P , which correspond to large values of η, the distribution p(m) monotonically decreases with m, while for large values of P , resp. small values of η, a peak at large cluster sizes emerges. A quantitative comparison between Eq.( 22) and individual-based simulations is beyond the scope of this work.

Conclusions

Clusters play a fundamental role in the macroscopic dynamics of self-propelled particle systems. Particularly, clusters are capable to transport orientation information for long distances. Orientational order and cluster dynamics are often closely linked.

Here we have addressed the problem of cluster dynamics in self-propelled particle systems in a systematic way. We have first studied single isolated clusters and learned that their dynamics can be understood in terms of simple mean-field arguments. We have derived expression for the diffusion coefficients of clusters and found that the spatial spreading of clusters is anisotropic, with the long axis of the cluster perpendicular to its moving direction. We have also shown that the mean-field approaches used to characterize the cluster dynamics can be applied to describe the behavior of self-propelled particles for small system sizes far away from the transition point.

Finally, we have analyzed the emergence of steady state cluster size distributions in systems of self-propelled particles that obey Eqs. ( 1) and [START_REF] Cavagna | Proc. Natl. Acad. Sci[END_REF]. Particularly, we have derived a set of equations of the Smoluchowski type to describe the cluster dynamics in the system. From direct numerical integration of these equations, we have shown that this approach leads to qualitatively similar distributions.

This study presents a first step towards the understanding of the complex clustering behavior exhibited by self-propelled particle systems. Several key aspects have been left out in this first approach. For instance, a system size analysis is required to check the validity of these results in large systems. Correlations among clusters have been completely ignored, though they may be crucial to the macroscopic dynamics. Finally, we suspect that clustering effects may induce giant number fluctuations as the ones predicted in [36]. However, a link between these two phenomena has not been still stablished. All these issues shall be the subject of future studies.

Fig. 1 .

 1 Fig.1. Spreading A of the swarm around it center of mass vs. time for different values of the noise intensity η. Symbols correspond to measurements of A(t) in simulations with N = 1000 self-propelled particles with nematic alignment, which initially were located at the origin and all pointing in direction + x. The dashed lines correspond to the approximation given by Eq.(18)

Fig. 2 .

 2 Fig. 2. (a) and (b): S LC and S F vs. time for a simulation performed with N = 1000 self-propelled particles with nematic alignment, initially located at the origin and all pointing in direction + x. The noise intensity corresponds to η = 0.3 in (a) and η = 1.5 in (b). The solid black line refers to the prediction for S F given by Eq. (19), while the dashed line is the approximation for S LC given by Eq. (20). (c) shows the trajectory of the center of mass of the particle ensemble corresponding to the numerical experiment with η = 0.3 (solid line) and η = 1.5 (dot-dashed line).

Fig. 3 .

 3 Fig.3. S F and S LC as function of η as predicted by Eq. (19) (solid curve) and (20) (dashed curve), respectively. Symbols correspond to S F (triangles) and S LC (squares) measured in simulations performed with N = 1000 self-propelled particles with nematic alignment initially located at the origin and all pointing in direction + x.

Fig. 4 .

 4 Fig. 4. S LC n vs. η. Symbols correspond to simulations with N = 2 12 self-propelled particles with nematic alignment at density ρ = 4 in a box with periodic boundary conditions and random initial conditions. The dot-dashed curve is a fitting of the first 8 data points to the left of ηc = 2 through S LC ∼ (ηcη) β, where β = 0.46 ± 0.03 (this may not be true for large system sizes). The solid curve corresponds to Eq. (20).

Fig. 5 .

 5 Fig. 5. Snapshots of simulations with self-propelled particles with nematic alignment at density ρ = 0.25 after 2.5 10 5 times units. The values of η are: 1.4 in (a), 1.3 in (b), and 0.6 in (c). The spatial arrangement of particles seen in the figures is representative of what is observed at the steady state.

Fig. 6 .

 6 Fig.6. Clustering at different values of the noise intensity η at low density. Simulation performed with N = 2 14 self-propelled particles with nematic alignment at density ρ = 0.25. Notice the transition from a monotonically decreasing distribution for large values of η to a non-monotonic distribution with a peak at large cluster sizes for small enough values of noise intensity η.

Fig. 7 .

 7 Fig. 7. Steady state cluster size distributions obtained from numerical integration of Eqs.(22) with N = 128 for various values of the dimensionless parameter P , where P = 2ǫd 2 v 0aD ef f(η) . Notice the transition from a monotonically decreasing distribution for small values of P to a non-monotonic distribution with a peak at large cluster sizes for large values of P .