Clustering and anisotropic correlated percolation in polar flocks
aa r X i v : . [ c ond - m a t . s o f t ] M a y Clustering and anisotropic correlated percolation in polar flocks
Nikos Kyriakopoulos, Hugues Chat´e,
2, 3 and Francesco Ginelli Department of Applied Physics, Aalto University, Helsinki, Finland Service de Physique de l’Etat Condens´e, CEA, CNRS,Universit´e Paris-Saclay, CEA-Saclay, 91191 Gif-sur-Yvette, France Beijing Computational Science Research Center, Beijing 100094, China Department of Physics and Institute for Complex Systems and Mathematical Biology,King’s College, University of Aberdeen, Aberdeen AB24 3UE, United Kingdom
We study clustering and percolation phenomena in the Vicsek model, taken here in its capacityof prototypical model for dry aligning active matter. Our results show that the order-disordertransition is not related in any way to a percolation transition, contrary to some earlier claims.We study geometric percolation in each of the phases at play, but we mostly focus on the orderedToner-Tu phase, where we find that the long-range correlations of density fluctuations give rise toa novel anisotropic percolation transition.
I. INTRODUCTION
Active matter [1] typically involves moving “particles”(such as social animals [2], cells [3, 4], biofilaments dis-placed by motor proteins [5], phoretic colloids [6], etc.).Energy, either stored internally or gathered from the en-vironment, is consumed locally to produce mechanicalwork. These systems display a wide range of collectivephenomena that are not possible in equilibrium. In par-ticular, Toner and Tu have shown that flocking systemssuch as the celebrated Vicsek model [7], where constant-speed particles locally align their velocities in the pres-ence of noise, can show true long-range orientational or-der even in 2 dimensions, in a strongly fluctuating phaseendowed by generic anisotropic long-range correlations[8–12].Active matter systems are also known to often showdense clusters that dynamically form, merge, shrink, andsplit. This has been observed experimentally in situ-ations as diverse as bacteria colonies [13], acto-myosinmotility assays [5, 14, 15], animal groups [16], and activecolloidal particles [17]. A wide, powerlaw-like, distribu-tion of cluster sizes has been reported in certain casessuch as gliding myxobacteria [18]. Simple models of self-propelled rods interacting solely via steric exclusion, putforward initially in the context of bacteria, have long beenknown to exhibit similarly broad distribution of clustersizes [19–21], a situation sometimes referred to as non-equilibrium clustering. In most of these systems, theseclusters are believed to be the consequence of arrested–or micro-phase separation, with size or mass distribu-tions bounded by a finite, albeit sometimes very large,intrinsic cut-off [22, 23].Clusters also appear in flocking models such as the Vic-sek model, where they are naturally and unambiguouslydefined by making use of the finite-range of interactions.Power-law distributions of cluster sizes have also beenreported [10, 24, 25]. Because these observations weremostly made in the region of parameter space where theorder-disorder transition takes place, some authors haveconjectured that, in active systems exhibiting collective motion, this transition from disorder to ordered collectivemotion could be somehow generically related to (or evenmediated by) non-equilibrium clustering[26]. This claim,at face value, may appear rather surprising: indeed, in anoisy model such as the Vicsek model, one expects thatat large enough density, particles would always form asingle, macroscopic, spanning cluster, irrespective of thedegree of orientational order present. Conversely, at lowenough densities, one has no chance to observe a perco-lating cluster. It is thus natural to expect a percolationtransition [27] separating these two regimes.Moreover, phase-separation has been recently shownto be at play in dry aligning active matter. It actuallyprovides the best framework to understand the phase di-agram of Vicsek-style models [28, 29], which contain 3phases, with a disordered gas separated from an orderedliquid by a coexistence phase. To the best of our knowl-edge, it remains unclear whether geometric percolationand the order-disorder transition can interfere in any wayin flocking models.In this work, we come back to this issue, and studyclustering phenomena in the Vicsek model, taken here asa prototypical model for dry aligning active matter. Ourresults show that the order-disorder transition is not re-lated in any way to a percolation transition. We studygeometric percolation in each of the phases at play, butwe mostly focus on the ordered Toner-Tu phase, wherewe find that the long-range correlations of density fluc-tuations give rise to a novel anisotropic percolation tran-sition.The remainder of this paper is organized as follows: inSection II, we summarize the phase diagram of the Vic-sek model and recall some of its basic properties. Sec-tions III,IV,V describe percolation and clustering in theToner-Tu liquid phase, while we briefly examine the dis-ordered and the coexistence phase in Section VI. A dis-cussion and some conclusions can be found in Section VII.
II. THE VICSEK MODEL FOR FLOCKINGAND ITS PHASE DIAGRAM
We consider the classic version of the Vicsek model[7] with metric interactions in two spatial dimensions.Particles are defined by an off-lattice position r i and anorientation θ i ∈ [0 , π ], with i = 1 , . . . , N . The discrete-time evolution is synchronous: orientations and positionsare updated at integer time steps according to the driven-overdamped dynamics θ i ( t ) = Arg N X j =1 A tij v i ( t ) + η ξ i ( t ) (1) r i ( t + 1) = r i ( t ) + v v i ( t + 1) , (2)where v i = (cos( θ i ) , sin( θ i )) is the unit vector pointingin the direction θ i , v is the speed of particles and ξ ti is a random angle drawn uniformly in [ − π, π ] with deltacorrelations in space and time. The alignment interactionis limited to a metric range with a radius r = 1 [30], andthe symmetric and time-dependent interaction matrix A t codes for the presence of neighbors within this interactionrange: A tij = || r i ( t ) − r j ( t ) || ≤
10 if || r i ( t ) − r j ( t ) || > . (3)Effectively, the sum in Eq. (1) runs over all particles inthe unit radius disk centered around particle i ( i itselfincluded). The finite interaction radius r allows for anatural and unambiguous definition of clusters: particleswithin distance r of each other belong to the same clus-ter. At any given time t clusters are then determined asthe connected components of the graph formed by theinteraction matrix A tij .We consider square domains of linear size L with peri-odic boundary conditions, corresponding to a global den-sity ρ = N/L . In the following we fix v = 0 . ρ and the noise amplitude η , the latter playing arole akin to that of temperature in equilibrium systems.For maximum noise, η = 1, particle orientations arecompletely random and decorrelated, so that at each timestep their spatial distribution is equivalent to one drawnfrom a Poisson point process [31]. As the noise is lowered,short range correlations initially build up (both in orien-tation and position) and, as a threshold η gas is passed,the system eventually undergoes a spontaneous symme-try breaking phase transition to long-ranged (polar) or-der, easily characterized by the mean particle orienta-tions order parameter, V ( t ) = N P Ni v i ( t ).Active particles move following the orientational de-grees of freedom that they themselves carry, linking localorder and local density in a simple but highly non-trivialway. As a result, the transition between the fluctuat-ing but homogeneous disordered and ordered phases is ρ η disordered gasToner-Tupolar liquidcoexistence phase FIG. 1. Phase diagram of the Vicsek model for v = 0 . ρ, η ) plane (from Ref. [35]). The binodal line η gas ( ρ ) sepa-rating the disordered gas from the coexistence phase made oftraveling high-density high-order bands is reported in black,while the red line marks the liquid binodal η liq ( ρ ) separatingthe coexistence region from the Toner-Tu polar liquid. Thegreen line links the diamonds locating the isotropic percola-tion threshold in the disordered gas phase. The 2 blue dia-monds linked by the solid line show the asymptotic location ofthe percolation transition in the Toner-Tu liquid phase deter-mined through finite size scaling (see Sec. IV). The η = 0 . not direct, like originally thought in analogy with mag-netic systems such as the XY model, but mediated bya coexistence phase where high-density ordered bandsmove in a low-density disordered background [10, 32–34]. Within the coexistence phase, increasing the globaldensity and/or the system size, the number of travelingbands increases linearly while the residual vapor densitybetween them remains constant [35]. We are thus in thepresence of a phase separation scenario: the disorderedgas (DG) is separated from the ordered Toner-Tu polarliquid (PL) by a coexistence region with a quantized liq-uid fraction (the traveling bands are microphases). Thecorresponding asymptotic phase diagram, following thenumerical results of [35], is reported in Fig. 1. Onehas thus two transitions, not one, marked by the twobinodal lines separating these different phases. Theyare non-decreasing functions of density, η gas = η gas ( ρ ), η liq = η liq ( ρ ), and in the limit of small densities onehas η gas ∼ √ ρ [36]. An inaccessible critical point ispushed towards infinite density [29]. The two transitionsare continuous (in the infinite-size limit) but not critical.At finite size, they appear discontinuous because of thelarge number of particle involved in nucleating a travelingband. III. CLUSTERING AND ANISOTROPICPERCOLATION IN THE TONER-TU LIQUIDPHASE
In the following 3 sections, we focus our attention onthe polarly-ordered Toner-Tu liquid phase. We initiallyfix the noise amplitude to η = 0 .
2, and study the clus-tering behavior for different particle densities ρ > . η liq in Fig. 1. Three typicalsnapshots, obtained in the stationary regime for increas-ing global densities ( ρ = 1 , ρ = 1 . ρ = 1 .
9) areshown in Fig. 2 for a system of linear size L = 256.At the lowest density value ρ = 1, the largest clus-ters are clearly smaller than system size. The largestof them contains less than 10% of the total number ofparticles. The transversal extension (with respect to thecurrent global order direction) of these largest clustersis much larger than the longitudinal one. Increasing thedensity, clusters remain clearly anisotropic and some ofthem are spanning across the system along the transver-sal direction. In the central panel of Fig. 2, one sees asingle spanning cluster, comprising less than one thirdof the total number of particles. It is only at densities ρ & . C ( r ) = | r ⊥ | χ f ( r k / | r ⊥ | ζ ) , (4)where k and ⊥ indices respectively refer to directionslongitudinal and transverse to the mean motion of theflock and the exponents χ and ζ as well as the function f are universal. A notable consequence of this fact isthat, in two spatial dimensions, the particles’ displace-ment transversal to the mean velocity is superdiffusive[10, 37], while it is simply diffusive in the longitudinaldirection (once substracted the mean motion).We now characterize the percolation transition and itsanisotropy from a more quantitative point of view. Indi-vidual clusters (labeled by k ) can be quantified by theirmass s k , that is, the number of particles in the clus-ter, and by their linear extension ℓ k , which we define astwice the in-cluster maximum distance between a clus-ter particle and the cluster center of mass[38]. The in-stantaneous maxima of these quantity are respectively s M = max k ( s k ) and ℓ M = max k ( ℓ k ) where the clusterindex k runs over all clusters of a given configuration. Two order parameters are routinely employed in theliterature [39] about isotropic percolation problems, the(normalized) mean largest cluster size n and the meancluster maximum linear extension d , where the averageis taken over many different realizations (e.g. samplinga long trajectory in the stationary state at regular timeintervals). The definitions of n and d and the associatedstandard deviations σ n and σ d read: n ≡ h s M i N , σ n ≡ p h ( s M − h s M i ) i N (5) d ≡ h ℓ M i√ L , σ d ≡ p h ( ℓ M − h ℓ M i ) i√ L . (6)In our anisotropic situation, n and d are expected to be-have differently as the density is increased, with d risingearlier to order 1 values and n following later. This isindeed observed in Fig. 3a where the two indicators arecompared for a system of size L = 256 and η = 0 . σ n and σ d peak at two differentdensity values (Fig. 3b).We also investigated directly the spanning probabil-ity S , i.e. the probability that a spanning cluster doesappear. While, in the thermodynamic limit, S ( ρ ) is astep function with the jump exactly located at the phasetransition, in finite systems S ( ρ ) is smoothed around the(finite size) percolation point [27]. To take into accountanisotropy, we consider both a transversal and a longi-tudinal spanning probability, S ⊥ and S k , defined respec-tively as the probability that a cluster wraps around the L × L torus in the transversal or longitudinal (w.r.t. theorder parameter V ) directions [40]. (We discuss the ac-curacy of these measures in finite systems – where fluc-tuations lead to the diffusion of the instantaneous meanorientation of motion V ( t ) – in the next section.) Forthe moment, we simply note that the transversal span-ning probability S ⊥ rises from zero towards one earlierthan the longitudinal probability S k , as shown in Fig. 3c.Before concluding this section, we note that a reliablenumerical evaluation of the above configuration averages– as the one presented in Fig. 3 – is considerably moredifficult to obtain than in standard percolation problems,where the probability distribution of the particles posi-tion is exactly known and systems configurations can begenerated from it. In our case, on the contrary, one has togenerate sufficiently uncorrelated configurations from thedynamics. This requires first to evolve the system fromsome initial condition into the stationary state (whichfor large systems may require a considerable number oftimesteps). Then, in order to obtain configuration aver-ages, one has to take averages over timescales T muchlarger than the typical autocorrelation time. Considerfor instance the timeseries of s M and ℓ M discussed above(an example of which for ρ = 1 . L = 256 and η = 0 . τ [41] in the low densityregime are of the order of 10 timesteps. Note howeverthat near the percolation transition τ drops suddenly by FIG. 2. Typical instantaneous snapshots in the Toner-Tu ordered phase at different densities ((a) ρ = 1, (b) ρ = 1 . ρ = 1 . v = 0 . η = 0 .
2. Colors correspond to connected clusters of particles, with the largest cluster in red. Notethat due to the large number of different clusters (in the order of thousands in all three panels), each color is used for severaldistinct clusters, hopefully sufficiently apart from each other to avoid confusion. The thick black arrow marks the instantaneousdirection of global order (i.e. the order parameter orientation) almost two orders of magnitude. This could seem coun-terintuitive, as phase transitions are typically associatedwith a slowing down of the dynamics. However, one hasto realize that near the percolation transition one has awide distribution of competing clusters with sizes closeto the spanning threshold, so that relatively small con-figuration changes may promote a different cluster to thelargest cluster status (either in total mass or linear ex-tension) thus resulting in a dramatic drop in the auto-correlation time for s M and ℓ M .Once the autocorrelation time has been estimated, theaccuracy of the empirical averages can be evaluated bythe standard error σ/ p T /τ , where σ is one standarddeviation and T /τ the number of independent configura-tions.
IV. FINITE SIZE SCALING ANALYSIS OFPERCOLATION IN THE TONER-TU LIQUIDPHASE
One of the best ways to numerically investigate crit-ical phase transitions is to perform a finite-size-scaling(FSS) study, measuring the lowest moments of suitableorder parameters as the system size is systematically in-creased. This is a classical approach in statistical physics,routinely applied to study both equilibrium and out-of-equilibrium critical phase transitions [42], and it has beenalready applied to the study of the percolation transition,for instance bond percolation on square lattices [27, 43].The main difficulty, generally, is to be sure to probesystem sizes large-enough so that one is in the scalingregime.
A. Percolation in the longitudinal direction
We first concentrate on the longitudinal percolationtransition, i.e. the point at which the spanning clusterbecomes two-dimensional and starts to span also in thebroken symmetry direction.The mean largest cluster size n measures the probabil-ity n/N that an arbitrary particle belongs to the largestcluster. In percolation theory, it is known to follow thefinite size scaling relation [43] n = L − β/ν f (( ρ − ρ ∞ c ) L /ν ) , (7)where f is a scaling function and β and ν two universalcritical exponents. At ρ = ρ ∞ c , the asymptotic criticalpoint, f (0) = const. and one obtains the power-law be-havior n ∼ L − β/ν . In two dimensional standard perco-lation, one has ν p = 4 / β p = 5 /
36 [44] (and thus β p /ν p = 5 / ρ and the sys-tem size between L = 64 and L = 1024, we find (seeFig. 4a) that for ρ ∞ c ≈ .
95 the mean largest clustersize indeed follows a power-law decay with an exponentcompatible (our best fit being β/ν = 0 . β p /ν p = 5 / χ n ≡ L σ n , χ Mn = L − γ/ν , (8)with – for d = 2 standard percolation – γ p = 43 /
18 [44] sothat γ p /ν p = 43 /
24. For each system size L we estimatethe peak susceptibility χ M by a quadratic fit of the peakregion of σ n ( ρ ). In Fig. 4b we show that once again ournumerical estimates are in very good agreement with thestandard percolation exponent. Indeed, our best fit of ρ n (a) ρ σ n (b) 1 2 ρ t ρ τ (d) d σ d S ⊥ S // l M s M FIG. 3. Anisotropic percolation transition in the Toner-Tuordered phase. (a) Normalized mean largest cluster size n (red dots) and normalized mean cluster maximum linear ex-tension d (black squares) as a function of the global density ρ .(b) Corresponding standard deviations σ n (red dots) and σ d (black squares). The blue dashed lines show quadratic fits ofthe peak regions (see text). (c) Transversal ( S ⊥ ) and longitu-dinal ( S k ) spanning probabilities as a function of total density.The dashed lines show a fit based on the error function (seetext). (d) Autocorrelation time τ of the timeseries of maximalcluster size s M (red dots) and maximal linear extension ℓ M (black squares) as a function of global density. Inset: typicalexcerpts from these timeseries for ρ = 1 .
5. Other parameters: L = 256, η = 0 .
2, and v = 0 .
5. Configurations averageshave been computed sampling every 100 timesteps a T = 10 timeseries in the stationary regime. The standard error (seetext) of the data shown in panels (a)-(c) is equal or smallerthan the symbol size. γ k /ν k = 1 . / ν that determines finite size corrections to thecritical point, ∆ ρ ≡ ρ ∞ c − ρ c ( L ) ∼ L − /ν . (9)Here we adopt and compare two different estimates forthe finite size critical density ρ c ( L ). We first estimate itas the location ρ M of the maximum of the largest clus-ter size standard deviation σ n ( ρ ) (once again, evaluatedthrough a quadratic fit of the peak region). Our results,illustrated in Fig. 4c (green squares) essentially confirmour previous estimate for the asymptotic critical point, ρ ∞ c = 1 . /ν k ≈ .
5. It has to benoted that this estimate is based on a second moment (the standard deviation), so that its reliability could bequestioned.A second, and perhaps more accurate estimate of thefinite size critical density can be obtained measuring thedensity value by which the finite size spanning probabil-ity crosses 1 /
2. We are here interested in the longitudi-nal spanning probability S k . Measuring it in relativelysmall systems, where the mean orientation V ( t ) stronglydiffuses in its anngular component, can be however a dif-ficult task. The central limit theorem implies that themean orientation should diffuse with an angular diffusionconstant proportional to η /N . For small enough systemsizes, thus, the mean orientation can change faster thanthe time needed by clusters to re-align transversally w.r.t. V ( t ). Therefore, in our FSS analysis we find preferableto consider, instead of the transversal and longitudinalspanning probabilities the one and two dimensional span-ning probabilities S and S . The former is the proba-bility that a cluster spanning along at least one spatialdirection (i.e. to join two opposite sides of the system)does exist. The latter probability, on the other hand,requires the spanning cluster to wrap along both spatialdirections, that is to join all four sides of our system.Numerical simulations show that – at least in the pa-rameter range we are interested in – for L '
256 we haveto a good accuracy S ⊥ ≈ S and S k ≈ S . For smallersystem sizes, however, we have S < S ⊥ < S k < S .In the following, we estimate the finite size critical den-sity as the density value by which an error function basedfit [45] of the finite size spanning probability S ( ρ ) crosses1 / ρ ∞ c = 1 . /ν k ≈ . ν k are clearly different from standard percolation in d = 2, ν p = 4 /
3. Combining our two different approaches weget ν k = 2 . ν k ≈ . /ν k = 0 . ν k = 2 . ν k , different fromthe one of standard percolation, is indeed confirmed byattempting a data collapse of the mean largest cluster size n according to the scaling relation (7). Our data clearlyrule out the value ν p = 4 /
3, and we obtain a satisfactorycollapse with β k /ν k = β p /ν p = 5 /
48 and 1 /ν k in therange 0 . ∼ . ρ ∞ c =1 . r randomly distributed accord-ing to a Poisson Point Process (PPP) corresponds to acritical area A c = πr ρ P P Pc = 1 . r = 1 /
2, we have ρ P P Pc = 1 . L n ρ=1.7ρ=1.85ρ=1.9ρ=1.92ρ=1.95ρ=2ρ=2.1 (a) L ∆ρ (c) 10 L χ n (b) ρ n -10 0 20 ∆ρ L ν n L β/ν (d) M FIG. 4. Finite size scaling analysis of the percolation tran-sition in the Toner-Tu ordered phase. (a) mean largest clus-ter size n vs. system size L for different densities (see leg-ends). The dashed red line marks the standard percolationcritical exponent ratio β p /ν p = 5 /
48. (b) susceptibility peakvalue (black dots) χ Mn vs. system size L . The dashed redline marks the standard percolation critical exponent ratio γ p /ν p = 43 /
24. (c) Critical point location finite size correc-tions ∆ ρ = ρ ∞ c − ρ c ( L ) evaluated either from the midpointof the spanning probability S (red full dots) or from thepeak location of the standard deviation σ n ( ρ ) (green emptysquares). Here we have used ρ ∞ c = 1 .
96. The dashed red linemarks a power law decay with an exponent − .
4, while thedashed green line falls off as L − . . (d) Data collapse of n according to the scaling relation (7) with ρ ∞ c = 1 .
96 and ex-ponents β = β p = 5 /
36, 1 /ν k = 0 . L = 64 and L = 1024. Inset: Non collapsed curves.From top to bottom: L = 64, L = 96, L = 128, L = 192, L = 256, L = 384, L = 512, L = 768 and L = 1024. Otherparameters are η = 0 . v = 0 .
5. As in Fig. 3, averageshave been computed sampling a 10 timesteps long timeseriesevery 100 timesteps. The standard error (see text) of the datashown in panels (a)-(c) is equal or smaller than symbol size. B. Harris criterion for percolation in correlateddensity fields
While a shift in the critical percolation point is notsurprising in the presence of activity, and indeed has beenobserved before in disordered active matter systems [39],the significant difference between our estimate for thecritical exponent ν and the standard percolation value ν p deserves a few more comments.It is indeed known that long-range correlations in theparticle density can change the value of the critical ex-ponent ν . In the percolation literature, this is known asthe Harris criterion [47, 48]: in the presence of sufficiently long-ranged density correlations C ρ ( r ) ∼ r − α with α < ν p , (10)finite size corrections are indeed stronger and the expo-nent ν takes larger values: ν = ν H = 2 α . (11)On the other hand, for correlations decaying faster, α > /ν p , correlations are not relevant and usual finite sizecorrections apply, ν = ν p .Applying the Harris criterion to our results suggeststhat the density field correlation should decay with apower-law with an exponent α = 2 /ν k in the range 0 . ∼
1. Using only the estimate derived from the spanningprobability distributions, we would have α = 0 . S ) could be well captured by measures of densitycorrelations averaged over all spatial directions. In thefollowing, therefore, we focus on isotropic correlations ofdensity fluctuations C ρ ( r, L ) = hh δρ ( x + r , t ) δρ ( x , t ) i S i t (12)where r = | r | and δρ ( x , t ) ≡ ˜ ρ ( x ) − ρ are the local den-sity fluctuations of a suitably coarse-grained density field˜ ρ ( x ), and h·i S indicates an average over the spatial co-ordinate x and the orientations of the displacement r .Isotropic correlations are then further averaged in time,with h·i t indicating an average over stationary state con-figurations. A more detailed analysis of anisotropic cor-relations will be reported in [49].The scaling of correlations is expected to be the samein the entire Toner-Tu phase. Here we focus on a pointclose to the percolation threshold, η = 0 . ρ = 1 . R d x δρ ( x , t ), and this implies that the cor-relation function C ( r, L ) should have at least one zero (asthere are surely anti-correlated regions). The smallestvalue of r for which correlations vanish can be taken as a L C( r ) (a)0 500 1000 L ξ
10 100 ξ -2 -1 h
10 100(c) 10 -2 r / ξ C( r ) ξ (d)0.01 1 r / ξ C( r ) (b) FIG. 5. (a) Isotropic density fluctuations correlationfunction in the Toner-Tu ordered phase ( v = 0 . η = 0 . ρ = 1 .
9) and increasing system sizes, L =64 , , , , , , , , ξ as a function of system size L .The dashed red line marks our best linear fit. (b) Same as (a)but after rescaling of the space variable. System size increasesalong the green arrow. (c) Finite size scaling of the (negative)slope h (see text). The dashed black line marks a power-lawdecay with exponent − α = − .
8. (d) Data collapse accordingto Eq. (16). Correlation functions have been averaged over10 different spatial configurations, sampled from the station-ary state dynamics every 10 timesteps. Standard errors areof the size of the symbols or smaller. measure of the correlation length ξ , that is C ( ξ, L ) = 0.Moreover, in systems where a continuous symmetry isspontaneously broken correlations are known to be scalefree, i.e. ξ ∼ L (see inset of Fig. 5a) and C ( r ) ∼ r − α inthe thermodynamic L → ∞ limit. Finite size correlationsare thus taken into account by [11] C ρ ( r, L ) = r − α g (cid:18) rξ (cid:19) , (13)where the scaling function obeys g ( u ) = 0 for u = 1 and g ( u ) → const. for u → α can be deter-mined by finite size analysis. Let us choose the rescaling y = r/ξ . From Eq. (13) we have (see Fig. 5b) C ρ ( y, L ) = y − α ξ − α g ( y ) . (14)From Eq. (14) it follows that a finite size analysis ofthe (negative) slope h of the rescaled correlation functionevaluated in y = 1 can be used to estimate the correlationexponent, h = − ddy C ρ ( y, L ) | y =1 = ξ − α | g ′ (1) | ∼ ξ − α . (15) Our best numerical estimates, reported in Fig. 5c, areindeed compatible with the correlation value suggestedby the Harris criterion, α ≈ . C Rρ ( y ) ≡ ξ − α C ρ (cid:18) rξ , L (cid:19) . (16)Our brief analysis of density correlations shows thatthe anomalous finite size corrections exponent ν wehave measured for our percolation transition (especiallythrough the more reliable spanning probabilities mea-sures) is fully compatible with the one expected by theHarris criterion for correlated percolation.We finally note that also the critical exponents β and γ may be modified by sufficiently strong correlations. How-ever, it has been however verified numerically [44] thatthe hyperscaling relation of standard percolation ν p d s = 2 β p + γ p (17)(with d s being the spatial dimension) is still verified bythe correlated Harris exponents ν H d s = 2 β H + γ H . (18)Interestigly, we find that this latter hyperscaling relationis also verified by our data: as we have seen, our two di-mensional estimates for the ratios β/ν and γ/ν are bothcompatible with the values expected by standard corre-lation theory, so that2 βν + γν = 2 . . (19) C. Percolation in the transversal direction
We finally discuss the percolation transition takingplace in the transversal direction. Repeating the proce-dure outlined for the longitudinal percolation transition,we evaluate the finite size transversal percolation thresh-old density ρ t ( L ) from the behavior of the one dimen-sional spanning probability S (see Fig. 6b). Our bestestimates are reported in the top panel of Fig. 6c (blacksquares) together with the ones for the finite-size lon-gitudinal percolation critical density ρ c ( L ) (red circles).Quite interestingly, our numerical results indicate thattheir difference ∆ c t ( L ) ≡ ρ c ( L ) − ρ t ( L ) (blue triangles inthe lower half of Fig. 6c) seems indeed to vanish in thelimit of large L , suggesting that in the thermodynamiclimit ρ t ( L ) → ρ ∞ c = 1 . ρ t ( L ) ≡ ρ t ( L ) − ρ ∞ c ∼ L − /ν k . (20) -4 0 4 ∆ρ L ν d (d)1 1.5 2 ρ S (b) 100 1000 L ∆ ct ρ S (a) 1.5 ρ t (c) ρ d σ d ρ c ∆ρ t ~ FIG. 6. (a)-(b) Two- and one-dimensional spanning probabil-ities S and S as a function of global density ρ for differentsystem sizes ( L = 64 , , , , / L . The horizontal green line marks our bestestimate for the asymptotic percolation point ρ ∞ c = 1 .
96. (c)Bottom panel: Transversal percolation point finite-size cor-rections ∆ ρ t = ρ ∞ c − ρ c ( L ) (black squares), longitudinal totransversal finite size difference ∆ c t (blue triangles, see text)and maximum variance of the largest cluster linear extension˜ σ d as a fuction of system size L in a double logarithmic scale.The dashed black line marks a power-law decay with an expo-nent 0 .
3, while the red one corresponds to constant behavior.(d) Data collapse of the mean maximum cluster extension d according to the scaling relation (21) with ρ ∞ c = 1 .
96 and1 /ν ⊥ = 0 . L = 64 and L = 1024 (color coded as in panels (a) and (b)). In the inset:Non collapsed curves for the mean largest cluster size L vs.density ρ . Along the cyan arrow: L = 64 , , , , This can be deduced from the lower panel of Fig. 6c),where ∆ ρ t ( L ) (black squares) exhibits a power-law de-cay compatible with an exponent 1 /ν ⊥ ≈ .
3, suggesting ν ⊥ ≈ . d seems to show rather anomalous scaling properties. Asit can be readily deduced from the inset of Fig. 6d, itsfinite size curves do cross near the critical density. Thisimplies that its scaling relation should take the form d = f ⊥ (( ρ − ρ ∞ c ) L /ν ⊥ ) (21)with f ⊥ a transversal scaling function. As we show inFig. 6d, one can indeed make use of this scaling rela- tion to achieve a satisfactory collapse of the d ( ρ ) curvesby only rescaling them along the abscissas. Comparisonwith the general scaling form (7) thus implies the rathersingular β ⊥ = 0.Finally we discuss the γ exponent, associated to themaximum linear extension susceptibility χ d ≡ L σ d , (22)whose peak value is expected to scale as χ Md ∼ L γ ⊥ /ν ⊥ . (23)Assuming that an hyperscaling relation analogous to (17)still holds between the transversal exponents, in two spa-tial dimensions we would get γ ⊥ /ν ⊥ ≈
2. By virtue ofEq. (22), this in turn implies that also the peak varianceof the largest cluster linear extension should not scalewith system size,˜ σ d ( L ) ≡ max ρ σ d ( ρ, L ) ∼ const. (24)This is indeed verified by our numerical data (see thebottom panel of Fig. 6c, where the maximum has beenevaluated by a quadratic fit of the peak region). We con-clude that the cluster maximum linear extension d showsno finite size scaling, apart from finite size corrections inits density dependence. D. Anisotropic percolation exponents
Our estimates for the Toner-Tu phase percolation ex-ponents in two spatial dimensions, as measured from sim-ulations of the Vicsek model, are summarized in Table I
TABLE I. Longitudinal and transversal percolation exponentscompared with the ones of standard percolation theory. Lon-gitudinal exponents are particularly difficult to evaluate dueto strong finite size effects and ours are rough estimates. Forthis reason we are not confident in providing precise uncer-tainty estimates. 1 /ν β/ν γ/ν
Standard percolation, d s = 2 3/4 5/48 43/24Longitudinal percolation 0.44(6) 0.108(5) 1.83(5)Transversal percolation 0.3 0 2 In the Toner-Tu theory, anisotropy is controlled by theexponent ξ [8], so that one should expect ν ⊥ = ν k /ξ or ξ = ν k ν ⊥ ≈ . ∼ . . (25)In two spatial dimensions, based on some renormalizationgroup conjectures, Toner and Tu suggested [37, 50] that ξ = 3 /
5, a value which coincides with the lower boundof our FSS measure. However, it should be noted thatthe more reliable spanning estimates rather support theupper bound of our estimates ξ ≈ .
75, thus suggestinga less severe anisotropy. More precise measures will berequired to shed light on this issue [49].Before concluding this section, we briefly comment onthe behavior of the percolation threshold as a functionof the Vicsek noise amplitude η . While a careful deter-mination of the full percolation line is beyond the scopeof this work, preliminary simulations indicates that for η = 0 . ρ ∞ c = 2 . V. CLUSTER SIZE DISTRIBUTION IN THETONER-TU LIQUID PHASE
We proceed to discuss cluster size distributions, awidely used quantity both in percolation theory and inthe literature on non-equilibrium clustering in active sys-tems.The cluster size distribution (CSD) is one of the sim-plest objects to be computed numerically in percolationtheory. One should notice though, that CSD correspondsto two different meanings in the literature. In a first ap-proach, the CSD P ( s ) measures the (properly normal-ized) number of clusters with size s one finds in givenconfigurations. This corresponds in practice to the prob-ability to find a cluster of size s when we pick at randomone of the many clusters we identify in our dynamics.Other authors, however, prefer to work with the proba-bility Q ( s ) that a particle picked at random belongs to acluster of size s . Obviously the two measures are related, Q ( s ) = sP ( s ), so that the choice between the two abovedefinitions is equivalent. In the following we consider P ( s ). We measure it by sampling a large number (typ-ically 10 ) of different steady-state configurations of ourdynamics, obtained from a single run (after a dynamicaltransient has been discarded), with 100 time units sep-arating consecutive configurations. In the following, wemay also find convenient to further rescale the clustersize s by the total number of particles N , so that we dealwith a normalized cluster size variable s/N ≤ η = 0 .
2. Our results, reportedin Fig. 7a, suggest that cluster size in the desity interval ρ ∈ [1 . , .
2] follows a power-like behavior over a widerange of scales (about four decades for the size consid-ered). This is in agreement with previous studies [10].Considering density values further out from the perco-lation point, see for instance Fig. 7b, one observes clearexponential cut-offs from power-law behavior. Note thatabove the percolation density, ρ > ρ ∞ c , where a single gi-ant connected cluster typically appears, the CSD showsan exponential cut-off, but also, beyond that, a finiteprobability of observing clusters of size s ≈ N .The proper way to discriminate a true power-law be-havior from an approximate one is, once again, finite- -6 -4 -2 s/N -10 -5 P( s ) (a) 10 -6 -4 -2 s/N -10 -5 P( s ) (b) FIG. 7. (a) Cluster size distribution P(s) vs. the rescaledcluster size s/N at different densities in the Toner-Tu orderedphase (from top to bottom: ρ = 1 . ρ = 1 . ρ = 1 . ρ = 2 . ρ = 2 . − .
9. (b) P ( s ) at the percolation threshold ( ρ = 1 .
95, fullred circles), compared with off-critical values ρ = 1 (blacksquares) and ρ = 4 (blue diamonds) for L = 1024. Themagenta dashed line marks a power-law decay with the Fisherexponent τ F = 187 / ≃ . timestepstrajectory every 100 timesteps. Other parameters: L = 1024, η = 0 . v = 0 . size analysis. We considered systems of different sizesbetween L = 64 and L = 1024. Off-critical CSDs, as theones shown in Fig. 8a, exhibit an exponential cut-off atsize Σ. While Σ may initially grow with system size, finitesize analysis of its estimated value [51] shows saturationeffects towards an asymptotic value Σ ∞ ( ρ ). As shownin Fig. 8c, this saturation seems to occur for all densi-tiy values different from the critical percolation density,with Σ ∞ ( ρ ) increasing as the percolation threshold is ap-proached from both sides. This implies that the powerlaws reported in Fig. 7a are not asymptotic. It is onlyat the anisotropic percolation point ρ ≈ ρ ∞ c = 1 . ∞ ( ρ ) diverges,and a truly asymptotic critical CSD appears. CSDs atthe percolation threshold at different system sizes are re-ported in Fig. 8b. They show a large size peak corre-sponding to the typical size of the percolating cluster,which is clearly scaling with the system size N , as it canalso be appreciated from Fig. 8b, where we have used thelocation of this peak to estimate the critical point typicalcluster size Σ (full red dots).We finally estimate the power law decay exponent atthe percolation point. At the critical point of standardpercolation, the cluster size distribution power law be-havior is controlled by the so-called Fisher exponent, τ F = 2 d s − β/νd s − β/ν , (26)which only depends on the spatial dimension d s and onthe critical exponents ratio β/ν [27]. In two spatial di-mensions we get τ F = ≈ .
05. We have seen that inthe longitudinal percolation transition, the scaling of thelargest cluster size n is still controlled by the standard0 s -10 -5 P ( s ) (a) s -10 -5 P ( s ) (b) N/ ρ Σ (c) LL FIG. 8. (a) Finite size variation of P ( s ) for ρ = 1 (follow-ing the cyan arrow L = 64 , , , , ρ = 1 .
95, where no size-dependent cut-off is present. The magenta dashed line marksa power-law decay with the Fisher exponent τ F = 187 / λ (see text) as a function of N/ρ for different density values: ρ = 1 (black circles), ρ = 1 . ρ = 1 . ρ = 1 .
95 (full redcircles), ρ = 2 . ρ = 2 . ∼ N . Alldistributions are log-binned, and have been computed sam-pling a 10 timesteps trajectory every 100 timesteps. Otherparameters: η = 0 . v = 0 .
5. All panels are in a doublelogarithmic scale. percolation exponent ratio β p /ν p , so that we also expectour cluster size distribution near the critical percolationdensity ρ ∞ c ≈ .
96 to behave as in standard percolation,that is P ( s ) ∼ s − τ F . (27)This is indeed verified by our data. For ρ = 1 . τ = 2 . ρ = ρ ∞ c , which have been foundtypically in the range [1 . ,
2] [10, 24].Altogether, our results show that, while truly criticalCSDs only appear at the percolation point, the Toner-Tuordered phase nevertheless displays an extended “quasi-critical” region, where cluster size distributions follow apower-law over several orders of magnitudes and for awide range of densities. This approximate critical regimehas also been reported in previous works [10, 24, 26]and – as we have discussed in the introduction – hasled some authors to speculate that the onset of collectivemotion should be accompanied by a percolation transi-tion. The analysis of the anisotropic percolation transi-tion carried on in the previous chapter however, clarifiesthat the Toner-Tu phase of finite size systems is charac-terized by a “double” percolation transition, with giant -6 -4 -2 s/N -12 -6 P( s ) (a) 0 0.4 0.8 η n (b)0 0.4 0.8 η n (c) 0.3 0.6 η n (d) ddd DG PL DGDG PLPL
FIG. 9. (a) Cluster size distribution P(s) vs. rescaled clustersize s/N for η = 0 . L = 1024); from top to bottom: ρ = 1 . ρ = 1 .
51 (the percolation point, full red circles),and ρ = 1 . τ F = 187 / n (green squares)and maximum cluster extension, d (orange circles) vs. noiseamplitude η across all 3 phases. Parameters: v = 0 . L = 1024, ρ = 0 .
5, (c) L = 1024, ρ = 1 (d) L = 512, ρ = 4.The black vertical dashed lines mark the finite-size onset oforder [52], η = η gas ( L ), while the red ones, at η = η liq separatephase coexistence (between the two vertical lines) from theToner-Tu polar liquid (PL). clusters first percolating transversally w.r.t. the mean di-rection of motion and, at higher densities, also spanningin the longitudinal direction. We conjecture that this ex-tended region of scaling is related to the two separatefinite size transitions at two clearly different densities.Note also that far away from this “extended region”, thecluster size distributions are clearly not scale free, see forinstance the case ρ = 4, η = 0 . VI. PERCOLATION AND CLUSTERING INTHE DISORDERED AND COEXISTENCEPHASES
A percolation transition is of course also found in thedisordered gas phase. It is a simple isotropic one withstandard exponents. Its transition line is reported inFig. 1, and for maximal noise culminates at the wellknown critical point for a Poisson Point Process, ρ P P Pc =1 . η < η = 0 . ρ = 0 . ρ = 1) and above ( ρ = 4) the percolation transition linesof both the disordered gas and the Toner-Tu phases, andvaried the noise amplitude as shown in Fig. 1, in order tocut across both binodals. We computed both the largestcluster size, n and the maximum cluster extension, d . Forlow densities, data shows that in the gas and Toner-Tuphases, clusters are small and do not reach a macroscopic,system spanning state (Fig. 9b-c). However, in the co-existence region [52], high and low local density patchesappear (signaling the presence of ordered liquid bandstravelling in a disordered gas) [35], and system spanningclusters suddenly appear. On the other hand, at largedensities (Fig. 9d), in both the gas and Toner-Tu phases,one has typically a single cluster encompassing almost allparticles, with n, d ≈
1. The appearance of lower densitydisordered patches, on the other hand, induces a drop inthe maximum cluster size in the coexistence region. Ithas been shown that due to these effects, also the clustersize distribution built by averaging over both phases inthe coexistence region show apparent power laws, albeitwith a decay exponent larger than the Fisher one [10].
VII. DISCUSSION AND CONCLUSIONS
Our numerical results show that nonequilibrium clus-tering effects in the two dimensional Vicsek model areessentially controlled by an underlying percolation point,and are therefore mainly geometrical in nature. Clusterdynamics and cluster size distributions behave differentlynot only in the different phases, but also within phases,as one always expects to cross a percolation transitionwhen the density is sufficiently large. Moreover, cross-ing one of the binodal lines delimiting the coexistencephase separating the disordered gas form the Toner-Tuordered liquid, sudden changes are typically observed inthe cluster dynamics and corresponding cluster size dis-tributions. These transitions, however, are dictated bythe overall phase-separation scenario of the phase dia-gram, and not vice-versa.In the disordered gas phase, a standard percolationtransition is observed, akin to that observed at maximalnoise (i.e. in a system fully equivalent to a Poisson point process), with standard percolation exponents [27] but aslight shift in the critical percolation density due to shortrange correlations.In the Toner-Tu symmetry-broken phase, on the otherhand, we have identified a novel anisotropic percolationtransition with clusters first spanning the transversal di-rection (w.r.t. the mean direction of motion) and onlylater, at higher densities, spanning also along the longi-tudinal direction. A careful finite size analysis revealedthat these two distinct percolation thresholds seem toconverge to the same density value in the thermodynamiclimit, albeit with two different correlation exponents ν ⊥ and ν k , which are are also clearly different from the well-known value of the standard percolation correlation ex-ponent ν p in two spatial dimensions.We have argued that the difference in the correla-tion exponents can be attributed to the long-range cor-relations which characterize density fluctuations in theToner-Tu phase. In particular, making use of the Har-ris criterion [47] for correlated percolation, we have beenable to link the value of the longitudinal correlation ex-ponent (the one controlling the onset of a cluster ofmacroscopic mass spanning in both directions) with theisotropic (i.e. averaged over all directions) density fluc-tuation correlations.The hyperscaling relation of standard percolationseems to hold also in the correlated Toner-Tu phase, withthe key exponents controlling the cluster size distribution(the Fisher exponent τ F ) and the first two momenta ofthe maximum cluster size ( β/ν and γ/ν ) compatible withtheir values from standard percolation theory.In general, it is only at the percolation point that thecluster size distribution is truly scale free ( P ( s ) ∼ s − τ F ).However, cluster size distributions resembling power-lawsover a wide range of scales occur for a finite range ofdensities in the Toner-Tu phase, presumably because ofthe “double-threshold” mechanism of anisotropic perco-lation. Only a careful finite size analysis can show thatthese power laws are not asymptotic but bounded by asize-independent cut-off.Here we have closely analyzed clustering and percola-tion in the classical Vicsek model for flocking, but weexpect our main conclusions to be generic and to hold inthe more general context of dry aligning active matter. ACKNOWLEDGMENTS
We have benefited from discussions with F. Perez-Reche. NK and FG acknowledge support from the MarieCurie Career Integration Grant (CIG) PCIG13-GA-2013-618399. FG also acknowledges support from EPSRCFirst Grant EP/K018450/1.2 [1] S. Ramaswamy, Ann. Rev. Cond. Matt. Phys. et al. Proc. Natl. Acad. Sci. USA ,1232(2008).[3] N. Sep´ulveda et al. , PLoS Comput Biol e1002944(2014).[4] X. Trepat et al. , Nature Phys.
426 (2209).[5] V. Schaller et al. , Nature , 73 (2010).[6] J. Palacci et al. , Science.
936 (2013).[7] T. Vicsek et al. , Phys. Rev. Lett. , 1226 (1995).[8] J. Toner and Y. Tu, Phys. Rev. Lett. , 4326 (1995).[9] S. Ramaswamy, R.A. Simha, and J. Toner, Europhys.Lett. , 196 (2003).[10] H. Chat´e, et al. , F. Ginelli, G. Gregoire, F. Raynaud,Phys. Rev. E , 046113 (2008).[11] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. San-tagati, F. Stefanini, and M. Viale, Proc. Natl. Acad. Sci.U.S.A. , 11865 (2010).[12] F. Giavazzi et al. J. Phys. D: Appl. Phys. , 13626 (2010).[14] S. K¨oler, V. Schaller, and A. Bausch, Nature Mater. ,462 (2011).[15] L. Huber, R. Suzuki, T. Kr¨uger, E. Frey, and A. R.Bausch, Science , 255 (2018).[16] F. Ginelli et al. , Proc. Natl. Acad. Sci. USA ,12729(2015).[17] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert,and L. Bocquet, Phys. Rev. Lett. et al. , Phys. Rev. Lett. ,030904(R) (2006).[20] Y. Yang, V. Marceau, and G. Gompper, Physical ReviewE , 031904 (2010).[21] M. Abkenar, K. Marx, T. Auth, and G. Gompper, Phys.Rev. E , 062314 (2013).[22] E. Tjhung, C. Nardini, and M. E. Cates, Phys. Rev. X , 031080 (2018)[23] X. Shi and H. Chat´e, arXiv:1807.00294 (2018).[24] C. Huepe and M. Aldana, Phys. Rev. Lett. , 168701(2004).[25] F. Peruani, L Schimansky-Geier, and M B¨ar, Eur. Phys.J. Spec. Top. , 111 (2008).[26] F. Peruani, and M. B¨ar, New J. Phys. Introduction to Percola-tion Theory, Revised 2nd ed. (Taylor & Francis, London1994).[28] For a recent review, see H. Chat´e, Ann. Rev. Cond. Mat.Phys, to appear.[29] A. P. Solon, H. Chat´e, and J. Tailleur Phys. Rev. Lett. , 068101 (2015).[30] Both the interaction range and the timestep can berescaled to unit values without any loss of generality.[31] A. Gabrielli, F. Sylos Labini, M. Joyce, and L.Pietronero,
Statistical Physics for Cosmic Structures ,Springer Verlag Inc. (New York - Berlin, 2005).[32] G. Gr´egoire and H. Chat´e, Phys. Rev. Lett. , 025702 (2004).[33] E. Bertin, M. Droz, and G. Gr´egoire, J. Phys. A ,445001 (2009).[34] T. Ihle, Phys. Rev. E , 030901(R) (2011).[35] A.P. Solon, J. B. Caussin, D. Bartolo, H. Chat´e, and J.Tailleur, Phys Rev E , 062111 (2015).[36] F. Ginelli, Eur. Phys. J. ST, , 4819(1998).[38] This is numerically much faster (order s vs order s ina cluster of s particles) to compute than the maximumdistance between any couple of particles in the cluster. Ofcourse periodic boundary conditions are properly takeninto account computing the cluster center of mass.[39] D. Levis, and L. Berthier, Phys. Rev. E V ( t ). This will bethe longitudinal axis, and the two opposite sides of oursquare system joined by the longitudinal axis are the thelongitudinal sides. Likewise, the remaining axis and op-posite sides are the transversal ones. Clusters spanningfrom one longitudinal (transversal) side to the other arelongitudinal (transversal) spanning clusters, and we eval-uate the longitudinal (transversal) spanning probabilityas the frequency of instantaneous system configurationswith at least one longitudinal (transversal) spanning clus-ter.[41] We measure the autocorrelation time τ as the timeneeded for the autocorrelation function to drop to 10%of its equal time value.[42] Finite Size Scaling and Numerical Simulations of Statis-tical Systems , edited by V. Privman (World Scientific,Singapore, 1990).[43] D.W. Heermann and D. Stauffer, Z. Physik B , 133(1980).[44] A.A. Saberi, Phys Rep. S ( ρ ) by F ( ρ ) ≡ [1 − erf ( a ( ρ − a ))] /
2, with the size-dependentfitting parameters a and a . Obviously, ρ c ( L ) = a .[46] S. Mertens, and C. Moore, Phys. Rev. E , 061109(2012).[47] A.B. Harris, J. Phys. C
387 (1984).[48] A. Weinrib Phys. Rev. B , 387 (1984).[49] B. Mahault, F. Ginelli, and H. Chat´e, in preparation.[50] J. Toner, Phys. Rev. E , 031918 (2012).[51] In order to estimate the off-critical cut-off size Σ from theCSD, we first subtract the apparent power law decay.Σ correspond to the cluster size at which this rescaledCSD drops below the threshold 1 /e/e