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

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–5] and hydrodynamic continuum theories [6–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 fluctuation–dissipation 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–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

$$\begin{eqnarray}&&{\chi }_{\parallel }\equiv \displaystyle \frac{\delta {\rm{\Phi }}(h)}{h},\end{eqnarray} \tag{ 1 }$$

where $\delta {\rm{\Phi }}(h)={\rm{\Phi }}(h)-{\rm{\Phi }}(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:

$$\begin{eqnarray}{\chi }_{{}_{\parallel }}={h}^{-\nu }f\left({{Lh}}^{\displaystyle \frac{1}{z}}\right)\propto \left\{\begin{array}{ll}{h}^{-\nu }, & L\gg {L}_{{\rm{c}}}(h),\\ {L}^{\gamma }, & L\ll {L}_{{\rm{c}}}(h),\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where ${L}_{{\rm{c}}}(h)\propto {h}^{-1/z}$ and, using a conjecture first put forward in [6],

$$\begin{eqnarray}&&\nu =\displaystyle \frac{4-d}{d+1},\;\;\;\;\;z=\displaystyle \frac{2(d+1)}{5},\;\;\;\;\;\gamma =\displaystyle \frac{2(4-d)}{5}\end{eqnarray} \tag{ 3 }$$

for any dimension $3/2\leqslant d\leqslant 4$, the upper critical dimension. In particular, we have $\nu =2/3$, $z=6/5$ and $\gamma =4/5$ in d = 2 and $\nu =1/4$, $z=8/5$ and $\gamma =2/5$ in d = 3. For $d\geqslant 4$, on the other hand, we predict $\delta {\rm{\Phi }}\propto 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.

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.

2.1. Hydrodynamic description

The hydrodynamic theory describes flocking by continuous, coarse grained number density $\rho ({\bf{r}},t)$ and velocity ${\bf{v}}({\bf{r}},t)$ fields. The hydrodynamic equations of motion governing these fields in the long-wavelength limit can be obtained either by symmetry arguments [6–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

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot ({\bf{v}}\rho )=0\end{eqnarray} \tag{ 4 }$$

and, in a schematic notation

$$\begin{eqnarray}&&{\partial }_{t}{\bf{v}}+{\boldsymbol{\Lambda }}[{\boldsymbol{\nabla }}{\bf{vv}}]=U(\rho ,| {\bf{v}}| ){\bf{v}}+{\bf{D}}[{\boldsymbol{\nabla }}{\boldsymbol{\nabla }}{\bf{v}}]+{{\bf{F}}}_{{\bf{P}}}+{\bf{f}}+{\bf{h}},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{\Lambda }}[{\boldsymbol{\nabla }}{\bf{vv}}]\equiv {\lambda }_{1}({\bf{v}}\cdot {\boldsymbol{\nabla }}){\bf{v}}+{\lambda }_{2}({\boldsymbol{\nabla }}\cdot {\bf{v}}){\bf{v}}+{\lambda }_{3}{\boldsymbol{\nabla }}(| {\bf{v}}{| }^{2})\end{eqnarray} \tag{ 6 }$$

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 ${\lambda }_{1}=1$ and ${\lambda }_{2}={\lambda }_{3}=0$, as in the usual Navier–Stokes equations.

The viscous terms

$$\begin{eqnarray}&&{\bf{D}}[{\boldsymbol{\nabla }}{\boldsymbol{\nabla }}{\bf{v}}]\equiv {D}_{1}{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\bf{v}})+{D}_{2}{({\bf{v}}\cdot {\boldsymbol{\nabla }})}^{2}{\bf{v}}+{D}_{3}{{\rm{\nabla }}}^{2}{\bf{v}}\end{eqnarray} \tag{ 7 }$$

reflect the tendency of localized fluctuations in the velocities to spread out because of local interactions.

The pressure term

$$\begin{eqnarray}&&{{\bf{F}}}_{{\bf{P}}}\equiv -{\boldsymbol{\nabla }}{P}_{1}-{\bf{v}}({\bf{v}}\cdot {\boldsymbol{\nabla }}{P}_{2})\end{eqnarray} \tag{ 8 }$$

is the sum of 'isotropic' (P₁) and 'anisotropic' (P₂) pressure terms, the latter being a genuinely non-equilibrium feature. Both terms tend to suppress local density fluctuations around the global mean value ${\rho }_{0}$. The pressures ${P}_{\mathrm{1,2}}$, and the convective and viscous parameters ${\lambda }_{k}$ and ${D}_{k}\gt 0$ ($k=1,2,3$) are functions of the local density ρ and the magnitude $| {\bf{v}}|$ of the local velocity.

Fluctuations are introduced through a Gaussian white noise ${\bf{f}}$ with correlations

$$\begin{eqnarray}&&\langle {f}_{i}({\bf{r}},t){f}_{j}({\bf{r}}^{\prime} ,t^{\prime} )\rangle ={\rm{\Delta }}{\delta }_{{ij}}{\delta }^{d}({\bf{r}}-{\bf{r}}^{\prime} )\delta (t-t^{\prime} ).\end{eqnarray} \tag{ 9 }$$

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 ${\bf{v}}$ have a nonzero magnitude ${v}_{0}(h)$ in the ordered phase. It satisfies the condition $U\gt 0$ for $| {\bf{v}}| \lt {v}_{0}(h=0)$, U = 0 for $| {\bf{v}}| ={v}_{0}(0)$, and $U\lt 0$ for $| {\bf{v}}| \gt {v}_{0}(0)$. This term thereby causes the flock to spontaneously break 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.

Equations (4) and (5) are identical to the unperturbed ones discussed in [12], except for the explicit addition of the coarse-grained constant field ${\bf{h}}$ in equation (5). By analyticity and rotational invariance, this field is linearly and isotropically proportional to the applied microscopic field when those fields are sufficiently small.

2.2. Mean-field analysis

We first discuss the system in the absence of fluctuations. Equations (4) and (5) admit a spatially uniform steady state solution

$$\begin{eqnarray}\rho ({\bf{r}},t) & = & {\rho }_{0}\\ {\bf{v}}({\bf{r}},t) & = & {{\bf{v}}}_{0}({\bf{h}}).\end{eqnarray} \tag{ 10 }$$

For any nonzero external field let ${\hat{{\bf{e}}}}_{\parallel }$ be the unit vector along ${\bf{h}}\equiv h\;{\hat{{\bf{e}}}}_{\parallel }$, while for strictly zero field ${\hat{{\bf{e}}}}_{\parallel }$ will be the direction of the spontaneous symmetry breaking. We have ${{\bf{v}}}_{0}({\bf{h}})={v}_{0}(h){\hat{{\bf{e}}}}_{\parallel }$, with the magnitude ${v}_{0}(h)$ of the homogeneous velocity ${{\bf{v}}}_{0}({\bf{h}})$ determined by the condition

$$\begin{eqnarray}&&U({v}_{0}(h)),{\rho }_{0}){v}_{0}(h)+h=0.\end{eqnarray} \tag{ 11 }$$

Since U is analytic in v, we have for small fields

$$\begin{eqnarray}&&{v}_{0}(h)-{v}_{0}(0)\propto h,\end{eqnarray} \tag{ 12 }$$

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

$$\begin{eqnarray}&&{\rm{\Phi }}(h)\equiv | \langle {\bf{v}}({\bf{r}},t)\rangle | \end{eqnarray} \tag{ 13 }$$

(here and hereafter $\langle \cdot \rangle$ denotes a global average in space and time) responds linearly in h.

2.3. 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 have just computed at mean field level.

In order to do so, we allow for small fluctuations around the homogeneous solution

$$\begin{eqnarray}\rho ({\bf{r}},t) & = & {\rho }_{0}+\delta \rho (h;{\bf{r}},t)\\ {\bf{v}}({\bf{r}},t) & = & {{\bf{v}}}_{0}(h)+\delta {\bf{v}}(h;{\bf{r}},t)\end{eqnarray} \tag{ 14 }$$

and distinguish between longitudinal and transverse velocity fluctuations, which are respectively parallel ($\parallel$) and perpendicular ($\perp$) to ${{\bf{v}}}_{0}(h)$

$$\begin{eqnarray}&&\delta {\bf{v}}(h;{\bf{r}},t)=\delta {v}_{\parallel }(h;{\bf{r}},t){{\bf{e}}}_{\parallel }+{{\bf{v}}}_{\perp }(h;{\bf{r}},t),\end{eqnarray} \tag{ 15 }$$

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, $\langle \delta \rho \rangle =0$, while symmetry considerations imply $\langle {{\bf{v}}}_{\perp }\rangle ={\bf{0}};$ that is, fluctuations can not steer the global average of $\langle {\bf{v}}({\bf{r}},t)\rangle$ away from the external field direction. This implies that corrections to the order parameter are linear in the longitudinal fluctuations: by making use of equations (12), (13) and (15) we have

$$\begin{eqnarray}&&{\rm{\Phi }}(h)={v}_{0}(h)+\langle \delta {\bf{v}}(h)\rangle ={v}_{0}(0)+\langle \delta {v}_{\parallel }(h)\rangle +O(h).\end{eqnarray} \tag{ 16 }$$

In order to compute longitudinal fluctuations, we have to expand the hydrodynamic equations (4) and (5) in the small fluctuations $\delta \rho$, $\delta {v}_{\parallel }$ and ${{\bf{v}}}_{\perp }$. 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 the appendix, but (fortunately) they are identical to order h with those for the zero-field case in [12]. The only difference at O(h) is that

$$\begin{eqnarray}&&U({v}_{0}(h),{\rho }_{0})\approx -\displaystyle \frac{h}{{v}_{0}}\ne 0.\end{eqnarray} \tag{ 17 }$$

Because slow modes dominate the long-distance behavior, and, it therefore proves, the small field response, we can eliminate the longitudinal fluctuations from equation (22), since they are a fast mode of the dynamics. The subtle details of this elimination are given in the appendix; 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

$$\begin{eqnarray}&&\delta {v}_{\parallel }\approx -\displaystyle \frac{| {{\bf{v}}}_{\perp }(h){| }^{2}}{2{v}_{0}(0)}+{\mu }_{1}\delta \rho +({\mu }_{2}{\partial }_{t}+{\mu }_{3}{\partial }_{\parallel })\;\delta \rho +{\mu }_{4}{{\rm{\nabla }}}_{\perp }\cdot {{\bf{v}}}_{\perp }+O(h).\end{eqnarray} \tag{ 18 }$$

Here ${\mu }_{1}$ is a constant which depends on the form of U. In typical flocking models with metric interactions ${\mu }_{1}\gt 0$ [11], so that density fluctuations are positively correlated with longitudinal fluctuations at the local level. In equation (18), ${{\rm{\nabla }}}_{\perp }$ denotes spatial derivatives in the transverse directions, and the constants ${\mu }_{2}$, ${\mu }_{3}$ and ${\mu }_{4}$ depend on the original parameters of the hydrodynamic equations (4) and (5). Full details, together with the derivation of equation (18), can be found in the appendix, but the exact form of these constants is unimportant here. Since these derivative terms are linear in $\delta \rho$ and ${{\bf{v}}}_{\perp }$, they vanish once averaged over space and time, so that from equation (18) we have

$$\begin{eqnarray}&&\langle \delta {v}_{\parallel }\rangle \approx -\displaystyle \frac{\langle | {{\bf{v}}}_{\perp }(h){| }^{2}\rangle }{2{v}_{0}(0)}+O(h)\end{eqnarray} \tag{ 19 }$$

which links the global average of transversal and longitudinal fluctuations and is the analog of the so-called principle of conservation of the modulus in an equilibrium ferromagnet [16]. From equation (16) we finally have

$$\begin{eqnarray}&&{\rm{\Phi }}(h)\approx {v}_{0}(0)-\displaystyle \frac{\langle | {{\bf{v}}}_{\perp }(h){| }^{2}\rangle }{2{v}_{0}(0)}+O(h)\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&\delta {\rm{\Phi }}(h)\equiv {\rm{\Phi }}(h)-{\rm{\Phi }}(0)\approx \displaystyle \frac{\langle | {{\bf{v}}}_{\perp }(0){| }^{2}\rangle -\langle | {{\bf{v}}}_{\perp }(h){| }^{2}\rangle }{2{v}_{0}(0)}+O(h).\end{eqnarray} \tag{ 21 }$$

We are then left with the problem of determining the fluctuations of the transverse velocity ${{\bf{v}}}_{\perp }$ in the presence of a non-zero field h. We will do so by analyzing the equations of motion for ${{\bf{v}}}_{\perp }$ and $\delta \rho$, 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 the appendix; the result is:

$$\begin{eqnarray}{\partial }_{t}\delta \rho & = & {[{\partial }_{t}\delta \rho ]}_{h=0}\\ {\partial }_{t}{{\bf{v}}}_{\perp } & = & {[{\partial }_{t}{{\bf{v}}}_{\perp }]}_{h=0}-{h}_{v}{{\bf{v}}}_{\perp },\end{eqnarray} \tag{ 22 }$$

where we have introduced the rescaled field

$$\begin{eqnarray}&&{h}_{v}\equiv \displaystyle \frac{h}{{v}_{0}(0)}\end{eqnarray} \tag{ 23 }$$

and ${[{\partial }_{t}\delta \rho ]}_{h=0}$ and ${[{\partial }_{t}{{\bf{v}}}_{\perp }]}_{h=0}$ are the terms originally given by equations (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 $\{{\mu }_{i}^{(0)}\}$. While the exact forms of equation (22) are not important for what follows, for completeness we also give them in the appendix.

2.4. Renormalization group

We have shown so far that the response $\delta {\rm{\Phi }}$ is determined by the global average of transversal fluctuations, equation (21). To compute this quantity, we proceed by a dynamical renormalization group (DRG) analysis [23] of equation (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}{\rm{\Lambda }}\leqslant | {{\bf{q}}}_{{}_{\perp }}| \leqslant {\rm{\Lambda }}$, where Λ is an 'ultra-violet cutoff', and b is an arbitrary rescaling factor. Then, one rescales lengths, time, $\delta \rho$ and ${{\bf{v}}}_{{}_{\perp }}$ in equation (22) according to ${{\bf{v}}}_{{}_{\perp }}={b}^{\alpha }{{\bf{v}}}_{{}_{\perp }}^{\prime }$, $\delta \rho ={b}^{\alpha }\delta {\rho }^{\prime }$, ${{\bf{r}}}_{{}_{\perp }}=b{{\bf{r}}}_{\perp }^{\prime }$, ${r}_{{}_{\parallel }}={b}^{\zeta }{r}_{{}_{\parallel }}^{\prime }$, and $t={b}^{z}{t}^{\prime }$ to restore the ultra-violet cutoff⁵ to Λ. 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, $\{{\mu }_{i}^{(0)}\}\to \{{\mu }_{i}^{(b)}\}$. For a suitable choice of the scaling exponents α, ζ, and z, these parameters flow to fixed, finite limits as $b\to \infty ;$ that is, $\{{\mu }_{i}^{(b\to \infty )}\}\to \{{\mu }_{i}^{* }\};$ 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 equation (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 equation (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

$$\begin{eqnarray}&&{h}_{v}={b}^{-z}{h}_{v}^{\prime },\end{eqnarray} \tag{ 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 original equations of motions can be related to those of the renormalized equations of motion via a simple scaling law. The example of interest for our problem is of course the correlation function

$$\begin{eqnarray}&&C({L}_{\perp },{L}_{\parallel },\{{\mu }_{i}^{0}\},{h}_{v})\equiv \langle | {{\bf{v}}}_{\perp }(0){| }^{2}\rangle -\langle | {{\bf{v}}}_{\perp }(h){| }^{2}\rangle .\end{eqnarray} \tag{ 25 }$$

Here L_⊥ and L_∥ are respectively the transverse and longitudinal system size. The DRG scaling law obeyed by C is thus

$$\begin{eqnarray}&&C({L}_{\perp },{L}_{\parallel },\{{\mu }_{i}^{0}\},{h}_{v})={b}^{2\alpha }C({b}^{-1}{L}_{\perp },{b}^{-\zeta }{L}_{\parallel },\{{\mu }_{i}^{b}\},{b}^{z}{h}_{v})\end{eqnarray} \tag{ 26 }$$

which follows simply from the fact that C involves two powers of ${{\bf{v}}}_{\perp }$, each of which gives a factor ${b}^{\alpha }$.

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={\left(\tfrac{{h}_{v}}{{h}_{v}^{* }}\right)}^{-1/z}$. Note that for small h—and thus small h_v—this choice implies $b\gg 1$, and that the parameters $\{{\mu }_{i}^{(b)}\}$ flow to $\{{\mu }_{i}^{* }\}$, their fixed point values. Hence, in the limit of small h, the scaling function (26) can be reduced to

$$\begin{eqnarray}&&C({L}_{\perp },{L}_{\parallel },\{{\mu }_{i}^{0}\},h)={h}^{-2\alpha /z}g({L}_{\perp }{h}^{1/z},{L}_{\parallel }{h}^{\zeta /z}),\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&g(x,y)\equiv {b}_{0}^{2\alpha /z}C({b}_{0}^{-1/z}x,{b}_{0}^{-\zeta /z}y,\{{\mu }_{i}^{* }\},{h}_{v}^{* })\;\;\;\;{\rm{with}}\;{b}_{0}\equiv {h}_{v}^{* }\;{v}_{0}(0),\end{eqnarray} \tag{ 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\to \infty$, and, hence $\{{\mu }_{i}^{(b)}\}\to \{{\mu }_{i}^{* }\}$. 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}_{\perp }\sim {L}_{\parallel }\sim L$. Assuming an anisotropy exponent $0\lt \zeta \lt 1$ (as expected [6–9, 12] for spatial dimensions $d\lt 4$), we have for small fields

$$\begin{eqnarray}&&{{Lh}}^{1/z}\propto {L}_{\perp }{h}^{1/z}\ll {L}_{\parallel }{h}^{\zeta /z}\end{eqnarray} \tag{ 29 }$$

so that finite-size scaling is controlled by the transverse flock extension ${L}_{\perp }$ and we can replace g with the universal scaling function $w(x)\equiv g(x,\infty )$. Above the upper critical dimension d_c = 4, where $\zeta =1$ [6], the transverse and longitudinal directions scale identically and we choose instead $w(x)\equiv g(x,x)$. Doing so, we finally obtain the scaling law

$$\begin{eqnarray}&&C(L,h)={h}^{1-\nu }\tilde{w}({{Lh}}^{1/z}),\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}&&\nu =1+2\alpha /z.\end{eqnarray} \tag{ 31 }$$

It is now straightforward to relate the order parameter change $\delta {\rm{\Phi }}(h)$ to this scaling law. From equation (21) we have

$$\begin{eqnarray}&&\delta {\rm{\Phi }}\propto C(L,h)+O(h).\end{eqnarray} \tag{ 32 }$$

Note that for $\nu \gt 0$ (again a condition we expect, as discussed below, to be satisfied for $d\lt 4$), we have $C(L,h)\gg O(h)$ and the corrections due to fluctuations dominate the mean field ones. To lowest order in h

$$\begin{eqnarray}&&\delta {\rm{\Phi }}=\displaystyle \frac{1}{2{v}_{0}(0)}[\langle | {{\bf{v}}}_{\perp }(0){| }^{2}\rangle -\langle | {{\bf{v}}}_{\perp }(h){| }^{2}\rangle ]={h}^{1-\nu }w({{Lh}}^{1/z}).\end{eqnarray} \tag{ 33 }$$

In the thermodynamic limit, ${{Lh}}^{1/z}\gg 1$ for any non-zero field and $w({{Lh}}^{1/z})\to w(\infty )={\rm{constant}}$, yielding the asymptotic result $\tfrac{\delta {\rm{\Phi }}}{h}\propto {h}^{-\nu }$. In practice, in any system large enough, the external field suppresses transverse fluctuations, thus increasing the scalar order parameter according to equation (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 $\alpha =1-d/2$ and z = 2, which implies $\nu =2-d/2\leqslant 0;$ therefore corrections due to fluctuations no longer dominate the mean field ones. Hence, ordinary linear response ${\chi }_{\parallel }\to {\rm{constant}}$ as $h\to 0$ is recovered for $d\gt 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 equation (22). To be more quantitative for $d\lt 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 equation (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}\leqslant {h}_{v}^{* }$.

In [6], it was argued that for any dimension $3/2\leqslant d\leqslant 4$ these zero-field exponents are

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{3-2d}{5},\;\;\;\;z=\displaystyle \frac{2(d+1)}{5},\;\;\;\;\ \zeta =\displaystyle \frac{d+1}{5}.\end{eqnarray} \tag{ 34 }$$

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 nonlineaarity at the fixed point is the term proportional to ${\lambda }_{1}$ in equations (5) and (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, $\zeta =1$ and $\alpha =1-d/2$.

2.5. Finite size effects and longitudinal response

We conclude this section the discussing finite size effects. The scaling form (33) implies that transverse fluctuations are suppressed by the field h on length scales

$$\begin{eqnarray}&&L\gg {L}_{{\rm{c}}}(h)\propto {h}^{-1/z}.\end{eqnarray} \tag{ 35 }$$

In small systems such that $L\ll {L}_{{\rm{c}}}$ (or equivalently for small field $h\ll {h}_{{\rm{c}}}(L)\propto {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(\infty )\propto w(\infty )=O(1)$ and $f(x)\propto {x}^{\nu z}$ for $x\ll 1$. This finally gives

$$\begin{eqnarray}\delta {\rm{\Phi }}={h}^{1-\nu }f(h\;{L}^{z})\propto \left\{\begin{array}{ll}{h}^{1-\nu }, & h\gg {L}^{-z},\\ {{hL}}^{\gamma }, & h\ll {L}^{-z},\end{array}\right.\end{eqnarray} \tag{ 36 }$$

with $\gamma =\nu z$. This scaling holds for external fields not too large. For $h\gt {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 ${\chi }_{\parallel }=\delta {\rm{\Phi }}/h$, our results imply equation (2), with the scaling exponents given by equation (3) (according to conjecture (34)).

We test our predictions (36) in two and three dimensions by simulating the well known Vicsek model [2] in an external homogeneous field ${\bf{h}}$. Each particle—labeled by $i=1,2,\ldots ,N$—is defined by a position ${{\bf{r}}}_{i}^{t}$ and a unit direction of motion ${{\bf{v}}}_{i}^{t}$. The model evolves with a synchronous discrete time dynamics

$$\begin{eqnarray}&&{{\bf{v}}}_{i}^{t+1}=({{ \mathcal R }}_{\omega }\;\circ \;\vartheta )\left(\displaystyle \sum _{j\in {S}_{i}}{{\bf{v}}}_{j}^{t}+{\bf{h}}\right),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{\bf{r}}}_{i}^{t+1}={{\bf{r}}}_{i}^{t}+{v}_{m}{{\bf{v}}}_{i}^{t},\end{eqnarray} \tag{ 38 }$$

where v_m is the particle speed⁷ , $\vartheta ({\bf{w}})={\bf{w}}/| {\bf{w}}|$ is a normalization operator and ${{ \mathcal R }}_{\omega }$ performs a random rotation (uncorrelated between different times t or particles i) uniformly distributed around the argument vector: ${{ \mathcal R }}_{\omega }{\bf{w}}$ is uniformly distributed around ${\bf{w}}$ inside an arc of amplitude $2\pi \omega$ (in d =  2) or in a spherical cap spanning a solid angle of amplitude $4\pi \omega$ (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, ${\rho }_{0}=N/{L}^{d}=1$ and $\omega =0.18$ (d =  2) or $\omega =0.11$ (d = 3). In both dimensions, this choice yields a zero-field order parameter ${\rm{\Phi }}(0)\approx 0.8$.

We perform simulations with different external field amplitudes and with different linear system sizes L. After discarding a transient T₀ sufficiently long for the system to settle into the stationary state, we estimate the mean global order parameter

$$\begin{eqnarray}&&{\rm{\Phi }}={\left\langle \displaystyle \frac{1}{N}\left|\displaystyle \sum _{i=1}^{N}{{\bf{v}}}_{i}^{t}\right|\right\rangle }_{t}\end{eqnarray} \tag{ 39 }$$

and its standard error, given by ${S}_{{\rm{E}}}=\sigma /\sqrt{n}$, with σ being the standard deviation and n the number of independent data points. We estimate n as the total number of stationary points T divided by the autocorrelation time τ of the mean order parameter timeseries, $n=T/\tau$. In equation (39), $\langle \cdot {\rangle }_{t}$ denotes time averages, performed over typically T = 10⁶ ∼ 10⁸ 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}\gg 1$, where our theory predicts $\delta {\rm{\Phi }}\sim {h}^{1-\nu }$. 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\ll {L}^{-z}$. We proceed by extrapolation, choosing a (somewhat arbitrary) h range of two decades and by measuring the effective power law exponent $1-{\nu }_{\mathrm{eff}}$ by linear regression. The resulting response $\delta {\rm{\Phi }}(h)$ is plotted in figure 1 for increasing system sizes L. As one expects, as system sizes increases, response curves approach the expected size-asymptotic behavior $\delta {\rm{\Phi }}\sim {h}^{1-\nu }$. In figure 1(a), for instance, the d = 2 response approaches the expected power-law $1-\nu =1/3$. In the inset of figure 1(a) we further quantify this convergence plotting $| {\nu }_{\mathrm{eff}}-\nu |$ versus the system size L. This shows that the effective exponent approaches the predicted one with corrections of order $1/\sqrt{L}$. We repeated the same procedure in d = 3. As shown in figure 1(b), the approach to the expected asymptotic exponent $1-\nu =3/4$ is faster, and the difference $| {\nu }_{\mathrm{eff}}-\nu |$ 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-\nu =3/4$ is already quite close to the value $\delta {\rm{\Phi }}\sim h$ expected at low values of h. A faster approach of ${\nu }_{\mathrm{eff}}$ to its asymptotic value is therefore not completely surprising.

Zoom In Zoom Out Reset image size

In d = 3, we can also easily compare the response behavior with the equilibrium prediction $\delta {\rm{\Phi }}\sim \sqrt{h}$ [16] (dashed blue line in figure 1(b)). 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 quasi-long 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

$$\begin{eqnarray}&&\nu =\displaystyle \frac{4-2\eta }{4-\eta },\end{eqnarray} \tag{ 40 }$$

where the order parameter correlation exponent η is bounded: $0\leqslant \eta \leqslant 1/4$, which implies that the susceptibility exponent ν varies over an extremely narrow range: $14/15\leqslant \nu \leqslant 1$. Our predicted value of $\nu =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 $\delta {\rm{\Phi }}\propto h$ predicted for small system sizes (or small fields), ${{hL}}^{z}\ll 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⁻⁴ or smaller. Nevertheless, our numerical simulations reveal in both d = 2 (figure 2(a)) and d = 3 (figure 2(b)) a linear growth of the response over more than one decade in the field amplitude h, especially for small system sizes.

Zoom In Zoom Out Reset image size

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 ,\delta {\rm{\Phi }}\sim h\;{L}^{\gamma }$. This is done in figures 2(c) and (d), where response values at different linear sizes L are compared to the predicted power-law with (respectively) $\gamma =4/5$ for d = 2 and $\gamma =2/5$ in d = 3. We obtain a good agreement in d = 2 (the best linear fit being $\gamma =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 $\gamma =0.48(7)$.

We finally consider the full range of accessible external fields values in figures 3(a) and (b), which shows data for the accessible range of external field values in both two and three dimensions. Fields h larger than ${h}_{s}\approx 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\approx {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 figures 3(c) and (d).

Zoom In Zoom Out Reset image size

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 equation (38). In [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

$$\begin{eqnarray}&&{{\bf{v}}}_{i}^{t+1}=({{ \mathcal R }}_{\omega }\;\circ \;\vartheta )\left[\vartheta \left(\displaystyle \sum _{j\in {S}_{i}}{{\bf{v}}}_{j}^{t}\right)+{\bf{h}}\right],\end{eqnarray} \tag{ 41 }$$

However, we do not expect these microscopic details to change the structure of the hydrodynamic equations (4) and (5), and thus we do not expect qualitative differences between the two microscopic external field implementations. Indeed, our preliminary simulations of equation (41) (not shown here) show no qualitative difference from the response extensively discussed above for equation (38).

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

$$\begin{eqnarray}&&{\chi }_{\perp }\sim \displaystyle \frac{{v}_{0}(0)}{h}\end{eqnarray} \tag{ 42 }$$

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\gt 3/2$. In particular, below the upper critical dimension d_c = 4, our results in the thermodynamic limit $L\to \infty$ show a diverging longitudinal response for $h\to 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 [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.

We have benefited from discussions with H Chaté and A Cavagna. We acknowledge support from the Marie Curie Career Integration Grant (CIG) PCIG13-GA-2013-618399. JT also acknowledges support from the SUPA distinguished visitor program and from the National Science Foundation through awards # EF-1137815 and 1006171, and thanks the University of Aberdeen for their hospitality while this work was underway. FG acknowledges support from EPSRC First Grant EP/K018450/1.

We first demonstrate that the longitudinal velocity $\delta {v}_{\parallel }$ is enslaved to the slow modes $\delta \rho$ and ${{\bf{v}}}_{\perp }$. We follow [12] and begin with the hydrodynamic equations (4) and (5) written out explicitly:

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot ({\bf{v}}\rho )=0,\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}{\partial }_{t}{\bf{v}} & = & -{\lambda }_{1}({\bf{v}}\cdot {\boldsymbol{\nabla }}){\bf{v}}-{\lambda }_{2}({\boldsymbol{\nabla }}\cdot {\bf{v}}){\bf{v}}-{\lambda }_{3}{\boldsymbol{\nabla }}(| {\bf{v}}{| }^{2})+U(\rho ,{\bf{v}}){\bf{v}}-{\boldsymbol{\nabla }}{P}_{1}\\ & & -{\bf{v}}({\bf{v}}\cdot {\boldsymbol{\nabla }}{P}_{2})+{D}_{1}{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\bf{v}})+{D}_{3}{{\rm{\nabla }}}^{2}{\bf{v}}+{D}_{2}{({\bf{v}}\cdot {\boldsymbol{\nabla }})}^{2}{\bf{v}}+{\bf{f}}+{\bf{h}},\end{eqnarray} \tag{ A.2 }$$

where, as noted in the main text, the parameters ${\lambda }_{i}(i=1\to 3)$, the local term U, the 'isotropic pressure' $P(\rho ,| {\bf{v}}| )$ and the 'anisotropic pressure' ${P}_{2}(\rho ,| {\bf{v}}| )$ are, in general, functions of the density ρ and the magnitude $| {\bf{v}}|$ of the local velocity. It is useful to Taylor expand P₁ and P₂ around the equilibrium density ${\rho }_{0}$:

$$\begin{eqnarray}&&{P}_{1}=\displaystyle \sum _{n=1}^{\infty }{\sigma }_{n}(| {\bf{v}}| ){(\rho -{\rho }_{0})}^{n},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{P}_{2}={P}_{2}(\rho ,| {\bf{v}}| )=\displaystyle \sum _{n=1}^{\infty }{{\rm{\Upsilon }}}_{n}(| {\bf{v}}| ){(\rho -{\rho }_{0})}^{n}.\end{eqnarray} \tag{ 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 be written as:

$$\begin{eqnarray}&&{\bf{v}}={v}_{0}{{\bf{e}}}_{{}_{\parallel }}+\delta {\bf{v}}=({v}_{0}+\delta {v}_{{}_{\parallel }}){{\bf{e}}}_{{}_{\parallel }}+{{\bf{v}}}_{\perp }\;\;\end{eqnarray} \tag{ A.5 }$$

(for simplicity, here and hereafter, we write ${v}_{0}\equiv {v}_{0}(h)$).

Taking the dot product of both sides of equation (A.2) with ${\bf{v}}$ itself, we obtain:

$$\begin{eqnarray}\displaystyle \frac{1}{2}({\partial }_{t}| {\bf{v}}{| }^{2}) & + & \displaystyle \frac{1}{2}({\lambda }_{1}+2{\lambda }_{3})({\bf{v}}\cdot {\boldsymbol{\nabla }})| {\bf{v}}{| }^{2})+{\lambda }_{2}({\boldsymbol{\nabla }}\cdot {\bf{v}})| {\bf{v}}{| }^{2}=U(| {\bf{v}})| {\bf{v}}{| }^{2}\\ & - & {\bf{v}}\cdot {\boldsymbol{\nabla }}{P}_{1}-| {\bf{v}}{| }^{2}{\bf{v}}\cdot {\boldsymbol{\nabla }}{P}_{2}+{D}_{1}{\bf{v}}\cdot {\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\bf{v}})+{D}_{3}{\bf{v}}\cdot {{\rm{\nabla }}}^{2}{\bf{v}}\\ & + & {D}_{2}{\bf{v}}\cdot ({({\bf{v}}\cdot {\boldsymbol{\nabla }})}^{2}{\bf{v}})+{\bf{v}}\cdot {\bf{f}}+{\bf{v}}\cdot {\bf{h}}.\end{eqnarray} \tag{ A.6 }$$

In this hydrodynamic approach, we are interested only in fluctuations $\delta {\bf{v}}({\bf{r}},t)$ and $\delta \rho ({\bf{r}},t)$ that vary slowly in space and time. (Indeed, the hydrodynamic equations (A.1) and (A.2) are only valid in this limit). Hence, terms involving space and time derivatives of $\delta {\bf{v}}({\bf{r}},t)$ and $\delta \rho ({\bf{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 $\delta {\bf{v}}({\bf{r}},t)$ and $\delta \rho ({\bf{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 $\delta {\bf{v}}({\bf{r}},t)$ and $\delta \rho ({\bf{r}},t)$. The ${\bf{v}}\cdot {\bf{f}}$ term can likewise be dropped, since it only leads to a term of order ${{\bf{v}}}_{{}_{\perp }}{f}_{{}_{\parallel }}$ in the ${{\bf{v}}}_{{}_{\perp }}$ equation of motion, which is negligible (since ${{\bf{v}}}_{{}_{\perp }}$ is small) relative to the ${{\bf{f}}}_{\perp }$ 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 $\delta {\bf{v}}$ and $\delta \rho$.

These observations can be used to eliminate many of the terms in equation (A.6), and solve for the quantity U; the solution is:

$$\begin{eqnarray}&&U=-\displaystyle \frac{h}{{v}_{0}}+{\lambda }_{2}{\boldsymbol{\nabla }}\cdot {\bf{v}}+{\bf{v}}\cdot {\boldsymbol{\nabla }}{P}_{2}+\displaystyle \frac{{\sigma }_{1}}{{v}_{0}}{\partial }_{{}_{\parallel }}\delta \rho +\displaystyle \frac{1}{2{v}_{0}}({\partial }_{t}+{\lambda }_{4}{\partial }_{{}_{\parallel }})\delta {v}_{{}_{\parallel }},\end{eqnarray} \tag{ A.7 }$$

where we have defined

$$\begin{eqnarray}&&{\lambda }_{4}\equiv ({\lambda }_{1}+2{\lambda }_{3}){v}_{0}\;.\end{eqnarray} \tag{ A.8 }$$

We can now express the longitudinal velocity $\delta {v}_{\parallel }$ in terms of the slow modes using equation (A.7) and the expansion

$$\begin{eqnarray}&&U\approx -{{\rm{\Gamma }}}_{1}\left(\delta {v}_{{}_{\parallel }}+\displaystyle \frac{| {{\bf{v}}}_{{}_{\perp }}{| }^{2}}{2{v}_{0}}\right)-{{\rm{\Gamma }}}_{2}\delta \rho \;\;,\end{eqnarray} \tag{ A.9 }$$

where we have defined

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1}\equiv -{\left(\displaystyle \frac{\partial U}{\partial | {\bf{v}}| }\right)}_{\rho }^{0},{{\rm{\Gamma }}}_{2}\equiv -{\left(\displaystyle \frac{\partial U}{\partial \rho }\right)}_{| {\bf{v}}| }^{0}\;,\end{eqnarray} \tag{ A.10 }$$

with, here and hereafter, super- or sub-scripts 0 denoting functions of ρ and $| {\bf{v}}|$ evaluated at $\rho ={\rho }_{0}$ and $| {\bf{v}}| ={v}_{0}$. We have also used the expansion (A.5) for the velocity in terms of the fluctuations $\delta {v}_{{}_{\parallel }}$ and ${\vec{v}}_{\perp }$ to write

$$\begin{eqnarray}&&| {\bf{v}}| ={v}_{0}+\delta {v}_{{}_{\parallel }}+\displaystyle \frac{| {{\bf{v}}}_{{}_{\perp }}{| }^{2}}{2{v}_{0}}+O(\delta {v}_{{}_{\parallel }}^{2},| {{\bf{v}}}_{{}_{\perp }}{| }^{4})\end{eqnarray} \tag{ A.11 }$$

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:

$$\begin{eqnarray}-{{\rm{\Gamma }}}_{1}\left(\delta {v}_{{}_{\parallel }}+\displaystyle \frac{| {{\bf{v}}}_{{}_{\perp }}{| }^{2}}{2{v}_{0}}\right)-{{\rm{\Gamma }}}_{2}\delta \rho & = & -\displaystyle \frac{h}{{v}_{0}}+{\lambda }_{2}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot {{\bf{v}}}_{{}_{\perp }}+{\lambda }_{2}{\partial }_{{}_{\parallel }}\delta {v}_{{}_{\parallel }}\\ & & +\displaystyle \frac{({{\rm{\Upsilon }}}_{1}{v}_{0}^{2}+{\sigma }_{1})}{{v}_{0}}{\partial }_{{}_{\parallel }}\delta \rho +\displaystyle \frac{1}{2{v}_{0}}({\partial }_{t}+{\lambda }_{4}{\partial }_{{}_{\parallel }})\delta {v}_{{}_{\parallel }},\end{eqnarray} \tag{ A.12 }$$

where we have kept only linear terms on the right-hand side of this equation, since the nonlinear terms are at least of order derivatives of $| {{\bf{v}}}_{{}_{\perp }}{| }^{2}$, and hence negligible, in the hydrodynamic limit, relative to the $| {{\bf{v}}}_{{}_{\perp }}{| }^{2}$ term explicitly displayed on the left-hand side.

This equation can be solved iteratively for $\delta {v}_{{}_{\parallel }}$ in terms of ${{\bf{v}}}_{{}_{\perp }}$, $\delta \rho$, and its derivatives. To lowest (zeroth) order in derivatives, $\delta {v}_{{}_{\parallel }}\approx {\mu }_{1}\delta \rho +\tfrac{h}{{v}_{0}{{\rm{\Gamma }}}_{1}}$, where we have defined

$$\begin{eqnarray}&&{\mu }_{1}=-\displaystyle \frac{{{\rm{\Gamma }}}_{2}}{{{\rm{\Gamma }}}_{1}}.\end{eqnarray} \tag{ A.13 }$$

Inserting this approximate expression for $\delta {v}_{{}_{\parallel }}$ into equation (A.12) everywhere $\delta {v}_{{}_{\parallel }}$ appears on the right-hand side of that equation gives $\delta {v}_{{}_{\parallel }}$ to first order in derivatives:

$$\begin{eqnarray}&&\delta {v}_{{}_{\parallel }}\approx -\displaystyle \frac{| {{\bf{v}}}_{{}_{\perp }}{| }^{2}}{2{v}_{0}}+{\mu }_{1}\;\delta \rho +\displaystyle \frac{{{\rm{\Gamma }}}_{2}}{{v}_{0}{{\rm{\Gamma }}}_{1}^{2}}{\partial }_{t}\delta \rho -\displaystyle \frac{{\lambda }_{5}}{{{\rm{\Gamma }}}_{1}}{\partial }_{{}_{\parallel }}\delta \rho -\displaystyle \frac{{\lambda }_{2}}{{{\rm{\Gamma }}}_{1}}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot {{\bf{v}}}_{{}_{\perp }}+\displaystyle \frac{h}{{v}_{0}{{\rm{\Gamma }}}_{1}}\;,\end{eqnarray} \tag{ A.14 }$$

where we have defined

$$\begin{eqnarray}&&{\lambda }_{5}\equiv \displaystyle \frac{({{\rm{\Upsilon }}}_{1}{v}_{0}^{2}+{\sigma }_{1})}{{v}_{0}}-\displaystyle \frac{{{\rm{\Gamma }}}_{2}}{{{\rm{\Gamma }}}_{1}}\left({\lambda }_{2}+\displaystyle \frac{{\lambda }_{4}}{{v}_{0}}\right)=\displaystyle \frac{({{\rm{\Upsilon }}}_{1}{v}_{0}^{2}+{\sigma }_{1})}{{v}_{0}}+{\mu }_{1}({\lambda }_{1}+{\lambda }_{2}+2{\lambda }_{3}).\end{eqnarray} \tag{ A.15 }$$

(In deriving the second equality in (A.15), we have used the definition (A.8) of ${\lambda }_{4}$.) Equation (A.14) coincides with equation (18) in the main text with

$$\begin{eqnarray}&&{\mu }_{2}=\displaystyle \frac{{{\rm{\Gamma }}}_{2}}{{v}_{0}{{\rm{\Gamma }}}_{1}^{2}},\end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray}&&{\mu }_{3}=-\displaystyle \frac{{\lambda }_{5}}{{{\rm{\Gamma }}}_{1}},\end{eqnarray} \tag{ A.17 }$$

$$\begin{eqnarray}&&{\mu }_{4}=-\displaystyle \frac{{\lambda }_{2}}{{{\rm{\Gamma }}}_{1}}\end{eqnarray} \tag{ A.18 }$$

and expresses the enslaving of $\delta {v}_{\parallel }$ to the slow modes $\delta \rho$ and ${{\bf{v}}}_{\perp }$.

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₂ and ${\lambda }_{2}$ cancel out of the ${\bf{v}}$ equation of motion, leaving

$$\begin{eqnarray}{\partial }_{t}{\bf{v}} & = & -{\lambda }_{1}({\bf{v}}\cdot {\boldsymbol{\nabla }}){\bf{v}}-{\lambda }_{3}{\boldsymbol{\nabla }}(| {\bf{v}}{| }^{2})-\displaystyle \frac{h}{{v}_{0}}{\bf{v}}+\displaystyle \frac{{\sigma }_{1}}{{v}_{0}}{\bf{v}}({\partial }_{{}_{\parallel }}\delta \rho )-{\boldsymbol{\nabla }}{P}_{1}+{D}_{1}{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\bf{v}})\\ & & +{D}_{3}{{\rm{\nabla }}}^{2}{\bf{v}}+{D}_{2}{({\bf{v}}\cdot {\boldsymbol{\nabla }})}^{2}{\bf{v}}+\left[\displaystyle \frac{1}{2{v}_{0}}({\partial }_{t}+{\lambda }_{4}{\partial }_{{}_{\parallel }})\delta {v}_{{}_{\parallel }}\right]{\bf{v}}+{\bf{f}}+{\bf{h}}.\end{eqnarray} \tag{ A.19 }$$

This can be made into an equation of motion for ${{\bf{v}}}_{{}_{\perp }}$ involving only ${{\bf{v}}}_{{}_{\perp }}({\bf{r}},t)$ and $\delta \rho ({\bf{r}},t)$ by (i) projecting (A.19) perpendicular to the direction of mean flock motion ${{\bf{e}}}_{{}_{\parallel }}$, and (ii) eliminating $\delta {v}_{{}_{\parallel }}$ by inserting equation (A.14) into the equation of motion (A.19) for ${\bf{v}}$. Using the expansions (A.5), (A.11) and neglecting 'irrelevant' terms we have:

$$\begin{eqnarray}{\partial }_{t}{{\bf{v}}}_{{}_{\perp }} & = & -{\lambda }_{1}^{0}{v}_{0}{\partial }_{{}_{\parallel }}{{\bf{v}}}_{{}_{\perp }}-{\lambda }_{1}^{0}({{\bf{v}}}_{{}_{\perp }}\cdot {{\boldsymbol{\nabla }}}_{{}_{\perp }}){{\bf{v}}}_{{}_{\perp }}-{g}_{1}\delta \rho {\partial }_{{}_{\parallel }}{{\bf{v}}}_{{}_{\perp }}-{g}_{2}{{\bf{v}}}_{{}_{\perp }}{\partial }_{{}_{\parallel }}\delta \rho -\displaystyle \frac{{c}_{0}^{2}}{{\rho }_{0}}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\delta \rho \\ & & -{g}_{3}{{\boldsymbol{\nabla }}}_{{}_{\perp }}(\delta {\rho }^{2})+{D}_{B}{{\boldsymbol{\nabla }}}_{{}_{\perp }}({{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot {{\bf{v}}}_{{}_{\perp }})+{D}_{3}{{\rm{\nabla }}}_{{}_{\perp }}^{2}{{\bf{v}}}_{{}_{\perp }}+{D}_{{}_{\parallel }}{\partial }_{{}_{\parallel }}^{2}{{\bf{v}}}_{{}_{\perp }}+{g}_{t}{\partial }_{t}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\delta \rho \\ & & +{g}_{{}_{\parallel }}{\partial }_{{}_{\parallel }}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\delta \rho +{{\bf{f}}}_{\perp }-{h}_{v}{{\bf{v}}}_{{}_{\perp }},\end{eqnarray} \tag{ A.20 }$$

where we have defined h_v as in the main text, and

$$\begin{eqnarray}&&{D}_{B}\equiv {D}_{1}+\displaystyle \frac{2{v}_{0}{\lambda }_{3}^{0}{\lambda }_{2}^{0}}{{{\rm{\Gamma }}}_{1}}\;,\end{eqnarray} \tag{ A.21 }$$

$$\begin{eqnarray}&&{D}_{{}_{\parallel }}\equiv {D}_{3}+{D}_{2}{v}_{0}^{2}\;,\end{eqnarray} \tag{ A.22 }$$

$$\begin{eqnarray}&&{g}_{1}\equiv {\mu }_{1}{\lambda }_{1}^{0}+{v}_{0}{\left(\displaystyle \frac{\partial {\lambda }_{1}}{\partial \rho }\right)}_{0}\;,\end{eqnarray} \tag{ A.23 }$$

$$\begin{eqnarray}&&{g}_{2}\equiv -{\mu }_{1}\displaystyle \frac{{\lambda }_{4}^{0}}{{v}_{0}}-\displaystyle \frac{{\sigma }_{1}^{0}}{{v}_{0}}\;,\end{eqnarray} \tag{ A.24 }$$

$$\begin{eqnarray}&&{g}_{3}\equiv {\sigma }_{2}^{0}+{\mu }_{1}^{2}{\lambda }_{3}^{0}+{v}_{0}{\mu }_{1}{\left(\displaystyle \frac{\partial {\lambda }_{3}}{\partial \rho }\right)}_{0}\;,\end{eqnarray} \tag{ A.25 }$$

$$\begin{eqnarray}&&{c}_{0}^{2}\equiv {\rho }_{0}{\sigma }_{1}^{0}+2{\rho }_{0}{v}_{0}{\lambda }_{3}^{0}{\mu }_{1}\;,\end{eqnarray} \tag{ A.26 }$$

$$\begin{eqnarray}&&{g}_{t}\equiv 2{\mu }_{1}{\lambda }_{3}^{0}\end{eqnarray} \tag{ A.27 }$$

and

$$\begin{eqnarray}&&{g}_{{}_{\parallel }}\equiv \displaystyle \frac{2{v}_{0}{\lambda }_{3}^{0}{\lambda }_{5}^{0}}{{{\rm{\Gamma }}}_{1}}-{\mu }_{1}{D}_{1}.\end{eqnarray} \tag{ A.28 }$$

Finally, using (A.5), (A.11) and (A.14) in the equation of motion (A.1) for ρ gives, again neglecting irrelevant terms:

$$\begin{eqnarray}{\partial }_{t}\delta \rho & = & -{\rho }_{0}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot {{\bf{v}}}_{{}_{\perp }}-{w}_{1}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot ({{\bf{v}}}_{{}_{\perp }}\delta \rho )-{v}_{2}{\partial }_{{}_{\parallel }}\delta \rho +{D}_{{\rho }_{\parallel }}{\partial }_{{}_{\parallel }}^{2}\delta \rho +{D}_{{\rho }_{\perp }}{{\rm{\nabla }}}_{{}_{\perp }}^{2}\delta \rho \\ & & +{D}_{\rho v}{\partial }_{{}_{\parallel }}({{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot {{\bf{v}}}_{{}_{\perp }})+{\rho }_{0}{\mu }_{2}{\partial }_{t}{\partial }_{{}_{\parallel }}\delta \rho -{\mu }_{1}{\partial }_{{}_{\parallel }}(\delta {\rho }^{2})+{\mu }_{5}{\partial }_{{}_{\parallel }}(| {{\bf{v}}}_{{}_{\perp }}{| }^{2}),\end{eqnarray} \tag{ A.29 }$$

where we have defined:

$$\begin{eqnarray}&&{v}_{2}\equiv {v}_{0}+{\mu }_{1}{\rho }_{0},\end{eqnarray} \tag{ A.30 }$$

$$\begin{eqnarray}&&{\mu }_{5}\equiv \displaystyle \frac{{\rho }_{0}}{2{v}_{0}}\;,\end{eqnarray} \tag{ A.31 }$$

$$\begin{eqnarray}&&{D}_{{\rho }_{\parallel }}\equiv \displaystyle \frac{{\rho }_{0}{\lambda }_{5}^{0}}{{{\rm{\Gamma }}}_{1}}=-{\rho }_{0}{\mu }_{3}\end{eqnarray} \tag{ A.32 }$$

and, last but by no means least

$$\begin{eqnarray}&&{D}_{\rho v}\equiv \displaystyle \frac{{\lambda }_{2}^{0}{\rho }_{0}}{{{\rm{\Gamma }}}_{1}}=-{\rho }_{0}\;{\mu }_{4}.\end{eqnarray} \tag{ A.33 }$$

The parameter ${D}_{{\rho }_{\perp }}$ is actually zero at this point in the calculation, but we have 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

$$\begin{eqnarray}{[{\partial }_{t}{{\bf{v}}}_{\perp }]}_{h=0} & \equiv & -{\lambda }_{1}^{0}{v}_{0}{\partial }_{{}_{\parallel }}{{\bf{v}}}_{{}_{\perp }}-{\lambda }_{1}^{0}({{\bf{v}}}_{{}_{\perp }}\cdot {{\boldsymbol{\nabla }}}_{{}_{\perp }}){{\bf{v}}}_{{}_{\perp }}-{g}_{1}\delta \rho {\partial }_{{}_{\parallel }}{{\bf{v}}}_{{}_{\perp }}-{g}_{2}{{\bf{v}}}_{{}_{\perp }}{\partial }_{{}_{\parallel }}\delta \rho \\ & & -\displaystyle \frac{{c}_{0}^{2}}{{\rho }_{0}}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\delta \rho -{g}_{3}{{\boldsymbol{\nabla }}}_{{}_{\perp }}(\delta {\rho }^{2})+{D}_{B}{{\boldsymbol{\nabla }}}_{{}_{\perp }}({{\boldsymbol{\nabla }}}_{{}_{\perp }}\cdot {{\bf{v}}}_{{}_{\perp }})+{D}_{3}{{\rm{\nabla }}}_{{}_{\perp }}^{2}{{\bf{v}}}_{{}_{\perp }}\\ & & +{D}_{{}_{\parallel }}{\partial }_{{}_{\parallel }}^{2}{{\bf{v}}}_{{}_{\perp }}+{g}_{t}{\partial }_{t}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\delta \rho +{g}_{{}_{\parallel }}{\partial }_{{}_{\parallel }}{{\boldsymbol{\nabla }}}_{{}_{\perp }}\delta \rho +{{\bf{f}}}_{\perp }.\end{eqnarray} \tag{ A.34 }$$