The Physics of the Vicsek Model

In these lecture notes, prepared for the Microswimmers Summer School 2015 at Forschungszentrum Juelich, I discuss the well known Vicsek model for collective motion and its main properties. In particular, I discuss its algorithmical implementation and the basic properties of its universality class. I present results from numerical simulations and insist on the role played by symmetries and conservation laws. Analytical arguments are presented in an accessible and simplified way, but ample references are given for more advanced readings.


Introduction
Collective motion (or flocking) is a ubiquitous phenomenon, observed in a wide array of different living systems and on an even wider range of scales, from fish schools [1] and mammal herds [2] to bacteria colonies [3] and cellular migrations [4], down to the cooperative behavior of molecular motors and biopolymers at the subcellular level [5]. The aerial displays of starling flocks and other social birds are of course among the most spectacular examples, and have attracted the interest of speculative observers for quite a long time [6].
To the physicist eye, these phenomena are also highly nontrivial because they occur far from equilibrium, as single constituent particles in a flock (whether they are birds, bacteria or cells) are active, i.e. they continuously dissipate free energy to perform systematic (i.e. non-thermal) motion. Also, collective motion often arises spontaneously, without any leader, external field or geometrical constraint guiding the process. In a more technical language, we may say that ordered motion follows from the spontaneous breaking of a continuous symmetry, viewing a collectively moving flock as an orientationally ordered phase of active matter [7].
Collective motion phenomena, of course, are not restricted to living matter, and in recent times they have been studied in various experimental systems, such as active colloids [8] and driven granular matter [10,9].
The ubiquity of collective motion phenomena at all scales, from groups of large vertebrates to subcellular collective dynamics, strongly hints at the existence of some universal features, possibly shared among the many different situations, regardless of many individual-level details. One way of approaching these problems is to construct and study minimal models of collective motion, that is models stripped of as many details as possible and only equipped with the basic features that we believe characterize the problem, typically its fundamental symmetries and conservation laws. This approach is fundamentally justified by hydrodynamics considerations, by which a great deal of microscopic details may be ignored, at least if we are interested in the large wavelength and long time behavior of our system [11].
In any case, even if one is interested in finer, non-asymptotic details, it is surely good practice, before starting toying with your favourite model, to first understand the underlying, long wavelength physics inevitably shared by all systems with the same fundamental features.
In these notes, I will introduce and discuss in details the properties of the Vicsek model -the simplest off-lattice model describing a flocking state -and of the related Vicsek class. Approaching the study of collective motion, it is important to understand that all physical systems and models sharing the same basic features with this class will also display the same asymptotic properties. The only way to escape this, is to alter some fundamental property of the system, like changing the broken symmetry (for instance from polar to nematic symmetry 1 ) or adding a further conservation law (for instance momentum conservation, which is relevant for most active suspensions). This is likely going to be the main message of this lecture.

The Vicsek model
The Vicsek model (VM) is perhaps the simplest model displaying a transition to collective motion; in the study of active matter plays a prototypical role, similar to the one played by the Ising model for equilibrium ferromagnetism. Its simple dynamical rule has been adopted as the starting point for many generalizations and variations which have been applied to a wide range of different problems. The Vicsek model has been originally introduced 20 years ago by the pioneering work of Vicsek and coworkers [12]. Subsequent numerical studies (see for instance Refs. [13] and [14]) greatly helped in clarifying its properties.

