Leading birds by their beaks: the response of flocks to external perturbations

We study the asymptotic response of polar ordered active fluids ("flocks") to small external aligning fields $h$. The longitudinal susceptibility $\chi_{_\parallel}$ diverges, in the thermodynamic limit, like $h^{-\nu}$ as $h \rightarrow 0$. In finite systems of linear size $L$, $\chi_{_\parallel}$ saturates to a value $\sim L^\gamma$. The universal exponents $\nu$ and $\gamma$ depend only on the spatial dimensionality $d$, and are related to the dynamical exponent $z$ and the"roughness exponent"$\alpha$ characterizing the unperturbed flock dynamics. Using a well supported conjecture for the values of these two exponents, we obtain $\nu = 2/3$, $\gamma = 4/5$ in $d = 2$ and $\nu = 1/4$, $\gamma = 2/5$ in $d = 3$. These values are confirmed by our simulations.


Introduction
Flocking -the collective motion of many active particles -is a ubiquitous emergent phenomenon that occurs in many living and synthetic systems over a wide range of scales. Examples range from mammal herds, fish schools and bird flocks to bacteria colonies 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 that the behavior of typical flocking systems is essentially determined by (i) the spontaneous breaking of continuous rotational symmetry and (ii) the far-from-equilibrium nature of locally interacting moving particles. While the former mechanism is common to many equilibrium systems (ranging from liquid crystals to magnetic systems and superfluid Helium-4 [13]) which spontaneously align a phase or orientational degree of freedom, the latter is unique to active matter systems. The self-propelled motion of active particles results in superdiffusive information propagation even in systems without momentum conservation, which in turn leads to many striking phenomena never found in equilibrium systems, such as long-range order in two spatial dimensions [6], and anomalously large number fluctuations [14].
However, little is known concerning the response of moving groups to external perturbations. This is an important question in statistical physics: symmetry breaking systems are often characterized by their response to a small external field, and studying response can also help answer the question of whether a generalized fluctuationdissipation relation (FDR) of some sort [15] holds in flocks. Ethologists, on the other hand, are interested in response to external threats and more generally in the biological significance of group response mechanisms. Finally, understanding response is essential for controlling flocking systems, either biological or artificial.
In equilibrium, the response of systems breaking a continuous symmetry to a small external field is a classic problem of statistical field theory, first solved in [16], where it was shown that fluctuations transversal to order couple to longitudinal ones, yielding a diverging longitudinal susceptibility in the entire ordered phase. This is a typical manifestation of symmetry-breaking, and it is a natural question to wonder how the far-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 been investigated numerically in the Vicsek model [5] and in the isotropic phase of an active dumbbells 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 with numerical simulations to characterize the static response of ordered flocks to a small homogeneous external field of amplitude h.
We are particularly interested in the asymptotic longitudinal response where δΦ(h) = Φ(h) − Φ(0) is the change in the magnitude of the time-averaged order parameter, which in our case is the mean velocity, due to the applied field. Our main result is the scaling law: where L c (h) ∝ h −1/z and, using a conjecture first put forward in [6], for any dimension 3/2 ≤ d ≤ 4, the upper critical dimension. In particular, we have ν = 2/3, z = 6/5 and γ = 4/5 in d = 2 and ν = 1/4, z = 8/5 and γ = 2/5 in 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.

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 the system. Obviously, Galilean invariance is broken, since the reference frame in which the static substrate or medium is at rest is preferred.

