Collective motion of self-propelled particles interacting without cohesion

We present a comprehensive study of Vicsek-style self-propelled particle models in two and three space dimensions. The onset of collective motion in such stochastic models with only local alignment interactions is studied in detail and shown to be discontinuous (first-order like). The properties of the ordered, collectively moving phase are investigated. In a large domain of parameter space including the transition region, well-defined high-density and high-order propagating solitary structures are shown to dominate the dynamics. Far enough from the transition region, on the other hand, these objects are not present. A statistically-homogeneous ordered phase is then observed, which is characterized by anomalously-strong density fluctuations, superdiffusion, and strong intermittency.


I. INTRODUCTION
Collective motion phenomena in nature have attracted the interest of scientists and other authors for quite a long time [1]. The question of the advantage of living and moving in groups, for instance, is a favorite one among evolutionary biologists [2]. In a different perspective, physicists are mostly concerned with the mechanisms at the origin of collective motion, especially when it manifests itself as a true, non-trivial, emerging phenomenon, i.e. in the absence of some obvious cause like the existence of a leader followed by the group, a strong geometrical constraint forcing the displacement, or some external field or gradient felt by the whole population. Moreover, the ubiquity of the phenomenon at all scales, from intra-cellular molecular cooperative motion to the displacement in group of large animals, raises, for physicists at least, the question of the existence of some universal features possibly shared among many different situations.
One way of approaching these problems is to construct and study minimal models of collective motion: if universal properties of collective motion do exist, then they should appear clearly within such models and thus could be efficiently determined there, before being tested for in more elaborate models and real-world experiments or observations. Such is the underlying motivation of recent studies of collective motion by a string of physicists [3,4,5,6,7,8]. Among them, the group of Tamas Vicsek has put forward what is probably the simplest possible model exhibiting collective motion in a non-trivial manner.
In the "Vicsek model" [9], point particles move off-lattice at constant speed v 0 , adjusting their direction of motion to that of the average velocity of their neighbors, up to some noise term accounting for external or internal perturbations (see below for a precise definition). For a finite density of particles in a finite box, perfect alignment is reached easily in the absence of noise: in this fluctuation-less collective motion, the macroscopic velocity equals the microscopic one. On the other hand, for strong noise particles are essentially non-interacting random walkers and their macroscopic velocity is zero, up to statistical fluctuations.
Vicsek et al. showed that the onset of collective motion occurs at a finite noise level. In other words, there exists, in the asymptotic limit, a fluctuating phase where the macroscopic velocity of the total population is, on average, finite. Working mostly in two space dimensions, they concluded, on the basis of numerical [5,9] simulations, that the onset of this ordered motion is well described as a novel non-equilibrium continuous phase transition leading to long range order, at odds with equilibrium where the continuous XY symmetry cannot be spontaneously broken in two space dimensions and below [10]. This brought support to the idea of universal properties of collective motion since the scaling exponents and functions associated to such phase transitions are expected to bear some degree of universality, even out of equilibrium.
Here, after a definition of the models involved (Section II), we come back, in Section III, to this central issue and present a rather comprehensive study of the onset of collective motion in Vicsek-style models. In Section IV, we describe the ordered, collective motion phase. Section V is devoted to a general discussion of our results together with some perspectives. Most of the numerical results shown were obtained in two space dimensions, but we also present three-dimensional results. Wherever no explicit mention is made, the default space dimension is two. Similarly, the default boundary conditions are periodic in a square or cubic domain.

A. Vicsek model: angular noise
Let us first recall the dynamical rule defining the Vicsek model [9]. Point particles labeled by an integer index i move off-lattice in a space of dimension d with a velocity v i of fixed modulus v 0 = | v i |. The direction of motion of particle i depends on the average velocity of all particles (including i) in the spherical neighbourhood S i of radius r 0 centered on i. The discrete-time dynamics is synchronous: the direction of motion and the position of all particles are updated at each timestep ∆t, in a driven, overdamped manner: where ϑ is a normalization operator (ϑ( w) = w/| w|) and R η performs a random rotation uniformely distributed around the argument vector: in d = 2, R η v is uniformely distributed around v inside an arc of amplitude 2π η; in d = 3, it lies in the solid angle subtended by a spherical cap of amplitude 4π η and centered around v. The particles positions r i are then simply updated by streaming along the chosen direction as in Note that the original updating scheme proposed by Vicsek et al. in [9] defined the speed as a backward difference, although we are using a forward difference. The simpler updating above, now adopted in most studies of Vicsek-style models, is not expected to yield different results in the asymptotic limit of infinite size and time.

B. A different noise term: vectorial noise
The "angular" noise term in the model defined above can be thought of as arising from the errors committed when particles try to follow the locally-averaged direction of motion. One could argue on the other hand that most of the randomness stems from the evaluation of each interaction between particle i and one of its neighbors, because, e.g., of perception errors or turbulent fluctuations in the medium. This suggests to replace Eq.(1) by: where ξ is a random unit vector and N i is the number of particles in S i . It is easy to realize that this "vectorial" noise acts differently on the system. While the intensity of angular noise is independent from the degree of local alignment, the influence of the vectorial noise decreases with increasing local order.