Definition and physical features
The model describes the overdamped dynamics of a collection of N selfpropelled particles (SPPs) characterized by their off-lattice position r t i and direction of motion (or heading) s t i , a unit vector, |s t i | = 1. Here i is the particle index, i = 1, . . . , N, and t labels time. All particles move with the same constant speed v 0 , according to the time-discrete dynamics so that orientation s and particle velocity v = v 0 s coincide but for a multiplicative constant (and often the term velocity is also used for the orientation s).
Particles tend to align their direction of motion with the one of their local neighbours, and s t i depends on the average direction of all particles (i included) in the spherical neighborhood S i of radius R 0 centered on i. Indeed, in the Vicsek algorithm the alignment with ones neighborhood is almost perfect, only hampered by a white noise term which plays a role analogous to the one of a temperature in equilibrium systems. In two spatial dimensions (d = 2), the direction of motion is defined by a single angle θ t i , with s = (cos θ, sin θ), and one may simply write the orientation dynamics as where ξ t i is a zero average, delta-correlated scalar noise uniformely distributed in [−π, π] 2 . Such a noise is often called white, since it has a flat Fourier spectrum.
In Eq. (2), the function Arg returns the angle defining the orientation of the average vector j n t ij s t j , and n t ij is the connectivity matrix, This way of chosing neighbours is sometimes defined as metric, being based on the metric notion of distance. The dynamics (1)-(2), depicted in Fig. 1A, is synchronous, meaning that all particles positions and headings are adjusted at the same time. In studying this model, one can always chose a convenient set of space and time units, such that ∆t = R 0 = 1 and the model behavior only depends on three control parameters: the noise amplitude η, the particles speed v 0 and the total density of particles ρ 0 = N/V , where V is the volume of the system. Being interested in the bulk properties of a system, one typically assumes periodic boundary conditions, so that V = L d , with L being the linear system size. In numerical simulations, periodic boundary conditions help to minimize finite size effects due to finite boundaries, and in the following we will implicitly assume them unless stated otherwise.
In the literature, one may find a number of slightly different flavours of the algorithm defined above. For instance, the noise in Eq. (2) may be distributed according to a Gaussian, a small, short ranged repulsion force between particles may be included to account for volume exclusion, or the position r i at time t + ∆t , as defined in Eq. (1), may be determined by the direction of motion at time t and not at time t + ∆t (indeed, this was actually the choice made in the original paper by Vicsek and coworkers [12]). However, typically all these differences do not matter much, and do not change the physical properties of the Vicsek model. On the other hand, there are some features which are essential, and define what we call the Vicsek class. It is worth discussing them explicitly: • Spontaneous symmetry breaking to polar order. Eqs. (1)-(2) are isotropic in space, as no preferred direction is given a priori. However, Eq. (2) contains an explicit polar (or ferromagnetic) alignment term. If this alignment term is strong enough to overcome the effect of the noise (or to put it differently, if the noise amplitude η is low enough), the system may develop global orientational order and thus collective motion, signaled by a finite polar order parameter (or center of mass velocity) Figure 1: (A) Cartoon depiction of the Vicsek dynamics. The particle in red aligns imperfectly with neigbours inside the metric range R 0 (dark blue) and then moves along its heading direction. Note that neighbours may change as a consequence of movement. In the cartoon only the red particle adjusts its position and heading, but in the full algorithm all particles update synchronously. (B) Molecular dynamics algorithm. Once the system is divided in boxes of linear size R 0 , we know immediately that the candidate neighbouring particles of the red particle are restricted to the ones laying in the nine adjacent boxes inside the red square. These particles are marked in light blue, but testing their distance from the red particle we find that only some of them (dark blue) are actually closer than R 0 (note that even a particle laying in the same box as the red one may occasionally be farther away than R 0 ).
an analogous of the total magnetization in spin systems. The value of the modulo of the polar order parameter ϕ(t) = |ϕ(t)| is essentially determined (minus fluctuations and finite size effects) by the three control parameters ρ 0 , η and v 0 . Its stationary time average φ = ϕ(t) t is typically used to describe the spontaneous symmetry breaking phenomenon, with φ > 0 in the ordered phase 3 .
Its orientation in the ordered phase, on the other hand, is not determined a priori, and all directions are equally likely (the one picked up at a given time being chosen by fluctuations). Since the orientation can change continuously in space 4 , in the transition to collective motion a continuous symmetry is spontaneously broken. This has a number of important consequences that will be explored in this notes.
• Self-propulsion and local alignment interactions. Particles are selfpropelled, that is, they move according to Eq. (1). In particular, they 3 There are of course finite size effects, and in the disordered phase the vectors s i do not cancel exactly one with each other, leading to |φ| ∼ 1 √ N . 4 As opposed, for instance, to the Ising model in which spins may only assume two values, ±1.
change their relative position according to their velocity fluctuations δs t i = s t i − ϕ(t). Thus, the connectivity matrix in Eq. (2) is not static, but it changes in time in a nontrivial way. This is exactly where the far from equilibrium nature manifests in the Vicsek model, as it will be discussed in section 3.3. Of course, the connectivity matrix will change as a consequence of particle motion only if the interactions are local, that is, if R 0 << L.
• Conservation laws. The only conservation law of the VM is the conservation of the total number of particles, that is, our birds do not die or get born on the fly. There are no other conservation laws, and in particular it should be noted that momentum is not conserved. Our self-propelled particles are thought to be moving over a dissipative substrate (or in a viscous medium) which acts as a momentum sink. This is of course not the case of a particle swimming in a three-dimensional suspensions 5 , where momentum is transferred from the swimmers (typically exerted as a force dipole) to the surrounding fluid, and long-range hydrodynamic interactions are probably relevant (and for man-made micro swimmers or self-propelled nano-rods they are typically the only interactions!).
As a consequence of the lack of momentum conservation, also Galileian invariance is broken. In fact, the VM is explicitly formulated in the reference frame in which the dissipative substrate is at rest, and it is not invariant under any arbitrary velocity shift.
All together, the features discussed above define the Vicsek class. Finally, we have to remark another obvious feature of the Vicsek model: all particles move with the same speed v 0 . However, to a certain extent it is possible to relax this conditions staying inside the Vicsek class. For instance, one can let the individual speeds fluctuate in some bounded interval without changing the model asymptotic properties.

