Hydrodynamic theory of flocking at a solid-liquid interface: long range order and giant number fluctuations
Niladri<EMAIL_ADDRESS>Instituut-Lorentz, Leiden University, P.O. Box 9506, 2300 RA Leiden, The Netherlands
Abhik<EMAIL_ADDRESS><EMAIL_ADDRESS>Condensed Matter Physics Division, Saha Institute of
Nuclear Physics, Calcutta 700064, West Bengal, India

John<EMAIL_ADDRESS>Department of Physics and Institute of Theoretical Science, University of Oregon, Eugene, Oregon 97403, USA

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 surface are 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 large area. Specifically, these number fluctuations grow as the $3/4$th power of the mean  number within the area.  Stable systems also exhibit anomalously rapid diffusion of tagged particles suspended in the passive fluid along any directions in a plane parallel to the solid-liquid interface,  whereas the diffusivity along the direction perpendicular to the plane is non-anomalous.
In other parameter regimes, the system becomes unstable.

§ INTRODUCTION

Active fluids  are fluids containing  self-propelled  units  (“active particles"),
each capable
of converting stored or ambient free energy into motion.    These intrinsically non-equilibrium systems often display orientationally ordered phases,  both in
living [1, 2, 3, 4, 5] and
nonliving [6, 7, 8, 9, 10]  systems. Such phases of active fluids exhibit many phenomena impossible in equilibrium orientationally ordered phases (e.g., nematics [11]), among them spontaneous breaking of continuous symmetries in two dimensions [12, 13, 14, 15], instability in the extreme Stokesian limit [16],
and giant number fluctuations [17, 18, 19].

As in  equilibrium systems, the  “hydrodynamic" (i.e., long distance, long-time) behavior of active fluids depends crucially on whether or not momentum is conserved.
In what has become the standard nomenclature of this field, systems lacking momentum conservation due to, e.g.,  friction with a substrate, are referred to as “dry" [22, 14, 19, 20, 21], while active fluids with momentum conservation are called “wet" [23, 23, 24, 25, 26].

In this paper, we treat a heretofore unconsidered  class of systems that is a natural hybrid of these two cases: a collection of
polar active particles at a solid-liquid interface. Our study is inspired by experiments on high-density actin motility
assays [27],
in which highly concentrated actin filaments on a solid-fluid interface are
propelled by  motor proteins. Another experimental realization is the work of Bricard et al  [28, 29], who studied the emergence of macroscopically directed motion in “Quincke rotators": a dilute collection of motile colloids, in contact with a solid substrate, and immersed in a passive bulk electrolytic  fluid. The colloids spontaneously roll, making them self propelled, when a sufficiently strong electric field is applied to the electrolytic fluid.

The hybrid system we consider here is in some ways reminiscent of  the phenomenon of “wet to dry crossover"  studied by other authors [30, 31]. However, there are important differences. The systems studied by [30, 31] are fully two dimensional, but had both viscous and frictional damping. In such a problem, the asymptotic, longest wavelength behavior is that of dry active matter. In our case, as long as our bulk fluid is semi-infinite, the asymptotic long-distance behavior is completely different, even in its scaling, from dry active matter. The mixture of wet and dry behavior in our system is not a crossover phenomenon, only occurring over an intermediate range of length scales. Rather, it is the true, asymptotic,  long-wavelength (i.e., hydrodynamic) behavior (again, for a semi-infinite bulk passive fluid).

The systems that we consider here differ  from both dry and  wet active matter, as defined above, by having both
friction from the underlying solid substrate and the long range
hydrodynamic interactions due to the overlying bulk passive fluid. Very
interestingly, the  experiments of  [27] found giant fluctuations of the active particle
number in a fixed volume,  with a standard deviation that scaled with $N^{0.8}$ in the
ordered phase, where $N$ is the mean number in the volume. Although the exponent of $0.8$ is the value predicted [17, 18] for dry polar active fluids, we will argue here that this is a coincidence, and that in fact polar fluids at a solid-liquid interface belong to an entirely different and heretofore unstudied universality class, with number fluctuations scaling like $N^{0.75}$, which probably lies within experimental error of the $0.8$ exponent found in the experiments.

In this  article, inspired by the experiments of Ref. [27, 28], we formulate
the hydrodynamic theory for a collection of polar, self propelled particles
at the interface between a solid
substrate and a bulk isotropic passive liquid that is semi-infinite in the
$z$-direction, as illustrated in Fig. <ref>.
We focus on the extreme Stokesian limit, in which the momentum of both the active particles and the bulk fluid is entirely negligible relative to viscous drag. The theory is then used to study the stability and fluctuations of a uniformly ordered active polar state of the particles.

(Color online) Schematic diagram of a layer of active polar particles
embedded on a solid substrate surrounded by passive ambient fluid from the top

The governing equations of the active particle number density and polarization fluctuations in this system, which are introduced here for the first time are very different from those of “dry" and “wet" active matter”.The breakdown of the Galilean invariance due to the friction from the substrate and the long ranged hydrodynamic interactions mediated by the passive bulk fluid lead to novel linear couplings between the density and polarization fluctuations, and  non-local damping. The resultant behavior constitutes  a new universality class, whose properties we describe in detail.

Our most surprising result is that this system can exhibit a long range ordered polar state, even in the presence of noise. This should be contrasted with both “wet" active systems, which are generically unstable[16] at very low Reynolds number, and with equilibrium systems, which cannot develop long ranged orientational order in two dimensions [32, 33, 34, 35].  Further, the number fluctuations in these systems are giant;  i.e.,
standard deviation $\sqrt{\langle(N-\langle N\rangle)^2\rangle}$ of the number $N$ of the active particles  contained in a fixed open area grows  with its average $\langle N\rangle$ much faster than its equilibrium off-critical counterpart. Specifically, we find
√(⟨(N-⟨N⟩)^2⟩)∝⟨N⟩^3/4  .

We also find that the diffusivities of tagged particle in any direction parallel to the active fluid layer is logarithmically  anomalous,  by which we mean the mean squared displacement $\langle|\brp(t)|^2\rangle$ grows according to  $<|r_\parallel(t)|^2>\propto t\ln (t)$, in contrast to the usual linear in $t$ scaling of normal diffusion. Diffusion normal to the plane remains conventional, with  $\langle|r_\perp(t)|^2\rangle\propto t$.

Unlike dry polar active matter, which can also exhibit long range order and giant number fluctuations [12, 13, 14, 15], these systems do so even in a linear theory. Indeed, such a linear theory proves to provide an asymptotically exact long wavelength description of these systems, again in contrast to dry active matter.

The principal condition for stability of our systems  is similar to the stability criterion for a simple compressible bulk fluid: the analog in our system of the bulk compressibility must be positive, as must our analog of the viscosities.

Thus, in contrast to “wet" active matter in the “Stokesian" limit[16, 15], our “mixed" system can be generically stable. Indeed, the requirements for stability are as easily met for these systems as for an equilibrium fluid.

These results are testable in standard experiments, e.g., in systems in which the active particles can be imaged, like those of [27, 28].

The remainder of this paper is organized as follows. We begin with a short summary of the technical results, which include the governing hydrodynamic equations themselves, in section <ref>.
In section (<ref>), we formulate the linear hydrodynamic theory, including the equations of motion, for our system, and calculate the eigenfrequencies of small fluctuations about the perfectly ordered state. We thereby demonstrate that the ordered state is dynamically stable. In section (<ref>), we calculate the correlation functions implied by our equations of motion, and use them to show that long ranged polar order is robust against noise in these systems. In section (<ref>), we show that the density-density correlation functions found in section (<ref>) imply giant number fluctuations (equation (<ref>)). We also show that those giant number fluctuations depend on the shape of the “counting volume", as well as the mean number of particles it contains, in contrast to equilibrium systems, in which the mean number alone determines the variance. In section (<ref>), we calculate the velocity fields induced in the passive fluid by the active particles, and the correlations of those fields. We also consider the motion of tracer particles in the passive fluid, and thereby derive the anomalous diffusion (equations (<ref>)-(<ref>).)  Finally, in section (<ref>), we use a dynamical renormalization group analysis   to show that these linear results are asymptotically exact at long wavelengths. That is, we show that all non-linear terms allowed in the hydrodynamic equations of motion are irrelevant, in the renormalization group sense, in the long-wavelength limit.
In appendix  (<ref>), we solve the bulk Navier-Stokes equation to relate the passive fluid velocity throughout the bulk to the fluid velocity on the liquid-solid interface. In (<ref>), we derive the stability limit for the ordered state. In appendix (<ref>), we evaluate an integral that arises in the calculation of the real-space correlation function.

§ SUMMARY OF THE PRINCIPAL RESULTS

The variables of our hydrodynamic theory are two time-dependent fields that are confined to the plane of a flat interface between a solid substrate and a three-dimensional (3D)
passive, isotropic, bulk fluid.  We choose our coordinates $(x,y,z)$ so that the interface sits at $z=0$, with the solid substrate underneath  (i.e., $z<0$), and the bulk fluid above (i.e., $z>0$). We will denote two dimensional position within that plane by a two component vector $\brp\equiv(x,y)$. One of the hydrodynamic fields is the conserved density fluctuation $\delta\rho(\brp,t)\equiv\rho(\brp,t)-\rho_0$ of the areal density $\rho(\brp,t)$ about its mean value $\rho_0$.

Our second hydrodynamic variable reflects the polar order. By “polar", we mean each particle (labeled by an index $\alpha$, where $\alpha$ runs from $1$ to the number of particles $N$), can have associated with it a unique unit vector $\hp_\alpha$ (imagine, e.g., a collection of arrows, with  $\hp_\alpha$ pointing from the tail to the head of arrow $\alpha$).

To formulate hydrodynamics, we coarse grain this particle-based unit vector  in to a position dependent vector field $\hp(\brp,t)$ defined on the surface at $z=0$.

Our theory studies small fluctuations around a uniform reference state with polarization ${\bf p}=\hx$ and density $\rho=\rho_0$.

The linear hydrodynamic equations for the transverse orientation fluctuation $p_y$ and density fluctuation $\delta\rho$  for a polar-ordered flock at a solid-liquid interface will be  systematically derived in section <ref> below. They are obtained by by linearizing around a uniform reference state with polarization ${\bf p}=\hat x$ and density $\rho=\rho_0$.    In Fourier space, these   equations  are
∂_tδρ=- iv_ρ[q_xδρ+ρ_c q_yp_y] + i·f_ρ ,

∂_tp_y =- iv_pq_xp_y -γ(q^2+q_y^2 q)p_y -(γ_ρρ_c)(q_xq_y q)
-iσ_t q_yδρ+ f_y  .

Here, $v_p$, $v_\rho$, $\gamma$, $\gamma_\rho$, and $\rho_c$ are phenomenological parameters of our model, the first four of which have the dimensions of speed, while $\rho_c$ has the dimensions of (areal) density. The parameters $\gamma$ and $\gamma_\rho$ are damping
coefficients that play roughly the same role in our problem as the viscosities play in a simple fluid[38]. The speeds $v_\rho$ and $v_p$ are fundamentally active parameters, with no analog in equilibrium fluids, which arise from the self-propulsion of the active particles.

The Gaussian noises ${\bf f}_\rho$ and $f_y$ have zero mean and variances

\begin{eqnarray}
&&\langle f_{\rho i} ({\bf r}_{_\perp},t)f_{\rho j} ({\bf r}_{_\perp}',t')\rangle=2 D_\rho\delta_{ij}^s\delta(
{\bf r}_{_\perp} - {\bf r}_{_\perp}') \delta(t-t')\,\nonumber\\\\
&&\langle f_{ y} ({\bf r}_{_\perp},t)f_{ y} ({\bf r}_{_\perp}',t')\rangle=2 D_p\delta(
{\bf r}_{_\perp}-{\bf r}_{_\perp}')
\delta(t-t')\,,
\end{eqnarray}

where ${\bf r}_{_\perp}=(x,y)$ denotes the position of a point on the plane of the active fluid layer, and the parameters $D_{p,\rho}$ are constants with the dimensions of diffusion coefficients.

Fluctuations about  the uniform ordered state are underdamped; that is, they both propagate and decay. The propagation is non-dispersive; that is, the wavespeed is independent of the wavenumber $q$, as in a simple compressible fluid. However, the wavespeed is direction dependent, as will be discussed in more detail below.

Even more strangely, the decay rate of these fluctuations scales completely differently from that of dry active matter; indeed, it scales linearly with $q$, just like the real part. This means that, unlike almost every other physical system that exhibits non-dispersive propagating modes [36, 37], the “quality factor" $Q$ of our system, defined as the ratio of the real part of $\omega(\bq)$ to its imaginary part, does not diverge as wavenumber $q$ goes to zero. Instead, it approaches a finite, system-dependent constant, independent of $q$, as $q\to0$.

This unusual damping
in turn leads to  the aforementioned  giant number fluctuations of the active particles,  as given by equation (<ref>).

We also find that, when the active particles at the interface are in their ordered state, they “stir" the passive fluid above them.
The components $\langle v_x^2(\brp,z)\rangle$ and $\langle v_y^2(\brp,z)\rangle$ of the mean squared velocity thereby induced in the passive fluid  are inversely proportional to the distance $z$ from the solid-fluid interface, whereas $\langle v_z^2(\brp,z)\rangle$ is inversely proportional to the cube of the distance $z$, i.e, it falls of as $1/z^3$.  Here, the $x$, $y$, and $z$ axes are respectively the direction of the polarization of the active particles, the in-plane direction orthogonal to that, and the normal to the interface, as illustrated in figure (<ref>).

These predictions for the passive fluid velocity correlations could be tested experimentally by tracking neutrally buoyant passive tracer particles in the passive fluid. Such particle tracking thus provides a probe of the fluctuations of the active particles in the ordered state.

These velocity fluctuations in the passive fluid also exhibit long ranged spatio-temporal correlations, which decay as $1/t$  for the projection of the velocity on the interface. The projection of the passive fluid velocity  along $z$ (perpendicular to the surface) decay as $1/t^3$.

These correlations in turn lead to anomalous diffusion of neutrally buoyant passive particles in the $x$- and $y$-direction, with variances of the displacements obeying

\begin{eqnarray}
&&\langle(x(t)-x(0))^2\rangle =2D_xt\left[\ln\left({v_0t\over z_0}\right)+O(1)\right]\,\,, \,\,\,\,\,\, t\ll{z_0^2\over D_z} \,,\nonumber\\
&&\langle(y(t)-y(0))^2\rangle =2D_yt\left[\ln\left({v_0t\over z_0}\right)+O(1)\right]\,\,, \,\,\,\,\,\, t\ll{z_0^2\over D_z} \,,\nonumber\\
&&\langle(x(t)-x(0))(y(t)-y(0)\rangle = 0\,\,, \,\,\,\,\,\, t\ll{z_0^2\over D_z} \,.
\label{anomdifshort}
\end{eqnarray}

The cross-correlation function in the third line of (<ref>) vanishes because the system is symmetric under inversion of $y$. Here, $v_0$ is a system-dependent characteristic speed associated with the active particles (roughly speaking, it is the speed at which they move over the surface due to their self-propulsion),  $z_0$ is the initial distance of the neutrally buoyant particle from the surface, and $D_x,\,D_y$ are   diffusion constants of the same order of magnitude as the diffusion constant in the $z$ direction.

The results (<ref>) hold for times $t$ short compared to the time it takes for the
neutrally buoyant particle to diffuse in the $z$ direction a distance comparable to $z_0$; that is,
for $t\ll{z_0^2\over D_z}$, where $D_z$ is the diffusion constant in 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\gg{z_0^2\over D_z}$, after which the neutrally
buoyant particle will have diffused in the $z$ direction a distance much larger than $z_0$, the
behavior of $\langle(x(t)-x(0))^2\rangle$, and $\langle(y(t)-y(0))^2\rangle$ can be obtained by replacing $z_0$ in (<ref>) with  $\sqrt{D_zt}$,
which  is  just the rms distance the particle will have diffused away from the surface in that time.
Doing so, we find that the motion on these much longer time scales is still superdiffusive, but with
precisely half the “superdiffusion" constant. That is,

\begin{eqnarray}
&&\langle(x(t)-x(0))^2\rangle =D_xt\left[\ln\left({v_0^2t\over D_z}\right)+O(1)\right]\,\,, \,\,\,\,\,\, t\gg{z_0^2\over D_z} \,,\label{anomdiflongx}\nonumber\\
\\
&&\langle(y(t)-y(0))^2\rangle =D_yt\left[\ln\left({v_0^2t\over D_z}\right)+O(1)\right]\,\,, \,\,\,\,\,\, t\gg{z_0^2\over D_z} \,,\label{anomdiflongy}\nonumber\\\\
&&\langle(x(t)-x(0))(y(t)-y(0)\rangle = 0 \,\,, \,\,\,\,\,\, t\gg{z_0^2\over D_z} \,,
\label{anomdiflongxy}
\end{eqnarray}

Again, as in (<ref>) above, the cross-correlation in the third line of (<ref>) vanishes because of symmetry under inversion of $y$. Diffusion in the $z$-direction remains conventional. Furthermore, the diffusion  constant for diffusion in the $z$-direction is independent of height $z$. This set of predictions could also be tested experimentally by particle tracking of neutrally buoyant    tracer particles in the passive fluid.

Particles denser than the passive fluid, which therefore sediment, will also be affected by this activity induced flow. We find that particles sedimenting from an initial height $z_0$ will, when they reach the surface, be spread over a region of dimensions $\sqrt{\langle (x(z=0)-x(z=z_0))^2\rangle}$ and $\sqrt{\langle (y(z=0)-y(z=z_0))^2\rangle}$ in the $x$ and $y$ directions, respectively, with
⟨(x(z=0)-x(z=z_0))^2⟩=2D_x(z_0v_sed)ln(v_0v_sed) v_sed≪v_0  ,

⟨(y(z=0)-y(z=z_0))^2⟩=2D_y(z_0v_sed)ln(v_0v_sed) v_sed≪v_0  ,

⟨(x(z=0)-x(z=z_0))(y(z=0)-y(z=z_0)⟩= 0 v_sed≪v_0   ,

where $v_0$ is roughly the mean speed of the active particles, and $v_{\rm sed}$ is  the speed at which the sedimenting particles sink. These results only hold when the sedimenting speed $v_{\rm sed}\ll v_0$; in the opposite limit, the lateral diffusion of a sedimenting particle becomes conventional.

Once again, these predictions should be readily testable in particle tracking experiments.

We turn now to correlations of the active particles themselves. These are determined by the  characteristic frequencies of the fluctuations. We find that these eigenfrequencies $\omega(\bq)$  are complex, and given by
ω()=Υ(θ_) q  .

Here $\Upsilon$ is a complex   (and complicated!) function of the angle $\theta_\bq$ between the wavevector $\bq$ and the direction of spontaneous polarization of the active particles, and is the solution of
a parameter and direction dependent quadratic equation   given in section <ref>.

Allowing the fluctuations of the stable modes to be driven by an additive white noise (which will arise from both noisy self-propulsion of the active particles, and thermal fluctuations), we find the orientational correlation function

C_pp(q,ω)≡⟨|p_y(, ω)|^2⟩=(1q^2)F_pp((ωq), θ_)  ,

where  $F_{pp}(\omega/q, \theta_\bq)$ is a scaling function whose precise form is given in Section <ref>. We also find that the peak positions of the scaling function, which we give in detail in Section <ref>, agree with those found in  the experiments of [29] on Quinke rotators.

The density-density correlation function and the orientation-density correlation function obey  similar scaling laws:
C_ρρ(q,ω)≡⟨|p_y(, ω)|^2⟩=1q^2F_ρρ((ωq), θ_)  ,

C_pρ(q,ω)≡⟨ρ(q,ω) p_y(-q,-ω)  ⟩
=1q^2F_pρ((ωq), θ_)  ,

with slightly different scaling functions $F_{\rho\rho}(\omega/q, \theta_\bq)$ and $ F_{p\rho}(\omega/q, \theta_\bq)$, whose precise forms   are  given by equations (<ref>) and (<ref>) of section <ref>.

Integrating these spatio-temporally Fourier-transformed correlation functions over all frequencies
$\omega$ gives  surprisingly simple expressions for the equal-time correlation functions.   Specifically,
C_pp()≡⟨|p_y(, t)|^2⟩=∫_-∞^∞dω2πC_pp(,ω)=f_pp(θ_)q  ,

where the function $f_{pp}$ depends only on $\theta_\bq$ and not $\bf q$, whose precise form is discussed later in Section <ref>.

For the equal-time density autocorrelation function, we obtain
C_ρρ()≡⟨|ρ(, t)|^2⟩=f_ρρ(θ_)q

where $f_{\rho\rho}$ is another function of $\theta_\bq$, to be given later.

Finally, the equal-time cross-correlation function is given by
C_pρ()≡⟨p_y(, t)ρ(-,t)⟩=f_pρ(θ_)q

where the direction dependence is given by the function $f_{p\rho}$ that depends only on $\theta_\bq$.

Fourier transforming these in space leads to equally simple expressions for the real space, equal-time correlation functions,  all of which scale as $1/r$ with anisotropic amplitudes.

That is:
C_pp()   ≡   ⟨δ(+,t)δ(,t)⟩≈⟨p_y(+,t)p_y(,t)⟩
=   f_pp(θ_)r  ,

where, the function $f_{pp}(\theta_\br)$ depends only on $\theta_\br$. the angle between $\br$ and the direction of the mean polarization, and not on the magnitude $r\equiv|\br|$. The  precise form of $f_{pp}(\theta_\br)$ is given in equation blah  of section blah.

For the real space equal-time density autocorrelation, we obtain

where $f_{\rho\rho}$ is another function of $\theta_\br$, to be given later.

Finally, the equal-time polarization-density cross-correlation function is given by

where the direction dependence is given by the function $f_{p\rho}$ that depends only on $\theta_\br$.

These predictions could also be tested experimentally in systems in which the active particles can be imaged, like those of [27, 28].

§ DERIVATION OF THE  HYDRODYNAMIC THEORY

§.§ Identifying the hydrodynamic variables

We consider active polar particles of  areal density $\rho(\brp,t)$

confined to  a flat interface between a solid substrate and a three-dimensional (3D)
passive, isotropic, bulk fluid.  We choose our coordinates $(x,y,z)$ so that the interface sits at $z=0$, with the solid substrate underneath  (i.e., $z<0$), and the bulk fluid above (i.e., $z>0$). By “polar", we mean each particle (labeled by an index $\alpha$, where $\alpha$ runs from $1$ to the number of particles $N$), can have associated with it a unique unit vector $\hp_\alpha$. (Imagine, e.g., a collection of arrows, with  $\hp_\alpha$ pointing from the tail to the head of arrow $\alpha$.)

To formulate hydrodynamics, we coarse grain this particle-based unit vector  in to a position dependent vector field $\hp(\brp,t)$ defined on the surface at $z=0$.
In principle, this coarse graining can lead to fluctuations in the magnitude of the coarse grained
$\hp(\brp,t)$ as well as its direction.
However, since  we are considering  a system   in an orientationally ordered state,  not at a continuous transition from an ordered to a disordered state, fluctuations of the  magnitude $|\hp|$ of $\hp$ relax fast; i.e., the relaxation time-scale remains finite in the long wavelength limit.  This is because fluctuations in the magnitude of $|\hp|$, unlike fluctuations in its direction, are not Goldstone modes associated with the spontaneous breaking of continuous rotation invariance. Therefore, there is no symmetry reason for them to be slow, so we  expect them to be fast generically. Hence,  in the spirit of hydrodynamics, we ignore these magnitude fluctuations, and set $\hp^2=1$ without any loss of generality.  A renormalization group analysis[41] confirms this heuristic argument.

We choose   our coordinate system such that $p_x=1$ is the reference polar
ordered state (see figure (<ref>).

In the presence of friction from the substrate, there is no momentum
conservation  on the surface, so  the only conserved
variable on the surface of relevance for the present problem is the active particle
number. The hydrodynamic theory that we develop below thus includes
the active particle number density $\rho (\brp,t)$ and the orientational order parameter ${\hat{\bf p}}(\brp,t)$ fluctuations
as the relevant slow  or “hydrodynamic" variables.  We also include the bulk fluid velocity $\bv(\br, t)$, which is defined throughout the semi-infinite three dimensional space above the surface, since in that space momentum (which is equivalent to velocity in the limit of an incompressible bulk fluid) is conserved. However, we will work in the Stokesian limit, in which viscous forces dominate inertial ones; this has the effect of “enslaving" the bulk velocity to the hydrodynamic variables $\rho (\brp,t)$ and  ${\hat{\bf p}}(\brp,t)$ on the surface, as we will show below.

§.§ Hydrodynamic equations of motion

We will formulate the hydrodynamic equations for these variables by expanding their equations of motion phenomenologically in powers of fluctuations of both fields $\hp$ and $\rho$ from their mean values, and in spatio-temporal gradients. In so doing, we must respect all symmetries and conservation laws of the underlying dynamics. 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 term consistent with the symmetries and conservation laws is allowed, and therefore will, in general, be present.

We will consider systems with underlying rotational invariance in the plane of the surface. That is, we will focus on the case in which there is no a priori preferred direction in that plane; the active particles are free to spontaneously break this symmetry, by choosing a direction in which to spontaneously point.

With these considerations in mind, we will now derive the hydrodynamic equations of motion for the fields $\rho(\brp,t)$, $\hp(\brp,t)$, and the bulk velocity velocity $\bv(\brp, z, t)$.

We begin with the active particle density $\rho(\brp,t)$.

§.§.§ Equation of motion for the active particle density $\rho$

The conservation of the active  particles implies that they obey a continuity equation:
∂_tρ+ ∇_s·J_ρ   =    0 ,
where ${\boldsymbol \nabla}_s\equiv {\hat{\bf x}}\partial/\partial x +
{\hat{\bf y}}\partial/\partial y$ is the 2D gradient operator, with ${\hat{\bf x}}$ and
${\hat{\bf y}}$ the unit vectors along the $x$ and $y$ axis respectively. The active particle
current ${\bf J}_\rho$ is generated by  friction with the  bulk fluid
velocity ${\bf v}(\br=\brp, z=0)$, plus additional diffusive and active currents.  We will phenomenologically expand all of these effects to leading order in powers of the bulk velocity  evaluated at the surface ${\bf v}(\brp, z=0)$, and gradients, keeping all terms consistent with the underlying rotation invariance. In practice, this means we can make the vector ${\bf J}_\rho$ only out of vectors the active particle configuration itself chooses, that is, out of gradients, the surface velocity $\bv_s(\brp, t)$,  by which we mean the bulk velocity  evaluated at the surface; i.e.,  $\bv_s(\brp, t)\equiv{\bf v}(\brp, z=0,t)$, and the polarization $\hp(\brp, t)$. These constraints force ${\bf J}_\rho $ to take the form:

J_ρ()    =    ρ_e(ρ, |_s|)  v_s(x,y,z=0) + κ(ρ,  |_s| ) -D_ρ∇_s ρ-D_ρ2(·∇_s)ρ-f_ρ(x,y)  .

to leading order in gradients, where $D_\rho>0$ is the diffusion constant, and $D_{\rho 2}$ is an anisotropic diffusion coefficient.  Any other symmetry-permitted vector terms involving additional gradients, e.g., $\nabla_s^2 {\bf \hat p}$ would be irrelevant in the long wavelength limit. Note also that the first two terms on the rhs of (<ref>) reflect the fact that the  the “active particle velocity” ${\bf J}_\rho/\rho$ will, in general, be different from the velocity ${\bf v}_s$  of the bulk fluid at the surface, due to the activity.  The existence of a non-zero $D_{\rho 2}$  reflects the fact that, once the underlying rotation invariance is broken by the development of a spontaneous polarization $\hp$, diffusion along the direction of $\hp$ need not, and, in general, will not, proceed at the same speed as diffusion perpendicular to that direction.  The  “effective density” $\rho_e(\rho, |\bv_s|)$ would simply be $\rho$ itself in a Galilean invariant model, but since the solid surface breaks Galilean invariance, it can in general depend non-linearly on $\rho$, and even on the magnitude $|\bv_s|$ of the surface velocity. All of these dependences disappear in the linear model; furthermore, our analysis in section (<ref>) of the non-linear model   shows that these dependences do not affect any of our predictions for the long distance behavior of the system. More precisely, they can all be absorbed into a suitable, finite “renormalization” of the parameters of the linear theory.  In renormalization group jargon, they are “irrelevant”. We have only included them in the current (<ref>) to make our starting model  more general.

The factor   $\kappa (\rho, |\bv_s|)$ is an active parameter reflecting the self-propulsion of the particles through interaction with the solid substrate. It can also in principle depend  on $\rho$ and  the magnitude $|\bv_s|$ of the surface velocity. As for $\rho_e$, all of these dependences disappear in the linear model, and prove to also be irrelevant in the RG sense just described.

For completeness,
we have included in this current a zero-mean Gaussian white noise  ${\boldsymbol f}_\rho$ that we take to be delta-correlated in space and time:

\begin{equation}
\langle f_{\rho i} (x,y,t)f_{\rho j} (x',y',t')\rangle=2 D_\rho\delta_{ij}^s\delta(
x- x') \delta(y- y') \delta(t-t')\,.
\label{rho noise}
\end{equation}

This noise is added to model the stochastic nature of the density dynamics. In the equilibrium limit, this would have been a thermal noise representing the effects of a finite temperature and with a variance proportional to the diffusivity, as required by the Fluctuation-Dissipation theorem [42]. In an active system  such as ours,  there is no constraint like the fluctuation-dissipation theorem or the Einstein relations relating the noise variance  (<ref>)  to the diffusivity. The white noise (<ref>) in (<ref>) reappears in (<ref>) as a conserved noise, consistent with the density being a conserved variable; see Eq. (<ref>) below.
However, as we will show later, this noise proves to be irrelevant in the long wavelength limit compared to the polarization noise we will introduce later. The diffusion constants $D_\rho$ and $D_{\rho2}$ prove to be irrelevant in the long-wavelength limit as well, as we shall see.

Using (<ref>) in (<ref>), the continuity
equation (<ref>)
can be written

∂_tρ=   -   ·∇_sκ(ρ) -κ(ρ)∇_s·
+    ∇_s·[D_ρ∇_sρ+D_ρ2(·∇_s)ρ]

-   ρ_e(ρ, |_s|)∇_s·_s-_s·∇_sρ_e(ρ, |_s|)+
∇_s·f_ρ .

§.§.§ The bulk velocity $\bv(\brp, z, t)$

In calculating the bulk velocity $\bv(\brp,z,t)$, we will assume  the bulk fluid is  in the extreme “Stokesian” limit, in which inertia is negligible relative to viscous drag. This should be appropriate for most systems in  which the active particles are microscopic, since the Reynolds' number will be extremely low for such particles.

The   3D incompressible bulk velocity field ${\bf
v}=(v_i,v_z),\,i=x,y$ satisfies
the 3D Stokes' equation

\begin{equation}
\eta\nabla^2_3v_\alpha (\brp,z)= \partial_\alpha \Pi (\brp,z)  \,,
\label{Stokes}
\end{equation}

where the incompressibility constraint

$\Pi$, the hydrostatic pressure.
In these expressions (<ref>), (<ref>), ${\boldsymbol \nabla}_3\equiv {\hat{\bf x}} \partial/\partial x +
{\hat{\bf y}}\partial/\partial y+
{\hat{\bf z}}\partial/\partial z$ is the  full three-dimensional  gradient operator, with ${\hat{\bf x}}$, ${\hat{\bf y}}$,  and
${\hat{\bf z}}$ the unit vectors along the $x$, $y$, and $z$  axes respectively.

We could, in principle, add a random Brownian noise to this Stokesian force balance equation (<ref>) to reflect thermal fluctuations in the non-zero temperature bulk fluid. However, since the bulk fluid is passive - that is, equilibrium, we know that the net effect of such a noise term on the surface polarization must be to introduce a purely equilibrium random noise in the polarization equation of motion. So we will instead take the simpler approach of simply adding such a thermal noise by hand to the active noises that will also appear in that equation.

In the absence of such a noise term, the noiseless Stokes equation (<ref>) can be solved exactly for the bulk velocity $\bv(\brp,z,t)$ in terms of the surface velocity $\bv_s(\brp,t)$. If we Fourier expand the surface velocity:

where $(L_x,L_y)$ are the linear dimensions of our (presumed rectangular) surface, then, as we show in Appendix <ref> , the bulk velocity $\bv(\brp,z,t)$  is given by

This is to be
supplemented by the boundary conditions at $z=\infty$, at which all stresses vanish, and at $z=0$, which we will now discuss.
For a bulk passive fluid with a free surface, the boundary condition is given
by the vanishing of the appropriate component of the shear stress.  On the other
hand, for a bulk fluid resting on a solid surface, one imposes the “no-slip”
boundary condition  that  the relative velocity between the solid surface and
the fluid layer in contact vanishes. More generally, if there is a slip
velocity between the solid surface and the 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 with the 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. Generalizing the partial-slip boundary condition [60] for an active surface
fluid, we impose the following active boundary condition at $z=0$:
(,z=0, t) ≡_s(,t)

=   v_a(ρ) (,z=0,t)+ μτ_∥(,z=0,t)  ,

where $v_a(\rho)$ is the spontaneous self-propulsion speed of the active particles relative to the solid substrate  that can depend upon the local active particle density, and the force density ${\bm \tau}_\parallel$ parallel to the surface is given by

τ_∥(,z=0,t)=·σ^bulk (,z=0,t)+ τ^surface(, t)  ,

with ${\bf \cO}$ the operator that projects any vector onto the plane of the surface, and
${\bf \sigma}^{{\rm bulk}}$ the bulk stress tensor. The projection operator ${\bf \cO}$  is given by
_ij=δ_ij-δ_izδ_jz  ,

while the bulk stress tensor is
given by

σ^bulk_ij=η(∂_iv_j+∂_jv_i) -Πδ_ij  .

The $\eta$ term in (<ref>)  is simply the viscous stress in the bulk fluid, while the pressure term $\Pi\delta_{ij}$ drops out of the force density ${\bm \tau}_\parallel$ parallel to the surface due to the projection operator ${\bf \cO}$.

The form of this boundary
condition is dictated by symmetry. Essentially, any vector in the two dimensional plane of the solid surface that is spontaneously chosen by the system (as opposed to specified a priori, is allowed, and, therefore, will be present. Because the projection  $\cO{\hat \bz}\cdot\bm{\sigma}^{{\rm bulk}}$ is such a vector, a term proportional to it is allowed, and, hence, will generically be present in the boundary condition (<ref>).

Physically, we expect contributions to the surface stress from the bulk fluid stress since this stress is also experienced at the interface where the active particles are situated. So we would expect this stress in general to influence the active particles, and, therefore, the boundary condition on the bulk fluid velocity at the interface on which these particles sit.
Note that if we set both the active surface stress $\tau^{\rm surface}$ and the self-propulsion speed $v_a$ to zero, equations II.8 to II.11 reduce to the conventional
partial slip boundary condition of passive fluid dynamics[60].

Requiring that the bulk stress tensor (<ref>) be divergenceless (as it must be in the Stokesian limit, since, by definition, in that limit the viscous stresses dominate the inertial ones) leads to
the Stokes' equation (<ref>) for the bulk flow. Note that $\tau^{\rm surface}$ exists only at the interface at $z=0$ and hence does not appear in the Stokes' equation for the bulk flow.  The  surface force ${\bm \tau}^{{\rm surface}} (\brp, t)$, when expanded to leading order in the polarization $\hp$ and gradients, must, by symmetry,  take the form[61, 62]
τ^surface_i=ζ_1(ρ)·∇_Sp_i +ζ_2(ρ)p_i∇_S·+p_i·∇_Sζ(ρ) - ∂_iP'_s(ρ)  .

In (<ref>), the $\zeta_{1,2}$ and $\zeta$ represent active stresses, and $P'_s$ is a surface osmotic pressure (not to be confused with the bulk pressure $\Pi$ in the Stokes equation (<ref>)), both of which depend on the density of the active particles.

Inserting (<ref>) and (<ref>) into our expression (<ref>) for the force density ${\bf \tau}_\parallel$ parallel to the surface, and using the result in our  active boundary condition (<ref>), gives,
for $i=(x,y)$,

\begin{equation}
v_{si}(\brp,t)=v_a(\rho) p_i(\brp,t) + \zeta_1(\rho)\hp\cdot\nabla_Sp_i +\zeta_2(\rho)p_i\nabla_S\cdot\hp +p_i\hp\cdot\nabla_S\zeta(\rho) + \mu\eta\bigg(\frac{\partial
v_i(\brp,z,t)}{\partial z}\bigg)_{z=0}- \partial_iP_s  \,,
\label{activebc}
\end{equation}

where we have defined $P_s\equiv\mu P'_s$, and used the fact that the velocity normal to the surface vanishes at the surface (i.e.,  $v_z(z=0)=0$) [61].  Note that the velocity $v_i(\brp,z,t)$ is $i$'th component  the full bulk velocity.

The term proportional to $\partial_zv_i$ can be given another interpretation: if we define a length $a\equiv\mu\eta$, $ \mu\eta\bigg(\frac{\partial
v_i}{\partial z}\bigg)_{z=0}=a\bigg(\frac{\partial
v_i}{\partial z}\bigg)_{z=0}$ is simply the velocity of a fluid element a distance $a$ above the surface. Hence, we can alternatively interpret this term as modeling  active particles of finite thickness $a$ on the surface   which are passively convected by the local fluid velocity at that height.
For a system in thermal  equilibrium, $v_a=0=\zeta_{1,2}(\rho) = \zeta(\rho)$, and (<ref>) reduces to the well-known equilibrium partial slip boundary condition.

§.§.§ Equation of motion for the polarization $\hp$

We now turn to the equation of motion for $\hp$. As the active particles are polar, the system is
not invariant under $\hp\rightarrow
-\hp$ symmetry. This means that terms even in $\hp$ are allowed in our phenomenological expression for $\partial_t\hp$.

The most general equation of motion for $p_k$  allowed by symmetry is therefore:
T_ki(αv_si -λ_pv(_s·∇_s)p_i +h_i γ_1+
(ν_1-1 2)p_j∂_iv_sj
+(ν_1+1 2)(·∇_s)v_si
-∂_iP_p(ρ)+f_i) ,

where we have defined the transverse projection operator

which projects any vector orthogonal to $\hp$. Its presence in (<ref>) insures that the fixed length condition $|\hp|=1$ on $\hp$ is preserved.

In (<ref>), the “molecular field” $\bh$ is given by
=-δFδ  ,

where the Frank free energy is
F=12∫d^2r_⊥ [K_1(∇_S·)^2+K_3|∇_S×|^2]  .
Here, $K_1$ and $K_3$ are respectively the Frank splay and bend elastic modulii.  This $h_i/\gamma_1$-term, together with the Frank free energy $F$, gives the damping terms in (<ref>) below. Notice that the more familiar   classical alignment-free energy term $K|\nabla_S\hp|^2$ is already implicitly included
in (<ref>) separately,  as can be seen by using the vector calculus identity
|∇_S×|^2=|∇_S|^2-(∇_S·)^2-∇_S·W  ,

where ${\bf W}\equiv(\hp\cdot\nabla_S)\hp-\hp(\nabla_S\cdot\hp)$, in (<ref>). Doing so, and dropping the surface term arising from the total divergence in (<ref>), we can rewrite (<ref>) in the form
F=12∫d^2r_⊥ [(K_1-K_3)(∇_S·)^2+K_3|∇_S|^2]  .
In this form, it is explicit that we have included the usual alignment-free energy term proportional to $|\nabla_S\hp|^2$ in (<ref>).

The reasoning that leads to (<ref>) is essentially identical to that leading to our earlier expression (<ref>) for the active boundary condition: we want to include every vector spontaneously picked by the system, provided that vector is also leading order in an expansion in powers of spatial gradients. It is straightforward to see that (<ref>) does so. It contains all of the vectors contained in (<ref>), except for an explicit term proportional to $\hp$ itself. We need not include such a term in (<ref>), since it would be “killed" by he projection operator $T_{ki}$. All the other terms in (<ref>) are vectors spontaneously picked by the system. Furthermore, they all have independent coefficients, since there is no symetry argument that would lock them ito any particular relation.

Note that terms proportional to some a priori specified vector, say $\hx$, are not allowed in (<ref>), since they violate rotation invariance by singling out a special direction that is not spontaneously chosen by the system (i.e., by any of the vectors $\hp$, $\bv$, etc., or by spatial gradients).

The fact that we are performing a gradient expansion allows us to drop terms proportional to, e.g., $\nabla_S^2 v_{si}$, $\nabla_S^4 v_{si}$, and so on, in (<ref>). While such terms are certainly allowed by symmetry, they are irrelevant at long distances relative to the $\alpha v_{si}$ and $p_j\partial_iv_{sj}$ terms already present  in (<ref>).

When all such symmetry forbidden and irrelevant in the gradient expansion terms are dropped, one is left with (<ref>) as the only possible equation of motion for $\hp$.

To summarize, the equation of motion (<ref>)
includes all of the relevant nonequilibrium terms allowed by symmetry, in addition to the usual
equilibrium terms and the convective covariant derivative of $\hp$. Here,
the presence of the solid substrate underneath breaks any Galilean
invariance, and hence a term directly proportional to $v_i$ in (<ref>) is
permissible: indeed the
$\alpha$-term in (<ref>) clearly breaks the Galilean invariance. The presence
of this terms implies flow alignment or antialignment of the polarization $\hp$ depending
upon the sign of $\alpha$  in
Eq. (<ref>). The “polarization pressure” $P_p(\rho)$ is an additional function of the density,
independent of the “osmotic pressure” $P_s(\rho)$ introduced earlier.
Furthermore,  the $\lambda$ term represents active self advection. Finally, the terms proportional to $\nu_1$ are “flow alignment terms”, identical in form to those found in nematic liquid crystals[63].

We have also added to the equation of motion (<ref>) a white noise ${\bf f}$ with statistics

\begin{equation}
\langle f_{ i} ({\bf r}_{_\perp},t)f_{ j} ({\bf r}_{_\perp}',t')\rangle=2 D_p\delta_{ij}\delta(
{\bf r}_{_\perp}-{\bf r}_{_\perp}')
\delta(t-t')\,.
\label{p noise}
\end{equation}

As discussed earlier, this noise also incorporates equilibrium contributions from thermal fluctuations of the bulk fluid. These must, by the fluctuation-dissipation theorem, be spatiotemporally white, as we have assumed in (<ref>), since the equilibrium dynamics of $\hp$ is local in space and time. The actual noise strength $D_p$ in (<ref>) is larger than the equilibrium value (which is proportional to $k_BT$)  due to active contributions to the noise.

§.§.§ Summary of the equations of motion

Our hydrodynamic model, then, is summarized by the equations of motion (<ref>) and (<ref>) for $\rho$  and $\hp$, respectively, and the solution (<ref>) of the Stokes equation (<ref>) for the bulk velocity field $\bv(x,y,z,t)$obtained with the boundary condition (<ref>). When considering fluctuations, we will also need the noise correlations (<ref>)  and (<ref>) .

§.§ Linearization of the equations of motion

We now show that linearizing Eqs. (<ref>) and (<ref>) about a uniform reference state with polarization ${\bf p}=\hat x$ and density $\rho=\rho_0$ produces Eqs. (<ref>) and (<ref>) in section <ref>.

The equations of motion and boundary conditions found in the previous section have an obvious spatially uniform, steady state solution:

where we have defined

and have chosen the $\hx$ axis of our coordinate system to be along the (spontaneously chosen) direction of
polarization, as illustrated in figure (<ref>).

As a first step towards understanding fluctuations about this steady state, we will  write

_s(,t)=(v_0+δv_sx(,t))+v_sy(,t)  ,

and expand the equations of motion <ref>) and (<ref>) for $\rho$  and $\hp$, and the boundary condition (<ref>), to linear order in $\delta\rho$ and $p_y$.
We will obtain the bulk velocity $\bv(\brp,z,t)$ from the surface velocity $\bv_s(\brp,t)$ using our solution (<ref>) of the Stokes equation.

We will expand the surface osmotic pressure of the particles $P_s(\rho)$  to linear order in $\delta\rho$:
P_s(ρ)=σδρ .

We will also expand the “surface polarization
pressure” $P_p(\rho)$ to linear order in $\delta\rho$ (and will show later that higher powers of $\delta\rho$ are irrelevant at long wavelengths):
P_p(ρ)= σ_p δρ ,

and likewise expand the $\zeta$'s in equation (<ref>) . To linear order, it is sufficient to take

ζ(ρ)=ζ(ρ_0)+ζ̅δρ        .

Finally, we expand $v_a(\rho)$, $\kappa(\rho)$ and $\rho_e(\rho)$ to linear order in $\delta\rho$:

ρ_e(ρ)=ρ_1+ρ_e'δρ        ,

where $\kappa_0$, $\bar\kappa$, $\rho_1$, and $\rho_e'$ are all constants to linear order. All of the other parameters (i.e., the diffusion constants $D_\rho$ and $D_{\rho2}$, $v_a$, $\alpha$, $\lambda$, $\lambda_{pv}$, $\gamma_1$, and $\nu_1$ can to linear order be replaced by their values at $\rho=\rho_0$, $|\bv_s|=v_0$, and treated as constants.

Now using our active boundary condition Eq. (<ref>), the fluctuating velocity components $\delta v_x$ and
$v_y$ can be expressed up to linear order as
δv_sx    =    v'_aδρ+ζ̅∂_xδρ+μη(∂_zv_x)_z=0 +ζ_20∂_yp_y
-σ∂_x δρ ,

v_sy    =    v_0 p_y +μη(∂_zv_y)_z=0  +ζ_10∂_xp_y
-σ∂_yδρ .

These expressions for $\delta v_{sx}$ and $v_{sy}$ are implicit equations, since the bulk velocities $v_x$ and $v_y$ on the right hand side also
depend on the surface velocity through (<ref>) .

We can make them explicit by solving them iteratively. This is most conveniently done in Fourier space, because we can then obtain the iterative solution to lowest order in powers of the wavenumber $q$.

Performing a two-dimensional Fourier transform on  (<ref>) and (<ref>) - that is, Fourier transforming over $\brp$, and using our solution (<ref>)  for the bulk velocity in terms of the surface velocity, we obtain
δv_sx(,t) =v'_aδρ(,t)+i ζ̅q_xδρ(,t) +μηW_x(,t) +iζ_20q_yp_y(,t)
-iσq_x δρ(,t) ,

δv_sy(,t) = v_0 p_y(,t) +μηW_y(,t)  +iζ_10q_xp_y(,t)
-iσq_yδρ(,t)  ,

where we have defined

W(,t)    ≡    (∂_z{[_s(,t)-z(·_s(,t))(+i)]e^-qz})_z=0

=    -q_s(,t)-(·_s(,t))+W_z(,t)

In the second equality of (<ref>), we have not bothered to evaluate $W_z(\bq,t)$, since it does not appear in (<ref>) or (<ref>).

If we now insert our expressions  (<ref>) and (<ref>) for $\bv_s(\bq,t)$ after the second equality in (<ref>), we see that the leading order in $q$ terms (which are  $O(q)$)  come from the $v'_a\delta\rho$ term   in  (<ref>) and the $v_0p_y$ term in  (<ref>). Keeping only those terms implies
_s(,t)≈v'_aδρ+v_0p_y(,t) .

Inserting this into (<ref>) gives
W_x    =    -(q^2+q_x^2 q)v'_aδρ-v_0(q_xq_yq)p_y

W_y    =     -v'_a(q_xq_yq)δρ-v_0(q^2+q_y^2 q) p_y

Inserting these expressions back into  (<ref>) and (<ref>) gives us our final closed form expression for the two components of the surface velocity, written entirely in terms of $p_y$ and $\delta\rho$:
δv_sx(,t)    =   [v'_a+i (ζ̅-σ) q_x -μηv'_a(q^2+q_x^2 q)]δρ(,t) +[iζ_20q_y-μηv_0 (q_xq_yq)]p_y(,t)

v_sy(,t)    =   [ v_0  -μη(q^2+q_y^2 q)   +iζ_10q_x]p_y(,t)
-[iσq_y +μηv'_a(q_xq_yq)]δρ ,

Note the non-analytic character of the $\mu\eta$ terms in these expressions; this reflects the long-ranged hydrodynamic interaction between active particles on the surface mediated by the bulk passive fluid. Indeed, it is only through these terms that the presence of the bulk fluid makes itself felt. Note also that these terms are real; we'll see in a moment that this makes them damping terms. They are also the same order in $q$ as the imaginary $\zeta$ and $\sigma$ terms, which are associated with propagation. This is the origin of the peculiar property of our system that damping and propagation are the same order in wavevector (or, equivalently, that the quality factor $Q$ of the normal modes of this system is finite and independent of wavenumber $q$ at small $q$). This should be contrasted with, e.g., a simple bulk equilibrium fluid, for which the propagating terms are $O(q)$, while the damping terms are $O(q^2)$.

Note that the parameter $\zeta_{20}$ has dropped out at this point in our calculation, because it only leads to terms of higher order in $q$ than we have kept here.

With these expressions (<ref>) and (<ref>) for the surface velocity  in terms of $p_y$ and $\delta\rho$
in hand, we can now derive closed equations of motion for $p_y$ and $\delta\rho$. First, we must linearize the general equations of motion (<ref>) and (<ref>), using the linearizations (<ref>)- (<ref>) for the pressures and $\zeta$'s, and    (<ref>) and (<ref>) for the density and the polarization. Doing so for the
continuity equation (<ref>), we obtain

∂_tδρ=-κ_0∂_yp_y -(κ̅+v_0ρ'_e)∂_xδρ-ρ_1∇_s·_s+ [D_ρ∇_s^2ρ+D_ρ2 ∂_x^2]δρ+

∇_s·f_ρ .

Fourier transforming equation (<ref>), and keeping only terms to leading order in $q$, we find
∂_tδρ(,t) =-iκ_0q_yp_y -i(κ̅+v_0ρ'_e)q_xδρ-iρ_1·_s+

i·f_ρ .

Setting $k=y$ in the equation of motion (<ref>) for $\hp$, and linearizing the resulting equation gives
∂_tp_y=α(v_sy-v_0p_y)-(λ_pvv_0+λ)∂_xp_y+(ν_1-1 2)∂_yδv_sx
+(ν_1+1 2)∂_xv_sy
-σ_p∂_yδρ+1γ_1(K_1∂_y^2+K_3∂_x^2)p_y+f_y  .

Note the appearance of the combination $v_{sy}-v_0p_y$ in the $\alpha$ term on the RHS of this equation. The fact that this particular combination appears is not an accident, but, rather, a consequence of rotation invariance: if we work to zeroeth order in gradients in our solution   for the surface velocity, this term vanishes, as it must, since $p_y$ is a Goldstone mode associated with the spontaneous breaking of the continuous rotational symmetry. It is because of this exact cancellation that we needed to evaluate $\bv_s$ to higher order in gradients of $p_y$ in (<ref>) and (<ref>), which is why we had to keep higher order gradient terms in that expression.  Note that no such cancellation happened when we calculated the vector ${\bf W}$ eqn. (<ref>) earlier, so there we could truncate the expansion for the surface velocity at zeroeth order in  the gradient term.

Fourier transforming equation (<ref>), and keeping only terms to leading order in $q$, we find that  the molecular field terms proportional to the Frank constants $K_{1,3}$ are higher order in $q$, and so can be dropped, leaving us with
∂_tp_y=α(v_sy-v_0p_y)-i(λ_pvv_0+λ)q_xp_y+i(ν_1-1 2)q_yδv_sx
+i(ν_1+1 2)q_xv_sy
-iσ_pq_yδρ+f_y)  .

Now inserting our expressions (<ref>) and (<ref>) for the surface velocity in terms of $p_y$ and $\delta\rho$ into the Fourier transformed equations of motion (<ref>) and (<ref>), and keeping only terms to leading order in $q$,  gives our final closed form for the linearized equations of motion:

∂_tδρ=- iv_ρ[q_xδρ+ρ_c q_yp_y] + i·f_ρ ,

∂_tp_y =- iv_pq_xp_y -γ(q^2+q_y^2 q)p_y -(γ_ρρ_c)(q_xq_y q)
-iσ_t q_yδρ+ f_y  ,

where we have defined the characteristic velocities
v_p≡λ_pv  v_0+λ-ζ_10α-(ν_1+12) v_0

v_ρ≡κ̅+ρ_e'v_0+ρ_1v'_a ,

the total inverse compressibility $\sigma_t\equiv\sigma_p+\sigma\alpha - \bigg(\frac{\nu_1-1}{2}\bigg)
the characteristic density
and the bulk fluid damping coefficients  $\gamma\equiv\alpha  v_0\mu\eta$ and $\gamma_\rho\equiv\alpha\mu\eta v'_a\rho_c$. Both $\gamma$ and $\gamma_\rho$ have the dimensions of speed. Further  more, all the contributions from $\delta v_x({\bf q},t)$ in (<ref>) in (<ref>) are subleading and hence do not appear in (<ref>), where only the terms leading order in $\bf q$ are kept.

We thus obtain the linearized hydrodynamic equations  (<ref>) and (<ref>),  which are nothing but  equations (<ref>) and (<ref>) mentioned earlier. These
are invariant under $q_y\rightarrow -q_y,\,p_y\rightarrow
-p_y,\rho\rightarrow\rho$, but not invariant under $q_x\rightarrow -q_x$. This
is equivalent to invariance in real space under  $y\rightarrow -y,
p_y\rightarrow -p_y,\rho\rightarrow \rho$, but  with no analogous invariance under
$x\rightarrow -x$.  In these two equations $v_\rho\neq v_p$ in general, due to the abscence of Galilean invariance because of
friction with the solid substrate.
The effective coefficient $\gamma$, when positive (which can be achieved by tuning the signs of the various original model parameters), serves as the effective damping coefficient in the model. Interestingly, the effective damping here is ${\cal O}(q)$, which is far stronger than the ${\cal O}(q^2)$ damping in the linearized Toner-Tu model for flocking [13, 14] in the hydrodynamic limit. Hydrodynamic interactions mediated by the passive, bulk fluid above are responsible for this ${\cal O}(q)$ damping.

An alert reader might wonder how hydrodynamic interactions mediated by the bulk fluid  can dominate friction from the underlying solid substrate,  which, after all, is $O(q^0)$. The reason is that friction with the substrate, while playing a very important role (in particular, it is the principal mechanism limiting the speed $v_0$ of the active particles), does not act to suppress fluctuations in the directions of motion of those particles (or, similarly, the polarization. The leading order damping of such fluctuations, which are the Goldstone modes of our problem, come from the hydrodynamic interactions.

§ MODE STRUCTURE OF THE LINEARIZED EQUATIONS AND STABILITY OF THE UNIFORM STATE

Having written down the equations of motion, we will now analyze the linear stability of the
We  work  in polar coordinates ${\bf
q}=(q\cos\theta_\bq,q\sin\theta_\bq)$,  and
set the noises to zero.
Assuming a time-dependence of the
form $\rho, p_y\sim  \exp (-i\omega t)$, this   leads to
the eigenvalue condition on $\omega$:

ω^2+iω[i(v_p+v_ρ)qcosθ_+γq (1+sin^2θ_)]-i[iv_pqcosθ_+γq (1+  φsin^2θ_)]v_ρqcosθ_-c_0^2q^2sin^2θ_=0  ,

where we have defined
c_0^2≡σ_tρ_c v_ρ
φ≡1-γ_ργ  .

The eigenfrequencies always scale as $q$, independent of all other parameters, as can be seen by inserting the ansatz
ω=Υ(θ_) q

into (<ref>). This leads to a $q$ independent condition on $\Upsilon$:
Υ^2+iΥ[i(v_p+v_ρ)cosθ_+γ(1+sin^2θ_)]-i[iv_pcosθ_+γ(1+ φsin^2θ_)]v_ρcosθ_-c_0^2sin^2θ_=0  ,

thereby proving that $\omega$ scales like $q$.

For stability, we must have ${\rm Im}(\omega)<0$ for all $\bq$, which means ${\rm Im}(\Upsilon(\theta_\bq))<0$ for all $\theta_\bq$. We show in appendix (<ref>)  that this condition is satisfied for sufficiently small $\gamma_\rho$   (which appears in (<ref>) through $\varphi\equiv1-{\gamma_\rho\over\gamma}$), provided $\gamma$ and $c_0^2$ are both $>0$. Note that this condition is similar to the stability criterion for a simple compressible bulk fluid: $c_0^2>0$ simply means the bulk compressibility is positive, while $\gamma>0$ is analogous to requiring positive shear and bulk viscosities.  The stability condition on $\gamma_\rho$ is
-|ϖ-1|-√((ϖ-1)^2+1m^2)<γ_ργ<-|ϖ-1|+√((ϖ-1)^2+1m^2)  ,

where we have defined the “Mach number”
m≡v_ρc_0  ,

and the speed ratio
ϖ≡v_pv_ρ  .

It is easy to see that this condition can always be satisfied for sufficiently small $\gamma_\rho$; in particular, the allowed region always includes $\gamma_\rho=0$. That a sufficiently large $\gamma_\rho$ relative to the effective damping $\gamma$ can lead to instabilities is not surprising. A non-zero $\gamma_\rho$ (or, equivalently, $v_a'$) implies that different patches of the system move at different speeds. Of course, damping tends to homogenize  the density by exchanging particles,  thereby reducing the  speed  differences. However, if $\gamma_\rho$ is too large, damping may not be sufficient to suppress these speed  fluctuations.  The eventual steady state  may be nonuniform, leading to a patterned state, though its actual nature cannot be ascertained from the  linearized equations of motion. Nonlinear amplitude equations would be necessary for a full-fledged analysis of the steady state. It will be interesting to explore the relation between these instabilities here and the banding instabilities reported in  [64, 65].

generic underdamped propagating waves with anisotropic, $\theta_\bq$-dependent
wavespeed proportional to $q$ are expected for the wide range of parameters satisfying the stability condition derived in appendix B.

§ CORRELATION FUNCTIONS AND ROBUSTNESS OF LONG-RANGED ORDER AGAINST NOISE

In the stable region of the parameter space,
correlation functions for the system in the steady state may be calculated
from the noise-driven equations of motion. Dropping for now the density force ${\bf f}_\rho$ (we will show later that it is irrelevant in the long-wavelength limit), and solving the linear equations of motion for the spatio-temporally Fourier transformed fields $p_y(\bq, \omega)$ and $\delta\rho(\bq, \omega)$ gives
p_y(, ω)=i(ω-v_ρqcosθ_)f_y(,ω)(ω-c_+(θ_)q)(ω-c_-(θ_)q)+ i(ωΨ(1,θ_)-v_ρq Ψ(φ,θ_)cosθ_)q ,

δρ(, ω)=iρ_cv_ρqsinθ_f_y(,ω)(ω-c_+(θ_)q)(ω-c_-(θ_)q)+ i(ωΨ(1,θ_)-v_ρ q Ψ(φ,θ_)cosθ_)q ,
where we have defined
Ψ(φ, θ)≡γ(1+φsin^2θ) ,

and with the “sound speeds” $c_\pm(\theta_\bq)$, defined by the positions of the peaks in the scaling function versus $u$  [39],  precisely those found for dry active matter in [13, 14, 15]; i.e.,

\begin{eqnarray}
c_{\pm}\left(\theta_\bq \right) =
\left({v_\rho + v_p \over 2}\right)\cos\theta_\bq
\pm \sqrt{{1 \over 4}\left(v_\rho -v_p\right)^2 \cos^2
\theta_\bq + c^2_0
\sin^2 \theta_\bq} \quad ,
\label{cplusminus}
\end{eqnarray}

The `sound speeds” $c_\pm(\theta_\bq)$ in (<ref>), defined by the positions of the peaks in the scaling function versus $u$ [39],  are precisely those found for dry active matter in [13, 14, 15]. Further, these peak positions agree with those found in  the experiments of [29] on Quinke rotators.
In Fig. (<ref>)  we show a  polar plot of these sound speeds.

Autocorrelating these fields with themselves then gives their spatio-temporally Fourier transformed correlations:

C_pp(q,ω)≡⟨|p_y(, ω)|^2⟩=2D_p(ω-v_ρqcosθ_)^2(ω-c_+(θ_)q)^2(ω-c_-(θ_)q)^2+ (ωΨ(1,θ_)-v_ρ q
Ψ(φ,θ_)cosθ_)^2q^2  ,

C_ρρ(q,ω)≡⟨|δρ(, ω)|^2 ⟩=2D_pρ_c^2v_ρ^2 q^2sin^2θ_(ω-c_+(θ_)q)^2(ω-c_-(θ_)q)^2+ (ωΨ(1,θ_)-v_ρ q Ψ(φ,θ_)cosθ_)^2q^2 ,

where we have used (<ref>) to obtain the autocorrelation $\langle|f_y(\bq,\omega)|^2\rangle=2D_p$.

Similar reasoning gives the cross-correlation function
C_pρ(q, ω)≡⟨ρ(q,ω) p_y(-q,-ω)⟩= 2D_p ρ_c v_ρqsinθ_(ω- v_ρqcosθ_)/(ω-c_+(θ_)q)^2(ω-c_-(θ_)q)^2+(ωΨ(1,θ_)-v_ρq Ψ(φ,θ_)cosθ_)^2q^2 .

Pulling a factor of $q^2$ out of the numerator of each of these expressions, and a factor of $q^4$ out of their denominators, gives the scaling forms  (<ref>), (<ref>),  and
(<ref>), with the scaling functions $F_{pp}(u_f)$, $F_{\rho\rho}(u_f)$ and $F_{p\rho}(u_f)$.

The scaling function $F_{pp}(u_f, \theta_\bq)$ is given by
F_pp(u_f, θ_)=2D_p(u_f-v_ρcosθ_)^2(u_f-c_+(θ_))^2(u_f-c_-(θ_))^2+(u_fΨ(1,θ_q)-v_ρΨ(φ,θ_q)cosθ_)^2 ,

where  $u_f\equiv \omega/q$.

Note that the existence of the scaling form (<ref>) implies that the ratio of the widths of the peaks
in  $C_{pp}(\bq, \omega)$, plotted versus $\omega$ for fixed $\bq$, to their positions does not change as $\bq\to{\bf 0}$; this is what we meant by our earlier cryptic comment that the “quality factor $Q$” becomes independent of $q$.

In Fig. (<ref>),  we show a  polar plot of the sound speeds $c_\pm(\theta_\bq)$, which are the positions of the peaks in the scaling function $F_{pp}$ when plotted versus the scaling argument $u_f$ for fixed direction of propagation $\theta_\bq$. As such, they are the analog in our system of the sound speeds in a simple compressible fluid.

(Color online) Polar plot of the sound speeds (<ref>) (black curve). Here we have taken $q=1$, $v_\rho=1$, $v_p=c_0=2$, and
$\gamma=.3$ (all in arbitrary units).The polarization points directly  to the right. The   red straight line indicates the direction of $\bq$ for which the scaling function $F_{pp}$ is plotted in figure (<ref>).

The density-density correlation function obeys a similar scaling law:
C_ρρ(q,ω)≡⟨|p_y(, ω)|^2⟩=1q^2F_ρρ((ωq), θ_)  ,

with a slightly different scaling function $F_{\rho\rho}(u_f, \theta_\bq)$  given by

F_ρρ(u_f, θ_)=2D_p  ρ_c^2v_ρ^2 sin^2θ_(u_f-c_+(θ_))^2(u_f-c_-(θ_))^2+(u_fΨ(1,θ_q)-v_ρΨ(φ,θ_q)cosθ_)^2  .

The orientation-density correlation function also obeys a similar scaling law:
C_pρ(q,ω)≡⟨ρ(q,ω) p_y(-q,-ω)  ⟩=1q^2F_pρ((ωq), θ_)  ,

where the scaling function for this correlation function is
F_pρ(u_f, θ_)
= 2D_p ρ_c v_ρsinθ_(u_f - v_ρcosθ_)/(u_f-c_+(θ_))^2(u_f-c_-(θ_))^2+(u_fΨ(1,θ_)-v_ρΨ(φ,θ_)cosθ_)^2 .

These three scaling functions are plotted versus the scaling argument $u$ for a fixed direction of propagation $\theta_\bq$ in Figure (<ref>). Note that the two peaks in the polarization scaling function $F_{pp}$ are exactly the same height, but have different widths. We have chosen to plot this figure for a fairly generic direction of propagation. In contrast, for $\theta=\pi/2$ (that is, propagation perpendicular to the polarization, which corresponds to a vertical line in figure (<ref>)), the scaling functions $F_{pp}$ and $F_{\rho\rho}$ both become even functions of the scaling argument $u$,  with two symmetrically placed peaks at $\pm c_0$. This implies that  the correlation functions both become even functions of $\omega$, with two symmetrically placed peaks at $\pm c_0q$.
Furthermore, for this direction of propagation, the two correlations functions are precisely proportional to each other; that is, their ratio is a constant, independent of both $q$ and $\omega$.

Another special direction is $\theta_\bq=0$; i.e., propagation along the polarization. In this case, the asymmetry is maximized: the polarization scaling function has only one peak, at scaling argument $u=v_p$, meaning the correlation function has a single peak at $\omega=v_pq$. The density-density scaling function  on the other hand, vanishes identically in this limit (except at the singular point $\omega=v_\rho q$, where its value depends on precisely how the limit $\omega\to v_\rho q$, $\theta_\bq$ is approached).

Both of these special directions are, indeed, special: for generic directions of propagation $\theta\ne0,\pm \pi/2$, the scaling functions, and, hence, the $\omega$-dependence of the correlation functions, look like figure (<ref>). That is, each scaling function has two in general asymmetrically placed peaks, in the same positions for both correlation functions.

The way the special limit $\theta=0$ is approached as $\theta\to0$ is that the first  (i.e., leftmost) peak of the polarization scaling function $F_{pp}$, while its height stays the same as that of the second (i.e., rightmost) peak, has its width continuously vanish as $\theta\to0$, while the density scaling function $F_{\rho\rho}$ has  the height of its first peak continuously vanish as $\theta\to0$. One can see the tendency towards this limit in the plot of figure (<ref>) at $\theta=27^\circ$.

(Color online) Plot of the scaling functions $F_{pp}(u_f,\theta_\bq)$  ( red curve), $F_{\rho\rho}((u_f,\theta_\bq)$ ( black curve), and  $F_{p\rho}(u_f,\theta_\bq)$ (green curve) versus scaling argument $u_f\equiv{\omega\over q}$  for fixed $\theta_\bq$. Here we have taken $v_\rho=1$, $v_p=c_0=2$,
$\gamma=.3$, $D_p=.5$,  $\rho_c=1$, ${\gamma_\rho/\gamma}={1/2}$  (all in arbitrary units), and $\theta=27^\circ$ (the direction indicated by the straight line in figure (<ref>).

We can use these expressions to calculate the equal time correlation functions; we find for the polarization $p_y$

\begin{eqnarray}
&&C_{pp}(\bq)\equiv\langle |p_y(\bq,t)|^2\rangle=\int_{-\infty}^\infty
\frac{d\omega}{2\pi}\langle |p_y(\bq, \omega)|^2\rangle
={D_p\over\pi} \int_{-\infty}^\infty
{(\omega-v_\rho q\cos\theta_\bq)^2d\omega\over(\omega-c_+(\theta_\bq)q)^2(\omega-c_-(\theta_\bq)q)^2+ (\omega\Psi(1,\theta_\bq)-v_\rho { q} \Psi(\varphi,\theta_\bq)\cos\theta_\bq)^2q^2}\nn\\
\label{pet1}
\end{eqnarray}

Introducing a new variable of integration
s≡ω-v_ρqcosθ_   Ξ(φ, θ_) Ψ(1,θ_)q

where we have defined
Ξ(Φ, θ)≡Ψ(Φ, θ)Ψ(1, θ)=1+Φsin^2θ1+sin^2θ

\begin{eqnarray}
&&C_{pp}(\bq)={D_p\over\pi{ \Psi(1, \theta_\bq)}{ q}} \int_{-\infty}^\infty
{{ (s+\psi(\varphi, \theta_\bq)}^2ds\over(s-\psi_+( \varphi,  \theta_\bq))^2(s-\psi_-({ \varphi}, \theta_\bq))^2+s^2} \,,
\label{pet2}
\end{eqnarray}

where we have defined
ψ_±(φ, θ_)≡c_±(θ_)-v_ρcosθ_   Ξ(φ, θ_)Ψ(1, θ_)  ,

ψ(φ, θ_)≡v_ρcosθ_ ( Ξ(φ, θ_)-1)Ψ(1,θ_)  .

We will show in appendix (<ref>)
that,  whenever the stability condition is satisfied, $\psi_+>0$ and $\psi_-<0$; we will make use of these facts later.

We also show in appendix <ref> that this integral is equal to

∫_-∞^∞(s+ψ(φ, θ_)^2ds(s-ψ_+(φ,  θ_))^2(s-ψ_-( φ, θ_))^2+s^2 = π[1+(γ_ργ)^2m^2sin^2θ_cos^2θ_/{(1+sin^2θ_)^2-m^2γ_ργ[ϖ-1+(ϖ-φ)sin^2θ_]
cos^2θ_}]  .

Using (<ref>), we obtain the expression for the equal-time correlation function:

\begin{align}
C_{pp}(\bq) &={D_p\over q\Psi(1, \theta_\bq)}\left[1+\left({\gamma_\rho\over\gamma}\right)^2\frac{m^2\sin^2\theta_\bq\cos^2\theta_\bq}{D(m^2,\varpi,\varphi,\theta_\bq)}\right] \nonumber \\
&={D_p\over\gamma q(1+\sin^2\theta_\bq)}\left[1+\left({\gamma_\rho\over\gamma}\right)^2\frac{m^2\sin^2\theta_\bq\cos^2\theta_\bq}{D(m^2,\varpi,\varphi,\theta_\bq)}\right]  \,,
\label{petfinal}
\end{align}

where we have defined

This expression is easily Fourier transformed in space to give the real space correlations $C_{pp}(\br)$ of the polarization fluctuations:
C_pp()≡⟨p_y(+,t)p_y(,t)⟩=∫d^2q(2π)^2 e^i·
C_pp()   .

Defining $\theta_\br$ as the angle between $\br$ and the $x$-axis, we can rewrite this as

C_pp() =D_p(2π)^2γ ∫_0^2π dθ_1+ sin^2θ_[1+(γ_ργ)^2m^2sin^2θ_cos^2θ_/D(m^2,ϖ,φ,θ_)]∫_0^∞dq  e^iqrcos(θ_-θ_)  .

The integral over $q$ is given by
∫_0^∞dq  e^iqrcos(θ_-θ_)=πδ(rcos(θ_-θ_))+𝒫[1rcos(θ_-θ_)]

where $\mathcal{P}$ denotes the Cauchy principal value.
Using this in (<ref>), we have

C_pp()=D_p2πγ ∫_0^2π dθ_1+ sin^2θ_[1+(γ_ργ)^2m^2sin^2θ_cos^2θ_/D(m^2,ϖ,φ,θ_)](δ(rcos(θ_-θ_))+𝒫[1πrcos(θ_-θ_)])  .

The Cauchy principal value term in this expression is odd under the operation $\theta_\bq\to\theta_\bq+\pi$, while ${1\over1+ \sin^2\theta_\bq}$ is even under this operation; hence, the contribution to integral from the Cauchy principal value term vanishes. We are therefore left with

C_pp()=D_p4πγ ∫_0^2π dθ_1+ sin^2θ_[1+(γ_ργ)^2m^2sin^2θ_cos^2θ_/D(m^2,ϖ,φ,θ_)]
δ(rcos(θ_-θ_))  .

δ(rcos(θ_-θ_))=δ(θ_-θ_+π/2)+δ(θ_-θ_-π/2)r|(cos(θ_-θ_)θ_)_θ_=θ_±π/2|=δ(θ_-θ_+π/2)+δ(θ_-θ_-π/2)r  ,

where the first equality follows from the familiar identity
δ(f(x))=∑_x_0δ(x-x_0)/|f'(x_0)|  ,
where $\{x_0\}$ are the roots of $f(x)$,
(<ref>) becomes

C_pp()    =    D_p4πγr ∫_0^2π dθ_1+ sin^2θ_[1+(γ_ργ)^2m^2sin^2θ_cos^2θ_/D(m^2,ϖ,φ,θ_)]

=D_pr2πγ(r^2+ x^2)[1+(γ_ργ)^2m^2y^2 x^2/{(r^2+x^2)^2-m^2γ_ργy^2[r^2(ϖ-1)+(ϖ-φ)x^2]}]  .

We now turn to the equal-time density correlations. Using (<ref>), we have

\begin{eqnarray}
&&C_{\rho\rho}(\bq)\equiv\langle |\rho(\bq,t)|^2\rangle\nn\\&&=\int_{-\infty}^\infty
\frac{d\omega}{2\pi}\langle |\rho(\bq, \omega)|^2\rangle
={D_p\rho_c^2v_\rho^2q^2\sin^2\theta_\bq\over\pi} \int_{-\infty}^\infty
{d\omega\over(\omega-c_+(\theta_\bq)q)^2(\omega-c_-(\theta_\bq)q)^2+(\omega\Psi(1,\theta_\bq)-v_\rho q \Psi(\varphi,\theta_\bq)\cos\theta_\bq)^2q^2} \,.\nn\\
\label{rhoet1}
\end{eqnarray}

Making the same change of variables of integration (<ref>) that we made earlier, and using partial fractions again, we find

\begin{eqnarray}
C_{\rho\rho}(\bq)={D_p\rho_c^2v_\rho^2q^2\sin^2\theta_\bq\over2i\pi[\gamma q(1+ \sin^2\theta_\bq)]^3} (I_\rho^--I_\rho^+) \,,
\label{rhoet2}
\end{eqnarray}

where we have defined

\begin{eqnarray}
&&I_\rho^\pm\equiv \lim_{\epsilon\to0} \int_{-\infty}^\infty
{ds\over(s-i\epsilon-\psi_+(\theta_\bq))(s-i\epsilon-\psi_-(\theta_\bq))\pm i(s-i\epsilon)}\frac{1}{s-i \epsilon} \,,
\label{irhopm}
\end{eqnarray}

where we have introduced the small parameter $\epsilon$ to regularize the integral.

We can now do both integrals by the usual complex contour techniques. Note that, unlike the integral for the polarization correlation function $C_{pp}$, the integrand in (<ref>) converges rapidly enough at infinity that the infinite semicircle needed to close the path of integration contributes nothing to the integral. Note also that the poles of $I^\pm_\rho$ are those of $I_\pm$ considered earlier, plus one more pole at $s=i\epsilon$.

If we choose $\epsilon>0$, then we can close the integral for $I_\rho^+$ in the upper half plane, picking up the pole at $s=i\epsilon$. This gives
I_ρ^+=2iπψ_+ψ_-  ,

where we have taken the limit $\epsilon\to0$.

On the other hand, we can close the contour for $I_\rho^-$ in the lower half plane, finding (for $\epsilon>0$) no poles at all. This implies
I_ρ^-=0  .

Inserting (<ref>) and (<ref>) into (<ref>), we obtain

\begin{eqnarray}
C_{\rho\rho}(\bq)=-{D_p\rho_c^2v_\rho^2q^2\sin^2\theta_\bq\over[\gamma q(1+ \sin^2\theta_\bq)]^3\psi_+\psi_-}  \,.
\label{rhoet3}
\end{eqnarray}

It is straightforward to check that, had we chosen $\epsilon<0$, we would have gotten the same answer, although now with $I_\rho^+=0$ and $I_\rho^-=-{2i\pi\over\psi_+\psi_-}$.

It is an equally straightforward algebraic exercise, using the definitions (<ref>) of $\psi_\pm$ and our expressions (<ref>)  for $c_\pm(\theta_\bq)$ to show that

ψ_+ψ_-=-D(m^2,ϖ,φ,θ_) v_ρ^2 γ^2sin^2θ_/m^2
[Ψ(1,θ_)]^4  .

Using this in (<ref>) gives

\begin{eqnarray}
C_{\rho\rho}(\bq)={D_p\rho_c^2 m^2 (1+ \sin^2\theta_\bq)\over \gamma q D(m^2,\varpi,\varphi,\theta_\bq)} \,.
\nonumber\\
\label{rhoet4}
\end{eqnarray}

The above expression can be easily Fourier transformed in space to obtain the real space density
correlations $C_{\rho\rho}(\br)$ just as we did above to obtain $C_{pp}(\br)$.  We thereby obtain:
D_p ρ_c^2m^2/2πγ r(r^2+x^2)/{(r^2+x^2)^2-m^2γ_ργy^2[r^2(ϖ-1)+(ϖ-φ)x^2]}  .

Finally we work out the expression for the crosscorrelation function

\begin{eqnarray}
C_{p\rho}({\bf q})\equiv \langle \rho({\bf q})p_y(-{\bf q})\rangle&=&\int_{-\infty}^{\infty} \frac{d\omega}{2\pi} \langle \rho({\bf q},\omega) p_y(-{\bf q},-\omega)\rangle \nonumber \\
&=& \int_{-\infty}^{\infty} \frac{d\omega}{2\pi} \frac{2D_p \rho_c v_\rho q\sin\theta_\bq (\omega - v_\rho q\cos\theta_\bq)} {(\omega-c_+(\theta_\bq)q)^2(\omega-c_-(\theta_\bq)q)^2+(\omega\Psi(1,\theta_\bq)-v_\rho q \Psi(\varphi,\theta_\bq)\cos\theta_\bq)^2q^2}\nonumber \\
&=&\frac{D_p\rho_cv_\rho q\sin\theta_\bq} {\pi\gamma^2 q^2 (1+ \sin^2\theta_\bq)^2}\int_{-\infty}^{\infty}\frac{(s+ \psi(\varphi,\theta_\bq))ds}{(s-\psi_+)^2(s -\psi_-)^2 +s^2}.
\label{crossprho}
\end{eqnarray}

From equations (<ref>) and (<ref>), we can easily  show that

\begin{eqnarray}
\int_{-\infty}^{\infty}\frac{(s+ \psi(\varphi,\theta_\bq))ds}{(s-\psi_+)^2(s -\psi_-)^2 +s^2} = I_2 +
\psi(\varphi,\theta_\bq) I_3
\end{eqnarray}

\begin{equation}
I_2=0, \; \mbox{and} \,\, I_3= -\frac{\pi}{\psi_+\psi_-}.
\label{i2i3}
\end{equation}

Using eqs. (<ref>) and (<ref>), we obtain

C_pρ()= - D_pρ_c sinθ_cosθ_m^2 γ_ρ/γ^2 q D(m^2,ϖ,φ,θ_)  .

Fourier transformation of the above expression gives the real space cross correlation, given by

C_pρ() = - D_pρ_c m^2 γ_ρ/2πγ^2 xyr/{(r^2+x^2)^2-m^2γ_ργy^2[r^2(ϖ-1)+(ϖ-φ)x^2]}.

To summarize, we have shown that  over a finite range of the model parameters, the linearized model displays long range orientational order. We argue below that this linear theory prediction is asymptotically exact, in contrast to the well-known Toner-Tu model for dry active matter,  for which nonlinear effects are essential for the existence of the long range order [13, 14].

§ GIANT NUMBER FLUCTUATIONS AND THEIR SHAPE DEPENDENCE

We will show in this section that the long-ranged density correlations found in equation (<ref>)   of the previous section lead to “giant number fluctuations”[66][15, 19, 17]. These can be defined as follows:

Consider a “counting box”, defined as
a rectangular area $A=L_x\times L_y$, and define the aspect ratio $\alpha_A=L_x/L_y$. Our experiment will consist of counting the number of active particles in this box. The mean number of particles $\overline N$ in the box is, of course, given by

\begin{equation}
\overline N= L_x L_y\rho_0 = \alpha_A L_y^2 \rho_0 \,.
\end{equation}

The total number fluctuation $\delta N\equiv N-\overline N$ in the area $A$ is

\begin{equation}
\delta N = \int d^2 r\,\delta\rho(\br,t).
\end{equation}

\begin{eqnarray}
&&\langle (\delta N)^2\rangle = \int_0^{L_x} dx\int_0^{L_y} dy \int_0^{L_x}dx'\int_0^{L_y} dy' \langle \rho ({\bf r},t)\rho({\bf r'},t)\rangle = \int_0^{L_x} dx\int_0^{L_y} dy \int_0^{L_x}dx'\int_0^{L_y} dy' C_{\rho\rho}({\bf r}- {\bf r'})\nonumber \\
&&=\frac{D_p\rho_c^2 m^2}{2\pi \gamma } \int_0^{L_x} dx\int_0^{L_y} dy \int_0^{L_y} dy'\int_0^{L_x}dx'\,\frac{|{\bf r-r'}|}{({\bf r-r'})^2 + (x-x')^2}\nonumber \\
&&\times \left[1-\frac{(v_p-v_\rho)(y-y')^2}{({\bf r-r'})^2 + (x-x')^2}\frac{\gamma_\rho v_\rho}{\gamma c_0^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2(x-x')^2} {\{({\bf r-r'})^2 + (x-x')^2\}^2}\right]^{-1}.
%\frac{|{\bf r-r'}|[(r-r')^2 + (x-x')^2]}{D_r(m^2,\varpi,\varphi,\br-\br')}.
\label{Nvar1}
\end{eqnarray}

We consider the limits $\alpha_A\gg 1$ and $\alpha_A \ll 1$.

Consider first the limit $L_x\gg L_y$, i.e., $\alpha_A \gg 1$. Now consider the integral over $x'$  for fixed ${\bf r}=(x,y)$ in this limit. Note first that in this limit, for most of the range of integration over $x$, $x\gg L_y$. For such values of $x$, we can split  the integral over $x'$  into  three parts:

\begin{eqnarray}
&&\int_0^{L_x}dx'\frac{|{\bf r-r'}|}{({\bf r-r'})^2 + (x-x')^2}\nonumber \\
&&\times \left[1-\frac{(v_p-v_\rho)(y-y')^2}{({\bf r-r'})^2 + (x-x')^2}\frac{\gamma_\rho v_\rho}{\gamma c_0^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2(x-x')^2} {\{({\bf r-r'})^2 + (x-x')^2\}^2}\right]^{-1}=I_1+I_2+I_3 \,,
\label{split1}
\end{eqnarray}

where we have defined

\begin{eqnarray}
I_1&\equiv&\int_0^{x-CL_y}dx'\frac{|{\bf r-r'}|}{({\bf r-r'})^2 + (x-x')^2}\nonumber \\
&&\times \left[1-\frac{(v_p-v_\rho)(y-y')^2}{({\bf r-r'})^2 + (x-x')^2}\frac{\gamma_\rho v_\rho}{\gamma c_0^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2(x-x')^2} {\{({\bf r-r'})^2 + (x-x')^2\}^2}\right]^{-1} \,,
\label{I1def}
\end{eqnarray}

\begin{eqnarray}
I_2&\equiv&\int_{x-CL_y}^{x+CL_y}dx'\frac{|{\bf r-r'}|}{({\bf r-r'})^2 + (x-x')^2}\nonumber \\
&&\times \left[1-\frac{(v_p-v_\rho)(y-y')^2}{({\bf r-r'})^2 + (x-x')^2}\frac{\gamma_\rho v_\rho}{\gamma c_0^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2(x-x')^2} {\{({\bf r-r'})^2 + (x-x')^2\}^2}\right]^{-1} \,,
\label{I2def}
\end{eqnarray}

\begin{eqnarray}
I_3&\equiv&\int_{x+CL_y}^{L_x}dx'\frac{|{\bf r-r'}|}{({\bf r-r'})^2 + (x-x')^2}\nonumber \\
&&\times \left[1-\frac{(v_p-v_\rho)(y-y')^2}{({\bf r-r'})^2 + (x-x')^2}\frac{\gamma_\rho v_\rho}{\gamma c_0^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2(x-x')^2} {\{({\bf r-r'})^2 + (x-x')^2\}^2}\right]^{-1} \,,
\label{I3def}
\end{eqnarray}

where $C$ is an arbitrary constant enough larger than $1$ that we can safely neglect terms of ${\cal O}({1\over C})$. With this choice, throughout the region of integration for $I_1$, $x-x'\gg L_y$. Therefore, throughout this region, $|y-y'|\ll x-x'$ (since $|y-y'|$ is always $<L_y$). This gives $\left(\frac{y-y'}{x-x'}\right)^2 \ll 1$.
Making this approximation for $I_1$, our expression for it reduces to
I_1≈12∫_0^x-CL_ydx'x-x'  .

Performing this elementary integral gives
I_1=12 ln(xCL_y)=12 ln(xL_y) -12 lnC  .

Virtually identical reasoning can be applied to $I_3$, giving the result
I_3=12 ln(L_x-x CL_y)  .

For $I_2$, a simple shift of variables of integration $x'=x+x''$ shows that $I_2$ is independent of $x$:

\begin{eqnarray}
I_2&=&\int_{-CL_y}^{CL_y}dx''\frac{\sqrt{x''^2 + (y-y')^2}}{(y-y')^2 + 2x''^2}\nonumber \\
&&\times \left[1-\frac{\gamma_\rho v_\rho}{\gamma c_0^2}\frac{(v_p-v_\rho)(y-y')^2} {(y-y')^2 + 2x''^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2{ x''}^2} {[(y-y')^2 + 2x''^2}]^2\right]^{-1}\,.
\label{I2result1}
\end{eqnarray}

Inserting these results (<ref>),  (<ref>), and  (<ref>) into our expression (<ref>), and using that result in our expression (<ref>) for $\langle (\delta N)^2\rangle$, we obtain

\begin{eqnarray}
&&\langle (\delta N)^2\rangle  %\frac{D_p\rho_c^2 m^2}{2\pi \gamma }\int_0^{L_x} dx'\frac{|{\bf r-r'}|}{({\bf r-r'})^2 + (x-x')^2}
%\nonumber\\ &&\times \left[1-\frac{(v_p-v_\rho)(y-y')^2}{({\bf r-r'})^2 + (x-x')^2}\frac{\gamma_\rho v_p}{\gamma c_0^2} - \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2(x-x')^2} {\{({\bf r-r'})^2 + (x-x')^2\}^2}\right]^{-1}\nonumber \\&&
\approx \frac{D_p\rho_c^2 m^2}{2\pi \gamma }\Bigg\{{1\over2 }\int_0^{L_x} dx\int_0^{L_y} dy \int_0^{L_y}dy'\left(\ln\left({x\over L_y}\right)+\ln\left({L_x-x\over L_y}\right) -2\ln C \right)\nonumber\\&&
+\int_0^{L_x} dx\int_0^{L_y} dy \int_0^{L_y}dy'\int_{-CL_y}^{CL_y}dx''\frac{ \sqrt{x''^2+(y-y')^2}}{ [(y-y')^2 +2x''^2]} %\nonumber \\
\times \left[1-\frac{\gamma_\rho v_\rho}{\gamma c_0^2}\frac{(v_p-v_\rho)(y-y')^2}
{ { [(y-y')^2 +2x''^2]} }
- \left(\frac{v_\rho\gamma_\rho}{c_0\gamma}\right)^2 \frac{(y-y')^2{ x''^2}} { [(y-y')^2 +2x''^2]^2}\right]^{-1}\Bigg\} \,,\nonumber\\
\label{delNfinal1}
\end{eqnarray}

The integrals over $y$ and $y'$ in the (triple) integral in the first line of this expression trivially give a factor of $L_y^2$,
since the integrand is independent of $y$ and $y'$. The remaining integral over $x$ is elementary. The net result is
12 ∫_0^L_x dx∫_0^L_y dy ∫_0^L_ydy'(ln(xL_y)+ln(L_x-xL_y) -2lnC )= L_xL_y^2(ln(L_xCL_y)-1) .

The $x$ integral in the remaining (quadruple) integral which appears on the second line in equation (<ref>) can be done immediately, since the integrand is independent of $x$, yielding a factor of $L_x$. The remaining triple integral over $y$, $y'$, and $x''$ can be done by changing variables of integration to new, rescaled variables $u_x$, $u_y$, and $u'_y$ via

x”≡L_yu_x  .

This gives
∫_0^L_x dx∫_0^L_y dy ∫_0^L_ydy'∫_-CL_y^CL_ydx”       √(x”^2+(y-y')^2)/ [(y-y')^2 +2x”^2]×[1-γ_ρv_ρ/γc_0^2(v_p-v_ρ)(y-y')^2/  [(y-y')^2 +2x”^2]
- (v_ργ_ρ/c_0γ)^2 (y-y')^2 x”^2/ [(y-y')^2 +2x”^2]^2]^-1

=C' L_xL_y^2  ,

C'   ≡    ∫_0^1 du_y∫_0^1 du_y'∫_-C^C du_x  √(u_x^2 + (u_y-u_y')^2)/(u_y - u_y')^2 + 2 u_x^2

[1-γ_ρv_ρ/γc_0^2(v_p-v_ρ)(u_y-u_y')^2/2u_x^2 + (u_y - u_y')^2 - (v_ργ_ρ/c_0γ)^2 (u_y-u_y')^2 u_x^2/[2u_x^2 + (u_y - u_y')^2]^2]^-1

is an ${\cal O}(1)$ constant. Comparing  (<ref>) and (<ref>), we see that the first line of
(<ref>)  actually dominates the second in the large aspect ratio limit $L_x\gg L_y$ that we are considering here. Therefore, we obtain, in the limit of large aspect ratio ($\alpha_A\gg1$):

\begin{equation}
\sqrt{\langle (\delta N)^2\rangle} \approx \frac{\rho_c m\sqrt{D_p}}{\sqrt{2\pi \gamma}}
\sqrt{L_xL_y^2\ln\left({L_x\over L_y}\right)}\approx\frac{\rho_c m\sqrt{D_p}}{\sqrt{2\pi \gamma}} L_y^{3/2}\sqrt{\alpha_A \ln \alpha_A}\,,
\label{Nvarlar}
\end{equation}

This expression (<ref>) can be rewritten in terms of
the mean particle number $\overline N$ in the same area $A$,  which is given by

\begin{equation}
\overline N= L_x L_y\rho_0 = \alpha_A L_y^2 \rho_0 \,.
\end{equation}

This gives
L_y=(N/(α_Aρ_0))^1/2  .

Using this result (<ref>) in our expression (<ref>)  for the variance $(\Delta N)^2$ of the number fluctuations gives

\begin{equation}
\Delta N= \sqrt{\langle(\delta N)^2\rangle} = \left({\rho_c m\over \rho_0^{3/4}}\right)\left(\sqrt {D_p\over
2\pi \gamma}\right) \left(\frac{{\overline N}^{3/4}\sqrt{\ln{\alpha_A}}}{\alpha_A^{1/4}}\right)    \,.
\label{bigalpha}
\end{equation}

Three points should be noted about this result:

1) The number fluctuations are giant; that is, they grow much more rapidly with $\bar N$ than the usual $\sqrt{\bar N}$ “law of large numbers” fluctuations, which are found in almost all systems, and, in particular, in most equilibrium systems away from fixed points[71].  Specifically, we have $\Delta N \propto  \overline N^{3/4}$. Note also  that the exponent $3/4$ is very close to the exponent $.8$ found in the experiments of  [27].

2) The size of the number fluctuations depend not only on the mean number $\bar N$, but also on the shape (i.e., on the aspect ratio $\alpha_A$).

3) The number fluctuations $\Delta N$ are a monotonically decreasing function of the aspect ratio $\alpha_A$ in this range of $\alpha_A\gg1$.

We analyze giant number fluctuations in the opposite regime of $L_y  \gg L_x$, i.e.,  small aspect ratio $\alpha_A\ll 1$ in Appendix <ref>.
As shown there, in this regime of small aspect ratio $\alpha_A$, $\Delta N$ is a monotonically increasing function of $\alpha_A$. Recall that in the opposite limit of  $\alpha_A\gg 1$, eqn (<ref>)    , we found that $\Delta N$ is a monotonically decreasing function of $\alpha_A$. Hence, the maximum value of $\Delta N$ for a given mean number of particles $\overline N$ will occur when $\alpha_A\sim 1$; i.e., for a roughly square counting box.

In contrast to our above results on giant number fluctuations, Ref. [28] did not report any giant number fluctuations in their experiment. We believe this is due to the fact that in their experiment, the density fluctuations were probed at length scales larger than the height of the passive fluid over the active layer, which is outside the regime of validity of our theory.

§ BULK VELOCITY FLUCTUATIONS

We can use the relations (<ref>) and (<ref>)  between the bulk fluid velocity field and the surface velocity, and our boundary condition (<ref>) for that surface velocity, to obtain expressions for the bulk velocity correlations  in terms of the polarization $\hp$ and density $\rho$ correlation functions. This gives, to the lowest order in gradients

⟨v_i(,z)v_j(, z)⟩=v_0^2 ∫d^2q(2π)^2 C_pp()e^-2qz[δ_iyδ_jy-q_yzq(δ_iyq_j+δ_jyq_i)+q_y^2z^2(q_iq_jq^2+δ_izδ_iz)]

+ v_0v_a' ∫d^2q(2π)^2 C_pρ()e^-2qz[δ_iyδ_jx+
δ_ixδ_jy - z(q_yqδ_ixq_j+ q_yqδ_jxq_i + q_xqδ_jyq_i
+ q_xqδ_iyq_j)+ 2q_xq_yz^2(q_iq_jq^2+δ_izδ_iz)]

+ v_a'^2 ∫d^2q(2π)^2 C_ρρ()e^-2qz[δ_ixδ_jx-q_xzq(δ_ixq_j+δ_jxq_i)+q_x^2z^2(q_iq_jq^2+δ_izδ_iz)]

The second ($\delta_{iy}q_j$) term in the first integral is odd in $q_x$ if $j=x$. Since $C_{pp}(\bq)e^{-2qz}$ and $q_y$ are even in $q_x$, the integral of this term vanishes when $i=x$. Hence, we can replace this term by $\delta_{iy}\delta_{jy}q_y$. Similar arguments imply that we can also replace the third ($\delta_{jy}q_i$) term with $\delta_{iy}\delta_{jy}q_y$. Similarly we can replace ($\delta_{ix}q_j$) with
$\delta_{ix}\delta_{jx}q_x$, and ($\delta_{jx}q_i$) with $\delta_{ix}\delta_{jx}q_x$ in the third integral. In the second integral ($\delta_{ix}q_j$) can be replaced with $\delta_{ix}\delta_{jy}q_y$, ($\delta_{jx}q_i$)
with $\delta_{jx}\delta_{iy}q_y$, ($\delta_{jy}q_i$) with $\delta_{ix}\delta_{jy}q_x$, and ($\delta_{iy}q_j$)
with $\delta_{iy}\delta_{jx}q_x$ respectively. Thus, we can rewrite (<ref>) as
⟨v_i(,z)v_j(, z)⟩   =   v_0^2 ∫d^2q(2π)^2 C_pp()e^-2qz[δ_iyδ_jy(1-2q_y^2zq)+q_y^2z^2(q_iq_jq^2+δ_izδ_jz)]

+ v_0v_a' ∫d^2q(2π)^2 C_pρ()e^-2qz[(δ_iyδ_jx+
δ_ixδ_jy)(1 - zq) + 2z^2q_xq_y(q_iq_jq^2+δ_izδ_jz)]

+ v_a'^2 ∫d^2q(2π)^2 C_ρρ()e^-2qz[δ_ixδ_jx(1-2q_x^2zq)+q_x^2z^2(q_iq_jq^2+δ_izδ_jz)].

We also note that the $q_iq_j\over q^2$ term in the first integral is also odd in $q_x$ unless $i=j$, because, if $i\ne j$, then one, and only one, of the indices $(i,j)$ must be $x$. Hence, its integral also vanishes if $i\ne j$.  Identical reasoning implies that the third integral  also vanishes for $i \ne j$. This makes the 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, since that term is proportional to $\delta_{iy}\delta_{jx}+
\delta_{ix}\delta_{jy}$. To see that it is also true for the second term (i.e., the $z^2q_xq_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 integrate to zero otherwise.

In light of these observations, we can rewrite our expression (<ref>) for the velocity correlations as
⟨v_i(,z)v_j(, z)⟩= v_0^2 ^(1)_ij(z)+v_0v_a'^(2)_ij(z)+v_a'^2^(3)_ij(z)   ,
where we have defined
^(1)_ij(z)≡∫d^2q(2π)^2 C_pp()e^-2qz[δ_iyδ_jy(1-2q_y^2zq)+q_y^2z^2(q_iq_jq^2+δ_izδ_jz)]

^(2)_ij(z)≡∫d^2q(2π)^2 C_pρ()e^-2qz[1 - zq + 2z^2q_x^2q_y^2q^2](δ_iyδ_jx+

^(3)_ij(z)≡∫d^2q(2π)^2 C_ρρ()e^-2qz[δ_ixδ_jx(1-2q_x^2zq)+q_x^2z^2(q_iq_jq^2+δ_izδ_jz)].

Using the fact that all three of the correlation functions $C_{\rho\rho}(\bq)$, $C_{p\rho}(\bq)$, and $C_{pp}(\bq)$ are proportional to ${1\over q}$ times functions that depend only on the direction $\theta_\bq$ of $\bq$, and making the  simple change of variables $\bq={\bQ\over z}$, shows that all three of these (tensor)  integrals $\varPi^{(1,2,3)}(z)$
are proportional to ${1\over z}$.

Thus, we conclude that
⟨v_i(,z)v_j(, z)⟩= (1z) ,

where we have defined the constant, parameter dependent matrix
≡v_0^2 ^(1)_ij+v_0v_a'^(2)_ij+v_a'^2^(3)_ij   ,
^(1)_ij(z)≡∫d^2Q(2π)^2 f_pp(θ_)Qe^-2Q[δ_iyδ_jy(1-2Q_y^2Q)+Q_y^2(Q_iQ_jQ^2+δ_izδ_jz)],

^(2)_ij(z)≡∫d^2Q(2π)^2 f_pρ(θ_)Qe^-2Q[1 - Q + 2Q_x^2Q_y^2Q^2](δ_iyδ_jx+

^(3)_ij(z)≡∫d^2Q(2π)^2 f_ρρ(θ_)Qe^-2Q[δ_ixδ_jx(1-2Q_x^2Q)+Q_x^2(Q_iQ_jQ^2+δ_izδ_jz)] .

In light of the above discussion, $\bM$ is a symmetric matrix, whose only non-zero off-diagonal components are $M_{xy}=M_{yx}$. Its diagonal entries are all, in general, different from each other, and from $M_{xy}$.

The important point about these velocity correlations is that they scale like ${1\over z}$; that is, inversely proportional to the distance $z$ from the solid surface. This ${1\over z}$ scaling will break down once $z$ becomes microscopic, as can be seen as follows:
our arguments above depended on our hydrodynamic theory, which breaks down for wavevectors $q$ comparable to an inverse microscopic length. Since the integrals over wavevector that we have done to derive (<ref>) were dominated by $q\sim{1\over z}$, the calculation clearly ceases to be valid once $z$ is a microscopic length, because then we'll need the correlation functions at wavevectors comparable to an inverse microscopic length, at which our hydrodynamic theory does not apply.

We now calculate the space and time-dependent velocity correlators

C^v_ij(-', z, z', t-t')≡⟨v_i(, z,t) v_j(', z',t')⟩

This can also be written in terms of the  in terms of the polarization $\hp$ and density $\rho$ correlation functions using the relation (<ref>)  between the bulk fluid velocity field and the surface velocity, and our boundary condition (<ref>) for that surface velocity.
We find

\begin{eqnarray}
&&C^v_{ij}(\brp-\brp', z, z', t-t')=\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i\omega (t-t')-i{\bf q}\cdot ({\brp-\brp'})]\langle v_i({\bf q},z,\omega)v_j(-{\bf q},z',-\omega) \rangle\nonumber \\
&&=v_0^2\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i\omega (t-t')-i{\bf q}\cdot ({\brp-\brp'})
-q(z+z')] C_{pp}({\bf q},\omega)\bigg[\delta_{iy}\delta_{jy}\left(1-{q_y^2( z+z') \over q}\right)+q_y^2z z'
\left({q_iq_j\over q^2}+\delta_{iz}\delta_{jz'}\right)\bigg] \nonumber \\
&&+ v_0 v_a'\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i\omega (t-t')-i{\bf q}\cdot ({\brp-\brp'})
-q(z+z')] C_{p\rho}({\bf q},\omega)\bigg[(\delta_{ix}\delta_{jy}+ \delta_{iy}\delta_{jx})\left(1-{q( z+z')
\over 2} + {2z z'q_x^2q_y^2 \over q^2} \right)\bigg] \nonumber \\
&&+v_a'^2\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i\omega (t-t')-i{\bf q}\cdot ({\brp-\brp'})
-q(z+z')] C_{\rho\rho}({\bf q},\omega)\bigg[\delta_{ix}\delta_{jx}\left(1-{q_x^2( z+z') \over q}\right)
+q_x^2z z'\left({q_iq_j\over q^2}+\delta_{iz}\delta_{jz'}\right)\bigg]
\label{vcorrreal}
\nonumber\\
\end{eqnarray}

where $C_{pp} ({\bf q},\omega)$ , $C_{p\rho}({\bf q},\omega)$, and $C_{\rho\rho}({\bf q},\omega)$ are given by equations (<ref>) and (<ref>), (<ref>), (<ref>),   (<ref>), and(<ref>),     of the introduction.

One of the best experimental probes of this correlation function is the diffusion of tracer particles. Consider first neutrally buoyant tracer particles, which, in the absence of diffusion, sit at a constant $z$, and are advected along $\hx$ at a speed $v_0$ by the mean motion of the fluid. If we wish to study the diffusion of these particles on a time scale small compared to the time required for them to diffuse a distance comparable to $z$, we therefore only need the velocity correlations (<ref>) at  $\br-\br'=v_0(t-t')\hx$, $z=z'$. In this limit, (<ref>) reduces to

\begin{eqnarray}
C^v_{ij}(v_0(t-t')\hx, z, z, t-t') = C^{pp}_{ij}(v_0(t-t')\hx, z, z, t-t') +
C^{p\rho}_{ij}(v_0(t-t')\hx, z, z, t-t') + C^{\rho\rho}_{ij}(v_0(t-t')\hx, z, z, t-t') \,,
\nonumber\\
\label{cvij}
\end{eqnarray}

\begin{eqnarray}
&& C^{pp}_{ij}(v_0(t-t')\hx, z, z, t-t') = \nonumber \\
&& v_0^2\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i(\omega-v_0q_x)(t-t')
-2qz] C_{pp}({\bf q},\omega)\bigg[\delta_{iy}\delta_{jy}\left(1-{2q_y^2z\over q}\right)+q_y^2z^2\left({q_iq_j\over q^2}+\delta_{iz}\delta_{jz}\right)\bigg] \nonumber \\
&&C^{p\rho}_{ij}(v_0(t-t')\hx, z, z, t-t') = \nonumber \\
&& v_a'v_0\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i(\omega-v_0q_x)(t-t')
-2qz] C_{p\rho}({\bf q},\omega)\bigg[(\delta_{ix}\delta_{jy}+ \delta_{iy}\delta_{jx})\left(1-qz + {2z^2 q_x^2q_y^2 \over q^2} \right)\bigg] \nonumber \\
&& C^{\rho\rho}_{ij}(v_0(t-t')\hx, z, z, t-t') = \nonumber \\
&& v_a'^2\int\frac{d\omega}{2\pi}\frac{d^2q}{(2\pi)^2} \exp [i(\omega-v_0q_x)(t-t')
-2qz] C_{\rho\rho}({\bf q},\omega)\bigg[\delta_{ix}\delta_{jx}\left(1-{2q_x^2z\over q}\right)+q_x^2z^2\left({q_iq_j\over q^2}+\delta_{iz}\delta_{jz}\right)\bigg]
\label{vcorrneutral}
\nonumber\\
\end{eqnarray}

It is straightforward to see from this expression that the tensor $C^{pp}_{ij}(v_0(t-t')\hx, z, z, t-t')$ is purely diagonal. Begin by noting that the only term in the integrand in (<ref>) that has an off-diagonal component is the $q_iq_j$ term. Since $\bq$ is a vector in the $xy$ plane, this term vanishes if either index $i$ or $j$ is $z$. Hence, the only off-diagonal terms are $i=x$, $j=y$, or vis-versa. In either case, the integrand then becomes odd in $q_y$ (recall that  both $C_{pp}(\bq, \omega)$ and $q=|\bq|$ are even in $\bq$. These are  the only $q_y$ dependent pieces of the rest of the integrand (it appears in the $2qz$ term in the argument of the first exponential)).

Hence, the integral in (<ref>), and, therefore, the correlation function $C^{pp}_{ij}(v_0(t-t')\hx, z, z, t-t')$, itself vanishes if $i=x$, $j=y$, or vice-versa. Since we have already established that the off-diagonal components of  $C^{pp}_{ij}(v_0(t-t')\hx, z, z, t-t')$ with one of the indices equal to $z$ also vanish, this completes the proof that $C^{pp}_{ij}(v_0(t-t')\hx, z, z, t-t')$ is, as claimed earlier,  purely diagonal. Similarly we can show from (<ref>) that $C^{\rho\rho}_{ij}(v_0(t-t')\hx, z, z, t-t')$ is also a
diagonal tensor.

However $C^{p\rho}_{ij}(v_0(t-t')\hx, z, z, t-t')$ is an off-diagonal symmetric tensor, as is easily seen from
(<ref>). Only the components $i=x$, $j=y$, and $i=y$, $j=x$ are non-zero, and

This shows that $C^v_{ij}(v_0(t-t')\hx, z, z, t-t')$ as a whole is symmetric tensor with four independent
components. The long-time scaling behavior of the four non-zero independent components $C^v_{xx}(v_0(t-t')\hx, z, z, t-t')$, $C^v_{yy}(v_0(t-t')\hx, z, z, t-t')$ and
$C^v_{xy}(v_0(t-t')\hx, z, z, t-t')$  are essentially same, as shown below .  However, $C^v_{zz}(v_0(t-t')\hx, z, z, t-t')$ has a different behavior in the long time limit. This is significant: diffusion along any direction in the $xy$-plane is anomalous, whereas it is normal (non-anomalous) along the vertical or $z$-direction.

To see this, let us consider each of these  four non-zero components in turn, starting with $C^v_{xx}(v_0(t-t')\hx, z, z, t-t')$. Using (<ref>), we see that this can be expressed as

C^v_xx(v_0(t-t'), z, z, t-t') = C^pp_xx(v_0(t-t'), z, z, t-t')+ C^ρρ_xx(v_0(t-t'), z, z, t-t')

as $C^{p\rho}_{xx}(v_0(t-t')\hx, z, z, t-t')$ is zero, being a non-diagonal tensor itself.

We evaluate each of the two terms in (<ref>). From (<ref>), we see that
$C^{pp}_{xx}(v_0(t-t')\hx, z, z, t-t')$ is given by

C^pp_xx(v_0(t-t'), z, z, t-t')

=v_0^2∫dω/2πd^2q/(2π)^2 exp[i(ω-v_0q_x)(t-t')
-2qz] C_pp(q,ω)q_x^2q_y^2z^2q^2 .

Making the change of variables of integration from $(\bq,\omega)$ to dimensionless variables ${\bf Q}, \Omega$ given by
≡Qz  ,

and recalling that $C_{pp}({\bf q},\omega)$ obeys the scaling form (<ref>), we find that

C^ pp_xx(v_0(t-t'), z, z, t-t')

=(v_0^3z)∫dΩ/2πd^2Q/(2π)^2 exp[i(Ω-Q_x)
v_0(t-t')/z-2Q] F_pp(v_0ΩQ, θ_Q)Q_x^2Q_y^2Q^4 ,

where $\theta_{\bf Q}$ is the angle between the rescaled, dimensionless vector ${\bf Q}$ and the direction of mean polarization $\hx$ (which is, of course, just the same as the angle $\theta_{\bf q}$between the original vector ${\bf q}$ and the direction of mean polarization $\hx$, since our rescaling (<ref>) was isotropic.

Using our expression (<ref>) for $F_{pp}$, we have

F_pp(v_0ΩQ, θ_Q)=2D_p(v_0ΩQ-v_ρcosθ_Q)^2(v_0ΩQ-c_+(θ_Q))^2(v_0ΩQ-c_-(θ_Q))^2+(v_0ΩQ
ψ(1,θ_Q)-v_ρψ(φ,θ_Q) cosθ_Q)^2  .

Factoring out $D_p v_0^2$ from the numerator of this expression, and $v_0^4$ from the denominator, we see that this can be rewritten in terms of a completely dimensionless scaling function of dimensionless arguments as
F_pp(v_0ΩQ, θ_Q)=(D_pv_0^2) H_pp(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})  ,

where we have defined

H_pp(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})=2(ΩQ-v_ρv_0 cosθ_Q)^2(ΩQ-c_+(θ_Q)v_0)^2( ΩQ-c_-(θ_Q)v_0)^2+(ΩQ ψ(1,θ_Q) v_0 -v_ρv_0ψ(φ,θ_Q) v_0 cosθ_Q)^2  .

Some of the dependence of this function on the dimensionless ratios ${v_\rho\over v_0}$, ${v_p\over v_0}$, ${\gamma\over v_0}$, and ${c_0\over v_0}$ is hidden in the ratios ${c_\pm(\theta_{\bf Q})\over v_0}$, since $c_\pm(\theta_{\bf Q})$ depend on $v_\rho$, $v_p$, $\gamma$, and $c_0$, as displayed in equation (<ref>).

Using (<ref>) and (<ref>) in (<ref>), we find that
$C^{pp}_{xx}(v_0(t-t')\hx, z, z, t-t')$
itself obeys a simple scaling law:
C^ pp_xx(v_0(t-t'), z, z, t-t')=(D_pv_0z)F^ pp_xx(v_0|t-t'|z)  ,

where the dimensionless scaling function

F^ pp_xx(u_r)

=∫dΩ/2πd^2Q/(2π)^2 exp[i(Ω-Q_x)u-2Q]H_pp(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})Q_x^2Q_y^2Q^4 ,

where the dimensionless scaling variable $u_r\equiv v_0|t-t'|/ z$.

The limiting behaviors of $F^{ pp}_{xx}(u_r)$ for small and large $u_r$ are easy to obtain;, and useful: $u_r\ll 1$ and $u_r\gg 1$ give information about the velocity  correlations, and from those, the displacement correlations, in the ballistic and diffusive limits, respectively, as we'll soon show. Intermediate values of $u_r\sim {\cal O}(1)$  correspond to the crossover between the two regimes. For  $u\ll1$, which corresponds to time differences obeying $v_0|t-t'|\ll z$,  the $e^{-2Q}$ factor in the integral (<ref>) kills the integrand before (that is, at smaller $Q$ than) the $(\Omega-Q_x)u$ term becomes important. Hence, the integral, and $F^{ pp}_{xx}(u_r)$ itself, become independent of $u$ in this limit; that is,
F^ pp_xx(u_r)=A_x      for      u_r≪1 ,

where the constant $A^{ pp}_{ xx}$ is given by

=∫dΩ/2πd^2Q/(2π)^2 e^-2QH_pp(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})Q_x^2Q_y^2Q^4 .

and is a function of all of the ratios ${v_\rho\over v_0}$, ${v_p\over v_0}$, ${\gamma\over v_0}$, and ${c_0\over v_0}$.  $A^{pp}_{xx}$  will be of ${\cal O}(1)$ when all of these ratios are of ${\cal O}(1)$.

In the opposite limit $u_r\gg1$, which corresponds to time differences obeying $v_0|t-t'|\gg z$,  the $e^{i(\Omega-Q_x)u}$ factor in the integral (<ref>) kills (by oscillation) all contributions to the integral coming from $Q\gtrsim{1\over u_r}\ll1$.  Hence, in the dominant region of the integral, $Q\ll1$, and so we can drop the $2Q$ term in the argument of the exponential. Doing so gives

F^ pp_xx(u_r)

=∫dΩ/2πd^2Q/(2π)^2 exp[i(Ω-Q_x)u]H_pp(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})Q_x^2Q_y^2Q^4 .

This integral can easily be done with one further change of variables:
Q≡Q'u_r  ,

which gives
F^ pp_xx(u_r)=B^pp_xxu_r^3       for      u_r≫1 ,

where the constant  $B^{pp}_{xx}$  is given by

=∫dΩ'/2πd^2Q'/(2π)^2 exp[i(Ω'-Q'_x)]H_pp(Ω'Q', θ_Q'; {v_ρv_0,v_pv_0,γv_0,c_0v_0})Q_x^'2Q_y^'2Q^'4 ,

like $A^{pp}_{xx}$,
again a function of all of the ratios ${v_\rho\over v_0}$, ${v_p\over v_0}$, ${\gamma\over v_0}$, and ${c_0\over v_0}$, and will again be of $O(1)$ when all of these ratios are of $O(1)$.

In summary,  the behavior of the scaling function $F^{ pp}_{xx}(u_r)$ is given by:

A^pp_xx   u_r≪1

B^pp_xxu_r^3   u_r≫1

We now evaluate $C^{\rho\rho}_{xx}(v_0(t-t')\hx, z, z, t-t')$, which from (<ref>) is expressed as
C^ρρ_xx(v_0(t-t'), z, z, t-t') = v_a'^2∫dω/2πd^2q/(2π)^2 exp[i(ω-v_0q_x)(t-t')
-2qz] C_ρρ(q,ω)[1-2q_x^2zq+ z^2q_x^4q^2]

Making  the same change of variables (<ref>)  as before, we obtain an scaling form for $C^{\rho\rho}_{ij}(v_0(t-t')\hx, z, z, t-t')$, given by
C^ρρ_xx(v_0(t-t'), z, z, t-t') = (D_pρ_c^2 v_a'^2v_0 z)F^ρρ_xx(v_0|t-t'|z)

where $F^{\rho\rho}_{xx}\left({v_0|t-t'|\over z}\right)$ is written as

=∫dΩ/2πd^2Q/(2π)^2 exp[i(Ω-Q_x)u_r-2Q]H_ρρ(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})
[1 Q^2 - 2 Q_x^2 Q^3 + Q_x^4 Q^4] .

where $H_{\rho\rho}\bigg({\Omega\over Q}, \theta_{\bf Q}; \bigg\{{v_\rho\over v_0},{v_p\over v_0},{\gamma\over v_0},{c_0\over v_0}\bigg\}\bigg)$ is defined as
H_ρρ(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0}) = 2(v_ρv_0)^2(sinθ_Q)^2 (ΩQ-c_+(θ_Q)v_0)^2( ΩQ-c_-(θ_Q)v_0)^2+(ΩQ ψ(1,θ_Q) v_0 -v_ρv_0ψ(φ,θ_Q) v_0 cosθ_Q)^2  .

The limiting behaviors of the scaling function $F^{\rho\rho}_{xx}(u)$ can be obtained by an almost identical analysis to that used for $F^{pp}_{xx}(u_r)$, giving:

A^ρρ_xx   u≪1

(B^ρρ_xx)_1u_r + (B^ρρ_xx)_2u_r^2 + (B^ρρ_xx)_3u_r^3   u_r≫1

=∫dΩ/2πd^2Q/(2π)^2 e^-2QH_ρρ(ΩQ, θ_Q; {v_ρv_0,v_pv_0,γv_0,c_0v_0})[1 Q^2 - 2 Q_x^2 Q^3 + Q_x^4 Q^4] .

(B^ρρ_xx)_1    =    ∫dΩ'/2πd^2Q'/(2π)^2 exp[i(Ω'-Q'_x)]H_ρρ(Ω'Q', θ_Q'; {v_ρv_0,v_pv_0,γv_0,c_0v_0})
1 Q^'2 ,

(B^ρρ_xx)_2    =    -2∫dΩ'/2πd^2Q'/(2π)^2 exp[i(Ω'-Q'_x)]H_ρρ(Ω'Q', θ_Q'; {v_ρv_0,v_pv_0,γv_0,c_0v_0})Q_x^'2Q^'3 ,
(B^ρρ_xx)_3    =    ∫dΩ'/2πd^2Q'/(2π)^2 exp[i(Ω'-Q'_x)]H_ρρ(Ω'Q', θ_Q'; {v_ρv_0,v_pv_0,γv_0,c_0v_0})
Q_x^'4Q^'4 ,

Like $A^{pp}_{xx}$ (<ref>) and $B^{pp}_{xx}$ (<ref>),  $A^{\rho\rho}_{xx}$ and $(B^{\rho\rho}_{xx})_1$, $(B^{\rho\rho}_{xx})_2$, and $(B^{\rho\rho}_{xx})_3$
are functions of all of the ratios ${v_\rho\over v_0}$, ${v_p\over v_0}$, ${\gamma\over v_0}$, and ${c_0\over v_0}$, and  will  all be of $O(1)$ when all of these ratios are of $O(1)$.

Taking (<ref>) and (<ref>) together, we find that $C^{v}_{xx}(v_0(t-t')\hx, z, z, t-t')$ is  given by
C^v_xx(v_0(t-t'), z, z, t-t') = (D_pv_0z)F^ pp_xx(v_0|t-t'|z) + (D_pρ_c^2 v_a'^2v_0 z)F^ρρ_xx(v_0|t-t'|z)

with $F^{ pp}_{xx}\left({v_0|t-t'|\over z}\right)$, and $F^{\rho\rho}_{xx}\left({v_0|t-t'|\over z}\right)$ given by (<ref>), and (<ref>) respectively.

Because of the ${1\over u_r}$ piece of the scaling function $F^{\rho\rho}_{xx}\left(u_r\right)$, the integral of this correlation function over $u$ for fixed $z$ diverges logarithmically at large $u_r$.
We will show  in section (<ref>)
below that this implies logarithmic superdiffusion (equations (<ref>) and (<ref>)  in the $x$ direction.

The derivation of similar results for the remain components of the velocity correlator tensor is given in Appendix <ref>.

§ ANOMALOUS DIFFUSION OF TRACERS

§.§ Neutrally buoyant tracers

A neutrally buoyant tracer is a particle that is simply carried around passively by the flows in the bulk fluid. That is, if the particle is at position $r_p(t)$ at time $t$, its  instantaneous velocity $\dot\br_p(t)$ is
_p(t)=(_p(t), t)  ,

where $\bv(\br,t)$ is the velocity field of the passive bulk fluid. Therefore, the instantaneous position of a tracer particle $\br_p(t)$ that starts at $\br_p(t=0)$ at some later time $t$ is given by

\begin{eqnarray}
\br_p(t)-\br_p(0) = \int^t_0
\bv(\br_p(t'), t') dt' \,.
\label{r tracer}
\end{eqnarray}

Since diffusion, and even the superdiffusion than we eventually will find, is a much slower process that ballistic transport, we can take, for a neutrally buoyant particle (i.e., one which has no net speed in the $z$-direction),
in (<ref>), where we have taken into account the fact that the passive fluid is, on average, flowing along the $\hx$ direction at a speed $v_0$. We have also simply written $z(t=0)$ as $z$. Doing so gives

\begin{eqnarray}
\br_p(t)-\br_p(0) = \int^t_0
\bv((v_0t'+x(t=0))\hx+y(0)\hy+z\hz, t') dt' \,.\nonumber\\
\label{r tracer 2}
\end{eqnarray}

Autocorrelating the $x$, $y$, and $z$ components of this equation with themselves gives
⟨(x(t)-x(0))^2⟩=∫_0^t dt'∫_0^t dt”   C^v_xx(t'-t”)  ,

⟨(y(t)-y(0))^2⟩=∫_0^t dt'∫_0^t dt”   C^v_yy(t'-t”)  ,

⟨(x(t)-x(0))(y(t)-y(0))⟩=∫_0^t dt'∫_0^t dt”   C^v_xy(t'-t”)=0  ,

⟨(z(t)-z(0))^2⟩=∫_0^t dt'∫_0^t dt”   C^v_zz(t'-t”)  ,

where the correlation functions $C^v_{xx}(t'-t'')$, $C^v_{yy}(t'-t'')$, $C^v_{xy}(t'-t'')$,and $C^v_{zz}(t'-t'')$ are precisely those we calculated in equations (<ref>), (<ref>), and (<ref>) respectively of
the preceding section.
We will now use these expressions to obtain the diffusive motion - that is, the small departure from the mean ballistic motion at velocity $v_0\hx$ of tracers in the $\hx$ direction - for all three directions.

We begin with the $x$-direction. Using our scaling form (<ref>)    for $C^v_{xx}(t'-t'')$, we can rewrite this as

⟨(x(t)-x(0))^2⟩=(D_pv_0z)∫_0^t dt'∫_0^t dt”   F^pp_xx(v_0|t'-t”|z) + (D_pρ_c^2v_a'^2v_0z)∫_0^t dt'∫_0^t dt”
F^ρρ_xx(v_0|t'-t”|z)   .

Now, changing variables of integration in the integrals above from $t''$ to $\delta t\equiv t''-t'$, we can rewrite this as

⟨(x(t)-x(0))^2⟩=(D_pv_0z)∫_0^t dt'∫_-t'^t-t' dδt   F^pp_xx(v_0|δt |z) + (D_pρ_c^2v_a'^2v_0z)∫_0^t dt'∫_-t'^t-t' dδt   F^ρρ_xx(v_0|δt |z)  .

Now  recalling, as shown by equations (<ref>)     and (<ref>),  that $F^{pp}_{xx}(u_r) \sim 1/u_r^3$, and $F^{\rho\rho}_{xx}(u_r)\sim 1/u_r$ for $u_r\gg 1$ (that is, for
$t\gg\left({z\over v_0}\right)$), we expect the first integral over $\delta t$ on the right hand side of (<ref>) to converge as its limits go to $\infty$, and the second integral to diverge. Thus
the second integral clearly dominates at long times and (<ref>) can be approximated by

⟨(x(t)-x(0))^2⟩≈(D_pρ_c^2v_a'^2v_0z)∫_0^t dt'∫_-t'^t-t' dδt   F^ρρ_xx(v_0|δt |z)  .

Changing variables of integration from $\delta t$ to $u\equiv {v_0\delta t\over z}$, we obtain for the integral over $\delta t$:
∫_-t'^t-t' dδt   F^ρρ_xx(v_0|δt |z)=(zv_0) ∫_-v_0t'z^v_0(t-t')z du_r   F^ρρ_xx(u_r) .

Using the fact that $F^{\rho\rho}_{xx}(u_r)$ is an even function of $u_r$, we can rewrite this as:
∫_-t'^t-t' dδt   F^ρρ_xx(v_0|δt |z)=(zv_0) {∫_0^v_0t'z du_r   F^ρρ_xx(u_r)+∫_0^v_0(t-t')z du_r   F^ρρ_xx(u_r)} .

It is convenient to break up the first integral in this expression into a part coming from $u_r<1$ and another part from $u_r>1$:
∫_0^v_0t'z du_r   F^ρρ_xx(u_r)=∫_0^1 du_r   F^ρρ_xx(u_r)+∫_1^v_0t'z du_r   F^ρρ_xx(u_r)

The integral from $0$ to $1$ in this expression is clearly $O(1)$, since, as we showed earlier, the integrand is. The second integral can be rewritten:
∫_1^v_0t'z du_r   F^ρρ_xx(u_r)
=    (B^ρρ_xx)_1∫_1^v_0t'z du_ru_r

+   ∫_1^v_0t'z du_r   (F^ρρ_xx(u_r)- (B^ρρ_xx)_1 u_r)  .

The second integral  on the right hand side of  this equation (<ref>)  is also ${\cal O}(1)$, as can be seen from the fact that the integrand is ${\cal O}(1)$, and the expression in parentheses falls off like ${1\over u_r^2}$ for large $u_r$, so the integral converges, even if the upper limit is taken to $\infty$. The first integral in (<ref>) is, of course, elementary (about as elementary as they come, in fact!). We thereby obtain:
∫_0^v_0t'z du_r   F^ρρ_xx(u_r)= (B^ρρ_xx)_1ln(v_0t'z)+O(1)  .

An almost identical analysis shows that
∫_0^v_0(t-t')z du_r   F^ρρ_xx(u_r)= (B^ρρ_xx)_1ln(v_0(t-t')z)+O(1)  .

Inserting (<ref>) and (<ref>) into  (<ref>) gives
∫_-t'^t-t' dδt   F^ρρ_xx(v_0|δt |z)=(zv_0) {
(B^ρρ_xx)_1 [ln(v_0t'z)+ln(v_0(t-t')z)]+O(1)} .

Inserting this in turn into (<ref>), we obtain:
⟨(x(t)-x(0))^2⟩=D_pρ_c^2v_a'^2 v_0^2∫_0^t dt'{ (B^ρρ_xx)_1 [ln(v_0t'z)+ln(v_0(t-t')z)]+O(1)} .

All of the integrals in this expression are elementary, yielding our final superdiffusive expression for the mean squared displacement in the $x$ direction:
⟨(x(t)-x(0))^2⟩= 2(B^ρρ_xx)_1 D_pρ_c^2v_a'^2 v_0^2 t (ln(v_0tz)+O(1)) ,

which the alert reader will recognize as equation (<ref>) of the introduction.

We have derived this result assuming that the particle has not moved appreciably in the $z$ direction from its original height $z$. This will  be true only for  $t\ll{z^2\over D_z}$, where $D_z$ is the ($z$-independent) diffusion constant we calculated  below in Eq. (<ref>). For much longer times, the particle will typically be a distance
$z(t)\sim\sqrt{D_z t}$ above the surface. Since the logarithm is quite insensitive to the precise position (i.e., to factors of $O(1)$ in this estimate of $z$, it will suffice, to leading logarithmic order, to replace
$z$ in (<ref>) with $z(t)\sim\sqrt{D_z t}$ in (<ref>). Doing so we obtain

⟨(x(t)-x(0))^2⟩= (B^ρρ_xx)_1 D_pρ_c^2v_a'^2 v_0^2 t (ln(v_0^2tD_z)+O(1)) ,

which that same  alert reader will recognize as equation (<ref>) of the introduction.

For short times $t\ll\left({z\over v_0}\right)$ (which, it should be noted, can actually get arbitrarily long as $z\to\infty$), the arguments of  $F^{pp}_{xx}$ and $F^{\rho\rho}_{xx}$ are always much less than $1$ throughout the region of integration over $t'$ and $t''$ in (<ref>). Therefore,
$F^{pp}_{xx}$ and $F^{\rho\rho}_{xx}$  can be replaced in that integral by their small $u_r$ limits
$A^{pp}_{xx}$ and $A^{\rho\rho}_{xx}$ respectively. And because they are constants, we obtain ballistic scaling in this regime:

⟨(x(t)-x(0))^2⟩= [A^pp_xx(D_pv_0z) + A^ρρ_xx
(D_pρ_c^2v_a'^2v_0z) ]t^2   t≪(zv_0) .

The interested reader can easily verify that this ballistic behavior matches smoothly onto the $z$-independent diffusive behavior (<ref>) for $t\sim\left({z\over v_0}\right)$, up to logarithmic factors, as it should.

Virtually identical reasoning applies to  diffusion in the $y$ direction, with the result:
⟨(y(t)-y(0))^2⟩= 2(B^pp_yy)_1 D_p t (ln(v_0tz)+O(1)) ,

and at very long times,
⟨(y(t)-y(0))^2⟩= (B^pp_yy)_1 D_p t (ln(v_0^2tD_z)+O(1)) .

The above equation is easy to recognize as (<ref>) of the introduction.

As for $x$-direction, the motion in the $y$-direction is also ballistic for short times $t\ll\left({z\over v_0}\right)$:
⟨(y(t)-y(0))^2⟩= [A^pp_yy(D_pv_0z)+ A^ρρ_yy(D_pρ_c^2v_a'^2v_0z) ] t^2 t≪(zv_0) .

Once again, the interested reader  can easily verify that this ballistic behavior matches smoothly onto the $z$-independent diffusive behavior (<ref>) for $t\sim\left({z\over v_0}\right)$, as it should.

In the $z$-direction, we see that both $F^{pp}_{zz}(u_r) \sim 1/u_r^3$, and $F^{\rho\rho}_{zz}(u_r)\sim 1/u_r^3$ for
$u\gg 1$ or $t\gg\left({z\over v_0}\right)$. This implies  diffusive behavior at long times; i.e.

⟨(z(t)-z(0))^2⟩= 2D_z t   ,

D_z=D_p∫_0^∞ du_r   F^pp_zz(u) + D_pρ_c^2 v_a'^2 v_0^2∫_0^∞ du_r
F^ρρ_zz(u_r) ,

where $D_z$ converges to a finite value.

As for motion in the $x$ and $y$ direction,  at short times $t\ll\left({z\over v_0}\right)$  we obtain ballistic scaling for $z$:
⟨(z(t)-z(0))^2⟩=[A^pp_zz(D_pv_0z)+ A^ρρ_zz(D_pρ_c^2v_a'^2v_0z) ] t^2 t≪(zv_0) .

As for $x$,  this ballistic behavior for $z$ matches smoothly onto the $z$-independent diffusive behavior (<ref>) for $t\sim\left({z\over v_0}\right)$, as it should.

We can also calculate the off diagonal correlation
⟨(x(t)-x(0))(y(t)-y(0))⟩=∫_0^t dt'∫_0^t dt”   C^v_xy(t'-t”)  ,

which turns out to be zero as $C^v_{xy}(t'-t'')$ consists of an integral over q of $C_{p\rho}(\bq)$ times an even function of $q_y$. As we have seen from Eq. () that $C_{p\rho}(\bq)$ is an odd in $q_y$, this integral turns out to be zero. So
⟨(x(t)-x(0))(y(t)-y(0))⟩=0  ,

The above equation is easy to recognize as (<ref>) of the introduction.

§.§ Sedimenting tracers

For sedimenting particles, we first note that, since we just showed that diffusion in the $z$  direction is normal, and homogeneous, we need to consider diffusive motion in the $x$, and $y$ direction as well as the
$x-y$ plane. We'll start by considering sedimenting particles whose sedimentation speed $v_{\rm sed}<<v_0$. Such particles will spend a time of $O({z_0\over v_{\rm sed}})$ at a height of order $z_0$. Due to the aforementioned insensitivity of the logarithmic factor in (<ref>), (<ref>), and the
(<ref>), we can therefore accurately estimate the total mean squared $y$ displacement
$\langle(x(z=z_0)-x(z=0))^2\rangle$, $\langle(y(z=z_0)-y(z=0))^2\rangle$, $\langle(x(z=z_0)-x(z=0))
(y(z=z_0)-y(z=0))\rangle$ of a sedimenting particle that starts at $z=z_0$ and sinks at speed $v_{\rm sed}$ all the way down to the surface by simply replacing $t$ in (<ref>), (<ref>), and (<ref>) equations by ${z_0\over v_{\rm sed}}$, which is the time it takes the sedimenting tracer to sink to the bottom. Doing so gives
⟨(x(z=z_0)-x(z=0))^2⟩   =    2(B^ρρ_xx)_1 D_pρ_c^2v_a'^2 v_0^2 (z_0v_sed) (ln(v_0v_sed)+O(1)) ,

⟨(y(z=z_0)-y(z=0))^2⟩   =    2 B^pp_yy D_p (z_0v_sed) (ln(v_0v_sed)+O(1)) ,

⟨(x(z=z_0)-x(z=0))(y(z=z_0)-y(z=0))⟩   =     0

which are just the equations (<ref>-<ref>) of the introduction.

For more rapidly sedimenting particles (i.e., denser ones), for which $v_{\rm sed}\gg v_0$, the time scale $t\sim{z_0\over v_{\rm sed}}$ of the sinking of the sedimenting particle is in the ballistic regime $t\ll\left({z\over v_0}\right)$, and so we need to use
(<ref>) and(<ref>) with $t$ replaced with the sedimenting time ${z_0\over v_{\rm sed}}$. This gives
⟨(x(z=z_0)-x(z=0))^2⟩   =    A^pp_xx(D_pv_0 z_0v_sed^2) +
A^ρρ_xx(D_pρ_c^2v_a'^2z_0v_0v_sed^2)   ,

⟨(y(z=z_0)-y(z=0))^2⟩   =    A^pp_yy(D_pv_0 z_0v_sed^2) +
A^ρρ_yy(D_pρ_c^2v_a'^2z_0v_0v_sed^2) ,

These are equations (<ref>) and (<ref>) of the introduction.

We close this Section by noting that, for both neutrally buoyant particles and sedimenting tracers, the aspect ratio $\sqrt{(y(t)-y(0))^2/(x(t)-x(0))^2}$ is not exactly 1, but is a model parameter dependent ${\cal O}(1)$ number, reflecting the geometric anisotropy due to the $x$-direction being the preferred direction of orientation for the polar particles. Nonetheless, the aspect ratio is independent of time $t$. This means an anisotropy exponent $\zeta_{ani}$ that describes the relative scaling between distances measured along the $x$- and $y$-directions is unity (i.e., $\zeta_{ani}=1$),  unlike in the Toner-Tu model of flocking, for which
$\zeta_{ani}\ne1$ [13, 14] due to strongly relevant anharmonic effects.

The anomalous diffusion we find here  has some features common with the phenomenon of
“Taylor diffusion”[75], which occurs in parallel plate shear flow. One could obtain the
Taylor diffusion geometry from ours by replacing the active fluid layer at the bottom  by a flat
surface, moving at constant velocity, and our free top surface   with a
stationary rigid wall with no-slip boundary conditions. However, our anomalous diffusion  is
fundamentally different  from Taylor diffusion. One principal difference is that, in our case,
anomalous diffusion arises purely due to the fluctuations away from uniform alignment
and motion in the active layer. There is obviously no analog of this if one replaces the active layer in our problem with a uniformly
moving rigid body, as in Taylor diffusion.

§ RENORMALIZATION GROUP ARGUMENT AND IRRELEVANCE OF NON-LINEARITIES

So far, we have worked strictly with the linear theory. For dry active matter, it is well-known[13, 14, 15]  that non-linear effects radically change the long-wavelength behavior (indeed, it is only the effects of non-linearities that even make the ordered state in two dimensions possible). It therefore clearly behooves us to ask whether non-linearities have such important effects in our problem.

In this section, we will use a simple renormalization group power counting argument to show that they do not. In fact, the linear theory presented earlier is asymptotically exact at long distances.

There are several sources of nonlinearities. For instance, treating the various speeds and  pressures in the problem as functions of $\rho$, and then expanding in powers of $\delta\rho$ produces these nonlinear terms. In addition there are nonlinear terms which originate from the fixed length constant on the polarization $\bf p$.

In order to ascertain the relevance or irrelevance of the nonlinear effects, we need to consider the lowest order nonlinear terms in the equations for $p_y$ and $\rho$,  and in the “active partial slip” boundary condition (<ref>). The RG analysis we are about to present will make it clear that the most important terms at long distances are those with the fewest possible spatial derivatives, and the smallest number of fields.

There are two types of additional terms. The first type arises very straightforwardly from expanding, e.g, the pressures $P_{p,s}(\rho)$ and the velocity $v_a(\rho)$, to higher order in $\rho$. The second come from including the $p_kp_i$ piece of the transverse projection operator $T_{ki}=\delta_{ki}-p_kp_i$ in the full equation of motion (<ref>). The latter gives rise to the following extra  quadratic  order in
the fields $p_y$ and $\delta\rho$
contribution to $\pp_tp_y$:
_tp_y(,t)_quad=-αp_yv_sx(,t)  .

Since this term already has one power of $p_y$ multiplying $v_{sx}$, it is  obviously sufficient, to quadratic  order in the fields, to use our linear solution  (<ref>) for $v_{sx}$ in (<ref>). Since it is most convenient for our RG analysis to work in real space, we rewrite (<ref>) in real space, where it reads
δv_sx(,t)=(v'_a+(ζ̅-σ) _x)δρ(,t)-μηv_a'∫d^2 r'_⊥K_ppρ(-') δρ(',t)+ζ_20_yp_y(,t)-μηv_0∫d^2 r'_⊥K_pρ(-') p_y (',t)  ,
where the kernels ${\cal K}_{pp\rho}(\brp)$ and ${\cal K}_{p\rho}(\brp)$ are the inverse  two-dimensional Fourier transforms of
$\left(\frac{q^2+q_x^2}{q}\right)$ and $\left(\frac{q_xq_y}{q}\right)$ respectively[72], and are given, at large distances [72],  by
K_ppρ() ≈-3 x^22πr_⊥^5

K_pρ() ≈-3 xy2πr_⊥^5

as we demonstrate in appendix <ref>.

Including such terms, the equations of motion take the form

\begin{eqnarray}
\frac{\partial p_y}{\partial t} &=& -v_p \partial_x p_y - \gamma\int d^2 r'_\perp {\cal K}_{pp}({\brp-\brp'}) p_y (\brp')  - \frac{\gamma_\rho}{\rho_c}\int d^2 r'_\perp {\cal K}_{p\rho}({\brp-\brp'}) \delta\rho (\brp')-
\sigma_t \partial_y\delta\rho+ f_y \nonumber \\
&-& { \gamma_{_{NL}}}\int d^2 r'_\perp p_y (\brp) {\cal K}_{p\rho}({\brp-\brp'}) p_y (\brp')-{ g_{pp\rho1}}\int d^2 r'_\perp p_y (\brp) {\cal K}_{pp\rho}({\brp-\brp'}) \delta\rho (\brp')
\nonumber \\
&-& { g_{pp\rho2}}\int d^2 r'_\perp p_y (\brp') {\cal K}_{pp}({\brp-\brp'}) \delta\rho (\brp') - { \sigma_{_{NL}}} p_y \partial_x\delta\rho \nonumber\\ &+&  \frac{\lambda_1}{2}\partial_y p_y^2  +
\lambda_2 \delta\rho\partial_x p_y+{\lambda_3} p_y f_x,
\label{py-real}
\end{eqnarray}

\begin{eqnarray}
\frac{\partial \delta\rho}{\partial t} = -v_\rho\partial_x\delta\rho - v_p\rho_c\partial_y p_y - \lambda_4\partial_y (p_y\delta \rho) + \lambda_5 \partial_x(\delta\rho)^2 + {\boldsymbol\nabla}\cdot {f}_\rho \,.\label{rho-gen}
\end{eqnarray}

In equation (<ref>), the non-linear terms appear on the last three lines. The “bare” values of the parameters in these equations of motion (hereafter denoted by a superscript “$0''$) are related as follows:

\begin{eqnarray}
\label{rr0a}\\
&&g_{pp\rho1}^0= \gamma^0_\rho/\rho^0_c \,,
\\
&&g_{pp\rho2}^0= \alpha  v'_a\mu\eta\,,
\\
&&\sigma^0_{_{NL}}= \sigma^0_t \,,
\\
&&\lambda_1^0=(\nu^0_1 - \lambda^0_{\rho v})v_0 - \lambda^0 \,,
\\
&&\lambda^0_2=\bigg(\frac{\nu^0_1+1}{2} \bigg) v'^0_a- \lambda^0_{\rho v}(v_a')^0 \,,
\\
&&\lambda^0_3=-1 \,,
\\
&&\lambda^0_4=v_\rho \,,
\\
&&\lambda^0_5=-\rho'^0_e(v_a')^0 \,.
\label{bare}
\end{eqnarray}

However, none of these parameters  will  continue to maintain these relations to the other parameters upon renormalization, which is why we have introduced them as independent parameters in the equations of motion (<ref>) and (<ref>).

Here, the kernel ${\cal K}_{pp}(\brp)$ is the inverse  two-dimensional Fourier transform of $\left(\frac{q^2+q_y^2}{q}\right)$,
and is given, at large distances [72],  by
K_pp() ≈-3 y^22πr_⊥^5  ,

as we demonstrate in appendix <ref>. Notice that we have restored  the  number-conserving  noise ${\bf f}_\rho$ in (<ref>). That noise  is assumed to be of zero-mean and Gaussian-distributed with  variance

\begin{equation}
\langle f_{\rho\,i}({\bf x},t) f_{\rho\,j}(({ {\bf x}',t'})\rangle = 2D_\rho { \delta_{ij}}\delta({\bf x}{ -{\bf x}'})\delta(t{ -t'}).
\label{frhocorr}
\end{equation}

We will now assess the importance of the non-linear terms in
these equations of motion using the dynamical renormalization group (DRG). Readers interested in a more complete and pedagogical  discussion of the DRG are referred to [73] for the details of this general approach.

This approach begins by decomposing the Fourier modes of the fields $p_y({\bf q},\omega)$ and $\rho({\bf q},\omega)$, and the noises $f_p$ and $f_\rho$ into a rapidly varying parts $p_y^>({\bf q},\omega)$,  $\rho^>({\bf q},\omega)$, $f_p^>({\bf q},\omega)$  and  $f_\rho^>({\bf q},\omega)$, and   slowly varying parts $p_y^<({\bf q},\omega)$ and  $\rho^><{\bf q},\omega)$ $f_p^<({\bf q},\omega)$  and  $f_\rho^<({\bf q},\omega)$. The rapidly varying parts are supported in the momentum shell $\Lambda\ee^{-\dd\ell}<q<\Lambda$, where $\dd\ell$ is an infinitesimal
and $\Lambda$ is the ultraviolet cutoff. The slowly varying part is supported in $q<\Lambda\ee^{-\dd\ell}$.

The  DRG procedure then consists of two steps. In step 1, we eliminate the rapidly varying parts $p_y^>({\bf q},\omega)$  and  $\rho^>({\bf q},\omega)$ from the equations of motion.
We do this by solving them iteratively for  $p_y^>({\bf q},\omega)$  and  $\rho^>({\bf q},\omega)$. This solution is a series in the non-linearities
which depends on the slow fields $p_y^<({\bf q},\omega)$ and  $\rho^<{\bf q},\omega)$. We substitute this solution into the equations of motion for the slow fields, and average over  the short wavelength components   $f^>_y({\bf q},\omega)$ and $\bff^>_\rho({\bf q},\omega)$ of the noises $f_y$ and $\bff_\rho$, which gives a closed EOM  for the slow fields  $p_y^<({\bf q},\omega)$ and  $\rho^<{\bf q},\omega)$.

Step 2 consists of rescaling  space and time as follows: ${\bf r}_{_\perp}'={\bf r}_\perp\ee^{-\dd\ell},\,t'=t\ee^{-z\dd\ell}$, where we will choose the “dynamical exponent” $z$ for our convenience. The rescaling of $\br$ has the effect of restoring the ultraviolet cutoff $\Lambda$ to its original value.

We simultaneously rescale the fields $p_y({\bf r}_{_\perp},t)$ and $\rho({\bf r}_{_\perp},t)$ according to

\begin{equation}
p_y({\bf r}_{_\perp},t)=e^{\chi_pd\ell}p_y({\bf x'},t'),\,\rho({\bf r}_{_\perp},t)=e^{\chi_\rho d\ell}\rho({\bf x'},t') \,
\end{equation}

where we will also choose the field rescaling exponents $\chi_p$ and $\chi_\rho$ for our convenience.
We then reorganize the resultant EOM so that it has the same form as  our original equations of motion, but with all of the coefficients changed, or, to use the standard jargon, “renormalized”.
This process is then repeated. The result is a set of differential recursion relations for the various parameters in the equations of motion (<ref>) and (<ref>), which are:

\begin{eqnarray}
&&{dD_p\over d\ell}=(z-2-2\chi_p)D_p +{\rm nonlinear \, corrections}\,,
\label{rr0}\\
&&{dD_\rho\over d\ell}= (z-4-2\chi_\rho)D_\rho+{\rm nonlinear \, corrections} \,,
\\
&&{d\sigma_t\over d\ell}= (z-1+\chi_\rho-\chi_p)\sigma_t+{\rm nonlinear \, corrections} \,,
\\
&&{d\gamma\over d\ell}=(z-1)\gamma +{\rm nonlinear  \, corrections} \,,
\\
&&{d\over d\ell}\left({\gamma_\rho\over\rho_c}\right)=(z+\chi_\rho-\chi_p-1)\left({\gamma_\rho\over\rho_c}\right)+{\rm nonlinear  \, corrections} \,,
\\
&&{d v_p\over d\ell}=(z-1)v_p+{\rm nonlinear  \, corrections} \,,
\\
&&{d v_\rho\over d\ell}=(z-1)v_\rho+{\rm nonlinear  \, corrections} \,,
\\
&& {d (v_p\rho_c)\over d\ell}=(\chi_p-\chi_\rho+z-1)v_p\rho_c+{\rm nonlinear  \, corrections} \,,
\\
&&{ {d\gamma_{_{NL}}\over d\ell}=(z+\chi_p-1)\gamma_{_{NL}}+{\rm nonlinear  \, corrections}} \,,
\\
&&{ {dg_{pp\rho1}\over d\ell}=(z+\chi_\rho-1)g_{pp\rho1}+{\rm nonlinear  \, corrections}} \,,
\\
&&{ {dg_{pp\rho2}\over d\ell}=(z+\chi_\rho-1)g_{pp\rho2}+{\rm nonlinear  \, corrections}} \,,
\\
&&{ {d\sigma_{_{NL}}\over d\ell}=(z+\chi_\rho-1)\sigma_{_{NL}}+{\rm nonlinear  \, corrections}} \,,
\\
&&{d\lambda_1\over d\ell}=(z+\chi_p-1)\lambda_1+{\rm nonlinear  \, corrections} \,,
\label{rr}
\end{eqnarray}

dλ_2dℓ=(z+χ_ρ-1)λ_2+nonlinear    corrections  ,

dλ_3dℓ=(χ_p)λ_3+nonlinear    corrections  ,

dλ_4dℓ=(z+χ_p-1)λ_4+nonlinear    corrections  ,

dλ_5dℓ=(z+χ_ρ-1)λ_4+nonlinear    corrections  .

The terms we have explicitly displayed in the above equations all come trivially from the rescaling of length, time, and fields. The “non-linear corrections” denote corrections arising due to the couplings between the “fast” and “slow” modes that arise in the equations of motion  (<ref>) and (<ref>) as a result of the non-linear terms in those equations. Since they arise from those non-linear terms, these terms must vanish when those non-linear terms do, and must, by continuity, be small when the coefficients $\lambda_i$ of those non-linear terms are small.

Of course, how small the $\lambda_i$ have to be to ensure that the non-linear corrections in the above recursion relations are negligible depends on the linear parameters in the equations of motion, since those determine the size of the fluctuations in the fields. This suggests a very simple, and very standard, way  to decide if the non-linear terms are important at long distances and times: simply choose the rescaling exponents $z$ and $\chi_{p,\rho}$ to keep the linear parameters that control the fluctuations fixed. If, with this choice of  $z$ and $\chi_{p,\rho}$, the $\lambda_i$ all then flow to zero as $\ell\to\infty$, then the non-linear terms are guaranteed to be unimportant at long distances. This is the approach we will now take.

We begin by noting that the linear parameters which control the size of the fluctuations in the linear theory are $D_p$, $v_{p,\rho}$, $\sigma_t$,$\gamma$, ${\gamma_\rho\over\rho_c}$, and $v_p\rho_c$, as we showed in section (<ref>)      above, where we calculated the fluctuations in the linear theory.

If we assume that all of the $\lambda_i$ are initially small enough that the non-linear corrections in the recursion relations (<ref>)-(<ref>) are negligible, then we can keep the above parameters fixed by choosing
χ_p=χ_ρ=-12,     z=1.

The choice $z=1$ implies that the scaling  argument of the scaling functions $F_{pp}, F_{\rho\rho}$ and $F_{p\rho}$ should be $\omega/q$  (since for general $z$, it would be $\omega/q^z$).

Using these choices (<ref>) in the recursion relation for $D_\rho$, and again neglecting the non-linear terms, we see that
dD_ρdℓ=-2D_ρ ,

which clearly shows that $D_\rho(\ell\to\infty)\to0$. Thus, $D_\rho$ is irrelevant at long distances and times, which justifies our neglect of it in our earlier, linear analysis.

We also find that all $\lambda_i$, $i=1\to5$ have same recursion relation:

dλ_idℓ=-(12)λ_i+nonlinear    corrections  ,
and that the other three non-linear coefficients $\gamma_{_{NL}}$, $g_{pp\rho}$, and $\sigma_{_{NL}}$ have the same RG eigenvalues:

\begin{eqnarray}
&&{ {d\gamma_{_{NL}}\over d\ell}=-\left({1\over2}\right)\gamma_{_{NL}}+{\rm nonlinear  \, corrections}} \,,
\nonumber\\
&&{ {dg_{pp\rho1}\over d\ell}=-\left({1\over2}\right)g_{pp\rho1}+{\rm nonlinear  \, corrections}} \,,
\nonumber\\
&&{ {dg_{pp\rho2}\over d\ell}=-\left({1\over2}\right)g_{pp\rho2}+{\rm nonlinear  \, corrections}} \,,
\nonumber\\
&&{ {d\sigma_{_{NL}}\over d\ell}=-\left({1\over2}\right)\sigma_{_{NL}}+{\rm nonlinear  \, corrections}}  \,.
\nonumber\\
\label{NLrg}
\end{eqnarray}

Therefore, all of these nonlinearities are also irrelevant at long distances, at least if they are initially small.

Since $\chi_{p, \rho}$ are both $<0$, any terms with more fields are less relevant. Likewise, any fields with more gradients are also less relevant. This includes all possible other nonlinearities. Therefore, all non-linearities are irrelevant. This implies that our linear results are asymptotically exact at long length and time scales.

§ SUMMARY AND CONCLUSIONS

In this paper, we have studied the  the stability and fluctuations of a large polar-ordered flock at a solid-liquid interface, which is a natural intermediate case between the two previously studied  cases of  wet and dry polar  active  fluids. Such a flock is affected by both the friction force from the solid substrate underneath and the long range hydrodynamic interaction mediated by the overlying passive, isotropic bulk fluid.   As a result, such a flock is simultaneously momentum nonconserving, but affected by the hydrodynamic interactions of the bulk surrounding fluid. Friction with the substrate also breaks Galilean invariance.