# Collective motion of self-propelled particles interacting without cohesion

Hugues Chaté, Francisco Ginelli, Guillaume Grégoire, Franck Raynaud

► To cite this version:

Hugues Chaté, Francisco Ginelli, Guillaume Grégoire, Franck Raynaud. Collective motion of self-propelled particles interacting without cohesion. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2008, 77 (4), pp.046113. 10.1103/PhysRevE.77.046113. hal-00517100

HAL Id: hal-00517100

https://hal.science/hal-00517100v1

Submitted on 31 Jul 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Collective motion of self-propelled particles interacting without cohesion

Hugues Chaté and Francesco Ginelli

CEA—Service de Physique de l'État Condensé, Centre d'Etudes de Saclay, 91191 Gif-sur-Yvette, France

Guillaume Grégoire

Matière et Systèmes Complexes, CNRS UMR 7057, Université Paris-Diderot, Paris, France

Franck Raynaud

CEA—Service de Physique de l'État Condensé, Centre d'Etudes de Saclay, 91191 Gif-sur-Yvette, France
and Matière et Systèmes Complexes, CNRS UMR 7057, Université Paris-Diderot, Paris, France

(Received 12 December 2007; published 18 April 2008)

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.

DOI: 10.1103/PhysRevE.77.046113

PACS number(s): 05.65.+b, 64.70.qj, 87.18.Gh

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, nontrivial, 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 intracellular molecular cooperative motion to the displacement in groups 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–8]. Among them, the group of Vicsek has put forward what is probably the simplest possible model exhibiting collective motion in a nontrivial 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 fluctuationless collective motion, the macroscopic velocity equals the microscopic one. On the other hand, for strong noise particles are essentially noninteracting 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 nonequilibrium 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 with such phase transitions are expected to bear some degree of universality, even out of equilibrium.

The above results caused a well-deserved stir and prompted a large number of studies at various levels [3,4,6–8,11–37]. In particular, two of us showed that the onset of collective motion is in fact discontinuous [38], and that the original conclusion of Vicsek et al. was based on numerical results obtained at too small sizes [5,9]. More recently, the discontinuous character of the transition was challenged in two publications, one by Vicsek and co-workers [39] and another by Aldana et al. [40].

Here, after a definition of the models involved (Sec. II), we come back, in Sec. III, to this central issue and present a rather comprehensive study of the onset of collective motion in Vicsek-style models. In Sec. 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 2. Similarly, the default boundary conditions are periodic in a square or cubic domain.

## II. THE MODELS

### 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 $\vec{v}_i$ of fixed modulus $v_0 = |\vec{v}_i|$. The direction of motion of particle *i* depends on the average velocity of all particles (including *i*) in the spherical neighborhood $S_i$ of radius $r_0$ centered on *i*. The discrete-time dynamics is synchronous: the direction of motion and the positions of all particles are updated at each time step $\Delta t$, in a driven, overdamped manner,

$$ \vec{v}_i(t + \Delta t) = v_0(\mathcal{R}_{\eta} \circ \vartheta)\left(\sum_{j \in S_i} \vec{v}_j(t)\right), \quad (1) $$

where $\vartheta$ is a normalization operator $[\vartheta(\vec{w}) = \vec{w}/|\vec{w}|]$ and $\mathcal{R}_{\eta}$ performs a random rotation uniformly distributed around the argument vector: in $d=2$, $\mathcal{R}_{\eta}\vec{v}$ is uniformly distributed around $\vec{v}$ inside an arc of amplitude $2\pi\eta$; in $d=3$, it lies in the solid angle subtended by a spherical cap of amplitude $4\pi\eta$ and centered around $\vec{v}$. The particle positions $\vec{r}_i$ are then simply updated by streaming along the chosen direction as in

$$ \vec{r}_i(t + \Delta t) = \vec{r}_i(t) + \Delta t \vec{v}_i(t + \Delta t). \quad (2) $$

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 the replacement of Eq. (1) by

$$ \vec{v}_i(t + \Delta t) = v_0 \vartheta \left( \sum_{j \in S_i} \vec{v}_j(t) + \eta \mathcal{N}_i \vec{\xi} \right), \quad (3) $$