Vectorial noise and Vicsek model in three space dimensions
In the literature, it is possible to find different ways of implementing the noise in the equation for the orientation dynamics. In the past, some attention has been given to the so called vectorial noise [13] (as opposed to the noise used in Eq. (2) which is sometimes called scalar or angular). One may replace Eq. (2) by where m i = j n ij is the number of interacting neighbours, and ξ is a random unit vector, delta-correlated in time and in the particle index. The denominator in the r.h.s. of Eq. (6) is a normalization term to ensure that |s| = 1. If one interprets the scalar noise of Eq. (2) as an error the SPP makes trying to take the (perfectly determined) mean direction of motions of his neighbours, the vectorial version of the noise can be thought as the sum of the errors made while trying to assess the direction of motion of the interacting neighbours 6 . Certain literature, also refers to these two noises implementation as (respectively) intrinsic and extrinsic, but it is important to stress that these two implementations do not yield different asymptotic properties, even if their finite size behavior may be slightly different (more later on this).
It is however worth remarking that Eq. (6) can be directly extended to any spatial dimension, while starting from Eq. (2) requires some more care. In order to write a Vicsek dynamics with scalar noise in d = 3, one has to introduce a rotation operator R η performing a random (and of course delta-correlated) rotation uniformly distributed around the argument vector, In d = 3, for instance, R η [s] will lay in the solid angle subtended by a spherical cap of amplitude 4πη and centered around s.

Limiting cases
It is instructing to consider the relations between the Vicsek model and some well-known models of equilibrium statistical physics. Obviously, the VM may be seen as an XY (or Heisenberg in d = 3) ferromagnet in which particles are not fixed in some lattice positions but can actually move along the spin direction. Indeed, the XY or Heisenberg equilibrium models can be formally recovered in the case v 0 → 0, where particles do not move at all and n ij is fixed once for all. If the local connectivity of the static connection network its dense enough, its dynamics converges to the equilibrium distribution of an XY or Heisenberg model, with a temperature T that is a monotonic function of the noise amplitude η. This is however a singular limit 7 .
Another way of looking at the Vicsek model is to see it as a persistent random walk in which particles may align their directions of motion one with each other by some local interaction rule. In continuous time and d = 2, the persistent random walk can be written aṡ (with ξ i being some white noise) which is just a Vicsek dynamics without the alignment interaction term (starting from Vicsek dynamics, Eq. (8) can be formally obtained by taking the limit R 0 → 0). Once again, this is a singular limit, and a collection of non-interacting persistent random walker has an equilibrium distribution with some temperature given by the noise term. The opposite limits, v 0 → ∞ and R 0 → ∞, also correspond to singular cases. As already mentioned, if the interactions are long ranged, the system is globally coupled and the connectivity matrix n ij is trivially static. In this way, motion is completely decoupled from long-ranged alignment, and most (if not all) of the fascinating Vicsek model properties are lost.
The infinite speed limit v 0 → ∞, on the other hand, just produces a random rewiring of the connectivity network: if v 0 >> L, any small fluctuation δs t i in the orientation will push nearby particles infinitely apart. In a system with periodic boundary conditions this is equivalent to random rewiring of interactions, another trivial case in which motion decouples from alignment.
The bottom line is that, while is interesting to understand the relations between the VM and its limiting cases, it is not possible in general to deduce properties of the former from the study of the latter (singular) limiting cases [15].

Algorithmic implementation
The Vicsek model is extremely simple and particularly well suited for numerical studies, as Eqs. (1)-(2) can be easily implemented on a computer. However, it should be noted that a straightforward implementation of the metric neighbouring condition (4) would require testing the distance of all i − j couples, an operation scaling with system size as order N 2 . This approach would quickly become unmanageable as the number of SPPs N grows, making practically impossible to run simulations with more than a few thousands particles.
There is of course a way around this problem, based on techniques originally developed for the study of molecular dynamics. The idea is rather simple, even if its algorithmic interpretation may not be so straightforward. One should ideally divide the system volume L d in boxes of linear size R 0 (remember that one can always rescale space so that R 0 = 1), assigning at each timestep each particle to a given box. Once this is done, it is clear that for any given particle i, all other particles laying outside the box containing i and its next neighbouring boxes cannot be closer than R 0 . Therefore, one immediately and effortlessly reduces its search to a handful of boxes per particle. In d = 2 one has to only look into 9 boxes (the general formula in spatial dimension d of course gives 3 d boxes). A sketch for this algorithm is depicted in Fig. 1B.
At any fixed total density 8 , the mean number of particles contained in these boxes does not grow with N, so that the number of operation needed to find all the interacting couples grows only linearly with N. Since also assigning particles to boxes is an order N operation, it is immediate to conclude that the entire molecular dynamics algorithm computational time is of order N rather than N 2 as the naive algorithm. A huge improvement if one is interested in asymptotic (i.e. long time and large N) properties.
Any serious numerical study should employ molecular dynamics algorithms. Current state of the art simulations of Vicsek model involve from a few millions to a few tens of millions of particles.