C. Repulsive force
In the original formulation of the Vicsek model as well as in the two variants defined above, the only interaction is alignment. In a separate work [22], we introduced a two-body repulsion/attraction force, to account for the possibility of maintaining the cohesion of a flock in an infinite space (something the Vicsek model does not allow). Here, we only study models without cohesion. Nevertheless, we have considered, in the following, the case of a pairwise repulsion "force", to estimate in particular the possible influence of the absence of volume exclusion effects in the basic model, which leaves the local density actually unbounded. We thus introduce the short ranged, purely repulsive interaction exerted by particle j on particle i: where e ij is the unit vector pointing from particle i to j and r c < r 0 is the typical repulsion range. Equations (1) and (3) are then respectively generalized to and where β measures the relative strength of repulsion with respect to alignment and noise strength.

D. Control and order parameters
The natural order parameter for our polar particles is simply the macroscopic mean velocity, conveniently nor- malized by the microscopic velocity v 0 where · i stands for the average over the whole population. Here, we mostly consider its modulus ϕ(t) = | ϕ(t)|, the scalar order parameter.
In the following, we set, without loss of generality, ∆t = 1 and r 0 = 1, and express all time-and lengthscales in terms of these units. Moreover repulsive force will be studied by fixing r c = 0.127 and β = 2.5.
This leaves us with two main parameters for these models: the noise amplitude η and the global density of particles ρ. Recently, the microscopic velocity v 0 has been argued to play a major role as well [39]. All three parameters (η, ρ, and v 0 ) are considered below.

III. ORDER-DISORDER TRANSITION AT THE ONSET OF COLLECTIVE MOTION
As mentioned above, the original Vicsek model attracted a lot of attention mostly because of the conclu- sions drawn from the early numerical studies [5,9]: the onset of collective motion was found to be a novel continuous phase transition spontaneously breaking rotational symmetry. However, it was later shown in [38] that beyond the typical sizes considered originally the discontinuous nature of the transition emerges, irrespective of the form of the noise term. Recently, the discontinuous character of the transition was argued to disappear in the limit of small v 0 [39]. We now address the problem of the nature of the transition in full detail.
Even though there is no rigorous theory for finitesize-scaling (FSS) for out-of-equilibrium phase transitions, there exists now ample evidence that one can safely rely on the knowledge gained in equilibrium systems [41,42,43]. The FSS approach [44,45] involves the numerical estimation of various moments of the order parameter distribution as the linear system size L is systematically varied. Of particular interest are the variance t and the so-called Binder cumulant where · t indicates time average. The Binder cumulant is especially useful in the case of continuous phase transitions, because it is one of the simplest ratio of moments which takes a universal value at the critical point η t , where all the curves G(η, L), obtained at different system sizes L, cross each other. At a first-order transition point, on the other hand, the Binder cumulant exhibits a sharp drop towards negative values [46]. This minimum is due to the simultaneous contributions of the two phases coexisting at threshold. Moreover, it is easy to compute that G(η, L) ≈ 2/3 in the ordered phase, while for a disordered state with a continuous rotational simmetry one has G(η, L) ≈ 1/3 in d = 2 and G(η, L) ≈ 4/9 in d = 3.

A. Overture
As an overture, we analyze systems of moderate size in two dimensions (N ≈ 10 4 particles) at the density ρ = 2, typical of the initial studies by Vicsek et al., but with the slightly modified update rule (2) and for both angular and vectorial noise. The microscopic velocity is set to v 0 = 0.5.
For angular noise, the transition looks indeed continuous, as found by Vicsek et al.. On the other hand, the time averaged scalar order parameter ϕ t displays a sharp drop for vectorial noise, and the Binder cumulant exibits a minimum at the transition point, indicating a discontinuous phase transition ( Fig. 1a-b). Simultaneously, the variance is almost delta-peaked. The difference betwen the two cases is also recorded in the probability distribution function (PDF) of ϕ which is bimodal (phase coexistence) in the vectorial noise case (Fig. 1c-d).
The qualitative difference observed upon changing the way noise is implemented in the dynamics is, however, only a finite-size effect. As shown in [38], the transition in the angular noise case reveals its asymptotic discontinuous character provided large-enough system sizes L are considered ( Fig. 2a-b). Remaining for now at a qualitative level, we show in Fig. 2c a typical time-series of the order parameter for the angular noise case in a large system in the transition region. The sudden jumps from the disordered phase to the ordered one and vice-versa are evidence for metastability and phase coexistence. Note that the system size beyond which the transition reveals its discontinuous character for the angular noise case at density ρ = 2 and velocity v 0 = 0.5 -the conditions of the original papers by Vicsek et al.-is of the order of L = 128, the maximum size then considered. It is clear also from Fig. 1 that the discontinuous nature of the transition appears earlier, when increasing system size, for vectorial noise than for angular noise. Thus, finite-size effects are stronger for angular noise. The same is true when one is in presence of the repulsive interactions (Fig.3). Finally, the same scenario holds in three space dimensions, with a discontinous phase transition separating the ordered from the disordered phases for both angular and vectorial noise (Fig. 4).
Before proceeding to a study of the complete phase diagram, we detail now how a comprehensive FSS study can be performed on a particular case.