where $\vec{\xi}$ is a random unit vector and $\mathcal{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 of 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 or 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*:

$$ \vec{f}_{ij} = -\vec{e}_{ij} \times [1 + \exp(|\vec{r}_j - \vec{r}_i|/r_c - 2)]^{-1}, \quad (4) $$

where $\vec{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

$$ \vec{v}_i(t + \Delta t) = v_0(\mathcal{R}_{\eta} \circ \vartheta)\left(\sum_{j \in S_i} \vec{v}_j(t) + \beta \sum_{j \in S_i} \vec{f}_{ij}\right) \quad (5) $$

and

$$ \vec{v}_i(t + \Delta t) = v_0 \vartheta \left( \sum_{j \in S_i} \vec{v}_j(t) + \beta \sum_{j \in S_i} \vec{f}_{ij} + \eta \mathcal{N}_i \vec{\xi} \right), \quad (6) $$

where $\beta$ 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 normalized by the microscopic velocity $v_0$,

$$ \vec{\varphi}(t) = \frac{1}{v_0} \langle \vec{v}_i(t) \rangle_i, \quad (7) $$

where $\langle \cdot \rangle_i$ stands for the average over the whole population. Here, we mostly consider its modulus $\varphi(t) = |\vec{\varphi}(t)|$, the scalar order parameter.

In the following, we set, without loss of generality, $\Delta t = 1$ and $r_0 = 1$, and express all time and length scales in terms of these units. Moreover, the repulsive force will be studied by fixing $r_c = 0.127$ and $\beta = 2.5$.

This leaves us with two main parameters for these models: the noise amplitude $\eta$ and the global density of particles $\rho$. Recently, the microscopic velocity $v_0$ has been argued to play a major role as well [39]. All three parameters ($\eta$, $\rho$, 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 conclusions 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,
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Typical behavior across the onset of collective motion for moderate-sized models ($\rho=2$, $v_0=0.5$, $L=64$) with angular (black circles) and vectorial noise (red triangles). The reduced noise amplitude $\epsilon=1-\eta/\eta_t$ is shown on the abscissa (transition points estimated at $\eta_t=0.6144(2)$—vectorial noise—and $\eta_t=0.478(5)$—angular noise). (a) Time-averaged order parameter $\langle\varphi(t)\rangle_t$. (b) Binder cumulant $G$. (c) Variance $\sigma$ of $\varphi$. (d) Order parameter distribution function $P$ at the transition point. Bimodal distribution for vectorial noise dynamics (red dashed line); unimodal shape for angular noise (black solid line). Time averages have been computed over $3 \times 10^5$ time steps.</figcaption></figure>

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 finite-size 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–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

$$\sigma(\eta, L) = L^d(\langle\varphi^2\rangle_t - \langle\varphi\rangle_t^2)$$

and the so-called Binder cumulant

$$G(\eta, L) = 1 - \frac{\langle\varphi^4\rangle_t}{3\langle\varphi^2\rangle_t^2}, \qquad (8)$$

where $\langle\cdot\rangle_t$ indicates time average. The Binder cumulant is especially useful in the case of continuous phase transitions, because it is one of the simplest ratios of moments which takes a universal value at the critical point $\eta_t$, where all the curves $G(\eta, 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 toward negative values [46]. This minimum is due to the simultaneous contributions of the two phases existing at threshold. Moreover, it is easy to compute that $G(\eta, L) \approx 2/3$ in the ordered phase, while for a disordered state with a continuous rota-

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) FSS analysis of angular noise dynamics ($\rho=2$, $v_0=0.5$, time averages computed over $2 \times 10^7$ time steps). Time-averaged order parameter (a) and Binder cumulant (b) as a function of noise for various system sizes $L$. (c) Piece of an order parameter time series close to the transition point ($L=256$, $\eta=0.476$).</figcaption></figure>

tional symmetry one has $G(\eta, L) \approx 1/3$ in $d=2$ and $G(\eta, L) \approx 4/9$ in $d=3$.

## A. Overture

As an overture, we analyze systems of moderate size in two dimensions ($N \approx 10^4$ particles) at the density $\rho=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 $\langle\varphi\rangle_t$ displays a sharp drop for vectorial noise, and the Binder cumulant exhibits a minimum at the transition point, indicating a discontinuous phase transition [Figs. 1(a) and 1(b)]. Simultaneously, the variance is almost $\delta$ peaked. The difference between the two cases is also recorded in the probability distribution function (PDF) of $\varphi$ which is bimodal (phase coexistence) in the vectorial noise case [Figs. 1(c) and 1(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 [Figs. 2(a) and 2(b)]. Remaining for now at a qualitative level, we show in Fig. 2(c) 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
<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color online) Transition to collective motion with short-range repulsive interactions. Left panels: angular noise. Right panels: vectorial noise. (a),(b) Order parameter vs noise amplitude at different system sizes. (c),(d) Binder cumulant G as a function of noise amplitude. ($\rho$=2, $v_0$=0.3, time averages carried over $10^7$ time steps.)</figcaption></figure>

density $\rho$=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 the system size is increased, for vectorial noise than for angular noise. Thus, finite-size effects are stronger for angular noise. The same is true when one is in the presence of repulsive interactions (Fig. 3). Finally, the same scenario holds in three space dimensions, with a discontinuous 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 $\tau(L)$, whose knowledge is needed to control the quality of time averaging: the duration $T$ of numerical simulations has been taken much larger than $\tau(L)$ ($T$=100$\tau$ in the largest systems, but typically 10 000$\tau$ for smaller sizes). Moreover, $\tau$ is also useful to correctly estimate the statistical errors on the various moments (as $\langle\phi\rangle_t$, $\sigma$, 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$

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Transition to collective motion in three spatial dimensions. Left panels: angular noise. Right panels: vectorial noise. (a),(b) Time-averaged order parameter vs noise amplitude at different system sizes. (c),(d) Binder cumulant G as a function of noise amplitude. ($\rho$=0.5, $v_0$=0.5, time averages carried over $10^5$ time steps.)</figcaption></figure>

measuring the exponential decay rate of the correlation function [Fig. 5(a)]

$$C(t) = \langle\phi(t_0)\phi(t_0+t)\rangle_t - \langle\phi(t_0)\rangle_{t_0}^2 \sim \exp\left(-\frac{t}{\tau}\right). \quad (9)$$

We found $\tau$ to vary roughly linearly with $L$ [see Fig. 5(b)]. It is interesting to observe that, at equilibrium, one would expect $\tau$ to scale as [48]

$$\tau = L^{d/2} \exp(\kappa L^{d-1}),$$

where $\kappa$ is the surface tension of the metastable state. Therefore, our result implies a very small or vanishing surface tension $\kappa \ll 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.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Correlation time $\tau$ of the order parameter near the transition point for vectorial noise dynamics with repulsion. System parameters are $\rho$=2, $v_0$=0.5, and $\eta \approx \eta_t$. (a) Time correlation function $C(t)$ at $L$=128. The lin-log inset shows the exponential decay. (b) Correlation time $\tau$ as a function of system size $L$. The dashed line marks linear growth with $L$. Correlation functions were computed on samples of $\approx 10^6$ realizations for typically $10^3$ time steps.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

FIG. 6. (Color online) FSS analysis of vectorial dynamics with short-range repulsive force ($\rho=2$, $v_0=0.3$). Convergence of the finite-size transition points measured from different moments of the order parameter FSS to the asymptotic transition point $\eta_t$ (see Fig. 3). Upper panels: finite-size transition points estimated from (a) time average, (b) variance, and (c) Binder cumulant. The horizontal dashed line marks the estimated asymptotic threshold $\eta_t = 0.5544(1)$. Lower panels: scaling of the finite-size reduced noise $\epsilon=1-\eta/\eta_t$ transition point. (d) Exponential convergence for the jump location in the time-averaged order parameter. (e) Power-law behavior of the variance peak position. (f) Power-law behavior of the Binder cumulant minimum. The dashed lines in (e) and (f) mark the estimated exponents $\gamma_\sigma=\gamma_G=2$.

Following Borgs and Kotecky [49], the asymptotic coexistence point $\eta_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 $\langle\varphi(t)\rangle_t$, located at $\eta_\varphi(L)$, is expected to converge exponentially to $\eta_t$ with $L$. Second, the location of the susceptibility peak $\eta_\sigma(L)$—which is the same as the peak in $\sigma$ provided some fluctuation-dissipation relation holds (see the Appendix)—also converges to $\eta_t$, albeit algebraically with an exponent $\gamma_\sigma$. Third, the location of the minimum of $G$, $\eta_G(L)$, is also expected to converge algebraically to $\eta_t$ with an exponent $\gamma_G=\gamma_\sigma$.

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 $\gamma_G=\gamma_\sigma=2d$, while for more than two phases $\gamma_G=\gamma_\sigma=d$. In Fig. 6, we show that our data are in good agreement with all these predictions. The three estimates of $\eta_t$ are consistent with each other within numerical accuracy. Moreover, $\eta_\varphi(L)$ is found to converge exponentially to the transitional noise amplitude, while both $\eta_\sigma(L)$ and $\eta_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: on 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. 7(a) 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 $\tau_\uparrow$ needed to jump from the disordered phase to the ordered one, as well as $\tau_\downarrow$, the decay time after which the ordered phase falls into the disordered one. Figure 7(b) shows, in three space dimensions, how these nucleation and decay rates vary with $\eta$ at two different sizes. A sharp divergence is observed, corresponding to the transition point. At a given time value $\tau$, one can read, from the distance between the “up” and the “down” curves, the average size of hysteresis loops for ramping rates of the order of $1/\tau$.

<figure><img src="image_8.png" /></figure>

FIG. 7. (Color online) Hysteresis in three spatial dimensions with vectorial noise. (a) Order parameter vs noise strength along the hysteresis loop observed with a ramp rate of $2 \times 10^{-6}$ per time step ($\rho=1/2$, $v_0=0.5$, $L=32$). Empty circles mark the path along the adiabatic increase of noise amplitude; full triangles for adiabatic decrease. (b) Nucleation times from the disordered phase to the ordered phase ($\tau_\uparrow$, left curves) and vice versa ($\tau_\downarrow$, right curves) for two system sizes [other parameters as in (a)]. Each point is averaged over 1000 realizations.

### D. Phase diagram

The above detailed FSS study would be very tedious to realize when the main parameters $\eta$, $\rho$, and $v_0$ are varied systematically, 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 $\eta$ (all other parameters being fixed). We call $L^*$ the crossover size marking the emergence of a minimum of $G(\eta)$.

We are now in the position to sketch the phase diagram in the $(\eta, \rho, 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-$\rho$ regime, one typically expects that the lower the density, the
<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) Asymptotic phase diagrams for the transition to collective motion. (a) Two space dimensions: threshold amplitude $\eta_t$ for angular noise as a function of density $\rho$ at $v_0 = 0.5$. Inset: Log-log plot to compare the low-density behavior with the mean-field predicted behavior $\eta_t \sim \sqrt{\rho}$ (dashed red line). (b) As in (a), but with vectorial noise dynamics. (c) Noise-density phase diagram in three dimensions for vectorial noise dynamics at fixed velocity $v_0=0.5$. In the log-log inset the transition line can be compared with the predicted behavior $\eta_t \sim \rho^{1/3}$ (dashed red line). (d) Two space dimensions: threshold amplitude $\eta_t$ for angular noise as a function of particle velocity $v_0$ at fixed density $\rho=1/2$ (black circles) and 1/8 (red triangles). The horizontal dashed line marks the noise amplitude considered in Ref. [39] (see Sec. III F).</figcaption></figure>

lower the transitional noise amplitude $\eta_t$. Indeed, for $\Delta v_0$ of the order of or not much smaller than the interaction range $r_0$ and in the low-density limit $\rho \ll 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 traveled before its velocity loses correlation with its initial direction of motion) varies as $v_0/\eta$. To allow for an ordered state, the noise amplitude should be small enough so that the persistence length remains larger than the average interparticle distance, i.e., $1/\rho^{1/d}$. Thus the transition noise amplitude is expected to behave as

$$ \eta_t \sim v_0 \rho^{1/d}. \qquad (10) $$

In [5], it was indeed found that $\eta_t \sim \rho^\alpha$ with $\alpha \approx 1/2$ in two dimensions. Our own data [Figs. 8(a)-8(c)] now confirm Eq. (10) for both the angular and vectorial noise in two and three spatial dimensions, down to very small $\rho$ values. The data deviate from the square-root behavior as the average interparticle distance becomes of the order of or smaller than the interaction range.

Finally, we also investigated the transition line when $v_0$ is varied [Fig. 8(d)]. For the vectorial noise case, at fixed density, the threshold noise value $\eta_t$ is almost constant (data obtained at $\rho=1/2$, not shown). For the angular noise, in the small-$v_0$ limit where the above mean-field argument does not

<figure><img src="image_10.png" /><figcaption>FIG. 9. (Color online) First-order transition for angular noise dynamics at high (left panels) and low (right panels) velocity $v_0$. Typical averaging time is $\approx 10^6$ time steps. (a) Time-averaged order parameter and (c) Binder cumulant at large particle velocity for angular noise in two spatial dimensions at increasing velocities and $L \approx L^*(\rho=2)$. (b) Time-averaged order parameter and (d) Binder cumulant for $v_0=0.05$ and two increasing system sizes ($\rho=1/2$).</figcaption></figure>

apply, we confirm the first-order character of the phase transition down to $v_0 \approx 0.05$ for both angular and vectorial noise [Figs. 9(b) and 9(d)]. For even smaller values of $v_0$, the investigation becomes numerically too costly (see the next subsection). Note that $\eta_t$ seems to be finite when $v_0 \to 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$ [Figs. 9(a) and 9(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.

Figure 10(a) shows that finite-size effects are stronger for angular noise than for vectorial noise for all densities $\rho$ at which we are able to perform these measurements. Note in particular that, at $\rho=2$, the density originally used by Vicsek et al., $L^* \sim 128$ for angular noise, while it is very small for vectorial noise, confirming the observation made in Sec. III A.

In the small-$\rho$ limit, the discontinuous character of the transition appears later and later, with $L^*$ roughly diverging
<figure><img src="image_11.png" /></figure>

FIG. 10. (Color online) Crossover system size $L^*$ above which the discontinuous character of the transition appears [as testified by the existence of a minimum in the $G(\eta)$ curve]. Black circles: angular noise. Red triangles: vectorial noise. (a) $L^*$ vs $\rho$ for $v_0=0.5$. Inset: the low-density behavior in log-log scales; the dashed line marks a power-law divergence proportional to $1/\sqrt{\rho}$. (b) $L^*$ vs $v_0$ at fixed density (open symbols, $\rho=1/2$; filled symbols, $\rho=2.0$).

as $1/\sqrt{\rho}$ [inset of Fig. 10(a)]. Note that this means that in the small-$\rho$ limit one needs approximately the same number of particles to start observing the discontinuity.

The large-$\rho$ limit reveals a difference between angular and vectorial noise: while $L^*$ remains small for vectorial noise, it seems to diverge for angular noise [Fig. 10(a)], making this case difficult to study numerically.

We also explored the role of 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. 10(b)]. In the small-$v_0$ limit, we record a strong increase of $L^*$ as $v_0 \to 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 order-disorder transition in Vicsek-like models is simple: for any finite density $\rho$, any finite velocity $v_0$, and both types of noise introduced, the transition is discontinuous. This was observed even in the numerically difficult limits of large or small $\rho$ 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, when the density and the noise intensity are kept fixed, a qualitative change is observed when $v_0$ is decreased: 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: by decreasing $v_0$ at fixed $\rho$ and $\eta$, 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 Sec. IV C) to the disordered phase. Our Fig. 8(d), obtained in the same conditions as in [39] (apart from harmless change of the time-updating rule), shows that if one keeps $\eta=0.1$ (as in [39]), one crosses the transition line at about $v_0 \approx 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 order-disorder 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 Vicsek-like models, they claim that the conclusions reached for the networks carry over to the transition to collective motion of the Vicsek-model-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 (Fig. 1 of [40]) to substantiate this claim are contradicted by our Figs. 9(a) and 9(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 \to \infty$ limit of Vicsek-like models.

## IV. NATURE OF THE ORDERED PHASE

We now turn our attention to the ordered, symmetry-broken 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, traveling, 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.

### A. Traveling in bands

Numerical simulations of the ordered phase dynamics ($\eta<\eta_t$), performed at large enough noise amplitudes, are characterized by the emergence of high-density moving bands ($d=2$) or sheets ($d=3$). Typical examples are given in Figs. 11 and 12. These moving structures appear for large enough systems after some transient. They extend 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 regions—are not ordered and perform random walks.

As shown in Fig. 11(a)–11(c) for the 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 \approx L^*$. It is therefore numerically
<figure><img src="image_12.png" /><figcaption>FIG. 11. (Color online) Typical snapshots in the ordered phase. Points represent the position of individual particles and the red arrow points along the global direction of motion. (a)–(c) Angular noise, $\rho=1/2$, $v_0=0.5$, $\eta=0.3$, and increasing system sizes, respectively, $L=64$, 256, and 1024. Sharp bands can be observed only if $L$ is larger than the typical bandwidth $w$. (d) Vectorial noise, $\rho=1/2$, $v_0=0.5$, $\eta=0.55$, and $L=64$: bands appear at relatively small system sizes for this type of noise. For clarity, only a representative sample of 10 000 particles is shown in (b) and (c). Boundary conditions are periodic.</figcaption></figure>

easier to observe bands in the ordered phase of vectorial noise dynamics (3), as in Fig. 11(d).

Bands may be observed asymptotically without and with a repulsive interaction [Fig. 12(c)] and for both kinds of noise. They appear for various choices of boundary conditions [see, for instance, Figs. 12(a) and 12(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 traveling parallel to one of the axes are favored 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 traveling in other directions may also appear, albeit with a smaller probability.

Bands can be described quantitatively through local quantities, such as the local density $\rho_{\ell}(\vec{x}, t)$, measured inside a domain $V(\vec{x})$ centered around $\vec{x}$, and the local order parameter

$$ \varphi_{\ell}(\vec{x}, t) = \frac{1}{v_0} |\langle \vec{v}_i(t) \rangle_{\vec{r}_i \in V(\vec{x})}|. \quad (11) $$

Further averaging these local quantities perpendicularly to the mean velocity (7), one has the density profile $\rho_{\perp}(x_{\parallel}, t) = \langle \rho_{\ell}(\vec{x}, t) \rangle_{\perp}$ and the order parameter profile $\varphi_{\perp}(x_{\parallel}, t) = \langle \varphi_{\ell}(\vec{x}, t) \rangle_{\perp}$, where $x_{\parallel}$ indicates the longitudinal direction with respect to the mean velocity. Bands are characterized by a sharp kink in both the density and the order parameter profiles [see Figs. 13(a), 13(c), and 13(d)]. They are typically

<figure><img src="image_13.png" /><figcaption>FIG. 12. (Color online) Same as Fig. 11 but in different geometries and boundary conditions or space dimensions. (a),(b) Vectorial noise ($\eta=0.325$, $\rho=1/8$, and $v_0=0.5$); boundary conditions are periodic along the $y$ (vertical) axis and reflecting in $x$. (a) A long single band travels along the periodic direction. (b) The domain size along the periodic direction is too small to accommodate bands, and a single band bouncing back and forth along the nonperiodic direction is observed. (c) Angular noise, repulsive force, and periodic boundary conditions ($\rho=2$, $\eta=0.23$, and $v_0=0.3$). (d) High-density sheet traveling in a three-dimensional box with periodic boundary conditions (angular noise with amplitude $\eta=0.355$, $\rho=1/2$, and $v_0=0.5$).</figcaption></figure>

asymmetric, as 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 $w$ [Fig. 13(b)].

Large systems may accommodate several bands at the same time, typically all moving in the same direction [see, for instance, Fig. 11(c) and the density profile in Fig. 14(e)]. 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 density-homogeneous, 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. Figures 14(a) and 14(b) show space-time plots of the density profile: initially flat, it develops structures with no well-defined wavelength [Fig. 14(c)]. 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. 14(d)]. This is at odds with a finite-wavelength supercritical instability, which would lead to a wave train of traveling bands. Furthermore, the asymptotic
<figure><img src="image_14.png" /><figcaption>FIG. 13. (Color online) (a) Typical density (black line) and order parameter (dashed red line) profiles for bands in two dimensions (vectorial noise, $\rho=2$, $\eta=0.6$, and $v_0=0.5$). (b) Tail of the density profile shown in (d) (black line) and its fit (blue dashed line) by the formula $\rho_{\perp}(x_{\parallel}, t) \approx a_0 + a_1(t) \exp(-x_{\parallel}/w)$, with $w \approx 6.3$ (lin-log scale). (c), (d) Traveling sheet in three dimensions (angular noise, $\rho=1/2$, $\eta=0.355$, and $v_0=0.5$). (e) Projection of particle positions on a plane containing the global direction of motion (marked by red arrow). (f) Density (black line) and order parameter (dashed red line) profiles along the direction of motion $x_{\parallel}$.</figcaption></figure>

(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 wave numbers [Fig. 14(f)]. 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. 14(e)].

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 a 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 variations of the density and order parameter profiles:

$$ \begin{aligned} \Delta\rho_{\perp}^2(t) &= \langle [\rho_{\perp}(x_{\parallel}, t) - \langle\rho_{\perp}(x_{\parallel}, t)\rangle_{\parallel}]^2 \rangle_{\parallel}, \\ \Delta\varphi_{\perp}^2(t) &= \langle [\varphi_{\perp}(x_{\parallel}, t) - \langle\varphi_{\perp}(x_{\parallel}, t)\rangle_{\parallel}]^2 \rangle_{\parallel}, \end{aligned} \quad (12) $$

where $\langle \cdot \rangle_{\parallel}$ indicates the average of the profile in the longitudinal direction with respect to mean velocity. Indeed, these

<figure><img src="image_15.png" /><figcaption>FIG. 14. (Color online) Emergence of high-density high-order traveling bands ($d=2$) from a spatially homogeneous (uniformly distributed random positions) initial condition with all particle velocities oriented along the major axis of a 196 × 1960 domain with periodic boundary conditions. Vectorial noise of amplitude $\eta=0.6$, density $\rho=2$, and $v_0=0.5$. (a) Space-time plot of the density profile. Time is running from left to right from $t=0$ to 12 000, while the longitudinal direction is represented on the ordinates. Color scale from blue (low values) to red (high values). (b) Same as (a) but at later times (from $t=148$ 000 to 160 000). (c) Spatial Fourier power spectrum $S$ of an early density profile ($t=12$ 000). (d) Early-time evolution of selected Fourier modes $k=10, 23, 28, 33, 41, 76$, and 121 (the black lowest curve is for $k=10$, the other curves are not significantly different). Inset: average over 50 different runs. (e) Density profile at a late time [$t=160$ 000, final configuration of (b)]. (f) Same as (c) but for the late-time density profile of (e).</figcaption></figure>

profile widths vanish in the infinite-size limit except if band or sheet structures are present.

In Figs. 15(a) and 15(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 and
<figure><img src="image_16.png" /><figcaption>FIG. 15. (Color online) (a),(b) Time-averaged profile width for both density (a) and order parameter (b) as a function of noise amplitude in the ordered phase (angular noise, 1024×256 domain, global motion along the major axis, $\rho=2$, and $v_0=0.5$). The dashed vertical blue line marks the order-disorder transition. (c),(d) Typical instantaneous profiles along the long dimension of the system described in (a) and (b) for intermediate noise value [(b) $\eta=0.4$] and in the bandless regime [(c) $\eta=0.15$].</figcaption></figure>

sheets stand less sharply out of the disordered background [Fig. 15(c)]. At some point [$\eta \approx 0.3$ for the parameter values considered in Figs. 15(a) and 15(b)], bands rather abruptly disappear and are no longer 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. 15(d)], but can no longer give rise to (meta)stable long-lived transverse structures.

Density fluctuations in the bandless regime are in fact anomalously strong: measuring number fluctuations in subsystems of linear size $\ell$, we find that their root mean square $\Delta n$ does not scale like the square root of $n=\rho\ell^d$, the mean number of particles they contain; rather we find $\sigma(n) \propto n^\alpha$ with $\alpha \approx 0.8$ both in $d=2$ and in $d=3$ [Fig. 16(a)]. 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 $\alpha$, although needing to be refined, does not seem to be compatible with the prediction $\alpha=\frac{1}{2}+\frac{1}{d}$ made in [53].) More work is needed to fully understand under what 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 the predictions of Toner and Tu [8,33,34], the dynamics of the symmetry-broken ordered phase of polar

<figure><img src="image_17.png" /><figcaption>FIG. 16. (Color online) Giant density fluctuation and transverse superdiffusion in the bandless ordered phase. (a) Anomalous density fluctuations (see text): $\Delta n$ scales approximately like $n^{0.8}$ (the dashed line has slope 0.8) both in two dimensions (black circles, $L=256$, $\rho=2$, $v_0=0.5$, angular noise amplitude $\eta=0.25$) and in three dimensions (red triangles, $L=64$, $\rho=1/2$, $v_0=0.5$, vectorial noise amplitude $\eta=0.1$, values shifted for clarity). (b) Average doubling time $\tau_2$ of the transverse (with respect to mean velocity) interparticle distance $\delta_\perp$. Black circles: ordered bandless regime ($\rho=4$, angular noise amplitude $\eta=0.2$ in a rectangular box of size $1024 \times 256$). The black dashed line marks the expected growth $\tau_2 \sim \delta_\perp^{3/2}$. Red squares: same but deep in the disordered phase ($\rho=4$, angular noise amplitude $\eta=1$, $L=512$). The dot-dashed red line shows normal diffusive behavior: $\tau_2 \sim \delta^2$.</figcaption></figure>

active particles should be characterized by a superdiffusive mean square displacement

$$ \Delta x_{\perp} = \sqrt{\langle [x_{\perp}(t) - x_{\perp}(0)]^2 \rangle_i} \quad (13) $$

in the direction(s) transversal to the mean velocity. In particular, in $d=2$ one has [34]

$$ \Delta x_{\perp}^2 \sim t^{\nu} \quad (14) $$

with $\nu=4/3$. While this analytical result has been successfully 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 subclusters of particles moving coherently (as discussed in Sec. IV D). As a consequence, an ensemble of test particles in a cohesionless model is exposed to different “transport” regimes (with respect to 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 $\tau_2$, the average time taken by two particles to double their transverse separation distance $\delta_\perp$. From Eq. (14) one immediately has

$$ \tau_2 \sim \delta_\perp^{2/\nu} \quad (15) $$

with $2/\nu=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. 16(b)] 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, subclusters emerge and propa-
gate ballistically and isotropically due to the absence of a
well-established global order. Particle 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 displace-
ment of particles exhibits the scaling $\Delta x^2 = \langle |\vec{x}(t) - \vec{x}(0)|^2 \rangle_i \propto t^{5/3}$ [21]. In view of our current understanding of the dis-
continuous 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 par-
ticles moving coherently. Even in the case where high-
density bands or sheets are present, these are in fact dynami-
cal 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 re-
gimes present at low noise intensities [Fig. 17(a) and 17(b)].

Thus, bands are truly complex, nontrivial structures
emerging out of the transverse dynamics of clusters with a
well-defined mean size (since μ > 2). It is only in the band-
less regime that one can speak, as do 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 consid-
ered 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 of [7] as soon as the system size is large enough.
Whether in the band and sheet regime or not, P(φ) shows
essentially Gaussian tails, is strongly peaked around its
mean, and its variance decreases with increasing system size
[Figs. 17(c) and 17(d)].

Although the picture of intermittent bursts between “lami-
nar” 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

<figure><img src="image_18.png" /></figure>

FIG. 17. (Color online) (a),(b) Cluster size distributions (arbi-
trary units) for domain sizes L=32 (black), 128 (cyan), and 512
(green) from left to right (d=2, ρ=2, v₀=0.5, angular noise). (c),(d)
Probability distribution functions of the order parameter φ (arbi-
trary units) for the same parameters and system sizes as in (a) and
(b). (The most peaked distributions are for the largest size L=512.)
Left panels (a),(c), η=0.1, bandless regime; Right panels (b),(d),
regime with bands at η=0.4.

intermittency of cluster dynamics in the bandless regime.
Again, these phenomena, reported also in the context of ac-
tive 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 and sheets
are the typical asymptotic structures appearing in *finite* dom-
mains 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 disor-
dered 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 corre-
lation function of both the density and velocity fields. These
fields are defined by a coarse-graining over a small length
scale ℓ (typically 4). These correlation functions have an
unusual shape [Fig. 19(a)]: after some rather fast initial de-
cay, they display an algebraic behavior whose effective ex-
ponent decreases with time, and finally display a near-
exponential cutoff. As a result, they cannot be easily
collapsed on a single curve using a simple, unique, rescaling
length scale. Nevertheless, using the late exponential cutoff,
a correlation length ξ can be extracted. Such a length scale ξ
grows roughly linearly with time [Fig. 19(b)]. Qualitatively
<figure><img src="image_19.png" /><figcaption>FIG. 18. (Color online) Phase ordering from disordered initial conditions ($d=2$, angular noise amplitude $\eta=0.08$, $\rho=1/8$, $v_0=0.5$, system size $L=4096$). Snapshots of the density field coarse grained on a scale $\ell=8$ at times $t=160, 320$, and $640$ from left to right.</figcaption></figure>

similar results are obtained whether or not the noise strength is in the range where bands and sheets appear in finite boxes.

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

## V. GENERAL DISCUSSION AND OUTLOOK

### 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 into 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. Role of bands and 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., with interactions restricted 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 needs to be further investigated, these results indicate that the 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 condi-

<figure><img src="image_20.png" /><figcaption>FIG. 19. (Color online) Phase ordering as in Fig. 18 ($L=4096$, $\rho=1/2$, $v_0=0.5$). (a) Two-point density correlation function $C(r,t) = \langle\rho_{\ell}(\vec{x},t)\rho_{\ell}(\vec{x}+\vec{r},t)\rangle_{\vec{r}}$ (coarse grained over a scale $\ell=4$) as a function of distance $r=|\vec{r}|$ at different time steps: from left to right $t=50, 50, 100, 200, 400, 800$, and $1600$. Noise amplitude is $\eta=0.25$; data have been further averaged over $\approx 40$ different realizations. Inset: log scales reveal the intermediate near-algebraic decay and the quasiexponential cutoff. (b) Length scale $\xi$, estimated from the exponential cutoff positions, as a function of time. Empty black circles: $\eta=0.25$ as in (a) (i.e., regime in which bands are observed asymptotically). Red full triangles: $\eta=0.1$ (i.e., in the bandless regime). The dashed black line marks linear growth.</figcaption></figure>
tions of the domain in which the particles are allowed to move. In rectangular domains with at least one periodic di- rection, these objects can form, span across the whole do- main, and move. But in, say, a circular domain with reflect- ing boundary conditions, they cannot develop freely, being repeatedly frustrated. Nevertheless, simulations performed in such a geometry indicate that the transition remains discon- tinuous, 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 toward 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 from developing 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 and sheets can de- velop 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 gen- eral picture. The key feature of the Vicsek-like models stud- ied 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 information about the order, advection, den- sity fluctuations, and order are intimately linked. High den- sity means strong local order (if the 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 or- der. In practice, such an “equation of state” approach can be justified by looking, e.g., at a scatter plot of local order pa- rameter vs local density. Figure 20 reveals that, in the or- dered bandless regime, such a scatter plot is characterized by a plateau over a large range of local density values corre- sponding 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 per- colate 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 fluctua- tions, the order parameter field is, on the other hand, rather constant; see Fig. 15(c)]. The noise intensity at which bands emerge roughly corresponds to the value where the low- density disordered regions percolate. The remaining discon- nected, dense patches then eventually self-organize into bands or sheets. The emergence of these elongated structures is rather natural: moving packets elongate spontaneously be-

cause they collect many particles; superdiffusion in the di- rections perpendicular to the mean motion endows these na- scent bands and sheets with some “rigidity.” At still stronger noise, the bands and 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 or 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 true only in the bandless regime.

### D. Outlook

The results presented here are almost entirely numerical. Although they were obtained with care, they need to be ul- timately backed up by more analytical results. A first step is the derivation of a continuous description in terms of a den- sity and a velocity field (or some combination of the two), which would allow one 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 fluc- tuations are crucial in the systems considered here if only because some of the effective noise terms will be multiplica-

<figure><img src="image_21.png" /></figure>

FIG. 20. (Color online) Scatter plots of local order parameter vs local density in the ordered phase [angular noise, $\rho=2$, $v_0=0.5$, in a domain of size $1024 \times 256$—the same parameters as in Fig. 15(a)]. The local quantities were measured in boxes of linear size $\ell=8$. Here, the local order is represented by the angle between the orien- tation $\Theta_{\text{local}}$ of the local order parameter and the global direction of motion $\Theta_{\text{global}}$. The black solid lines are running averages of the scatter plots. The red solid lines indicate the global density $\rho=2$ (and thus mark the percolation threshold in a two-dimensional square lattice). (a),(b) $\eta=0.1$ and 0.2: in the bandless regime, the ordered plateau starts below $\rho=2$, i.e., ordered regions percolate. (c) Approximately at the limit of existence of bands: the start of the plateau is near $\rho=2$. (d) At higher noise amplitude in the presence of bands.
tive 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.

## ACKNOWLEDGMENTS

Most of this research has been funded by the European Union via the FP6 StarFLAG project. Partial support from the French ANR Morphoscale project is acknowledged. We thank A. Vulpiani and M. Cencini for introducing us to the method outlined in Ref. [54].

## APPENDIX: FLUCTUATION-DISSIPATION RELATION

In [5], Vicsek et al. also studied the validity of the fluctuation-dissipation theorem and concluded, from numerical analysis, that it is violated. Here we approach this question again, using the “dynamical” approach put forward in [60], rather than the equilibrium used in [5]. The fluctuation-dissipation relation is expressed as

$$R(t - t_0) = \frac{1}{T_{\text{eff}}} \frac{\partial C(t - t_0)}{\partial t}, \qquad (\text{A1})$$

where R is the response, C the associated correlation function, and $T_{\text{eff}}$ some “effective temperature.”

In both cases, at any rate, the key point is to investigate the effect of an external field on the ordering process. Since here one cannot rely on any Hamiltonian structure, the external field remains somewhat arbitrary, as it cannot be unambiguously defined as the conjugate variable of the order parameter.

In [5], the field $\vec{h}$ was directly confronted with the local average velocity, and the governing equation was replaced by

$$\vec{v}_i(t + \Delta t) = v_0(\mathcal{R}_{\eta} \circ \vartheta)\left(\sum_{j \in \mathcal{V}_i} v_j(t) + \vec{h}\right). \qquad (\text{A2})$$

This manner of introducing $\vec{h}$ leads to an effective intensity which depends on the local ordering: $\vec{h}$ is comparatively stronger in disordered regions (an important effect at the early stages of ordering) than in ordered regions. This could in fact prevent the necessary linear regime from occurring, even at very low field values (see, for 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:

$$\vec{v}_i(t + \Delta t) = v_0 \vartheta \left( \sum_{j \in \mathcal{V}_i} \vec{v}_j(t) + \left| \sum_{j \in \mathcal{V}_i} \vec{v}_j(t) \right| \vec{h} + \eta \mathcal{N}_i \vec{\xi} \right), \qquad (\text{A3})$$

where the vectorial noise was chosen because, as shown above, it leads more easily to the asymptotic regime. The

<figure><img src="image_22.png" /><figcaption>FIG. 21. (Color online) Test of the fluctuation-dissipation relation on the vectorial model with repulsive force ($\rho=2$, $v_0=0.3$, $L=128$). (a) Susceptibility vs time at reduced noise amplitude $\epsilon=1 - \eta/\eta_i=0.005$, $|\vec{h}|=10^{-2}$ (dashed red line), and $|\vec{h}|=10^{-3}$ (plain black line). (b) Susceptibility vs correlation at $|\vec{h}|=10^{-3}$ and $\epsilon=0.08, 0.26$, and 0.44 from top to bottom. (c) Effective temperature vs noise amplitude. The vertical dashed line marks the transition point.</figcaption></figure>

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_\perp$. We consider the former. In practice, we quenched, at time $t_0$, a strong-noise, highly disordered system $\varphi(t_0) \approx 0$ to a smaller noise value and started applying the constant homogeneous field $\vec{h}$ immediately. We then followed the response of the system by monitoring the growth of the order parameter. We measured the susceptibility $\chi_\|$, which is nothing but the integrated response function: $\chi_\|(t, t_0) = \int_{t_0}^t R_\|(t, t') dt'$. In practice, we have

$$\chi_\|(t - t_0) = \frac{1}{|\vec{h}|} \vec{\varphi}(t) \cdot \vartheta[\vec{h}]. \qquad (\text{A4})$$

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. 21(a).

Correspondingly, the correlation function is defined as

$$C(t - t_0) = \frac{1}{v_0^2} \langle \vec{v}_i(t_0) \cdot \vec{v}_i(t) \rangle_i. \qquad (\text{A5})$$

The fluctuation-dissipation relation (A1) can then be written in its integrated form:

$$\chi_\|(t - t_0) = \frac{1}{T_{\text{eff}}} [C(0) - C(t - t_0)]. \qquad (\text{A6})$$

In Fig. 21(b), we show that $\chi_\|$ and C are related linearly in time, confirming the validity of this relation and allowing an estimation of $T_{\text{eff}}$.

This well-defined—although not uniquely defined— effective temperature varies as expected in parameter space. In particular, it increases systematically with the noise strength $\eta$ [Fig. 21(c)], although this variation is not linear. Note also, that, intriguingly, there is a small jump of $T_{\text{eff}}$ at the noise value corresponding to the transition in this case ($\eta_t \approx 0.55$).
[1] Pliny, *Natural History*, translated by H. Rackham (Harvard University Press, 1968), Vol. III, book 10, xxxii, p. 63.

[2] *Animal Groups in Three Dimensions*, edited by J. K. Parrish and W. M. Hamner (Cambridge University Press, Cambridge, U.K., 1997).

[3] E. V. Albano, Phys. Rev. Lett. **77**, 2129 (1996).

[4] I. D. Couzin and J. Krause, Adv. Stud. Behav. **32**, 1 (2003).

[5] A. Czirók, H. E. Stanley, and T. Vicsek, J. Phys. A **30**, 1375 (1997).

[6] G. Grégoire, H. Chaté, and Y. Tu, Phys. Rev. E **64**, 011902 (2001).

[7] C. Huepe and M. Aldana, Phys. Rev. Lett. **92**, 168701 (2004).

[8] J. Toner and Y. Tu, Phys. Rev. Lett. **75**, 4326 (1995).

[9] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. **75**, 1226 (1995).

[10] N. D. Mermin and H. Wagner, Phys. Rev. Lett. **17**, 1133 (1966).

[11] E. Bertin, M. Droz, and G. Grégoire, Phys. Rev. E **74**, 022101 (2006).

[12] B. Birnir, J. Stat. Phys. **128**, 535 (2007).

[13] H. J. Bussemaker, A. Deutsch, and E. Geigant, Phys. Rev. Lett. **78**, 5018 (1997).

[14] H. Chaté, F. Ginelli, and R. Montagne, Phys. Rev. Lett. **96**, 180602 (2006).

[15] I. D. Couzin, J. Theor. Biol. **218**, 1 (2002).

[16] I. D. Couzin, J. Krause, N. Franks, and S. Levin, Nature (London) **433**, 513 (2005).

[17] Z. Csahók and T. Vicsek, Phys. Rev. E **52**, 5297 (1995).

[18] A. Czirók, A.-L. Barabási, and T. Vicsek, Phys. Rev. Lett. **82**, 209 (1999).

[19] A. Czirók, M. Vicsek, and T. Vicsek, Physica A **264**, 299 (1999).

[20] Y. L. Duparcmeur, H. Herrmann, and J. P. Troadec, J. Phys. I **5**, 1119 (1995).

[21] G. Grégoire, H. Chaté, and Y. Tu, Phys. Rev. Lett. **86**, 556 (2001).

[22] G. Grégoire, H. Chaté, and Y. Tu, Physica D **181**, 157 (2003).

[23] J. Hemmingsson, J. Phys. A **28**, 4245 (1995).

[24] H. Levine, W.-J. Rappel, and I. Cohen, Phys. Rev. E **63**, 017101 (2000).

[25] A. S. Mikhailov and D. H. Zanette, Phys. Rev. E **60**, 4571 (1999).

[26] A. Mogilner and L. Edelstein-Keshet, J. Math. Biol. **38**, 534 (1999).

[27] O. J. O'Loan and M. R. Evans, J. Phys. A **32**, L99 (1999).

[28] M. R. D'Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Phys. Rev. Lett. **96**, 104302 (2006).

[29] N. Shimoyama, K. Sugawara, T. Mizuguchi, Y. Hayakawa, and M. Sano, Phys. Rev. Lett. **76**, 3870 (1996).

[30] R. A. Simha and S. Ramaswamy, Physica A **306**, 262 (2002).

[31] R. A. Simha and S. Ramaswamy, Phys. Rev. Lett. **89**, 058101 (2002).

[32] B. Szabó, G. J. Szolosi, B. Gonci, Z. Juranyi, D. Selmeczi, and T. Vicsek, Phys. Rev. E **74**, 061908 (2006).

[33] J. Toner and Y. Tu, Phys. Rev. E **58**, 4828 (1998).

[34] Y. Tu, J. Toner, and M. Ulm, Phys. Rev. Lett. **80**, 4819 (1998).

[35] C. M. Topaz and A. L. Bertozzi, SIAM J. Appl. Math. **65**, 152 (2004).

[36] C. Topaz, A. Bertozzi, and L. M.A., Bull. Math. Biol. **68**, 1601 (2006).

[37] T. Vicsek, A. Czirók, I. J. Farkas, and D. Helbing, Physica A **274**, 182 (1999).

[38] G. Grégoire and H. Chaté, Phys. Rev. Lett. **92**, 025702 (2004).

[39] M. Nagy, I. Daruka, and T. Vicsek, Physica A **373**, 445 (2007).

[40] M. Aldana, V. Dossetti, C. Huepe, V. M. Kenkre, and H. Laralde, Phys. Rev. Lett. **98**, 095702 (2007).

[41] S. Lübeck, Int. J. Mod. Phys. B **18**, 3977 (2004).

[42] P. Marcq, H. Chaté, and P. Manneville, Phys. Rev. Lett. **77**, 4003 (1996).

[43] P. Marcq, H. Chaté, and P. Manneville, Phys. Rev. E **55**, 2606 (1997).

[44] K. Binder, in *Phase Transitions and Critical Phenomena*, edited by C. Domb and M. S. Green (Academic Press, 1976).

[45] *Finite Size Scaling and Numerical Simulations of Statistical Systems*, edited by V. Privman (World Scientific, Singapore, 1990).

[46] K. Binder, Rep. Prog. Phys. **60**, 487 (1997).

[47] B. Efron, *The Jackknife, The Bootstrap and Other Resampling Plans* (SIAM, Philadelphia, 1982).

[48] J. C. Niel and J. Zinn-Justin, Nucl. Phys. B **280**, 355 (1987).

[49] C. Borgs and R. Kotecky, J. Stat. Phys. **61**, 79 (1990).

[50] H. Chaté, F. Ginelli, and G. Grégoire, Phys. Rev. Lett. **99**, 229601 (2007).

[51] S. Mishra and S. Ramaswamy, Phys. Rev. Lett. **97**, 090602 (2006).

[52] V. Narayan, S. Ramaswamy, and N. Menon, Science **317**, 105 (2007).

[53] S. Ramaswamy, R. A. Simha, and J. Toner, Europhys. Lett. **62**, 196 (2003).

[54] G. Boffetta, A. Celani, M. Cencini, G. Lacorata, and A. Vulpiani, Chaos **10**, 50 (2000).

[55] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. **49**, 435 (1977).

[56] J. Toner, Y. Tu, and S. Ramaswamy, Ann. Phys. (N.Y.) **318**, 170 (2005).

[57] E. Bertin, M. Droz, and G. Grégoire, (unpublished).

[58] Z. Csahók and A. Czirók, Physica A **243**, 304 (1997).

[59] E. Ben-Jacob, I. Cohen, and H. Levine, Adv. Phys. **49**, 395 (2000).

[60] L. F. Cugliandolo, J. Kurchan, and L. Peliti, Phys. Rev. E **55**, 3898 (1997).

[61] E. Bertin (private communication).