Physical properties
We now proceed to discuss the main physical properties exhibited by the Vicsek class. As we shall see, they mostly emerge from the intriguing interplay between particles self propulsion and the spontaneous symmetry breaking characterizing the ordered state.

Transition to collective motion and phase separation
Numerical simulations easily show that the Vicsek model display a transition from disorder to ordered collective motion. For instance, as the noise amplitude η is decreased below a certain threshold (and both ρ 0 and v 0 kept fixed), particles start to synchronize their heading and to move together. Starting from disordered initial condistions, this coarsening process is relatively fast, and the size ℓ d of ordered domains grows linearly in time, ℓ d ∼ t [14].
The easiest way to capture the transition to collective motion is to monitor the order parameter ϕ (the center of mass velocity) defined in Eq. (5). At high noise amplitudes, SPPs are unable to synchronize their headings, which tend to cancel out in the sum i s i . It can be shown that the sum of N randomly oriented unit vectors has a modulo of order √ N, so that in the disordered phase the scalar order parameter ϕ ∼ 1 √ N , or essentially zero for any large number N of SPPs.
At lower noise amplitudes, below a certain threshold η c , the system undergoes a spontanous symmetry breaking phase transition as SPPs synchronize their heading. The scalar order parameter becomes finite and roughly of order one (note that perfect order exactly implies ϕ = 1). This is resumed in Fig. 2a, where the long-time (or stationary) average φ = ϕ(t) t is shown for different noise amplitudes. The parameter that is varied as the system goes through the symmetry breaking is referred to as control parameter. The threshold noise amplitude value for the onset of collective motion is, of course, not independent from the other model parameters, and one has η c = η c (ρ 0 , v 0 ).
One simple way to understand the onset of collective motion is to consider that, in order to synchronize the heading of all SPPs, information should be able to propagate through the entire system. While alignment interactions between particles produce such information, noise clearly destroys it. A simple mean-field like argument can then be put forward for low densities. To simplify things, lets rescale our units so that the interaction range is one, R 0 = 1. If ρ 0 << 1 particles are often isolated, and their relatively rare interactions can be treated as instantaneous collisions 9 from which particles emerge agreeing on their headings. The distance ℓ that a particle travels between collisions, the mean free path, scales as ℓ ∼ 1/ρ 0 . Information can propagate through the system only if the mean free path is larger than the SPP persistence length ℓ p , that is the distance a particle can travel before losing its out-of-collision heading. At the onset of order one expects these two quantities to have the same magnitude, ℓ ∼ ℓ p . Given that the persistence length is inversely proportional to noise variance, a relation that has been numerically verified for ρ 0 << 1 (at least in d = 2), 9 To be more precise, we want the mean inter-collision time to be much larger than the time two nearby particles spend at a mutual distance shorter than R 0 . For more details, see reference [14].
and that defines a critical line in the (η, ρ 0 ) plane. This implies that one can also use the total density as a control parameter, keeping the noise amplitude fixed. In this case, one crosses to collective motion as the density is increased. At a first glance, one may think that the symmetry breaking transition to collective motion should be similar to the transition to order in an equilibrium spin system, leading directly to some homogeneously ordered state. This is however not the case, due to the interplay between local order and local density induced by motion. Moving particles, indeed, may gather in high density patches, increasing in turn the number of interacting neighbors, i.e. particles with a mutual distance smaller than R 0 . Locally high density has a positive feedback on the efficiency of the alignment interaction, so that high density patches may be able to locally align while the rest of the systems does not; this is something that cannot happen in an equilibrium spin system! One can indeed show that this feedback mechanism inevitably leads to a long wavelength instability near the onset of order [16,17], that destabilizes the homogeneous ordered phase and leads to (spontaneous) phase separation. For the polar symmetry of the Vicsek class, these phase separation takes the form of high-density ordered bands that travel in a low-density sea of dis-ordered particles [13,14] (see Fig. 2c). Bands extend transversally to the direction of motion and are characterized by a well-defined width, so that it is possible to accomodate several in the same system. Indeed, on very large timescales they seem to settle in a regularly spaced pattern, leading to a smectic arrangement of traveling ordered bands [18]. In d = 3, simple symmetry considerations imply that these structures manifest as sheets, again extending transversally to the direction of motion (Fig. 2d).
At lower noise values or larger densities, the long wavelength instability disappears, and a second transition leads to a homogeneous ordered phase. The resulting Vicsek class phase diagram, sketched qualitatively in Fig. 2e, is thus composed of three phases. A disordered one, akin to a collection of persistent random walkers, a phase-separated ordered regime, characterized by high density ordered bands, and finally an homogeneous ordered phase. In the latter two phases, the rotational symmetry is spontaneously broken and the system exhibit collective motion.
As a consequence of phase separation, the symmetry breaking transition to polar order is a first order one, rather than a second order critical one [13,14]. At the onset of order, the system is bistable, and alternates between the disordered regime and the appearance of a single ordered band, with the order parameter ϕ(t) showing the corresponding jumps characteristics of phase coexistence and of first order phase transtions (Fig. 2b).
The transition to collective motion can also be interpreted as a liquidgas transition, albeit in a non-equilibrium context and with no accessible supercritical region [18].
This phase diagram, with phase separation and a first order transition characterizing the onset of order, is rather generic; it is indeed common not only to the entire Vicsek class, but also to systems whose broken state is characterized by a different symmetry (although details of the phase separated regime may change with different symmetries). However, the existence of this phase separated regime has proven rather elusive, and it took a decade from the first introduction of the Vicsek model to discover it. In fact, the long-wavelength instability leading to phase separation is characterized by a rather large instability wavelength Λ c , so that in systems not too large, where L < Λ c , phase separation cannot be observed and the transition may be mistakenly thought to be continuous and critical. It is only when L is sufficiently larger than Λ c that the true asymptotic behavior of the Vicsek model emerge.
The instability wavelength -of course -depends on model parameters and, to make things worse, also on non-universal details such as the noise implementation. In particular, it is rather larger in systems with a scalar noise (as in Eq. 2) than in systems with a vectorial one (as in Eq. 6), so that it is not uncommon to be able to observe bands only in systems with several hundred thousands of particles or more. Moreover, Λ c seem to diverge both in the low density and in the low speed limits [14]. These difficulties (one can say that the VM is characterized by very strong finite size effects) have fueled a long debate on the order of the phase transition, on the genericity of the phase separated band regime and on the eventual difference between scalar and vectorial noise models in the thermodynamic limit. Careful finite size analisyis and large scale simulations, together with the study of hydrodynamic theories [16,17] for the Vicsek class however, have gathered convincing evidence over the last decade. Nowadays there is a general consensus for the scenario detailed above: no asymptotic difference between scalar and vectorial models, first order transition to collective motion and genericity of the phase separation scenario. 10 We conclude this section noting that moving bands quite similar to VM ones have been observed in in vitro experiments with motility assays, i.e. in a mixture of molecular motors and actin filaments which are among the constituents of cellular cytoskeleton [5]. Such a systems is of course much more complicated than the VM, but is still characterized by self-propulsion (due to the molecular motors) and may undergo a spontaneous symmetry breaking thanks to filament interactions which are effectively aligning. These experimental result demonstrate the power of the minimal model approach.