B. Complete FSS analysis
For historical reasons, the following study has been performed on the model with vectorial noise and repulsive force (Eq. 6). It has not been repeated in the simpler case of the "pure" Vicsek model because its already high numerical cost would have been prohibitive due to the strong finite size effects.
As a first step, we estimated the correlation time τ (L), whose knowledge is needed to control the quality of timeaveraging: the duration T of numerical simulations has been taken much larger than τ (L) (T = 100τ in the largest systems, but typically 10000τ for smaller sizes). Moreover, τ is also useful to correctly estimate the statistical errors on the various moments (as ϕ t , σ, and G ) of the PDF of the order parameter, for which we used the Jackknife procedure [47]. The correlation time was estimated near the transition (where it is expected to be largest) as function of system size L measuring the exponential decay rate of the correlation function ( Fig. 5a) We found τ to vary roughly linearly with L (see Fig. 5b).
It is interesting to observe that at equilibrium, one would expect τ to scale as [48] where κ is the surface tension of the metastable state. Therefore, our result implies a very small or vanishing surface tension κ ≪ 1/L, a situation reminiscent of observations made in the cohesive case [22], where the surface tension of a cohesive droplet was found to vanish near the onset.
Following Borgs and Koteckỳ [49], the asymptotic coexistence point η t (i.e. the first order transition point) can be determined from the asymptotic convergence of various moments of the order parameter PDF. First, the observed discontinuity in ϕ(t) t , located at η ϕ (L), is expected to converge exponentially to η t with L. Second, the location of the susceptibility peak η σ (L) -which is the same as the peak in σ provided some fluctuationdissipation relation holds, see the Appendix -also converges to η t , albeit algebraically with an exponent γ σ . Third, the location of the minimum of G, η G (L), is also expected to converge algebraically to η t with with an ex- Interestingly, the value taken by these exponents actually depends on the number of phases and of the dimension d of the system: for two-phase coexistence one has γ G = γ σ = 2d, while for more than two phases γ G = γ σ = d. In Figure 6, we show that our data are in good agreement with all these predictions. The three estimates of η t are consistent with each other within numerical accuracy. Moreover η ϕ (L) is found to converge exponentially to the transitional noise amplitude, while both η σ (L) and η G (L) show algebraic convergence with an exponent close to 2. This agrees with the fact that due to the continuous rotational symmetry, the ordered phase is degenerate and amounts to an infinite number of possible phases.

C. Hysteresis
One of the classical hallmarks of discontinuous phase transitions is the presence, near the transition, of the hysteresis phenomenon: ramping the control parameter at a fixed (slow) rate up and down through the transition point, a hysteresis loop is formed, inside which phase coexistence is manifest (see Fig. 7a for the d = 3 case with vectorial noise). The size of such hysteresis loops varies with the ramping rate. An intrinsic way of assessing phase coexistence and hysteresis is to study systematically the nucleation time τ ↑ needed to jump from the disordered phase to the ordered one, as well as τ ↓ the decay time after which the ordered phase falls into the disordered one. Fig. 7b shows, in three space dimensions, how these nucleation and decay vary with η at two different sizes. A sharp divergence is observed, corresponding to the transition point. At a given time value τ , one can read, from the distance between the "up" and the "down" curve, the average size of hysteresis loops for ramping rates of the order of 1/τ .