Hydrodynamic description
The hydrodynamic theory describes flocking by continuous, coarse grained number density ρ(r, t) and velocity v(r, t) fields. The hydrodynamic equations of motion governing these fields in the long-wavelength limit can be obtained either by symmetry arguments [6,7,8,9], or by kinetic theory [10,11] and describe the asymptotic dynamics of polar flocks regardless of the precise nature of the interactions, provided only that they 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 where are all the convective-like terms permitted by the symmetries and conservation laws of the system. Here, all three coefficients are, in general, neither zero nor one, as opposed to systems with Galilean invariance where one has simply λ 1 = 1 and λ 2 = λ 3 = 0 as in usual Navier-Stokes equations.
The viscous terms reflect the tendency of localized fluctuations in the velocities to spread out because of local interactions. The pressure term is the sum of an isotropic and anisotropic pressure terms, the latter being a genuinely non-equilibrium feature. Both terms tend to suppress local density fluctuations around the global mean value ρ 0 . The pressures P 1,2 , and the convective and viscous parameters λ k and D k > 0 (k = 1, 2, 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 This accounts in a simple way for any source of microscopic fluctuations, such as the microscopic noise term opposing order in simple SPP models ‡. Finally, the local term U simply makes the local v have a nonzero magnitude v 0 (h) in the ordered phase. It satisfies the condition U > 0 for |v| < v 0 (h = 0), U = 0 for |v| = v 0 (0), and U < 0 for |v| > v 0 (0). This term thereby spontaneously breaks rotational symmetry even in the absence of an external field. Small departures of the statistics of the noise from these assumptions, 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 of the long distance scaling properties of the flock. Eqs. (4)-(5) are identical to the unperturbed ones discussed in [12], except for the explicit addition of the coarse-grained constant field h in Eq. (5). By analyticity and rotational invariance, this field is linearly and isotropically proportional to the applied microscopic field when those fields are sufficiently small.

Mean-field analysis
We first discuss the system in the absence of fluctuations. Eqs. (4)-(5) admit a spatially uniform steady state solution ρ(r, t) = ρ 0 v(r, t) = v 0 (h) (10) ‡ One can argue that the fluctuating term arising from direct coarse-graining of such models is typically multiplicative (i.e., with correlations proportional to the density) rather than addictive [21]. This difference, however, is irrelevant for the asymptotic properties discussed here, because the local density fluctuations (not to be confused with the giant number fluctuations) in the TT phase are small compared to the mean density.
For any nonzero external field letê be the unit vector along h ≡ hê , while for strictly zero fieldê will be the direction of the spontaneous symmetry breaking. We have v 0 (h) = v 0 (h)ê , with the magnitude v 0 (h) of the homogeneous velocity v 0 (h) determined by the condition Since U is analytic in v, we have for small fields where v 0 (0) is the zero field symmetry broken solution. It is well known that sufficiently deep 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 (here and hereafter · denotes a global average in space and time) responds linearly in h.

Fluctuations
We now move beyond mean field to consider the effect of fluctuations; we will show that the corrections to the order parameter Φ due to fluctuations are much larger than linear ones we've just computed at mean field level. In order to do so, we allow for small fluctuations around the homogeneous solution, and distinguish between longitudinal and transverse velocity fluctuations, which are respectively parallel ( ) and perpendicular (⊥) to v 0 (h), where we have made explicit the field dependence of fluctuations. For simplicity, we will hereafter often not explicitly display the space, time and field dependence of the fluctuations. Note that, due to number conservation, δρ = 0, while symmetry considerations imply v ⊥ = 0; that is, fluctuations can't steer the global average of v(r, t) away from the external field direction. This implies that corrections to the order parameter are linear in the longitudinal fluctuations: By making use of Eqs. (12), (13) and (15) we have In order to compute longitudinal fluctuations, we have to expand the hydrodynamic equations (4)-(5) in the small fluctuations δρ, δv and v ⊥ . We are interested only in fluctuations that vary slowly in space and time (indeed the hydrodynamic equations are only valid in this limit), so that space and time derivatives of the fluctuations are always of higher order than the fluctuations themselves. The details of this surprisingly subtle calculation are given in Appendix A, but (fortunately) they are identical to order h with those for the zero-field case in Ref. [12]. The only difference at O(h) is that Because slow modes dominate the long-distance behavior, and, it therefore proves, the small 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 are given 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 by the instantaneous values of those slow modes) via the relation Here µ 1 is a constant which depends on the form of U. In typical flocking models with metric interactions µ 1 > 0 [11], so that density fluctuations are positively correlated with longitudinal fluctuations at the local level. In equation (18), ∇ ⊥ denotes spatial derivatives in the transverse directions, and the constants µ 2 , µ 3 and µ 4 depend on the original parameters of the hydrodynamic equations (4)- (5). Full details, together with the derivation of Eq. (18), can be found in Appendix A, but the exact form of these constants 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 which links the global average of transversal and longitudinal fluctuations and is the analogous of the so-called principle of conservation of the modulus in an equilibrium ferromagnet [16]. From Eq. (16) we finally have and We are then left with the problem of determining the fluctuations of the transverse velocity v ⊥ in the presence of a non-zero field h. We will do so by analyzing the equations of motion for v ⊥ and δρ, which follow from inserting (18) into the velocity equation of motion (5) projected transverse to the direction of mean motion, and into the density equation of motion (4), and expanding in fields and derivatives. Again, details are relegated to Appendix A; the result is: and 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 appearing in those terms as {µ While the exact forms of Eqs. (22) are not important for what follows, for completeness we also give them in Appendix A.

Renormalization Group
We have shown so far that the response δΦ is determined by the global average of transversal fluctuations, Eq. (21). To compute this quantity, we proceed by a dynamical renormalization group (DRG) analysis [23] of Eqs. (22). Once again, this standard analysis is almost identical to that carried out in [12] for the zero field case. We start by averaging the equations of motion over the short-wavelength fluctuations: i.e., those with support in the "shell" of Fourier space b −1 Λ ≤ |q ⊥ | ≤ Λ, where Λ is an "ultra-violet cutoff", and b is an arbitrary rescaling factor. Then, one rescales lengths, time, δρ and v ⊥ in equations (22) The scaling exponents α, ζ, and z, known respectively as the "roughness", "anisotropy", and "dynamical" exponents, are at this point arbitrary.
This DRG process leads to a new, "renormalized" pair of equations of motion of the same form as (22), but with "renormalized" values of the parameters, {µ For a suitable choice of the scaling exponents α, ζ, and z, these parameters flow to fixed, finite limits as b → ∞; that is, {µ this is referred to as a "renormalization group 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 only generate 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 in the 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 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 original equations 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 α used for v ⊥ . However, since fluctuations of δρ and v ⊥ have the same scaling with distance and time, they prove 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 χ.
a simple scaling law. The example of interest for our problem is of course the correlation function Here L ⊥ and L are respectively the transverse and longitudinal system size. The DRG scaling law obeyed by C is thus which follows simply from the fact that C involves two powers of v ⊥ , each of which gives a factor b α . In order to examine the scaling of C with field amplitude h, we use the completely standard [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 to have the same value regardless of the bare value of h v . This implies that b = hv Note that for small h -and thus small h v -this choice implies b ≫ 1, and that the parameters {µ Hence, in the limit of small h, the scaling function (26) can be reduced to where 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 {µ Hence, we expect this scaling law to break down for large fields, and, in fact, it does.
We now focus our attention on roughly square systems, with L ⊥ ∼ L ∼ L. Assuming an anisotropy exponent 0 < ζ < 1 (as expected [6,7,8,9,12] for spatial dimensions d < 4), we have for small fields so that finite-size scaling is controlled by the transverse flock extension L ⊥ and we can replace g with the universal scaling function w(x) ≡ g(x, ∞). Above the upper critical dimension d c = 4, where ζ = 1 [6], scaling is the same in both the transversal and longitudinal directions and we choose instead w(x) ≡ g(x, x). Doing so, we finally obtain the scaling law where It is now straightforward to relate the order parameter change δΦ(h) to this scaling law. From Eq. (21) we have Note that for ν > 0 (again a condition we expect, as discussed below, to be satisfied for d < 4), we have C(L, h) ≫ O(h) and the corrections due to fluctuations dominate the mean field ones. To lowest order in h In the thermodynamic limit, Lh 1/z ≫ 1 for any non-zero field and w(Lh 1/z ) → w(∞) = constant, yielding the asymptotic result δΦ h ∝ h −ν . In practice, in any system large enough, the external field suppresses transverse fluctuations, thus increasing the scalar order 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 other hand, it is known [6] that α = 1−d/2 and z = 2, which implies ν = 2−d/2 ≤ 0; therefore corrections due to fluctuations no longer dominate the mean field ones. Hence, ordinary linear response χ → constant as h → 0 is recovered for d > 4.
So far, we have kept our discussion of DRG at a qualitative level, independent of the precise form of the zero-field terms in Eqs. (22). To be more quantitative for d < 4, we need the actual values of the scaling exponents α, ζ, and z for which the DRG flows to a fixed point in those dimensions. These values actually do depend on the form of Eqs. (22) (and, in particular, to the nature of their relevant nonlinear terms), but luckily for small fields h they have to coincide with their zero-field values. Indeed, there is no reason for which these zero-field values should be affected by a sufficiently small rescaled field, h v ≤ h * v . In [6], it was argued that for any dimension 3/2 ≤ d ≤ 4 these zero-field exponents are It has since been since realized [12] that the original arguments leading to these values are flawed. However, the simple conjecture that the only relevant non-linearity at the fixed point is the term proportional to λ 1 in Eqs. (5)-(6) leads to precisely these values for 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 values for d = 2 and, to a lesser extent, d = 3. In the following we will assume this conjecture holds, verifying it a-posteriori by numerical measures of asymptotic response in the Vicsek [2] model. Above the upper critical dimension, d c = 4, finally, the scaling exponents take the exact linear values z = 2, ζ = 1 and α = 1 − d/2.

Finite size effects and longitudinal response
We conclude this section discussing finite size effects. The scaling form (33) implies that transverse fluctuations are suppressed by the field h on length scales In small systems such that 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 parameter should revert to linear order in the external field h. We can include this behavior in a single universal scaling function f by requiring that f (∞) ∝ w(∞) = O(1) and f (x) ∝ x νz for x ≪ 1. This finally gives, with γ = νz. This scaling holds for external fields not too large. For h > h * v , on the other hand, the small field approximation discussed here is no longer valid, and saturation effects change the scaling (36). Once expressed in terms of the longitudinal susceptibility χ = δΦ/h, our results imply Eq. (2), with the scaling exponents given by Eq. (3) (according to conjecture (34)).

Numerical simulations
We test our predictions (36) in two and three dimensions by simulating the well known Vicsek model [2] in an external homogeneous field h. Each particle -labeled by i = 1, 2, . . . , N -is defined by a position r t i and a unit direction of motion v t i . The model evolves with a synchronous discrete time dynamics where v m is the particle speed , ϑ(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 angle of amplitude 4πω (d = 3). The interaction is "metric": that is, each particle i interacts with all of its neighbors within unit distance. In the following, we adopt periodic boundary conditions and choose typical microscopic parameters so that the system lies within the TT phase [5]: v m = 0.5, ρ 0 = N/L d = 1 and ω = 0.18 (d = 2) or ω = 0.11 (d = 3). In both dimensions, this choice yields a zero-field order parameter Φ(0) ≈ 0.8. We perform simulations with different external field amplitudes and with different linear system sizes L. After discarding a transient T 0 sufficiently long for the system to settle into the stationary state, we estimate the mean global order parameter and its standard error, given by S E = σ/ √ n, with σ being the standard deviation and n the number of independent data points. We estimate n as the total number of Note that v m = 0 is the equilibrium limit of this model [26], which is singular. stationary points T divided by the autocorrelation time τ of the mean order parameter timeseries, n = T /τ . In Eq. (39), · t denotes time averages, performed over typically T = 10 6 ∼ 10 8 time steps. In particular, as the precision of the zero-field order parameter affects all response and the autocorrelation time decreases with h, in our numerical simulations we take care to estimate the zero-field order parameter over times as large as possible.
We begin addressing response in the large system size (or large field) regime hL z ≫ 1, where our theory predicts δΦ ∼ h 1−ν . Measuring this power law is a particularly difficult task, as it is sandwiched between saturation effects at larger values of 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 effective power law exponent 1 − ν ef f by linear regression. The resulting response δΦ(h) is plotted in Fig. 1 for increasing system sizes L. As one expects, as system sizes increases, response curves approach the expected size-asymptotic behavior δΦ ∼ h 1−ν . In Fig. 1a, for instance, the d = 2 response approaches the expected power-law 1 − ν = 1/3. In the inset of Fig. 1a we further quantify this convergence plotting |ν ef f − ν| vs. the system size L. This shows that the effective exponent approaches the predicted one with corrections of order 1/ √ L. We repeated the same procedure in d = 3. As shown in Fig. 1b, the approach to the expected asymptotic exponent 1 − ν = 3/4 is faster, and the difference |ν ef f − ν| vanishes faster, as L −1.5 . In d = 3 our simulations are obviously limited to a much smaller range of linear size values, but it should be noted that in d = 3 finite size effects vanish quicker (being the exponent z larger) while the asymptotic exponent 1 − ν = 3/4 is already quite close to the value δΦ ∼ h expected at low values of h. A faster approach of ν ef f to its asymptotic value is therefore not completely surprising.
In d = 3, we can also easily compare the response behavior with the equilibrium prediction δΦ ∼ √ h [16] (dashed blue line in Fig. 1b). This clearly shows that the far-from-equilibrium nature of the Vicsek model makes the susceptibility exponent very different from that in equilibrium ferromagnets. In d = 2 equilibrium systems with a continuous symmetry cannot develop long range order, but, rather, exhibit only a quasilong range ordered phase, characterized by scaling exponents that vary continuously with temperature [13,27]. The equilibrium susceptibility exponent [27,28] in d = 2 is given by where the order parameter correlation exponent η is bounded: 0 ≤ η ≤ 1/4, which implies that the susceptibility exponent ν varies over an extremely narrow range: 14/15 ≤ ν ≤ 1. Our predicted value of ν = 2/3 for 2d flocks lies well outside this range; far enough, in fact, that our simulations both support the theory presented here, and rule out any equilibrium interpretation. Next, we consider the linear behavior δΦ ∝ h predicted for small system sizes (or small fields), hL z ≪ 1. This inequality imposes a (severe) upper limit on the range of h, while there is a lower limit set by our numerical precision in evaluating responses of the order of 10 −4 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 decade in the field amplitude h, especially for small system sizes.
By selecting a single h value lying in the linear regime for all accessible system sizes L, one is also able to test the saturation exponent gamma, δΦ ∼ h L γ . This is done in Figs. 2c-d, where response values at different linear sizes L are compared to the predicted power-law with (respectively) γ = 4/5 for d = 2 and γ = 2/5 in d = 3. We obtain a good agreement in d = 2 (the best linear fit being γ = 0.79 (5), while data in d = 3 is less clear, in rough agreement with the expected power-law behavior only for sizes L ≥ 60, with a best linear fit of γ = 0.48 (7).
We finally consider the full range of accessible external fields values in Fig. 3ab, which shows data for the accessible range of external field values in both two and three dimensions. Fields h larger than h s ≈ 0.1, of course, are out of the small field regime and show saturation effects, while due to statistical fluctuations we have been unable to obtain reliable estimates for external fields smaller than h ≈ 10 −4 . Within this range, comparison with the predicted scaling (2) (as given by the dashed lines) is overall satisfactory, especially in d = 2. By a proper rescaling, making use of the three scaling exponents (3), we can also collapse our data at different sizes on roughly a single  curve, as shown in Fig. 3c-d. To summarize, numerical simulations are in good agreement with our theoretical predictions, at least in d = 2. Results in d = 3 are prone to larger errors and obviously explore a more limited range of linear sizes, but nevertheless are still compatible with our 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 of the 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 the microscopic Vicsek equations (38). In Ref. [5] it was argued that different microscopic implementations could lead to different response, and in particular it was recommended to choose one by which the external field was normalized by the local order parameter  value, such as in However, we do not expect these microscopic details to change the structure of the hydrodynamic equations (4)- (5), and thus we do not expect qualitative differences between the two microscopic external field implementations. Indeed, our preliminary simulations of Eq. (41) (not shown here) show no qualitative difference from the response extensively discussed above for Eq. (38).

Conclusions
So far, we have only considered the longitudinal susceptibility, which characterizes the response of the magnitude of the order parameter to a small external field. Simple considerations, based on symmetry, imply that the flock polarization vector will eventually align with any non-zero stationary external field, including one applied transversal to the initial polarization. Thus, for the transvere susceptibility we have, trivially, as in equilibrium systems.
In this paper, we have fully characterized the static response of homogeneous ordered flocks to small external fields for any dimension d > 3/2. In particular, below the 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. This is ultimately a consequence of the spontaneous symmetry breaking of the continuous rotation symmetry, albeit the far-from equilibrium nature of flocks yields different results from, say, equilibrium ferromagnets in d = 2, 3 [13]. We have also fully characterized finite size effects -typically of great importance in biological applications of collective motion -and verified our results via numerical simulations. We believe that the finite numerical values reported in Ref. [3] for the longitudinal susceptibility are entirely due to finite size effects.
Incidentally, our numerics thereby also provide further evidence supporting the conjecture (34) for the scaling exponent values [6].
Our results are expected to hold generically for all collective motion systems showing a bona fide TT phase. This class encompasses both systems with metric interactions and those with topological interactions. It also includes the inertial spin model recently put forward in [30] to account for the turning dynamics measured experimentally in starling flocks. This is because the long time hydrodynamic theory of the inertial spin theory relaxes to the TT theory [31]; hence, the static response will be unchanged. Dynamical response (i.e. how quickly the flock turns towards the field direction), however, could be different at short times in inertial spin models, while the long-time behavior should be the same.
In future work, we will explore more thoroughly the phase diagram of Vicsek-like models beyond the TT phase, investigating the disordered and phase separated regimes Other future directions include the study of the finite-time, dynamical reponse [32] in both overdamped Vicsek-like models and inertial spin ones, and the study of spatially localized perturbations. fluctuations We first demonstrate that the longitudinal velocity δv is enslaved to the slow modes δρ and v ⊥ . We follow Ref. [12] and begin with the hydrodynamic Eqs. (4) and (5) written out explicitly: where, as noted in the main text, the parameters λ i (i = 1 → 3), the local term U, the "isotropic Pressure" P (ρ, |v|) and the "anisotropic Pressure"P 2 (ρ, |v|) are, in general, functions of the density ρ and the magnitude |v| of the local velocity. It is useful to Taylor expand P 1 and P 2 around the equilibrium density ρ 0 : Here D 1 , D 2 and D 3 are all positive in the ordered state. Again as discussed in the main text, in the ordered phase, the velocity field can be written as: (for simplicity, here and hereafter, we write v 0 ≡ v 0 (h)).
Taking the dot product of both sides of equation (A.2) with v itself, we obtain: In this hydrodynamic approach, we are interested only in fluctuations δ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 derivatives of δv(r, t) and δρ(r, t) are always negligible, in the hydrodynamic limit, compared to terms involving the same number of powers of fields without any time or space derivatives.
Furthermore, the fluctuations δv(r, t) and δρ(r, t) can themselves be shown to be small 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 in the v ⊥ equation of motion, which is negligible (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, we need only keep terms involving h that are proportional to h and independent of the small 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: where we've defined We can now express the longitudinal velocity δv in terms of the slow modes using equation (A.7) and the expansion where we've defined with, here and hereafter, super-or sub-scripts 0 denoting functions of ρ and |v| evaluated at ρ = ρ 0 and |v| = v 0 . We've also used the expansion (A.5) for the velocity in terms of the fluctuations δv and v ⊥ to write and kept only terms that an DRG analysis shows to be relevant in the long wavelength limit [12]. Inserting (A.9) into (A.7) gives: where we've kept only linear terms on the right hand side of this equation, since the non-linear terms are at least of order derivatives of |v ⊥ | 2 , and hence negligible, in the hydrodynamic limit, relative to the |v ⊥ | 2 term explicitly displayed on the left-hand side. This equation can be solved iteratively for δv in terms of v ⊥ , δρ, and its derivatives. To lowest (zeroth) order in derivatives, δv ≈ µ 1 δρ + h v 0 Γ 1 where we have defined Inserting this approximate expression for δv into equation (A.12) everywhere δv appears on the right hand side of that equation gives δv to first order in derivatives: where we've defined . and expresses the enslaving of δv to the slow modes δρ and v ⊥ .
We finally derive the equations of motion for the slow modes. Inserting the expression (A.7) for U back into equation (A.2), we find that P 2 and λ 2 cancel out of the v equation of motion, leaving 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 , and ii) eliminating δv 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: where we've defined h v as in the main text, and where we've defined: and, last but by no means least, The parameter D ρ ⊥ is actually zero at this point in the calculation, but we've included it in equation (A.29) anyway, because it is generated by the nonlinear terms under the Renormalization Group. Likewise, the parameter w 1 = 1, but will change from that value upon renormalization.
The equation of motion (A.29) is, as claimed in the main text, exactly the same as that in the absence of the external field h, while the equation of motion (A.20) is of the form (22), with