Topological Vicsek model
A relevant change of the Vicsek rule (2) is given by topological interactions [19]. In topological models, one choses interacting neighbours not as the SPPs lying inside a metric range R 0 , but on the basis of some local topological (or metric-free) rule, such as the n c nearest neighbours or the Voronoi neighbours 11 It is important to stress that this is still a local interaction rule, albeit in the topological rather than metric sense. These choices are motivated by experimental evidence, gathered in starling flocks [20] and in other social vertebrates [21], that individual do not interact with neighbours chosen inside a certain fixed range, but rather with a more or less fixed number of neighbours regardless of local density. One can think that visual perception, limited by occlusions to the first shell of neighbours, is better modeled by topological rather than metric interactions.
In topological models, fluctuations in local density do not affect the interaction frequency or the number of interacting neighbours, so that there is no positive feedback on the efficiency of the alignment interaction. In the absence of an interaction range, it is indeed possible to rescale lengths in order to always have a unit total density, ρ 0 = 1.
Moreover, the long-wavelength instability which destabilizes the homogeneous ordered phase at the onset of order is not present in models with topological interactions, and phase separation is removed [22]. The corresponding phase diagram is much simpler and independent from total density. As the noise is lowered, one directly crosses from disorder to an homogeneously ordered, collectively moving, phase. In the absence of phase separation, the transition is a second order, continuous one, characterized by a novel set of critical exponents.