D. Phase diagram
The above detailed FSS study would be very tedious to realize when varying systematically the main parameters η, ρ, and v 0 , as well as the nature of the noise and the presence or not of repulsive interactions. From now on, to characterize the discontinuous nature of the transition, we rely mainly on the presence, at large-enough system sizes L, of a minimum in the variation of the Binder cumulant G with η (all other parameters being fixed). We call L * the crossover size marking the emergence of a minimum of G(η). We are now in the position to sketch the phase diagram in the (η, ρ, v 0 ) parameter space. The numerical protocol used is, at given parameter values, to run a large-enough system so that the discontinuous character of the transition is seen (i.e. L > L * ). For larger sizes, the location of the transition point typically varies very little, so that for most practical purposes, locating the (asymptotic) transition point from systems sizes around L * is satisfactory.
The results presented below are in agreement with simple mean-field-like arguments in the diluted limit: in the small-ρ regime, one typically expects that the lower the density, the lower the transitional noise amplitude η t . Indeed, for ∆t v 0 of the order of or not much smaller than the interaction range r 0 and in the low-density limit ρ ≪ 1/r 0 d , the system can be seen as a dilute gas in which particles interact by short range ordering forces only. In this regime, the persistence length of an isolated particle (i.e. the distance travelled before its velocity loses correlation with its initial direction of motion) varies like v 0 /η. To allow for an ordered state, the noise amplitude should be small enough so that the persistence length remains larger than the average inter-particle distance, i.e. 1/ρ 1/d . Thus the transition noise amplitude is expected to behave as In [5], it was indeed found that η t ∼ ρ α with α ≃ 1 2 in two dimensions. Our own data ( Fig. 8a-c) now confirm Eq. (10) for both the angular and vectorial noise in two and three spatial dimensions, down to very small ρ values. The data deviate from the square-root behavior as the average inter-particle distance gets of the order of or smaller than the interaction range. Finally, we also investigated the transition line when v 0 is varied (Fig. 8d). For the vectorial noise case, at fixed density, the threshold noise value η t is almost constant (data obtained at ρ = 1 2 , not shown). For the angular noise, in the small v 0 limit where the above mean-field argument does not apply, we confirm the first-order character of the phase transition down to v 0 ≈ 0.05 for both angular and vectorial noise (Fig. 9b,d). For even smaller values of v 0 , the investigation becomes numerically too costly (see section below). Note that η t seems to be finite when v 0 → 0 + , a limit corresponding to the XY model on a randomly connected graph. Still for angular noise, the large velocity limit is also difficult to study numerically. Again, we observe that the transition is discontinuous as far as we can probe it, i.e. v 0 = 20 (Fig. 9a,b).

E. Special limits and strength of finite-size effects
We now discuss particular limits of the models above together with the relative importance of finite-size effects. Recall that these are quantified by the estimated value of the crossover size L * beyond which the transition appears discontinuous. All the following results have been obtained for d = 2. Partial results in three dimensions indicate that the same conclusions should hold there. Keep in mind that in all cases reported, the transition is discontinuous. We are just interested here in how large a system one should use in order to reach the asymptotic regime. Fig. 10a shows that finite-size effects are stronger for angular noise than for vectorial noise for all densities ρ at which we are able to perform these measurements. Note in particular that at ρ = 2, the density originally used by Vicsek et al., L * ∼ 128 for angular noise, while it is very small for vectorial noise, confirming the observation made in Section III A.
In the small-ρ limit, the discontinuous character of the transition appears later and later, with L * roughly diverging as 1/ √ ρ (inset of Fig. 10a). Note that this means that in the small-ρ limit one needs approximately the same number of particles to start observing the discontinuity.
The large-ρ limit reveals a difference between angular and vectorial noise: while L * remains small for vectorial noise, it seems to diverge for angular noise (Fig. 10a), making this case difficult to study numerically.
We also explored the role the microscopic velocity v 0 in the strength of finite-size effects. Qualitatively, the effects observed are similar to those just reported when the density is varied (Fig. 10b). In the small v 0 limit, we record a strong increase of L * as v 0 → 0 for both types of noise. In the large velocity limit, L * decreases for vectorial noise, whereas it increases for angular noise.

F. Summary and discussion
The summary of the above lengthy study of the orderdisorder transition in Vicsek-like models is simple: for any finite density ρ, any finite velocity v 0 , and for both types of noise introduced, the transition is discontinuous. This was observed even in the numerically-difficult limits of large or small ρ or v 0 . These results contradict recent claims made about the angular noise case (original Vicsek model). We now comment on these claims.
Vicsek and co-workers [39], showed that keeping the density and the noise intensity fixed, a qualitative change is observed when decreasing v 0 : for not too small v 0 values, in the ordered phase, particles diffuse anisotropically (and the transition is discontinuous), while diffusion becomes isotropic at small v 0 , something interpreted as a sign of a continuous transition in this region. Rather than the convoluted arguments presented there, what happens is in fact rather simple: decreasing v 0 at fixed ρ and η, one can in fact cross the transition line, passing from the ordered phase (where particles obviously diffuse anisotropically due to the transverse superdiffusive effects discussed in Section IV C) to the disordered phase. Our Figure 8d, obtained in the same conditions as in [39] (apart from harmless change of time updating rule), shows that if one keeps η = 0.1 (as in [39]), one crosses the transition line at about v 0 ≃ 0.1, the value invoked by Vicsek and co-workers to mark a crossover from discontinuous to "continuous" transitions.
In a recent letter [40], Aldana et al. study orderdisorder phase transitions in random network models and show that the nature of these transitions may change with the way noise is implemented in the dynamics (they consider the angular and vectorial noises defined here). Arguing that these networks are limiting cases of Vicseklike models, they claim that the conclusions reached for the networks carry over to the transition to collective motion of the VM-like systems. They conclude in particular that in the case of "angular" noise the transition to collective motion is continuous. We agree with the analysis of the network models, but the claim that they are relevant as limits of Vicsek-like models is just wrong: the data presented there (Figure 1 of [40]) to substantiate this claim is contradicted by our Figures 9a,c (see also [50]) obtained at larger system sizes. Again, for large-enough system sizes, the transition is indeed discontinuous. Thus, at best, the network models of Aldana et al. constitute a singular v 0 → ∞ limit of Vicsek-like models.

