Hydrodynamics of Turning Flocks
HHydrodynamics of Turning Flocks
Xingbo Yang and M. Cristina Marchetti
1, 2 Physics Department, Syracuse University, Syracuse NY 13244, USA Syracuse Biomaterials Institute, Syracuse University, Syracuse NY 13244, USA
We present a hydrodynamic model of flocking that generalizes the familiar Toner-Tu equations to in-corporate turning inertia of well-polarized flocks. The continuum equations controlled by only two di-mensionless parameters, orientational inertia and alignment strength, are derived by coarse graining theinertial spin model recently proposed by Cavagna et al. [1]. The interplay between orientational inertiaand bend elasticity of the flock yields anisotropic spin waves that mediate the propagation of turning infor-mation throughout the flock. The coupling between spin current density to the local vorticity field througha nonlinear friction gives rise to a hydrodynamic mode with angular-dependent propagation speed at longwavelength. This mode goes unstable as a result of the growth of bend and splay deformations augmentedby the spin wave, signaling the transition to complex spatio-temporal patterns of continuously turning andswirling flocks.
PACS numbers: 87.18.Gh, 05.65.+b, 47.54.-r, 87.18.Hf
The Vicsek model [2, 3] and related continuous-timevariations [4] have been used to model flocking in a va-riety of systems, from birds [5] to cells [6] to in vitro cellular components [7, 8] and synthetic swimmers [9].These are examples of active systems, consisting of indi-vidually driven, dissipative units that exhibit coordinatedmotion (flocking) at large scales [10, 11]. In the Vicsekmodel the active units are described as point particles withoverdamped dynamics carrying a velocity vector of fixedmagnitude, hence “flying spins”. Each spin tends to alignwith its neighbors, but makes errors, modeled as angularnoise [3]. The system exhibits a liquid-gas phase transitionfrom a disordered gas state to a polar liquid state as thenoise is decreased or the number density is increased, withmicrophase separation in the coexistence region [12]. Theexistence of the transition has been put on firm groundsby a large number of numerical studies [3, 13]. Tonerand Tu also proposed a continuum version of the modelinspired by dynamical field theories of condensed mattersystems [14, 15].Recent work [16] has suggested that the description ofthe observed collective turning of bird flocks requires amodification of the Vicsek model to include angular inertiain the dynamics. This allows propagation of angular cor-relations through the flock on large scales via spin-wave-like excitations [1]. In this paper we derive the continuumequations for such an “inertial spin model” by explicitlycoarse-graining the microscopic dynamics. The resultingequations (Eqs. 6-8) generalize the Toner-Tu model to ac-count for turning modes by incorporating the dynamics ofthe spin angular momentum of the flock. These equations,governed by only two dimensionless parameters, are thefirst important result of our work. They contain new termsas compared to the phenomenological model of Ref. [17],most importantly a nonlinear friction that couples spin anddensity fluctuations to bend and splay deformations of theorder parameter. This new coupling transforms the prop-agating density bands ubiquitously observed in flocking models into turning bands of spin currents, driving the tran-sition to a novel state of continuously swirling and rotat-ing flocks, where turning information are transmitted byanisotropic propagating spin waves. The predicted soundspeeds could in principle be measured in experiments.Our starting point is the continuous-time model of iner-tial spins proposed by Cavagna et al. [1], where N pointparticles in a two-dimensional box of area L , with aver-age number density ρ = N/L , interact via a pairwisealigning interaction. Each particle is described by its posi-tion r i and the direction of its velocity, identified by a unitvector ˆ e θ i = (cos θ i , sin θ i ) in 2D. The dynamics of the i -th spin is described by d r i dt = v ˆ e θ i , dθ i dt = 1 χ s i , (1) ds i dt = γ (cid:88) j ˜ F ( θ j − θ i , r ji ) − ηχ s i + √ (cid:15) ξ i ( t ) , (2)with r ji = r j − r i , v the self-propulsion speed, s i thespin angular momentum and χ the spin moment of inertia.The spin is an internal angular momentum that generatesthe self-rotation, and is distinct from the angular momen-tum of the center of mass. The polar aligning interaction ofstrength γ is given by ˜ F ( θ, r ) = sin( θ ) / ( πR ) if | r | ≤ R and zero otherwise, with R the range of interaction. Thisform of the interaction used before in the literature [18]allows us to make analytical progress in the derivation ofthe continuum equations. Finally, η is a friction and (cid:15) de-scribes the strength of the angular noise, with ξ i ( t ) a Gaus-sian white noise with zero mean and unit variance.On time scales large compared to the relaxation time τ η = χ/η , one can neglect the time derivative on theleft hand side of Eq. (2) and eliminate the spin angularmomentum, s i , from the angular dynamics. This yieldsa continuous-time version of the Vicsek model, with ef-fective alignment strength γ/η and effective angular noise (cid:15)/η . Two additional time scales govern the dynamics of a r X i v : . [ c ond - m a t . s o f t ] J a n − − − − a b c d e f FIG. 1: Flocking patterns obtained via numerical solutions ofEqs. 6 -8 (top & middle rows) and via particle simulations usingEqs.1-2 (bottom row). The left/right columns show the patternsobtained at points
A/B of the phase diagram Fig.2a, correspond-ing to the spinodal portion of the coexistence region of travelingbands of polar liquid and disordered gas and to the region of thespin-wave instability, respectively. The snapshots (a-d) are ob-tained with ˜ γ = 2 . , ˜ χ = 0 . (a,c) and ˜ γ = 7 . , ˜ χ = 1 . (b,d) ona 300 by 300 grid lattice with grid size . , integration time step . and periodic boundary conditions. The particle simulationsare performed with particles in a box of size L = 10 withperiodic boundary conditions. Simulation parameters are R = 1 , (cid:15) = 2 . , v = 2 . , χ = 1 . , η = 1 . and γ = 0 . (e) and γ = 0 . (f). The integration time step is . . a-d: the arrowsrepresent the local polarization, with length proportional to thepolarization strength. The color indicates spin current density ina-b and number density in c-d. e-f: the arrows represent polariza-tion of individual particles (see Supplementary Movies). the system: the effective rotational diffusion time, τ (cid:15) = η /(cid:15) , and the alignment time, τ γ = η/ ( γρ ) .Following standard methods [19, 20], one obtains thenoise-averaged Fokker-Planck equation associated with themicroscopic dynamics described by Eqs. (1) and (2), as (cid:18) D t + sχ ∂ θ (cid:19) P = ∂ s (cid:20) ( η sχ + T [ P ]) P (cid:21) + (cid:15)∂ s P , (3)where D t = ∂ t + v e θ · ∇ is the material derivative, P ( r , θ, s, t ) is the probability density of particles at po-sition r , with velocity in direction θ and spin s at time t , and T [ P ] is the aligning torque T [ P ] = − γ (cid:90) θ (cid:48) (cid:90) s (cid:48) F ( θ (cid:48) − θ ) P ( r , θ (cid:48) , s (cid:48) , t ) . (4)For simplicity we have assumed ˜ F ( θ, r ) = δ ( r ) F ( θ ) ,with F ( θ ) = sin( θ ) , neglecting interaction between pairsat different positions.We describe the large-scale dynamics in terms of a fewcoarse-grained fields that vary slowly relative to micro-scopic time scales. For polarized flocks in addition to thenumber density, ρ ( r , t ) , of active units and their polariza-tion current density, w ( r , t ) , we include the spin angularmomentum density, S ( r , t ) . These are obtained from theprobability density P as ρ ( r , t ) w ( r , t ) S ( r , t ) = (cid:90) θ (cid:90) s ˆ e θ s P ( r , θ, s, t ) . (5)To obtain a closed set of hydrodynamic equations for ρ , w and S = S ˆ z , we combine moment techniquesused to approximate the velocity-dependent part of theFokker-Planck equation [21] with the closure developedin Ref. 22, 23 to handle kinetic equations of active sys-tems (see Supplementary Material). To minimize the num-ber of parameters, we nondimensionalize the equations byscaling time with τ (cid:15) = η /(cid:15) , length with v τ (cid:15) and den-sity with ρ = N/L . The resulting equations are con-trolled by only two dimensionless parameters ˜ χ = τ η /τ (cid:15) and ˜ γ = τ (cid:15) /τ γ [29]. For simplicity, we drop the tildes andall parameters are dimensionless in the following discus-sion unless otherwise noted. The continuum equations aregiven by ∂ρ∂t = − ∇ · w , (6) D wt w = − (cid:2) α ( ρ ) + β | w | (cid:3) w − ∇ ρ + λ w ( ∇ · w )+Ω S × w + Ω ∇ × S + D w ∇ w , (7) D st S = − ∇ × (cid:2)(cid:0) α ( ρ ) + β | w | (cid:1) w (cid:3) + Ω w × ∇ w − λ s ( ∇ · w ) S − ξ S + D s ∇ S , (8)where D wt = ∂ t + λ w · ∇ and D st = ∂ t + λ s w · ∇ are convective derivatives, α ( ρ ) = (1 − γρ ) / (1 + χ ) , β = γ / [8(1 + χ )] and ξ = 1 /χ . Explicit expressions forall other dimensionless parameters are given in the supple-mentary material. A pressure-type term λ ∇| w | has beenneglected in Eq.(7) because this term is known to lead toa spurious instability even in the overdamped limit when λ is evaluated with the closure used here [23, 24]. Thisinstability has not been observed in particle simulations ofVicsek models. We have also verified that it is not obtainedin particle simulations of the inertial spin model.Equations (6-8) augment the flocking model of Tonerand Tu [14] by incorporating the dynamics of the spin cur-rent. When S is neglected, these equations reduce to theToner-Tu equations as derived by Farrell et al. [18] (butin the case of constant self-propulsion speed). As in theToner-Tu model, the vector field w plays the dual role ofpolarization density and flow velocity. In equilibrium sys-tems of rotors both the equations for the spin and the veloc-ity field v w would contain dissipative couplings describ-ing friction with the substrate proportional to the combi-nation S /χ − v ∇ × w , guaranteeing that the angularvelocity S /χ and the vorticity v ∇ × w be equal whenthe whole system is rotating as a rigid body [25, 26]. In thenonequilibrium system considered here, in contrast, fric-tional terms proportional to angular velocity and vorticitywill in general appear with different coefficients. The firstterm on the right hand side of Eq. (8) was not includedin previous phenomenological model [17] and has a natu-ral interpretation of a nonlinear, velocity-dependent vorti-cal friction. The “self-spinning” term S × w couples thecenter-of-mass motion to the turning dynamics. In contrastto systems of passive rotors [25, 27], in the self-propelledparticle model considered here, these two degrees of free-dom are coupled because the spinning angle also controlsthe direction of translational motion [1]. We expect theseequations will provide useful to describe a number of activesystems where collective turning controls the large-scaledynamics.The homogeneous steady states of the continuum equa-tions have uniform density, ρ = 1 , and zero mean valueof the spin, S = 0 . As in the Toner-Tu model withno angular inertia, there are two such states: an isotropicgas state, with w = 0 , and a polarized liquid or flock-ing state, with w = w ˆ x and w = (cid:112) − α /β , where α = α ( ρ = 1) . We have chosen the ˆ x axis along the di-rection of spontaneously broken symmetry. The isotropicstate is always linearly stable for γ < . We examine be-low the linear stability of the polarized state by consider-ing the dynamics of fluctuations. We let w = w ˆ x + δ w , ρ = 1 + δρ , S = ˆ z δs and introduce Fourier amplitudes ( δρ, δ w , δs ) = (cid:80) q ( ρ q , w q , s q ) e i q · r + σt to obtain a set oflinearized equations in Fourier space (see SupplementaryMaterial).For spatial variations along the direction of broken sym-metry ( q = q ˆ x ), w y q and s q decouple from ρ q and w x q . Thecoupled linear dynamics of fluctuations in the density andthe magnitude of polarization ( w x q ) is unaffected by angu-lar inertia and is controlled by a longitudinal propagatingmode, with propagation speed c ρ = | α ρ | / (2 βw ) , where α ρ = ∂ ρ α . This mode goes unstable when γ < / ,corresponding to region A in Fig.2a. This instability isknown in Vicsek and Toner-Tu models as banding instabil-ity, but has recently been identified as the spinodal bound-ary within the liquid-gas coexistence region (Fig.1 left col-umn) [12, 22, 24, 28]. The coupled dynamics of spinand bending fluctuations ( w y q ) gives rise to overdamped,finite-wavelength spin waves that mediate the propaga-tion of turning information throughout the flock with wavespeed c s = w √ Ω Ω that increases with alignmentstrength. The existence of such propagating spin waves has been demonstrated on the basis of general arguments [1]and phenomenological continuum models [17], where theywere dubbed “second sound”.For wavevectors along any directions other than the di-rection of broken symmetry, all four equations are cou-pled and the analysis of the modes is rather cumbersome.For small wavevectors, we find two stable and relaxationalmodes that will not be discussed further and two hydrody-namic propagating modes, with dispersion relation σ ± ( q, θ ) = ic ± ( θ ) q − D sw ( θ ) q + O ( q ) , (9)and wave velocity c ± ( θ ) = α ρ cos( θ ) ± (cid:113) α ρ cos ( θ ) + 8 β w sin ( θ )4 βw , (10)where θ is the angle between the direction of q and thedirection of broken symmetry. The full expression for thedamping D sw ( θ ) is not instructive thus is not given here.For θ = 0 , c − (0) = α ρ / (2 βw ) and the mode σ − yieldsthe banding instability that delimits the spinodal region ofmicrophase separation [12]. For arbitrary angle θ , however,both modes are propagating with anisotropic speed and de-scribe coupled fluctuations of density, spin, and bend/splaydeformations of the polarization field. The angular depen-dence of the instability is shown in Fig.2b that displays theregions where D sw < . At small angles the instability isdriven by density fluctuations, as in the Toner-Tu model.At large angles the instability is dominated by spin fluc-tuations. At θ = π/ the longitudinal banding instabilityis suppressed and the dynamics is controlled by transversespin wave propagating at speed | c ± ( π/ | = 1 / √ . Interms of our dimensionless parameters, this transverse spinwave is unstable for γ > (1 + 4 χ )(1 + χ ) / (8 χ ) + 4 , cor-responding to region B in Fig.2a. The instability is drivenby the growth of bend − ∇ × [( α ( ρ ) + β | w | ) w ] andsplay λ w ( ∇ · w ) deformations augmented by the spinwave through the self-rotation term Ω S × w . This long-wavelength instability of the ordered state is a new resultof our work and will be referred to as spin-wave instabil-ity. It leads to a complex spatio-temporal dynamics withlarge density and spin fluctuations characterized by contin-uously turning and swirling flocks as confirmed by numer-ical solutions of the hydrodynamic equations and particlesimulations (see Fig.1 right column).By carrying out the small wavevector expansion of thedispersion relation Eq (9) up to fourth order in q we canidentify the wavector q c of the fastest growing mode corre-sponding to the maximum of Re [ σ ± t ( q )] shown in Fig. 2cand d for various values of γ and χ . This defines the char-acteristic length scale λ c ∼ /q c that can be thought of ascontrolling the size of the turning flock at the linear level.To gain more insight on the complex spatio-temporalstructures that emerge in the unstable regions of parame-ters and to confirm the results of the linear stability anal-ysis, we have solved numerically Eqs. (6-8) with periodic (cid:97) (cid:114) Γ Θ ! ! q Σ t R ! q " ! ! q Σ t R ! q " c (cid:1) d (cid:1) Increase)) (cid:1) Decrease)) (cid:1) a (cid:1) A (cid:1) B (cid:1) b (cid:1) A (cid:1) B (cid:1) FIG. 2: a. Phase diagram in the plane of dimensionless γ and χ . b. Phase diagram in the plane of dimensionless γ and θ at χ = 1 . . The shaded region A is the spinodal portion of the re- gion of coexistence of disordered gas (existing for γ < ) andtraveling bands of polar liquid (existing for γ > / ). The co-existence region is delimited by the binodals (not shown) and ex-tends inside the white regions, both to the right and to the left ofregion A, as verified via particle simulations. The shaded regionB corresponds to the region where the homogeneous polar liquidis linearly unstable to spin-waves as shown in Fig.1. c. Real partof the dispersion relation of the transverse mode σ ± t = σ ± ( π/ at χ = 1 and γ = 7 . , . , . , . . d. Real part of the disper-sion relation of the transverse mode σ ± t = σ ± ( π/ at γ = 9 and χ = 0 . , . , . , . . (cid:97) c s (cid:97) c s (cid:97) c s a b FIG. 3: a: Snapshot of the anisotropic spin wave in the polarizedstate at γ = 7 . and χ = 2 . . Color indicates the spin currentdensity. b: Speed of spin waves in the polarized state as a functionof alignment strength γ for χ = 1 . , . , . (red, blue, black) inthe directions longitudinal (circles) and transverse (squares) tothat of mean polarization. The dashed line is the transverse speed | c ± ( π/ | in Eqn.10. The system is evolved for 5000 time steps. boundary conditions starting from the homogeneous polarstate with small perturbations. The results are summarizedin the phase diagram of Fig.2a. The shaded region A isbounded to the left by the line γ = 2 where the disor-dered gas is linearly unstable and to the right by the line γ = 8 / where the homogeneous polar liquid is linearlyunstable to longitudinal fluctuations (the banding instabil-ity). These instability lines delimit the spinodal portionof the gas/liquid coexistence region and are distinct from the binodal lines that mark the boundaries of such a re-gion [12]. In fact particle simulations reveal that the co-existence region extends to the left and right of region A.The squares in Fig.2a correspond to mean density fluctu-ations ∆ ρ = (cid:113) N (cid:80) r < ( ρ ( r ) − ρ ) >/ρ = 0 . ,with N the number of grid points, evaluated in the contin-uum model, starting in a uniform polar state. The shadedregion B is the region where the homogenous polar liquidis linearly unstable to spin wave fluctuations. Again, par-ticle simulations show that the inhomogeneous spinningbands are found beyond the linear stability boundary thatdelimits region B, suggesting that this region is also a spin-odal region. The diamonds correspond to spin fluctuations ∆ S = (cid:113) N (cid:80) r < ( S ( r ) − < S > ) > = 0 . . In theoverdamped limit χ → , the spin-wave instability van-ishes due to the rapid decay of spin current fluctuationsover time χ/η , and the dynamics of the system is con-trolled solely by a rescaled alignment strength γ , with ageneric banding instability close to the flocking transition,as in the Vicsek and the Toner-Tu models [22, 24, 28]. Ourresult, together with Ref. [17], highlights for the first timethe importance of inertia in controlling dynamics of activepolar systems at large length scales.To understand the nature of the spin waves that medi-ate the transfer of turning information within the flock, westudy the propagation of the spin waves numerically withEqs. (6-8) by initializing the system in the uniformly po-larized state, with a concentrated spin current at the cen-ter (Fig. 3a). We measure the longitudinal and transversespeed as a function of alignment strength γ for various χ and plot the results in Fig. 3b. The longitudinal speed (cir-cles) increases with the strength of alignment interactionwhile the transverse speed (squares) stays approximatelyconstant over the range of parameters.In the longitudinal direction, where δw y and δs de-couple from δw x and δρ , the spin wave is governed bya damped wave equation at finite wavelength with wavespeed c s = w √ Ω Ω proportional to alignment strength.In the transverse direction, all fluctuations are coupled andthe dynamics is governed at long wavelength by the hy-drodynamic mode (Eq.9) with an angular-dependent prop-agating speed that reduces to | c ± ( π/ | = 1 / √ in thetransverse direction as given in Eqn.10, and fits the dataquantitatively in Fig.3b.We have derived continuum equations that generalize theToner-Tu model of flocking to incorporate turning inertiaby coarse-graining the active inertial spin model proposedrecently by Cavagna et al. [1]. The coarse-graining sim-plifies the analysis by shrinking the number of independentparameters to two. The interplay between rotational inertiaand bending elasticity of a polarized flock provides a mech-anism for the propagation of turning information throughthe flock in the form of collective spin-wave excitations.By studying the continuum equations analytically and nu-merically, we predict a new instability of the polarized stateassociated with large density and spin current fluctuationsthat leads to complex spatio-temporal patterns of continu-ously swirling and rotating flocks. This long-wavelengthinstability is associated with the growth of anisotropic spinwaves and is referred to as spin-wave instability .We thank Sriram Ramaswamy and Andrea Cavagna foruseful discussions. The research leading to this workwas supported by the National Science Foundation (NSF)awards DMR-1305184 and DGE-1068780 at SyracuseUniversity and NSF award PHY11-25915 and the Gordonand Betty Moore Foundation Grant No. 2919 at the KITPat the University of California, Santa Barbara. MCM alsoacknowledges support from the Simons Foundation. APPENDIX A: DIMENSIONLESS PARAMETERS
It is useful to clarify our choice of dimensionless param-eters by making contact with the non-inertial Vicsek modelfamiliar from the literature. The continuous-time Vicsekmodel can be obtained from Eqs. (1) and (2) of the maintext by letting χ = 0 and eliminating s i , with the result, d r i dt = v ˆ e θ i ; , (11) dθ i dt = γη (cid:88) j ˜ F ( θ j − θ i , r ji ) + (cid:115) (cid:15)η ξ i ( t ) , (12)where ˜ F ( θ ) = F ( θ ) / ( πR ) for | r ij | ≤ R and zero oth-erwise, with F ( θ ) = sin θ . The non-inertial limit cor-responds to a Vicsek model with alignment strength γ/η and noise amplitude (cid:15)/η . The additional parameters inEqs. (11) and (12) are v , the mean density ρ = N/L ,and the radius R of the interaction. If we scale lengthswith R and times with R/v the continuous time Vicsekmodel described by Eqs. (11) and (12) contains three di-mensionless parameters: the scaled noise, (cid:15)R/ ( η v ) , thescaled alignment strength, γ/ ( ηRv ) , and the mean den-sity ρ R . The discrete time Vicsek model can be re-covered by assuming that the alignment is instantaneous, i.e., ( γ/ηR ) − is short compared to all other times scales(specifically τ (cid:15) = η /(cid:15) and the time step for updatingthe dynamics). The resulting model contain two dimen-sionless parameters: the mean density ρ R and the noise (cid:15)R/ ( η v ) , as expected.Alternatively, in Eqs. (11) and (12) we can scale timeswith τ (cid:15) and lengths with v τ (cid:15) . The microdynamics thentakes the form d r i dt = ˆ e θ i ; , (13) dθ i dt = γη(cid:15) (cid:88) j ∈ C ˜ R F ( θ j − θ i ) + √ ξ i ( t ) , (14)where R i and time are now all dimensionless, F ( θ ) =sin θ and we have made explicit the dependence on the in-teraction range, with C ˜ R a circle of radius ˜ R = Rτ (cid:15) /v .The mean field limit of these equations will only depend onthe dimensionless parameter γηρ /(cid:15) = τ (cid:15) /τ γ . There are,however, two additional parameters that provide cutoffs tothe mean-field theory: the interaction range at small scalesand the system size at large scales, both scale with v τ (cid:15) .In other words, although seemingly magically rewritten interms of a single parameter, the model still contains threeindependent parameters. When comparing to solution ofthe nonlinear PDE’s obtained in mean-field to the resultsof particle simulations where we set R = 1 one shouldthink of the density ρ R and γη/ ( (cid:15)R ) as independentparameters.For the inertial continuous time model described by Eqs.(1) and (2) of the main text, the same transformation yieldsa mean field theory that contains only two dimensionlessparameters, defined as ˜ γ and ˜ χ in the main text. To these,however, we must add the two cutoffs at large and smallscales. The non-inertial limit is recovered for ˜ χ = 0 . APPENDIX B: DERIVATION OF THE HYDRODYNAMICEQUATIONS
The Fokker-Planck equation for the one-particle proba-bility density P ( r , θ, s, t ) associated with Eqs. (1) and (2)of the main text is given by ˙ P ( r , θ, s, t ) + v θ · ∇ P = − ∂∂θ ( 1 χ sP ) + T ( θ, r , t ) ∂P∂s + ∂∂s ( ηχ sP ) + (cid:15) ∂ P∂s , (15)where T ( θ, r , t ) = − γ (cid:82) π − π dθ (cid:48) F ( θ (cid:48) − θ ) p ( r , θ (cid:48) , t ) is the torque. We have assumed local interaction F ( θ, r ) = δ ( r ) sin ( θ ) and defined p ( r , θ (cid:48) , t ) = (cid:82) s P ( r , θ (cid:48) , s, t ) [30]. To make the notation more compact, we define the Fokker-Planck operator as L k = L rev + L ir , (16) L rev = − sχ ∂∂θ + T ( r , θ, t ) ∂∂s − v θ · ∇ , (17) L ir = ηχ ∂∂s ( s + s ∂∂s ) , (18)where L rev and L ir represent the reversible and irriversiblepart of the Fokker-Planck operator respectively, and wehave introduced the steady state value of the spin s = (cid:15)χ/η . In the absence of interaction and activity, the steadystate distribution of the spin, obtained by setting the timederivative to zero, has a Maxwell-like form, given by P ( s ) = 1 (cid:112) πs exp( − s s ) . (19)Following standard methods [21], we transform theFokker-Planck operator by multiplying it from the rightand the left by φ ( s ) = P ( s ) and φ − ( s ) = P − ( s ) ,respectively, with the result ¯ L k = φ − ( s ) L k φ ( s ) = ¯ L rev + ¯ L ir , (20) ¯ L ir = − ηχ b + b , ¯ L rev = − bD − b + ˆ D − v θ · ∇ , (21)where b + and b are creation and annihilation operators, re-spectively. b + φ n ( s ) = √ n + 1 φ n +1 ( s ) , (22) bφ n ( s ) = √ nφ n − ( s ) . (23) D and ˆ D are the differential operators, with the latter con-taining the information of the interaction, b + = − s ∂∂s + 12 ss , b = s ∂∂s + 12 ss , (24) D = s χ ∂∂θ , ˆ D = s χ ∂∂θ + T ( θ, r , t ) s . (25)The normalized eigenfunctions φ n ( s ) of the operator ¯ L ir = − ηχ b + b are defined by the eigenvalue equation ¯ L ir φ n ( s ) = − ηχ nφ n ( s ) , (26) with φ n ( s ) = ( b + ) n φ ( s ) / √ n ! , (27) φ ( s ) = exp( − s s ) / (cid:113) s √ π. (28)Finally, φ n ( s ) are related to the physicists’ Hermite poly-nomials H n ( x ) = (cid:0) x − ddx (cid:1) n · as φ n ( s ) = H n ( s √ s ) exp( − s s ) / (cid:113) n !2 n s √ π. (29)We now expand the probability distribution function interms of φ n ( s ) , P ( r , θ, s, t ) = φ ( s ) ∞ (cid:88) n =0 p n ( r , θ, t ) φ n ( s ) , (30)and we insert the expansion into the Fokker-Planck equa-tion, ∂ t P ( r , θ, s, t ) = L k P ( r , θ, s, t ) , (31)where the Fokker-Planck operator is obtained after an in-verse transformation, as L k = φ ( s )( − ηχ b + b − bD − b + ˆ D − v · ∇ ) φ − ( s ) . (32)Using the properties of the operators and the orthogonalityof the Hermite polynomials, we obtain a hierachy of equa-tions for the moments p n ( r , θ, t ) , D t p n ( r , θ, t ) = − ηχ np n ( r , θ, t ) − √ n + 1 Dp n +1 ( r , θ, t ) − √ n ˆ Dp n − ( r , θ, t ) , (33)where D t = ∂ t + v θ · ∇ is the material derivative. Ex-plicitly, the equations for the first three moments are givenby D t p = − Dp , (34) D t p = − ηχ p − √ Dp − ˆ Dp , (35) D t p = − ηχ p − √ Dp − √ Dp . (36)The first two moments are related to the probability density c ( r , θ, t ) of finding a particle at r , with velocity directed along θ at time t and the spin current j ( r , θ, t ) as c ( r , θ, t ) = p = (cid:90) ∞−∞ P ( r , θ, s, t ) ds , (37) j ( r , θ, t ) = s p = (cid:90) ∞−∞ sP ( r , θ, s, t ) ds. (38)To obtain closed equations for c and j , we set D t p = 0 fortimes long compared to χ/ η , and let p n = 0 for n ≥ .We then eliminate p in favor of p and p to obtain closedequations. The equations for density and current are thengiven by D t c ( r , θ, t ) = − χ ∂j∂θ , (39) D t j ( r , θ, t ) = − ηχ j + (cid:15)η ∂ j∂θ + 1 η ∂ [ T ( r , θ, t ) j ] ∂θ (40) − (cid:15)η ∂c∂θ − T ( r , θ, t ) c. The goal is to obtain closed equations for the number den-sity ρ ( r , t ) , polarization density w ( r , t ) and spin current S ( r , t ) , which are the conserved, symmetry-breaking andrelevant dynamic variables in the flocking system, respec-tively. Generalizing the method described in Ref.[22], we introduce the angular Fourier transform of c and j as c k ( r , t ) = (cid:90) π − π c ( r , θ, t ) e ikθ dθ , (41) j k ( r , t ) = (cid:90) π − π j ( r , θ, t ) e ikθ dθ , (42)which are related to ρ ( r , t ) , w ( r , t ) and S ( r , t ) by ρ ( r , t ) = c ( r , t ) , S ( r , t ) = j ( r , t ) , (43) w x ( r , t ) = Re [ c ( r , t )] , w y ( r , t ) = Im [ c ( r , t )] , (44)whose dynamic equations are ∂ t c k ( r , t ) + v ∇ ∗ c k +1 + v ∇ c k − = ikχ j k , (45) ∂ t j k ( r , t ) + v ∇ ∗ j k +1 + v ∇ j k − = − η k χ j k + ik(cid:15)η c k + ikγ πη (cid:88) m j k − m F − m c m + γ π (cid:88) m c k − m F − m c m , (46)where ∇ = ∂ x + i∂ y , ∇ ∗ = ∂ x − i∂ y and F ± = ± iπ . Wehave introduced an effective friction η k = η + k (cid:15)χ/η .Explicity, the equations for c , c and j are given by ∂ t c + v ∇ ∗ c + v ∇ c ∗ = 0 , (47) ∂ t c + v ∇ ∗ c + v ∇ c = iχ j , (48) ∂ t j + v ∇ ∗ j + v ∇ j ∗ = − ηχ j . (49) To close these equations, we need to express j and c interms of c , c and j . To do so, we consider the equationsfor j , j and c , ∂ t j + v ∇ ∗ j + v ∇ j = − η χ j + i(cid:15)η c − γ η ( j c ∗ − j c ) + iγ c c ∗ − c c ) , (50) ∂ t j + v ∇ ∗ j + v ∇ j = − η χ j + 2 i(cid:15)η c − γη ( j c ∗ − j c ) + iγ c c ∗ − c c ) , (51) ∂ t c + v ∇ ∗ c + v ∇ c = 2 iχ j . (52)For times long compared to χ/η , we set ∂ t j = ∂ t j =0 . Retaining terms up to first order in χ/η we obtain the expression for j and j , j = χη (cid:20) i(cid:15)η c + γ η j c + iγ c c ∗ − c c ) − v ∇ j (cid:21) + O ( χ ) , (53) j = χη ( 2 i(cid:15)η c − iγ c ) + O ( χ ) . (54)Inserting Eq. 54 into the equation for c , we obtain, ∂ t c + v ∇ ∗ c + v ∇ c = γη c − (cid:15)ηη c . (55)For times long compared to ηη / (4 (cid:15) ) , we follow themethod of Ref. [22, 23] and set ∂ t c = 0 and c n = 0 for n ≥ to obtain the expression for c , c = γη (cid:15) c − v ηη (cid:15) ∇ c . (56)Using the expressions for j and c , we obtain the closedequations, ∂c ∂t + v ∇ c ∗ + v ∇ c = 0 , (57) ∂c ∂t + v ∇ c + v γη (cid:15) ∇ ∗ c = (cid:18) γ η c − (cid:15)ηη − γ η (cid:15)η | c | (cid:19) c + iγ ηη j c − iv η ∇ j + γv ηη (cid:15)η c ∗ ∇ c + v ηη (cid:15) ∇ c , (58) ∂j ∂t + iχγv ηη (cid:15)η ( ∇ [( ∇ ∗ c ∗ ) c ]) − ∇ ∗ [( ∇ c ) c ∗ ]) = − v γχ η [ i ∇ ( c c ∗ ) − i ∇ ∗ ( c c )] − v (cid:15)χ ηη ( − i ∇ c ∗ + i ∇ ∗ c ) (59) − v γ ηχ (cid:15)η (cid:2) − i ∇ ( | c | c ∗ ) + i ∇ ∗ ( | c | c ) (cid:3) − v γχ ηη [ ∇ ( c ∗ j ) + ∇ ∗ ( c j )] + χv η ∇ j − ηχ j . (60)Using the following identities, ∇ ∗ c = (cid:2) w · ∇ ) w + 2 w ( ∇ · w ) − ∇| w | (cid:3) , (61) ic j = S × w , (62) i ∇ j = −∇ × S (63) i ∇ ( c c ∗ ) − i ∇ ∗ ( c c ) = − ∇ × ( ρ w ) , (64) − i ∇ c ∗ + i ∇ ∗ c = 2 ∇ × w , (65) − i ∇ ( | c | c ∗ ) + i ∇ ∗ ( | c | c ) = 2 ∇ × ( | w | w ) , (66) ∇ ( c ∗ j ) + ∇ ∗ ( c j ) = 2 S ∇ · w + 2( w · ∇ ) S , (67) c ∗ ∇ c = ( w · ∇ ) w − w ( ∇ · w ) + 12 ∇| w | , (68) i ( ∇ [( ∇ ∗ c ∗ ) c ]) − ∇ ∗ [( ∇ c ) c ∗ ]) = − w × ∇ w , (69) we finally obtain the hydrodynamic equations [31], ∂ρ∂t = −∇ · ( v w ) , (70) ∂ w ∂t + v ∇ ρ + v γη (cid:15) (cid:2) w · ∇ ) w + 2 w ( ∇ · w ) − ∇| w | (cid:3) = (cid:18) γ η ρ − (cid:15)ηη − γ η (cid:15)η | w | (cid:19) w (71) + γ ηη S × w + v η ∇ × S + γv ηη (cid:15)η (cid:20) ( w · ∇ ) w − w ( ∇ · w ) + 12 ∇| w | (cid:21) + v ηη (cid:15) ∇ w , (72) ∂ S ∂t = v γχ η ∇ × ( ρ w ) − v χ(cid:15)ηη ∇ × w − v γ χη (cid:15)η ∇ × ( | w | w ) − v γχ ηη [ S ( ∇ · w ) + ( w · ∇ ) S ]+ χγv ηη (cid:15)η w × ∇ w + χv η ∇ S − ηχ S . (73) APPENDIX C: MODE ANALYSIS
We start with the dimensionless hydrodynamic equa-tions. Time is scaled by the rotational diffusion time τ (cid:15) = η /(cid:15) and length by the persistence length v τ (cid:15) . w and ρ are scaled by the average number density ρ and S by ρ χ/τ (cid:15) , leading to [32] ∂ρ∂t = −∇ · w , (74) ∂ w ∂t + λ ( w · ∇ ) w = − (cid:2) α ( ρ ) + β | w | (cid:3) w − ∇ ρ + λ w ( ∇ · w ) + Ω S × w + Ω ∇ × S + D w ∇ w , (75) ∂ S ∂t + λ s ( w · ∇ ) S = − ∇ × (cid:2) ( α ( ρ ) + β | w | ) w (cid:3) + Ω w × ∇ w − λ s S ( ∇ · w ) − ξ S + D s ∇ S . (76)All parameters are related to two microscopic dimension-less variables: the scaled alignment strength ˜ γ = τ (cid:15) /τ γ and inertia ˜ χ = τ η /τ (cid:15) , where τ (cid:15) = η /(cid:15) , τ η = χ/η and τ γ = η/ ( γρ ) are the three natural timescales in the sys- tem corresponding to rotational diffusion, frictional dissi-pation and alignment interaction.We drop the tilde in the following discussion for simplic-ity of notation. α ( ρ ) = 11 + χ (1 − γρ , β = 11 + χ γ , Ω = χγ χ ) , Ω = χ χ ) , Ω = γ
16 ( 1 + 4 χ χ ) ,λ = γ − γ
16 ( 1 + 4 χ χ ) , λ = − (cid:20) γ γ
16 ( 1 + 4 χ χ ) (cid:21) λ s = χγ χ ) ,ξ = 1 χ , D w = 1 + 4 χ , D s = χ χ ) . To perform linear mode analysis, we restrict ourselves tothe 2D planar case. The isotropic state is always linearlystable therefore trivial, and we focus on the uniformly po-larized state for γ > with the direction of spontaneous broken symmetry along ˆ x . Perturbing around the polarizedstate ρ = 1 + δρ , w = w ˆ x + δ w and S = δS ˆ z , we arriveat the linearized equations ∂δρ∂t = −∇ · δ w , (77) ∂δ w ∂t + λ w ∂ x δ w = ( µ δρ + µ δw x ) w ˆ x − ∇ δρ + λ w ˆ x ∇ · δ w + Ω δS ˆ z × w ˆ x + Ω ∇ × δS ˆ z + D w ∇ δ w , (78) ∂δS ˆ z∂t + λ s w ∂ x ( δS z ˆ z ) = ∇ × [( µ δρ + µ δw x ) w ˆ x ] + Ω w ˆ x × ∇ δw − ξδS ˆ z + D s ∇ ( δS ˆ z ) , (79)where µ = − ∂ ρ α = γ χ ) , µ = − βw = − w γ χ ) and w = (cid:104) ( γ − γ (cid:105) / . Longitudinal mode q = q x Considering mode along the direction of broken symme-try, we obtain0 σδρ = − iqδw x , (80) σδw x = µ w δρ + µ w δw x − iq δρ + iq ( λ − λ ) w δw x − D w q δw x , (81) σδw y = − iqλ w δw y + Ω w δS − iq Ω δS − D w q δw y , (82) σδS = − q Ω w δw y − ξδS − iqw λ s δS − D s q δS (83) “Banding Instability” Notice that δρ and δw x decou-ple from δw y and δS , leading to the dispersion relation σ l ( q ) = iµ µ q + 1 µ w (cid:20) ( λ − λ ) w µ µ + 12 − µ µ (cid:21) q + O ( q ) . (84)Fluctuations in density and magnitude of polarization leadto the “banding instability” close to the isotropic-polarphase transition as generally observed in polar active fluid,the condition of which is given by ( λ − λ ) w µ µ + 12 < µ µ . (85)In terms of the microscopic parameters, it reads γ < . (86) Spin wave
Dynamics of δw y and δS gives rise to thespin wave, carrying the information of turning. Neglecting convections and diffusion, the dispersion relation for thespin wave is σ ± s = − ξ ± c s q (cid:115) [ ξ/ (2 c s )] q − , (87)where c s = w (cid:112) Ω Ω = (cid:20) ( γ − χ (1 + 4 χ )4(1 + χ ) (cid:21) (88)is the wave speed. Transverse mode q = q y Transverse instability
Transverse mode is governedby the full coupled equations: − iq µ w µ w − D w q iqλ w iq Ω − iq iqλ s w − D w q Ω w − iqw µ − iqw µ − q Ω w − ξ − D s q δρδw x δw y δS = σ t ( q ) δρδw x δw y δS , which leads to the dispersion relation once treated pertur-batively in the long wavelength limit: σ ± t ( q ) = ± √ i q + 1 ξ (cid:18) w µ Ω µ − Ω Ω w − λ w Ω − D w ξ (cid:19) q + O ( q ) , (89)from which the condition for transverse instability is ob-tained as w Ω µ µ − λ w Ω > D w ξ + Ω Ω w , (90) or in terms of microscopic parameters γ > (1 + 4 χ )(1 + χ )8 χ + 4 . (91)The phase diagram is plotted in Fig.2a in the main text,with quantitative agreement between the numerical and an-alytical phase boundaries. This transverse instability ren-1ders the system spatially inhomogeneous with large densityand spin fluctuations characterized by continuously turningand swirling flocks with propagating spin waves. There-fore, we term it the spin-wave instability. The spatial-temporal patterns have been observed from both the numer-ical simulations of the hydrodynamic equations and parti-cle simulations [33].To understand the origin of the instability, we write downthe minimal equations that yield this instability. For clarity,we write down the dimensionful form. ∂ t w = − β | w | w + λ w ( ∇ · w ) + Ω S × w , (92) ∂ t S = − v χ ∇ × (cid:0) β | w | w (cid:1) − ξ S , (93)where λ = − v γη/ (4 (cid:15) ) − γv ηη / (16 (cid:15)η ) < , β = ηγ / (8 (cid:15)η ) , Ω = γ/ (2 ηη ) and ξ = η/χ . The lin-earized equations are ∂ t δw x = µ w δw x + λ w ∂ y δw y , (94) ∂ t δw y = Ω w δs z , (95) ∂ t δs z = − µ s w ∂ y δw x − ξδs z , (96)where µ = − βw < and µ s = v χµ < . Theylead to the dispersion relation σ t ( q ) = − w λ µ s Ω ξµ q + O ( q ) , (97)which yields the instability condition w λ µ s Ω ξµ < . (98)This condition can be interpreted as the growth of bendand splay deformations augmented by the spin wave. If weinclude the density-dependent alignment interaction, rota-tional diffusion and spin elasticity, all of which serve asstabalization factors, we recover the full condition 90. Thecompetition among these effects yields the spin-wave in-stability, which is model-dependent. [1] A. Cavagna, L. D. Castello, I. Giardina, T. Grigera, A. Jelic,S. Melillo, T. Mora, L. Parisi, E. Silvestri, M. Viale, et al.,J. Stat. Phys. , 601 (2014).[2] C. W. Reynolds, Computer Graphics , 25 (1987).[3] T. Vicsek, A. Czir´ok, E. Ben-Jacob, I. Cohen, andO. Shochet, Phys. Rev. Lett. , 1226 (1995).[4] P. Romanczuk, M. B¨ar, W. Ebeling, B. Lindner, andL. Schimansky-Geier, Eur. Phys. J Special Topics , 1(2012).[5] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cis- bani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Pro-caccini, et al., PNAS , 1232 (2008).[6] B. Szab´o, G. J. Szolosi, B. Gonci, Z. Juranyi, D. Selmeczi,and T. Vicsek, Phys. Rev. E , 061908 (2006).[7] V. Schaller, C. Weber, C. Semmrich, E. Frey, and A. R.Bausch, Nature , 73 (2010). [8] Y. Sumino, K. H. Nagai, Y. Shitaka, D. Tanaka,K. Yoshikawa, H. Chat´e, and K. Oiwa, Nature , 448(2012).[9] A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot, andD. Bartolo, Nature , 95 (2013).[10] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liver-pool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. , 1143 (2013).[11] S. Ramaswamy, Annual Review of Condensed MatterPhysics , 323 (2010).[12] A. P. Solon, H. Chat´e, and J. Tailleur, Phys. Rev. Lett. ,068101 (2015).[13] G. Gr´egoire and H. Chat´e, Phys. Rev. Lett. , 025702(2004).[14] J. Toner and Y. Tu, Phys. Rev. Lett. , 4326 (1995).[15] J. Toner, Y. Tu, and S. Ramaswamy, Ann. Phys. , 170(2005).[16] A. Attanasi, A. Cavagna, L. D. Castello, I. Giardina, T. S.Grigera, A. Jeli´c, S. Melillo, L. Parisi, O. Pohl, E. Shen,et al., Nature Physics , 691 (2014). [17] A. Cavagna, I. Giardina, T. S. Grigera, A. Jelic, D. Levine,S. Ramaswamy, and M. Viale, Phys. Rev. Lett. , 218101(2015).[18] F. Farrell, M. C. Marchetti, D. Marenduzzo, and J. Tailleur,Phys. Rev. Lett. , 248101 (2012).[19] D. S. Dean, J. Phys. A: Math. Theor. , L613 (1996).[20] R. Zwanzig, Nonequilibrium Statistical Mechanics (OxfordUniversity Press, 2001).[21] H. Risken,
The Fokker-Planck Equation (Springer-Verlag,1988), 2nd ed.[22] E. Bertin, M. Droz, and G. Gr´egoire, J. Phys. A: Math.Theor. , 445001 (2009).[23] A. Peshkov, E. Bertin, F. Ginelli, and H. Chat´e, Eur. Phys. JSpecial Topics , 1315 (2014).[24] S. Mishra, A. Baskaran, and M. C. Marchetti, Phys. Rev. E , 061916 (2010).[25] J. C. Tsai, F. Ye, J. Rodriguez, J. P. Gollub, and T. C. Luben-sky, Phys. Rev. Lett. , 214301 (2005).[26] E. Braun, O. L. Fuchs, and S. Godoy, Chemical Physics Let-ters , 434 (1997).[27] M. W. James T. Hynes, Raymond Kapral, Physica A , 427(1977).[28] E. Bertin, M. Droz, and G. Gr´egoire, Physical Review E ,022101 (2006).[29] See the SI for a description of our choice of parameters.[30] This interaction does not describe the mean polarizationdeep in the ordered state. We have verified that better be-haved models such as F ( θ ) ∼ sin( θ/ , give equations ofthe same structure and do not affect the qualitative behaviorand instabilities.[31] An alternative closure proposed in [9] yields continuumequations with the same structure as those obtained here,but with different coefficients.[32] If we neglect ∂ t S in Eq. (76) and use the resulting equationsto eliminate S in favor of ρ and w , the resulting continuumequations have the same structure as those obtained in [18], with O ( τ η /τ (cid:15) ) corrections to various coefficients.[33] We have performed extensive particle simulations that con-firm the existence of a region of turning flocks and largespin density fluctuations at large γ . Typical snapshots fromsimulations are shown in Fig. 1 in the main text, but the full2