Long range order in d=2
We now turn our attention to the homogeneously ordered phase.
One interesting and, to a certain extent, surprising property emerging from numerical simulations of the Vicsek models, is its ability to display true collective motion in d = 2, that is, to have a true long range ordered phase in which the order parameter ϕ t is finite for any system size. This is in apparent contradiction with a well-known theorem due to Mermin and Wagner (MW) [23], stating that no system breaking a continuous symmetry in two spatial dimensions may achieve long range order (LRO). A classical example of this theorem is given by the XY model in d = 2. In this case, the system may only achieve a lesser kind of order, called quasi long range order (QLRO), where the order parameter decays algebraically with the number of spins N, albeit with a very small exponent 12 . While this means that, strictly speaking, no order is present asymptotically, a trace of order can still be found in the algebraically decaying spin-spin correlation function, and it is thus possible to formally define a phase transition -the Kosterlitz-Thouless (KT) transition [24].
There is of course a caveat, since the Mermin-Wagner theorem only applyes to equilibrium systems, and the Vicsek model is of course out-ofequilibrium. It is however interesting to understand why the ability of Vicsek particles to move can beat the MW theorem. It is instructive to consider a simplified argument -first introduced by John Toner in his lecture notes on flocking - [25] showing that the VM does so thanks to more efficient information transfer mechanisms.
First consider XY spins on a d dimensional lattice. Suppose they all point in the same direction s , with the exception of a single "mistaken" spin, that lies at an angle δθ 0 from all the others. How this mistake will evolve in time on our lattice? Ferromagnetic alignment cannot simply reset δθ 0 to zero: all it can do is to "iron it out", spreading it to nearby lattice sites. On a lattice, this propagation mechanism is purely diffusive, ∂ t δθ ∼ ∇ 2 δθ, and in a time τ the original error will spread out over a distance r ∼ √ τ , or a volume V e ∼ τ d/2 . Since the total error inside the volume is conserved, the error per spin decays as δθ ∼ δθ 0 /V e ∼ δθ 0 τ −d/2 . This is what happens to a single mistake. However, noise fluctuations constantly produces local errors with a number of errors per spin proportional to time. In a propagation volume V e one has n d ∼ τ V e ∼ τ 1+d/2 errors. Their combined root mean square, according to central limit theorem, is Ω e ∼ √ n e ∼ √ τ V e . We are finally in the position to compute the total error amplitude per spin, that is If d > 2, Eq. (10) predicts that fluctuation errors per spin should decay algebraically in space. This means that order is resistant to fluctuations, and the system displays long range order. On the other hand, if d < 2, fluctuations grows algebraically in space, so that no global order is possible. The case d = 2 is marginal, with a zero algebraic exponent but a logarithmic divergence in the system size L 13 . In this case, fluctuations are still unbounded, but only logarithmically, so that the order is destroyed extremely slowly and the equilibrium system displays QLRO. Note that the fact that we are breaking a continuous symmetry is essential to this argument. Only in this case, in fact, arbitrary small fluctuations can induce an arbitrarily small mistake δθ 0 in spin orientation.
In the VM, however, orientation fluctuations are coupled to motion. Indeed, fluctuations induce a separation between particles of order δx ⊥ ∼ v 0 τ sin ∆θ ∼ τ ∆θ in the directions transversal to the mean direction of motion, and δx ∼ v 0 τ (1 − cos ∆θ) ∼ τ ∆θ 2 in the longitudinal direction, so that two different mechanisms compete to transport orientation information: particle motion and standard diffusion. The propagation volume is readily decomposed in its transversal and longitudinal directions V e ∼ w d−1 ⊥ w (see Fig. 3a), where we have so that the error per spin in the Vicsek model is given by The three equations (11)-(13), where we have introduced the three unknown exponents γ, γ ⊥ and γ , should be solved simultaneously. They yield a system of three linear equations in the three unknown exponents so that above the upper critical dimension d c = 4 the sistem is fully diffusive and γ < 0. For 7/3 ≤ d < 4 transversal propagation is superdiffusive and we have γ = 3 − 2d 2(d + 1) , γ ⊥ = 5 2(d + 1) and γ = 1 2 (16) with again a negative γ. Finally, for d < 7/3 our simple argument also predicts superdiffusion propagation also in the longitudinal direction: which gives γ < 0 for any d > 1, so that orientation fluctuations are suppressed on large scales and the VM can attain long ranger order in any d > 1, thanks to the non-equilibrium, self propelled nature of its particles 14 The fact that, below the upper critical dimension d c = 4, particle motion dominates over simple diffusion -resulting in a superdiffusive propagation -is related to the so called breakdown of linearized hydrodynamics. This phenomenon can be studied more rigorously by a dynamical renormalization group (DRG) study of the hydrodynamic equations for the Vicsek universality class, first obtaines by Toner & Tu by by symmetry arguments [26,27,28,29]. Their detailed analysis clearly lies out of the scope of this notes, but it is worth mentioning that DRG calculations suggest that it is only in the transversal direction that particle motion dominates over simple diffusion. This consideration forces γ = 1/2 in the above argument. This invalidates Eq. (17) and extends Eq. (16) below d = 7/3, yelding γ < 0 and thus LRO in any dimension larger than d = 3/2. Finally, we also note that generically w ⊥ >> w so that fluctuations propagate much slower in the longitudinal directions than in the transversal ones (see Fig. 3a). This spatial anisotropy is of course due to the symmetry breaking process. Once a direction of motion is picked up, spatial isotropy is broken and the longitudinal direction can have different scaling properties from the transversal ones.