IV. NATURE OF THE ORDERED PHASE
We now turn our attention to the ordered, symmetrybroken phase. In previous analytical studies, it has often been assumed that the density in the ordered phase is spatially homogeneous, albeit with possibly large fluctuations (see, e.g. [8]). This is indeed what has been reported in early numerical studies, in particular by Vicsek et al. [9]. In the following, we show that this is not true in large enough systems, where, for a wide range of noise amplitudes near the transition point, density fluctuations lead to the formation of localized, travelling, high-density and high-order structures. At low enough noise strength, though, a spatially-homogeneous ordered phase is found, albeit with unusually strong density fluctuations. transversally with respect to the mean direction of motion, and have a center of mass velocity close to v 0 . While particles inside bands are ordered and, in the asymptotic regime, move coherently with the global mean velocity, particles lying outside bands -in low density regionsare not ordered and perform random walks.
As shown in Fig. 11a-c for angular noise dynamics (1), there exists a typical system size L b , below which the bands or sheets cannot be observed. Numerical simulations indicate that L b depends only weakly on the noise amplitude and is of the same order of magnitude as the crossover size marking the appearance of the discontinuous character of the transition: L b ≈ L * . It is therefore numerically easier to observe bands in the ordered phase of vectorial noise dynamics (3), as in Fig 11d. Bands may be observed asymptotically without and with a repulsive interaction (Fig 12c) and for both kinds of noise. They appear for various choices of boundary conditions (see for instance Fig. 12a-b, where reflecting boundary conditions have been employed), which may play a role in determining the symmetry-broken mean direction of motion. For instance, bands travelling parallel to one of the axis are favoured when periodic boundary conditions are employed in a rectangular box (they represent the simplest way in which an extended structure can wrap around a torus, and are thus reached more easily from disordered initial conditions), but bands travelling in other directions may also appear, albeit with a smaller probability.
Bands can be described quantitatively through local quantities, such as the local density ρ ℓ ( x, t), measured inside a domain V( x) centered around x, and the local order parameter Further averaging these local quantities perpendicularly to the mean velocity (7) one has the density profile ρ ⊥ (x , t) = ρ ℓ ( x, t) ⊥ and the order parameter profile ϕ ⊥ (x , t) = ϕ ℓ ( x, t) ⊥ , where x indicates the longitudinal direction w.r.t. mean velocity. Bands are characterized by a sharp kink in both the density and the order parameter profiles (see Fig. 13a and 13c-d). They are typically asymmetric, as it can be expected for moving structures, with a rather sharp front edge, a well defined mid-height width w -which typically is of the same order as L b -and an exponentially decaying tail with a characteristic decay length of the order of the w (Fig. 13b). Large systems may accomodate several bands at the same time, typically all moving in the same direction (see for instance Fig. 11c and the density profile in Fig. 14e). However, they do not form well-defined wave trains, but rather a collection of solitary objects, as hinted by the following numerical experiments. We investigated the instability of the densityhomogeneous, ordered state in a series of numerical simulations starting from particles uniformly distributed in space but strictly oriented along the major axis in a large rectangular domain. Figure 14a,b show space-time plots of the density profile: initially flat, it develops structures with no well-defined wavelength (Fig. 14c). Density fluctuations destroy the initially-ordered state in a rather unusual way: a dynamical Fourier analysis of the density profile show a weakly-peaked, wide band of wavelengths growing subexponentially (Fig. 14d). This is at odds with a finite-wavelength supercritical instability, which would lead to a wavetrain of traveling bands. Furthermore, the asymptotic (late time) power spectra of the density profiles are not peaked around a single frequency either, but rather broadly distributed over a large range of wavenumbers (Fig. 14f). In the asymptotic regime, bands are extremely long-lived metastable (or possibly stable) objects, which are never equally-spaced (a typical late-time configuration is shown in Fig. 14e).
To summarize, the emerging band or sheet structure in the asymptotic regime is not a regular wave train characterized by a single wavelength, but rather a collection of irregularly-spaced localized traveling objects, probably weakly interacting through their exponentially decaying tails.

