Surprising mappings of 2D polar active fluids to 2D soap and 1D sandblasting

Active fluids and growing interfaces are two well-studied but very different non-equilibrium systems. Each exhibits non-equilibrium behavior quite different from that of their equilibrium counterparts. Here we demonstrate a surprising connection between these two: the ordered phase of incompressible polar active fluids in two spatial dimensions without momentum conservation, and growing one-dimensional interfaces (that is, the 1+1-dimensional Kardar-Parisi-Zhang equation), in fact belong to the same universality class. This universality class also includes two equilibrium systems : two-dimensional smectic liquid crystals, and a peculiar kind of constrained two-dimensional ferromagnet. We use these connections to show that two-dimensional incompressible flocks are robust against fluctuations, and exhibit universal long-ranged, anisotropic spatio-temporal correlations of those fluctuations. We also thereby determine the exact values of the anisotropy exponent $\zeta$ and the roughness exponents $\chi_{_{x,y}}$ that characterize these correlations.

Since many fluids flow much slower than the speed of sound, a great deal of the work done over the past two centuries on equilibrium fluids has focused on incompressible fluids [23,24]. In this paper, we consider 2D active incompressible fluids; more specifically, we consider them in rotation invariant, but non-Galilean-invariant situations in which momentum is not conserved (e.g., active fluids moving over an isotropic frictional substrate such as cells crawling on a substrate). Such an active system contains rich physics: it has recently been shown that their static-moving transition belongs to a new universality class [25]. Here, we focus on the long-range properties of the system in the moving phase.
We note that the incompressibility condition is not merely a theoretical contrivance; not only can it be readily simulated [26,27] but it can arise in a variety of real experimental situations, including systems with longranged repulsive interactions [28], and dense systems of active particles with strong repulsive short-ranged interactions, such as bacteria [26]. In addition, incompressibility plays an important role in the motile colloidal systems in fluid-filled microfluidic channels recently studied [29], although these systems differ in detail from those we study here in being two component (background fluid plus colloids).
In this paper, we formulate a hydrodynamic (i.e., longwavelength and long-time) theory of the ordered, moving phase of a 2D incompressible polar active fluid. We find that the equal-time velocity correlation functions of the type of incompressible polar active fluids we study here can be mapped exactly onto those of two equilibrium problems: a divergence-free 2D XY model (a peculiar FIG. 1: | Visual representation of the mappings. The flow lines of the ordered phase of a 2D incompressible polar active fluid, the magnetization lines of the ordered phase of divergence-free 2D XY magnets, dislocation-free 2D smectic layers, and the surfaces of a growing one-dimensional crystal (which can be obtained by taking equal-time-interval snap shots), undulate in exactly the same way over space; their fluctuations share exactly the same asymptotic scaling behavior at large length scales. Note that the vertical axis is time for KPZ surface growth and the y Cartesian coordinate for the other three systems.
type of ferromagnet different from ordinary ferromagnets, which are divergenceful) and a dislocation-free 2D smectic A liquid crystal [21,22,[30][31][32], as well as onto the time dependent correlation functions of the non-equilibrium 1+1-dimensional KPZ equation [8]. The mapping of the 2D smectic onto the 1 + 1-dimensional KPZ equation was discovered by Golubovic and Wang [21,22]; the other two mappings are new (although 2D ferromagnets with 2D dipolar interactions, which are similar but not identical systems, have also been mapped onto 2D smectics [30]). This series of mappings is illustrated in Fig. 1.
Our results imply in particular that incompressible polar active fluids can develop long-ranged orientational order (by developing a non-zero mean velocity v ) in two dimensions, just as found previously for compressible polar active fluids, but in complete contrast to their equilibrium counterparts (i.e., ordinary divergenceful ferromagnets) with underlying rotation invariance, which cannot so order. However, the scaling behavior of the velocity correlation functions is very different from those for compressible polar active fluids studied in Ref. [11,12]. Specifically, we find that the equal-time velocity correlation function in the ordered phase has the following limiting behaviours: where X ≡ |x − x |/ξ x and Y ≡ |y − y |/ξ y are rescaled lengths in the x and y directions, and we define the scaling ratio κ ≡ X Y 2/3 . Here the function Φ(κ 1) ≈ cκ 3 and the constant c ≈ 0.295 are both universal (i.e., system-independent), while C 0 and A are non-universal (i.e., system-dependent), positive, finite constants, and ξ x,y are non-universal lengths. Note that the fact that |v(r, t) − v(r , t)| 2 goes to a finite value in the large separation limit |r − r | → ∞ implies long-ranged orientational order.

Results
Model. We start with the hydrodynamic model for compressible polar active fluids without momentum conservation [9,11,12]: where v(r, t), and ρ(r, t) are respectively the coarse grained continuous velocity and density fields. All of the parameters λ i (i = 1 → 3), U , the "damping coefficients" µ B,T,2 , the "isotropic pressure" P (ρ, v) and the "anisotropic Pressure" P 2 (ρ, v) are, in general, functions of the density ρ and the magnitude v ≡ |v| of the local velocity. Note that we omit higher order damping terms because, as our analysis will show later, they are irrelevant. In addition, because we focus here on the ordered phase, µ T,B,2 is taken to be positive, as required for the stability of the ordered phase. The U term makes the local v have a nonzero magnitude v 0 in the ordered phase, by the simple expedient of having U > 0 for v < v 0 , U = 0 for v = v 0 , and U < 0 for v > v 0 . The f term is a random driving force. It is assumed to be Gaussian with white noise correlations: where the "noise strength" D is a constant parameter of the system, and i, j denote Cartesian components. Note that in contrast to thermal fluids (e.g., Model A in [24]), we are concerned with active systems that are not momentum conserving. As a result, the leading contribution to the noise correlations is of the form depicted in (4). We now take the incompressible limit by taking the isotropic pressure P only to be extremely sensitive to departures from the mean density ρ 0 . One could alternatively consider making U (ρ, v) and P 2 (ρ, v) extremely sensitive to changes in ρ as well. This would be appropriate for an active fluid near its "active jamming" [33] transition, since in that regime a small change in the local density can change the speed from a non-zero value for ρ < ρ jam to zero for ρ > ρ jam . We will discuss this case in a future publication.
Focusing here on the case in which only the isotropic pressure P becomes extremely sensitive to changes in the density, we see that, in this limit, in which the isotropic pressure suppresses density fluctuations extremely effec-tively, changes in the density are too small to affect U (ρ, v), λ 1,2,3 (ρ, v), µ B,T,2 (ρ, v), and P 2 (ρ, v). As a result, in the incompressible limit taken this way, all of them effectively become functions only of the speed v; their ρ-dependence drops out since ρ is essentially constant.
Another consequence of the suppression of density fluctuations by the isotropic pressure P is that the continuity equation (2) reduces to the familiar condition for incompressible flow, which can, as in simple fluid mechanics, be used to determine the isotropic pressure P . All of the above discussion taken together leads to the following equation of motion in tensor notation for an incompressible polar active fluid, ignoring irrelevant terms: dv , and the λ 2 and µ B terms vanish due to the incompressibility constraint ∇ · v = 0 on v. In writing (6), we absorb a term W (v) into the pressure We now analyze the implications of equation (6) for the ordered state.
Linear theory. In the ordered phase, the system spontaneously breaks rotational symmetry by moving on average along some spontaneously chosen direction which we callx; we call the direction orthogonal to thisŷ. In the absence of fluctuations (i.e., if we set the noise f in (6) to zero), the velocity will be the same everywhere in space and time, and have magnitude v 0 , which we remind the reader is defined by U (v 0 ) = 0. We treat fluctuations by expanding v around v 0x , defining u(r, t) as the small fluctuation in the velocity field about this mean: Plugging Eq. (7) into Eq. (6) and expanding to linear order in u, leads to a linear stochastic partial differential equation with constant coefficients. Like all such equations, this can be solved simply by spatiotemporally Fourier transforming, and solving the resultant linear algebraic equations for the Fourier transformed field u(q, ω) in terms of the Fourier transformed noise f (q, ω). We can thereby relate the two point correlation function |u y (q, ω)| 2 to the known correlations (4) of the random force f . Integrating the result over all frequencies ω, and dividing by 2π, gives the equal time, spatially Fourier transformed velocity autocorrelation |u y (q, t)| 2 . Details of this straightforward calculation are given in "Methods"; the result is where Γ(q) ≡ µq 2 T,2 are µ T,2 (v) evaluated at v = v 0 , and the second, approximate equality applies for all q → 0. This can be seen by noting that, for q y q 2 x and q → 0, q 2 y Γ(q)q 2 , while for q y < ∼ q 2 x and q → 0, Γ(q)q 2 ≈ µq 4 x . Hence, in both cases, (which together cover all possible ranges of q for q → 0), the approximation 2αq 2 y + Γ(q)q 2 ≈ 2αq 2 y + µq 4 x is valid. We can now obtain the real space transverse fluctuations where L is the lateral extent of the system in the xdirection (its extent in the y-direction is taken for the purposes of this argument to be infinite). Note that the longitudinal fluctuations u 2 x (r, t) are negligibale compared to u 2 y (r, t) . Using (8), the integral in (9) is readily seen to converge in the infra-red, and, hence, as system size L → ∞. Since the integral is finite, and proportional to the noise strength D, it is clear that, for sufficiently small D, the transverse fluctuations u 2 y (r, t) can be made small enough that long-ranged orientational order -i.e., a non-zero v(r, t) -is preserved in the presence of fluctuations; therefore, the ordered state is stable against fluctuations for sufficiently small noise strength D.
We show in the next section that this conclusion remains valid when nonlinear effects are taken into account (even though those nonlinearities change the scaling laws from those predicted by the linear theory).
Nonlinear Theory. We begin by expanding the full equation of motion (6) to higher order in u. This gives where the superscript "0" means that the v-dependent coefficients are evaluated at v = v 0 , and we define the "longitudinal mass" α ≡ − v0 .
The first line of equation (10) contains the linear terms, including the noise f ; the first three terms on the second line are the relevant non-linearities, while the fourth term proves to be irrelevant, as we'll soon show.
In writing (10), we have neglected "obviously irrelevant" terms, by which we mean terms that differ from those explicitly displayed in (10) by having more powers of the small fluctuations u, or more spatial derivatives of a given type. For more discussion of these "obviously irrelevant" terms, see "Methods". Note that only one of the non-linearities associated with the λ 1,2,3 terms, namely, λ 0 1 u y ∂ y u y actually remains at this point.
To proceed further, we must power count more carefully.
We only need to calculate one of the two fields u x,y , since they are related by the incompressibility condition ∇ · v = 0. We choose to solve for u y ; its Fourier transformed equation of motion can be obtained by Fourier transforming (10) and acting on both sides of the resultant equation with the transverse projection operator P lm (q) = δ lm − q l q m /q 2 which projects orthogonal to the spatial wavevector q. This eliminates the pressure term. Taking the l = y component of the resulting equation gives: where F q represents the Fourier component at wavevector q, i.e., F q [g(r)] ≡ d 2 r g(r)e −iq·r ; the "bare" value of the speed v 1 , before rescaling and renormalization, is v 1 = λ 0 1 v 0 , and Γ(q) is given after equation (8). We now rescale co-ordinates (x, y), time t, and the components of the real space velocity field u x,y (r, t) according to x → e x , y → e ζ y , t → e z t (12) where the scalings of u x (r, t) and u y (r, t) are related by the incompressibility condition. Note that our convention for the anisotropy exponent here is exactly the opposite of that used in references [9][10][11][12][13]; that is, we define ζ by q y ∼ q ζ x being the dominant regime of wavevector, while [9][10][11][12][13] defines this regime as q x ∼ q ζ y . Upon this rescaling, the form of Eq. (11) remains unchanged, but the various coefficients become dependent on the rescaling parameter .
Details of this simple power counting (including the slightly subtle question of how to rescale the projection operators) are given in "Methods". The results for the three parameters (damping coefficient µ, "longitudinal mass" α, and noise strength D) that control the size of the fluctuations in the linear theory are: µ → e (z−2) µ, α → e (z−2ζ+2) α, and D → e (z−2χy−ζ−1) D.
We now use the standard renormalization group logic to assess the importance of the non-linear terms in (11). This logic is to choose the rescaling exponents z, ζ, and χ y so as to keep the size of the fluctuations in the field u fixed upon rescaling. This is clearly accomplished by keeping α, µ, and D fixed. From the rescalings just found, this leads to three simple linear equations in the three unknown exponents z, ζ, and χ y ; solving these, we find the values of these exponents in the linearized theory: ζ lin = z lin = 2, χ ylin = −1. With these exponents in hand, we can now assess the importance of the nonlinear terms in (11) at long length scales, simply by looking at how their coefficients rescale. (We don't have to worry about the size of the actual non-linear terms themselves changing upon rescaling, because we have chosen the rescalings to keep them constant in the linear theory.) We find that all of the non-linearities whose coefficients are proportional to α are "relevant" (i.e., grow upon rescaling), while those associated with the last remaining non-linearity, λ 0 1 , associated with the λ terms get smaller upon rescaling: λ 0 1 → e − 2 λ 0 1 . Hence, this term will not affect the long-distance behavior, and can be dropped from the problem. This is very different from the compressible problem, in which the α non-linearities are unimportant, while the λ ones dominate; the reasons for this difference are discussed in "Methods".
Dropping the λ 0 1 term in (10), and making a Galilean transformation to a "pseudo-co-moving" co-ordinate system moving in the directionx of mean flock motion at speed v 1 ≡ λ 0 1 v 0 to eliminate the "convective term" v 1 ∂ x u m from the right hand side of (10), leaves us with our final simplified form for the equation of motion: We now show that Eq. (15) also describes an equilibrium system: the ordered phase of the 2D XY model subject to the divergence-free constraint ∇ · M = 0, where M is the magnetization. This connection enables us to use purely equilibrium statistical mechanics (in particular, the Boltzmann distribution) to determine the equal-time correlations of 2D incompressible polar active fluids.
Divergence-free 2D XY model. The 2D XY model describes a 2D ferromagnet whose magnetization field M(r) and position r both have two components. The Hamiltonian for this model can be written, ignoring ir-relevant terms, as [34] where µ is the "spin wave stiffness". In the ordered phase, the "potential" V (|M|) has a circle of global minima at a non-zero value of |M|, which we will take to be v 0 .
Expanding in small fluctuations about this minimum by writing M = (v 0 + u x )x + u yŷ , we obtain, keeping only "relevant" terms, where we define the "longitudinal mass" 2α ≡ .
We now add to this model the divergence-free constraint ∇ · M = 0, which obviously implies ∇ · u = 0. To enforce this constraint, we introduce to the Hamiltonian a Lagrange multiplier P (r): The simplest dynamical model that relaxes back to the equilibrium Boltzmann distribution e −βH (u) for the Hamiltonian H is [34,35] where f is the thermal noise whose statistics can also be described by Eq. (4) with D = k B T = 1/β. This TDGL equation is readily seen to be exactly Eq. (15) with µ 0 T = µ. Therefore, we conclude that the ordered phase of 2D incompressible polar active fluids has the same static (i.e., equal-time) scaling behaviors as the ordered phase of the 2D XY model subject to the constraint ∇ · M = 0.
This mapping between a nonequilibrium active fluid model and a "divergence-free" XY model allows us to investigate the fluctuations in our original active fluid model by studying the partition function of the equilibrium model.
To deal with the exact identity ∇·u = 0, we use a trick familiar from the study of incompressible fluid mechanics: we introduce a "streaming function"; i.e., a new scalar field h(r) such that Because this construction guarantees that the incompressibility condition ∇ · u = 0 is automatically satisfied, there is no constraint on the field h(r).
The field h(r) has a simple interpretation as the displacement of the fluid flow lines from set of parallel lines alongx that would occur in the absence of fluctuations, as illustrated in Fig. 2 for pointing out this pictorial interpretation to us.) This fact, which is explained in more detail in "Methods", is a consequence of the fact that, as in conventional 2D fluid mechanics, contours of the streaming function h(r) are flow lines.

. (We thank Pawel Romanczuk
This picture of a set of lines that "wants" to be parallel being displaced by a fluctuation h(r) looks very much like a 2D smectic liquid crystal (i.e., "soap"), for which the layers are actually 1D fluid stripes.
2D smectic and KPZ models. This resemblance between our system and a 2D smectic is not purely visual.
Indeed, making the substitution (19), the Hamiltonian (17) becomes (ignoring irrelevant terms like (∂ x ∂ y h) 2 , which is irrelevant compared to (∂ 2 x h) 2 because y-derivatives are less relevant than x-derivatives): where B = 2αv 2 0 and K = µv 2 0 . This Hamiltonian is exactly the Hamiltonian for the dislocation-free 2D smectic model with h(r) in Eq. (20) interpreted as the displacement field of the smectic layers, as also illustrated in Fig.  2.
The scaling behaviours of the dislocation-free 2D smectic model are extremely non-trivial, since the "critical dimension" d c below which a purely harmonic description of these systems breaks down is d c = 3 [36]. Fortunately, these non-trivial scaling behaviours are known, thanks to an ingenious further mapping [21,22] of this problem onto the 1+1-dimensional KPZ equation [8], which is a model for interface growth or erosion (e.g., "sandblasting"). In this mapping, which connects the equal-time correlation functions of the 2D smectic to the KPZ equation, the y-coordinate in the smectic is mapped onto time t in the KPZ equation with h(x, t) the height of the "surface" at position x and time t above some reference height. As a result, the dynamical exponent z KPZ of the 1+1-dimensional KPZ equation becomes the anisotropy exponent ζ of the 2D smectic. Since the scaling laws of the 1+1-dimensional KPZ equation are known exactly [8], those of the equal-time correlations of the 2D smectic can be obtained as well.
This gives [21,22] ζ = 3/2 and χ h = 1/2 as the exponents for the 2D smectic, where χ h gives the scaling of the smectic layer displacement field h(r) with spatial coordinate x. Given the streaming function relation (19) between h(r) and u(r), we see that the scaling exponent χ y for u y is just χ y = χ h − 1 = −1/2 and that the scaling exponent χ x for u x is just χ x = χ y + 1 − ζ = −1. Note that these exponents are different from those for compressible polar active fluids [9][10][11][12][13] where ζ = 5/3 and χ y = −1/5. (Note that our convention here (q y ∼ q ζ x ) is the inverse of that (q x ∼ q ζ y ) used in references [9][10][11][12][13].) The fact that both of the scaling exponents χ y and χ x are less than zero implies that both u y and u x fluctuations remain finite as system size L → ∞; this, in turn, implies that the system has long-ranged orientational order since |v(r, t) − v(r , t)| 2 remains finite as |r − r | → ∞. That is, the ordered state is stable against fluctuations, at least for sufficiently small noise D.
The velocity correlation function can be calculated through the connection between u and h. Using the aforementioned connection between 2D smectics and the 1+1-dimensional KPZ equation, the equal time layer displacement correlation function takes the form [21,22]: where we define the scaling variable κ ≡ X Y 2/3 , with X ≡ |x − x |/ξ x , and Y ≡ |y − y |/ξ y , and the non-universal constant B is an overall multiplicative factor; estimates of the non-universal nonlinear lengths ξ x,y are given in "Methods".
Rewriting the velocity correlation function (1) in terms of the fluctuation u using (7) gives (24) where C 0 = 2 u 2 (r, t) is finite, and the two correlation functions on the right hand side of the equality are just the derivatives of the layer displacement correlation function: To derive (25,26) we use (19) and the definition of C h (i.e., the first equality of formula (21)).
Inserting (25,26) into (24) and using the asymptotic forms (21,22) for C h , we obtain (as explained in more detail in "Methods") the asymptotic form of the velocity correlation function given by (1). We can also obtain the Fourier transformed equal time correlation functions; these are given in "Methods".

Discussion
We formulate a universal equation of motion describing the ordered phase of 2D incompressible polar active fluids. After using renormalization group analysis to identify the relevant non-linearities of this model, we perform a series of mathematical transformation which map our model to three other interesting, but seemingly unrelated, models. Specifically, we make heretofore unanticipated connections between four seemingly unrelated systems: the ordered phase of 2D incompressible polar active fluids, the ordered phase of the divergencefree 2D XY model, dislocation-free 2D smectics, and growing one-dimensional interfaces.
Through this connection, we show that 2D incompressible polar active fluids spontaneously break continuous rotational invariance (which their equilibrium counterparts (i.e., ordinary divergenceful ferromagnets) cannot do), and obtain the exact scaling behavior of the equal-time velocity correlation function of the original model. Because this mapping only involves equal-time correlations, the dynamical scaling of the original model is currently unknown. We hope to determine this scaling in further work.

Methods
Linear theory. In this section we give the details of the derivation of the linearized theory of incompressible polar active fluids. We begin with the linearized equation of motion, obtained by expanding Eq. (6) of the main text to linear order in the fluctuation u of the velocity around its mean value v 0x : where the superscript "0" means that the v-dependent coefficients are evaluated at v = v 0 , and we define the "longitudinal mass" α ≡ − v0 Our goal now is to determine the scaling of the fluctuations u of the velocity with length and time scales, and to determine the relative scaling of the two Cartesian components x and y of position with each other, and with time t. That is, in the language of hydrodynamics, we seek the "roughness exponents" χ x,y , the anisotropy exponent ζ, and the dynamical exponent z characterizing respectively the scaling of: velocity fluctuations u x,y , "transverse" (i.e., perpendicular to the direction of flock motion) position y, and time t with "longitudinal" (i.e., parallel to the direction of flock motion) position x. Knowing this scaling (in particular, χ x,y ) allows us to answer the most important question about this system: is the ordered state actually stable against fluctuations?
To obtain this scaling in the linear theory, we begin by calculating the fluctuations of u predicted by that theory. Since the two components of u are not independent, but, rather, locked to each other by the incompressibility condition ∇ · v = 0, it is only necessary to calculate one of them. We choose to focus on the y-component, which can be calculated by first spatio-temporally Fourier transforming (27), and then acting on both sides with the transverse projection operator P lm (q) = δ lm − q l q m /q 2 which projects orthogonal to the spatial wavevector q.

The component = y of the resultant equation then gives
where we define with µ ≡ µ 0 T + µ 0 2 v 2 0 .
We can eliminate u x from (28) using the incompressibility condition ∇·v = 0, which implies, in Fourier space, q x u x = −q y u y . Solving the resultant linear algebraic equation for u y (q, ω) in terms of f m (q, ω) gives where we define the direction-dependent "sound speed" Using Eq. (28), we can obtain |u y (q, ω)| 2 from the known correlations of the random force f (i.e., formula (4) in the main text). Integrating the result over all frequencies ω, and dividing by 2π, gives the equal time, spatially Fourier transformed velocity autocorrelation: where the second, approximate equality applies for all q → 0. This can be seen by noting that, for q y q 2 x and q → 0, q 2 y Γ(q)q 2 , while for q y < ∼ q 2 x and q → 0, Γ(q)q 2 ≈ µq 4 x . Hence, in both cases, (which together cover all possible ranges of q for q → 0), the approximation 2αq 2 y + Γ(q)q 2 ≈ 2αq 2 y + µq 4 x is valid. Equation (32) implies that fluctuations diverge most rapidly as q → 0 if q is taken to zero along a locus in the q plane that obeys q y < ∼ q 2 x ; along such a locus, asymptotically, |u y (q, t)| 2 ∝ 1 q 2 . In contrast, along all other locii, i.e., those for which q y q 2 x , |u y (q, t)| 2 ∝ q 2 x q 2 y 1 q 2 . In this sense, one can say that the regime q y < ∼ q 2 x shows the largest fluctuations at small q; this implies the anisotropy exponent ζ = 2.
We can get the dynamical exponent z predicted by the linear theory by inspection of (30), although some care is required. The form of the first term in the denominator might suggest ω ∝ q, which would imply z = 1. However, the propagating c(q)q term in this expression does not appear in our final expression (32) for the fluctuations; rather, these are controlled entirely by the damping term Balancing ω against that term in the dominant regime of wavevector q y ∼ q 2 x gives ω ∝ q 2 x , which implies z = 2. Now we seek χ y , which determines whether or not the ordered state is stable against fluctuations in an arbitrarily large system. This can be obtained by looking at the real space fluctuations u 2 where L is the lateral extent of the system in the x-direction (its extent in the y-direction is taken for the purposes of this argument to be infinite). Using (32), this integral is readily seen to converge in the infra-red, and, hence, as system size L → ∞. Since the integral is finite, and proportional to the noise strength D, it is clear that, for sufficiently small D, the transverse u y fluctuations in real space can be made small enough that long-ranged orientational order, and, hence, a nonzero v(r, t) , is preserved in the presence of fluctuations; the ordered state is stable against fluctuations for sufficiently small noise strength D.
The exponent χ y can be obtained by looking at the departure δu 2 y of the u y fluctuations from their infinite system limit: (2π) 2 |u y (q, t)| 2 ; we define the "roughness exponent" χ y by the way this quantity scales with system size L: δu 2 y ∝ L 2χy . Note that this definition of χ y requires χ y < 0, since it depends on the existence of an ordered state, which necessarily implies that the velocity fluctuations δu 2 y do not diverge as L → ∞. If u 2 y (r, t) | L=∞ is not finite, one can obtain χ y by performing exactly the type of scaling argument outlined here directly on u 2 y (r, t) | L itself. Approximating (32) for the dominant regime of wavevector q y ∼ q 2 x , and changing variables in the integral from q x,y to Q x,y according to q x ≡ Qx L , q y ≡ Qy L 2 shows that δu 2 y ∝ L −1 , and hence χ y = − 1 2 . Note also that the fluctuations of u x are much smaller than those of u y . This can be seen by using the incompressibility condition, which implies, in Fourier space, which is clearly finite as q → 0 along any locus; indeed, it is bounded above by D 2α . We can calculate a roughness exponent χ x for u x for the linear theory from this result exactly as we calculate the roughness exponent χ y for u y ; we find χ x = 1 − ζ + χ y = − 3 2 . We shall see in the next section that the first line of this equality also holds in the full non-linear theory, even though the values of the exponents χ x , ζ, and χ y all change.
The fact that u x has much smaller fluctuations than u y means that we have to work to higher order in u y than in u x when we treat the non-linear theory, as we do in next section.
Mapping to an equilibrium "incompressible" magnet. We now go beyond the linear theory, and expand the full equation of motion (6) of the main text to higher order in u. We obtain We keep terms that might naively appear to be higher order in the small fluctuations (e.g., the u 3 y δ my term relative to the u x u y δ my term) because, as we saw in the linearized theory, the two different components u x,y of u scale differently at long length scales. Hence, it is not immediately obvious, e.g., which of the two terms just mentioned is actually most important at long distances. We therefore, for now, keep them both. For essentially the same reason, it is not obvious whether u 2 y δ mx or u 3 y δ my is more important, so we shall for now keep both of these terms as well.
On the other hand, it is immediately obvious that a term like, e.g., u x u 2 y δ mx is less relevant than u 2 y δ mx , since, whatever the relative scaling of u x and u y , u x u 2 y δ mx is much smaller at large distances than u 2 y δ mx , since u x is. Likewise, we drop the term 1 2 dλ1 dv v=v0 u 2 y ∂ x u y δ my , since it is manifestly smaller, by one ∂ x , than the u 3 y δ my term already displayed explicitly in (34).
This sort of reasoning guides us very quickly to the reduced model (34). As explained in the main text, acting on both sides of (34) with the transverse projection operator P lm (q) = δ lm − q l q m /q 2 which projects orthogonal to the spatial wavevector q eliminates the pressure term. Then taking the l = y component of the resulting equation gives (11) of the main text, which we now use to calculate the rescaled coefficients.
To do this, we must also determine how the projection operators P yx and P yy rescale upon the rescalings (i.e., (12) of the main text). Since in the linear theory (see, e.g., the u y -u y correlation function (32)) fluctuations are dominated by the regime q y < ∼ q 2 x , it follows that P yx (q) = − qxqy q 2 ≈ −q y /q x 1 and P yy (q) = 1 − q 2 y q 2 ≈ 1. This implies that these rescale according to P yx (q) → e (1−ζ) P yx (q) , P yy (q) → P yy (q) . (35) Performing the rescalings (12)(13)(14) of the main text, and (35) above on the equation of motion (11) of the main text, we obtain, from the rescalings of first three (i.e., the linear) terms on the right hand side the following rescalings of the parameters: and α → e (z−2ζ+2) α .
Note that the Γ(q) term in (11) of the main text involves two parameters (µ and µ 0 T ); hence, we get the rescalings of both of these parameters from this term.
Similarly, looking at the rescaling of the non-linear terms proportional to u 2 y and u 3 y , respectively, we obtain the rescalings: We recover the first of these by looking at the rescaling of the non-linear term proportional to u x u y as well. We note that the two rescalings (38) are both consistent with (37) By power counting on the u y ∂ y u y term, we obtain the rescaling of λ 0 1 : Finally, by looking at the rescaling of the noise correlations (i.e., (4) of the main text), we obtain the scaling of the noise strength D: We now use the standard renormalization group logic to assess the importance of the non-linear terms in (11) of the main text. This logic is to choose the rescaling exponents z, ζ, and χ y so as to keep the size of the fluctuations in the field u fixed upon rescaling. Since, as we saw in our treatment of the linearized theory (in particular, Eq. (32)), that size is controlled by three parameters: the "longitudinal mass" α, the damping coefficient µ, and the noise strength D, the choice of z, ζ, and χ y that keeps these fixed will clearly accomplish this. From the rescalings (36), (37), and (41), this leads to three simple linear equations in the three unknown exponents z, ζ, and χ y ; solving these, we find the values of these exponents in the linearized theory: which, unsurprisingly, are the linearized exponents we found earlier.
With these exponents in hand, we can now assess the importance of the non-linear terms in (11) of the main text at long length scales, simply by looking at how their coefficients rescale. (We don't have to worry about the size of the actual non-linear terms themselves changing upon rescaling, because we have chosen the rescalings to keep them constant in the linear theory.) The mass α, of course, is kept fixed. Inserting the linearized exponents (42) into the rescaling relation (39) for v 0 , we see that Since v 0 appears in the denominator of all three of the non-linear terms associated with α, and α itself is fixed, this implies that all three of those terms are "relevant", in the renormalization group sense of growing larger as we go to longer wavelengths (i.e., as grows). As usual in the RG, this implies that these terms ultimately alter the scaling behavior of the system at sufficiently long distances. In particular, the exponents z, ζ, and χ x,y change from their values (42) predicted by the linear theory. The same is not true of the λ 0 1 non-linearity, however, because it is irrelevant; that is, it gets smaller upon renormalization. This follows from inserting the linearized exponents (42) into the rescaling relation (40) for λ 0 1 , which gives which shows clearly that λ 0 1 vanishes as → ∞; that is, in the long-wavelength limit.
Since λ 0 1 was the only remaining non-linearity associated with the λ terms in our original equation of motion (34), we can accurately treat the full, long distance behavior of this problem by leaving out all of those nonlinear terms.
Doing so reduces the equation of motion (34) to Before proceeding to analyze this equation, we note the differences between the structure of this problem and that for the compressible case. In the compressible problem, there is no constraint analogous to the incompressibility condition relating u x and u y . Hence, u x is free to relax quickly (to be precise, on a time scale 1 2α ) to its local "optimal" value, which is readily seen to be − u 2 y 2v0 . Once this relaxation has occurred, all of the non-linearities associated with α drop out of that compressible problem, leaving the λ non-linearities as the dominant ones. For a detailed discussion of the rather tricky analysis of the compressible problem that leads to this conclusion, see Ref. [41]. Here, in the incompressible problem, u x is, because of the incompressibility constraint, not free to relax in such a way as to cancel out the α nonlinearities, which, because they involve no spatial derivatives, wind up dominating the λ non-linearities, which do involve spatial derivatives. In addition, the suppression of fluctuations by the incompressibility condition, which as we've already seen in the linear theory, makes the λ non-linearities not only less relevant than the α ones, but actually irrelevant. Hence, we can drop them in this incompressible problem, leaving us with Eq. (45) as our equation of motion.
As one final simplification, we make a Galilean transformation to a "pseudo-co-moving" co-ordinate system moving in the directionx of mean flock motion at speed λ 0 1 v 0 . Note that if the parameter λ 0 1 had been equal to 1, this would be precisely the frame co-moving with the flock. The fact that it is not is a consequence of the lack of Galilean invariance in our problem.
This boost eliminates the "convective" term λ 0 1 v 0 ∂ x u m from the right hand side of (45), leaving us with our final simplified form for the equation of motion: which is just equation (15) of the main text.
Mapping of equilibrium "incompressible" magnet to 2D smectic. We begin by demonstrating the pictorial interpretation of the "streaming function" introduced in the main text via This implies that the streaming function φ for the full velocity field v(r) = v 0x + u, defined via v x = ∂ y φ, v y = −∂ x φ, is given by As in conventional 2d fluid mechanics, contours of the streaming function φ are flow lines. When the system is in its uniform steady state (i.e., v = v 0x ), these contour lines, defined via φ = nC, n = 0, 1, 2, 3...
where C is some arbitrary constant, are a set of parallel, uniformly spaced lines given by y n = nC/v 0 . Now let's ask what the flow lines are if there are fluctuations in the velocity field: v = v 0x + u. Combining our expression for φ (48) and the expression (49) for the flow lines, we see that the positions of the flow lines are now given by y n = nC/v 0 + h, n = 0, 1, 2, 3... , which shows that h(r) can be interpreted as the local displacement of the flow lines from their positions in the ground state configuration. This picture of a set of lines that "wants" to be parallel being displaced by a fluctuation h(r) looks very much like a 2D smectic liquid crystal (i.e., "soap"), for which the layers are actually one-dimensional fluid stripes.
This resemblance between our system and a 2D smectic is not purely visual. Indeed, making the substitution (47), the Hamiltonian (17) of the main text becomes (ignoring irrelevant terms like (∂ x ∂ y h) 2 , which is irrelevant compared to (∂ 2 x h) 2 because y-derivatives are less relevant than x-derivatives) where B = 2αv 2 0 and K = µv 2 0 . This Hamiltonian is exactly the Hamiltonian for the equilibrium 2D smectic model with h in Eq. (51) interpreted as the displacement field of the smectic layers. For the equilibrium 2D smectic the partition function is where it should be noted that there is no constraint on the functional integral over h(r) in this expression, since, as noted earlier, h(r) is unconstrained.
Since the variable transformation Eq. (47) is linear, the partition functions for the smectic: Z s (Eq. (52)) and that for the constrained XY model: are the same up to a constant Jacobian factor, which changes none of the statistics.
To summarize what we have learned so far: we have successfully mapped the model for the ordered phase of an incompressible polar active fluid onto the ordered phase of the equilibrium 2D XY model with the constraint ∇ · u = 0, which in turn we have mapped onto the standard equilibrium 2D smectic model [30] The scaling behaviors of the former can therefore be obtained by studying the latter. Note that the connection between our problem and the dipolar magnet, which was studied in [30], is that the long-ranged dipolar interaction in magnetic systems couples to, and therefore suppresses, the longitudinal component of the magnetization. See [30] for more details.
Mapping the 2D smectic to the one-dimensional KPZ equation. Fortunately, the scaling behaviours of the equilibrium 2D smectic model are known, thanks to an ingenious further mapping [21,22] of this problem onto the 1+1-dimensional KPZ equation [8], which is a model for interface growth or erosion (e.g., "sandblasting"). In this mapping, which connects the equaltime correlation functions of the 2D smectic to the 1+1-dimensional KPZ equation, the y-coordinate in the smectic is mapped onto time t in the 1+1-dimensional KPZ equation with h(x, t) the height of the "surface" at position x and time t above some reference height. As a result, the dynamical exponent z KPZ of the 1+1dimensional KPZ equation becomes the anisotropy exponent ζ of the 2d smectic. Since the scaling laws of the 1+1-dimensional KPZ equation are known exactly [8], those of the equal-time correlations of the 2D smectic can be obtained as well.
This gives [21,22] ζ = 3/2 and χ h = 1/2 as the exponents for the 2D smectic, where χ h gives the scaling of the smectic layer displacement field h with spatial coordinate x. Given the streaming function relation (47) between h and u, we see that the scaling exponent χ y for u y is just χ y = χ h − 1 = −1/2 and that the scaling exponent χ x for u x is just χ x = χ y + 1 − ζ = −1.
The fact that both of the scaling exponents χ y and χ x are less than zero implies that both u y and u x fluctuations remain finite as system size L → ∞; this, in turn, implies that the system has long-ranged orientational order. That is, the ordered state is stable against fluctuations, at least for sufficiently small noise D.
The velocity correlation function can be calculated through the connection between u and h. Using the aforementioned connection between 2D smectics and the 1 + 1-dimensional KPZ equation, the layer displacement correlation function takes the form [21,22]: where the scaling variable κ ≡ X Y 2/3 , X = |x − x |/ξ x , and Y = |y − y |/ξ y , B is a non-universal overall multiplicative factor extracted from the scaling function, and the non-universal nonlinear lengths ξ x,y are calculated in the next section. The universal scaling function Ψ has been numerically estimated (www-m5.ma.tum.de/KPZ) [37][38][39][40]: where for κ 1, Here, the constants c and c 1,2 are all universal and are given by c ≈ 0.295, c 1 ≈ 1.843465, and c 2 ≈ 1.060... (www-m5.ma.tum.de/KPZ) [39,40]. Only the lengths ξ x,y , and the overall multiplicative factor of B in (54) are non-universal (i.e., system-dependent).
Rewriting the velocity correlation function (Eq. (1) of the main text) in terms of the fluctuations u of the velocity from its mean value (as defined in Eq. (7) of the main text), we find |v(r, t) − v(r , t)| 2 = C 0 − 2 u y (r, t)u y (r , t) − 2 u x (r, t)u x (r , t) , (57) where C 0 = 2 |u(r, t)| 2 is finite. The two correlation functions on the right hand side of the equality are just the derivatives of the layer displacement correlation function: The velocity component scaling functions Ψ x,y can both be expressed in terms of the height scaling function Ψ, via Ψ y (κ) = κ (2Ψ + κΨ ) , Ψ x (κ) = κ 4 (5Ψ + 2κΨ ) .
Using the asymptotic forms (55) for the height scaling function Ψ in (60) and (61), we obtain the asymptotic behaviors: Using these expressions (62, 63) for the scaling functions in the scaling expressions (58, 59) for the u x and u y correlation functions, and in turn using those in our expression (57) for the velocity correlation function, we obtain x−x y−y 2 , κ 1 where the non-universal constant A is given by We've also defined a new universal function which has the same limiting behaviour as Φ h , namely, since the additive logarithm in (67) is sub-dominant to the leading κ 3 term.
Calculation of the nonlinear lengths. The nonlinear lengths ξ x,y can be calculated most conveniently from the equilibrium 2D smectic model (51). By definition, ξ x and ξ y are the lengths along x and y beyond which the anharmonic terms in Eq. (51) become important. To determine these lengths, we treat the anharmonic terms perturbatively, and calculate the lowest order correction to the harmonic terms. In a finite system of linear dimensions L x,y , this "perturbative" correction will indeed be perturbative (i.e., small) compared to the "bare" values of the harmonic terms. However, they grow without bound with increasing L x,y , and, hence, eventually cease to be small; that is, the perturbation theory breaks down at large L x,y . The values of L x,y above which the perturbation theory breaks down are the nonlinear lengths ξ x,y . Calculating the lowest order correction to the compression modulus B (i.e., the coefficient of (∂ y h) 2 in the smectic Hamiltonian (51)) can be graphically represented by the Feynman diagram in Fig. 3. This leads to a correction to compression modulus: In this calculation, we have taken L y , the system size along y, to be infinite. By the definition of ξ x , |δB| = B for L x = ξ x , which gives where in the second equality we have used the relations B = 2αv 2 0 , K = µv 2 0 , and k B T = D between the parameters of the smectic and those of the original incompressible active fluid.
Likewise, doing the same calculation for L y = ξ y , L x = ∞, we find