The Toner & Tu phase: scale free correlations and anomalous density fluctuations
The homogeneous ordered phase of the Vicsek class is sometimes referred to as the Toner & Tu phase, after the authors of the pioneering papers that first discussed its hydrodynamic behavior [26,27,28]. In this section we briefly discuss its most important properties, which hold for both metric and topological interactions. It is well known that in systems where a continuous symmetry is spontaneously broken, the entire ordered phase is characterized by an algebraic decay of its connected correlation functions (i.e. the corresponding fluctuations correlation function) [11,30]. This is also true for the Vicsek model; moreover, by virtue of the coupling between orientation and local particle density, both the density-density and the orientation-orientation connected correlation functions show an algebraic decay. In particular, it is instructive to consider orientation fluctuations δs i = s i − 1 N i s i . Their equal time, two points correlation function is defined as where r ij is the distance between particle i and j and · is an average over  [14].
realizations (or time in a stationary states). It can be shown that one has C s (r) ∼ r −χ .
In systems of finite linear size L, due to the global constraint i δs i = 0, the correlation function has a zero, which can be used as a finite-size definition of the correlation length ξ, C s (r = ξ) = 0. As a consequence of the spontaneous symmetry breaking one has ξ ∼ L, i.e. the correlation length scales with the system size. In finite systems one can thus write C s (r, L) = r −χ g r ξ (19) where g(r) is a universal scaling function with g(1) = 0.
We have just shown that in the Vicsek class orientation (or velocity) fluctuations are scale free. While a rigorous demonstration is beyond the scope of these notes, it is important to remark that this is just a consequence of the spontaneous breaking of a continuous symmetry. The concept of scale free correlations in collective motion, in fact, received a certain attention after they have been measured in starling flocks observed in the wild [31].
The algebraic nature of correlation functions has a number of other nontrivial consequences. The so called giant particles number fluctuations are one of the most relevant. We begin giving an operative, computational definition. Define a box of linear size ℓ inside your system, containing n t particles at time t. One can then measure the mean number of particles n contained in the box by taking a mean in time over different countings. In the homogeneous phase, this will be simply given by n = ρ 0 ℓ d . Together with the mean, one can also measure root mean square fluctuations ∆n = ( (n t − n ) 2 ) 1/2 . By considering boxes of different size ℓ (see Fig. 3b), one can then explore numerically the relation between the mean and its fluctuations ∆n ∼ n α (20) In equilibrium systems, away from critical point one, has generally α = 1/2 in agreement with the central limit theorem, but numerical simulations [14] show that in the entire Toner & Tu phase one has α ≈ 0.8 in both two and three spatial dimensions, as shown in Fig. 3c. Fluctuations in number density are anomalously large in the Vicsek class! This is indeed another manifestation of the slow power-law decay of correlations. A slow enough decay in space of local density 15 fluctuations δρ(r, t) correlations, corresponds indeed to an algebric divergence behavior at small frequencies q in Fourier space (where q = |q|), as opposed to ordinary equilibrium systems where σ = 0. The small frequency behavior of the stationary density structure factor S(q) gives indeed the fluctuations to mean ratio in the limit of a large particle number, Remembering that n = ℓ d ρ 0 , and that transforming back into Fourier space one has S(q → 0) ∼ 1 q σ ∼ ℓ σ (the box linear extension ℓ being a small frequency cutoff), one obtains [32] ∆n ∼ n 1/2+σ/(2d) (23) or As anticipated, the equilibrium result α = 1/2 is recovered when the structure factor is finite for q → 0, that is for σ = 0. The argument given above is slightly simplified in implicitly assuming spatial isotropy of correlation functions and of the corresponding structure factor. We indeed know that this is not the case: due to symmetry breaking spatial isotropy is broken, and correlation functions show different algebraic behaviors in the transversal and longitudinal directions. In fact, by measuring means and fluctuations in square boxes, we are taking an average over the different directions. Correspondingly, in the above argument, one should average S(q → 0) over all directions. In the Appendix we carry on this procedure in detail making use of Toner & Tu theory predictions for S(q) [27]. It yields an estimate of α = 4/5 for d = 2 and α = 23/30 for d = 3. Note that the estimates for d = 2 and d = 3 are very close to each other and in substantial agreement with current numerical data as shown in Fig. 3c.