B. Low-noise regime and giant density fluctuations
As the noise amplitude is decreased away from the transition point, bands are less sharp, and eventually disappear, giving way to an ordered state characterized by an homogeneous local order parameter and large fluctuations of the local density.
A quantitative measure of the presence, in the ordered phase, of structures spanning the dimension transverse to the mean motion (i.e. bands or sheets) is provided by the variances of the density and order parameter profiles: where · indicates the average of the profile in the longitudinal direction with respect to mean velocity. Indeed, these profile widths vanish in the infinite-size limit except if band/sheet structures are present. In Fig 15a,b, we plot these profile widths averaged over time as a function of noise amplitude. Both quantities present a maximum close to the transition point in the ordered phase, and drop drastically as soon as the disordered phase is entered. Lowering the noise away from the transition point, these profiles decrease steadily: bands/sheets stand less sharply out of the disordered background (Fig 15c). At some point (η ≈ 0.3 for the parameters values considered in Fig 15a,b), bands rather abruptly disappear and are no more well-defined transversal objects. It is difficult to define this point accurately, but it is clear that for lower noise intensities the local order parameter is strongly homogeneous in space. Nevertheless, fluctuations in the density field are strong (Fig 15d), but can no more give rise to (meta)stable longlived transverse structures. Density fluctuations in the bandless regime are in fact anomalously strong: measuring number fluctuations in sub-systems of linear size ℓ, we find that their root meansquare ∆n does not scale like the square root of n = ρℓ d , the mean number of particles they contain; rather we find σ(n) ∝ n α with α ≈ 0.8 both in d = 2 and in d = 3 (Fig. 16a). This is reminiscent of the recent discovery of "giant density fluctuations" in active nematics [14,51,52]. However, the theoretical argument which initially predicted such fluctuations [53] cannot be invoked directly in the present case. (Indeed, the above value of α, although needing to be refined, does not seem to be compatible with the prediction α = 1 2 + 1 d made in [53].) More work is needed to fully understand under which circumstances the coupling between density and order in systems of "active", self-propelled particles gives rise to such anomalous density fluctuations.

C. Transverse superdiffusion
According to Toner and Tu predictions [8,33,34], the dynamics of the symmetry-broken ordered phase of polar active particles should be characterized by a superdiffusive mean square displacement in the direction(s) transversal to mean velocity. In particular, in d = 2 one has [34] with ν = 4/3. While this analytical result has been succesfully tested by numerical simulations of models with cohesive interactions [34,38], numerical simulations in models without cohesion present substantial difficulties, mainly due to the presence of continuously merging and splitting sub-clusters of particles moving coherently (as discussed in Section IV D). As a consequence, an ensemble of test particles in a cohesionless model is exposed to different "transport" regimes (w.r.t. center of mass motion) which are not well separated in time. When the mean displacement is averaged at fixed time, this tends to mask the transverse superdiffusion.
To overcome this problem, we chose to follow [54] and to measure τ 2 , the average time taken by two particles to double their transverse separation distance δ ⊥ . From Eq. (14) one immediately has with 2/ν = 3/2 in d = 2. In order to easily separate the transverse from the parallel component, we considered an ordered system in a large rectangular domain with periodic boundary conditions and the mean velocity initially oriented along the long side. The mean direction of motion then stays oriented along this major axis, so that we can identify the transverse direction with the minor axis. Furthermore, a high density and a small (angular) noise amplitude (corresponding to the bandless regime) have been chosen to avoid the appearance of large, locally disordered patches. Our results (Fig. 16b) confirm the prediction of Toner and Tu: transverse superdiffusion holds at low-enough noise, while normal diffusion is observed in the disordered, high-noise phase. Note that the systematic deviation appearing in our data at some large scale is induced by large fluctuations in the orientation of the global mean velocity during our numerical simulations (not shown here).
We take the opportunity of this discussion to come back to the superdiffusive behavior of particles observed in the transition region [21]. There, sub-clusters emerge propagate ballistically and isotropically due to the absence of a well-established global order. Particles trajectories consist in "ballistic flights", occurring when a particle is caught in one of these coherently moving clusters, alternated with ordinary diffusion in disordered regions. The mean square displacement of particles exhibits the scaling ∆x 2 = | x(t) − x(0))| 2 i ∝ t 5/3 [21]. In view of our current understanding of the discontinuous nature of the transition, we now tend to believe that this isotropic superdiffusion is probably not asymptotic.

