Hydrodynamic theory of flocking at a solid-liquid interface: long range order and giant number fluctuations
HHydrodynamic theory of flocking at a solid-liquid interface: long range order andgiant number fluctuations
Niladri Sarkar, ∗ Abhik Basu, † and John Toner ‡ Instituut-Lorentz, Leiden University, P.O. Box 9506, 2300 RA Leiden, The Netherlands Condensed Matter Physics Division, Saha Institute of Nuclear Physics, Calcutta 700064, West Bengal, India Department of Physics and Institute of Theoretical Science,University of Oregon, Eugene, Oregon 97403, USA (Dated: February 5, 2021)We construct the hydrodynamic theory of coherent collective motion (“flocking”) at a solid-liquid interface. The polar order parameter and concentration of a collection of “active” (self-propelled) particles at a planar interface between a passive, isotropic bulk fluid and a solid surfaceare dynamically coupled to the bulk fluid. We find that such systems are stable, and have long-range orientational order, over a wide range of parameters. When stable, these systems exhibit“giant number fluctuations”, i.e., large fluctuations of the number of active particles in a fixed largearea. Specifically, these number fluctuations grow as the 3 / I. INTRODUCTION AND SUMMARY OFRESULTS
Active fluids are fluids containing self-propelled units(” active particles ”), each capable of converting stored orambient free energy into motion. These intrinsically non-equilibrium systems often display orientationally orderedphases, both in living [1–5] and nonliving [6–10] systems.Such phases of active fluids exhibit many phenomenaimpossible in equilibrium orientationally ordered phases(e.g., nematics [11]), among them spontaneous breakingof continuous symmetries in two dimensions [12–15], in-stability in the extreme Stokesian limit [16], and giantnumber fluctuations [17–19].As in equilibrium systems, the ”hydrodynamic” (i.e.,long distance, long-time) behavior of active fluids de-pends crucially on whether or not momentum is con-served. In what has become the standard nomencla-ture of this field, systems lacking momentum conserva-tion due to, e.g., friction with a substrate, are referred toas ”dry” [14, 19, 20], while active fluids with momentumconservation are called ”wet” [21].In this paper, we treat a heretofore unconsidered classof systems that is a natural hybrid of these two cases:a collection of polar active particles at a solid-liquid in-terface. Our study is inspired by experiments on high-density actin motility assays [22], in which highly con-centrated actin filaments on a solid-fluid interface arepropelled by motor proteins. Another experimental real-ization is the work of Bricard et al [23, 24], who stud-ied the emergence of macroscopically directed motion ∗ [email protected] † [email protected], [email protected] ‡ [email protected] in ”Quincke rotators”: a dilute collection of motile col-loids, in contact with a solid substrate, and immersedin a passive bulk electrolytic fluid. The colloids sponta-neously roll, making them self propelled, when a suffi-ciently strong electric field is applied to the electrolyticfluid.These systems differ from both dry and wet active mat-ter, as defined above, by having both friction from theunderlying solid substrate and the long range hydrody-namic interactions due to the overlying bulk passive fluid.Very interestingly, the experiments of [22] found giantfluctuations of the active particle number in a fixed vol-ume, with a standard deviation that scaled with N . inthe ordered phase, where N is the mean number in thevolume. Although the exponent of 0 . N . , which probablylies within experimental error of the 0 . z -direction, as illustratedin Fig. 1. We focus on the extreme Stokesian limit, inwhich the momentum of both the active particles and thebulk fluid is entirely negligible relative to viscous drag.The theory is then used to study the stability and fluc-tuations of a uniformly ordered active polar state of theparticles.Our most surprising result is that this system can ex-hibit a long range ordered polar state, even in the pres-ence of noise. This should be contrasted with both ”wet”active systems, which are generically unstable[16] at very a r X i v : . [ c ond - m a t . s o f t ] F e b Solid Substrate Polar particlesBulk fluid z xy
FIG. 1. (Color online) Schematic diagram of a layer of activepolar particles embedded on a solid substrate surrounded bypassive ambient fluid from the top low Reynolds number, and with equilibrium systems,which cannot develop long ranged orientational order intwo dimensions [25–28].Unlike dry polar active matter, which can also exhibitsuch order[12–15], these systems do so even in a lineartheory . Indeed, such a linear theory proves to providean asymptotically exact long wavelength description ofthese systems, again in contrast to dry active matter.The principal condition for stability of our systems issimilar to the stability criterion for a simple compressiblebulk fluid: the analog in our system of the bulk com-pressibility must be positive, as must our analog of theviscosities.Thus, in contrast to ”wet” active matter in the ”Stoke-sian” limit[15, 16], our ”mixed” system can be genericallystable. Indeed, the requirements for stability are as easilymet for these systems as for an equilibrium fluid.Fluctuations about this uniform ordered state are un-derdamped; that is, they both propagate and decay. Thepropagation is non-dispersive; that is, the wavespeed isindependent of the wavenumber q , as in a simple com-pressible fluid. However, the wavespeed is direction de-pendent, as will be discussed in more detail below.However, the decay rate of these fluctuations scalescompletely differently from that of dry active matter; in-deed, it scales linearly with q , just like the real part. Thismeans that, unlike almost every other physical systemthat exhibits non-dispersive propagating modes[29, 30],the ”quality factor” Q of our system does not diverge aswavenumber q (not to be confused with the quality factor Q ) goes to zero. Instead, it approaches a finite, system-dependent constant, independent of q , as q →
0. Thisunusual damping in turn leads to giant number fluctua-tions of the active particles, such that with the standarddeviation (cid:112) (cid:104) ( N − (cid:104) N (cid:105) ) (cid:105) of the number N of the activeparticles contained in a fixed open area scaling with its average (cid:104) N (cid:105) according to (cid:112) (cid:104) ( N − (cid:104) N (cid:105) ) (cid:105) ∝ (cid:104) N (cid:105) / . (I.1)We also find that, when the active particles at the in-terface are in their ordered state, they ”stir” the pas-sive fluid above them. The components (cid:104) v x ( r ⊥ , z ) (cid:105) and (cid:104) v y ( r ⊥ , z ) (cid:105) of the mean squared velocity thereby inducedin the passive fluid are inversely proportional to the dis-tance z from the solid-fluid interface, whereas (cid:104) v z ( r ⊥ , z ) (cid:105) is inversely proportional to the cube of the distance z ,i.e, it falls of as 1 /z . Here, the x , y , and z axes are re-spectively the direction of the polarization of the activeparticles, the in-plane direction orthogonal to that, andthe normal to the interface, as illustrated in figure (1).Here and throughout this paper, we will use r ⊥ to denotethe two component position vector in real space in theplane of the surface; i.e., r ⊥ ≡ x ˆ x + y ˆ y .These predictions for the passive fluid velocity corre-lations could be tested experimentally by tracking neu-trally buoyant passive tracer particles in the passive fluid.Such particle tracking thus provides a probe of the fluc-tuations of the active particles in the ordered state.These velocity fluctuations in the passive fluid also ex-hibit long ranged spatio-temporal correlations, which de-cay as 1 /t for the projection of the velocity on the inter-face. The projection of the passive fluid velocity along z (perpendicular to the surface) decay as 1 /t .These correlations in turn lead to anomalous diffusionof neutrally buoyant passive particles in the x - and y -direction, with variances of the displacements obeying (cid:104) ( x ( t ) − x (0)) (cid:105) = 2 D x t (cid:20) ln (cid:18) v tz (cid:19) + O (1) (cid:21) , t (cid:28) z D z , (cid:104) ( y ( t ) − y (0)) (cid:105) = 2 D y t (cid:20) ln (cid:18) v tz (cid:19) + O (1) (cid:21) , t (cid:28) z D z , (cid:104) ( x ( t ) − x (0))( y ( t ) − y (0) (cid:105) = 0 , t (cid:28) z D z . (I.2)The cross-correlation function in the third line of (I.2)vanishes because the system is symmetric under inversionof y . Here, v is a system-dependent characteristic speedassociated with the active particles (roughly speaking, itis the speed at which they move over the surface dueto their self-propulsion), z is the initial distance of theneutrally buoyant particle from the surface, and D x , D y are diffusion constants of the same order of magnitude asthe diffusion constant in the z direction.The results (I.2) hold for times t short compared tothe time it takes for the neutrally buoyant particle todiffuse in the z direction a distance comparable to z ;that is, for t (cid:28) z D z , where D z is the diffusion constantin the z -direction. This is a genuine diffusion constant;that is, the motion along z is simple diffusion, with a z -independent diffusion constant. For t (cid:29) z D z , after whichthe neutrally buoyant particle will have diffused in the z direction a distance much larger than z , the behaviorof (cid:104) ( x ( t ) − x (0)) (cid:105) , and (cid:104) ( y ( t ) − y (0)) (cid:105) can be obtainedby replacing z in (I.2) with √ D z t , which is just the rmsdistance the particle will have diffused away from thesurface in that time. Doing so, we find that the motion onthese much longer time scales is still superdiffusive, butwith precisely half the ”superdiffusion” constant. That is, (cid:104) ( x ( t ) − x (0)) (cid:105) = D x t (cid:20) ln (cid:18) v tD z (cid:19) + O (1) (cid:21) , t (cid:29) z D z , (I.3) (cid:104) ( y ( t ) − y (0)) (cid:105) = D y t (cid:20) ln (cid:18) v tD z (cid:19) + O (1) (cid:21) , t (cid:29) z D z , (I.4) (cid:104) ( x ( t ) − x (0))( y ( t ) − y (0) (cid:105) = 0 , t (cid:29) z D z , (I.5)Again, as in (I.2) above, the cross-correlation in the thirdline of (I.5) vanishes because of symmetry under inver-sion of y . Diffusion in the z -direction remains conven-tional. Furthermore, the diffusion constant for diffusionin the z -direction is independent of height z . This set ofpredictions could also be tested experimentally by parti-cle tracking of neutrally buoyant tracer particles in thepassive fluid.Particles denser than the passive fluid, which thereforesediment, will also be affected by this activity inducedflow. We find that particles sedimenting from an initialheight z will, when they reach the surface, be spreadover a region of dimensions (cid:112) (cid:104) ( x ( z = 0) − x ( z = z )) (cid:105) and (cid:112) (cid:104) ( y ( z = 0) − y ( z = z )) (cid:105) in the x and y directions,respectively, with (cid:104) ( x ( z = 0) − x ( z = z )) (cid:105) = 2 D x (cid:18) z v sed (cid:19) ln (cid:18) v v sed (cid:19) , v sed (cid:28) v , (I.6) (cid:104) ( y ( z = 0) − y ( z = z )) (cid:105) = 2 D y (cid:18) z v sed (cid:19) ln (cid:18) v v sed (cid:19) , v sed (cid:28) v , (I.7) (cid:104) ( x ( z = 0) − x ( z = z ))( y ( z = 0) − y ( z = z ) (cid:105) = 0 , v sed (cid:28) v , (I.8)where v is roughly the mean speed of the active parti-cles, and v sed is the speed at which the sedimenting par-ticles sink. This results holds only when the sedimentingspeed v sed (cid:28) v ; in the opposite limit, the lateral dif-fusion of a sedimenting particle becomes conventional;specifically, we find that (cid:104) ( x ( z = 0) − x ( z = z )) (cid:105) = 2 D x v (cid:18) z v (cid:19) , v sed (cid:29) v , (I.9) (cid:104) ( y ( z = 0) − y ( z = z )) (cid:105) = 2 D y v (cid:18) z v (cid:19) , v sed (cid:29) v . (I.10) Once again, these predictions should be readilytestable in particle tracking experiments.We turn now to correlations of the active particlesthemselves. These are determined by the characteristicfrequencies of the fluctuations, we find that their eigen-frequencies ω are complex, and given by ω = Υ( θ q ) q , (I.11)where Υ is a complex function of the angle θ q betweenthe wavevector q and the direction of spontaneous po-larization of the active particles, and is the solution of aparameter and direction dependent quadratic equation:Υ + i Υ[ i ( v p + v ρ ) cos θ q + γ (1 + sin θ q )] − i [ iv p cos θ q + γ (1 + ϕ sin θ q )] v ρ cos θ q − c sin θ q = 0 , (I.12)where v p , v ρ , c , γ , and ϕ are phenomenological param-eters of our model, the first four of which have the di-mensions of speed. The speed c is related to the com-pressibility of the active particles in the same way as thesound speed in a simple equilibrium liquid is, while γ isa damping coefficient that plays roughly the same role inour problem as the viscosities play in a simple fluid[31].The speeds v ρ and v p are fundamentally active param-eters, with no analog in equilibrium fluids, which arisefrom the self-propulsion of the active particles.For stability, we must have Im( ω ) < q , whichmeans Im(Υ( θ q )) < θ q . We show in section(III) that this condition is satisfied provided γ and c areboth >
0, and that φ is sufficiently close to 1. Note thatthe first two conditions are similar to the stability crite-rion for a simple compressible bulk fluid: c > γ > φ has no analog in equilibrium systems,because φ itself has no analog in equilibrium either.Thus, in contrast to ”wet” active matter in the ”Stoke-sian” limit[15, 16], our ”mixed” system can be genericallystable. Indeed, the requirements for stability are almostas easily met for these systems as for an equilibrium fluid.Allowing these fluctuations to be driven by an addi-tive white noise (which will arise from both noisy self-propulsion of the active particles, and thermal fluctua-tions), we find that the polarization, defined as a unitvector pointing along the direction of self-propulsion ofthe active particles, has spatio-temporally Fourier trans-formed correlation functions given by C pp ( q , ω ) ≡ (cid:104)| p y ( q , ω ) | (cid:105) = (cid:18) q (cid:19) F pp (cid:18) (cid:18) ωq (cid:19) , θ q (cid:19) , (I.13)where the scaling function F pp ( u, θ q ) is given by F pp ( u, θ q ) = 2 D p ( u − v ρ cos θ q ) ( u − c + ( θ q )) ( u − c − ( θ q )) + ( u Ψ(1 , θ q ) − v ρ Ψ( ϕ, θ q ) cos θ q ) , (I.14)where the parameter D p , defined more precisely by equa-tion (II.17), is proportional to the variance of the afore-mentioned active particle noise Ψ( ϕ, θ ) ≡ γ (1 + ϕ sin θ ), and the “sound speeds” c ± ( θ q ), defined by the positionsof the peaks in the scaling function versus u [32], are pre-cisely those found for dry active matter in [13–15]; i.e., c ± ( θ q ) = (cid:18) v ρ + v p (cid:19) cos θ q ± (cid:114)
14 ( v ρ − v p ) cos θ q + c sin θ q . (I.15)These peak positions agree with those found in theexperiments of [24] on Quinke rotators.Note that the existence of the scaling form (I.13) im-plies that the ratio of the widths of the peaks in C pp ( q , ω ),plotted versus ω for fixed q , to their positions does notchange as q → ; this is what we meant by our earliercryptic comment that the ”quality factor Q ” becomesindependent of q .In Fig. (2), we show a polar plot of the sound speeds c ± ( θ q ), which are the positions of the peaks in the scaling function F pp when plotted versus the scaling argument u ≡ ωq for fixed direction of propagation θ q . As such,they are the analog in our system of the sound speeds ina simple compressible fluid.The density-density correlation function obeys a simi-lar scaling law: C ρρ ( q , ω ) ≡ (cid:104)| p y ( q , ω ) | (cid:105) = 1 q F ρρ (cid:18) (cid:18) ωq (cid:19) , θ q (cid:19) , (I.16)with a slightly different scaling function F ρρ ( u, θ q ) givenby F ρρ ( u, θ q ) = 2 D p ρ c v ρ sin θ q ( u − c + ( θ q )) ( u − c − ( θ q )) + ( u Ψ(1 , θ q ) − v ρ Ψ( ϕ, θ q ) cos θ q ) . (I.17)The orientation-density correlation function also obeys a similar scaling law: C pρ ( q , ω ) ≡ (cid:104) ρ ( q , ω ) p y ( − q , − ω ) (cid:105) = 1 q F pρ (cid:18) (cid:18) ωq (cid:19) , θ q (cid:19) , (I.18) FIG. 2. (Color online) Polar plot of the sound speeds (IV.4)(black curve). Here we have taken q = 1, v ρ = 1, v p = c = 2, and γ = . q for which the scaling function F pp is plotted infigure (3). where the scaling function for this correlation function is F pρ ( u, θ q ) = 2 D p ρ c v ρ sin θ q ( u − v ρ cos θ q )( u − c + ( θ q )) ( u − c − ( θ q )) + ( u Ψ(1 , θ q ) − v ρ Ψ( ϕ, θ q ) cos θ q ) . (I.19)These three scaling functions are plotted versus thescaling argument u for a fixed direction of propagation θ q in Figure (3). Note that the two peaks in the polar-ization scaling function F pp are exactly the same height,but have different widths. We have chosen to plot thisfigure for a fairly generic direction of propagation. Incontrast, for θ = π/ F pp and F ρρ bothbecome even functions of the scaling argument u , withtwo symmetrically placed peaks at ± c . This impliesthat the correlation functions both become even func-tions of ω , with two symmetrically placed peaks at ± c q .Furthermore, for this direction of propagation, the twocorrelations functions are precisely proportional to eachother; that is, their ratio is a constant, independent ofboth q and ω .Another special direction is θ = 0; i.e., propagation along the polarization. In this case, the asymmetry ismaximized: the polarization scaling function has only one peak, at scaling argument u = v p , meaning the cor-relation function has a single peak at ω = v p q . Thedensity-density scaling function also has a single peak, ata different position, namely at scaling argument u = v ρ ,meaning the correlation function has a single peak at ω = v ρ q .Both of these special directions are, indeed, special: for generic directions of propagation θ (cid:54) = 0 , ± π/
2, the scalingfunctions, and, hence, the ω -dependence of the correla-tion functions, look like figure (3). That is, each scalingfunction has two in general asymmetrically placed peaks,in the same positions for both correlation functions.The way the special limit θ = 0 is approached as θ → F pp , while its height stays the same asthat of the second (i.e., rightmost) peak, has its widthcontinuously vanish as θ →
0, while the density scalingfunction F ρρ has the height of its first peak continuouslyvanish as θ →
0. One can see the tendency towards thislimit in the plot of figure (3) at θ = 27 ◦ .Integrating these spatio-temporally Fourier-transformed correlation functions over all frequencies ω gives surprisingly simple expressions for the equal-timecorrelation functions. For the equal-time polarizationautocorrelation, we obtain C pp ( q ) ≡ (cid:104)| p y ( q , t ) | (cid:105) = (cid:90) ∞−∞ dω π C pp ( q , ω ) = f pp ( θ q ) q , (I.20)where f pp ( θ q ) = D p γ (1 + sin θ q ) × (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) . (I.21) FIG. 3. (Color online) Plot of the scaling functions F pp ( u, θ q )( red curve), F ρρ (( u, θ q ) ( black curve), and F pρ ( u, θ q ) (greencurve) versus scaling argument u ≡ ωq for fixed θ q . Here wehave taken v ρ = 1, v p = c = 2, γ = . D p = . ρ c = 1, γ ρ /γ = 1 / θ = 27 ◦ (the directionindicated by the straight line in figure (2). Here we have defined D ( m , (cid:36), ϕ, θ q ) = (cid:20)(cid:0) θ q (cid:1) − m γ ρ γ [ (cid:36) − (cid:36) − ϕ ) sin θ q ] cos θ q (cid:21) , (I.22)with (cid:36) ≡ v p /v ρ and m ≡ v p /c .For the equal-time density autocorrelation, we obtain C ρρ ( q ) ≡ (cid:104)| ρ ( q , t ) | (cid:105) = f ρρ ( θ q ) q (I.23)where f ρρ ( θ q ) = D p ρ c m (1 + sin θ q ) γD ( m , (cid:36), ϕ, θ q ) . (I.24)Finally, the equal-time cross-correlation function isgiven by C pρ ( q ) ≡ (cid:104) p y ( q , t ) ρ ( − q , t ) (cid:105) = f pρ ( θ q ) q (I.25) where the direction dependence is given by f pρ ( θ q ) = − D p ρ c sin θ q cos θ q m γ ρ γ D ( m , (cid:36), ϕ, θ q ) . (I.26)Fourier transforming these in space leads to equallysimple expressions for the real space, equal-time correla-tion functions: C pp ( r ) ≡ (cid:104) δ ˆ p ( r + R , t ) δ ˆ p ( R , t ) (cid:105) ≈ (cid:104) p y ( r + R , t ) p y ( R , t ) (cid:105) = D p r πγ ( r + x ) (cid:18) γ ρ γ (cid:19) m y x (cid:110) ( r + x ) − m γ ρ γ y [ r ( (cid:36) −
1) + ( (cid:36) − ϕ ) x ] (cid:111) , (I.27) C ρρ ( r ) ≡ (cid:104) δρ ( r + R , t ) δρ ( R , t ) (cid:105) = D p ρ c m πγ r ( r + x ) (cid:110) ( r + x ) − m γ ρ γ y [ r ( (cid:36) −
1) + ( (cid:36) − ϕ ) x ] (cid:111) , (I.28) C pρ ( r ) = − D p ρ c m γ ρ πγ xyr (cid:110) ( r + x ) − m γ ρ γ y [ r ( (cid:36) −
1) + ( (cid:36) − ϕ ) x ] (cid:111) . (I.29)In the above, we have defined δ ˆ p ( r , t ) ≡ ˆ p ( r , t ) − (cid:104) ˆ p ( r , t ) (cid:105) , δρ ( r , t ) ≡ ρ ( r , t ) − ρ , where (cid:104) ˆ p ( r , t ) (cid:105) and ρ arethe spatial averages of the polarization ˆ p and the arealnumber density of the active particles respectively.These predictions could also be tested experimentallyin systems in which the active particles can be imaged,like those of [22, 23].The remainder of this paper is organized as follows. Insection (II), we formulate the linear hydrodynamic the-ory, including the equations of motion, for our system,and calculate the eigenfrequencies of small fluctuationsabout the perfectly ordered state. We thereby demon-strate that the ordered state is dynamically stable. Insection (IV), we calculate the correlation functions im-plied by our equations of motion, and use them to showthat long ranged polar order is robust against noise inthese systems. In section (V), we show that the density-density correlation functions found in section (IV) implygiant number fluctuations (equation (I.1)). We also showthat those giant number fluctuations depend on the shapeof the ”counting volume”, as well as the mean numberof particles it contains, in contrast to equilibrium sys-tems, in which the mean number alone determines thevariance. In section (VI), we calculate the velocity fieldsinduced in the passive fluid by the active particles, andthe correlations of those fields. We also consider the mo-tion of tracer particles in the passive fluid, and therebyderive the anomalous diffusion (equations (I.2)-(I.8).) Fi-nally, in section (VIII), we use a dynamical renormaliza-tion group analysis to show that these linear results areasymptotically exact at long wavelengths. That is, weshow that all non-linear terms allowed in the hydrody-namic equations of motion are irrelevant, in the renor-malization group sense, in the long-wavelength limit. Inappendix (A), we solve the bulk Navier-Stokes equationto relate the passive fluid velocity throughout the bulkto the fluid velocity on the liquid-solid interface. In (B),we derive the stability limit for the ordered state. Inappendix (C), we evaluate an integral that arises in thecalculation of the real-space correlation function. II. DERIVATION OF THE HYDRODYNAMICTHEORYA. Identifying the hydrodynamic variables
We consider active polar particles of areal density ρ ( r ⊥ , t ) confined to a flat interface between a solid sub-strate and a three-dimensional (3D) passive, isotropic,bulk fluid. We choose our coordinates ( x, y, z ) so thatthe interface sits at z = 0, with the solid substrate un-derneath (i.e., z < z > α , where α runs from 1 to the number of parti-cles N ), can have associated with it a unique unit vectorˆ p α . (Imagine, e.g., a collection of arrows, with ˆ p α point-ing from the tail to the head of arrow α .)To formulate hydrodynamics, we coarse grain thisparticle-based unit vector in to a position dependent vec- tor field ˆ p ( r ⊥ , t ) defined on the surface at z = 0. Whilein principle this coarse graining can lead to fluctuationsin the magnitude of the coarse grained ˆ p ( r ⊥ , t ) as well asits direction, such amplitude fluctuations can be shownto be irrelevant[33] in the ordered state, and so we willignore such amplitude fluctuations of p and set p = 1.We choose our coordinate system such that p x = 1 is thereference polar ordered state (see figure (1).In the presence of friction from the substrate, there isno momentum conservation on the surface, so the onlyconserved variable on the surface of relevance for thepresent problem is the active particle number. The hy-drodynamic theory that we develop below thus includesthe active particle number density ρ ( r ⊥ , t ) and the orien-tational order parameter ˆ p ( r ⊥ , t ) fluctuations as the rel-evant slow or ”hydrodynamic” variables. We also includethe bulk fluid velocity v ( r , t ), which is defined throughoutthe semi-infinite three dimensional space above the sur-face, since in that space momentum (which is equivalentto velocity in the limit of an incompressible bulk fluid) isconserved. However, we will work in the Stokesian limit,in which viscous forces dominate inertial ones; this hasthe effect of ”enslaving” the bulk velocity to the hydro-dynamic variables ρ ( r ⊥ , t ) and ˆ p ( r ⊥ , t ) on the surface, aswe will show below. B. Hydrodynamic equations of motion
We will formulate the hydrodynamic equations forthese variables by expanding their equations of motionphenomenologically in powers of fluctuations of bothfields ˆ p and ρ from their mean values, and in spatio-temporal gradients. In so doing, we must respect all sym-metries and conservation laws of the underlying dynam-ics. Because this is a non-equilibrium system, however,additional constraints that would apply in equilibrium,such as detailed balance, do not hold here, and any termconsistent with the symmetries and conservation laws isallowed, and therefore will, in general, be present.We will consider systems with underlying rotationalinvariance in the plane of the surface. That is, we willfocus on the case in which there is no a priori preferreddirection in that plane; the active particles are free to spontaneously break this symmetry, by choosing a direc-tion in which to spontaneously point.With these considerations in mind, we will now de-rive the hydrodynamic equations of motion for thefields ρ ( r ⊥ , t ), ˆ p ( r ⊥ , t ), and the bulk velocity velocity v ( r ⊥ , z, t ).We begin with the active particle density ρ ( r ⊥ , t ).
1. Equation of motion for the active particle density ρ The conservation of the active particles implies thatthey obey a continuity equation: ∂ t ρ + ∇ s · J ρ = 0 , (II.1)where ∇ s ≡ ˆ x ∂/∂x + ˆ y ∂/∂y is the 2D gradient operator,with ˆ x and ˆ y the unit vectors along the x and y axisrespectively. The active particle current J ρ is generatedby friction with the bulk fluid velocity v ( r = r ⊥ , z = 0),plus additional diffusive and active currents. We will phenomenologically expand all of these effects to lead-ing order in powers of the bulk velocity evaluated at thesurface v ( r ⊥ , z = 0), and gradients, keeping all termsconsistent with the underlying rotation invariance. Inpractice, this means we can make the vector J ρ only outof vectors the active particle configuration itself chooses,that is, out of gradients, the surface velocity v s ( r ⊥ , t ),by which we mean the bulk velocity evaluated at the sur-face; i.e., v s ( r ⊥ , t ) ≡ v ( r ⊥ , z = 0 , t ), and the polarizationˆ p ( r ⊥ , t ). These constraints force J ρ to take the form: J ρ ( r ⊥ ) = ρ e ( ρ, | v s | ) v s ( x, y, z = 0) + κ ( ρ, | v s | )ˆ p − D ρ ∇ s ρ − D ρ ˆ p (ˆ p · ∇ s ) ρ − f ρ ( x, y ) . (II.2)to leading order in gradients, where D ρ > D ρ is an anisotropic diffusion co-efficient. The existence of a non-zero D ρ reflects thefact that, once the underlying rotation invariance is bro-ken by the development of a spontaneous polarizationˆ p , diffusion along the direction of ˆ p need not, and, ingeneral, will not, proceed at the same speed as diffusionperpendicular to that direction. The ”effective density” ρ e ( ρ, | v s | ) would simply be ρ itself in a Galilean invari-ant model, but since the solid surface breaks Galileaninvariance, it can in general depend non-linearly on ρ ,and even on the magnitude | v s | of the surface velocity.All of these dependences disappear in the linear model;furthermore, our analysis in section (VIII) of the non-linear model shows that these dependences do not affectany of our predictions for the long distance behavior ofthe system. More precisely, they can all be absorbed intoa suitable, finite ”renormalization” of the parameters ofthe linear theory. In renormalization group jargon, theyare ”irrelevant”. We have only included them in the cur-rent (II.2) to make our starting model more general.The factor κ ( ρ, | v s | ) is an active parameter reflectingthe self-propulsion of the particles through interactionwith the solid substrate. It can also in principle dependon ρ and the magnitude | v s | of the surface velocity. Asfor ρ e , all of these dependences disappear in the linearmodel, and prove to also be irrelevant in the RG sensejust described.For completeness, we have included in this current azero-mean Gaussian white noise f ρ that we take to bedelta-correlated in space and time: (cid:104) f ρi ( x, y, t ) f ρj ( x (cid:48) , y (cid:48) , t (cid:48) ) (cid:105) = 2 D ρ δ sij δ ( x − x (cid:48) ) δ ( y − y (cid:48) ) δ ( t − t (cid:48) ) . (II.3)However, as we will show later, this noise proves to beirrelevant in the long wavelength limit compared to thepolarization noise we will introduce later. The diffusionconstants D ρ and D ρ prove to be irrelevant in the long-wavelength limit as well, as we shall see. Using (II.2) in (II.1), the continuity equation (II.1) canbe written ∂ t ρ = − ˆ p · ∇ s κ ( ρ ) − κ ( ρ ) ∇ s · ˆ p + ∇ s · [ D ρ ∇ s ρ + D ρ ˆ p (ˆ p · ∇ s ) ρ ] − ρ e ( ρ, | v s | ) ∇ s · v s − v s · ∇ s ρ e ( ρ, | v s | ) + ∇ s · f ρ . (II.4)
2. The bulk velocity v ( r ⊥ , z, t ) In calculating the bulk velocity v ( r ⊥ , z, t ), we will as-sume the bulk fluid is in the extreme ”Stokesian” limit, inwhich inertia is negligible relative to viscous drag. Thisshould be appropriate for most systems in which the ac-tive particles are microscopic, since the Reynolds’ num-ber will be extremely low for such particles.The three-dimensional (3D) incompressible bulk veloc-ity field v = ( v i , v z ) , i = x, y satisfies the 3D Stokes’equation η ∇ v α ( r ⊥ , z ) = ∂ α Π( r ⊥ , z ) , (II.5)together with the incompressibility constraint ∇ · v = 0.Here ∇ ≡ ˆ x ∂/∂x + ˆ y ∂/∂y + ˆ z ∂/∂z is the full three-dimensional gradient operator, with ˆ x , ˆ y , and ˆ z as theunit vectors along the x , y , and z axes respectively.We could, in principle, add a random Brownian noiseto this Stokesian force balance equation (II.5) to reflectthermal fluctuations in the non-zero temperature bulkfluid. However, since the bulk fluid is passive - that is, equilibrium , we know that the net effect of such a noiseterm on the surface polarization must be to introducea purely equilibrium random noise in the polarizationequation of motion. So we will instead take the simplerapproach of simply adding such a thermal noise by handto the active noises that will also appear in that equation.In the absence of such a noise term, the noiseless Stokesequation (II.5) can be solved exactly for the bulk velocity v ( r ⊥ , z, t ) in terms of the surface velocity v s ( r ⊥ , t ). If weFourier expand the surface velocity: v s ( r ⊥ , t ) = 1 (cid:112) L x L y (cid:88) q v s ( q , t ) e i q · r ⊥ (II.6)where ( L x , L y ) are the linear dimensions of our (pre-sumed rectangular) surface, then, as we show in Ap-pendix A , the bulk velocity v ( r ⊥ , z, t ) is given by v ( r ⊥ , z, t ) = 1 (cid:112) L x L y (cid:88) q [ v s − z ( q · v s )(ˆ q + i ˆ z )] e − qz + i q · r ⊥ (II.7)This is to be supplemented by the boundary condi-tions at z = ∞ , at which all stresses vanish, and at z = 0, which we will now discuss. For a bulk passive fluidwith a free surface, the boundary condition is given bythe vanishing of the appropriate component of the shearstress. On the other hand, for a bulk fluid resting on asolid surface, one imposes the ”no-slip” boundary condi-tion that the relative velocity between the solid surfaceand the fluid layer in contact vanishes. More generally,if there is a slip velocity between the solid surface andthe liquid layer in contact with the solid surface, then one must impose what in fluid dynamics is known as a partial slip condition on the fluid layer in contact withthe solid substrate. For an active fluid layer at the solid-liquid interface, the fluid layer in contact is active ; hence,there should be active forces acting on this layer. Gener-alizing the partial-slip boundary condition for an activesurface fluid, we impose the following active boundarycondition at z = 0: v ( r ⊥ , z = 0 , t ) ≡ v s ( r ⊥ , t )= v ˆ p ( r ⊥ , z = 0 , t ) + µ τ (cid:107) ( r ⊥ , z = 0 , t ) , (II.8)where v is the spontaneous self-propulsion speed of theactive particles relative to the solid substrate, and theforce density τ (cid:107) parallel to the surface is given by τ (cid:107) ( r ⊥ , z = 0 , t ) = ˆ z · σ bulk ( r ⊥ , z = 0 , t ) + τ surface ( r ⊥ , t ) , (II.9)with σ bulk the bulk stress tensor. This is given by σ bulk ij = η ( ∂ i v j + ∂ j v i ) (II.10)which is simply the viscous stress in the bulk fluid. Thesurface force τ surface ( r ⊥ , t ), when expanded to leadingorder in the polarization ˆ p and gradients, must take theform[34, 35] τ surface i = ζ ( ρ )ˆ p · ∇ S p i + ζ ( ρ ) p i ∇ S · ˆ p + p i ˆ p · ∇ S ζ ( ρ ) − ∂ i P (cid:48) s ( ρ ) . (II.11)In (II.11), the ζ , and ζ represent active stresses, and P (cid:48) s is a surface osmotic pressure (not to be confused withthe bulk pressure Π in the Stokes equation (II.5)), bothof which depend on the density of the active particles. Inserting (II.10) and (II.11) into our expression (II.9)for the for the force density τ (cid:107) parallel to the surface, andusing the result in our partial slip boundary condition(II.8), gives, for i = ( x, y ), v si ( r ⊥ , t ) = v a ( ρ ) p i ( r ⊥ , t ) + ζ ( ρ )ˆ p · ∇ S p i + ζ ( ρ ) p i ∇ S · ˆ p + p i ˆ p · ∇ S ζ ( ρ ) + µη (cid:18) ∂v i ( r ⊥ , z, t ) ∂z (cid:19) z =0 − ∂ i P s , (II.12)where we have defined P s ≡ µP (cid:48) s , and used the fact thatthe velocity normal to the surface vanishes at the sur-face (i.e., v z ( z = 0) = 0) [34]. Note that the velocity v i ( r ⊥ , z, t ) is i ’th component the full bulk velocity.The term proportional to ∂ z v i can be given another in-terpretation: if we define a length a ≡ µη , µη (cid:18) ∂v i ∂z (cid:19) z =0 = a (cid:18) ∂v i ∂z (cid:19) z =0 is simply the velocity of a fluid element a dis-tance a above the surface. Hence, we can alternativelyinterpret this term as modeling active particles of finite thickness a on the surface which are passively convectedby the local fluid velocity at that height. For a systemin thermal equilibrium, v a = 0 = ζ , ( ρ ) = ζ ( ρ ), and(II.12) reduces to the well-known equilibrium partial slipboundary condition.
3. Equation of motion for the polarization ˆ p We now turn to the equation of motion for ˆ p . Asthe active particles are polar, the system is not invariantunder ˆ p → − ˆ p symmetry. This means that terms even0in ˆ p are allowed in our phenomenological expression for ∂ t ˆ p . The most general equation of motion for p k allowedby symmetry is therefore: ∂ t p k = T ki (cid:18) αv si − λ pv ( v s · ∇ s ) p i + h i γ + (cid:18) ν − (cid:19) p j ∂ i v sj + (cid:18) ν + 12 (cid:19) (ˆ p · ∇ s ) v si − λ (ˆ p · ∇ s ) p i − ∂ i P p ( ρ ) + f i (cid:19) , (II.13)where we have defined the transverse projection operator T ki ≡ δ ski − p k p i (II.14)which projects any vector orthogonal to ˆ p . Its presencein (II.13) insures that the fixed length condition | ˆ p | = 1on ˆ p is preserved.The equation of motion (II.13) includes all of the rel-evant nonequilibrium terms allowed by symmetry, in ad-dition to the usual equilibrium terms and the convectivecovariant derivative of ˆ p . Here, the presence of the solidsubstrate underneath breaks any Galilean invariance, andhence a term directly proportional to v i in (II.13) is per-missible: indeed the α -term in (II.13) clearly breaks theGalilean invariance. The presence of this terms implies flow alignment or antialignment of the polarization ˆ p de-pending upon the sign of α in Eq. (II.13). The ”polar-ization pressure” P p ( ρ ) is an additional function of thedensity, independent of the ”osmotic pressure” P s ( ρ ) in-troduced earlier. Furthermore, the λ term represents ac-tive self advection. Finally, the terms proportional to ν are ”flow alignment terms”, identical in form to thosefound in nematic liquid crystals[36].In (II.13), the ”molecular field” h is given by h = − δFδ ˆ p , (II.15)where the Frank free energy is F = 12 (cid:90) d r ⊥ [ K ( ∇ S · ˆ p ) + K | ( ∇ S × ˆ p ) | . (II.16)Here, K and K are respectively the Frank splay andbend elastic modulii.We have also added to the equation of motion (II.13)a white noise f with statistics (cid:104) f i ( r ⊥ , t ) f j ( r (cid:48) ⊥ , t (cid:48) ) (cid:105) = 2 D p δ ij δ ( r ⊥ − r (cid:48) ⊥ ) δ ( t − t (cid:48) ) . (II.17)As discussed earlier, this noise also incorporates equilib-rium contributions from thermal fluctuations of the bulkfluid. These must, by the fluctuation-dissipation theo-rem, be spatiotemporally white, as we have assumed in(II.17), since the equilibrium dynamics of ˆ p is local inspace and time. The actual noise strength D p in (II.17)is larger than the equilibrium value (which is proportionalto k B T ) due to active contributions to the noise.
4. Summary of the equations of motion
Our hydrodynamic model, then, is summarized by theequations of motion (II.4) and (II.13) for ρ and ˆ p , re-spectively, and the solution (II.7) of the Stokes equation(II.5) for the bulk velocity field v ( x, y, z, t )obtained withthe boundary condition (II.12). When considering fluc-tuations, we will also need the noise correlations (II.3)and (II.17) . C. Linearization of the equations of motion
The equations of motion and boundary conditionsfound in the previous section have an obvious spatiallyuniform, steady state solution: ρ ( r ⊥ , t ) = ρ (II.18)ˆ p ( r ⊥ , t ) = ˆ x (II.19) v s ( r ⊥ , t ) = v ˆ x (II.20) v ( r ⊥ , z, t ) = v ˆ x (II.21)where we have defined v ≡ v a ( ρ ) (II.22)and have chosen the ˆ x axis of our coordinate system tobe along the (spontaneously chosen) direction of polar-ization, as illustrated in figure (1).As a first step towards understanding fluctuationsabout this steady state, we will write ρ ( r ⊥ , t ) = ρ + δρ ( r ⊥ , t ) , (II.23)ˆ p ( r ⊥ , t ) = ˆ x (cid:113) − p y ( r ⊥ , t ) + p y ( r ⊥ , t )ˆ y , (II.24) v s ( r ⊥ , t ) = ( v + δv sx ( r ⊥ , t ))ˆ x + v sy ( r ⊥ , t )ˆ y , (II.25)and expand the equations of motion II.4) and (II.13) for ρ and ˆ p , and the boundary condition (II.12), to linearorder in δρ and p y . We will obtain the bulk velocity v ( r ⊥ , z, t ) from the surface velocity v s ( r ⊥ , t ) using oursolution (II.7) of the Stokes equation.We will expand the surface osmotic pressure of theparticles P s ( ρ ) to linear order in δρ : P s ( ρ ) = σδρ . (II.26)1We will also expand the ”surface polarization pressure” P p ( ρ ) to linear order in δρ (and will show later that higherpowers of δρ are irrelevant at long wavelengths): P p ( ρ ) = σ p δρ , (II.27)and likewise expand the ζ ’s in equation (II.12) . To linearorder, it is sufficient to take ζ ( ρ ) = ζ ( ρ ) ≡ ζ , (II.28) ζ ( ρ ) = ζ ( ρ ) ≡ ζ , (II.29) ζ ( ρ ) = ζ ( ρ ) + ¯ ζδρ . (II.30)Finally, we expand v a ( ρ ), κ ( ρ ) and ρ e ( ρ ) to linear orderin δρ : v a ( ρ ) ≡ v + v (cid:48) a δρ , (II.31) κ ( ρ ) ≡ κ + ¯ κδρ , (II.32) ρ e ( ρ ) = ρ + ρ (cid:48) e δρ , (II.33)where κ , ¯ κ , ρ , and ρ (cid:48) e are all constants to linear order.All of the other parameters (i.e., the diffusion constants D ρ and D ρ , v a , α , λ , λ pv , γ , and ν can to linear orderbe replaced by their values at ρ = ρ , | v s | = v , andtreated as constants. Now using our active boundary condition Eq. (II.12),the fluctuating velocity components δv x and v y can beexpressed up to linear order as δv sx = v (cid:48) a δρ + ¯ ζ∂ x δρ + µη ( ∂ z v x ) z =0 + ζ ∂ y p y − σ∂ x δρ , (II.34) v sy = v p y + µη ( ∂ z v y ) z =0 + ζ ∂ x p y − σ∂ y δρ . (II.35)These expressions for δv sx and v sy are implicit equa-tions, since the bulk velocities v x and v y on the righthand side also depend on the surface velocity through(II.7) .We can make them explicit by solving them iteratively.This is most conveniently done in Fourier space, becausewe can then obtain the iterative solution to lowest orderin powers of the wavenumber q .Performing a two-dimensional Fourier transform on(II.34) and (II.35) - that is, Fourier transforming over r ⊥ , and using our solution (II.7) for the bulk velocity interms of the surface velocity, we obtain δv sx ( q , t ) = v (cid:48) a δρ ( q , t ) + i ¯ ζq x δρ ( q , t ) + µηW x ( q , t ) + iζ q y p y ( q , t ) − iσq x δρ ( q , t ) , (II.36) δv sy ( q , t ) = v p y ( q , t ) + µηW y ( q , t ) + iζ q x p y ( q , t ) − iσq y δρ ( q , t ) , (II.37)where we have defined W ( q , t ) ≡ (cid:0) ∂ z (cid:8) [ v s ( q , t ) − z ( q · v s ( q , t ))(ˆ q + i ˆ z )] e − qz (cid:9)(cid:1) z =0 = − q v s ( q , t ) − ( q · v s ( q , t ))ˆ q + W z ( q , t )ˆ z (II.38)In the second equality of (II.38), we have not botheredto evaluate W z ( q , t ), since it does not appear in (II.36)or (II.37).If we now insert our expressions (II.36) and (II.37) for v s ( q , t ) after the second equality in (II.38), we see thatthe leading order in q terms (which are O ( q )) come fromthe v (cid:48) a δρ term in (II.36) and the v p y term in (II.37).Keeping only those terms implies v s ( q , t ) ≈ v (cid:48) a δρ ˆ x + v p y ( q , t )ˆ y . (II.39) Inserting this into (II.38) gives W x = − (cid:18) q + q x q (cid:19) v (cid:48) a δρ − v (cid:18) q x q y q (cid:19) p y , (II.40) W y = − v (cid:48) a (cid:18) q x q y q (cid:19) δρ − v (cid:32) q + q y q (cid:33) p y . (II.41)Inserting these expressions back into (II.36) and (II.37)gives us our final closed form expression for the two com-ponents of the surface velocity, written entirely in termsof p y and δρ :2 δv sx ( q , t ) = (cid:20) v (cid:48) a + i (¯ ζ − σ ) q x − µηv (cid:48) a (cid:18) q + q x q (cid:19)(cid:21) δρ ( q , t ) + (cid:20) iζ q y − µηv (cid:18) q x q y q (cid:19)(cid:21) p y ( q , t ) , (II.42) v sy ( q , t ) = (cid:34) v − µη (cid:32) q + q y q (cid:33) + iζ q x (cid:35) p y ( q , t ) − (cid:20) iσq y + µηv (cid:48) a (cid:18) q x q y q (cid:19)(cid:21) δρ , (II.43)Note the non-analytic character of the µη terms in theseexpressions; this reflects the long-ranged hydrodynamicinteraction between active particles on the surface medi-ated by the bulk passive fluid. Indeed, it is only throughthese terms that the presence of the bulk fluid makes it-self felt. Note also that these terms are real; we’ll see in amoment that this makes them damping terms. They arealso the same order in q as the imaginary ζ and σ terms,which are associated with propagation. This is the originof the peculiar property of our system that damping andpropagation are the same order in wavevector (or, equiv-alently, that the quality factor Q of the normal modes of this system is finite and independent of wavenumber q atsmall q ). This should be contrasted with, e.g., a simplebulk equilibrium fluid, for which the propagating termsare O ( q ), while the damping terms are O ( q ).With these expressions (II.42) and (II.43) for the sur-face velocity in terms of p y and δρ in hand, we can nowderive closed equations of motion for p y and δρ . First,we must linearize the general equations of motion (II.4)and (II.13), using the linearizations (II.27)- (II.33) forthe pressures and ζ ’s, and (II.23) and (II.24) for the den-sity and the polarization. Doing so for the continuityequation (II.4), we obtain ∂ t δρ = − κ ∂ y p y − (¯ κ + v ρ (cid:48) e ) ∂ x δρ − ρ ∇ s · v s + [ D ρ ∇ s ρ + D ρ ∂ x ] δρ + ∇ s · f ρ . (II.44)Fourier transforming equation (II.44), and keeping only terms to leading order in q , we find ∂ t δρ ( q , t ) = − iκ q y p y − i (¯ κ + v ρ (cid:48) e ) q x δρ − iρ q · v s + i q · f ρ . (II.45)Setting k = y in the equation of motion (II.13) for ˆ p ,and linearizing the resulting equation gives ∂ t p y = α ( v sy − v p y ) − ( λ pv v + λ ) ∂ x p y + (cid:18) ν − (cid:19) ∂ y δv sx + (cid:18) ν + 12 (cid:19) ∂ x v sy − σ p ∂ y δρ + 1 γ ( K ∂ y + K ∂ x ) p y + f y . (II.46)Note the appearance of the combination v sy − v p y in the α term on the RHS of this equation. The fact that thisparticular combination appears is not an accident, but,rather, a consequence of rotation invariance: if we workto zeroeth order in gradients in our solution for the sur-face velocity, this term vanishes, as it must, since p y is aGoldstone mode associated with the spontaneous break-ing of the continuous rotational symmetry. It is becauseof this exact cancellation that we needed to evaluate v s to higher order in gradients of p y in (II.36) and (II.37), which is why we had to keep higher order gradient termsin that expression. Note that no such cancellation hap-pened when we calculated the vector W eqn. (II.38)earlier, so there we could truncate the expansion for thesurface velocity at zeroeth order in the gradient term.Fourier transforming equation (II.46), and keepingonly terms to leading order in q , we find that the molec-ular field terms proportional to the Frank constants K , are higher order in q , and so can be dropped, leaving uswith ∂ t p y = α ( v sy − v p y ) − i ( λ pv v + λ ) q x p y + i (cid:18) ν − (cid:19) q y δv sx + i (cid:18) ν + 12 (cid:19) q x v sy − iσ p q y δρ + f y (cid:19) . Now inserting our expressions (II.42) and (II.43) for the surface velocity in terms of p y and δρ into the Fourier3transformed equations of motion (II.45) and (II.47), andkeeping only terms to leading order in q , gives our finalclosed form for the linearized equations of motion: ∂ t δρ = − iv ρ [ q x δρ + ρ c q y p y ] + i q · f ρ , (II.47) ∂ t p y = − iv p q x p y − γ (cid:32) q + q y q (cid:33) p y − (cid:18) γ ρ ρ c (cid:19) (cid:18) q x q y q (cid:19) δρ − iσ t q y δρ + f y , (II.48)where we have defined the characteristic velocities v p ≡ λ pv v + λ − ζ α − (cid:18) ν + 12 (cid:19) v (II.49)and v ρ ≡ ¯ κ + ρ (cid:48) e v + ρ v (cid:48) a , (II.50)the total inverse compressibility σ t ≡ σ p + σα − (cid:18) ν − (cid:19) v (cid:48) a , the characteristic density ρ c ≡ ρ v + κ v ρ , (II.51)and the bulk fluid damping coefficients γ ≡ αv µη and γ ρ ≡ αµηv (cid:48) a ρ c . Both γ and γ ρ have the dimensions ofspeed.Clearly, the linearized equations (II.47) and (II.48) areinvariant under q y → − q y , p y → − p y , ρ → ρ , but notinvariant under q x → − q x . This is equivalent to invari-ance in real space under y → − y, p y → − p y , ρ → ρ , butwith no analogous invariance under x → − x . The effec-tive coefficient γ , when positive (which can be achievedby tuning the signs of the various original model param-eters), serves as the effective damping coefficient in themodel. Interestingly, the effective damping here is O ( q ),which is far stronger than the O ( q ) damping in the lin-earized Toner-Tu model for flocking [13, 14] in the hy-drodynamic limit. Hydrodynamic interactions mediated by the passive, bulk fluid above are responsible for this O ( q ) damping.An alert reader might wonder how hydrodynamic in-teractions mediated by the bulk fluid can dominate fric-tion from the underlying solid substrate, which, after all,is O ( q ). The reason is that friction with the substrate,while playing a very important role (in particular, it isthe principal mechanism limiting the speed v of the ac-tive particles), does not act to suppress fluctuations inthe directions of motion of those particles (or, similarly,the polarization. The leading order damping of such fluc-tuations, which are the Goldstone modes of our problem,come from the hydrodynamic interactions. III. MODE STRUCTURE OF THE LINEARIZEDEQUATIONS AND STABILITY OF THEUNIFORM STATE
Having written down the equations of motion, we willnow analyze the linear stability of the system. We workin polar coordinates q = ( q cos θ q , q sin θ q ), and set thenoises to zero. Assuming a time-dependence of the form ρ, p y ∼ exp( − iωt ), this leads to the eigenvalue conditionon ω : ω + iω [ i ( v p + v ρ ) q cos θ q + γq (1 + sin θ q )] − i [ iv p q cos θ q + γq (1 + ϕ sin θ q )] v ρ q cos θ q − c q sin θ q = 0 , (III.1)where we have defined c ≡ σ t ρ c v ρ (III.2)and ϕ ≡ − γ ρ γ . (III.3) . The eigenfrequencies always scale as q , independentof all other parameters, as can be seen by inserting theansatz ω = Υ( θ q ) q (III.4)into (III.1). This leads to a q independent condition onΥ:4Υ + i Υ[ i ( v p + v ρ ) cos θ q + γ (1 + sin θ q )] − i [ iv p cos θ q + γ (1 + ϕ sin θ q )] v ρ cos θ q − c sin θ q = 0 , (III.5)thereby proving that ω scales like q .For stability, we must have Im( ω ) < q , whichmeans Im(Υ( θ q )) < θ q . We show in appendix(B) that this condition is satisfied for sufficiently small γ ρ (which appears in (III.5) through ϕ ≡ − γ ρ γ ), provided γ and c are both >
0. Note that this condition is similar to the stability criterion for a simple compressible bulkfluid: c > γ > γ ρ is − | (cid:36) − | − (cid:114) ( (cid:36) − + 1 m < γ ρ γ < −| (cid:36) − | + (cid:114) ( (cid:36) − + 1 m , (III.6)where we have defined the ”Mach number” m ≡ v ρ c , (III.7)and the speed ratio (cid:36) ≡ v p v ρ . (III.8)It is easy to see that this condition can always be satisfiedfor sufficiently small γ ρ ; in particular, the allowed regionalways includes γ ρ = 0. That a sufficiently large γ ρ rel-ative to the effective damping γ can lead to instabilitiesis not surprising. A non-zero γ ρ (or, equivalently, v (cid:48) a )implies that different patches of the system move at dif-ferent speeds. Of course, damping tends to homogenizethe density by exchanging particles, thereby reducing thespeed differences. However, if γ ρ is too large, dampingmay not be sufficient to suppress these speed fluctua-tions. The eventual steady state may be nonuniform,leading to a patterned state, though its actual naturecannot be ascertained from the linearized equations ofmotion. Nonlinear amplitude equations would be neces-sary for a full-fledged analysis of the steady state. It will be interesting to explore the relation between these in-stabilities here and the banding instabilities reported in[37, 38].Thus generic underdamped propagating waves withanisotropic, θ q -dependent wavespeed proportional to q are expected for the wide range of parameters satisfyingthe stability condition derived in appendix B. IV. CORRELATION FUNCTIONS ANDROBUSTNESS OF LONG-RANGED ORDERAGAINST NOISE
In the stable region of the parameter space, the corre-lation functions for the system in the steady state maybe calculated from the noise-driven equations of mo-tion. Dropping for now the density force f ρ (we willshow later that it is irrelevant in the long-wavelengthlimit), and solving the linear equations of motion for thespatio-temporally Fourier transformed fields p y ( q , ω ) and δρ ( q , ω ) gives p y ( q , ω ) = i ( ω − v ρ q cos θ q ) f y ( q , ω )( ω − c + ( θ q ) q )( ω − c − ( θ q ) q ) + i ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q , (IV.1) δρ ( q , ω ) = iρ c v ρ q sin θ q f y ( q , ω )( ω − c + ( θ q ) q )( ω − c − ( θ q ) q ) + i ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q , (IV.2)where we have definedΨ( ϕ, θ ) ≡ γ (1 + ϕ sin θ ) , (IV.3) and with the “sound speeds” c ± ( θ q ), defined by the po-5sitions of the peaks in the scaling function versus u [32], precisely those found for dry active matter in [13–15];i.e., c ± ( θ q ) = (cid:18) v ρ + v p (cid:19) cos θ q ± (cid:114)
14 ( v ρ − v p ) cos θ q + c sin θ q . (IV.4)In Fig. (2) we show a polar plot of these sound speeds.Autocorrelating these fields with themselves then givestheir spatio-temporally Fourier transformed correlations: C pp ( q , ω ) ≡ (cid:104)| p y ( q , ω ) | (cid:105) = 2 D p ( ω − v ρ q cos θ q ) ( ω − c + ( θ q ) q ) ( ω − c − ( θ q ) q ) + ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q , (IV.5) C ρρ ( q , ω ) ≡ (cid:104)| δρ ( q , ω ) | (cid:105) = 2 D p ρ c v ρ q sin θ q ( ω − c + ( θ q ) q ) ( ω − c − ( θ q ) q ) + ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q , (IV.6)where we have used (II.17) to obtain the autocorrelation (cid:104)| f y ( q , ω ) | (cid:105) = 2 D p .Similar reasoning gives the cross-correlation function C pρ ( q , ω ) ≡ (cid:104) ρ ( q , ω ) p y ( − q , − ω ) (cid:105) = 2 D p ρ c v ρ q sin θ q ( ω − v ρ q cos θ q )( ω − c + ( θ q ) q ) ( ω − c − ( θ q ) q ) + ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q . (IV.7)Pulling a factor of q out of the numerator of each ofthese expressions, and a factor of q out of their denom-inators, gives the scaling forms (I.13), (I.16), and (I.18), with the scaling functions (I.14) and ((I.17)) for thesecorrelation functions given in the introduction, and plot- ted there in figure (3).We can use these expressions to calculate the equaltime correlation functions; we find for the polarization p y C pp ( q ) ≡ (cid:104)| p y ( q , t ) | (cid:105) = (cid:90) ∞−∞ dω π (cid:104)| p y ( q , ω ) | (cid:105) = D p π (cid:90) ∞−∞ ( ω − v ρ q cos θ q ) dω ( ω − c + ( θ q ) q ) ( ω − c − ( θ q ) q ) + ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q (IV.8)Introducing a new variable of integration s ≡ ω − v ρ q cos θ q Ξ( ϕ, θ q )Ψ(1 , θ q ) q (IV.9)6where we have definedΞ(Φ , θ ) ≡ Ψ(Φ , θ )Ψ(1 , θ ) = 1 + Φ sin θ θ (IV.10)gives C pp ( q ) = D p π Ψ(1 , θ q ) q (cid:90) ∞−∞ ( s + ψ ( ϕ, θ q ) ds ( s − ψ + ( ϕ, θ q )) ( s − ψ − ( ϕ, θ q )) + s , (IV.11)where we have defined ψ ± ( ϕ, θ q ) ≡ c ± ( θ q ) − v ρ cos θ q Ξ( ϕ, θ q )Ψ(1 , θ q ) , (IV.12) and ψ ( ϕ, θ q ) ≡ v ρ cos θ q (Ξ( ϕ, θ q ) − , θ q ) . (IV.13)We will show in appendix (C) that, whenever the stabilitycondition is satisfied, ψ + > ψ − <
0; we will makeuse of these facts later.We also show in appendix C that this integral is equalto (cid:90) ∞−∞ ( s + ψ ( ϕ, θ q ) ds ( s − ψ + ( ϕ, θ q )) ( s − ψ − ( ϕ, θ q )) + s = π (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q (cid:110)(cid:0) θ q (cid:1) − m γ ρ γ [ (cid:36) − (cid:36) − ϕ ) sin θ q ] cos θ q (cid:111) . (IV.14)Using (IV.14), we obtain the expression for the equal- time correlation function: C pp ( q ) = D p q Ψ(1 , θ q ) (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) = D p γq (1 + sin θ q ) (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) , (IV.15)where we have defined D ( m , (cid:36), ϕ, θ q ) = (cid:20)(cid:0) θ q (cid:1) − m γ ρ γ [ (cid:36) − (cid:36) − ϕ ) sin θ q ] cos θ q (cid:21) (IV.16)This expression is easily Fourier transformed in spaceto give the real space correlations C pp ( r ) of the polariza-tion fluctuations: C pp ( r ) ≡ (cid:104) p y ( r + R , t ) p y ( R , t ) (cid:105) = (cid:90) d q (2 π ) e i q · r q C pp ( q ) . (IV.17) Defining θ r as the angle between r and the x -axis, we canrewrite this as7 C pp ( r ) = D p (2 π ) γ (cid:90) π dθ q θ q (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) (cid:90) ∞ dq e iqr cos( θ q − θ r ) . (IV.18)The integral over q is given by (cid:90) ∞ dq e iqr cos( θ q − θ r ) = πδ ( r cos( θ q − θ r ))+ P (cid:20) r cos( θ q − θ r ) (cid:21) where P denotes the Cauchy principal value. Using thisin (IV.18), we have C pp ( r ) = D p πγ (cid:90) π dθ q θ q (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) (cid:18) δ ( r cos( θ q − θ r )) + P (cid:20) πr cos( θ q − θ r ) (cid:21)(cid:19) . (IV.19)The Cauchy principal value term in this expression is oddunder the operation θ q → θ q + π , while θ q is evenunder this operation; hence, the contribution to integral from the Cauchy principal value term vanishes. We aretherefore left with C pp ( r ) = D p πγ (cid:90) π dθ q θ q (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) δ ( r cos( θ q − θ r )) . (IV.20)Using δ ( r cos( θ q − θ r )) = δ ( θ q − θ r + π/
2) + δ ( θ q − θ r − π/ r (cid:12)(cid:12)(cid:12)(cid:12)(cid:16) ∂ cos( θ q − θ r ) ∂θ q (cid:17) θ q = θ r ± π/ (cid:12)(cid:12)(cid:12)(cid:12) = δ ( θ q − θ r + π/
2) + δ ( θ q − θ r − π/ r , (IV.21)where the first equality follows from the familiar identity δ ( f ( x )) = (cid:88) x δ ( x − x ) / | f (cid:48) ( x ) | , (IV.22) where { x } are the roots of f ( x ), (IV.20) becomes C pp ( r ) = D p πγr (cid:90) π dθ q θ q (cid:34) (cid:18) γ ρ γ (cid:19) m sin θ q cos θ q D ( m , (cid:36), ϕ, θ q ) (cid:35) ( δ ( θ q − θ r + π/
2) + δ ( θ q − θ r − π/ D p r πγ ( r + x ) (cid:18) γ ρ γ (cid:19) m y x (cid:110) ( r + x ) − m γ ρ γ y [ r ( (cid:36) −
1) + ( (cid:36) − ϕ ) x ] (cid:111) . (IV.23)We now turn to the equal-time density correlations. Using (IV.6), we have8 C ρρ ( q ) ≡ (cid:104)| ρ ( q , t ) | (cid:105) = (cid:90) ∞−∞ dω π (cid:104)| ρ ( q , ω ) | (cid:105) = D p ρ c v ρ q sin θ q π (cid:90) ∞−∞ dω ( ω − c + ( θ q ) q ) ( ω − c − ( θ q ) q ) + ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q . (IV.24)Making the same change of variables of integration (IV.9) that we made earlier, and using partial fractionsagain, we find C ρρ ( q ) = D p ρ c v ρ q sin θ q iπ [ γq (1 + sin θ q )] ( I − ρ − I + ρ ) , (IV.25)where we have defined I ± ρ ≡ lim (cid:15) → (cid:90) ∞−∞ ds ( s − i(cid:15) − ψ + ( θ q ))( s − i(cid:15) − ψ − ( θ q )) ± i ( s − i(cid:15) ) 1 s − i(cid:15) , (IV.26)where we have introduced the small parameter (cid:15) to reg-ularize the integral.We can now do both integrals by the usual complexcontour techniques. Note that, unlike the integral forthe polarization correlation function C pp , the integrandin (IV.26) converges rapidly enough at infinity that theinfinite semicircle needed to close the path of integrationcontributes nothing to the integral. Note also that thepoles of I ± ρ are those of I ± considered earlier, plus onemore pole at s = i(cid:15) .If we choose (cid:15) >
0, then we can close the integral for I + ρ in the upper half plane, picking up the pole at s = i(cid:15) .This gives I + ρ = 2 iπψ + ψ − , (IV.27)where we have taken the limit (cid:15) → I − ρ inthe lower half plane, finding (for (cid:15) >
0) no poles at all.This implies I − ρ = 0 . (IV.28) Inserting (IV.27) and (IV.28) into (IV.25), we obtain C ρρ ( q ) = − D p ρ c v ρ q sin θ q [ γq (1 + sin θ q )] ψ + ψ − . (IV.29)It is straightforward to check that, had we chosen (cid:15) < I + ρ = 0 and I − ρ = − iπψ + ψ − .It is an equally straightforward algebraic exercise, us-ing the definitions (IV.12) of ψ ± and our expressions(IV.4) for c ± ( θ q ) to show that ψ + ψ − = − D ( m , (cid:36), ϕ, θ q ) v ρ γ sin θ q m [Ψ(1 , θ q )] . (IV.30)Using this in (IV.29) gives C ρρ ( q ) = D p ρ c m (1 + sin θ q ) γqD ( m , (cid:36), ϕ, θ q ) . (IV.31)The above expression can be easily Fourier transformedin space to obtain the real space density correlations C ρρ ( r ) just as we did above to obtain C pp ( r ). We therebyobtain: C ρρ ( r ) ≡ (cid:104) δρ ( r + R , t ) δρ ( R , t ) (cid:105) = D p ρ c m πγ r ( r + x ) (cid:110) ( r + x ) − m γ ρ γ y [ r ( (cid:36) −
1) + ( (cid:36) − ϕ ) x ] (cid:111) . (IV.32)9Finally we work out the expression for the crosscorre- lation function C pρ ( q ) ≡ (cid:104) ρ ( q ) p y ( − q ) (cid:105) = (cid:90) ∞−∞ dω π (cid:104) ρ ( q , ω ) p y ( − q , − ω ) (cid:105) = (cid:90) ∞−∞ dω π D p ρ c v ρ q sin θ q ( ω − v ρ q cos θ q )( ω − c + ( θ q ) q ) ( ω − c − ( θ q ) q ) + ( ω Ψ(1 , θ q ) − v ρ q Ψ( ϕ, θ q ) cos θ q ) q = D p ρ c v ρ q sin θ q πγ q (1 + sin θ q ) (cid:90) ∞−∞ ( s + ψ ( ϕ, θ q )) ds ( s − ψ + ) ( s − ψ − ) + s . (IV.33)From equations (C9) and (C10), we can easily show that (cid:90) ∞−∞ ( s + ψ ( ϕ, θ q )) ds ( s − ψ + ) ( s − ψ − ) + s = I + ψ ( ϕ, θ q ) I (IV.34)where I = 0 , and I = − πψ + ψ − . (IV.35)Using eqs. (IV.35) and (IV.33), we obtain C pρ ( q ) = − D p ρ c sin θ q cos θ q m γ ρ γ qD ( m , (cid:36), ϕ, θ q ) . (IV.36)Fourier transformation of the above expression givesthe real space cross correlation, given by C pρ ( r ) = − D p ρ c m γ ρ πγ xyr (cid:110) ( r + x ) − m γ ρ γ y [ r ( (cid:36) −
1) + ( (cid:36) − ϕ ) x ] (cid:111) . (IV.37) V. GIANT NUMBER FLUCTUATIONS ANDTHEIR SHAPE DEPENDENCE
We will show in this section that the long-ranged den-sity correlations found in equation (IV.32) of the previ-ous section lead to ”giant number fluctuations”[39][15,17, 19]. These can be defined as follows:Consider a ”counting box”, defined as a rectangu-lar area A = L x × L y , and define the aspect ratio α A = L x /L y . Our experiment will consist of counting the number of active particles in this box. The mean number of particles N in the box is, of course, given by N = L x L y ρ = α A L y ρ . (V.1)The total number fluctuation δN ≡ N − N in the area A is δN = (cid:90) d r δρ ( r , t ) . (V.2)Thence, (cid:104) ( δN ) (cid:105) = (cid:90) L x dx (cid:90) L y dy (cid:90) L x dx (cid:48) (cid:90) L y dy (cid:48) (cid:104) ρ ( r , t ) ρ ( r (cid:48) , t ) (cid:105) = (cid:90) L x dx (cid:90) L y dy (cid:90) L x dx (cid:48) (cid:90) L y dy (cid:48) C ρρ ( r − r (cid:48) )= D p ρ c m πγ (cid:90) L x dx (cid:90) L y dy (cid:90) L y dy (cid:48) (cid:90) L x dx (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − . (V.3)We consider the limits α A (cid:29) α A (cid:28)
1. Consider first the limit L x (cid:29) L y , i.e., α A (cid:29)
1. Now0consider the integral over x (cid:48) for fixed r = ( x, y ) in thislimit. Note first that in this limit, for most of the range of integration over x , x (cid:29) L y . For such values of x , wecan split the integral over x (cid:48) into three parts: (cid:90) L x dx (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − = I + I + I , (V.4)where we have defined I ≡ (cid:90) x − CL y dx (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − , (V.5) I ≡ (cid:90) x + CL y x − CL y dx (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − , (V.6) I ≡ (cid:90) L x x + CL y dx (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − , (V.7)where C is an arbitrary constant enough larger than1 that we can safely neglect terms of O ( C ). Withthis choice, throughout the region of integration for I , x − x (cid:48) (cid:29) L y . Therefore, throughout this region, | y − y (cid:48) | (cid:28) x − x (cid:48) (since | y − y (cid:48) | is always < L y ). Thisgives (cid:16) y − y (cid:48) x − x (cid:48) (cid:17) (cid:28)
1. Making this approximation for I ,our expression for it reduces to I ≈ (cid:90) x − CL y dx (cid:48) x − x (cid:48) . (V.8)Performing this elementary integral gives I = 12 ln (cid:18) xCL y (cid:19) = 12 ln (cid:18) xL y (cid:19) −
12 ln
C . (V.9) Virtually identical reasoning can be applied to I , givingthe result I = 12 ln (cid:18) L x − xCL y (cid:19) . (V.10)For I , a simple shift of variables of integration x (cid:48) = x + x (cid:48)(cid:48) shows that I is independent of x : I = (cid:90) CL y − CL y dx (cid:48)(cid:48) (cid:112) x (cid:48)(cid:48) + ( y − y (cid:48) ) ( y − y (cid:48) ) + 2 x (cid:48)(cid:48) × (cid:34) − γ ρ v ρ γc ( v p − v ρ )( y − y (cid:48) ) ( y − y (cid:48) ) + 2 x (cid:48)(cid:48) − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) x (cid:48)(cid:48) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] (cid:35) − . (V.11)1Inserting these results (V.9), (V.11), and (V.10) into our expression (V.4), and using that result in our expression(V.3) for (cid:104) ( δN ) (cid:105) , we obtain (cid:104) ( δN ) (cid:105) ≈ D p ρ c m πγ (cid:40) (cid:90) L x dx (cid:90) L y dy (cid:90) L y dy (cid:48) (cid:18) ln (cid:18) xL y (cid:19) + ln (cid:18) L x − xL y (cid:19) − C (cid:19) + (cid:90) L x dx (cid:90) L y dy (cid:90) L y dy (cid:48) (cid:90) CL y − CL y dx (cid:48)(cid:48) (cid:112) x (cid:48)(cid:48) + ( y − y (cid:48) ) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] × (cid:34) − γ ρ v ρ γc ( v p − v ρ )( y − y (cid:48) ) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) x (cid:48)(cid:48) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] (cid:35) − (cid:41) , (V.12)The integrals over y and y (cid:48) in the (triple) integral inthe first line of this expression trivially give a factor of L y , since the integrand is independent of y and y (cid:48) . Theremaining integral over x is elementary. The net resultis12 (cid:90) L x dx (cid:90) L y dy (cid:90) L y dy (cid:48) (cid:18) ln (cid:18) xL y (cid:19) + ln (cid:18) L x − xL y (cid:19) − C (cid:19) = L x L y (cid:18) ln (cid:18) L x CL y (cid:19) − (cid:19) . (V.13)The x integral in the remaining (quadruple) integralwhich appears on the second line in equation (V.12) canbe done immediately, since the integrand is independentof x , yielding a factor of L x . The remaining triple inte-gral over y , y (cid:48) , and x (cid:48)(cid:48) can be done by changing variables of integration to new, rescaled variables u x , u y , and u (cid:48) y via y ≡ L y u y , y (cid:48) ≡ L y u (cid:48) y , x (cid:48)(cid:48) ≡ L y u x . (V.14)This gives (cid:90) L x dx (cid:90) L y dy (cid:90) L y dy (cid:48) (cid:90) CL y − CL y dx (cid:48)(cid:48) (cid:112) x (cid:48)(cid:48) + ( y − y (cid:48) ) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] × (cid:34) − γ ρ v ρ γc ( v p − v ρ )( y − y (cid:48) ) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) x (cid:48)(cid:48) [( y − y (cid:48) ) + 2 x (cid:48)(cid:48) ] (cid:35) − = C (cid:48) L x L y , (V.15)where C (cid:48) ≡ (cid:90) du y (cid:90) du (cid:48) y (cid:90) C − C du x (cid:113) u x + ( u y − u (cid:48) y ) ( u y − u (cid:48) y ) + 2 u x (cid:34) − γ ρ v ρ γc ( v p − v ρ )( u y − u (cid:48) y ) u x + ( u y − u (cid:48) y ) − (cid:18) v ρ γ ρ c γ (cid:19) ( u y − u (cid:48) y ) u x [2 u x + ( u y − u (cid:48) y ) ] (cid:35) − (V.16)is an O (1) constant. Comparing (V.13) and (V.15), wesee that the first line of (V.12) actually dominates the second in the large aspect ratio limit L x (cid:29) L y that weare considering here. Therefore, we obtain, in the limitof large aspect ratio ( α A (cid:29) (cid:112) (cid:104) ( δN ) (cid:105) ≈ ρ c m (cid:112) D p √ πγ (cid:115) L x L y ln (cid:18) L x L y (cid:19) ≈ ρ c m (cid:112) D p √ πγ L / y (cid:112) α A ln α A , (V.17)This expression (V.17) can be rewritten in terms ofthe mean particle number N in the same area A , whichis given by N = L x L y ρ = α A L y ρ . (V.18)This gives L y = ( N / ( α A ρ )) / . (V.19)Using this result (V.19) in our expression (V.17) for thevariance (∆ N ) of the number fluctuations gives∆ N = (cid:112) (cid:104) ( δN ) (cid:105) = (cid:32) ρ c mρ / (cid:33) (cid:32)(cid:115) D p πγ (cid:33) (cid:32) N / √ ln α A α / A (cid:33) . (V.20)Three points should be noted about this result: 1) The number fluctuations are giant ; that is, they growmuch more rapidly with ¯ N than the usual √ ¯ N ”law oflarge numbers” fluctuations, which are found in almostall systems, and, in particular, in most equilibrium sys-tems away from fixed points[40]. Specifically, we have∆ N ∝ N / . Note also that the exponent 3 / . N , but also on the shape (i.e., onthe aspect ratio α A ).3) The number fluctuations ∆ N are a monotonically de-creasing function of the aspect ratio α A in this range of α A (cid:29) L y (cid:29) L x , i.e., aspect ratio α A (cid:28)
1. We again begin with Eq. (V.3). Now considerthe integral over y (cid:48) for fixed r = ( x, y ) in this limit. Notefirst that in this limit, for most of the range of integrationover | y − y | (cid:29) L x . For such values of y , in analogy withour approach for α A (cid:29)
1, we can split the integral over y (cid:48) into three parts (cid:90) L y dy (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ ) ( y − y (cid:48) )( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − = I + I + I , (V.21)where we have defined I ≡ (cid:90) y − ˜ CL x dy (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − , (V.22) I ≡ (cid:90) y + ˜ CL x y − ˜ CL x dy (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − , (V.23)and I ≡ (cid:90) L y y + ˜ CL x dy (cid:48) | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) × (cid:34) − ( v p − v ρ )( y − y (cid:48) ) ( r − r (cid:48) ) + ( x − x (cid:48) ) γ ρ v ρ γc − (cid:18) v ρ γ ρ c γ (cid:19) ( y − y (cid:48) ) ( x − x (cid:48) ) { ( r − r (cid:48) ) + ( x − x (cid:48) ) } (cid:35) − . (V.24)3Once again, we have chosen the constant ˜ C to be suffi-ciently large compared to 1 that, throughout the regionsof integration of I and I , | y − y (cid:48) | (cid:29) L x . Hence, in thoseregions of integration, | r − r (cid:48) | ( r − r (cid:48) ) + ( x − x (cid:48) ) ≈ y − y (cid:48) , (V.25)( y − y (cid:48) ) ( x − x (cid:48) ) + ( r − r (cid:48) ) ≈ , (V.26)( y − y (cid:48) ) ( x − x (cid:48) ) [( x − x (cid:48) ) + ( r − r (cid:48) ) ] (cid:28) . (V.27) This implies that I ≈ (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − (cid:90) y − ˜ CL x dy (cid:48) y − y (cid:48) = ln (cid:18) y ˜ CL x (cid:19) (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − = (cid:26) ln (cid:18) yL x (cid:19) − ln ˜ C (cid:27) (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − . (V.28)Virtually identical reasoning can be applied to I , giv- ing the result I = (cid:26) ln (cid:18) L y − yL x (cid:19) − ln ˜ C (cid:27) (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − . (V.29)For I , a simple shift of variables of integration y (cid:48) = y + y (cid:48)(cid:48) shows that I is independent of y : I = (cid:90) ˜ CL x − ˜ CL x dy (cid:48)(cid:48) (cid:112) y (cid:48)(cid:48) + ( x − x (cid:48) ) x − x (cid:48) ) + y (cid:48)(cid:48) × (cid:34) − γ ρ v ρ γc ( v p − v ρ ) y (cid:48)(cid:48) x − x (cid:48) ) + y (cid:48)(cid:48) − (cid:18) v ρ γ ρ c γ (cid:19) ( x − x (cid:48) ) y (cid:48)(cid:48) [2( x − x (cid:48) ) + y (cid:48)(cid:48) ] (cid:35) − . (V.30)Inserting these results (V.28), (V.30), and (V.29) into our expression (V.21), and using that result in our ex-pression (V.3) for (cid:104) ( δN ) (cid:105) , we obtain (cid:104) ( δN ) (cid:105) ≈ D p ρ c m πγ (cid:40)(cid:90) L x dx (cid:90) L x dx (cid:48) (cid:90) L y dy (cid:18) ln (cid:18) yL x (cid:19) + ln (cid:18) L y − yL x (cid:19) − C (cid:19) (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − + (cid:90) L x dx (cid:90) L x dx (cid:48) (cid:90) L y dy (cid:90) ˜ CL x − ˜ CL x dy (cid:48)(cid:48) (cid:112) y (cid:48)(cid:48) + ( x − x (cid:48) ) x − x (cid:48) ) + y (cid:48)(cid:48) × (cid:34) − γ ρ v ρ γc ( v p − v ρ ) y (cid:48)(cid:48) [2( x − x (cid:48) ) + y (cid:48)(cid:48) ] − (cid:18) v ρ γ ρ c γ (cid:19) ( x − x (cid:48) ) y (cid:48)(cid:48) [2( x − x (cid:48) ) + y (cid:48)(cid:48) ] (cid:35) − (V.31)The integrals over x and x (cid:48) in the (triple) integral inthe first line of this expression trivially give a factor of L x , since the integrand is independent of x and x (cid:48) . Theremaining integral over y is elementary. The net result is4 (cid:90) L x dx (cid:90) L x dx (cid:48) (cid:90) L y dy (cid:18) ln (cid:18) yL x (cid:19) + ln (cid:18) L y − yL x (cid:19) − C (cid:19) (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − = 2 L x L y (cid:18) ln (cid:18) L y ˜ CL x (cid:19) − (cid:19) (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − . (V.32)The y integral in the remaining (quadruple) integralwhich appears on the second line in equation (V.31) canbe done immediately, since the integrand is independentof y , yielding a factor of L y . The remaining triple integralover x , x (cid:48) , and y (cid:48)(cid:48) can be done by changing variables ofintegration to new, rescaled variables u x , u y , and u (cid:48) y via x ≡ L x u x , x (cid:48) ≡ L x u (cid:48) x , y (cid:48)(cid:48) ≡ L x u y . (V.33)This gives (cid:90) L x dx (cid:90) L x dx (cid:48) (cid:90) L y dy (cid:90) ˜ CL x − ˜ CL x dy (cid:48)(cid:48) (cid:112) y (cid:48)(cid:48) + ( x − x (cid:48) ) x − x (cid:48) ) + y (cid:48)(cid:48) × (cid:34) − γ ρ v ρ γc ( v p − v ρ ) y (cid:48)(cid:48) [2( x − x (cid:48) ) + y (cid:48)(cid:48) ] − (cid:18) v ρ γ ρ c γ (cid:19) ( x − x (cid:48) ) y (cid:48)(cid:48) [2( x − x (cid:48) ) + y (cid:48)(cid:48) ] (cid:35) − = ˜ C (cid:48) L x L y , (V.34)where˜ C (cid:48) ≡ (cid:90) du x (cid:90) du (cid:48) x (cid:90) C − C du y (cid:113) ( u x − u (cid:48) x ) + u y u y + 2( u x − u (cid:48) x ) (cid:34) − γ ρ v ρ γc ( v p − v ρ ) u y u x − u (cid:48) x ) + u y − (cid:18) v ρ γ ρ c γ (cid:19) u y ( u x − u (cid:48) x ) [2( u x − u (cid:48) x ) + u y ] (cid:35) − (V.35)is an O (1) constant.Comparing (V.32) and (V.34), we see that the firstline of (V.31) actually dominates the second in the small aspect ratio limit L x (cid:28) L y that we are considering here.Therefore, we obtain, in the limit of small aspect ratio( α A (cid:28) (cid:112) (cid:104) ( δN ) (cid:105) ≈ ρ c m (cid:113) ˜ C (cid:48)(cid:48) D p √ πγ (cid:115) L x L y ln (cid:18) L y L x (cid:19) ≈ ρ c m (cid:113) ˜ C (cid:48)(cid:48) D p √ πγ L / y α A (cid:115) ln (cid:18) α A (cid:19) , (V.36)where we have defined˜ C (cid:48)(cid:48) ≡ (cid:20) − ( v p − v ρ ) γ ρ v ρ γc (cid:21) − (V.37)Note that C (cid:48) is an O (1), parameter-dependent but aspectratio and box size independent constant.This can be rewritten in terms of the mean particle number N in the same area A using (V.19), which gives∆ N = (cid:112) (cid:104) ( δN ) (cid:105) = (cid:32) ρ c mρ / (cid:33) (cid:32)(cid:115) C (cid:48) D p πγ (cid:33) N / α / A (cid:115) ln (cid:18) α A (cid:19) . (V.38)Note that, in this regime of small aspect ratio α A , ∆ N isa monotonically increasing function of α A . Recall that inthe opposite limit of α A (cid:28)
1, eqn (V.20) , we found that∆ N is a monotonically decreasing function of α A . Hence,5the maximum value of ∆ N for a given mean number ofparticles N will occur when α A ∼
1; i.e., for a roughlysquare counting box.In contrast to our above results on giant number fluc-tuations, Ref. [23] did not report any giant number fluc-tuations in their experiment. We believe this is due tothe fact that in their experiment, the density fluctuationswere probed at length scales larger than the height of thepassive fluid over the active layer, which is outside theregime of validity of our theory.
VI. BULK VELOCITY FLUCTUATIONS
We can use the relation (A20) between the bulk fluidvelocity field and the surface velocity, and our boundarycondition (II.12) for that surface velocity, to obtain ex-pressions for the bulk velocity correlations in terms of thepolarization ˆ p and density ρ correlation functions. Thisgives (cid:104) v i ( r ⊥ , z ) v j ( r ⊥ , z ) (cid:105) = v (cid:90) d q (2 π ) C pp ( q ) e − qz (cid:20) δ iy δ jy − q y zq ( δ iy q j + δ jy q i ) + q y z (cid:18) q i q j q + δ iz δ iz (cid:19) (cid:21) + v v (cid:48) a (cid:90) d q (2 π ) C pρ ( q ) e − qz (cid:20) δ iy δ jx + δ ix δ jy − z ( q y q δ ix q j + q y q δ jx q i + q x q δ jy q i + q x q δ iy q j ) + 2 q x q y z (cid:18) q i q j q + δ iz δ iz (cid:19) (cid:21) + v (cid:48) a (cid:90) d q (2 π ) C ρρ ( q ) e − qz (cid:20) δ ix δ jx − q x zq ( δ ix q j + δ jx q i ) + q x z (cid:18) q i q j q + δ iz δ iz (cid:19) (cid:21) (VI.1)The second ( δ iy q j ) term in the first integral is odd in q x if j = x . Since C pp ( q ) e − qz and q y are even in q x , the inte-gral of this term vanishes when i = x . Hence, we can re-place this term by δ iy δ jy q y . Similar arguments imply thatwe can also replace the third ( δ jy q i ) term with δ iy δ jy q y . Similarly we can replace ( δ ix q j ) with δ ix δ jx q x , and ( δ jx q i )with δ ix δ jx q x in the third integral. In the second inte-gral ( δ ix q j ) can be replaced with δ ix δ jy q y , ( δ jx q i ) with δ jx δ iy q y , ( δ jy q i ) with δ ix δ jy q x , and ( δ iy q j ) with δ iy δ jx q x respectively. Thus, we can rewrite (VI.1) as (cid:104) v i ( r ⊥ , z ) v j ( r ⊥ , z ) (cid:105) = v (cid:90) d q (2 π ) C pp ( q ) e − qz (cid:20) δ iy δ jy (cid:32) − q y zq (cid:33) + q y z (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) + v v (cid:48) a (cid:90) d q (2 π ) C pρ ( q ) e − qz (cid:20)(cid:18) δ iy δ jx + δ ix δ jy (cid:19) (1 − zq ) + 2 z q x q y (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) + v (cid:48) a (cid:90) d q (2 π ) C ρρ ( q ) e − qz (cid:20) δ ix δ jx (cid:18) − q x zq (cid:19) + q x z (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) . (VI.2)We also note that the q i q j q term in the first integral isalso odd in q x unless i = j , because, if i (cid:54) = j , then one,and only one, of the indices ( i, j ) must be x . Hence, itsintegral also vanishes if i (cid:54) = j . Identical reasoning impliesthat the third integral also vanishes for i (cid:54) = j . This makesthe first and third integrals diagonal.The second integral is non-zero only when i = x , j = y ,or i = y , j = x . This is explicit for the first term, sincethat term is proportional to δ iy δ jx + δ ix δ jy . To see that it is also true for the second term (i.e., the z q x q y term),note that term is odd in at least one of q x or q y unless i = x , j = y , or i = y , j = x , and will therefore integrateto zero otherwise.In light of these observations, we can rewrite our ex-pression (VI.2) for the velocity correlations as (cid:104) v i ( r ⊥ , z ) v j ( r ⊥ , z ) (cid:105) = v Π (1) ij ( z )+ v v (cid:48) a Π (2) ij ( z )+ v (cid:48) a Π (3) ij ( z ) , (VI.3)where we have defined Π (1) ij ( z ) ≡ (cid:90) d q (2 π ) C pp ( q ) e − qz (cid:20) δ iy δ jy (cid:32) − q y zq (cid:33) + q y z (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) (VI.4) Π (2) ij ( z ) ≡ (cid:90) d q (2 π ) C pρ ( q ) e − qz (cid:20) − zq + 2 z q x q y q (cid:21)(cid:18) δ iy δ jx + δ ix δ jy (cid:19) (VI.5)6 Π (3) ij ( z ) ≡ (cid:90) d q (2 π ) C ρρ ( q ) e − qz (cid:20) δ ix δ jx (cid:18) − q x zq (cid:19) + q x z (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) . (VI.6)Using the fact that all three of the correlation functions C ρρ ( q ), C pρ ( q ), and C pp ( q ) are proportional to q timesfunctions that depend only on the direction θ q of q , andmaking the simple change of variables q = Q z , showsthat all three of these (tensor) integrals Π (1 , , ( z )areproportional to z .Thus, we conclude that (cid:104) v i ( r ⊥ , z ) v j ( r ⊥ , z ) (cid:105) = (cid:18) z (cid:19) M , (VI.7) where we have defined the constant, parameter depen-dent matrix M ≡ v M (1) ij + v v (cid:48) a M (2) ij + v (cid:48) a M (3) ij , (VI.8)with M (1) ij ( z ) ≡ (cid:90) d Q (2 π ) f pp ( θ Q ) Q e − Q (cid:20) δ iy δ jy (cid:32) − Q y Q (cid:33) + Q y (cid:18) Q i Q j Q + δ iz δ jz (cid:19) (cid:21) , (VI.9) M (2) ij ( z ) ≡ (cid:90) d Q (2 π ) f pρ ( θ Q ) Q e − Q (cid:20) − Q + 2 Q x Q y Q (cid:21)(cid:18) δ iy δ jx + δ ix δ jy (cid:19) , (VI.10) M (3) ij ( z ) ≡ (cid:90) d Q (2 π ) f ρρ ( θ Q ) Q e − Q (cid:20) δ ix δ jx (cid:18) − Q x Q (cid:19) + Q x (cid:18) Q i Q j Q + δ iz δ jz (cid:19) (cid:21) . (VI.11)In light of the above discussion, M is a symmetric matrix,whose only non-zero off-diagonal components are M xy = M yx . Its diagonal entries are all, in general, differentfrom each other, and from M xy .The important point about these velocity correlationsis that they scale like z ; that is, inversely proportional tothe distance z from the solid surface. This z scaling willbreak down once z becomes microscopic, as can be seenas follows: our arguments above depended on our hy-drodynamic theory, which breaks down for wavevectors q comparable to an inverse microscopic length. Since theintegrals over wavevector that we have done to derive(VI.7) were dominated by q ∼ z , the calculation clearlyceases to be valid once z is a microscopic length, because then we’ll need the correlation functions at wavevectorscomparable to an inverse microscopic length, at whichour hydrodynamic theory does not apply.We now calculate the space and time-dependent veloc-ity correlators C vij ( r ⊥ − r (cid:48) ⊥ , z, z (cid:48) , t − t (cid:48) ) ≡ (cid:104) v i ( r ⊥ , z, t ) v j ( r (cid:48) ⊥ , z (cid:48) , t (cid:48) ) (cid:105) (VI.12)This can also be written in terms of the in terms ofthe polarization ˆ p and density ρ correlation functionsusing the relation (A20) between the bulk fluid velocityfield and the surface velocity, and our boundary condition(II.12) for that surface velocity. We find C vij ( r ⊥ − r (cid:48) ⊥ , z, z (cid:48) , t − t (cid:48) ) = (cid:90) dω π d q (2 π ) exp[ iω ( t − t (cid:48) ) − i q · ( r ⊥ − r (cid:48) ⊥ )] (cid:104) v i ( q , z, ω ) v i ( − q , z (cid:48) , − ω ) (cid:105) = v (cid:90) dω π d q (2 π ) exp[ iω ( t − t (cid:48) ) − i q · ( r ⊥ − r (cid:48) ⊥ ) − q ( z + z (cid:48) )] C pp ( q , ω ) (cid:20) δ iy δ jy (cid:32) − q y ( z + z (cid:48) ) q (cid:33) + q y zz (cid:48) (cid:18) q i q j q + δ iz δ jz (cid:48) (cid:19) (cid:21) + v v (cid:48) a (cid:90) dω π d q (2 π ) exp[ iω ( t − t (cid:48) ) − i q · ( r ⊥ − r (cid:48) ⊥ ) − q ( z + z (cid:48) )] C pρ ( q , ω ) (cid:20) ( δ ix δ jy + δ iy δ jx ) (cid:32) − q ( z + z (cid:48) )2 + 2 zz (cid:48) q x q y q (cid:33) (cid:21) + v (cid:48) a (cid:90) dω π d q (2 π ) exp[ iω ( t − t (cid:48) ) − i q · ( r ⊥ − r (cid:48) ⊥ ) − q ( z + z (cid:48) )] C ρρ ( q , ω ) (cid:20) δ ix δ jx (cid:18) − q x ( z + z (cid:48) ) q (cid:19) + q x zz (cid:48) (cid:18) q i q j q + δ iz δ jz (cid:48) (cid:19) (cid:21) (VI.13)where C pp ( q , ω ) , C pρ ( q , ω ), and C ρρ ( q , ω ) are given by equations (I.13) and (I.14), (I.16), (I.17), (I.18),7and(I.19), of the introduction.One of the best experimental probes of this correlationfunction is the diffusion of tracer particles. Consider firstneutrally buoyant tracer particles, which, in the absenceof diffusion, sit at a constant z , and are advected along ˆ x at a speed v by the mean motion of the fluid. If we wishto study the diffusion of these particles on a time scalesmall compared to the time required for them to diffusea distance comparable to z , we therefore only need thevelocity correlations (VI.13) at r − r (cid:48) = v ( t − t (cid:48) )ˆ x , z = z (cid:48) .In this limit, (VI.13) reduces to C vij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = C ppij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) + C pρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) + C ρρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) , (VI.14)where C ppij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = v (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C pp ( q , ω ) (cid:20) δ iy δ jy (cid:32) − q y zq (cid:33) + q y z (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) C pρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = v (cid:48) a v (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C pρ ( q , ω ) (cid:20) ( δ ix δ jy + δ iy δ jx ) (cid:32) − qz + 2 z q x q y q (cid:33) (cid:21) C ρρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = v (cid:48) a (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C ρρ ( q , ω ) (cid:20) δ ix δ jx (cid:18) − q x zq (cid:19) + q x z (cid:18) q i q j q + δ iz δ jz (cid:19) (cid:21) (VI.15)It is straightforward to see from this expression thatthe tensor C ppij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is purely diagonal.Begin by noting that the only term in the integrand in(VI.15) that has an off-diagonal component is the q i q j term. Since q is a vector in the xy plane, this termvanishes if either index i or j is z . Hence, the only off-diagonal terms are i = x , j = y , or vis-versa. In eithercase, the integrand then becomes odd in q y (recall thatboth C pp ( q , ω ) and q = | q | are even in q . These arethe only q y dependent pieces of the rest of the integrand(it appears in the 2 qz term in the argument of the firstexponential)).Hence, the integral in (VI.15), and, therefore, thecorrelation function C ppij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ), itselfvanishes if i = x , j = y , or vice-versa. Since wehave already established that the off-diagonal compo-nents of C ppij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) with one of the in-dices equal to z also vanish, this completes the proofthat C ppij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is, as claimed earlier, purely diagonal. Similarly we can show from (VI.15) that C ρρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is also a diagonal tensor.However C pρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is an off-diagonalsymmetric tensor, as is easily seen from (VI.15). Onlythe components i = x , j = y , and i = y , j = x arenon-zero, and equal.This shows that C vij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) as a wholeis symmetric tensor with four independent components.The long-time scaling behavior of the four non-zeroindependent components C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ), C vyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) and C vxy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) are essentially same, as shown below . However, C vzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) has a different behavior in thelong time limit. This is significant: diffusion along anydirection in the xy -plane is anomalous, whereas it is nor-mal (non-anomalous) along the vertical or z -direction.To see this, let us consider each of these four non-zerocomponents in turn, starting with C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ). Using (VI.14), we see that this can be expressed as C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = C ppxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) + C ρρxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) (VI.16)as C pρxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is zero, being a non-diagonal tensor itself.8We evaluate each of the two terms in (VI.16). From (VI.15), we see that C ppxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is givenby C ppxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = v (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C pp ( q , ω ) q x q y z q . (VI.17)Making the change of variables of integration from ( q , ω ) to dimensionless variables Q , Ω given by ω ≡ v Ω z , q ≡ Q z , (VI.18)and recalling that C pp ( q , ω ) obeys the scaling form (I.13),we find that C ppxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = (cid:18) v z (cid:19) (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) v ( t − t (cid:48) ) /z − Q ] F pp ( v Ω Q , θ Q ) Q x Q y Q , (VI.19)where θ Q is the angle between the rescaled, dimensionlessvector Q and the direction of mean polarization ˆ x (whichis, of course, just the same as the angle θ q between the original vector q and the direction of mean polarizationˆ x , since our rescaling (VI.18) was isotropic.Using our expression (I.14) for F pp , we have F pp (cid:18) v Ω Q , θ Q (cid:19) = 2 D p ( v Q − v ρ cos θ Q ) ( v Q − c + ( θ Q )) ( v Q − c − ( θ Q )) + ( v Q ψ (1 , θ Q ) − v ρ ψ ( ϕ, θ Q ) cos θ Q ) . (VI.20)Factoring out D p v from the numerator of this expres-sion, and v from the denominator, we see that this canbe rewritten in terms of a completely dimensionless scal- ing function of dimensionless arguments as F pp (cid:18) v Ω Q , θ Q (cid:19) = (cid:18) D p v (cid:19) H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) , (VI.21)where we have defined H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) = 2 (cid:18) Ω Q − v ρ v cos θ Q (cid:19) (cid:18) Ω Q − c + ( θ Q ) v (cid:19) (cid:18) Ω Q − c − ( θ Q ) v (cid:19) + (cid:18) Ω Q ψ (1 ,θ Q ) v − v ρ v ψ ( ϕ,θ Q ) v cos θ Q (cid:19) . (VI.22)Some of the dependence of this function on the dimen-sionless ratios v ρ v , v p v , γv , and c v is hidden in the ratios c ± ( θ Q ) v , since c ± ( θ Q ) depend on v ρ , v p , γ , and c , as dis-played in equation (IV.4).Using (VI.21) and (VI.22) in (VI.19), we find that C ppxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) itself obeys a simple scalinglaw: C ppxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = (cid:18) D p v z (cid:19) F ppxx (cid:18) v | t − t (cid:48) | z (cid:19) , (VI.23)9where the dimensionless scaling function F ppxx ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u − Q ] H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q y Q . The limiting behaviors of F ppxx ( u ) for small and large u are easy to obtain. For u (cid:28)
1, which corresponds totime differences obeying v | t − t (cid:48) | (cid:28) z , the e − Q factorin the integral (VI.24) kills the integrand before (thatis, at smaller Q than) the (Ω − Q x ) u term becomes im-portant. Hence, the integral, and F ppxx ( u ) itself, become independent of u in this limit; that is, F ppxx ( u ) = A x for u (cid:28) , (VI.24)where the constant A ppxx is given by A ppxx = (cid:90) d Ω2 π d Q (2 π ) e − Q H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q y Q . and is a function of all of the ratios v ρ v , v p v , γv , and c v . A ppxx will be of O (1) when all of these ratios are of O (1).In the opposite limit u (cid:29)
1, which corresponds to timedifferences obeying v | t − t (cid:48) | (cid:29) z , the e i (Ω − Q x ) u factor in the integral (VI.24) kills (by oscillation) all contributionsto the integral coming from Q (cid:38) u (cid:28)
1. Hence, in thedominant region of the integral, Q (cid:28)
1, and so we candrop the 2 Q term in the argument of the exponential.Doing so gives F ppxx ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u ] H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q y Q . This integral can easily be done with one further changeof variables: Ω ≡ Ω (cid:48) u , Q ≡ Q (cid:48) u , (VI.25)which gives F ppxx ( u ) = B ppxx u for u (cid:29) , (VI.26) where the constant B ppxx is given by B ppxx = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H pp (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) x Q (cid:48) y Q (cid:48) , and, like A ppxx , is again a function of all of the ratios v ρ v , v p v , γv , and c v , and will again be of O (1) when all ofthese ratios are of O (1).In summary, the behavior of the scaling function F ppxx ( u ) is given by: F ppxx ( u ) = A ppxx , u (cid:28) B ppxx u , u (cid:29) , (VI.27)We now evaluate C ρρxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ), which from0(VI.15) is expressed as C ρρxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = v (cid:48) a (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C ρρ ( q , ω ) (cid:20) − q x zq + z q x q (cid:21) (VI.28)Making the same change of variables (VI.18) as before,we obtain an scaling form for C ρρij ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ),given by C ρρxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = (cid:18) D p ρ c v (cid:48) a v z (cid:19) F ρρxx (cid:18) v | t − t (cid:48) | z (cid:19) (VI.29) where F ρρxx (cid:16) v | t − t (cid:48) | z (cid:17) is written as F ρρxx ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u − Q ] H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19)(cid:20) Q − Q x Q + Q x Q (cid:21) . where H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) is defined as H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) = 2 (cid:18) v ρ v (cid:19) (sin θ Q ) (cid:18) Ω Q − c + ( θ Q ) v (cid:19) (cid:18) Ω Q − c − ( θ Q ) v (cid:19) + (cid:18) Ω Q ψ (1 ,θ Q ) v − v ρ v ψ ( ϕ,θ Q ) v cos θ Q (cid:19) . (VI.30)The limiting behaviors of the scaling function F ρρxx ( u )can be obtained by an almost identical analysis to that used for F ppxx ( u ), giving: F ρρxx ( u ) = A ρρxx , u (cid:28) ( B ρρxx ) u + ( B ρρxx ) u + ( B ρρxx ) u , u (cid:29) A ρρxx = (cid:90) d Ω2 π d Q (2 π ) e − Q H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19)(cid:20) Q − Q x Q + Q x Q (cid:21) . and( B ρρxx ) = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H ρρ (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) , (VI.32)( B ρρxx ) = − (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H ρρ (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) x Q (cid:48) , (VI.33)1( B ρρxx ) = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H ρρ (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) x Q (cid:48) , (VI.34)Like A ppxx (VI.25) and B ppxx (VI.27), A ρρxx and ( B ρρxx ) ,( B ρρxx ) , and ( B ρρxx ) are functions of all of the ratios v ρ v , v p v , γv , and c v , and will all be of O (1) when all of theseratios are of O (1).Taking (VI.23) and (VI.29) together, we find that C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is given by C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = (cid:18) D p v z (cid:19) F ppxx (cid:18) v | t − t (cid:48) | z (cid:19) + (cid:18) D p ρ c v (cid:48) a v z (cid:19) F ρρxx (cid:18) v | t − t (cid:48) | z (cid:19) (VI.35)with F ppxx (cid:16) v | t − t (cid:48) | z (cid:17) , and F ρρxx (cid:16) v | t − t (cid:48) | z (cid:17) given by (VI.27),and (VI.31) respectively.Because of the u piece of the scaling function F ρρxx ( u ),the integral of this correlation function over u for fixed z diverges logarithmically at large u . We will show in sec-tion (VII) below that this implies logarithmic superdif-fusion (equations (I.2) and (I.3) in the x direction.We now turn to velocity correlations in the y -direction.Using (VI.14) and (VI.15), we see that the velocitycorrelation in the y direction is given by C vyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = C ppyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) + C ρρyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) (VI.36)= v (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C pp ( q , ω ) (cid:32) − q y zq + q y z q (cid:33) + v (cid:48) a (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C ρρ ( q , ω ) z q x q y q . (VI.37)Making the change of variables (VI.18) as before, we find that C vyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) also obeys a scalinglaw: C vyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = (cid:18) D p v z (cid:19) F ppyy (cid:18) v | t − t (cid:48) | z (cid:19) + (cid:18) D p ρ c v (cid:48) a v z (cid:19) F ρρyy (cid:18) v | t − t (cid:48) | z (cid:19) , (VI.38)where now the dimensionless scaling functions are given by F ppyy ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u − Q ] H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) (cid:32) − Q y Q + Q y Q (cid:33) (cid:18) Q (cid:19) . (VI.39) F ρρyy ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u − Q ] H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q y Q . (VI.40)The limiting behaviors of these scaling functions can beobtained by an almost identical analysis to that used for F ppxx ( u ), with the result: F ppyy ( u ) = A ppyy , u (cid:28) ( B ppyy ) u + ( B ppyy ) u + ( B ppyy ) u , u (cid:29) F ρρyy ( u ) = A ρρyy , u (cid:28) B ρρyy u , u (cid:29) A ppyy , A ρρyy , ( B ppyy ) , ( B ppyy ) , ( B ppyy ) ,and B ρρyy are given by A ppyy = (cid:90) d Ω2 π d Q (2 π ) e − Q H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) (cid:32) − Q y Q + Q y Q (cid:33) (cid:18) Q (cid:19) , (VI.43) A ρρyy = (cid:90) d Ω2 π d Q (2 π ) e − Q H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q y Q . (VI.44)and( B ppyy ) = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H pp (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) , (VI.45)( B ppyy ) = − (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H pp (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) y Q (cid:48) , (VI.46)( B ppyy ) = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H pp (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) y Q (cid:48) . (VI.47) B ρρyy = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H ρρ (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) x Q (cid:48) y Q (cid:48) . (VI.48)As we found earlier, all of these constants are again func-tions of all of the ratios v ρ v , v p v , γv , and c v , and will againbe of O (1) when all of these ratios are of O (1).However, the important point is that the dominantterm in (VI.47) at large u , which is obviously the ( B ppyy ) term, is clearly not integrable over u all the way out to u = ∞ ; its integral diverges logarithmically. We will show in the next subsection that this leads to superdiffu-sive motion in the y direction. Thus, taking this togetherwith our result for velocity fluctuations in the x direction,we see that there will be superdiffusion for both directionsparallel to the surface.Interestingly, the superdiffusion in the x -direction is not correlated with that in the y -direction. To see this,we calculate the cross-correlation function C vxy . Using(VI.14) and (VI.15), we see that this is given by C vxy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = C pρyy ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) )= v v (cid:48) a (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C pρ ( q , ω ) (cid:32) − qz + 2 q x q y z q (cid:33) . (VI.49)This vanishes due to the fact that C pρ ( q , ω ) is odd in q y ; see (IV.7) above. Thus, there is no cross-correlation between the superdiffusion in the x -direction and that inthe y -direction.3We now turn to the motion of the tracer particles in the z -direction. The relevant velocity correlation using(VI.14), can be expressed as C vzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = C ppzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) + C ρρzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) (VI.50)We evaluate each of the two terms in (VI.16). From (VI.15), we see that C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) is givenby C vzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) )= v (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C pp ( q , ω ) q y z + v (cid:48) a (cid:90) dω π d q (2 π ) exp[ i ( ω − v q x )( t − t (cid:48) ) − qz ] C ρρ ( q , ω ) q x z . (VI.51)Making the same change of variables (VI.18) that wemade when we were analyzing C vxx ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ), we find that C vzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) obeys an almostidentical scaling law: C vzz ( v ( t − t (cid:48) )ˆ x , z, z, t − t (cid:48) ) = (cid:18) D p v z (cid:19) F ppzz (cid:18) v | t − t (cid:48) | z (cid:19) + (cid:18) D p ρ c v (cid:48) a v z (cid:19) F ρρzz (cid:18) v | t − t (cid:48) | z (cid:19) (VI.52)where now the dimensionless scaling functions are given by F ppzz ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u − Q ] H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q y Q , (VI.53) F ρρzz ( u ) = (cid:90) d Ω2 π d Q (2 π ) exp[ i (Ω − Q x ) u − Q ] H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q . (VI.54)The limiting behaviors of these scaling functions canbe obtained by an almost identical analysis to that usedfor F ppxx ( u ), with almost identical results: F ppzz ( u ) = A ppzz , u (cid:28) B ppzz u , u (cid:29) F ρρzz ( u ) = A ρρzz , u (cid:28) B ρρzz u , u (cid:29) A ppzz = (cid:90) d Ω2 π d Q (2 π ) e − Q H pp (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q y Q , (VI.57) A ρρzz = (cid:90) d Ω2 π d Q (2 π ) e − Q H ρρ (cid:18) Ω Q , θ Q ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q x Q , (VI.58)4and B ppzz = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H pp (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) y Q (cid:48) , (VI.59) B ρρzz = (cid:90) d Ω (cid:48) π d Q (cid:48) (2 π ) exp[ i (Ω (cid:48) − Q (cid:48) x )] H ρρ (cid:18) Ω (cid:48) Q (cid:48) , θ Q (cid:48) ; (cid:26) v ρ v , v p v , γv , c v (cid:27)(cid:19) Q (cid:48) x Q (cid:48) , (VI.60)Note that, as we found for the xx correlations, A ppzz , A ρρzz , B ppzz and B ρρzz are functions of all of the ratios v ρ v , v p v , γv , and c v , and will be of O (1) when all of theseratios are of O (1).Note also that the scaling function for the z directionfalls off sufficiently rapidly as a function of its argument u that their integral over all u converges. We will showbelow that this implies that the motion of tracer particlesin the z direction is simply conventional diffusion. Fur-thermore, we will also show that the diffusion constantfor the z -direction is independent of the height z of theparticles above the surface. VII. ANOMALOUS DIFFUSION OF TRACERSA. Neutrally buoyant tracers
A neutrally buoyant tracer is a particle that is simplycarried around passively by the flows in the bulk fluid.That is, if the particle is at position r p ( t ) at time t , itsinstantaneous velocity ˙ r p ( t ) is˙ r p ( t ) = v ( r p ( t ) , t ) , (VII.1)where v ( r , t ) is the velocity field of the passive bulk fluid.Therefore, the instantaneous position of a tracer particle r p ( t ) that starts at r p ( t = 0) at some later time t is givenby r p ( t ) − r p (0) = (cid:90) t v ( r p ( t (cid:48) ) , t (cid:48) ) dt (cid:48) . (VII.2)Since diffusion, and even the superdiffusion than weeventually will find, is a much slower process that ballistictransport, we can take, for a neutrally buoyant particle(i.e., one which has no net speed in the z -direction), r p ( t (cid:48) ) = ( v t (cid:48) + x ( t = 0))ˆ x + y (0)ˆ y + z ˆ z (VII.3) in (VII.2), where we have taken into account the factthat the passive fluid is, on average, flowing along theˆ x direction at a speed v . We have also simply written z ( t = 0) as z . Doing so gives r p ( t ) − r p (0) = (cid:90) t v (( v t (cid:48) + x ( t = 0))ˆ x + y (0)ˆ y + z ˆ z , t (cid:48) ) dt (cid:48) . (VII.4)Autocorrelating the x , y , and z components of this equa-tion with themselves gives (cid:104) ( x ( t ) − x (0)) (cid:105) = (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) C vxx ( t (cid:48) − t (cid:48)(cid:48) ) , (VII.5) (cid:104) ( y ( t ) − y (0)) (cid:105) = (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) C vyy ( t (cid:48) − t (cid:48)(cid:48) ) , (VII.6) (cid:104) ( x ( t ) − x (0))( y ( t ) − y (0)) (cid:105) = (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) C vxy ( t (cid:48) − t (cid:48)(cid:48) ) = 0 , (VII.7) (cid:104) ( z ( t ) − z (0)) (cid:105) = (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) C vzz ( t (cid:48) − t (cid:48)(cid:48) ) , (VII.8)where the correlation functions C vxx ( t (cid:48) − t (cid:48)(cid:48) ), C vyy ( t (cid:48) − t (cid:48)(cid:48) ), C vxy ( t (cid:48) − t (cid:48)(cid:48) ),and C vzz ( t (cid:48) − t (cid:48)(cid:48) ) are precisely those we calcu-lated in equations (VI.35), (VI.38), and (VI.52) respec-tively of the preceding section. We will now use theseexpressions to obtain the diffusive motion - that is, thesmall departure from the mean ballistic motion at ve-locity v ˆ x of tracers in the ˆ x direction - for all threedirections.We begin with the x -direction. Using our scaling form(VI.35) for C vxx ( t (cid:48) − t (cid:48)(cid:48) ), we can rewrite this as (cid:104) ( x ( t ) − x (0)) (cid:105) = (cid:18) D p v z (cid:19) (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) F ppxx (cid:18) v | t (cid:48) − t (cid:48)(cid:48) | z (cid:19) + (cid:18) D p ρ c v (cid:48) a v z (cid:19) (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) F ρρxx (cid:18) v | t (cid:48) − t (cid:48)(cid:48) | z (cid:19) . (VII.9)Now, changing variables of integration in the second integral from t (cid:48)(cid:48) to δt ≡ t (cid:48)(cid:48) − t (cid:48) , we can rewrite this as5 (cid:104) ( x ( t ) − x (0)) (cid:105) = (cid:18) D p v z (cid:19) (cid:90) t dt (cid:48) (cid:90) t − t (cid:48) − t (cid:48) dδt F ppxx (cid:18) v | δt | z (cid:19) + (cid:18) D p ρ c v (cid:48) a v z (cid:19) (cid:90) t dt (cid:48) (cid:90) t − t (cid:48) − t (cid:48) dδt F ρρxx (cid:18) v | δt | z (cid:19) . (VII.10)Now recalling, as shown by equations (VI.26) and(VI.31), that F ppxx ( u ) ∼ /u , and F ρρxx ( u ) ∼ /u for u (cid:29) t (cid:29) (cid:16) zv (cid:17) ), we expect the first integralover δt on the right hand side of (VII.10) to converge asits limits go to ∞ , and the second integral to diverge.Thus the second integral clearly dominates at long timesand (VII.10) can be approximated by (cid:104) ( x ( t ) − x (0)) (cid:105) ≈ (cid:18) D p ρ c v (cid:48) a v z (cid:19) (cid:90) t dt (cid:48) (cid:90) t − t (cid:48) − t (cid:48) dδt F ρρxx (cid:18) v | δt | z (cid:19) . (VII.11) Changing variables of integration from δt to u ≡ v δtz ,we obtain for the integral over δt : (cid:90) t − t (cid:48) − t (cid:48) dδt F ρρxx (cid:18) v | δt | z (cid:19) = (cid:18) zv (cid:19) (cid:90) v t − t (cid:48) ) z − v t (cid:48) z du F ρρxx ( u ) . (VII.12)Using the fact that F ρρxx ( u ) is an even function of u , wecan rewrite this as: (cid:90) t − t (cid:48) − t (cid:48) dδt F ρρxx (cid:18) v | δt | z (cid:19) = (cid:18) zv (cid:19) (cid:90) v t (cid:48) z du F ρρxx ( u ) + (cid:90) v t − t (cid:48) ) z du F ρρxx ( u ) . (VII.13)It is convenient to break up the first integral in this ex-pression into a part coming from u < u > (cid:90) v t (cid:48) z du F ρρxx ( u ) = (cid:90) du F ρρxx ( u ) + (cid:90) v t (cid:48) z du F ρρxx ( u )(VII.14)The integral from 0 to 1 in this expression is clearly O (1),since, as we showed earlier, the integrand is. The secondintegral can be rewritten: (cid:90) v t (cid:48) z du F ρρxx ( u ) = ( B ρρxx ) (cid:90) v t (cid:48) z duu + (cid:90) v t (cid:48) z du (cid:18) F ρρxx ( u ) − ( B ρρxx ) u (cid:19) . (VII.15) The second integral on the right hand side of this equa-tion (VII.15) is also O (1), as can be seen from the factthat the integrand is O (1), and the expression in paren-theses falls off like u for large u , so the integral con-verges, even if the upper limit is taken to ∞ . The firstintegral in (VII.15) is, of course, elementary (about aselementary as they come, in fact!). We thereby obtain: (cid:90) v t (cid:48) z du F ρρxx ( u ) = ( B ρρxx ) ln (cid:18) v t (cid:48) z (cid:19) + O (1) . (VII.16)An almost identical analysis shows that (cid:90) v t − t (cid:48) ) z du F ρρxx ( u ) = ( B ρρxx ) ln (cid:18) v ( t − t (cid:48) ) z (cid:19) + O (1) . (VII.17)Inserting (VII.16) and (VII.17) into (VII.13) gives (cid:90) t − t (cid:48) − t (cid:48) dδt F ρρxx (cid:18) v | δt | z (cid:19) = (cid:18) zv (cid:19) (cid:26) ( B ρρxx ) (cid:20) ln (cid:18) v t (cid:48) z (cid:19) + ln (cid:18) v ( t − t (cid:48) ) z (cid:19)(cid:21) + O (1) (cid:27) . (VII.18)Inserting this in turn into (VII.11), we obtain: (cid:104) ( x ( t ) − x (0)) (cid:105) = D p ρ c v (cid:48) a v (cid:90) t dt (cid:48) (cid:26) ( B ρρxx ) (cid:20) ln (cid:18) v t (cid:48) z (cid:19) + ln (cid:18) v ( t − t (cid:48) ) z (cid:19)(cid:21) + O (1) (cid:27) . (VII.19)All of the integrals in this expression are elementary, yielding our final superdiffusive expression for the mean6squared displacement in the x direction: (cid:104) ( x ( t ) − x (0)) (cid:105) = 2( B ρρxx ) D p ρ c v (cid:48) a v t (cid:18) ln (cid:18) v tz (cid:19) + O (1) (cid:19) , (VII.20)which the alert reader will recognize as equation (I.2) ofthe introduction.We have derived this result assuming that the parti-cle has not moved appreciably in the z direction fromits original height z . This will be true only for t (cid:28) z D z ,where D z is the ( z -independent) diffusion constant we calculated above. For much longer times, the particlewill typically be a distance z ( t ) ∼ √ D z t above the sur-face. Since the logarithm is quite insensitive to the pre-cise position (i.e., to factors of O (1) in this estimate of z , it will suffice, to leading logarithmic order, to replace z in (VII.20) with z ( t ) ∼ √ D z t in (VII.20). Doing so weobtain (cid:104) ( x ( t ) − x (0)) (cid:105) = ( B ρρxx ) D p ρ c v (cid:48) a v t (cid:18) ln (cid:18) v tD z (cid:19) + O (1) (cid:19) , (VII.21)which that same alert reader will recognize as equation(I.3) of the introduction.For short times t (cid:28) (cid:16) zv (cid:17) (which, it should be noted,can actually get arbitrarily long as z → ∞ ), the ar-guments of F ppxx and F ρρxx are always much less than 1 throughout the region of integration over t (cid:48) and t (cid:48)(cid:48) in(VII.9). Therefore, F ppxx and F ρρxx can be replaced in thatintegral by their small u limits A ppxx and A ρρxx respectively.And because they are constants, we obtain ballistic scal-ing in this regime: (cid:104) ( x ( t ) − x (0)) (cid:105) = (cid:20) A ppxx (cid:18) D p v z (cid:19) + A ρρxx (cid:18) D p ρ c v (cid:48) a v z (cid:19)(cid:21) t , t (cid:28) (cid:18) zv (cid:19) . (VII.22)The interested reader can easily verify that this ballisticbehavior matches smoothly onto the z -independent dif-fusive behavior (VII.12) for t ∼ (cid:16) zv (cid:17) , up to logarithmic factors, as it should.Virtually identical reasoning applies to diffusion in the y direction, with the result: (cid:104) ( y ( t ) − y (0)) (cid:105) = 2( B ppyy ) D p t (cid:18) ln (cid:18) v tz (cid:19) + O (1) (cid:19) , (VII.23)and at very long times, (cid:104) ( y ( t ) − y (0)) (cid:105) = ( B ppyy ) D p t (cid:18) ln (cid:18) v tD z (cid:19) + O (1) (cid:19) . (VII.24)The above equation is easy to recognize as (I.4) of theintroduction. As for x -direction, the motion in the y -direction is also7ballistic for short times t (cid:28) (cid:16) zv (cid:17) : (cid:104) ( y ( t ) − y (0)) (cid:105) = (cid:20) A ppyy (cid:18) D p v z (cid:19) + A ρρyy (cid:18) D p ρ c v (cid:48) a v z (cid:19)(cid:21) t , t (cid:28) (cid:18) zv (cid:19) . (VII.25)Once again, the interested reader can easily verify thatthis ballistic behavior matches smoothly onto the z -independent diffusive behavior (VII.23) for t ∼ (cid:16) zv (cid:17) ,as it should.In the z -direction, we see that both F ppzz ( u ) ∼ /u ,and F ρρzz ( u ) ∼ /u for u (cid:29) t (cid:29) (cid:16) zv (cid:17) . This implies diffusive behavior at long times; i.e. (cid:104) ( z ( t ) − z (0)) (cid:105) = 2 D z t , (VII.26)with D z = D p (cid:90) ∞ du F ppzz ( u ) + D p ρ c v (cid:48) a v (cid:90) ∞ du F ρρzz ( u ) , (VII.27)where D z converges to a finite value.As for motion in the x and y direction, at short times t (cid:28) (cid:16) zv (cid:17) we obtain ballistic scaling for z : (cid:104) ( z ( t ) − z (0)) (cid:105) = (cid:20) A ppzz (cid:18) D p v z (cid:19) + A ρρzz (cid:18) D p ρ c v (cid:48) a v z (cid:19)(cid:21) t , t (cid:28) (cid:18) zv (cid:19) . (VII.28)As for x , this ballistic behavior for z matches smoothlyonto the z -independent diffusive behavior (VII.26) for t ∼ (cid:16) zv (cid:17) , as it should.We can also calculate the off diagonal correlation (cid:104) ( x ( t ) − x (0))( y ( t ) − y (0)) (cid:105) = (cid:90) t dt (cid:48) (cid:90) t dt (cid:48)(cid:48) C vxy ( t (cid:48) − t (cid:48)(cid:48) ) , (VII.29)which turns out to be zero as C vxy ( t (cid:48) − t (cid:48)(cid:48) ) consists of anintegral over q of C pρ ( q ) times an even function of q y .As we have seen from Eq. () that C pρ ( q ) is an odd in q y ,this integral turns out to be zero. So (cid:104) ( x ( t ) − x (0))( y ( t ) − y (0)) (cid:105) = 0 , (VII.30)The above equation is easy to recognize as (I.5) of theintroduction. B. Sedimenting tracers
For sedimenting particles, we first note that, since wejust showed that diffusion in the z direction is normal,and homogeneous, we need to consider diffusive motionin the x , and y direction as well as the x − y plane. We’llstart by considering sedimenting particles whose sedi-mentation speed v sed << v . Such particles will spenda time of O ( z v sed ) at a height of order z . Due to theaforementioned insensitivity of the logarithmic factor in(VII.20), (VII.23), and the (VII.30), we can therefore ac-curately estimate the total mean squared y displacement (cid:104) ( x ( z = z ) − x ( z = 0)) (cid:105) , (cid:104) ( y ( z = z ) − y ( z = 0)) (cid:105) , (cid:104) ( x ( z = z ) − x ( z = 0))( y ( z = z ) − y ( z = 0)) (cid:105) of a sedi-menting particle that starts at z = z and sinks at speed v sed all the way down to the surface by simply replacing t in (VII.20), (VII.23), and (VII.30) equations by z v sed ,which is the time it takes the sedimenting tracer to sinkto the bottom. Doing so gives (cid:104) ( x ( z = z ) − x ( z = 0)) (cid:105) = 2( B ρρxx ) D p ρ c v (cid:48) a v (cid:18) z v sed (cid:19) (cid:18) ln (cid:18) v v sed (cid:19) + O (1) (cid:19) , (VII.31) (cid:104) ( y ( z = z ) − y ( z = 0)) (cid:105) = 2 B ppyy D p (cid:18) z v sed (cid:19) (cid:18) ln (cid:18) v v sed (cid:19) + O (1) (cid:19) , (VII.32) (cid:104) ( x ( z = z ) − x ( z = 0))( y ( z = z ) − y ( z = 0)) (cid:105) = 0 (VII.33)8which are just the equations (I.6-I.8) of the introduction.For more rapidly sedimenting particles (i.e., denserones), for which v sed (cid:29) v , the time scale t ∼ z v sed ofthe sinking of the sedimenting particle is in the ballis- tic regime t (cid:28) (cid:16) zv (cid:17) , and so we need to use (VII.22)and(VII.25) with t replaced with the sedimenting time z v sed . This gives (cid:104) ( x ( z = z ) − x ( z = 0)) (cid:105) = A ppxx (cid:18) D p v z v (cid:19) + A ρρxx (cid:18) D p ρ c v (cid:48) a z v v (cid:19) , (VII.34) (cid:104) ( y ( z = z ) − y ( z = 0)) (cid:105) = A ppyy (cid:18) D p v z v (cid:19) + A ρρyy (cid:18) D p ρ c v (cid:48) a z v v (cid:19) , (VII.35)(VII.36)These are equations (I.6) and (I.7) of the introduction.We close this Section by noting that, for both neutrallybuoyant particles and sedimenting tracers, the aspect ra-tio (cid:112) ( y ( t ) − y (0)) / ( x ( t ) − x (0)) is not exactly 1, butis a model parameter dependent O (1) number, reflectingthe geometric anisotropy due to the x -direction being thepreferred direction of orientation for the polar particles.Nonetheless, the aspect ratio is independent of time t .This means an anisotropy exponent ζ ani that describesthe relative scaling between distances measured along the x - and y -directions is unity (i.e., ζ ani = 1), unlike in theToner-Tu model of flocking, for which ζ ani (cid:54) = 1 [13, 14]due to strongly relevant anharmonic effects. VIII. RENORMALIZATION GROUPARGUMENT AND IRRELEVANCE OFNON-LINEARITIES
So far, we have worked strictly with the linear the-ory. For dry active matter, it is well-known[13–15] thatnon-linear effects radically change the long-wavelengthbehavior (indeed, it is only the effects of non-linearitiesthat even make the ordered state in two dimensions pos-sible). It therefore clearly behooves us to ask whethernon-linearities have such important effects in our prob-lem.In this section, we will use a simple renormalizationgroup power counting argument to show that they donot. In fact, the linear theory presented earlier is asymp- totically exact at long distances.There are several sources of nonlinearities. For in-stance, treating the various speeds and pressures in theproblem as functions of ρ , and then expanding in powersof δρ produces these nonlinear terms. In addition thereare nonlinear terms which originate from the fixed lengthconstant on the polarization p .In order to ascertain the relevance or irrelevance of thenonlinear effects, we need to consider the lowest ordernonlinear terms in the equations for p y and ρ , and in the”active partial slip” boundary condition (II.12). The RGanalysis we are about to present will make it clear thatthe most important terms at long distances are those withthe fewest possible spatial derivatives, and the smallestnumber of fields.There are two types of additional terms. The first typearises very straightforwardly from expanding, e.g, thepressures P p,s ( ρ ) and the velocity v a ( ρ ), to higher or-der in ρ . The second come from including the p k p i pieceof the transverse projection operator T ki = δ ki − p k p i inthe full equation of motion (II.13). The latter gives riseto the following extra quadratic order in the fields p y and δρ contribution to ∂ t p y : ∂ t p y ( r ⊥ , t ) quad = − αp y v sx ( r ⊥ , t ) . (VIII.1)Since this term already has one power of p y multiplying v sx , it is obviously sufficient, to quadratic order in thefields, to use our linear solution (II.42) for v sx in (VIII.1).Since it is most convenient for our RG analysis to workin real space, we rewrite (VIII.1) in real space, where itreads δv sx ( r ⊥ , t ) = ( v (cid:48) a +(¯ ζ − σ ) ∂ x ) δρ ( r ⊥ , t ) − µηv (cid:48) a (cid:90) d r (cid:48)⊥ K ppρ ( r ⊥ − r (cid:48) ⊥ ) δρ ( r (cid:48) ⊥ , t )+ ζ ∂ y p y ( r ⊥ , t ) − µηv (cid:90) d r (cid:48)⊥ K pρ ( r ⊥ − r (cid:48) ⊥ ) p y ( r (cid:48) ⊥ , t ) , (VIII.2)where the kernels K ppρ ( r ⊥ ) and K pρ ( r ⊥ ) are the in-verse two-dimensional Fourier transforms of (cid:16) q + q x q (cid:17) and (cid:16) q x q y q (cid:17) respectively[41], and are given, at large dis-9tances [41], by K ppρ ( r ⊥ ) ≈ − x πr ⊥ (VIII.3)and K pρ ( r ⊥ ) ≈ − xy πr ⊥ (VIII.4) as we demonstrate in appendix C. Including such terms,the equations of motion take the form ∂p y ∂t = − v p ∂ x p y − γ (cid:90) d r (cid:48)⊥ K pp ( r ⊥ − r (cid:48) ⊥ ) p y ( r (cid:48) ⊥ ) − γ ρ ρ c (cid:90) d r (cid:48)⊥ K pρ ( r ⊥ − r (cid:48) ⊥ ) δρ ( r (cid:48) ⊥ ) − σ t ∂ y δρ + f y − γ NL (cid:90) d r (cid:48)⊥ p y ( r ⊥ ) K pρ ( r ⊥ − r (cid:48) ⊥ ) p y ( r (cid:48) ⊥ ) − g ppρ (cid:90) d r (cid:48)⊥ p y ( r ⊥ ) K ppρ ( r ⊥ − r (cid:48) ⊥ ) δρ ( r (cid:48) ⊥ ) − g ppρ (cid:90) d r (cid:48)⊥ p y ( r (cid:48) ⊥ ) K pp ( r ⊥ − r (cid:48) ⊥ ) δρ ( r (cid:48) ⊥ ) − σ NL p y ∂ x δρ + λ ∂ y p y + λ δρ∂ x p y + λ p y f x , (VIII.5)and ∂δρ∂t = − v ρ ∂ x δρ − v p ρ c ∂ y p y − λ ∂ y ( p y δρ ) + λ ∂ x ( δρ ) + ∇ · f ρ . (VIII.6)In equation (VIII.5), the non-linear terms appear on thelast three lines. The ”bare” values of the parameters inthese equations of motion (hereafter denoted by a super-script ”0”) are related as follows: γ NL = γ , (VIII.7) g ppρ = γ ρ /ρ c , (VIII.8) g ppρ = αv (cid:48) a µη , (VIII.9) σ NL = σ t , (VIII.10) λ = ( ν − λ ρv ) v − λ , (VIII.11) λ = (cid:18) ν + 12 (cid:19) v (cid:48) a − λ ρv ( v (cid:48) a ) , (VIII.12) λ = − , (VIII.13) λ = v ρ , (VIII.14) λ = − ρ (cid:48) e ( v (cid:48) a ) . (VIII.15)However, none of these parameters will continue to main-tain these relations to the other parameters upon renor-malization, which is why we have introduced them as in-dependent parameters in the equations of motion (VIII.5)and (VIII.6).Here, the kernel K pp ( r ⊥ ) is the inverse two-dimensionalFourier transform of (cid:16) q + q y q (cid:17) , and is given, at large dis- tances [41], by K pp ( r ⊥ ) ≈ − y πr ⊥ , (VIII.16)as we demonstrate in appendix C. Notice that we haverestored the number-conserving noise f ρ in (VIII.6).That noise is assumed to be of zero-mean and Gaussian-distributed with variance (cid:104) f ρ i ( x , t ) f ρ j (( x (cid:48) , t (cid:48) ) (cid:105) = 2 D ρ δ ij δ ( x − x (cid:48) ) δ ( t − t (cid:48) ) . (VIII.17)We will now assess the importance of the non-linearterms in these equations of motion using the dynamicalrenormalization group (DRG). Readers interested in amore complete and pedagogical discussion of the DRGare referred to [42] for the details of this general approach.This approach begins by decomposing the Fouriermodes of the fields p y ( q , ω ) and ρ ( q , ω ), and the noises f p and f ρ into a rapidly varying parts p >y ( q , ω ), ρ > ( q , ω ), f >p ( q , ω ) and f >ρ ( q , ω ), and slowly varying parts p y ( q , ω ) and0 ρ > ( q , ω ) from the equations of motion. We do this bysolving them iteratively for p >y ( q , ω ) and ρ > ( q , ω ). Thissolution is a series in the non-linearities which dependson the slow fields p 12 (VIII.36)Using this choice in the recursion relation for D ρ , andagain neglecting the non-linear terms, we see that dD ρ d(cid:96) = − D ρ , (VIII.37)which clearly shows that D ρ ( (cid:96) → ∞ ) → 0. Thus, D ρ isirrelevant at long distances and times, which justifies ourneglect of it in our earlier, linear analysis.We also find that all λ i , i = 1 → dλ i d(cid:96) = − (cid:18) (cid:19) λ i + nonlinear corrections , (VIII.38)and that the other three non-linear coefficients γ NL , g ppρ , and σ NL have the same RG eigenvalues: dγ NL d(cid:96) = − (cid:18) (cid:19) γ NL + nonlinear corrections ,dg ppρ d(cid:96) = − (cid:18) (cid:19) g ppρ + nonlinear corrections ,dg ppρ d(cid:96) = − (cid:18) (cid:19) g ppρ + nonlinear corrections ,dσ NL d(cid:96) = − (cid:18) (cid:19) σ NL + nonlinear corrections . (VIII.39)Therefore, all of these nonlinearities are also irrelevant atlong distances, at least if they are initially small.Since χ p,ρ are both < 0, any terms with more fields areless relevant. Likewise, any fields with more gradients arealso less relevant. This includes all possible other nonlin-earities. Therefore, all non-linearities are irrelevant. Thisimplies that our linear results are asymptotically exact atlong length and time scales. IX. SUMMARY AND CONCLUSIONS In this paper, we have studied the the stability andfluctuations of a large polar-ordered flock at a solid-liquidinterface, which is a natural intermediate case betweenthe two previously studied cases of wet and dry polar ac-tive fluids. Such a flock is affected by both the frictionforce from the solid substrate underneath and the longrange hydrodynamic interaction mediated by the overly-ing passive, isotropic bulk fluid. As a result, such a flockis simultaneously momentum nonconserving, but affectedby the hydrodynamic interactions of the bulk surround-ing fluid. Friction with the substrate also breaks Galileaninvariance .These features lead to novel behavior at long lengthand time scales, radically different from both dry andwet polar flocks. First of all, a flock of arbitrarily largesize with long range polar order is stable for a range ofthe model parameters. Although this prediction is qual-itatively same as that of the original Toner-Tu modelfor a dry polar-ordered flock [13, 14], there are signif-icant differences. For instance, the effective dampingin the present theory is O ( q ) in contrast to the O ( q )damping in the linearized Toner-Tu model [13, 14] in thelong wavelength limit. As a result, fluctuations are sig-nificantly smaller at a solid-fluid interface than in theToner-Tu model. Therefore, the predictions of thelinear theory are asymptotically exact in the long wave-length limit. This is contrast to the larger fluctuations inthe Toner-Tu model, in which the predictions from thelinear theory break down at long length and time scalesdue to the fluctuations. Furthermore, although both thesolid-fluid interface problem treated here and the Toner-Tu model are anisotropic, we find isotropic scaling in ourproblem, while the Toner-Tu model exhibits anisotropicscaling [13, 14].We have shown the existence of giant number fluctua-tion in our model, with the variance of the number scalingas the 3/4th power of the mean. In addition to this un-usual scaling, we also find that the number fluctuationsin a given area depend on the shape , as well as the size,of the area.In addition, we find that the bulk fluid is ”stirred” bythe active particles on the interface, giving rise to longranged fluctuations in the bulk fluid velocity. These inturn lead to anomalous diffusion of tagged particles inthe bulk fluid. Specifically, we find that the displacementvariances (cid:104) ( x ( t ) − x (0)) (cid:105) and (cid:104) ( y ( t ) − y (0)) (cid:105) in the planebecome anomalous, scaling as t ln t in the large time limit.In contrast, diffusion in the z -direction remains normal,i.e., (cid:104) ( z ( t ) − z (0)) (cid:105) scales linearly with t .Non-interacting particles sedimenting through the bulkfluid, as a result, will land on the solid substrate in a re-gion whose typical dimensions exhibit an anomalous log-arithmic dependence on the sedimenting speed, as sum-marized in equations (VII.31)-(VII.33).Finally, we have shown that non-linearities are irrel-evant to the long-distance, long-time behavior of thesesystems, in contrast to dry active matter.There are many possible extensions of the work re-ported here. One could, for example, consider replacingeither the bulk fluid or the bulk solid of our problem witha liquid crystal (e.g., nematic or smectic). One could alsoask how the presence of multiple species, instead of oneas here, might affect the macroscopic properties.It would also be interesting to study the order-disordertransition in the present model. Will a linear theory suf-fice, as we have found it does for the ordered phase?It would also be interesting to study a variant of ourmodel, in which the passive fluid layer has finite height h .This is precisely the experimental geometry of Ref. [23].In the limit of lateral length scales L ⊥ (cid:29) h , such a systemreduces to the system of polar-ordered flocks suspendedin a fluid on a substrate studied in Ref. [43]. Lastly, onemight consider another variant of our system, in whichthere is a bulk fluid of finite thickness resting on a solidsubstrate, and is covered by a fluid membrane at thetop surface containing self-propelled particles attachedto the membrane. How the order of the polar flock cou-ples with the membrane undulations, and how this de-pends on the bulk fluid thickness, is an interesting, andcompletely open, question.Acknowledgements: One of us (AB) thanks the SERB,DST (India) for partial financial support through theMATRICS scheme [file no.: MTR/2020/000406]. NSis partially supported by Netherlands Organization forScientific Research (NWO), through the Vidi grant No.2016/N/00075794. We thank S. Ramaswamy for shar-ing reference [35] with us. NS thanks Institut Curie andMPIPKS for their support through postdoctoral fellow-ships while some of this work was being done. AB thanksthe MPIPKS, Dresden for their hospitality, and their sup-port through their Visitors’ Program, while a portion ofthis work was underway. JT likewise thanks the MPIPKS for their hospitality, and their support through the Mar-tin Gutzwiller Fellowship, and the Higgs Center of theUniversity of Edinburgh for their support with a HiggsFellowship. Appendix A: Relating the 3D bulk fluid flow to thesurface velocity field The 3D bulk passive fluid satisfies the Stoke’s equation η ∇ v = ∇ Π , (A1)where ∇ is the 3D gradient operator, and Π is the hy-drostatic pressure associated with the bulk fluid in the z > ∇ · v = 0. Therefore, taking the divergenceof (A1) implies ∇ Π = 0 . (A2)Taking the Laplacian of (A1) and using (A2) then implies ∇ v α = 0 . (A3)We now look for plane wave solutions of this equation;that is, solutions of the form v ( r ⊥ , z ) = V ( z ) exp( i q · r ⊥ ) , (A4)where q is a two-dimensional in-plane Fourier wavevectorand r ⊥ ≡ ( x, y, 0) is the projection of r onto the plane ofthe surface.We likewise assume that the velocity v s ( x, y ) on thesurface (which by assumption is parallel to the surface)also takes a plane wave form: v s ( r ⊥ ) = v s ( q ) exp( i q · r ⊥ ) , (A5)The general solution of this equation (A3) of the form(A4) that vanishes as z → ∞ is v ( r ⊥ , z ) = ( w ( q ) + w ( q ) z ) exp( − qz ) exp( i q · r ⊥ ) . (A6)The boundary condition we impose on the 3D passivefluid at the fluid-solid interface is v z ( z = 0) = 0 and v ⊥ ( x, y, z = 0) = v s ( x, y ), where ⊥ denotes componentsof the velocity parallel to the fluid-solid interface, andThese boundary conditions entirely determine w inequation(A6): w ( q ) = v s ( q ) , (A7)which implicitly fixes w z = 0.To determine w , we take the curl of (A1) to eliminatethe pressure Π . This gives: ∇ × ∇ v = . (A8)The Laplacian of our velocity field (A6) is ∇ v = w ∇ (cid:20) exp( − qz ) exp( i q · r ⊥ ) (cid:21) + w ∇ (cid:20) z exp( − qz ) exp( i q · r ⊥ ) (cid:21) = − q w [exp( − qz ) exp( i q · r ⊥ ) . (A9)Inserting this into (A8) gives( − q ˆ z + i q ) × w = . (A10)Taking the dot product of this with q (which we remindthe reader lies in the plane of the surface, and, hence,perpendicular to ˆ z ) implies q · (ˆ z × w ) = 0 . (A11)Using the cyclic properties of triple products, this in turn implies ˆ z · ( q × w ) = 0 . (A12)Since ˆ z · q × ˆ z = 0, this only imposes a condition on thein-plane components w ⊥ of w : q × w ⊥ = , (A13)the most general solution of which is w ⊥ = q φ . (A14)To determine the remaining unknown quantities, whichare φ and w z , we insert our expression (A6) for v intothe incompressibility condition: ∇ · v = [ i ( q · v s ) + w z ] exp( − qz ) exp( i q · r ⊥ ) + ( i q · w − q w z ) z exp( − qz ) exp( i q · r ⊥ ) = 0 , (A15)where we have replaced w with v s everywhere it appears.To satisfy this equation, the coefficients of bothexp( − qz ) exp( i q · r ⊥ ) and z exp( − qz ) exp( i q · r ⊥ ) mustvanish. The former condition impliesw z = − i q · v s , (A16)which can be used in the latter to give i q · w = − iq ( q · v s ) . (A17)Using (A14) to rewrite this equation in terms of φ , andsolving for φ , gives φ = − q · v s q , (A18)which in turn implies via (A14) that w ⊥ = − ˆ q q · v s . (A19)Now using (A7), (A16), and (A19) in (A6), we obtain: v ( r ⊥ , z, t ) = [ v s − z ( q · v s )(ˆ q + i ˆ z )] e − qz + i q · r ⊥ . (A20) Appendix B: Derivation of the stability condition We will prove that the aforementioned conditions im-ply that Im(Υ( θ )) < θ in two steps. First, wewill show that Im(Υ( θ = π/ < 0. Then we will show that Im(Υ( θ )) cannot equal θ . SinceIm(Υ( θ )) is obviously a continuous function of θ , this im-plies that it can never become positive, since to get fromits negative value at θ = π/ θ = π/ < 0. Setting θ = π/ π/ + 2 iγ Υ( π/ − c = 0 , (B1)which is readily solved to giveΥ( π/ 2) = − iγ ± (cid:113) c − γ . (B2)If c and γ are both > 0, then the magnitude of theargument of the square root in this expression is clearlyless than γ . Hence, even if γ > c , so that the argumentof the square root is negative, making the square rootitself purely imaginary, it can only add to the imaginarypart of Υ a contribution smaller in magnitude than − iγ ;hence, it cannot make the imaginary part positive. If γ < c , things are even simpler: the square root in (B2)is real, and Im(Υ( θ = π/ − γ , which is < γ > c > γ > θ = π/ < 0. We now complete our proof thatIm(Υ( θ )) < all θ by showing that Im(Υ( θ )) (cid:54) = 0for any value of θ . We will prove this by contradiction:assume that Im(Υ( θ )) = 0 at some value of θ . Then atthat value of θ , Υ( θ )) is real. Therefore, the imaginarypart of (III.5) readsΥ γ (1 + sin θ ) − γv ρ cos θ (1 + ϕ sin θ ) = 0 . (B3)which can be solved for Υ:Υ = v ρ cos θ (cid:18) ϕ sin θ θ (cid:19) , (B4) Inserting (B4) into the real part of (III.5) then implies v ρ cos θ (cid:18) ϕ sin θ θ (cid:19) − v ρ ( v ρ + v p ) cos θ (cid:18) ϕ sin θ θ (cid:19) + v p v ρ cos θ − c sin θ = 0 . (B5)After a bit of algebra, and using ϕ = 1 − γ ρ γ , this can be rewritten (cid:20) (cid:18) v ρ γ ρ γ (cid:19) cos θ (cid:2) v p − v ρ + ( v p − ϕv ρ ) sin θ (cid:3) − c (1 + sin θ ) (cid:21) sin θ = 0 . (B6)Our original assumption that Im(Υ( θ )) = 0 can there-fore only be satisfied if (B6) has a solution for some real θ between − π and π . To say this another way, our systemwill be stable if (B6) has no solution for any real θ . Since, in that range of θ , sin θ (cid:54) = 0, we can divide (B6)by sin θ . Doing so, and reorganizing a bit further, wecan rewrite (B6) in dimensionless form as (cid:18) m γ ρ γ (cid:19) cos θ (cid:2) (cid:36) − (cid:36) − ϕ ) sin θ (cid:3) = (1 + sin θ ) , (B7)where we have defined the ”Mach number” m ≡ v ρ c , (B8)and the speed ratio (cid:36) ≡ v p v ρ . (B9)This equation clearly has no solution for real θ , imply-ing that the system is stable, if the left hand side of the equation is bounded above by 1, since the right hand sideis clearly bounded below by 1.We can derive such a bound by using cos θ ≤ θ ≤ (cid:36) − ≤ | (cid:36) − | , and ϕ = 1 − γ ρ γ to showthat the left hand side (LHS) obeysLHS < (cid:18) m γ ρ γ (cid:19) (cid:20) | (cid:36) − | + γ ρ γ (cid:21) . (B10)The right hand side of equation (B10) is a quadraticfunction of the ratio γ ρ γ . Requiring that it be less than 1leads to the bounds on that ratio: − | (cid:36) − | − (cid:114) ( (cid:36) − + 1 m < γ ρ γ < −| (cid:36) − | + (cid:114) ( (cid:36) − + 1 m . (B11)It is easy to see that this condition can always be satisfiedfor sufficiently small γ ρ ; in particular, the allowed regionalways includes γ ρ = 0.Note that we have shown that the condition (B11) is sufficient for stability; however, it is by no means neces-sary.Nonetheless, having such a sufficient condition proves that the ordered state in this system, unlike that of ”wet”systems, can be stable. More specifically, as long as c and γ are positive, there is always a window of stabilityat sufficiently small γ ρ .At θ = 0 , π , Im Υ( θ ) = 0; however, there are stabiliz-ing O ( q ) damping term that has been neglected, whichprovide damping at θ = 0 , π .Q.E.D. Appendix C: Evaluation of the integral for the equaltime correlation function In this appendix, we evaluate the integral I ≡ (cid:90) ∞−∞ ( s + ψ ( ϕ, θ q ) ds ( s − ψ + ( ϕ, θ q )) ( s − ψ − ( ϕ, θ q )) + s , (C1)in equation (IV.11), which arises in the calculation of theequal time correlation function C pp ( q ) ≡ (cid:104)| p y ( q , t ) | (cid:105) , is given by I = π , independent of all parameters, and θ q .This proves to be true if the two quantities ψ ± ( θ q ) ≡ c ± ( θ q ) − v ρ cos θ q Ξ( ϕ, θ q )Ψ(1 , θ q ) (C2)obey ψ + ( θ q ) > 0, and ψ − ( θ q ) < 0. We will now showthat ψ + ( θ q ), and ψ − ( θ q ) do indeed obey these inequal-ities whenever the system is stable. We will then showthat when these inequalities are satisfied, I = π .We begin by showing that ψ + ( θ q ) > 0, and ψ − ( θ q ) < 0. This can be shown by first rewriting (C2) using ourexpression (I.15) for the sound speeds c ± ( θ q ): ψ ± ( θ q ) = (cid:40)(cid:18) v ρ (1 − v p (cid:19) cos θ q ± (cid:114) 14 ( v ρ − v p ) cos θ q + c sin θ q (cid:41) / Ψ(1 , θ q ) . (C3)Since Ψ(1 , θ q ) = γ (1 + sin θ ) is always positive if γ > ψ + ( θ q ) > 0, and ψ − ( θ q ) < (cid:16) v ρ (1 − v p (cid:17) cos θ q in that equation. This isequivalent to the argument of that square root being big-ger than the square of (cid:16) v ρ (1 − v p (cid:17) cos θ q . This leads to the condition14 ( v ρ − v p ) cos θ q + c sin θ q > (cid:20)(cid:18) v ρ (1 − v p (cid:19) cos θ q (cid:21) (C4)as a necessary and sufficient condition for making ψ + ( θ q ) > 0, and ψ − ( θ q ) < , θ ) ≡ Ψ(Φ , θ )Ψ(1 , θ ) = 1 + Φ sin θ θ , (C5)this can be rewritten as (cid:18) v ρ γ ρ γ (cid:19) cos θ [ v p − v ρ + ( v p − ϕv ρ )] < c (1 + sin θ ) , (C6)which the alert reader will note is precisely the same asour stability condition (III.6). So, if our system is stable, ψ + ( θ q ) > 0, and ψ − ( θ q ) < I equa-tion (C1). This can clearly be rewritten as I = I + 2 ψI + ψ I , (C7) with I = (cid:90) ∞−∞ s ds ( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s , (C8) I = (cid:90) ∞−∞ sds ( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s , (C9)and I = (cid:90) ∞−∞ ds ( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s . (C10)We will now evaluate each of these using the methodof partial fractions. This begins with the identity1( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s = 12 is (cid:20) s − ψ + ( θ q ))( s − ψ − ( θ q )) − is − s − ψ + ( θ q ))( s − ψ − ( θ q )) + is (cid:21) . (C11)Using this, we can rewrite I (C8) as I = (cid:90) ∞−∞ s ds ( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s = 12 i ( I − − I ) , (C12)where we have defined I ± ≡ (cid:90) ∞−∞ sds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) ± is . (C13)We will now evaluate I ± ; in particular, we will showthat both are independent of ψ ± .Consider first I , and the poles of its integrand in thecomplex s plane. There are obviously two of these, whichwe can write as s , = r , e iφ , , (C14)with r , both real and positive, and φ , real. Inspectionof the denominator in (C13) shows that these must obey r e iφ + r e iφ = ψ + ( θ q ) + ψ − ( θ q ) − i (C15)and r r e i ( φ + φ ) = ψ + ( θ q ) ψ − ( θ q ) . (C16)Since ψ + > ψ − < ψ + ψ − < 0. Using this factin (C16), and recalling that r and r are both real andpositive, we see that φ + φ = π , (C17) from which it immediately follows that s = r e i ( π − φ ) = − r e − iφ . (C18)Note that (C18) implies that s and s lie on the sameside of the real axis. Using (C18) in (C15) gives r e iφ − r e − iφ = ψ + ( θ q ) + ψ − ( θ q ) − i . (C19)Taking the imaginary part of this expression gives( r + r ) sin φ = − , (C20)which implies that sin φ < 0, which in turn implies that s , and, hence, s , both lie in the lower half plane.Then if we consider the contour C in the complex s plane consisting of the real axis plus the infinite semicircle above the real axis connecting its ends, it follows from theabsence of poles in the upper half plane that (cid:73) C sds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is = 0 . (C21)The integral in this expression can be written as (cid:73) C sds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is = I + (cid:90) S sds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is , (C22)where the contour S is the aforementioned infinite semi-circle. Since that semicircle is infinite, ψ ± and is can beneglected in the denominator of the integrand relative to s . Thus we have (cid:90) S sds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is = (cid:90) S sdss = (cid:90) S dss = iπ . (C23)Using this in (C22) and using the result in (C21) gives I + = − iπ . (C24) It is straightforward to repeat this reasoning for I − ; inthat case, both poles lie in the upper half plane, and thesemicircle must lie in the lower half plane. We obtain I − = iπ . (C25)Using (C24) and (C25) in (C12), we obtain I = (cid:90) ∞−∞ s ds ( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s = π . (C26)as claimed earlier. Similar reasoning can be applied to I I 3. For I I = (cid:90) ∞−∞ sds ( s − ψ + ( θ q )) ( s − ψ − ( θ q )) + s = 12 i ( I − − I ) , (C27)where we have defined I ± ≡ (cid:90) ∞−∞ ds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) ± is . (C28)Proceeding as we did above for I , we find (cid:73) C ds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is = I + (cid:90) S ds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is = 0 . (C29)Again using the fact that semicircle is infinite, we can again neglect ψ ± and is in the denominator of the inte-grand relative to s . Thus we have (cid:90) S ds ( s − ψ + ( θ q ))( s − ψ − ( θ q )) + is = (cid:90) S dss = lim R →∞ (cid:90) S dss = lim R →∞ (1 − e − iπ ) /R = 0 . (C30)Thus we find I = 0 . (C31)Similar reasoning shows that I − = 0 , (C32)as well. Taking these two results C31 and C31 togetherin C27 gives I = 0 . (C33)Finally, applying this partial fraction approach to I ,we obtain I = (cid:18) i (cid:19) ( I − − I ) , where we have defined I ± ≡ (cid:90) ∞−∞ ds ( s − i(cid:15) ) [( s − ψ + ( θ q ))( s − ψ − ( θ q )) ± is ] , (C34)and we have shifted the pole of s above the real axis bya small amount (cid:15) , which we will take to zero at the endof our calculation. It is straightforward to check that choosing to move the pole below the real axis by a smallamount (cid:15) leads to exactly the same final answer for I .Proceeding with the choice of moving the pole abovethe axis, we note that, for this choice, all of the poles inthe integrand for I − lie in the upper half plane. There-fore, evaluating the integral by closing the contour in thelower half plane (which we can do with impunity, sincethe integrand vanishes like s as | s | → ∞ , which impliesthat the semi-infinite semi-circle with we close the con-tour contributes nothing to the integral), we find that I − = 0.On the other hand, for I , there is a single pole at s = i(cid:15) and the upper half plane, and two poles (as discussedearlier) in the lower half plane. Making the easy choiceof closing the contour in the upper half plane, and takingthe limit (cid:15) → 0, we easily find I = 2 πiψ + ψ − . (C35)Using this in (C34), we obtain our final result for I : I = − πψ + ψ − . (C36)Using this result C36 and our earlier results I = π and I = 0 in our expression C7 for the integral I gives I = π (cid:18) − ψ ψ + ψ − (cid:19) . (C37)Using our expressions (C3) and (IV.13) for ψ ± and ψ in this result, we obtain, after a little (!) algebra, ourfinal expression I = π (cid:18) (cid:18) γ ρ γ (cid:19) m sin θ cos θ (cid:110)(cid:0) θ (cid:1) − m [ (cid:36) − (cid:36) − ϕ ) sin θ ] cos θ (cid:111) (cid:19) . (C38) Appendix D: Evaluation of the non-local kernels Introducing an exponential UV cutoff, we can write K pp ( r ⊥ ) = (cid:90) d q (2 π ) (cid:32) q + q y q (cid:33) e i q · r ⊥ − aq = − ( ∂ x + 2 ∂ y ) φ ( r ⊥ ) (C1)where a ≡ π/ Λ, with Λ being the ultraviolet cutoff, and φ ( r ⊥ ) = (cid:90) d q (2 π ) (cid:18) q (cid:19) e i q · r ⊥ − aq = 14 π (cid:90) π dθ (cid:90) ∞ exp([ ir ⊥ cos θ − a ] q ) dq = 14 π (cid:90) π dθa − ir ⊥ cos θ = 14 π (cid:90) π ( a + ir ⊥ cos θ ) dθa + r ⊥ cos θ . (C2)The cos θ term in the numerator is odd under θ → π − θ ,so its integral vanishes. The remainder of the integral iseven under θ → π − θ , so we can replace its integral withtwice the integral from − π/ π/ 2. Doing so, we get φ ( r ⊥ ) = a π (cid:90) π − π/ dθa + r ⊥ cos θ . (C3)The integral can be done straightforwardly with the sub-stitution u = tan θ , with the result φ ( r ⊥ ) = 12 π (cid:112) a + r ⊥ . (C4)Inserting (C4) into the last equality of (C1) and takingthe derivatives gives K pp ( r ⊥ ) = (cid:18) a − y π ( a + r ⊥ ) / (cid:19) . (C5) Likewise, K ppρ ( r ⊥ ) = (cid:18) a − x π ( a + r ⊥ ) / (cid:19) . (C6)and, K pρ ( r ⊥ ) = − ∂ x ∂ y φ ( r ⊥ ) = − xy π ( a + r ⊥ ) / . (C7)The limiting behaviors of these for r (cid:29) a are given byequations (VIII.16) and (VIII.4) respectively. Note alsothat both are well-behaved as r → [1] K. Kruse, J.F. Joanny, F. J¨ulicher, J. Prost, and K. Seki-moto, Asters, vortices, and rotating spirals in active gels of polar filaments, Phys. Rev. Lett. , 078101(2004). [2] K. Kruse, J.F. Joanny, F. J¨ulicher, J. Prost, and K. Seki-moto, Generic theory of active polar gels: a paradigm forcytoskeletal dynamics, Eur. Phys. J. E, , 5, (2005).[3] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler,and R.E. Goldstein, Confinement stabilizes a bacterialsuspension into a spiral vortex, Phys. Rev. Lett., ,268102, (2013).[4] D. Saintillan, and M. J. Shelley, Instabilities, pattern for-mation, and mixing in active suspensions, Phys. Fluids , 123304, (2008).[5] Y. Hatwalne, S. Ramaswamy, M. Rao, S. Madan andA. Simha, Rheology of active-particle suspensions, Phys.Rev. Lett. , 118101 (2004).[6] S. Saha, R. Golestanian, and S. Ramaswamy, Clusters,asters, and collective oscillations in chemotactic colloids,Phys. Rev. E, , 062316, (2014).[7] B. Liebchen, D. Marenduzzo, I. Pagonabarraga, andM.E. Cates, Clustering and Pattern Formation inChemorepulsive Active Colloids, Phys. Rev. Lett. ,258301(2015).[8] V. Narayan, S. Ramaswamy, and N. Menon, Long-LivedGiant Number Fluctuations in a Swarming Granular Ne-matic, Science , 105, (2007).[9] L.J. Daniels, Y. Park, T. C. Lubensky, and D. J. Durian,Dynamics of gas-fluidized granular rods, Phys. Rev. E ,041301(2009).[10] A. Baskaran, and M.C. Marchetti, Hydrodynamics ofself-propelled hard rods, Phys. Rev. E , 011920 (2008).[11] P.G. de Gennes and J. Prost, The Physics of Liquid Crys-tals (Oxford University Press, Oxford, 1995)[12] T. Vicsek, A. Czir´ok, E. Ben-Jacob, I. Cohen, andO. Shochet, Novel type of phase transition in a system ofself-driven particles, Phys. Rev. Lett. , 1226, (1995).[13] J. Toner and Y. Tu, Long-Range Order in a Two-Dimensional Dynamical XY Model: How Birds Fly To-gether, Phys. Rev. Lett. , 4326 (1995).[14] J. Toner and Y. Tu, Flocks, herds, and schools: A quan-titative theory of flocking, Phys. Rev. E , (1998).[15] J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics andphases of flocks, Annals of Physics , 170 (2005).[16] R.A. Simha and S. Ramaswamy, Hydrodynamic fluctu-ations and instabilities in ordered suspensions of self-propelled particles, Phys. Rev. Lett. , 058101(2002).[17] H. Chat´e, F. Ginelli, G. Gregoire and F. Raynaud, Col-lective motion of self-propelled particles interacting with-out cohesion, Phys Rev E , 046113 (2008); F. Ginelli,The Physics of the Vicsek model, Eur. Phys. J. SpecialTopics , 2099 (2016).[18] J. Toner, Giant number fluctuations in dry active polarfluids: A shocking analogy with lightning rods, J. Chem.Phys. , 154120 (2019).[19] S. Ramaswamy, R.A. Simha, and J. Toner, Active ne-matics on a substrate: Giant number fluctuations andlong-time tails, EPL (Europhys. Lett.) , 196 (2003).[20] C. Wolgemuth, E. Hoiczyk, D. Kaiser, and G. Oster, Howmyxobacteria glide, Current Biology, , 369 (2002).[21] E. Lushi, and H. Wioland, and R. Goldstein, Fluid flowscreated by swimming bacteria drive self-organizationin confined suspensions, Proceedings of the NationalAcademy of Sciences , 9733 (2014).[22] V. Schaller, and A.R. Bausch, Topological defects anddensity fluctuations in collectively moving systems, Pro-ceedings of the National Academy of Sciences , 4488(2013). [23] A. Bricard,J.-B. Caussin, N. Desreumaux, O. Dauchot,and D. Bartolo, Emergence of macroscopic directed mo-tion in populations of motile colloids, Nature , 95(2013).[24] D. Geyer, A. Morin and D. Bartolo, Sounds and hydro-dynamics of polar active fluids, Nature Materials , 789(2018).[25] N. D. Mermin and H. Wagner, Absence of Ferro-magnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev.Lett. , 1133 (1966); P. C. Hohenberg, Existence ofLong-Range Order in One and Two Dimensions, Phys.Rev. , 383 (1967); N. D. Mermin, Absence of Order-ing in Certain Classical Systems, J. Math. Phys. , 1061(1967).[26] Exceptions to this result are two dimensional crystals(see, e.g., [27]) and fluctuating tethered membranes (see,e.g., [28]).[27] B. I. Halperin and D. R. Nelson, Theory of Two-Dimensional Melting, Phys. Rev. Lett. , 121 (1978).[28] Y. Kantor, M. Kardar, and D. R. Nelson, StatisticalMechanics of Tethered Surfaces, Phys. Rev. Lett. , 791(1986).[29] The only example of a finite Q as q → /ω in smectic-A liquid crystals, Phys.Rev. Lett. , 51 (1982); G. F. Mazenko, S. Ramaswamyand J. Toner, Breakdown of conventional hydrodynamicsfor smectic-A, hexatic-B, and cholesteric liquid crystals,Phys. Rev. A , 1618 (1983).[31] It is a peculiarity of our problem, and in particular ofthe fact that damping in this system scales as q , ratherthan the usual q , that a damping coefficient can have thedimensions of a speed, rather than a diffusion constant.[32] These are not actually the sound speeds one would obtainfrom the real part of the eigenfrequencies ω .[33] J. Toner, Harvard University Ph. D. Thesis, unpublished(1981).[34] In most of the literature, the ζ terms in this expressionare taken to be the divergence of a tensor: that is, insteadof those terms, one writes ∂ j ( ζp i p j ) = ζ ( ρ )(ˆ p · ∇ s ) p i + ζ ( ρ ) p i ( ∇ s · ˆ p ) + p i (ˆ p · ∇ s ) ζ ( ρ ), which is a special case(specifically, ζ = ζ = ζ ) of the more general expressionwe have used here. That our more general expression isactually allowed - i.e., that these terms need not add upto the divergence of a tensor in non-equilibrium systems-was first recognized by [35].[35] A. Maitra, P. Srivastava, M.C. Marchetti, J.S. Lintu-vuori, S. Ramaswamy, and M. Lenz, A nonequilibriumforce can stabilize 2D active nematics, Proceedings ofthe National Academy of Sciences , 6934 (2018).[36] P.C. Martin, O. Parodi, and P.S. Pershan, Unified hydro-dynamic theory for crystals, liquid crystals, and normalfluids, Phys. Rev. A , 2401 (1972).[37] E. Bertin, M. Droz, G. Gregoire, Hydrodynamic equa-tions for self-propelled particles: microscopic derivationand stability analysis, J. Phys. A: Math. Theor. ,445001 (2009).[38] S. Mishra, A. Baskaran, and M. C. Marchetti, Fluctu-ations and pattern formation in self-propelled particles,Phys. Rev. E 81, 061916 (2010).[39] S. Shankar, S. Ramaswamy, and M. C. Marchetti, Low-noise phase of a two-dimensional active nematic system, Phys. Rev. E , 012707 (2018); S. Mishra, A. Baskaran,and M. C. Marchetti, Fluctuations and pattern forma-tion in self-propelled particles, Phys. Rev. E 81, 061916(2010).[40] The only equilibrium systems we know of that shows gi-ant number fluctuations away from a critical point aresuperfluids. See, e.g., P. M. Chaikin and T. C. Luben-sky, Principles of condensed matter physics (CambridgeUniversity Press, Cambridge UK, 1995).[41] The apparent divergence of these kernels as r ⊥ → willactually be cutoff at small r ⊥ if we impose an ultravi-olet cutoff on our wavevector integrals, as we should,since our hydrodynamic theory does not apply out toarbitrarily large wavelength. Hence, the forms (VIII.16) and (VIII.4) should be considered as large distance ex-pansions of the kernel, asymptotically valid for r ⊥ (cid:29) a ,where a is a short-distance cutoff. An explicit demon-stration that these kernels do not divegre as r → , 732 (1977).[43] A. Maitra, P. Srivastava, M. Cristina Marchetti, S. Ra-maswamy, and M. Lenz, Swimmer Suspensions on Sub-strates: Anomalous Stability and Long-Range Order,Phys. Rev. Lett.124