Microscopic theory for the phase separation of self-propelled repulsive disks
eepl draft
Microscopic theory for the phase separation of self-propelled re-pulsive disks
Julian Bialk´e, Hartmut L¨owen, and Thomas Speck
Institut f¨ur Theoretische Physik II, Heinrich-Heine-Universit¨at, D-40225 D¨usseldorf, Germany
PACS – Fluctuation phenomena, random processes, noise, and Brownian motion
PACS – Phase separation and segregation in colloidal systems
Abstract –Motivated by recent experiments on colloidal suspensions, we study analytically andnumerically a microscopic model for self-propelled particles lacking alignment interactions. In thismodel, even for purely repulsive interactions, a dynamical instability leading to phase separationhas been reported. Starting from the many-body Smoluchowski equation, we develop a mean-fielddescription based on a novel closure scheme and derive the effective hydrodynamic equations. Wedemonstrate that the microscopic origin of the instability is a force imbalance due to an anisotropicpair distribution leading to self-trapping. The phase diagram can be understood in terms of twoquantities: a minimal drive and the force imbalance. At sufficiently high propulsion speeds thereis a reentrance into the disordered fluid.
Introduction. –
Living “active matter” [1] rangingfrom bacterial suspensions [2] to flocks of birds [3] is anemerging paradigm on the interface of physics, chemistryand biology. These systems are composed of identical sub-units, which are able to show collective dynamical behav-ior such as swarming [4], active clustering [5], and evenmicro-bacterial turbulence [6]. The fact that every indi-vidual by itself is far from equilibrium allows a virtuallyunlimited number of propagation rules to be devised and,moreover, the emergence of purely dynamical phases. Thisis in contrast to the phenomenon of equilibrium phase sep-aration [7] such as the liquid-vapor transition, which canbe understood on thermodynamic grounds going back tothe seminal work of van der Waals [8]. For one-componentsystems as considered here, equilibrium phase separationrequires attractive interactions.Theoretical progress in statistical physics is often basedon insights gained from the study of models that are sim-ple but still capture an essential feature of real, complexsystems. Self-propelled particles constitute one such classof models that has been very successful in describing non-equilibrium collective behavior. For example, a large classof models [9,10]–the most famous of which is the Vicsek [4]model–describe particles that move with constant velocityand align their orientations with the average orientationof neighboring particles. Such an alignment can also arisedue to volume exclusion as has been observed for granularrods [11]. In the Vicsek model, a dynamical phase transi- tion from an ordered into an unordered phase takes placeas a function of the orientational noise strength.Following earlier work [12, 13], suspensions of self-propelled spherical colloidal particles have recently beenrealized experimentally. The propulsion mechanism is dif-fusiophoresis: either the catalytic decomposition of hydro-gen peroxide on a platinum hemisphere [14] or hematite inconjunction with light [15], or the reversible local demix-ing of a near-critical water-lutidine mixture [16, 17]. Thenon-equilibrium behavior of these driven suspensions isunexpectedly rich and includes a stable fluid of “living”clusters that form, break, and merge but with a steadymean size [15]; and a phase separation, i.e., the largestcluster grows until it is composed of a finite fraction ofthe particles [17]. In a model of run-and-tumble bacteria,Tailleur and Cates have shown theoretically that such adynamical instability resembling liquid-vapor phase sepa-ration can be attributed to a density-dependent mobility(or propulsion speed) [18–20]. Non-equilibrium phase sep-aration has been shown to also occur in computer simula-tions of a minimal model for self-propelled repulsive disks,where particle orientations are independent and can bemodelled as free rotational diffusion [17, 21, 22]. Althoughrecently a link of run-and-tumble motion to active Brown-ian motion has been made [23], a connection between themicroscopic details and the large-scale dynamical instabil-ity leading to phase separation has been missing.In this Letter, we close this gap between numericalp-1 a r X i v : . [ c ond - m a t . s o f t ] J u l ulian Bialk´e, Hartmut L¨owen, and Thomas Specksimulations and experiments of colloidal suspensions onone side and phenomenological models positing a density-dependent mobility on the other side by deriving an ex-plicit expression for the effective swimming speed fromfirst principles. To this end, we start from the full many-body Smoluchowski equation of the minimal model andobtain the effective evolution equation for a single taggedparticle. We employ a novel scheme to close the ensu-ing hierarchy of evolution equations and discuss the ap-proximations involved. The effective swimming speed ofthe tagged particle is reduced compared to its free swim-ming speed due to a force imbalance that is quantified bythe pair distribution function. We finally test our resultsemploying Brownian dynamics computer simulations fordifferent repulsive pair potentials. Model. –
We consider a suspension of N colloidalparticles with number density ¯ ρ moving in two dimensions.The coupled equations of motion are˙ r k = −∇ k U + v e k + ξ k , (1)where the Gaussian noise ξ k models the interactions withthe solvent. The noise has zero mean and temporal short-ranged correlations (cid:104) ξ k ( t ) ξ Tk (cid:48) ( t (cid:48) ) (cid:105) = 2 δ kk (cid:48) δ ( t − t (cid:48) ) . (2)Particles interact through a pair potential u ( r ) with to-tal potential energy U = (cid:80) k The time evolution of the normalizedjoint probability distribution ψ N ( { r k , ϕ k } ; t ) is governedby the Smoluchowski equation ∂ t ψ N = N (cid:88) k =1 ∇ k · [( ∇ k U ) − v e k + ∇ k ] ψ N + D r N (cid:88) k =1 ∂ ψ N ∂ϕ k . (5)The first step is to derive an approximate equation of mo-tion for the projected density ψ ( r , ϕ ; t ) = (cid:90) d r · · · d r N (cid:90) d ϕ · · · d ϕ N N ψ N (6)of a single particle. Since particles are identical, in thefollowing we simply designate particle 1 as the tagged par-ticle and drop the subscript.Performing the integration, Eq. (5) becomes ∂ t ψ = −∇ · [ F + v e ψ − ∇ ψ ] + D r ∂ ϕ ψ (7)with mean force F ( r , ϕ ; t ) ≡ (cid:90) d r · · · d r N (cid:90) d ϕ · · · d ϕ N ( −∇ U ) N ψ N = − (cid:90) d r (cid:48) u (cid:48) ( | r − r (cid:48) | ) r − r (cid:48) | r − r (cid:48) | ψ ( r , ϕ, r (cid:48) ; t ) (8)exerted by the surrounding particles onto the tagged par-ticle. Here, ψ ( r , ϕ, r (cid:48) ; t ) is the two-body density distribu-tion function to find another particle at r (cid:48) (with arbitraryorientation) together with the tagged particle at r withorientation ϕ . Following this scheme leads to the well-known BBGKY hierarchy of coupled equations [24].In order to proceed, we need to find an approximate clo-sure. The two-body density ψ can be decomposed intothe product of the conditional probability g ( | r − r (cid:48) | , θ | r , ϕ )to find another particle at position r (cid:48) given there is a self-propelled particle at r with orientation e times the prob-ability to find a particle at r , ψ ( r , ϕ, r (cid:48) ; t ) = ψ ( r , ϕ ; t )¯ ρg ( | r − r (cid:48) | , θ | r , ϕ ; t ) . (9)Here, θ is the angle enclosed by the displacement vector r (cid:48) − r (pointing away from the tagged particle) and theorientation e .Inserting the decomposition Eq. (9) into Eq. (8), theprojection of the force onto the orientation reads e · F = − ¯ ρζψ with coefficient ζ ≡ (cid:90) ∞ d r r [ − u (cid:48) ( r )] (cid:90) π d θ cos θg ( r, θ ) , (10)where u ( r ) is the pair potential. So far, we have made noapproximations. To proceed, we assume that the coeffi-cient ζ is independent of the position of the tagged par-ticle, which is equivalent to assuming that the system ishomogeneous. Moreover, we neglect the time-dependenceof the pair distribution function, g ( r, θ | r , ϕ ; t ) ≈ g ( r, θ ).These are sufficient conditions to study the onset of ap-2icroscopic theory for the phase separation of self-propelled repulsive diskspossible dynamical instability of an initially homogeneoussuspension. We decompose the mean force in the (non-orthogonal) basis spanned by the orientation e and thedensity gradient ∇ ψ , and expand to linear order of |∇ ψ | with result F = ( e · F ) e + (1 − D ) ∇ ψ + O ( |∇ ψ | ) , (11)where formally D = 1 − F · [ ∇ ψ − e ( e · ∇ ψ )] / |∇ ψ | .Hence, Eq. (7) now reads ∂ t ψ = −∇ · [ v e − D ∇ ] ψ + D r ∂ ϕ ψ , (12)which is the desired closed equation for the temporal evo-lution of the tagged particle density. Here, the effectiveswimming speed v ≡ v − ¯ ρζ (13)enters. In the following, we assume that D is a constantcoefficient that does not depend on the position of thetagged particle. It thus corresponds to the long-time diffu-sion coefficient of the passive suspension ( v = 0). Eq. (12)together with Eq. (13) constitutes our first main result.The physical picture we have in mind is the following:As shown in Fig. 1(a), the propulsion of particles leads toa pair distribution g ( r, θ ) that is anisotropic, i.e., there isa depletion zone behind the tagged particle and an excesszone in front of the particle. The effective swimming speed v is reduced due to this anisotropy. This is intuitivelyclear since particles block each other in a dense suspensiondue to volume exclusion or repulsive interactions. Theprojection coefficient ζ quantifies how strongly the taggedparticle is slowed down by its neighbors as a function ofthe particle density ¯ ρ . Note that for passive particles with v = 0 the pair distribution g ( r ) becomes isotropic and,therefore, ζ = 0 vanishes. Dynamical instability. – To make the connection toprevious work on mobility-induced phase separation [23],we now cast Eq. (12) into the form of the more familiareffective hydrodynamic equations (see also Refs. 25, 26).To this end we introduce the particle density ρ ( r , t ) ≡ (cid:90) π d ϕ ψ ( r , ϕ, t ) (14)and the orientational field p ( r , t ) ≡ (cid:90) π d ϕ e ψ ( r , ϕ, t ) (15)as the first moment of the tagged particle orientation.Plugging in Eq. (12), we find ∂ t ρ ( r , t ) = −∇ · [ v p ( r , t ) − D ∇ ρ ( r , t )] , (16)which couples the temporal evolution of the density to theorientational field. The orientational field in turn couplesto an integral involving ee T = 12 + 12 (cid:18) cos 2 ϕ sin 2 ϕ sin 2 ϕ − cos 2 ϕ (cid:19) . (17) y x (a) (b) -2-1 0 1 2 -2 -1 0 1 2 0 1 2 3 4 5 6 7 8 ρ ζ / v * v /v * Fig. 1: (a) Anisotropic pair distribution function g ( r, θ ) for atagged particle (white circle) at the origin with its orientationpointing along the y -axis (arrow). There is a larger probabilityto find another particle in front of the tagged particle comparedto finding it behind. (Shown is simulation data for hard disksat area fraction φ = 0 . v = 20.) (b) Instabilityregion (shaded) as a function of propulsion speed v and forcecoefficient ζ . Along the dashed line v = 0, the upper half planewould correspond to a negative effective swimming speed v < For a closure, we drop the second term corresponding tothe second moment. Again, this approximation is onlyjustified at the onset of a large-scale instability close tothe homogenous state. We thus arrive at ∂ t p ( r , t ) = − ∇ ( vρ ) + D ∇ p − D r p . (18)The stationary solution of the two effective hydrodynamicequations (16) and (18) is ρ ( r , t ) = ¯ ρ and p ( r , t ) = 0.Clearly, for a constant effective speed v this solution isalways stable.Now suppose that ρ ( r , t ) corresponds to a slowly varyingdensity profile such that within the short range of the pairpotential the density is approximately constant. We canthen replace the global density ¯ ρ in the expression for theeffective swimming speed Eq. (13) by the local density with v ( ρ ) = v − ρζ . Following Ref. 23, we consider the limitof large length scales and time scales much longer than1 /D r . Eq. (18) then implies the quasi-stationary solution p ≈ − D r ∇ ( vρ ) = − v − ρζ D r ∇ ρ, (19)i.e., the orientational field is adiabatically enslaved to thedensity and points along the density gradient. Pluggingthis result back into Eq. (16) leads to the diffusion equa-tion ∂ t ρ = ∇ · (cid:20) D + ( v − ρζ )( v − ρζ )2 D r (cid:21) ∇ ρ ≡ ∇ · D∇ ρ (20)with collective diffusion coefficient D ( ρ ) governing thedensity relaxation. Instability of the homogenous suspen-sion is signaled by D (¯ ρ ) < 0. From this criterion we findthe instability region ζ − (cid:54) ζ (cid:54) ζ + bounded by¯ ρζ ± v ∗ = 34 ( v /v ∗ ) ± (cid:112) ( v /v ∗ ) − P v (a) (b) (c) (d) ϕ = 0.4 ϕ = 0.5 0 4 8 12 16 0 4 8 12 16 ρ ζ / v * v /v * ϕ = 0.4 ϕ = 0.510 -3 -2 -1 -80 -40 0 40 80 L / P L(v /v c -1)L = 70L = 50L = 40 04080120160200 0 10 20 30 D l t v ϕ = 0.4 ϕ = 0.5 Fig. 2: Simulation results for (almost) hard disks interacting via Eq. (24): (a) Size of largest cluster P as a function of propulsionspeed v for two area fractions φ using N = 4900 particles. (b) Finite size scaling using the equilibrium 2D-Ising exponentsfor densities φ = 0 . (cid:4) , v c (cid:39) 50) and φ = 0 . (cid:72) , v c (cid:39) 38) for three different system sizes N = L . (c) Long-time diffusioncoefficient D lt vs . the effective speed v = v − ¯ ρζ . The dashed line shows D lt = v / (2 D r ). (d) Reduced force coefficient ¯ ρζ/v ∗ as a function of reduced speed v /v ∗ , cf. Fig. 1(b). The solid line encloses the instability region. The vertical dashed linesindicate the corresponding critical speeds v c . with minimal speed v ∗ ≡ (cid:112) DD r . (22)Eq. (21) constitutes our second main result. Fig. 1(b)shows the instability region assuming that both thepropulsion speed v and the coefficient ζ can be controlledindependently. Of course, keeping all other parametersfixed, ζ is a function of v and thus describes a curve in thisplot. Our theory is only able to predict the onset of thisdynamical instability. Below we will demonstrate usingcomputer simulations that phase separation ensues fromthis instability. Hence, phase separation is predicted tooccur within the shaded region if ζ ( v ) crosses the bound-ary. A curve ζ ( v ) can enter the instability region for v > v ∗ and, as we will demonstrate below, leave it forlarger speeds corresponding to a reentrance into the dis-ordered fluid phase.From Eq. (12) we can also determine the self-diffusioncoefficient D lt ≡ lim t →∞ t (cid:104) [ r ( t ) − r (0)] (cid:105) = D + v D r , (23)i.e., the long-time diffusion coefficient of the interactingsuspension can be calculated as the diffusion coefficient ofa free particle but employing the reduced swimming speed v instead of v . The swimming speed v can be calculatedvia Eq. (13) and is not a fit parameter as in Ref. 21. Brownian dynamics simulations. – We test ourtheoretical predictions with Brownian dynamics simula-tions. We consider up to N = 4900 particles in a quadraticbox and employ periodic boundary conditions. The equa-tions of motion Eq. (1) are integrated with time step5 · − –10 − depending on the speed v . Motivated byrecent experiments [17], we first study a suspension of (al-most) hard disks employing the WCA potential [27] u ( r ) = (cid:40) ε (cid:110)(cid:0) δr (cid:1) − (cid:0) δr (cid:1) (cid:111) + ε ( r (cid:54) / δ )0 ( r > / δ ) , (24) where we set the potential strength ε = 100. Particleshave diameter unity in reduced units. We set 2 / δ = 1 sothat interactions are only present for overlapping particles.We study two area fractions φ = 0 . φ = 0 . 5, where φ = ¯ ρπ/ 4. Moreover, we assume that the no-slip boundarycondition holds, which determines the rotational diffusioncoefficient as D r = 3 in reduced units.We perform a cluster analysis based on a simple over-lap criterion. Two particles are considered as “bonded” iftheir distance is smaller than unity (they overlap), and acluster is the set of all mutually bonded particles. As a ge-ometrical order parameter, in Fig. 2(a) we plot the meansize P of the largest cluster as a function of the propulsionspeed v . For both densities there is a transition from theunordered, homogeneous phase to an ordered phase, wherethe largest cluster occupies a finite fraction of the system.This drive-induced transition occurs as we increase thespeed v and resembles a second order transition, i.e., theorder parameter P does not jump but increases continu-ously.In order to gain a first insight into the universalityclass of the transition, we have performed finite-size scal-ing [28]. To this end, we have simulated smaller systemswith particles N = L for L = 40 , , 70. We employ L as the relevant system size and not the actual length ofthe simulation box since the order parameter P is basedon the fixed particle size. Defining (cid:15) ≡ v /v c − 1, finite-size scaling predicts the order parameter to behave like P L ( v ) = L − β/ν ˜ P ( L /ν (cid:15) ) with universal scaling function˜ P ( x ). Fig. 2(b) shows that the data for different systemsizes indeed collapses onto a single curve using the equi-librium Ising exponents ν = 1 and β = 1 / v (cid:54) v c . From the collapsewe have estimated the critical speeds v c (cid:39) 50 ( φ = 0 . v c (cid:39) 38 ( φ = 0 . ζ by evaluating the integral Eq. (10). In Fig. 2(c),p-4icroscopic theory for the phase separation of self-propelled repulsive disksthe long-time self-diffusion coefficient D lt is plotted as afunction of the effective speed v = v − ¯ ρζ . Below thetransition v < v c , the diffusion coefficient is indeed verywell described by D lt = D + v / (2 D r ) ≈ v / (2 D r ) as pre-dicted and confirms the validity of Eq. (12). For v > v c we observe that the effective speed v stays approximatelyconstant while the diffusion coefficient D lt grows further.Note that in the phase-separated suspension the global speed v and diffusion coefficient D lt are still calculatedfrom all particles. Although particles are slower while partof a large cluster, there is a constant particle exchange be-tween clusters and dilute phase leading to a monotonouslyincreasing D lt as a function of v . Finally, in Fig. 2(d) weplot the force coefficient as a function of reduced speed v /v ∗ , where v ∗ follows from Eq. (22) with the long-timediffusion coefficient D ( (cid:39) . , . 36 for φ = 0 . , . 5) mea-sured at equilibrium ( v = 0). In agreement with Fig. 2(a),for low speeds the value of ζ corresponds to the homoge-neous suspension. The speeds where ζ crosses the stabilityboundary agree well with the finite-size scaling estimatesof the critical speeds v c .The observed phase separation is a robust phenomenonthat does not depend on the employed pair potential. Fora demonstration we have gathered numerical data for threefurther potentials: (H) soft spheres with a harmonic re-pulsion u ( r ) = ε ( r − for r (cid:54) u ( r ) = 0 other-wise [21], (GCM) the Gaussian core model u ( r ) = εe − r ,and (Y) the screened Coulomb or Yukawa potential u ( r ) = εe − κ ( r − − /r with κ = 5. For the last two potentials weassume the effective particle diameter to be unity. In thissection, we relax the coupling of rotational and transla-tional diffusion by setting D r = 3 · − and neglectingthe translational noise altogether, assuming that the ef-fective noise induced by the propulsion is dominant. Tobe as general as possible, we absorb the strength of thepotential in the modified time unit τ /ε .We study the three potentials for different speeds v at densities: (H) φ = 0 . 7, (GCM) φ = 0 . φ = 0 . 2. Fig. 3(a) shows the structure factor S ( q ) ≡ (cid:80) kk (cid:48) (cid:104) e i q · ( r k − r k (cid:48) ) (cid:105) /N for different speeds v . For all threepotentials, we see an abrupt change for small q from ap-proaching a finite value S (0) to a power law, whereby thelatter signals phase separation. In Fig. 3(b), the force coef-ficient ζ is plotted as a function of the speed v , where opensymbols correspond to the homogeneous, and closed sym-bols to the phase separated suspension as determined fromthe behavior of the structure factor at small q . Since trans-lational noise is neglected, v ∗ (cid:39) v , we observe a reentrance into the fluid phase.Finally, Fig. 3(c) shows a snapshot for the harmonic po-tential (H) at v = 0 . -2 -1 S ( q ) q (a)(b)(c) (d) (H) 10 -1 q(GCM) 10 -1 q (Y) ρ ζ v (H)(Y)(GCM) 0 0.1 0.2 0 0.1 0.2 0.3 0.4 0.5 -0.2 0 0.2 0.4 0.6 0.8 1 Fig. 3: Simulation results neglecting translational diffusion for:(H) soft spheres with harmonic repulsion, (GCM) the Gaussiancore model, and (Y) the Yukawa potential. (a) Structure fac-tors S ( q ) for different speeds v increasing from bottom to top.(b) Force coefficient ζ as a function of the propulsion speed v (note that the unit of time compared to Fig. 2 is 1 / v = 0 . v = 0 . t = 25 quantifying the persistence ofparticle motion with respect to the initial particle orientation. disordered fluid at higher propulsion speeds can be ob-served. The color code corresponds to the quantity α k ≡ e k ( t ) · [ r k ( t + ∆ t ) − r k ( t )] / ( v ∆ t ) (25)measuring the persistence of particle motion over the timeinterval ∆ t with respect to the initial orientation. For∆ t (cid:28) /D r we expect for a free particle α k ∼ α k < Discussion. – What is the origin of the dynamicalinstability? The emerging microscopic picture is the fol-lowing: particles collide and, due to the persistence oftheir swimming motion, they block each other forming asmall cluster [14, 15, 17]. Collisions with other particleslead to a growth of the cluster. For a particle situated inthe rim of the cluster to become free, it has to wait a time ∼ /D r for its orientation to point outward. Depending onp-5ulian Bialk´e, Hartmut L¨owen, and Thomas Speckthe ratio between this waiting time and the collision timecontrolled by the propulsion speed, there is either a stablesteady state corresponding to a homogeneous suspensionwith dynamically evolving clusters, or clusters grow untilthe largest cluster occupies a finite fraction of the system.This microscopic picture is confirmed by the effectivelarge scale evolution equations. The relaxation of a densityfluctuation is dominated by the flux p along the local meanorientation, which is given by p ∼ −∇ ( vρ ). Note thatthe “potential” is the product of swimming speed v anddensity ρ . Now image a small region in which the densityis increased but at the same time the swimming speedis decreased. The density of the surrounding region issmaller but the speed of particles is higher such that theproduct vρ might be larger than that of the dense region.In this case there is a flux towards the denser region andphase separation sets in. Conclusions. – We have studied a dynamical insta-bility in a minimal model for self-propelled repulsive disksfor densities below the freezing transition [29]. The insta-bility results in phase separation and intrinsically requiresthe system to be driven away from thermal equilibrium.It cannot be modelled as an effective isotropic attractionthat shifts an equilibrium coexistence line [30]. We haveshown that instead the relevant mechanism is a force im-balance due to the propulsion, which implies an effectiveswimming speed that depends on the density.Starting from the Smoluchowski equation, we have de-rived a closed equation of motion for the tagged parti-cle density by decomposing the mean force due to thesurrounding particles into two independent contributions:along the particle orientation and along the density gradi-ent. We have employed Brownian dynamics simulations tocalculate the two free parameters entering this equation:the force coefficient ζ and the passive long-time diffusioncoefficient D . Using these values as input, our theory isable to reproduce the single-particle long-time diffusioncoefficient in the homogeneous suspension and to predictthe onset of the instability. For different potentials wehave confirmed numerically that there is also a reentranceinto the disordered fluid at sufficiently high propulsion.Here we have studied a two dimensional system but ourapproach also extends to three dimensions. For the sake ofclarity and brevity, we have neglected hydrodynamic inter-actions. At least pair-wise hydrodynamic interactions areeasily incorporated in our framework and will change thecoefficients D , D r , and ζ . However, the qualitative predic-tions, and in particular the instability diagram Fig. 1(b),will remain unchanged. Models that, e.g., imply a reorien-tation of particles during collisions instead of a persistentmotion (see, e.g., Ref. 31) might move out of the instabilityregion either because v ∗ is increased or the force coefficient ζ is reduced. While we have used computer simulations inorder to calculate ζ , as a next step it will be desirable togo to the next level of the BBGKY hierarchy in order tocalculate the anisotropic pair distribution, and therefore the force coefficient ζ , from first principles. ∗ ∗ ∗ This work has been supported financially by the DFGwithin SFB TR6 (project C3) and by the ERC (advancedgrant INTERCOCOS under project number 267499). REFERENCES[1] S. Ramaswamy, Annu. Rev. Cond. Mat. Phys. , 323(2010).[2] K. Drescher et al. , Proc. Natl. Acad. Sci. U.S.A. ,10940 (2011).[3] A. Cavagna et al. , Proc. Natl. Acad. Sci. U.S.A. ,11865 (2010).[4] T. Vicsek et al. , Phys. Rev. Lett. , 1226 (1995).[5] H. H. Wensink and H. L¨owen, Phys. Rev. E , 031409(2008).[6] H. H. Wensink et al. , Proc. Natl. Acad. Sci. U.S.A. ,14308 (2012).[7] H. Tanaka, J. Phys.: Condens. Matter , R207 (2000).[8] J. van der Waals, Verhandel. Konink. Akad. Weten. Am-sterdam (Sect. 1) , (1893).[9] H. Levine, W.-J. Rappel, and I. Cohen, Phys. Rev. E ,017101 (2000).[10] T. Vicsek and A. Zafeiris, Phys. Rep. , 71 (2012).[11] V. Narayan, S. Ramaswamy, and N. Menon, Science ,105 (2007).[12] W. F. Paxton et al. , J. Am. Chem. Soc. , 13424 (2004).[13] R. Golestanian, T. B. Liverpool, and A. Ajdari, Phys.Rev. Lett. , 220801 (2005).[14] I. Theurkauff et al. , Phys. Rev. Lett. , 268303 (2012).[15] J. Palacci et al. , Science , 936 (2013).[16] I. Buttinoni et al. , J. Phys.: Cond. Matter , 284129(2012).[17] I. Buttinoni et al. , Phys. Rev. Lett. , 238301 (2013).[18] J. Tailleur and M. E. Cates, Phys. Rev. Lett. , 218103(2008).[19] M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J.Tailleur, Proc. Natl. Acad. Sci. U.S.A. , 11715 (2010).[20] A. G. Thompson, J. Tailleur, M. E. Cates, and R. A.Blythe, J. Stat. Mech. P02029 (2011).[21] Y. Fily and M. C. Marchetti, Phys. Rev. Lett. ,235702 (2012).[22] G. S. Redner, M. F. Hagan, and A. Baskaran, Phys. Rev.Lett. , 055701 (2013).[23] M. E. Cates and J. Tailleur, EPL , 20010 (2013).[24] J. Hansen and I. McDonald, Theory of Simple Liquids ,3rd ed. (Academic Press, Amsterdam, 2006).[25] J. Toner and Y. Tu, Phys. Rev. Lett. , 4326 (1995).[26] A. Baskaran and M. C. Marchetti, Proc. Natl. Acad. Sci.U.S.A. , 15567 (2009).[27] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem.Phys. , 5237 (1971).[28] K. Binder, Z. Phys. B , 119 (1981).[29] J. Bialk´e, T. Speck, and H. L¨owen, Phys. Rev. Lett. ,168301 (2012).[30] J. Schwarz-Linek et al. , Proc. Natl. Acad. Sci. U.S.A. ,4052 (2012).[31] S. M. Fielding, arXiv:1210.5464 (2013).,4052 (2012).[31] S. M. Fielding, arXiv:1210.5464 (2013).