Models with attractive/repulsive interactions and surface tension
It is finally worth noticing that the alignment rule alone is not able to maintain the cohesion a of a finite flock in open space. Fluctuations, in fact, will inevitably pull apart particles one from each other, finally disintegrating the flock. As already mentioned, in numerical simulations this problem is usually solved by introducing periodic boundary conditions, an appropriate choice when one is interested in the bulk, asymptotic properties of the Vicsek class. However, if one wants to simulate a finite group in open space, some attractive interaction should be added to introduce a surface tension and stabilize the finite flock. Attraction (together with short range repulsion) was already present in Ref. [33], where a pioneering flocking model has been proposed in the context of computer graphics, but the first study of a VM model with cohesion in a statistical physics context has been performed in Ref. [34], where Eq. (6) has been modified by adding an attraction/repulsion term. One has where e ij is the unit vector going from particle i to j r ij = |r i − r j | is the reciprocal distance and in [34] topological interaction were used. Here f is a two body force, repulsive at short range and attracting further away. For instance, one can chose f (r) = min(1, r − r e ), where r e is the equilibrium distance. By increasing the cohesion parameter β, it has been shown that the finite flock can pass from a gas phase -where the group disintegrates in open space -to a (moving) liquid one and eventually to a (moving) crystal phase. The effect of strong repulsion alone added to alignment has been discussed in [35].

Concluding remarks
In these notes, we have discussed the Vicsek model and its relative "universality class" by making use of numerical experiments and of a number of illustrative but somehow simplified arguments. A more rigorous analytical treatment of the VM asymptotic properties is given by hydrodynamic theories, but their detailed discussions clearly lies out of the scope of this lecture notes. The interested reader should consult the original work of Toner & Tu on phenomenological hydrodynamics [26,27,28,29], where an RG approach to the study of the homogeneous ordered phase is carried on, and the Boltzmann-Ginzburg-Landau approach developed in [16,17].
While the Vicsek universality class is robust to many variations, such as changes in the way the noise is implementes (as long as no long-range correlations are introduced) or the details of the local alignment interaction (but relevant changes can be introduced switching the interaction from metric to topological as discussed in Section 3.2), changes in some fundamental features are typically relevant. Modifying the nature of broken symmetry, for instance, is a typical example of such a change. For example, one may consider nematic rather than ferromagnetic alignment, without altering the polar self-propelled nature of particles (the so called self propelled rods model) [36], or consider altogether completely nematic particles (which have a preferred axis of motion but not a well defined direction) such as in active nematics [37]. These models are relevant to the modelling of elongated active particles interacting by volume exclusion forces, which typically induce an effective nematic interaction. In general, these so called Vicsek-like models constitute different universality classes, but share a very similar phase diagram structure with the Vicsek class: the phase diagram of all metric Vicsek-like models, for isntance, exhibit a phase separated regime (possibly with different symmetries/properties w.r.t. the VM) taking place at the onset of order and separating the disordered from the homogeneously ordered phase.
Other relevant changes include violation of particles number conservation, as discussed in Ref. [38], or -as previously discussed -the inclusion of momentum conservation and long-ranged hydrodynamic interactions.

A Appendix -The anomalous density fluctuations exponent
In this appendix, we compute explicitly the anomalous density fluctuations exponent making use of the results of Toner & Tu theory. The density structure factor has an anysotropic structure and it is given by (in units of the interaction distance R 0 ) [27] S(q) ∼ where q and q ⊥ are (respectively) the projection of the Fourier space vector q in the longitudinal and transversal directions w.r.t. the direction of motion. The two exponents ζ and χ are scaling exponents for which the DRG flows to a fixed point. According to a conjecture first put forward in [27], in any dimension 3/2 ≤ d ≤ 4 they are While this conjecture has never been proven rigorously, there is a reasonable numerical [39,14] evidence supporting the above scaling exponent values for d = 2 and, to a lesser extent, d = 3. In the following we will assume the above values hold. We can visualize the three different sectors which in Eq. (26) determines the scaling of the density structure factor as in Fig. 4. In particular, we are interested in the scaling behavior as one approaches q = 0 along different paths in the (q ⊥ , q ) plane (or moves towards infinity in the real axis representation (1/q ⊥ , 1/q ). It is easy to see that moving towards infinity along the line q ∼ q ⊥ ∼ q the structure factor picks up a divergence S(q) ∼ q 1−d−ζ−2χ ∼ q −2(d+1)/5 (28) This is actually the strongest possible divergence in any d < 4. Moving to infinity along the line q ∼ q ζ ⊥ , for instance, gives while chosing other paths towards q = 0 lying in the sector I, II or III Fig. 4 also produces weaker divergences or no divergences at all. This can be checked by chosing a family of paths q ∼ q ν ⊥ . The value of the exponent ν determines the chosen sector for our path, with ν > 1 corresponding to sector I, 1 > ν > ζ to sector II and ν < ζ to sector III. To summarize, the structure factor is dominated by divergences along the q ∼ q ⊥ line,