PRL 113, 258104 (2014)
PHYSICAL REVIEW LETTERS
week ending 19 DECEMBER 2014
Vortex Arrays and Mesoscale Turbulence of Self-Propelled Particles Robert Großmann,1,* Pawel Romanczuk,1 Markus Bär,1 and Lutz Schimansky-Geier2 1
Physikalisch-Technische Bundesanstalt Berlin, Abbestraße 2-12, 10587 Berlin, Germany Department of Physics, Humboldt-Universität zu Berlin, Newtonstraße 15, 12489 Berlin, Germany (Received 26 May 2014; revised manuscript received 11 September 2014; published 19 December 2014) 2
Inspired by the Turing mechanism for pattern formation, we propose a simple self-propelled particle model with short-range alignment and antialignment at larger distances. It is able to produce orientationally ordered states, periodic vortex patterns, and mesoscale turbulence, which resembles observations in dense suspensions of swimming bacteria. The model allows a systematic derivation and analysis of a kinetic theory as well as hydrodynamic equations for density and momentum fields. A phase diagram with regions of pattern formation as well as orientational order is obtained from a linear stability analysis of these continuum equations. Microscopic Langevin simulations of self-propelled particles are in agreement with these findings. DOI: 10.1103/PhysRevLett.113.258104
PACS numbers: 47.63.Gd, 05.10.Gg, 87.10.Mn, 87.18.Gh
The term “active matter” refers to nonequilibrium systems of interacting, self-propelled entities which are able to take up energy from their environment and convert it into motion [1,2]. Examples, such as cytoskeletal filaments [3], chemically driven colloids [4], or flocks of birds [5] have recently received a lot of attention in physics, chemistry, and biology. They exhibit a wide range of collective phenomena which are absent in systems at thermodynamic equilibrium, for example large-scale traveling bands and polar clusters [6,7] as well as arrays of vortices [8,9]. In this context, bacteria represent important model systems which have been used to investigate such different aspects as clustering [10] and rheological properties [11] of active matter systems. Recently, irregular vortex structures were experimentally observed in dense suspensions of swimming bacteria [11–17]. In addition, a phenomenological model with a Swift-Hohenberg coupling was proposed which describes the observed behavior including the power spectrum of the bacterial dynamics [16,18]. The spectra reveal pattern formation on finite length scales. The dynamic vortex patterns observed in experiments and simulations are reminiscent of a turbulent state which led the authors to denominate this new phenomenon mesoscale turbulence. In this Letter, we propose a self-propelled particle model which is capable of producing such mesoscopic spatiotemporal patterns based on competing interactions: local alignment interaction at short length scales and antialignment at larger distances. Our model is inspired by the Turing mechanism of short-range activation and long-range inhibition in reaction-diffusion systems sufficient to explain instabilities at finite length scales [19,20]. In this context, activation corresponds to creation of local polar order, whereas antialignment inhibits the emergence of long-range polar order. The effects of related competing interactions were also studied in the Ising model with 0031-9007=14=113(25)=258104(5)
ferromagnetic nearest neighbor and antiferromagnetic nextnearest neighbor coupling [21–25], as well as in the Kuramoto model [26]. We abstain from considering a detailed model of individual swimmers immersed and interacting through a surrounding fluid. However, we note that hydrodynamic interactions of bacteria have been shown to destabilize polar order at high densities [27], and antialignment is the simplest interaction inhibiting local polar order. Thus, our model has to be viewed as a generalization of the seminal Vicsek model [28], with its focus on generic behavior emerging from simple, effective interactions instead of detailed modeling of specific interactions derived from first principles. The original Vicsek model displays surprisingly complex spatiotemporal behavior [29,30], but it does not exhibit mesoscale turbulence or stationary vortex patterns. In contrast, our generalized model including antialignment produces both orientationally ordered states and periodic vortex arrays. Moreover, a region with spatiotemporal irregular behavior of the order parameter field is found between both states. The simulated patterns reproduce most features found in the turbulent state in experiments with dense suspensions of swimming bacteria and corresponding phenomenological macroscopic models [16,18]. The simplicity of our microscopic model makes it analytically tractable and allows us to derive kinetic equations as well as approximated equations for the coarse-grained density and momentum fields. Thus, we can establish a link between microscopic parameters and the hydrodynamic transport coefficients. The resulting hydrodynamic model is a modification of the Toner-Tu equations [31,32] similar to the phenomenological model of bacterial turbulence [16,33]. We show that antialignment interaction may lead to a negative viscosity of the order parameter field. The latter was proposed in [16,33] on a
258104-1
© 2014 American Physical Society
PRL 113, 258104 (2014)
PHYSICAL REVIEW LETTERS
phenomenological basis. Moreover, we analyze the stability of the spatially homogeneous, isotropic state in both frameworks and compare our results to direct simulations of the microscopic dynamics. Model.—We consider N particles moving with constant speed v0 in a two-dimensional system with periodic boundary conditions. Each particle interacts with all neighbors within a finite interaction range lint . We work in natural units such that the particle mass, speed, and interaction range equal one (m ¼ 1, v0 ¼ 1, lint ¼ 1). The stochastic equations of motion for the individual particles read r_ i ¼ vðφi Þ;
φ_ i ¼
N X
T φ ðrji ;φi ;φj Þ þ
pffiffiffiffiffiffiffiffiffi 2Dφ ζi ðtÞ;
ð1Þ
j≠i
with ri and vðφi Þ ¼ ðcos φi ; sin φi ÞT denoting the position and velocity vector of the ith particle, respectively. Because of the constant speed jvðφi Þj ¼ 1, the velocity vector is fully determined by the polar angle φi representing the direction of motion of the particle. The particles reorient according to pair interactions T φ ðrji ; φi ; φj Þ, which depend on the distance vector rji ¼ rj − ri and the orientation angles of the interaction partners. The last term on the right-hand side of the angle equation (1) represents angular noise with intensity Dφ. Here, ζ i ðtÞ describes Gaussian random processes with zero mean and hζ i ðtÞζ j ðt0 Þi ¼ δij δðt − t0 Þ. The pair-wise interaction is modeled as follows: T φ ðrji ; φi ; φj Þ ¼ þ μðrji Þ sin ðφj − φi Þ − κ sin ðαji − φi ÞΘðξr − rji Þ:
ð2Þ
The first term aligns the interacting particles either parallel or antiparallel depending on the sign of the interaction strength μðrji Þ, which is a function of the scalar distance rji ¼ jrji j. For μðrji Þ > 0, the velocity vectors align, whereas μðrji Þ < 0 implies an antiparallel orientation of these vectors (antialignment). For simplicity, we model μðrji Þ as a piecewise parabolic function ( μðrji Þ ¼
þμþ ð1 − ðrji =ξa Þ2 Þ for 0 ≤ rji ≤ ξa −μ−
4ðrji −ξa Þð1−rji Þ ð1−ξa Þ2
for ξa < rji ≤ 1
;
(a)
week ending 19 DECEMBER 2014 (b)
FIG. 1 (color online). (a) Interaction scheme in the frame of the focal particle. (b) The alignment strength μðrji Þ versus distance, cf. Eq. (3).
κ ≥ 0 is assumed below a distance rji ≤ ξr and no interaction for rji > ξr, indicated by the Heaviside function ΘðxÞ. The repulsive interaction is motivated by steric interaction of finite-sized particles [34]. We assume that short-range alignment and repulsion act on comparable length scale ξr ≲ ξa and use ξr ¼ ξa =2. Langevin simulations.—We performed numerical simulations of the microscopic particle dynamics (1) for high densities ρ0 > 1, which correspond to the situation encountered in dense bacterial suspensions. Furthermore, the local mean-field approximation used below yields a good approximation in this regime [35,36]. Examples of the observed patterns are shown in Fig. 2 (see Supplemental Material (SM) for movies [36]). For low μþ and μ− , we observe a spatially homogeneous, disordered state (region I in the phase diagram Fig. 4, Movie 1 in [36]). For low μ−, an increase in μþ eventually leads to the onset of longrange orientational order (region II in Fig. 4, Movie 2 in [36]). In contrast, for large μþ and large μ− , we observe periodic arrays of vortices in the velocity field (Fig. 2(a), region III in Fig. 4, Movie 3 in [36]). These vortices are typically accompanied by a periodic modulation of the density. However, a strong short-range repulsion (large κ) enforces a quasihomogeneous density profile with only weak fluctuations. Interestingly, for parameters in the
ð3Þ
with μ > 0. Thus, the interaction favors alignment at short distances and antialignment at larger distances within the interaction range (see also Fig. 1). The maximal strengths of alignment and antialignment are denoted by μþ and μ− , respectively. The second term on the right-hand side of Eq. (2) describes a repulsive soft-core interaction where αji ¼ argðrji Þ is the polar positional angle of particle j in the reference frame of particle i. A constant repulsion strength
FIG. 2 (color online). Snapshots of the coarse-grained velocity field v (coarse graining length δl ¼ ξa ¼ 0.2) obtained from numerical solution of the microscopic model (1): (a) spatially periodic vortex state (μþ ¼ 3.2, μ− ¼ 3.2 × 10−2 ). (b) Turbulent state close to the transition from the vortex phase to onset of polar order (μþ ¼ 1.6, μ− ¼ 8 × 10−3 ). Colors indicate the local vorticity ∇∧v ¼ ∂ x vy − ∂ y vx. Simulation parameters: L ¼ 30, N ¼ 135 000, Dφ ¼ 1.0, ξa ¼ 0.2, ξr ¼ 0.1, κ ¼ 10, Δt ¼ 10−2 .
258104-2
PRL 113, 258104 (2014) (a)
(b)
PHYSICAL REVIEW LETTERS
week ending 19 DECEMBER 2014
(c)
FIG. 3 (color online). (a) Radial velocity correlation function ~ vv ðrÞ for a periodic vortex array (μþ ¼ 3.2, μ− ¼ 3.2 × 10−2 ) C and (b) in the turbulent state close to the onset of global polar order (μþ ¼ 1.6, μ− ¼ 8 × 10−3 ). The insets show the Fourier transform E2 ðkÞ of the corresponding two-dimensional velocity correlation function Cvv ðrÞ. (c) The energy spectrum E~ 2 ðkÞ in the turbulent state. The solid lines indicate the scalings observed in dense Bacillus subtilis suspensions (a ≈ 1.7, b ≈ 2.7) and the phenomenological theory [16]. The inset shows the most unstable wave number kmax versus μ− predicted from kinetic theory (solid line) and the position of the maximum of E~ 2 ðkÞ in the stationary state of microscopic simulations (dots). Simulation parameters: L ¼ 30, N ¼ 135 000, Dφ ¼ 1.0, ξa ¼ 0.2, ξr ¼ 0.1, κ ¼ 10, Δt ¼ 10−2 .
vortex phase (III) close to the emergence of polar order (II), the short-range alignment induces mesoscopic convective flows, whereas μ− is just sufficiently large to break global orientational order. These flows destroy the spatial periodicity of the vortex array, and the emergent, irregular pattern is best described as mesoscale turbulence (Fig. 2(b), Movie 4 in [36]), discussed in more detail below. The patterns observed in Langevin simulations of the microscopic model correspond to the ones observed by Dunkel et al. [33] in their phenomenological model. ~ vv ðrÞ of the normalized, Figure 3 shows the radial part C equal time velocity correlation function Cvv ðrÞ ¼ vðr0 ; tÞ · vðr0 þ r; tÞ (overlines indicate temporal averages) as well as the 2 D Fourier transforms E2 ðkÞ of Cvv ðrÞ for a periodic vortex array and in the turbulent state. In the vortex state, we typically observe discrete peaks in E2 ðkÞ, indicating a spatially periodic vortex lattice. The periodicity ~ vv ðrÞ in this is also reflected by damped oscillations of C ~ vv ðrÞ shows only a single regime. In the turbulent state, C minimum and E2 ðkÞ has a ringlike shape indicating a rotational invariant pattern with a characteristic wave number, but without spatial periodicity. In the turbulent regime, the energy spectrum E~ 2 ðkÞ [average of E2 ðkÞ over the wave number angle] is in good agreement with observations made in dense bacterial suspensions [Fig. 3(c)] [11–13,16,17]. However, based on the Langevin simulations, we infer that the observed scalings in the energy spectra are not universal, but depend on the specific choice of parameters. Moreover, our results suggest that mesoscale turbulence emerges due to a secondary instability of the vortex arrays induced by mesoscopic convective flows.
FIG. 4 (color online). Left: Phase diagram as predicted by the kinetic and hydrodynamic theory, based on linear stability analysis of the spatially homogeneous, isotropic state. From kinetic theory, we distinguish three regions: stability of the disordered state (I, white), instability of the spatially homogeneous system (II, blue), and instability at a finite length scale (III, red). Lines represent instabilities predicted by the linearized hydrodynamic theory: solid lines indicate boundaries of region (III); the black dashed line represents boundary for the onset of polar order (α ¼ 0). The labeled points indicate the parameter values corresponding to the snapshots in Fig. 2. Numerical simulations of the model reveal the emergence of polarly ordered states in (II) as well as vortex arrays and mesoscale turbulence in (III). Parameters: Dφ ¼ 1, ξa ¼ 0.2, ξr ¼ 0.1, κ ¼ 10. Right: Exemplary dispersion relations Re½σðkÞ for region II (top, ρ0 μþ ¼ 60, μ− =μþ ¼ 10−3 ) and region III (bottom, ρ0 μþ ¼ 70, μ− =μþ ¼ 5 × 10−2 ).
Kinetic theory.—Following Dean [37], we derived the approximate nonlinear Fokker-Planck equation ∂ t pðr; φ; tÞ ¼ −vðφÞ · ∇pðr; φ; tÞ þ Dφ ∂ 2φ pðr; φ; tÞ ZZ 2 0 0 0 0 0 0 d r dφ T φ ðr ; φ; φ Þpðr; φ; tÞpðr þ r ; φ ; tÞ ; − ∂φ ð4Þ for P the one-particle density function pðr; φ; tÞ ¼ h Ni¼1 δ(r − ri ðtÞ)δ(φ − φi ðtÞ)i. The derivation relies on a mean-field assumption which offers a good approximation at high densities. Details of the derivation are discussed in the SM [36]. Equation (4) is solved by the spatially homogeneous, isotropic state p0 ¼ ρ0 =ð2πÞ where ρ0 is the spatially homogeneous particle density. We study the stability of this solution by linearizing (4) around this stationary state. We obtain the linearized dynamics ∂ t δpðr; φ; tÞ ¼ Lp0 ½δpðr; φ; tÞ for the deviation of the one-particle density from the stationary state δpðr; φ; tÞ ¼ pðr; φ; tÞ − p0 . The eigenvalues σðkÞ, with k ¼ jkj, of the linear operator Lp0 were calculated numerically in Fourier space [36]. They determine the stability of the ground state. The sign convention is chosen such that instabilities are related to positive real parts ReðσÞ. The different behavior of the largest eigenvalues σðkÞ reveals the existence of two distinct instability mechanisms: We observe an instability
258104-3
PHYSICAL REVIEW LETTERS
PRL 113, 258104 (2014)
of the spatially homogeneous system indicated by a maximum of Re½σðkÞ at k ¼ 0 (cf. upper right panel Fig. 4). For other parameter values, we find an additional instability with a maximum of Re½σðkÞ at a finite wave number k ≠ 0 (Turing instability) predicting the emergence of spatial patterns with a characteristic length scale (cf. lower right panel Fig. 4). This particular instability cannot be found in a system with alignment only, because it crucially depends on the presence of sufficiently strong antialignment interaction [42]. In Fig. 4, the parameter regions where the instabilities occur are illustrated by color. Langevin simulations of the microscopic dynamics in the various parameter regimes reveal the existence of global polar order as well as vortex arrays and mesoscale turbulence, respectively. Linearized hydrodynamic theory.—In order to obtain theoretical insights in the mechanism leading to the instabilities described above, we derived approximated equationsR[36] for the relevant observables—the density ρðr; R tÞ ¼ dφpðr; φ; tÞ and the momentum field wðr; tÞ ¼ dφvðφÞpðr; φ; tÞ—via a moment expansion of Eq. (4). At first glance, it is not obvious that other observables, like the nematic tensor [43], are less relevant. However, it is important to note that we consider all observables locally in space. Because of the short-range alignment interaction, the system will exhibit either local disorder or local polar order, irrespective of the global features of the emergent pattern. This is also apparent from the numerical Langevin simulations of the microscopic dynamics (cf. Fig. 2). For the analysis of the linear stability, we consider only the linearized hydrodynamic equations ∂ρ ¼ − ∇ · w; ∂t
ð5aÞ
∂w Δw ¼ − c∇ρ − Dφ w þ πρ0 μˆ Δ w þ ; ∂t 16Dφ
ð5bÞ
where the positive coefficient c ¼ ð1=2 þ πρ0 κξ3r =6Þ is proportional to the compressibility of the active particle gas. The repulsive interaction leads to an increase of the local pressure. The alignment interaction is encoded in the pseudodifferential operator μˆ Δ ¼
∞ X n¼0
Z μn Δn ;
μn ¼
0
∞
drμðrÞ
r2nþ1 : ð6Þ 4n ðn!Þ2
The nonlocality of the alignment interaction is reflected by the infinite series over Laplacians [44]. A comparison of the kinetic theory with results on linear stability of the spatially homogeneous, isotropic state fρðr; tÞ ¼ ρ0 ; w ¼ 0g using Eqs. (5) shows that the hydrodynamic theory qualitatively reproduces the relevant instabilities (see Fig. 4). In particular, the transition lines predicted by the hydrodynamic theory (black lines in
week ending 19 DECEMBER 2014
Fig. 4) are in good quantitative agreement with the results of the kinetic theory for large noise strength, whereas deviations are found for small noise values [45]. So far, we have kept the full interaction operator μˆ Δ defined in Eq. (6). Now, we connect our theory to phenomenological theories as the one by Toner and Tu [31,32] or the one by Dunkel et al. [16,33] truncating P by c the infinite series in Eq. (6) as μˆ Δ;N c ¼ Nn¼0 μn Δn . The stability of the linearized dynamics at short length scales (large wave numbers) requires ð−1ÞN c μN c < 0 to hold. In region III, where antialignment dominates, the lowest order of truncation is given by N c ¼ 2. In this case, the linearized hydrodynamic equation of the momentum field reads ∂w ≈ − c∇ρ þ αw þ ΓΔw þ πρ0 μ2 Δ2 w; ∂t
ð7Þ
where α ¼ πρ0 μ0 − Dφ and Γ ¼ πρ0 μ1 þ ð16Dφ Þ−1 denotes the effective viscosity of the momentum field w. Important insights into the mechanism leading to pattern formation are gained from this equation. First, polar order arises if the coefficient α is positive (right of black dashed line in Fig. 4). Note that this statement only holds locally, since the existence of local polar order does not imply the emergence of global orientational order. Furthermore, Eq. (7) allows us to relate the emergence of patterns with a characteristic length scale [Fig. 3(c)], such as vortex arrays and mesoscale turbulence, to a Turing instability in the velocity field. This instability only occurs if the effective viscosity Γ becomes negative due to the antialignment interaction (necessary condition). This instability mechanism was proposed earlier within a phenomenological theory for vortex arrays and mesoscale turbulence in [33], which is structurally equivalent to (7). From the linear stability analysis of the kinetic or hydrodynamic equations, we cannot predict the actual nature of the emergent pattern. The results of the Langevin simulations, however, allow us to distinguish between vortex arrays and mesoscale turbulence by analyzing correlation functions and energy spectra from the full dynamics of the self-propelled particles (see Fig. 3). Summary.—We propose a simple self-propelled particle model with purely short-range interactions: alignment of close-by neighbors and antialignment with particles at larger distances. It exhibits not only a polar ordered phase, but also periodic vortex arrays and mesoscale turbulence, although the competing interactions do not a priori favor vortex structures. Since our model captures the phenomenology of mesoscale turbulence, we conjecture that this state may emerge due to a competition of short-range alignment and an opposing interaction at larger distances suppressing the emergence of long-range orientational order. On a more general level, the emergence of spatially periodic patterns of collective dynamics in our model is analogous to the Turing mechanism. However, there is a major difference to classical reaction-diffusion systems:
258104-4
PRL 113, 258104 (2014)
PHYSICAL REVIEW LETTERS
Here, the order parameter is the vectorial flow velocity which is related to convective mass transport. This goes beyond the original Turing idea and introduces a novel coupling between order and convection leading to complex spatiotemporal dynamics: vortex arrays and mesoscale turbulence. We thank S. Heidenreich for enlightening discussions and gratefully acknowledge the support by the German Research Foundation via Research Training Group Grant No. GRK 1558 and International Research Training Group Grant No. IRTG 1740. Likewise, we would like to thank three anonymous referees for constructive comments. P. R. acknowledges the hospitality of the Kavli Institute for Theoretical Physics (UCSB) and support in part by the National Science Foundation under Grant No. NSF PHY11-25915.
*
[1] [2]
[3] [4] [5]
[6] [7] [8] [9]
[10] [11] [12] [13] [14] [15]
Corresponding author. [email protected] P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Eur. Phys. J. Spec. Top. 202, 1 (2012). M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013). F. Jülicher, K. Kruse, J. Prost, and J.-F. Joanny, Phys. Rep. 449, 3 (2007). I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 108, 268303 (2012). M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini et al., Proc. Natl. Acad. Sci. U.S.A. 105, 1232 (2008). V. Schaller, C. Weber, C. Semmrich, E. Frey, and A. R. Bausch, Nature (London) 467, 73 (2010). A. M. Menzel and H. Löwen, Phys. Rev. Lett. 110, 055702 (2013). I. H. Riedel, K. Kruse, and J. Howard, Science 309, 300 (2005). Y. Sumino, K. H. Nagai, Y. Shitaka, D. Tanaka, K. Yoshikawa, H. Chaté, and K. Oiwa, Nature (London) 483, 448 (2012). F. Peruani, J. Starruß, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch, and M. Bär, Phys. Rev. Lett. 108, 098102 (2012). A. Sokolov, R. E. Goldstein, F. I. Feldchtein, and I. S. Aranson, Phys. Rev. E 80, 031903 (2009). C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein, and J. O. Kessler, Phys. Rev. Lett. 93, 098103 (2004). H. P. Zhang, A. Be’er, R. S. Smith, E.-L. Florin, and H. L. Swinney, Europhys. Lett. 87, 48011 (2009). T. Ishikawa, N. Yoshida, H. Ueno, M. Wiedeman, Y. Imai, and T. Yamaguchi, Phys. Rev. Lett. 107, 028102 (2011). A. Sokolov and I. S. Aranson, Phys. Rev. Lett. 109, 248109 (2012).
week ending 19 DECEMBER 2014
[16] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Proc. Natl. Acad. Sci. U.S.A. 109, 14308 (2012). [17] K.-A. Liu and L. I, Phys. Rev. E 86, 011924 (2012). [18] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Phys. Rev. Lett. 110, 228102 (2013). [19] A. M. Turing, Phil. Trans. R. Soc. B 237, 37 (1952). [20] S. Kondo and T. Miura, Science 329, 1616 (2010). [21] S. Krinsky and D. Mukamel, Phys. Rev. B 16, 2313 (1977). [22] M. N. Barber, J. Phys. A 12, 679 (1979). [23] J. Oitmaa, J. Phys. A 14, 1159 (1981). [24] J. D. Shore and J. P. Sethna, Phys. Rev. B 43, 3782 (1991). [25] J. D. Shore, M. Holzer, and J. P. Sethna, Phys. Rev. B 46, 11376 (1992). [26] H. Hong and S. H. Strogatz, Phys. Rev. Lett. 106, 054102 (2011). [27] G. Subramanian and D. L. Koch, J. Fluid Mech. 632, 359 (2009). [28] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995). [29] H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud, Phys. Rev. E 77, 046113 (2008). [30] J.-B. Caussin, A. Solon, A. Peshkov, H. Chaté, T. Dauxois, J. Tailleur, V. Vitelli, and D. Bartolo, Phys. Rev. Lett. 112, 148102 (2014). [31] J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998). [32] J. Toner, Phys. Rev. E 86, 031918 (2012). [33] J. Dunkel, S. Heidenreich, M. Bär, and R. E. Goldstein, New J. Phys. 15, 045016 (2013). [34] We have confirmed that the qualitative results do not depend on how the repulsive interactions are implemented in the model. [35] R. Großmann, L. Schimansky-Geier, and P. Romanczuk, New J. Phys. 15, 085014 (2013). [36] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.113.258104 for movies and complementary technical details, which includes Refs. [37–41]. [37] D. S. Dean, J. Phys. A 29, L613 (1996). [38] C. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences, Springer Series in Synergetics (Springer, New York, 2010). [39] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer Verlag, Berlin, 1996). [40] B. N. Belintsev, M. A. Livshits, and M. V. Volkenstein, Z. Phys. B 44, 345 (1981). [41] E. Hernández-García and C. López, Phys. Rev. E 70, 016216 (2004). [42] One may argue that repulsion can induce instabilities at finite wavelengths. However, the instability we report does emerge for κ ¼ 0 as well, i.e., it is not due to repulsion, but rather antialignment interaction. [43] P.-G. De Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon Press, Oxford, 1993), Vol. 23. [44] Any truncation of the infinite series at finite n yields a local approximation of the interaction integral. [45] Additional phase diagrams for other noise values are shown in the Supplemental Material [36].
258104-5
PHYSICAL REVIEW LETTERS
week ending 19 DECEMBER 2014
Vortex Arrays and Mesoscale Turbulence of Self-Propelled Particles Robert Großmann,1,* Pawel Romanczuk,1 Markus Bär,1 and Lutz Schimansky-Geier2 1
Physikalisch-Technische Bundesanstalt Berlin, Abbestraße 2-12, 10587 Berlin, Germany Department of Physics, Humboldt-Universität zu Berlin, Newtonstraße 15, 12489 Berlin, Germany (Received 26 May 2014; revised manuscript received 11 September 2014; published 19 December 2014) 2
Inspired by the Turing mechanism for pattern formation, we propose a simple self-propelled particle model with short-range alignment and antialignment at larger distances. It is able to produce orientationally ordered states, periodic vortex patterns, and mesoscale turbulence, which resembles observations in dense suspensions of swimming bacteria. The model allows a systematic derivation and analysis of a kinetic theory as well as hydrodynamic equations for density and momentum fields. A phase diagram with regions of pattern formation as well as orientational order is obtained from a linear stability analysis of these continuum equations. Microscopic Langevin simulations of self-propelled particles are in agreement with these findings. DOI: 10.1103/PhysRevLett.113.258104
PACS numbers: 47.63.Gd, 05.10.Gg, 87.10.Mn, 87.18.Gh
The term “active matter” refers to nonequilibrium systems of interacting, self-propelled entities which are able to take up energy from their environment and convert it into motion [1,2]. Examples, such as cytoskeletal filaments [3], chemically driven colloids [4], or flocks of birds [5] have recently received a lot of attention in physics, chemistry, and biology. They exhibit a wide range of collective phenomena which are absent in systems at thermodynamic equilibrium, for example large-scale traveling bands and polar clusters [6,7] as well as arrays of vortices [8,9]. In this context, bacteria represent important model systems which have been used to investigate such different aspects as clustering [10] and rheological properties [11] of active matter systems. Recently, irregular vortex structures were experimentally observed in dense suspensions of swimming bacteria [11–17]. In addition, a phenomenological model with a Swift-Hohenberg coupling was proposed which describes the observed behavior including the power spectrum of the bacterial dynamics [16,18]. The spectra reveal pattern formation on finite length scales. The dynamic vortex patterns observed in experiments and simulations are reminiscent of a turbulent state which led the authors to denominate this new phenomenon mesoscale turbulence. In this Letter, we propose a self-propelled particle model which is capable of producing such mesoscopic spatiotemporal patterns based on competing interactions: local alignment interaction at short length scales and antialignment at larger distances. Our model is inspired by the Turing mechanism of short-range activation and long-range inhibition in reaction-diffusion systems sufficient to explain instabilities at finite length scales [19,20]. In this context, activation corresponds to creation of local polar order, whereas antialignment inhibits the emergence of long-range polar order. The effects of related competing interactions were also studied in the Ising model with 0031-9007=14=113(25)=258104(5)
ferromagnetic nearest neighbor and antiferromagnetic nextnearest neighbor coupling [21–25], as well as in the Kuramoto model [26]. We abstain from considering a detailed model of individual swimmers immersed and interacting through a surrounding fluid. However, we note that hydrodynamic interactions of bacteria have been shown to destabilize polar order at high densities [27], and antialignment is the simplest interaction inhibiting local polar order. Thus, our model has to be viewed as a generalization of the seminal Vicsek model [28], with its focus on generic behavior emerging from simple, effective interactions instead of detailed modeling of specific interactions derived from first principles. The original Vicsek model displays surprisingly complex spatiotemporal behavior [29,30], but it does not exhibit mesoscale turbulence or stationary vortex patterns. In contrast, our generalized model including antialignment produces both orientationally ordered states and periodic vortex arrays. Moreover, a region with spatiotemporal irregular behavior of the order parameter field is found between both states. The simulated patterns reproduce most features found in the turbulent state in experiments with dense suspensions of swimming bacteria and corresponding phenomenological macroscopic models [16,18]. The simplicity of our microscopic model makes it analytically tractable and allows us to derive kinetic equations as well as approximated equations for the coarse-grained density and momentum fields. Thus, we can establish a link between microscopic parameters and the hydrodynamic transport coefficients. The resulting hydrodynamic model is a modification of the Toner-Tu equations [31,32] similar to the phenomenological model of bacterial turbulence [16,33]. We show that antialignment interaction may lead to a negative viscosity of the order parameter field. The latter was proposed in [16,33] on a
258104-1
© 2014 American Physical Society
PRL 113, 258104 (2014)
PHYSICAL REVIEW LETTERS
phenomenological basis. Moreover, we analyze the stability of the spatially homogeneous, isotropic state in both frameworks and compare our results to direct simulations of the microscopic dynamics. Model.—We consider N particles moving with constant speed v0 in a two-dimensional system with periodic boundary conditions. Each particle interacts with all neighbors within a finite interaction range lint . We work in natural units such that the particle mass, speed, and interaction range equal one (m ¼ 1, v0 ¼ 1, lint ¼ 1). The stochastic equations of motion for the individual particles read r_ i ¼ vðφi Þ;
φ_ i ¼
N X
T φ ðrji ;φi ;φj Þ þ
pffiffiffiffiffiffiffiffiffi 2Dφ ζi ðtÞ;
ð1Þ
j≠i
with ri and vðφi Þ ¼ ðcos φi ; sin φi ÞT denoting the position and velocity vector of the ith particle, respectively. Because of the constant speed jvðφi Þj ¼ 1, the velocity vector is fully determined by the polar angle φi representing the direction of motion of the particle. The particles reorient according to pair interactions T φ ðrji ; φi ; φj Þ, which depend on the distance vector rji ¼ rj − ri and the orientation angles of the interaction partners. The last term on the right-hand side of the angle equation (1) represents angular noise with intensity Dφ. Here, ζ i ðtÞ describes Gaussian random processes with zero mean and hζ i ðtÞζ j ðt0 Þi ¼ δij δðt − t0 Þ. The pair-wise interaction is modeled as follows: T φ ðrji ; φi ; φj Þ ¼ þ μðrji Þ sin ðφj − φi Þ − κ sin ðαji − φi ÞΘðξr − rji Þ:
ð2Þ
The first term aligns the interacting particles either parallel or antiparallel depending on the sign of the interaction strength μðrji Þ, which is a function of the scalar distance rji ¼ jrji j. For μðrji Þ > 0, the velocity vectors align, whereas μðrji Þ < 0 implies an antiparallel orientation of these vectors (antialignment). For simplicity, we model μðrji Þ as a piecewise parabolic function ( μðrji Þ ¼
þμþ ð1 − ðrji =ξa Þ2 Þ for 0 ≤ rji ≤ ξa −μ−
4ðrji −ξa Þð1−rji Þ ð1−ξa Þ2
for ξa < rji ≤ 1
;
(a)
week ending 19 DECEMBER 2014 (b)
FIG. 1 (color online). (a) Interaction scheme in the frame of the focal particle. (b) The alignment strength μðrji Þ versus distance, cf. Eq. (3).
κ ≥ 0 is assumed below a distance rji ≤ ξr and no interaction for rji > ξr, indicated by the Heaviside function ΘðxÞ. The repulsive interaction is motivated by steric interaction of finite-sized particles [34]. We assume that short-range alignment and repulsion act on comparable length scale ξr ≲ ξa and use ξr ¼ ξa =2. Langevin simulations.—We performed numerical simulations of the microscopic particle dynamics (1) for high densities ρ0 > 1, which correspond to the situation encountered in dense bacterial suspensions. Furthermore, the local mean-field approximation used below yields a good approximation in this regime [35,36]. Examples of the observed patterns are shown in Fig. 2 (see Supplemental Material (SM) for movies [36]). For low μþ and μ− , we observe a spatially homogeneous, disordered state (region I in the phase diagram Fig. 4, Movie 1 in [36]). For low μ−, an increase in μþ eventually leads to the onset of longrange orientational order (region II in Fig. 4, Movie 2 in [36]). In contrast, for large μþ and large μ− , we observe periodic arrays of vortices in the velocity field (Fig. 2(a), region III in Fig. 4, Movie 3 in [36]). These vortices are typically accompanied by a periodic modulation of the density. However, a strong short-range repulsion (large κ) enforces a quasihomogeneous density profile with only weak fluctuations. Interestingly, for parameters in the
ð3Þ
with μ > 0. Thus, the interaction favors alignment at short distances and antialignment at larger distances within the interaction range (see also Fig. 1). The maximal strengths of alignment and antialignment are denoted by μþ and μ− , respectively. The second term on the right-hand side of Eq. (2) describes a repulsive soft-core interaction where αji ¼ argðrji Þ is the polar positional angle of particle j in the reference frame of particle i. A constant repulsion strength
FIG. 2 (color online). Snapshots of the coarse-grained velocity field v (coarse graining length δl ¼ ξa ¼ 0.2) obtained from numerical solution of the microscopic model (1): (a) spatially periodic vortex state (μþ ¼ 3.2, μ− ¼ 3.2 × 10−2 ). (b) Turbulent state close to the transition from the vortex phase to onset of polar order (μþ ¼ 1.6, μ− ¼ 8 × 10−3 ). Colors indicate the local vorticity ∇∧v ¼ ∂ x vy − ∂ y vx. Simulation parameters: L ¼ 30, N ¼ 135 000, Dφ ¼ 1.0, ξa ¼ 0.2, ξr ¼ 0.1, κ ¼ 10, Δt ¼ 10−2 .
258104-2
PRL 113, 258104 (2014) (a)
(b)
PHYSICAL REVIEW LETTERS
week ending 19 DECEMBER 2014
(c)
FIG. 3 (color online). (a) Radial velocity correlation function ~ vv ðrÞ for a periodic vortex array (μþ ¼ 3.2, μ− ¼ 3.2 × 10−2 ) C and (b) in the turbulent state close to the onset of global polar order (μþ ¼ 1.6, μ− ¼ 8 × 10−3 ). The insets show the Fourier transform E2 ðkÞ of the corresponding two-dimensional velocity correlation function Cvv ðrÞ. (c) The energy spectrum E~ 2 ðkÞ in the turbulent state. The solid lines indicate the scalings observed in dense Bacillus subtilis suspensions (a ≈ 1.7, b ≈ 2.7) and the phenomenological theory [16]. The inset shows the most unstable wave number kmax versus μ− predicted from kinetic theory (solid line) and the position of the maximum of E~ 2 ðkÞ in the stationary state of microscopic simulations (dots). Simulation parameters: L ¼ 30, N ¼ 135 000, Dφ ¼ 1.0, ξa ¼ 0.2, ξr ¼ 0.1, κ ¼ 10, Δt ¼ 10−2 .
vortex phase (III) close to the emergence of polar order (II), the short-range alignment induces mesoscopic convective flows, whereas μ− is just sufficiently large to break global orientational order. These flows destroy the spatial periodicity of the vortex array, and the emergent, irregular pattern is best described as mesoscale turbulence (Fig. 2(b), Movie 4 in [36]), discussed in more detail below. The patterns observed in Langevin simulations of the microscopic model correspond to the ones observed by Dunkel et al. [33] in their phenomenological model. ~ vv ðrÞ of the normalized, Figure 3 shows the radial part C equal time velocity correlation function Cvv ðrÞ ¼ vðr0 ; tÞ · vðr0 þ r; tÞ (overlines indicate temporal averages) as well as the 2 D Fourier transforms E2 ðkÞ of Cvv ðrÞ for a periodic vortex array and in the turbulent state. In the vortex state, we typically observe discrete peaks in E2 ðkÞ, indicating a spatially periodic vortex lattice. The periodicity ~ vv ðrÞ in this is also reflected by damped oscillations of C ~ vv ðrÞ shows only a single regime. In the turbulent state, C minimum and E2 ðkÞ has a ringlike shape indicating a rotational invariant pattern with a characteristic wave number, but without spatial periodicity. In the turbulent regime, the energy spectrum E~ 2 ðkÞ [average of E2 ðkÞ over the wave number angle] is in good agreement with observations made in dense bacterial suspensions [Fig. 3(c)] [11–13,16,17]. However, based on the Langevin simulations, we infer that the observed scalings in the energy spectra are not universal, but depend on the specific choice of parameters. Moreover, our results suggest that mesoscale turbulence emerges due to a secondary instability of the vortex arrays induced by mesoscopic convective flows.
FIG. 4 (color online). Left: Phase diagram as predicted by the kinetic and hydrodynamic theory, based on linear stability analysis of the spatially homogeneous, isotropic state. From kinetic theory, we distinguish three regions: stability of the disordered state (I, white), instability of the spatially homogeneous system (II, blue), and instability at a finite length scale (III, red). Lines represent instabilities predicted by the linearized hydrodynamic theory: solid lines indicate boundaries of region (III); the black dashed line represents boundary for the onset of polar order (α ¼ 0). The labeled points indicate the parameter values corresponding to the snapshots in Fig. 2. Numerical simulations of the model reveal the emergence of polarly ordered states in (II) as well as vortex arrays and mesoscale turbulence in (III). Parameters: Dφ ¼ 1, ξa ¼ 0.2, ξr ¼ 0.1, κ ¼ 10. Right: Exemplary dispersion relations Re½σðkÞ for region II (top, ρ0 μþ ¼ 60, μ− =μþ ¼ 10−3 ) and region III (bottom, ρ0 μþ ¼ 70, μ− =μþ ¼ 5 × 10−2 ).
Kinetic theory.—Following Dean [37], we derived the approximate nonlinear Fokker-Planck equation ∂ t pðr; φ; tÞ ¼ −vðφÞ · ∇pðr; φ; tÞ þ Dφ ∂ 2φ pðr; φ; tÞ ZZ 2 0 0 0 0 0 0 d r dφ T φ ðr ; φ; φ Þpðr; φ; tÞpðr þ r ; φ ; tÞ ; − ∂φ ð4Þ for P the one-particle density function pðr; φ; tÞ ¼ h Ni¼1 δ(r − ri ðtÞ)δ(φ − φi ðtÞ)i. The derivation relies on a mean-field assumption which offers a good approximation at high densities. Details of the derivation are discussed in the SM [36]. Equation (4) is solved by the spatially homogeneous, isotropic state p0 ¼ ρ0 =ð2πÞ where ρ0 is the spatially homogeneous particle density. We study the stability of this solution by linearizing (4) around this stationary state. We obtain the linearized dynamics ∂ t δpðr; φ; tÞ ¼ Lp0 ½δpðr; φ; tÞ for the deviation of the one-particle density from the stationary state δpðr; φ; tÞ ¼ pðr; φ; tÞ − p0 . The eigenvalues σðkÞ, with k ¼ jkj, of the linear operator Lp0 were calculated numerically in Fourier space [36]. They determine the stability of the ground state. The sign convention is chosen such that instabilities are related to positive real parts ReðσÞ. The different behavior of the largest eigenvalues σðkÞ reveals the existence of two distinct instability mechanisms: We observe an instability
258104-3
PHYSICAL REVIEW LETTERS
PRL 113, 258104 (2014)
of the spatially homogeneous system indicated by a maximum of Re½σðkÞ at k ¼ 0 (cf. upper right panel Fig. 4). For other parameter values, we find an additional instability with a maximum of Re½σðkÞ at a finite wave number k ≠ 0 (Turing instability) predicting the emergence of spatial patterns with a characteristic length scale (cf. lower right panel Fig. 4). This particular instability cannot be found in a system with alignment only, because it crucially depends on the presence of sufficiently strong antialignment interaction [42]. In Fig. 4, the parameter regions where the instabilities occur are illustrated by color. Langevin simulations of the microscopic dynamics in the various parameter regimes reveal the existence of global polar order as well as vortex arrays and mesoscale turbulence, respectively. Linearized hydrodynamic theory.—In order to obtain theoretical insights in the mechanism leading to the instabilities described above, we derived approximated equationsR[36] for the relevant observables—the density ρðr; R tÞ ¼ dφpðr; φ; tÞ and the momentum field wðr; tÞ ¼ dφvðφÞpðr; φ; tÞ—via a moment expansion of Eq. (4). At first glance, it is not obvious that other observables, like the nematic tensor [43], are less relevant. However, it is important to note that we consider all observables locally in space. Because of the short-range alignment interaction, the system will exhibit either local disorder or local polar order, irrespective of the global features of the emergent pattern. This is also apparent from the numerical Langevin simulations of the microscopic dynamics (cf. Fig. 2). For the analysis of the linear stability, we consider only the linearized hydrodynamic equations ∂ρ ¼ − ∇ · w; ∂t
ð5aÞ
∂w Δw ¼ − c∇ρ − Dφ w þ πρ0 μˆ Δ w þ ; ∂t 16Dφ
ð5bÞ
where the positive coefficient c ¼ ð1=2 þ πρ0 κξ3r =6Þ is proportional to the compressibility of the active particle gas. The repulsive interaction leads to an increase of the local pressure. The alignment interaction is encoded in the pseudodifferential operator μˆ Δ ¼
∞ X n¼0
Z μn Δn ;
μn ¼
0
∞
drμðrÞ
r2nþ1 : ð6Þ 4n ðn!Þ2
The nonlocality of the alignment interaction is reflected by the infinite series over Laplacians [44]. A comparison of the kinetic theory with results on linear stability of the spatially homogeneous, isotropic state fρðr; tÞ ¼ ρ0 ; w ¼ 0g using Eqs. (5) shows that the hydrodynamic theory qualitatively reproduces the relevant instabilities (see Fig. 4). In particular, the transition lines predicted by the hydrodynamic theory (black lines in
week ending 19 DECEMBER 2014
Fig. 4) are in good quantitative agreement with the results of the kinetic theory for large noise strength, whereas deviations are found for small noise values [45]. So far, we have kept the full interaction operator μˆ Δ defined in Eq. (6). Now, we connect our theory to phenomenological theories as the one by Toner and Tu [31,32] or the one by Dunkel et al. [16,33] truncating P by c the infinite series in Eq. (6) as μˆ Δ;N c ¼ Nn¼0 μn Δn . The stability of the linearized dynamics at short length scales (large wave numbers) requires ð−1ÞN c μN c < 0 to hold. In region III, where antialignment dominates, the lowest order of truncation is given by N c ¼ 2. In this case, the linearized hydrodynamic equation of the momentum field reads ∂w ≈ − c∇ρ þ αw þ ΓΔw þ πρ0 μ2 Δ2 w; ∂t
ð7Þ
where α ¼ πρ0 μ0 − Dφ and Γ ¼ πρ0 μ1 þ ð16Dφ Þ−1 denotes the effective viscosity of the momentum field w. Important insights into the mechanism leading to pattern formation are gained from this equation. First, polar order arises if the coefficient α is positive (right of black dashed line in Fig. 4). Note that this statement only holds locally, since the existence of local polar order does not imply the emergence of global orientational order. Furthermore, Eq. (7) allows us to relate the emergence of patterns with a characteristic length scale [Fig. 3(c)], such as vortex arrays and mesoscale turbulence, to a Turing instability in the velocity field. This instability only occurs if the effective viscosity Γ becomes negative due to the antialignment interaction (necessary condition). This instability mechanism was proposed earlier within a phenomenological theory for vortex arrays and mesoscale turbulence in [33], which is structurally equivalent to (7). From the linear stability analysis of the kinetic or hydrodynamic equations, we cannot predict the actual nature of the emergent pattern. The results of the Langevin simulations, however, allow us to distinguish between vortex arrays and mesoscale turbulence by analyzing correlation functions and energy spectra from the full dynamics of the self-propelled particles (see Fig. 3). Summary.—We propose a simple self-propelled particle model with purely short-range interactions: alignment of close-by neighbors and antialignment with particles at larger distances. It exhibits not only a polar ordered phase, but also periodic vortex arrays and mesoscale turbulence, although the competing interactions do not a priori favor vortex structures. Since our model captures the phenomenology of mesoscale turbulence, we conjecture that this state may emerge due to a competition of short-range alignment and an opposing interaction at larger distances suppressing the emergence of long-range orientational order. On a more general level, the emergence of spatially periodic patterns of collective dynamics in our model is analogous to the Turing mechanism. However, there is a major difference to classical reaction-diffusion systems:
258104-4
PRL 113, 258104 (2014)
PHYSICAL REVIEW LETTERS
Here, the order parameter is the vectorial flow velocity which is related to convective mass transport. This goes beyond the original Turing idea and introduces a novel coupling between order and convection leading to complex spatiotemporal dynamics: vortex arrays and mesoscale turbulence. We thank S. Heidenreich for enlightening discussions and gratefully acknowledge the support by the German Research Foundation via Research Training Group Grant No. GRK 1558 and International Research Training Group Grant No. IRTG 1740. Likewise, we would like to thank three anonymous referees for constructive comments. P. R. acknowledges the hospitality of the Kavli Institute for Theoretical Physics (UCSB) and support in part by the National Science Foundation under Grant No. NSF PHY11-25915.
*
[1] [2]
[3] [4] [5]
[6] [7] [8] [9]
[10] [11] [12] [13] [14] [15]
Corresponding author. [email protected] P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Eur. Phys. J. Spec. Top. 202, 1 (2012). M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013). F. Jülicher, K. Kruse, J. Prost, and J.-F. Joanny, Phys. Rep. 449, 3 (2007). I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 108, 268303 (2012). M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini et al., Proc. Natl. Acad. Sci. U.S.A. 105, 1232 (2008). V. Schaller, C. Weber, C. Semmrich, E. Frey, and A. R. Bausch, Nature (London) 467, 73 (2010). A. M. Menzel and H. Löwen, Phys. Rev. Lett. 110, 055702 (2013). I. H. Riedel, K. Kruse, and J. Howard, Science 309, 300 (2005). Y. Sumino, K. H. Nagai, Y. Shitaka, D. Tanaka, K. Yoshikawa, H. Chaté, and K. Oiwa, Nature (London) 483, 448 (2012). F. Peruani, J. Starruß, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch, and M. Bär, Phys. Rev. Lett. 108, 098102 (2012). A. Sokolov, R. E. Goldstein, F. I. Feldchtein, and I. S. Aranson, Phys. Rev. E 80, 031903 (2009). C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein, and J. O. Kessler, Phys. Rev. Lett. 93, 098103 (2004). H. P. Zhang, A. Be’er, R. S. Smith, E.-L. Florin, and H. L. Swinney, Europhys. Lett. 87, 48011 (2009). T. Ishikawa, N. Yoshida, H. Ueno, M. Wiedeman, Y. Imai, and T. Yamaguchi, Phys. Rev. Lett. 107, 028102 (2011). A. Sokolov and I. S. Aranson, Phys. Rev. Lett. 109, 248109 (2012).
week ending 19 DECEMBER 2014
[16] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Proc. Natl. Acad. Sci. U.S.A. 109, 14308 (2012). [17] K.-A. Liu and L. I, Phys. Rev. E 86, 011924 (2012). [18] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Phys. Rev. Lett. 110, 228102 (2013). [19] A. M. Turing, Phil. Trans. R. Soc. B 237, 37 (1952). [20] S. Kondo and T. Miura, Science 329, 1616 (2010). [21] S. Krinsky and D. Mukamel, Phys. Rev. B 16, 2313 (1977). [22] M. N. Barber, J. Phys. A 12, 679 (1979). [23] J. Oitmaa, J. Phys. A 14, 1159 (1981). [24] J. D. Shore and J. P. Sethna, Phys. Rev. B 43, 3782 (1991). [25] J. D. Shore, M. Holzer, and J. P. Sethna, Phys. Rev. B 46, 11376 (1992). [26] H. Hong and S. H. Strogatz, Phys. Rev. Lett. 106, 054102 (2011). [27] G. Subramanian and D. L. Koch, J. Fluid Mech. 632, 359 (2009). [28] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995). [29] H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud, Phys. Rev. E 77, 046113 (2008). [30] J.-B. Caussin, A. Solon, A. Peshkov, H. Chaté, T. Dauxois, J. Tailleur, V. Vitelli, and D. Bartolo, Phys. Rev. Lett. 112, 148102 (2014). [31] J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998). [32] J. Toner, Phys. Rev. E 86, 031918 (2012). [33] J. Dunkel, S. Heidenreich, M. Bär, and R. E. Goldstein, New J. Phys. 15, 045016 (2013). [34] We have confirmed that the qualitative results do not depend on how the repulsive interactions are implemented in the model. [35] R. Großmann, L. Schimansky-Geier, and P. Romanczuk, New J. Phys. 15, 085014 (2013). [36] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.113.258104 for movies and complementary technical details, which includes Refs. [37–41]. [37] D. S. Dean, J. Phys. A 29, L613 (1996). [38] C. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences, Springer Series in Synergetics (Springer, New York, 2010). [39] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer Verlag, Berlin, 1996). [40] B. N. Belintsev, M. A. Livshits, and M. V. Volkenstein, Z. Phys. B 44, 345 (1981). [41] E. Hernández-García and C. López, Phys. Rev. E 70, 016216 (2004). [42] One may argue that repulsion can induce instabilities at finite wavelengths. However, the instability we report does emerge for κ ¼ 0 as well, i.e., it is not due to repulsion, but rather antialignment interaction. [43] P.-G. De Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon Press, Oxford, 1993), Vol. 23. [44] Any truncation of the infinite series at finite n yields a local approximation of the interaction integral. [45] Additional phase diagrams for other noise values are shown in the Supplemental Material [36].
258104-5