D. Internal structure of the ordered region
We now turn our attention to the internal structure of the ordered regimes. As we noted in the previous section, these regimes do not consist of a single cluster of interacting particles moving coherently. Even in the case where high density bands/sheets are present, these are in fact dynamical objects made of splitting and merging clusters. Note that for the models considered here, clusters are unambiguously defined thanks to the strictly finite interaction range r 0 .
As noticed first by Aldana and Huepe [7], clusters of size n are distributed algebraically in the ordered region, i.e. P (n) ∼ n −µ . But a closer look reveals that the exponent µ characterizing the distribution of cluster sizes changes with the distance to the transition point. For noise intensities not too far from the threshold, when bands are observed, we find µ values larger than 2, whereas µ < 2 in the bandless regimes present at low noise intensities (Fig. 17a,b).
Thus, bands are truly complex, non-trivial structures emerging out of the transverse dynamics of clusters with a well defined mean size (since µ > 2). It is only in the bandless regime that one can speak, as Aldana and Huepe, of "strong intermittency". We note in passing that the parameter values they considered correspond in fact to a case where bands are easily observed (at larger sizes than those considered in [7]). Thus, clusters do have a well-defined mean size in their case. Consequently, the probability distribution P (ϕ) of the order parameter ϕ does not show the behavior reported in Fig. 1   as soon as the system size is large enough. Whether in the band/sheet regime or not, P (ϕ) shows essentially Gaussian tails, is strongly peaked around its mean, and its variance decreases with system size (Fig. 17c,d).
Although the picture of intermittent bursts between "laminar"intervals proposed by Aldana and Huepe has thus to be abandoned, the anomalous density fluctuations reported in the previous section are probably tantamount to the strong intermittency of cluster dynamics in the bandless regime. Again, these phenomena, reported also in the context of "active" nematics [14,51,53], deserve further investigation.

E. Phase ordering
The ordered regimes presented above are the result of some transient evolution. In particular, the bands/sheets are the typical asymptotic structures appearing in finite domains with appropriate boundary conditions. In an infinite system, the phase ordering process is, on the other hand, infinite, and worth studying for its own sake.
Numerically, we have chosen to start from highly disordered initial conditions which have a homogeneous density and vanishing local order parameter. In practice, we quench a system "thermalized" at strong noise to a smaller, subcritical, η value. Typical snapshots show the emergence of structures whose typical scale seems to increase fast (Fig. 18). During this domain growth, we monitor the two-point spatial correlation function of both the density and velocity fields. These fields are defined by a coarse-graining over a small lengthscale ℓ (typically 4). These correlation functions have an unusual shape (Fig. 19a): after some rather fast initial decay, they display an algebraic behavior whose effective exponent decreases with time, and finally display a near-exponential cut-off. As a result, they cannot be easily collapsed on a single curve using a simple, unique, rescaling length-scale. Nevertheless, using the late exponential cut-off, a correlation length ξ can be extracted. Such a lengthscale ξ grows roughly linearly with time (Fig. 19b). Qualitatively similar results are obtained whether the noise strength is in the range where bands/sheets appear in finite boxes or not.
We note that the above growth law is reminiscent of that of the so-called model H of the classification of Halperin and Hohenberg [55]. Since this model describes, in principle, the phase separation in a viscous binary fluid, the fast growth observed could thus be linked to the hydrodynamic modes expected in any continuous description of Vicsek-like models [8,56].

A. Summary of main results
We now summarize our main results before discussing them at a somewhat more general level.
We have provided ample evidence that the onset of collective motion in Vicsek-style models is a discontinuous (first-order) phase transition, with all expected hallmarks, in agreement with [38]. We have made the (numerical) effort of showing this in the limits of small and large velocity and/or density.
We have shown that the ordered phase is divided in two regions: near the transition and down to rather low noise intensities, solitary structures spanning the directions transverse to the global collective velocity (the bands or sheets) appear, leading to an inhomogeneous density field. For weaker noise, on the other hand, no such structures appear, but strong, anomalous density fluctuations exist and particles undergoes superdiffusive motion transverse to the mean velocity direction.
Finally, we have reported a linear growth (with time) of ordered domains when a disordered configuration is quenched in the ordered phase. This fast growth can probably be linked to the expected emergence of long wavelength hydrodynamic modes in the ordered phase of active polar particles models.

B. On the role of bands/sheets
The high-density high-order traveling bands or sheets described here appear central to our main findings. They seem to be intimately linked to the discontinuous character of the transition which can, to some extent, be considered as the stability limit of these objects. In the range of noise values where they are observed, the anomalous density fluctuations present at lower noise intensities are suppressed.
One may then wonder about the universality of these objects. Simple variants of the Vicsek-style models studied here (e.g. restricting the interactions to binary ones involving only the nearest neighbor) do exhibit bands and sheets [11]. Moreover, the continuous deterministic description derived by Bertin et al. [11] does possess localized, propagating solitary solutions rather similar to bands [57]. Although the stability of these solutions need to be further investigated, these results indicate that these objects are robust and that their existence is guaranteed beyond microscopic "details". However, the emergence of regular, stable, bands and sheets is obviously conditioned to the shape and the boundary conditions of the domain in which the particles are allowed to move. In rectangular domains with at least one periodic direction, these objects can form, span across the whole domain, and move. But in, say, a circular domain with reflecting boundary conditions, they cannot develop freely, being repeatedly "frustrated". Nevertheless, simulations performed in such a geometry indicate that the transition remains discontinuous, with the ordered phase consisting of one or several dense packets traveling along the circular boundary. Note, though, that these packets intermittently emit elongated structures (bands) traveling towards the interior of the disk before colliding on the boundary. To sum up, bands appear as the "natural" objects in the transition region, but they may be prevented by the boundaries to develop into full-size straight objects.
At any rate, time-series of the order parameter such as the one presented in Fig. 2 clearly show that the transition is discontinuous irrespective of the geometry and boundaries of the domain, and thus of whether bands/sheets can develop into stable regular structures or not: the sudden, abrupt, jumps from the disordered state to some ordered structure are tantamount to a nucleation phenomenon characteristic of a discontinuous transition.

