Leading birds by their beaks: the response of flocks to external perturbations
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Leading birds by their beaks: the response of flocksto external perturbations
Nikos Kyriakopoulos and Francesco Ginelli
SUPA, Institute for Complex Systems and Mathematical Biology, Kings College,University of Aberdeen, Aberdeen AB24 3UE, United Kingdom
John Toner
Department of Physics and Institute of Theoretical Science, University of Oregon,Eugene, OR 97403, USA
Abstract.
We study the asymptotic response of polar ordered active fluids (“flocks”)to small external aligning fields h . The longitudinal susceptibility χ k diverges, in thethermodynamic limit, like h − ν as h →
0. In finite systems of linear size L , χ k saturatesto a value ∼ L γ . The universal exponents ν and γ depend only on the spatialdimensionality d , and are related to the dynamical exponent z and the “roughnessexponent” α characterizing the unperturbed flock dynamics. Using a well supportedconjecture for the values of these two exponents, we obtain ν = 2 / γ = 4 / d = 2and ν = 1 / γ = 2 / d = 3. These values are confirmed by our simulations. eading birds by their beaks: the response of flocks to external perturbations
1. Introduction
Flocking – the collective motion of many active particles – is a ubiquitous emergentphenomenon that occurs in many living and synthetic systems over a wide range ofscales. Examples range from mammal herds, fish schools and bird flocks to bacteriacolonies and cellular migrations, down to subcellular molecular motors and biopolymers[1]. Over the last 20 years, studies of minimal models of self-propelled particles (SPP)[2, 3, 4, 5] and hydrodynamic continuum theories [6, 7, 8, 9, 10, 11, 12] have shown thatthe behavior of typical flocking systems is essentially determined by (i) the spontaneousbreaking of continuous rotational symmetry and (ii) the far-from-equilibrium nature oflocally interacting moving particles. While the former mechanism is common to manyequilibrium systems (ranging from liquid crystals to magnetic systems and superfluidHelium-4 [13]) which spontaneously align a phase or orientational degree of freedom, thelatter is unique to active matter systems. The self-propelled motion of active particlesresults in superdiffusive information propagation even in systems without momentumconservation, which in turn leads to many striking phenomena never found in equilibriumsystems, such as long-range order in two spatial dimensions [6], and anomalously largenumber fluctuations [14].However, little is known concerning the response of moving groups to externalperturbations. This is an important question in statistical physics: symmetry breakingsystems are often characterized by their response to a small external field, and studyingresponse can also help answer the question of whether a generalized fluctuation-dissipation relation (FDR) of some sort [15] holds in flocks. Ethologists, on the otherhand, are interested in response to external threats and more generally in the biologicalsignificance of group response mechanisms. Finally, understanding response is essentialfor controlling flocking systems, either biological or artificial.In equilibrium, the response of systems breaking a continuous symmetry to a smallexternal field is a classic problem of statistical field theory, first solved in [16], where itwas shown that fluctuations transversal to order couple to longitudinal ones, yieldinga diverging longitudinal susceptibility in the entire ordered phase. This is a typicalmanifestation of symmetry-breaking, and it is a natural question to wonder how thefar-from-equilibrium nature of flocks may change this fundamental result.Until now, only a few studies, mostly numerical, have addressed these questions.Asymptotic response has been first studied numerically in the well-known Vicsek model[3], but that work focused on the behavior of the susceptibility near the transition,rather than in the ordered phase. Short time response and the dynamic FDR has beeninvestigated numerically in the Vicsek model [5] and in the isotropic phase of an activedumbbells system [17]. The response to finite and/or localized perturbations, finally,has also been studied in [18, 19, 20].Here we provide a different approach, combining hydrodynamic theory results withnumerical simulations to characterize the static response of ordered flocks to a small homogeneous external field of amplitude h . eading birds by their beaks: the response of flocks to external perturbations longitudinal response χ k ≡ δ Φ( h ) h (1)where δ Φ( h ) = Φ( h ) − Φ(0) is the change in the magnitude of the time-averaged orderparameter, which in our case is the mean velocity, due to the applied field. Our mainresult is the scaling law: χ k = h − ν f (cid:16) Lh z (cid:17) ∝ ( h − ν , L ≫ L c ( h ) L γ , L ≪ L c ( h ) , (2)where L c ( h ) ∝ h − /z and, using a conjecture first put forward in [6], ν = 4 − dd + 1 , z = 2( d + 1)5 , γ = 2(4 − d )5 (3)for any dimension 3 / ≤ d ≤
4, the upper critical dimension. In particular, we have ν = 2 / z = 6 / γ = 4 / d = 2 and ν = 1 / z = 8 / γ = 2 / d = 3. For d >
4, on the other hand, we predict δ Φ ∝ h .In the remainder of this paper, we will first derive the results (2) and (3) analytically,and then present numerical simulations that confirm them.
2. Response theory
We consider “dry” flocks, by which we mean flocks which move over a or through a static dissipative substrate or medium that acts as a momentum sink. Total momentum,thus, is not conserved, and no long ranged hydrodynamic interactions are present in thesystem. Obviously, Galilean invariance is broken, since the reference frame in which thestatic substrate or medium is at rest is preferred.
The hydrodynamic theory describes flocking by continuous, coarse grained numberdensity ρ ( r , t ) and velocity v ( r , t ) fields. The hydrodynamic equations of motiongoverning these fields in the long-wavelength limit can be obtained either by symmetryarguments [6, 7, 8, 9], or by kinetic theory [10, 11] and describe the asymptotic dynamicsof polar flocks regardless of the precise nature of the interactions, provided only thatthey are local; in particular, the same hydrodynamic equations apply for both ”metric”and ”topological” interactions[21, 22]. They are ∂ t ρ + ∇ · ( v ρ ) = 0 (4)and, in a schematic notation ∂ t v + Λ [ ∇ vv ] = U ( ρ, | v | ) v + D [ ∇∇ v ] + F P + f + h (5)where Λ [ ∇ vv ] ≡ λ ( v · ∇ ) v + λ ( ∇ · v ) v + λ ∇ ( | v | ) (6) eading birds by their beaks: the response of flocks to external perturbations λ = 1 and λ = λ = 0 as inusual Navier-Stokes equations.The viscous terms D [ ∇∇ v ] ≡ D ∇ ( ∇ · v ) + D ( v · ∇ ) v + D ∇ v (7)reflect the tendency of localized fluctuations in the velocities to spread out because oflocal interactions.The pressure term F P ≡ −∇ P − v ( v · ∇ P ) (8)is the sum of an isotropic and anisotropic pressure terms, the latter being a genuinelynon-equilibrium feature. Both terms tend to suppress local density fluctuations aroundthe global mean value ρ . The pressures P , , and the convective and viscous parameters λ k and D k > k = 1 , ,
3) are functions of the local density ρ and the magnitude | v | of the local velocity.Fluctuations are introduced through a Gaussian white noise f with correlations h f i ( r , t ) f j ( r ′ , t ′ ) i = ∆ δ ij δ d ( r − r ′ ) δ ( t − t ′ ) (9)This accounts in a simple way for any source of microscopic fluctuations, such as themicroscopic noise term opposing order in simple SPP models ‡ . Finally, the local term U simply makes the local v have a nonzero magnitude v ( h ) in the ordered phase. Itsatisfies the condition U > | v | < v ( h = 0), U = 0 for | v | = v (0), and U < | v | > v (0). This term thereby spontaneously breaks rotational symmetry even in theabsence of an external field. Small departures of the statistics of the noise from theseassumptions, e.g., slightly non-Gaussian statistics, or the introduction of “local color”in the sense of short-ranged spatio-temporal correlations of the noise, change none ofthe long distance scaling properties of the flock.Eqs. (4)-(5) are identical to the unperturbed ones discussed in [12], except for theexplicit addition of the coarse-grained constant field h in Eq. (5). By analyticity androtational invariance, this field is linearly and isotropically proportional to the appliedmicroscopic field when those fields are sufficiently small. We first discuss the system in the absence of fluctuations. Eqs. (4)-(5) admit a spatiallyuniform steady state solution ρ ( r , t ) = ρ v ( r , t ) = v ( h ) (10) ‡ One can argue that the fluctuating term arising from direct coarse-graining of such models is typicallymultiplicative (i.e., with correlations proportional to the density) rather than addictive[21]. Thisdifference, however, is irrelevant for the asymptotic properties discussed here, because the local densityfluctuations (not to be confused with the giant number fluctuations) in the TT phase are small comparedto the mean density. eading birds by their beaks: the response of flocks to external perturbations e k be the unit vector along h ≡ h ˆ e k , while forstrictly zero field ˆ e k will be the direction of the spontaneous symmetry breaking. Wehave v ( h ) = v ( h )ˆ e k , with the magnitude v ( h ) of the homogeneous velocity v ( h )determined by the condition U ( v ( h )) , ρ ) v ( h ) + h = 0 . (11)Since U is analytic in v , we have for small fields v ( h ) − v (0) ∝ h , (12)where v (0) is the zero field symmetry broken solution. It is well known that sufficientlydeep in the ordered phase such a zero-field solution is stable against spatial perturbations[11]. In the following we will restrict our analysis to this so-called Toner-Tu (TT) phase.To summarize: in mean field theory, the magnitude of the order parameterΦ( h ) ≡ |h v ( r , t ) i| (13)(here and hereafter h·i denotes a global average in space and time) responds linearly in h . We now move beyond mean field to consider the effect of fluctuations; we will show thatthe corrections to the order parameter Φ due to fluctuations are much larger than linearones we’ve just computed at mean field level.In order to do so, we allow for small fluctuations around the homogeneous solution, ρ ( r , t ) = ρ + δρ ( h ; r , t ) v ( r , t ) = v ( h ) + δ v ( h ; r , t ) , (14)and distinguish between longitudinal and transverse velocity fluctuations, which arerespectively parallel ( k ) and perpendicular ( ⊥ ) to v ( h ), δ v ( h ; r , t ) = δv k ( h ; r , t ) e k + v ⊥ ( h ; r , t ) (15)where we have made explicit the field dependence of fluctuations. For simplicity, wewill hereafter often not explicitly display the space, time and field dependence of thefluctuations.Note that, due to number conservation, h δρ i = 0, while symmetry considerationsimply h v ⊥ i = ; that is, fluctuations can’t steer the global average of h v ( r , t ) i awayfrom the external field direction. This implies that corrections to the order parameterare linear in the longitudinal fluctuations: By making use of Eqs. (12), (13) and (15)we have Φ( h ) = v ( h ) + h δ v ( h ) i = v (0) + h δv k ( h ) i + O ( h ) (16)In order to compute longitudinal fluctuations, we have to expand the hydrodynamicequations (4)-(5) in the small fluctuations δρ , δv k and v ⊥ . We are interested only influctuations that vary slowly in space and time (indeed the hydrodynamic equations eading birds by their beaks: the response of flocks to external perturbations h with those for the zero-field case in Ref. [12]. The only difference at O ( h ) is that U ( v ( h ) , ρ ) ≈ − hv = 0 . (17)Because slow modes dominate the long-distance behavior, and, it therefore proves, thesmall field response, we can eliminate the longitudinal fluctuations from Eqs. (22),since they are a fast mode of the dynamics. The subtle details of this elimination aregiven in Appendix A; the result is that the longitudinal velocity fluctuation becomes“enslaved” to the slow modes (that is, its instantaneous value is entirely determined bythe instantaneous values of those slow modes) via the relation δv k ≈ − | v ⊥ ( h ) | v (0) + µ δρ + ( µ ∂ t + µ ∂ k ) δρ + µ ∇ ⊥ · v ⊥ + O ( h ) . (18)Here µ is a constant which depends on the form of U . In typical flocking models withmetric interactions µ > ∇ ⊥ denotes spatialderivatives in the transverse directions, and the constants µ , µ and µ depend on theoriginal parameters of the hydrodynamic equations (4)-(5). Full details, together withthe derivation of Eq. (18), can be found in Appendix A, but the exact form of theseconstants is unimportant here. Since these derivative terms are linear in δρ and v ⊥ ,they vanish once averaged over space and time, so that from Eq. (18) we have h δv k i ≈ − h| v ⊥ ( h ) | i v (0) + O ( h ) (19)which links the global average of transversal and longitudinal fluctuations and is theanalogous of the so-called principle of conservation of the modulus in an equilibriumferromagnet [16]. From Eq. (16) we finally haveΦ( h ) ≈ v (0) − h| v ⊥ ( h ) | i v (0) + O ( h ) (20)and δ Φ( h ) ≡ Φ( h ) − Φ(0) ≈ h| v ⊥ (0) | i − h| v ⊥ ( h ) | i v (0) + O ( h ) (21)We are then left with the problem of determining the fluctuations of the transversevelocity v ⊥ in the presence of a non-zero field h . We will do so by analyzing the equationsof motion for v ⊥ and δρ , which follow from inserting (18) into the velocity equation ofmotion (5) projected transverse to the direction of mean motion, and into the densityequation of motion (4), and expanding in fields and derivatives. Again, details arerelegated to Appendix A; the result is: ∂ t δρ = [ ∂ t δρ ] h =0 (22) ∂ t v ⊥ = [ ∂ t v ⊥ ] h =0 − h v v ⊥ eading birds by their beaks: the response of flocks to external perturbations h v ≡ hv (0) (23)and [ ∂ t δρ ] h =0 and [ ∂ t v ⊥ ] h =0 are the terms originally given by Eqs. (2.18) and (2.28) of[12] for the zero field case. For later use, we denote collectively the parameters appearingin those terms as { µ (0) i } . While the exact forms of Eqs. (22) are not important for whatfollows, for completeness we also give them in Appendix A. We have shown so far that the response δ Φ is determined by the global average oftransversal fluctuations, Eq. (21). To compute this quantity, we proceed by a dynamicalrenormalization group (DRG) analysis [23] of Eqs. (22). Once again, this standardanalysis is almost identical to that carried out in [12] for the zero field case. We startby averaging the equations of motion over the short-wavelength fluctuations: i.e., thosewith support in the “shell” of Fourier space b − Λ ≤ | q ⊥ | ≤ Λ, where Λ is an “ultra-violetcutoff”, and b is an arbitrary rescaling factor. Then, one rescales lengths, time, δρ and v ⊥ in equations (22) according to v ⊥ = b α v ′ ⊥ , δρ = b α δρ ′ , r ⊥ = b r ′⊥ , r k = b ζ r ′ k , and t = b z t ′ to restore the ultra-violet cutoff to Λ § . The scaling exponents α , ζ , and z ,known respectively as the “roughness”, “anisotropy”, and “dynamical” exponents, areat this point arbitrary.This DRG process leads to a new, “renormalized” pair of equations of motion of thesame form as (22), but with “renormalized” values of the parameters, { µ (0) i } → { µ ( b ) i } .For a suitable choice of the scaling exponents α , ζ , and z , these parameters flow to fixed,finite limits as b → ∞ ; that is, { µ ( b →∞ ) i } → { µ ∗ i } ; this is referred to as a “renormalizationgroup fixed point”. The utility of this choice will be discussed in a moment.Since all terms except the h term in Eqs. (22) are rotation invariant, they can onlygenerate other rotation invariant terms in the first (averaging) step of the DRG. Hence,they cannot renormalize h , which breaks rotation invariance. Thus, the only change inthe h term in Eqs. (22) occurs in the second (rescaling) step. Since the coefficient h v scales as the inverse of time, this is easily seen to lead to the recursion relation h v = b − z h ′ v , (24)which – for the reasons just given – is exact to linear order in h .By construction, the DRG has the property that correlation functions in the originalequations of motions can be related to those of the renormalized equations of motion via § One could more generally rescale δρ with a different rescaling exponent α ρ from the exponent α usedfor v ⊥ . However, since fluctuations of δρ and v ⊥ have the same scaling with distance and time, theyprove to rescale with the same exponent α [12]. Note also that the exponent we call α here is called χ in most of the literature; we have broken this convention here to avoid confusion with the susceptibility χ . eading birds by their beaks: the response of flocks to external perturbations
8a simple scaling law. The example of interest for our problem is of course the correlationfunction C (cid:0) L ⊥ , L k , { µ i } , h v (cid:1) ≡ h| v ⊥ (0) | i − h| v ⊥ ( h ) | i (25)Here L ⊥ and L k are respectively the transverse and longitudinal system size. The DRGscaling law obeyed by C is thus C (cid:0) L ⊥ , L k , { µ i } , h v (cid:1) = b α C (cid:0) b − L ⊥ , b − ζ L k , { µ bi } , b z h v (cid:1) (26)which follows simply from the fact that C involves two powers of v ⊥ , each of which givesa factor b α .In order to examine the scaling of C with field amplitude h , we use the completelystandard[23, 24] renormalization group trick of choosing the scaling factor b such that b z h v is equal to some constant reference field strength h ∗ v , which we will always choose tohave the same value regardless of the bare value of h v . This implies that b = (cid:16) h v h ∗ v (cid:17) − /z .Note that for small h – and thus small h v – this choice implies b ≫
1, and that theparameters { µ ( b ) i } flow to { µ ∗ i } , their fixed point values. Hence, in the limit of small h ,the scaling function (26) can be reduced to C ( L ⊥ , L k , { µ i } , h ) = h − α/z g ( L ⊥ h /z , L k h ζ/z ) (27)where g ( x, y ) ≡ b α/z C ( b − /z x, b − ζ/z y, { µ ∗ i } , h ∗ v ) with b ≡ h ∗ v v (0) , (28)is a universal scaling function (since we always make the same choice of h ∗ v ).Note that this expression only applies for small h , since it is only in that limit that b → ∞ , and, hence { µ ( b ) i } → { µ ∗ i } . Hence, we expect this scaling law to break down forlarge fields, and, in fact, it does.We now focus our attention on roughly square systems, with L ⊥ ∼ L k ∼ L .Assuming an anisotropy exponent 0 < ζ < d < Lh /z ∝ L ⊥ h /z ≪ L k h ζ/z (29)so that finite-size scaling is controlled by the transverse flock extension L ⊥ and we canreplace g with the universal scaling function w ( x ) ≡ g ( x, ∞ ). Above the upper criticaldimension d c = 4, where ζ = 1 [6], scaling is the same in both the transversal andlongitudinal directions and we choose instead w ( x ) ≡ g ( x, x ). Doing so, we finallyobtain the scaling law C ( L, h ) = h − ν ˜ w ( Lh /z ) (30)where ν = 1 + 2 α/z . (31)It is now straightforward to relate the order parameter change δ Φ( h ) to this scaling law.From Eq. (21) we have δ Φ ∝ C ( L, h ) + O ( h ) . (32) eading birds by their beaks: the response of flocks to external perturbations ν > d < C ( L, h ) ≫ O ( h ) and the corrections due to fluctuations dominate themean field ones. To lowest order in hδ Φ = 12 v (0) (cid:2) h| v ⊥ (0) | i − h| v ⊥ ( h ) | i (cid:3) = h − ν w ( Lh /z ) (33)In the thermodynamic limit, Lh /z ≫ w ( Lh /z ) → w ( ∞ ) =constant, yielding the asymptotic result δ Φ h ∝ h − ν . In practice, in any system largeenough, the external field suppresses transverse fluctuations, thus increasing the scalarorder parameter according to Eq. (20), an effect that below the upper critical dimension d c = 4 proves stronger than mean field corrections linear in h . Above d c , on the otherhand, it is known [6] that α = 1 − d/ z = 2, which implies ν = 2 − d/ ≤
0; thereforecorrections due to fluctuations no longer dominate the mean field ones. Hence, ordinarylinear response χ k → constant as h → d > d < α , ζ , and z for which the DRG flowsto a fixed point in those dimensions. These values actually do depend on the formof Eqs. (22) (and, in particular, to the nature of their relevant nonlinear terms), butluckily for small fields h they have to coincide with their zero-field values. Indeed, thereis no reason for which these zero-field values should be affected by a sufficiently smallrescaled field, h v ≤ h ∗ v .In [6], it was argued that for any dimension 3 / ≤ d ≤ α = 3 − d , z = 2( d + 1)5 , ζ = d + 15 . (34)It has since been since realized [12] that the original arguments leading to these valuesare flawed. However, the simple conjecture that the only relevant non-linearity at thefixed point is the term proportional to λ in Eqs. (5)-(6) leads to precisely these valuesfor the exponents. While this conjecture has never been proven, there is solid numerical[5, 7] and even experimental [25] evidence supporting the above scaling exponent valuesfor d = 2 and, to a lesser extent, d = 3. In the following we will assume this conjectureholds, verifying it a-posteriori by numerical measures of asymptotic response in theVicsek [2] model.Above the upper critical dimension, d c = 4, finally, the scaling exponents take theexact linear values z = 2, ζ = 1 and α = 1 − d/ We conclude this section discussing finite size effects. The scaling form (33) implies thattransverse fluctuations are suppressed by the field h on length scales L ≫ L c ( h ) ∝ h − /z . (35) eading birds by their beaks: the response of flocks to external perturbations L ≪ L c (or equivalently for small field h ≪ h c ( L ) ∝ L − z ),however, this suppression is ineffective, and leading corrections to the order parametershould revert to linear order in the external field h . We can include this behaviorin a single universal scaling function f by requiring that f ( ∞ ) ∝ w ( ∞ ) = O (1) and f ( x ) ∝ x νz for x ≪
1. This finally gives, δ Φ = h − ν f ( h L z ) ∝ ( h − ν , h ≫ L − z hL γ , h ≪ L − z . (36)with γ = νz . This scaling holds for external fields not too large. For h > h ∗ v , onthe other hand, the small field approximation discussed here is no longer valid, andsaturation effects change the scaling (36). Once expressed in terms of the longitudinalsusceptibility χ k = δ Φ /h , our results imply Eq. (2), with the scaling exponents givenby Eq. (3) (according to conjecture (34)).
3. Numerical simulations
We test our predictions (36) in two and three dimensions by simulating the well knownVicsek model [2] in an external homogeneous field h . Each particle – labeled by i = 1 , , . . . , N – is defined by a position r ti and a unit direction of motion v ti . Themodel evolves with a synchronous discrete time dynamics v t +1 i = ( R ω ◦ ϑ ) X j ∈ S i v tj + h ! (37) r t +1 i = r ti + v m v ti (38)where v m is the particle speed k , ϑ ( w ) = w / | w | is a normalization operator and R ω performs a random rotation (uncorrelated between different times t or particles i )uniformly distributed around the argument vector: R ω w is uniformly distributed around w inside an arc of amplitude 2 πω (in d = 2) or in a spherical cap spanning a solid angleof amplitude 4 πω ( d = 3). The interaction is “metric”: that is, each particle i interactswith all of its neighbors within unit distance. In the following, we adopt periodicboundary conditions and choose typical microscopic parameters so that the system lieswithin the TT phase [5]: v m = 0 . ρ = N/L d = 1 and ω = 0 .
18 ( d = 2) or ω = 0 . d = 3). In both dimensions, this choice yields a zero-field order parameter Φ(0) ≈ . L . After discarding a transient T sufficiently long for the system tosettle into the stationary state, we estimate the mean global order parameterΦ = * N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N X i =1 v ti (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)+ t (39)and its standard error, given by S E = σ/ √ n , with σ being the standard deviationand n the number of independent data points. We estimate n as the total number of k Note that v m = 0 is the equilibrium limit of this model [26], which is singular . eading birds by their beaks: the response of flocks to external perturbations T divided by the autocorrelation time τ of the mean order parametertimeseries, n = T /τ . In Eq. (39), h·i t denotes time averages, performed over typically T = 10 ∼ time steps. In particular, as the precision of the zero-field order parameteraffects all response and the autocorrelation time decreases with h , in our numericalsimulations we take care to estimate the zero-field order parameter over times as largeas possible.We begin addressing response in the large system size (or large field) regime hL z ≫
1, where our theory predicts δ Φ ∼ h − ν . Measuring this power law is aparticularly difficult task, as it is sandwiched between saturation effects at larger valuesof h and the crossover to linear behavior at h ≪ L − z . We proceed by extrapolation,choosing a (somewhat arbitrary) h range of two decades and by measuring the effectivepower law exponent 1 − ν eff by linear regression. The resulting response δ Φ( h ) isplotted in Fig. 1 for increasing system sizes L . As one expects, as system sizes increases,response curves approach the expected size-asymptotic behavior δ Φ ∼ h − ν . In Fig. 1a,for instance, the d = 2 response approaches the expected power-law 1 − ν = 1 /
3. Inthe inset of Fig. 1a we further quantify this convergence plotting | ν eff − ν | vs. thesystem size L . This shows that the effective exponent approaches the predicted onewith corrections of order 1 / √ L . We repeated the same procedure in d = 3. As shownin Fig. 1b, the approach to the expected asymptotic exponent 1 − ν = 3 / | ν eff − ν | vanishes faster, as L − . . In d = 3 our simulations areobviously limited to a much smaller range of linear size values, but it should be notedthat in d = 3 finite size effects vanish quicker (being the exponent z larger) while theasymptotic exponent 1 − ν = 3 / δ Φ ∼ h expectedat low values of h . A faster approach of ν eff to its asymptotic value is therefore notcompletely surprising.In d = 3, we can also easily compare the response behavior with the equilibriumprediction δ Φ ∼ √ h [16] (dashed blue line in Fig. 1b). This clearly shows that thefar-from-equilibrium nature of the Vicsek model makes the susceptibility exponent verydifferent from that in equilibrium ferromagnets. In d = 2 equilibrium systems with acontinuous symmetry cannot develop long range order, but, rather, exhibit only a quasi-long range ordered phase, characterized by scaling exponents that vary continuously withtemperature [13, 27]. The equilibrium susceptibility exponent [27, 28] in d = 2 is givenby ν = 4 − η − η , (40)where the order parameter correlation exponent η is bounded: 0 ≤ η ≤ /
4, whichimplies that the susceptibility exponent ν varies over an extremely narrow range:14 / ≤ ν ≤
1. Our predicted value of ν = 2 / δ Φ ∝ h predicted for small system sizes (orsmall fields), hL z ≪
1. This inequality imposes a (severe) upper limit on the range of h , eading birds by their beaks: the response of flocks to external perturbations -3 -2 -1 h -3 -2 -1 δΦ -3 -2 -1 -3 -2 -1 (a) -3 -2 h -4 -3 -2 δΦ -3 -2 -4 -3 -2 (b) L -1 | ν e ff − ν | L -1 | ν e ff − ν | Figure 1. (color online) Size-asymptotic regime – Order parameter change vs. theapplied field amplitude for different system sizes. a) For d = 2, from bottom to top, L = 32 , , , , , δ Φ ∼ h / , while the dashed blue line marks the upper bound forthe d = 2 equilibrium response δ Φ ∼ h / . In the inset: Absolute difference betweenthe measured effective exponent ν eff (see text) and its expected asymptotic value ν asa function of system size. The dashed black lines marks a power law decay as 1 / √ L .(b) For d = 3, from bottom to top, L = 40 , , , , δ Φ ∼ h / , while the dashed blue linecorrespond to the d = 3 equilibrium response δ Φ ∼ h / . Inset: d = 3 data as in theinset of panel (a). The dashed black lines marks a power law decay as L − . . Errorbars measure standard errors (see text). All graphs are in a double logarithmic scale. while there is a lower limit set by our numerical precision in evaluating responses of theorder of 10 − or smaller. Nevertheless, our numerical simulations reveal in both d = 2(Fig. 2a) and d = 3 (Fig. 2b) a linear growth of the response over more than one decadein the field amplitude h , especially for small system sizes.By selecting a single h value lying in the linear regime for all accessible systemsizes L , one is also able to test the saturation exponent gamma, δ Φ ∼ h L γ . This isdone in Figs. 2c-d, where response values at different linear sizes L are compared to thepredicted power-law with (respectively) γ = 4 / d = 2 and γ = 2 / d = 3. Weobtain a good agreement in d = 2 (the best linear fit being γ = 0 . d = 3 is less clear, in rough agreement with the expected power-law behavior only forsizes L ≥
60, with a best linear fit of γ = 0 . h larger than h s ≈ .
1, of course, are out of the small fieldregime and show saturation effects, while due to statistical fluctuations we have beenunable to obtain reliable estimates for external fields smaller than h ≈ − . Withinthis range, comparison with the predicted scaling (2) (as given by the dashed lines) isoverall satisfactory, especially in d = 2. By a proper rescaling, making use of the threescaling exponents (3), we can also collapse our data at different sizes on roughly a single eading birds by their beaks: the response of flocks to external perturbations -4 -3 h -3 -2 δΦ (a) L -3 -2 δΦ (c) -3 -2 h -4 δΦ (b)
40 100 L -4 -3 δΦ (d) Figure 2. (color online) Linear regime. (a)-(b) Order parameter change vs. theapplied field amplitude in the linear regime for different system sizes (the cyan arrowindicates increasing system sizes): a) For d = 2 L = 32 , , , , d = 3 L = 40 , , , , δ Φ ∼ h . Error barsmeasure data standard errors (see text). (c) d = 2 Response at fixed h – as shownby the red arrows in panel (a) – and different system sizes in the linear regime. Thedashed black line marks a power law with the predicted slope 0 .
8. (d) Same as in (c),but for d = 3. The dashed black line marks a power law with the predicted slope 0 . curve, as shown in Fig. 3c-d.To summarize, numerical simulations are in good agreement with our theoreticalpredictions, at least in d = 2. Results in d = 3 are prone to larger errors and obviouslyexplore a more limited range of linear sizes, but nevertheless are still compatible withour predictions.We also performed a few additional numerical studies of response (not shown here)with different parameter values (but still in the TT phase), and in the ordered phase ofthe so-called topological Vicsek model [29], confirming the generality of these results.It is also worth commenting on the way the external field is implemented in themicroscopic Vicsek equations (38). In Ref. [5] it was argued that different microscopicimplementations could lead to different response, and in particular it was recommendedto choose one by which the external field was normalized by the local order parameter eading birds by their beaks: the response of flocks to external perturbations -4 -2 h -2 δΦ (a) -2 h L -2 χ // L -0.8 (c) -4 -3 -2 -1 h -4 -2 δΦ (b) h L -2 -1 χ // L -0.4 (d) Figure 3. (color online) (a)-(b) Order parameter change vs. the applied fieldamplitude in both the linear and size-asymptotic regimes for different system sizes (thearrow indicates increasing system sizes): a) For d = 2 L = 32 , , , , , d = 3 L = 40 , , , , d = 2 and (d) d = 3. In all panels, the dashed black lines markthe linear response ( δ Φ ∝ h or χ k ∼ L γ expected by our theory for h ≪ h c ∼ L − z .Dashed red lines, on the other hand, mark the nonlinear regime predicted for h ≫ h c , δ Φ ∝ h − ν or χ k ∝ h − ν . We have ν ( d = 2) = 2 /
3, and ν ( d = 3) = 1 /
4. Error barsmeasure standard errors (see text). All graphs are in a double logarithmic scale. value, such as in v t +1 i = ( R ω ◦ ϑ ) " ϑ X j ∈ S i v tj ! + h (41)However, we do not expect these microscopic details to change the structure of thehydrodynamic equations (4)-(5), and thus we do not expect qualitative differencesbetween the two microscopic external field implementations. Indeed, our preliminarysimulations of Eq. (41) (not shown here) show no qualitative difference from the responseextensively discussed above for Eq. (38). eading birds by their beaks: the response of flocks to external perturbations
4. Conclusions
So far, we have only considered the longitudinal susceptibility, which characterizesthe response of the magnitude of the order parameter to a small external field.Simple considerations, based on symmetry, imply that the flock polarization vectorwill eventually align with any non-zero stationary external field, including one appliedtransversal to the initial polarization. Thus, for the transvere susceptibility we have,trivially, χ ⊥ ∼ v (0) h (42)as in equilibrium systems.In this paper, we have fully characterized the static response of homogeneousordered flocks to small external fields for any dimension d > /
2. In particular, belowthe upper critical dimension d c = 4, our results in the thermodynamic limit L → ∞ show a diverging longitudinal response for h →
0, i.e. a diverging susceptibility. Thisis ultimately a consequence of the spontaneous symmetry breaking of the continuousrotation symmetry, albeit the far-from equilibrium nature of flocks yields different resultsfrom, say, equilibrium ferromagnets in d = 2 , bona fide TT phase. This class encompasses both systems with metric interactions andthose with topological interactions. It also includes the inertial spin model recently putforward in [30] to account for the turning dynamics measured experimentally in starlingflocks. This is because the long time hydrodynamic theory of the inertial spin theoryrelaxes to the TT theory [31]; hence, the static response will be unchanged. Dynamicalresponse (i.e. how quickly the flock turns towards the field direction), however, couldbe different at short times in inertial spin models, while the long-time behavior shouldbe the same.In future work, we will explore more thoroughly the phase diagram of Vicsek-likemodels beyond the TT phase, investigating the disordered and phase separated regimesOther future directions include the study of the finite-time, dynamical reponse [32] inboth overdamped Vicsek-like models and inertial spin ones, and the study of spatiallylocalized perturbations. eading birds by their beaks: the response of flocks to external perturbations Acknowledgments
We have benefited from discussions with H. Chat´e and A. Cavagna. We acknowledgesupport from the Marie Curie Career Integration Grant (CIG) PCIG13-GA-2013-618399. JT also acknowledges support from the SUPA distinguished visitor program andfrom the National Science Foundation through awards
Appendix A. Expansion of the Hydrodynamic equations for smallfluctuations
We first demonstrate that the longitudinal velocity δv k is enslaved to the slow modes δρ and v ⊥ . We follow Ref. [12] and begin with the hydrodynamic Eqs. (4) and (5) writtenout explicitly: ∂ t ρ + ∇ · ( v ρ ) = 0 (A.1) ∂ t v = − λ ( v · ∇ ) v − λ ( ∇ · v ) v − λ ∇ ( | v | ) + U ( ρ, v ) v − ∇ P − v ( v · ∇ P ) + D ∇ ( ∇ · v ) + D ∇ v + D ( v · ∇ ) v + f + h , (A.2)where, as noted in the main text, the parameters λ i ( i = 1 → U , the“isotropic Pressure” P ( ρ, | v | ) and the “anisotropic Pressure” P ( ρ, | v | ) are, in general,functions of the density ρ and the magnitude | v | of the local velocity. It is useful toTaylor expand P and P around the equilibrium density ρ : P = ∞ X n =1 σ n ( | v | )( ρ − ρ ) n (A.3) P = P ( ρ, | v | ) = ∞ X n =1 Υ n ( | v | )( ρ − ρ ) n . (A.4)Here D , D and D are all positive in the ordered state.Again as discussed in the main text, in the ordered phase, the velocity field can bewritten as: v = v e k + δ v = ( v + δv k ) e k + v ⊥ (A.5)(for simplicity, here and hereafter, we write v ≡ v ( h )).Taking the dot product of both sides of equation (A.2) with v itself, we obtain:12 (cid:0) ∂ t | v | (cid:1) + 12 (cid:0) λ + 2 λ )( v · ∇ ) | v | (cid:1) + λ ( ∇ · v ) | v | = U ( | v ) | v | − v · ∇ P − | v | v · ∇ P + D v · ∇ ( ∇ · v ) + D v · ∇ v + D v · (cid:0) ( v · ∇ ) v (cid:1) + v · f + v · h . (A.6) eading birds by their beaks: the response of flocks to external perturbations δ v ( r , t ) and δρ ( r , t ) that vary slowly in space and time. (Indeed, the hydrodynamic equations (A.2)and (A.1) are only valid in this limit). Hence, terms involving space and time derivativesof δ v ( r , t ) and δρ ( r , t ) are always negligible, in the hydrodynamic limit, comparedto terms involving the same number of powers of fields without any time or spacederivatives.Furthermore, the fluctuations δ v ( r , t ) and δρ ( r , t ) can themselves be shown to besmall in the long-wavelength limit [12]. Hence, we need only keep terms in equation(A.6) up to linear order in δ v ( r , t ) and δρ ( r , t ). The v · f term can likewise be dropped,since it only leads to a term of order v ⊥ f k in the v ⊥ equation of motion, which isnegligible (since v ⊥ is small) relative to the f ⊥ term already there.In addition, treating the magnitude h of the applied field as a small quantity, weneed only keep terms involving h that are proportional to h and independent of thesmall fluctuating quantities δ v and δρ .These observations can be used to eliminate many of the terms in equation (A.6),and solve for the quantity U ; the solution is: U = − hv + λ ∇ · v + v · ∇ P + σ v ∂ k δρ + 12 v (cid:16) ∂ t + λ ∂ k (cid:17) δv k , (A.7)where we’ve defined λ ≡ ( λ + 2 λ ) v . (A.8)We can now express the longitudinal velocity δv k in terms of the slow modes usingequation (A.7) and the expansion U ≈ − Γ (cid:18) δv k + | v ⊥ | v (cid:19) − Γ δρ , (A.9)where we’ve definedΓ ≡ − (cid:18) ∂U∂ | v | (cid:19) ρ , Γ ≡ − (cid:18) ∂U∂ρ (cid:19) | v | , (A.10)with, here and hereafter, super- or sub-scripts 0 denoting functions of ρ and | v | evaluatedat ρ = ρ and | v | = v . We’ve also used the expansion (A.5) for the velocity in terms ofthe fluctuations δv k and ~v ⊥ to write | v | = v + δv k + | v ⊥ | v + O ( δv k , | v ⊥ | ) , (A.11)and kept only terms that an DRG analysis shows to be relevant in the long wavelengthlimit [12]. Inserting (A.9) into (A.7) gives: − Γ (cid:18) δv k + | v ⊥ | v (cid:19) − Γ δρ = − hv + λ ∇ ⊥ · v ⊥ + λ ∂ k δv k + (Υ v + σ ) v ∂ k δρ + 12 v (cid:16) ∂ t + λ ∂ k (cid:17) δv k , (A.12) eading birds by their beaks: the response of flocks to external perturbations | v ⊥ | , and hence negligible, in thehydrodynamic limit, relative to the | v ⊥ | term explicitly displayed on the left-hand side.This equation can be solved iteratively for δv k in terms of v ⊥ , δρ , and its derivatives.To lowest (zeroth) order in derivatives, δv k ≈ µ δρ + hv Γ where we have defined µ = − Γ Γ . (A.13)Inserting this approximate expression for δv k into equation (A.12) everywhere δv k appears on the right hand side of that equation gives δv k to first order in derivatives: δv k ≈ − | v ⊥ | v + µ δρ + Γ v Γ ∂ t δρ − λ Γ ∂ k δρ − λ Γ ∇ ⊥ · v ⊥ + hv Γ , (A.14)where we’ve defined λ ≡ (Υ v + σ ) v − Γ Γ (cid:18) λ + λ v (cid:19) = (Υ v + σ ) v + µ ( λ + λ + 2 λ ) . (A.15)(In deriving the second equality in (A.15), we’ve used the definition (A.8) of λ .)Eq. (A.14) coincides with Eq. (18) in the main text with µ = Γ v Γ (A.16) µ = − λ Γ (A.17) µ = − λ Γ (A.18)and expresses the enslaving of δv k to the slow modes δρ and v ⊥ .We finally derive the equations of motion for the slow modes. Inserting theexpression (A.7) for U back into equation (A.2), we find that P and λ cancel outof the v equation of motion, leaving ∂ t v = − λ ( v · ∇ ) v − λ ∇ ( | v | ) − hv v + σ v v ( ∂ k δρ ) − ∇ P + D ∇ ( ∇ · v )+ D ∇ v + D ( v · ∇ ) v + (cid:20) v (cid:16) ∂ t + λ ∂ k (cid:17) δv k (cid:21) v + f + h . (A.19)This can be made into an equation of motion for v ⊥ involving only v ⊥ ( r , t ) and δρ ( r , t ) by i) projecting (A.19) perpendicular to the direction of mean flock motion e k ,and ii) eliminating δv k by inserting equation (A.14) into the equation of motion (A.19)for v . Using the expansions (A.5), (A.11) and neglecting “irrelevant” terms we have: ∂ t v ⊥ = − λ v ∂ k v ⊥ − λ ( v ⊥ · ∇ ⊥ ) v ⊥ − g δρ∂ k v ⊥ − g v ⊥ ∂ k δρ − c ρ ∇ ⊥ δρ eading birds by their beaks: the response of flocks to external perturbations − g ∇ ⊥ ( δρ ) + D B ∇ ⊥ ( ∇ ⊥ · v ⊥ ) + D ∇ ⊥ v ⊥ + D k ∂ k v ⊥ + g t ∂ t ∇ ⊥ δρ + g k ∂ k ∇ ⊥ δρ + f ⊥ − h v v ⊥ (A.20)where we’ve defined h v as in the main text, and D B ≡ D + 2 v λ λ Γ , (A.21) D k ≡ D + D v , (A.22) g ≡ µ λ + v (cid:18) ∂λ ∂ρ (cid:19) , (A.23) g ≡ − µ λ v − σ v , (A.24) g ≡ σ + µ λ + v µ (cid:18) ∂λ ∂ρ (cid:19) , (A.25) c ≡ ρ σ + 2 ρ v λ µ , (A.26) g t ≡ µ λ (A.27)and g k ≡ v λ λ Γ − µ D . (A.28)Finally, using (A.5), (A.11) and (A.14) in the equation of motion (A.1) for ρ gives, againneglecting irrelevant terms: ∂ t δρ = − ρ ∇ ⊥ · v ⊥ − w ∇ ⊥ · ( v ⊥ δρ ) − v ∂ k δρ + D ρ k ∂ k δρ + D ρ ⊥ ∇ ⊥ δρ + D ρv ∂ k ( ∇ ⊥ · v ⊥ ) + ρ µ ∂ t ∂ k δρ − µ ∂ k ( δρ ) + µ ∂ k ( | v ⊥ | ) , (A.29)where we’ve defined: v ≡ v + µ ρ (A.30) µ ≡ ρ v , (A.31) D ρ k ≡ ρ λ Γ = − ρ µ (A.32)and, last but by no means least, D ρv ≡ λ ρ Γ = − ρ µ , (A.33) eading birds by their beaks: the response of flocks to external perturbations D ρ ⊥ is actually zero at this point in the calculation, but we’veincluded it in equation (A.29) anyway, because it is generated by the nonlinear termsunder the Renormalization Group. Likewise, the parameter w = 1, but will changefrom that value upon renormalization.The equation of motion (A.29) is, as claimed in the main text, exactly the same asthat in the absence of the external field h , while the equation of motion (A.20) is of theform (22), with[ ∂ t v ⊥ ] h =0 ≡ − λ v ∂ k v ⊥ − λ ( v ⊥ · ∇ ⊥ ) v ⊥ − g δρ∂ k v ⊥ − g v ⊥ ∂ k δρ − c ρ ∇ ⊥ δρ − g ∇ ⊥ ( δρ ) + D B ∇ ⊥ ( ∇ ⊥ · v ⊥ ) + D ∇ ⊥ v ⊥ + D k ∂ k v ⊥ + g t ∂ t ∇ ⊥ δρ + g k ∂ k ∇ ⊥ δρ + f ⊥ (A.34) [1] S. Ramaswamy, Ann. Rev. Cond. Matt. Phy.
323 (2010).[2] T. Vicsek et al. , Phys. Rev. Lett. , 1226 (1995).[3] A. Czirok, H. E. Stanley, and T. Vicsek, J. Phys. A , 1375 (1997).[4] G. Gr´egoire and H. Chat´e, Phys. Rev. Lett. , 025702 (2004).[5] H. Chat´e, F. Ginelli, G. Gr´egoire, and F. Reynaud, Phys. Rev. E , 046113 (2008).[6] J. Toner and Y.-h. Tu, Phys. Rev. Lett. ,4326 (1995).[7] Y.-h. Tu, M. Ulm and J. Toner, Phys. Rev. Lett. , 4819 (1998).[8] J. Toner and Y.-h. Tu, Phys. Rev. E , 4828(1998).[9] J. Toner, Y.-h. Tu, and S. Ramaswamy, Ann. Phys. , 170(2005).[10] E. Bertin, M. Droz, and G. Gr´egoire, Phys. Rev. E , 022101 (2006).[11] E. Bertin, M. Droz, and G. Gr´egoire, J. Phys. A , 445001 (2009).[12] J. Toner, Phys. Rev. E , 031918 (2012).[13] V. L. Pokrovskii and A. Z. Pata˘shinskii, Fluctuation Theory of Phase Transitions, 2nd ed.(Pergamon, Oxford, 1979; Nauka, Moscow, 1982).[14] S. Ramaswamy, R. A. Simha and J. Toner, Eur. Phys. Lett.
196 (2003).[15] U. M. B. Marconi, A. Puglisi, L. Rondoni, A. Vulpiani, Phys. Rep.
111 (2008).[16] A. Z. Pata˘shinskii and V. L. Pokrovskii, Zh. Eksp. Teor. Fiz. , 1445 (1973).[17] D. Loi, S. Mossa, L. E. Cugliandolo, Phys. Rev. E , 051111 (2008); Soft Matter , 3726 (2011).[18] A. Cavagna, I. Giardina, F. Ginelli, Phys. Rev. Lett. , 168107 (2013)[19] I. D. Couzin, J. Krause, N. R. Franks, S. A. Levin, Nature , 513 (2005).[20] D. J. G. Pierce, L. Giomi, arXiv:1511.03652 (2015).[21] A. Peshkov, E. Bertin, F. Ginelli, H. Chat´e, Eur. Phys. J. Special Topics , 1315 (2014).[22] A. Peshkov, S. Ngo, E. Bertin, H. Chat´e, F. Ginelli, Phys. Rev. Lett , 098101 (2012).[23] See, e.g., D. Forster, D. R. Nelson, and M. J. Stephen, Phys. Rev. A , 732 (1977).[24] For an excellent discussion of precisely the logic we use here, see, e.g., S. K. Ma, Modern Theoryof Critical Phenomena, 1st. ed., (Benjamin, Reading, Mass., 1977).[25] K. Gowrishankar et al. , Cell , 1353 (2012).[26] F. Ginelli, arXiv:1511.01451 (2015).[27] J. M. Kosterlitz and D. J. Thouless, J. Phys. C , 1181 (1973); J. M. Kosterlitz, J. Phys. C ,1046 (1974).[28] J. V. Jos´e, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Phys. Rev. B , 1217 (1977).[29] F. Ginelli, H. Chat´e, Phys. Rev. Lett., et al. , Nat. Phys.
691 (2014).[31] A. Cavagna et al. , Phys. Rev. Lett., eading birds by their beaks: the response of flocks to external perturbations [32] L. F. Cugliandolo, J. Kurchan, and L. Peliti, Phys. Rev. E55