C. A speculative picture
We would now like to offer the following speculative general picture: the key feature of the Vicsek-like models studied here -as well as of other models for active media made of self-propelled particles [11,35,56,58,59]is the coupling between density and order. Particles are forced to move, and, since they carry informations about the order, advection, density fluctuations and order are intimately linked. High-density means strong local order (if noise is low enough) because the many particles in a given neighborhood will adopt roughly the same orientation. The reverse is also true: in a highly-ordered region, particles will remain together for a long time and thus will sweep many other particles leading to a denser and denser group.
At a given noise level, one can thus relate, in the spirit of some local equilibrium hypothesis, local density to local order. In practice, such an "equation of state" approach can be justified by looking, e.g., at a scatter plot of local order parameter vs local density. Fig. 20 reveals that, in the ordered bandless regime, such a scatter plot is characterized by a plateau over a large range of local density values corresponding to order, followed, below some crossover density, by more disordered local patches. The regions in space where local density is below this crossover level do not percolate in the bandless regime, and order can be maintained very steadily in the whole domain (this is corroborated by the fact that in spite of the large, anomalous density fluctuations, the order parameter field is, on the other hand, rather constant, see Fig. 15c). The noise intensity at which bands emerge roughly corresponds to the value where the low-density disordered regions percolate. The reamining disconnected, dense patches then eventually self-organize into bands/sheets. The emergence of these elongated structures is rather natural: moving packets elongate spontaneously because they collect many particles; superdiffusion in the directions perpendicular to the mean motion endow these nascent bands/sheets with some "rigidity". At still stronger noise, the bands/sheets are destroyed and global order disappears.
The above features are at the root of the approach by Toner and Tu [8,33,34]. Their predictions of strong density fluctuations, transverse superdiffusion, and peculiar sound propagation properties are correct as long as bands/sheets do not exist, i.e. for not too strong noise intensities. This is indeed in agreement with their assumption that the density field is statistically homogeneous in the frame moving at the global velocity (albeit with strong fluctuations), which is only true in the bandless regime.

D. Outlook
The results presented here are almost entirely numerical. Although they were obtained with care, they need to be ultimately backed up by more analytical results. A first step is the derivation of a continuous description in terms of a density and a velocity field (or some combination of the two), which would allow to go beyond microscopic "details". In that respect the deterministic equation derived by Bertin et al. from a Boltzmann description in the dilute limit [11] is encouraging. However, one may suspect that intrinsic fluctuations are crucial in the systems considered here if only because some of the effective noise terms will be multiplicative in the density. A mesoscopic, stochastic equation description is thus a priori preferable. This is especially true in view of the "giant" anomalous fluctuations present in the bandless ordered phase. These fluctuations clearly deserve further investigation, all the more so as they seem to be generic features of active particle models [14,51].
Ongoing work is devoted to both these general issues.
instance, Fig. 6 of [5]). A mean-field analysis has confirmed this view, showing a logarithmic variation of the response with the field intensity [61].
To bypass this problem, we have preferred to use the following equation: where the vectorial noise was chosen because, as shown above, it leads more easily to the asymptotic regime. The effective intensity of the field is now proportional to the local order.
Two (scalar) response functions can be defined in our problem: the longitudinal response R along the field direction, and the transverse response R ⊥ . We consider the former. In practice, we quenched, at time t 0 , a strongnoise, highly-disordered system ϕ(t 0 ) ≈ 0 to a smaller noise value and started applying the constant, homogeneous field h immediately. We then followed the response of the system by monitoring the growth of the order parameter. We measured the susceptibility χ which is nothing but the integrated response function: χ (t, t 0 ) = t t0 R (t, t ′ )dt ′ . In practice, we have: In a well-behaved system, the susceptibility should be independent of the amplitude of the field, at least at small enough values ("linear" regime). This is what we observed, as shown in Fig. 21a.
Correspondingly, the correlation function is defined as: The fluctuation-dissipation relation (A1) can then be written in its integrated form: In Fig. 21b, we show that χ and C are related linearly in time, confirming the validity of this relation and allowing an estimation of T eff .
This well-defined -although not uniquely definedeffective temperature varies as expected in parameter space. In particular, it increases systematically with the noise strength η (Fig. 21c), although this variation is not linear. Note also, that, intriguingly, there is a small jump of T eff at the noise value corresponding to the transition in this case (η t ≃ 0.55).