# Spinodal decomposition and phase separation in polar active matter

Maxx Miller<EMAIL_ADDRESS>Department of Physics and Institute for
Fundamental Science, University of Oregon, Eugene, OR $97403$ John Toner
<EMAIL_ADDRESS>Department of Physics and Institute for Fundamental Science,
University of Oregon, Eugene, OR $97403$

###### Abstract

We develop and study the hydrodynamic theory of flocking with autochemotaxis.
This describes large collections of self-propelled entities all spontaneously
moving in the same direction, each emitting a substance which attracts the
others (e.g., ants). The theory combines features of the Keller-Segel model
for autochemotaxis with the Toner-Tu theory of flocking. We find that
sufficiently strong autochemotaxis leads to an instability of the uniformly
moving state (the “flock”), in which bands of different density form moving
parallel to the mean flock velocity with different speeds. These bands, which
are reminiscent of ant trails, coarsen over time to reach a phase-separated
state, in which one high density and one low density band fill the entire
system. The same instability, described by the same hydrodynamic theory, can
occur in flocks phase separating due to any microscopic mechanism (e.g.,
sufficiently strong attractive interactions). Although in many ways analogous
to equilibrium phase separation via spinodal decomposition, the two steady
state densities here are determined not by a common tangent construction, as
in equilibrium, but by an uncommon tangent construction very similar to that
found for motility induced phase separation (MIPS) of disordered active
particles. Our analytic theory agrees well with our numerical simulations of
our equations of motion.

##  I Introduction and Motivation

Two of the principle mechanisms of ordered biological motion are “taxis” and
“flocking”Vicsek  _et al._ (1995); Toner and Tu (1995); Tu  _et al._ (1998);
Toner and Tu (1998); Toner  _et al._ (2005). “Taxis” is motion in response to
an external stimulus. One extremely large class of this type is movement
directed preferentially along the gradient of some quantity. Examples include:
“thermotaxis” (following temperature gradients)Pline  _et al._ (1988) and
“phototaxis” (light intensity) Anderson and Mackie (1977). In this paper, we
focus on the extremely common example of “chemotaxis”, in which motion is
directed along the gradient of the concentration of a chemical substance.
Almost all motile organisms perform chemotaxis Stock and Baker . Examples
include: Escherichia coliMesibov and Adler (1972), spongy mothsDoane and
McManus , spermHildebrand and Kaupp (2005), locusts Topaz  _et al._ (2012),
and antsCouzin and Franks (2003); Dussutour  _et al._ (2009).

In “autochemotaxis”, the chemical substance in question is secreted by the
moving creatures themselves. Clearly, this can lead to aggregation of those
creaturesKeller and Segel (1971a, 1970, b); Horstmann (2003); Stark (2018).

Another ubiquitous, and familiar, form of motion is the phenomenon of
flocking. Flocking is the collective coherent motion of a large number of
organisms. Examples include: a school of fish, a swarm of bees, or a flock of
birds. Such motion can occur any time creatures exhibit a tendency to follow
their neighbors. Strikingly, even purely short ranged interactions prove to be
sufficient to generate long-ranged orientational order, even in an arbitrarily
large flock, even in two dimensionsToner and Tu (1995); Tu  _et al._ (1998);
Toner and Tu (1998); Toner  _et al._ (2005), in contrast to equilibrium
systems, in which the development of long-ranged orientational order in a
rotation invariant system with purely short-ranged interactions is forbidden
by the “Mermin-Wagner-Hohenberg Theorem”Mermin and Wagner (1966).

Like chemotaxis, flocking can be observed on a large range of length scales
from micrometers, as in, (e.g., slime molds Dormann and Weijer (2006), to
kilometers (e.g., locust plagues) Topaz  _et al._ (2012)). Locust plagues can
contain billions of locusts, can occupy up to a thousand square kilometers,
and are capable of causing billions of dollars of damageTopaz  _et al._
(2012).

Hydrodynamic theories of autochemotaxis without flockingKeller and Segel
(1971a, 1970, b); Horstmann (2003); Stark (2018) and flocking without
autochemotaxisToner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998);
Toner  _et al._ (2005); Toner (2012a) have been successful in quantitatively
explaining each of their respective phenomena in isolation. In this paper, we
develop and study a hydrodynamic theory for systems in which both effects are
present. We will show the presence of both effects can lead to formation of
trails (e.g., ant trailsCouzin and Franks (2003)).

Specifically, we consider active self-propelled particles (hereafter called
“boids”) that move over a frictional substrate, so that momentum is not
conserved and time reversal symmetry is broken. The underlying dynamical rules
are assumed to be rotation and translation invariant.

We investigate the stability in the presence of autochemotaxis of an “active
homogeneous polar ordered state” which spontaneously breaks the underlying
rotation invariance by having all of the boids move, on average, in the same
direction (i.e., “flock”). The homogeneity of the state means that translation
invariance is not broken.

The chemotaxic element of the hydrodynamic theory is the same as that of the
Keller-Segel modelKeller and Segel (1971a, 1970, b); Horstmann (2003).
Specifically, boids release a chemical signal. This chemical signal can be a
substance that either repels neighboring boids (i.e., is a “chemo-repellent”)
or attracts them (i.e., is a “chemo-attractant”). In this paper, we study the
chemo-attractant case, in which, in addition to following their neighbors,
each boid also tends to follow the local chemo-attractant concentration
gradient. The chemo-attractant diffuses, and decays with a finite lifetime
$\tau$.

The “strength” of the autochemotaxis can therefore be increased by increasing
any one of three parameters: the chemo attractant lifetime $\tau$, the rate of
production of the chemo attractant by the boids $\gamma$, or the coupling of
their motion to the gradient of the concentration of the chemo-attractant.

The linear stability analysis of the hydrodynamic theory that we perform in
this paper shows that, like disordered active systems undergoing only
autochemotaxis Stark (2018), the spatially homogeneous state of ordered flocks
(or, to use the technical term, “polar ordered active fluids”) becomes
unstable when the autochemotaxis is sufficiently strong. However, in contrast
to this instability in disordered systems, in ordered flocks the instability,
near threshold, is extremely anisotropic. Specifically, plane wave modulations
of the density with wavevector $\mathbf{q}$ are only unstable in a narrow
region of $\mathbf{q}$ space, extended in the direction perpendicular to the
direction of flock motion, as illustrated in figure 1.

Figure 1:  The region of instability in $\mathbf{q}$ space, for autochemotaxic
systems very close to the instability threshold on the unstable side. Note
that the instability only occurs at wavevectors very nearly perpendicular to
the mean direction of flock motion, which is the horizontal ($x$)- axis in
this figure. Specifically, the width along $q_{x}$ of the unstable region
scales like $\epsilon^{2}$, while its length perpendicular to $x$ scales like
${\epsilon}$, where $\epsilon$ is a measure of the distance from the
instability threshold (defined more precisely below). As a result, the
instability is towards forming “bands” in real space running parallel to the
direction of mean flock motion, as illustrated in figures 2 and 3.

As a result, the instability is towards forming “bands” in real space running
parallel to the direction of mean flock motion. This is illustrated in figures
2 and 3.

The instability we find here is completely different from the “banding
instability” Bertrand and Lee (2022); Bertin  _et al._ (2009); Mishra  _et
al._ (2010); Toner (2012a). The banding instability occurs in a completely
different region of parameter spaceToner (2012a), is driven by completely
different physicsBertrand and Lee (2022); Bertin  _et al._ (2009); Mishra  _et
al._ (2010); Toner (2012a), and is characterized by the formation of bands
perpendicular to the direction of mean flock motion.

Figure 2: The “band” structure of the instability at intermediate times. The
density is only modulated along one of the directions (which we call $y$, and
which is indicated in the figure) perpendicular to the direction $\hat{x}$ of
mean flock motion. A plot of the modulation of the density along the
$y$-direction (e.g., along the path of the $y$-axis shown in the figure) is
given on 3. Figure 3: A plot of the density for a typical numerical solution
of the hydrodynamic equations for 1d configurations, as described in the text,
at intermediate times, with weak non-linear terms. The dashed blue curve is
the initial condition. The orange curve is the density at the last time step
our simulation reached. The density is only modulated along one of the
directions (which we call $y$, and which is indicated in the figure)
perpendicular to the direction $\hat{x}$ of mean flock motion. Here $\ell$ is
the linear spatial extent of our system in the $y$-direction (which is
periodic). The parameters used for this solution were $\Gamma,\Lambda=0.25$,
and the solution was iterated for 3000 units of time. The choice of $\rho_{0}$
was arbitrary.

To investigate what happens once the density and velocity perturbations
induced by this instability became too large to be treated by our linear
stability analysis, we found both analytic steady-state and numerical time-
dependent solutions to our equations of motion. We did so for “one-
dimensional” (1d) solutions: that is, solutions in which the fields only
depend on time and a single Cartesian coordinate perpendicular to the
direction of mean flock motion. This one-dimensional restriction is justified,
at least in the early stages, by the aforementioned fact that the instability
only occurs for modes with their wavevectors nearly perpendicular to the
direction of mean flock motion.

Our numerical solutions show that the these bands fairly quickly evolve into a
set of well-separated parallel “plateaus” of almost constant density, as
illustrated in figure 3. The evolution then continues via the merger of high
density bands with other high density bands and low density bands with other
low density ones. This suggests that bands of like density are attracted to
each other; however, it appears that this interaction falls off quickly with
distance, because the merging process gets very slow. Indeed, the mergers
become so slow that our 1d simulations simply cannot run long enough (at least
in the time we’re willing to wait!) to reach what we believe will be the
ultimate steady state of the system. That state is predicted by our analytic
solution to be phase separation into one high density band and one low density
band, moving parallel to each other (and in the direction of the original
flock motion) at slightly different speeds. This final geometry is illustrated
in figure 4, in which we compare our analytic solution of the steady state
equations to the numerical, time-dependent solution at a sufficiently late
time.

Figure 4: The density profile of the final phase separated steady state. The
long-time limit of our numerical solution, and our fit curve are both plotted.
The solid orange curve is the numerical solution. The dashed blue curve is the
fit curve. The fit curve is $\tanh{(y/\sqrt{1.0263})}+0.0355$. Although the
form of the fit curve is not derived analytically, its large and small $y$
limits are taken from our analytic theory. As can be seen, the agreement
between the analytic theory and the numerical solution is good. The parameters
used are
$\frac{1}{\rho_{0}}\sqrt{\frac{m}{u}},\frac{\lambda}{D_{{}_{L\perp}}\rho_{0}}\sqrt{\frac{m}{u}}=0.25$.

Our argument for this final state is based on our analytic steady state
solution of our equations of motion, which is illustrated in figure 4. This
solution, which is asymptotically valid sufficiently close to the analog in
our problem of the critical point in equilibrium phase separation, has all of
the features of the density profile of an equilibrium system in vapor-liquid
coexistence: on the left we have a low density “vapor” phase. As we move to
the right, the density increases smoothly through an interface with a well-
defined width, until it turns over and plateaus at a higher “liquid” density.
The approach to the liquid density as one moves to the right of the interface,
and the approach to the vapor density as one moves to the left, is
exponential.

It is this exponential tail - or, more precisely, the overlap of the
exponential tail of one high density region with that of one on the other side
of an intervening low density region- that leads to the attractive interaction
between bands. It is the fact that this interaction falls off exponentially
that makes the merger process so slow. Indeed, we expect that the time to
reach the fully phase separated state will, in this 1d picture, grow
exponentially with the system size. We believe this is the reason the system
cannot reach the final phase separated state in our numerical solutions.

However, we also believe that this exponential growth of equilibration time
with length scale is an artifact of the 1d nature of the analysis just
presented. While we believe this approach is accurate for early times - i.e.,
times up to and including the formation of the plateaus- we also believe that
once the plateaus are formed, their merger will be dominated by processes that
are missed in a 1d picture.

Specifically, we expect that what will really happen, once the plateaus form,
is that the bands will begin to undulate due to noise in the system, as
illustrated in figure 5. These undulations will grow with time until bands
begin bumping into their neighbors, at which point the bands can start to
merge. The merged region will then rapidly “zipper” along in both directions,
merging the two bands. This process will then repeat with the larger bands
formed as a result of such mergers, until full phase separation is achieved.

Figure 5: Our conjectured picture of the evolution of the density bands at
later times. The boundaries fluctuate due to noise in ways that are not
captured in the noiseless one dimensional approach. The central blue band is
depicted in our conjectured transient “zipper” state.

The time scale needed for undulations of the bands to grow large enough to
reach neighboring bands probably grows only algebraically with the separation
between bands - much as the interface fluctuations in the KPZ equation
doKardar  _et al._ (1986) \- rather than exponentially. We therefore expect
that the time required for our autochemotaxic instability to completely phase
separate a large system will grow only algebraically with system size.

However, to make this idea more precise, a theory of the dynamics of
fluctuations of an interface between two flocks of different densities and
speeds is needed. No such theory currently exists; its development remains a
topic for future research.

Those familiar with the dynamics of equilibrium phase separation will see the
similarities with our system, although ours differs from the equilibrium case
because of its strong anisotropy, and non-equilibrium effects.

One manifestation of these non-equilibrium effects is the determination of the
values of the two co-existing steady-state densities. In equilibrium, these
can be determined by the well-known “common tangent construction”Chaikin and
Lubensky (1995). Our aforementioned analytic steady-state solution shows that
in flocks that phase separate, the two coexisting densities are determined by
an uncommon tangent construction. This is very similar to the behavior found
in mobility induced phase separation (MIPS) in disordered active systems Cates
and Tailleur (2015).

More specifically, we find that as one approaches the dynamical critical
point, the uncommon tangent construction approaches a common tangent
construction - but only in this limit. In general, the uncommon tangent
construction applies.

Our analytic steady-state solution implies the phase diagram shown in figure 6
and 13, in which the vertical axis is a model parameter defined later which
decreases with increasing autochemotaxis , and the horizontal axis is the mean
density of the system. This analytic solution is asymptotically valid
sufficiently close to the critical point in this figure. In the region above
the curve labelled “Binodal line”, the uniform state is the only stable steady
state of the flock. Below the curve labelled “Spinodal line”, only the phase
separated state is stable. Between these two curves, either state can be
stable.

Figure 6: Phase diagram of flock phase-separation. Solid blue and dashed
orange lines are the binodal and spinodal lines respectively. The orange
filled region, labeled “unstable”, corresponds to the coexistence of the
‘liquid’ (high density) and ‘gas’ (low density) phases. See figure 13 and
section V for more details.

When the instability is weak, and when the system is close to the dynamical
analog of an equilibrium liquid-gas critical pointChaikin and Lubensky (1995),
we can analytically determine the two densities in terms of parameters of the
model. This analytic theory agrees extremely well with numerical solutions of
our equations of motion, even when the instability is not weak and the system
is fairly far from a critical point.

This said, our treatment of this phase separation is entirely “mean-field”:
that is, it ignores fluctuations in the local density and velocity.
Fluctuations in the density are well-knownChaikin and Lubensky (1995); Ma
(1995) to radically change the scaling behavior of the density near the
critical point in equilibrium systems. Furthermore, even stable flocks
experience enormous (indeed, divergent) fluctuation corrections to their
dynamicsToner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998); Toner
_et al._ (2005). Among other effects, these make the diffusion “constants” no
longer constant; they become functions of the wave vector $\mathbf{q}$ under
consideration, diverging as $\mathbf{q}\to{\bf 0}$. Since both types of
effects - i.e., critical and non-critical fluctuations- are present in our
system, we expect the scaling exponents for, e.g., the density difference
between the two phase separated regions, to be quantitatively different from
those we find here. Therefore, the effect of fluctuations on the non-
equilibrium phase separation we study here is a promising open topic for
future research.

Although we have focused in the above discussion on autochemotaxis, the
ultimate, longest length and time scale limit of our model, and the results we
obtain for it herein, also apply to any dry polar ordered active matter system
in which any microscopic mechanism drives the inverse compressibility
negative. There are many other mechanisms one could imagine that could do
this: direct attractive interactions between the active agents, for one, would
be sufficient. Any system in which any microscopic mechanism leads to a
negative inverse compressibility will be described, at the longest length and
time scales, by a Toner-Tu model with a negative inverse compressibility,
which is the basis for all of our predictions here. Therefore, , the
phenomenology we predict here - in particular, that the phase separation is
into bands of alternating higher and lower density form running parallel to
the direction of mean flock motion, as illustrated in figure 2\- will hold for
any mechanism of phase separation (e.g., sufficiently strong attractive
interactions between the boids).

The remainder of this paper is organized as follows: in section (II), we
review the Keller-Segel model of autochemotaxisKeller and Segel (1971a, 1970,
b); Horstmann (2003); Stark (2018)and the hydrodynamic (Toner-Tu) theory of
flockingToner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998); Toner
_et al._ (2005); Toner (2012a). In section III, we develop the equations of
motion for autochemotaxic flocks. In section IV, we linearize the equations of
motion, examine their mode structure, and identify the instability. In section
V, we analytically and numerically determine the final state to which the
instability involves, using a one-dimensional analysis for systems very near a
“critical point”, analogous to the critical point of phase separating
equilibrium systems, that we identify. This proves to be described by a common
tangent construction very similar to that for equilibrium systems. In section
VI, we demonstrate that further from the critical point, this phase separation
is described by an “uncommon tangent construction” quite similar to that found
in disordered active systems (e.g., motility-induced phase separation
(MIPS)Cates and Tailleur (2015). In section VII, we summarize our results, and
discuss possible directions for future research.

##  II Review of prior theories of (auto)chemotaxis and flocking

###  II.1 (Auto)chemotaxis

In the 1970s, Keller and Segel developed their model of chemotaxisKeller and
Segel (1971a, 1970, b); Horstmann (2003); Stark (2018). The general model
includes autochemotaxis, which describes a collection of creatures (which
we’ll call “boids”) emitting “chemo-attractant”; i.e., a substance to which
they themselves are attracted. The chemo-attractant then diffuses, and decays
with a finite lifetime. The boids diffuse as well, but with a bias in the
direction of the gradient of the chemo-attractant concentration. In addition,
Keller and Segel in their most general formulation allow for birth and death
of the boids, making their system what we would now describe as a “Malthusian
flock”Toner (2012b); Chen  _et al._ (2020).

All of these processes (diffusion, chemo-attractant generation and decay) will
in general depend on the local number density $\rho(\mathbf{r},t)$ of boids at
position $\mathbf{r}$ and time $t$, and the local concentration
$\eta(\mathbf{r},t)$ of the chemo-attractant at $(\mathbf{r},t)$. This
reasoning lead Keller and Segel to the following hydrodynamic equations in
their most general form:

$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle\nabla\cdot\bigg{(}k_{1}(\rho,\eta)\nabla\rho-
k_{2}(\rho,\eta)\rho\nabla\eta\bigg{)}+k_{3}(\rho,\eta)\,,$
$\displaystyle\partial_{t}\eta$ $\displaystyle=$ $\displaystyle
D_{\eta}\nabla^{2}\eta+k_{4}(\rho,\eta)-k_{5}(\rho,\eta)\eta\,.$ (II.2)

In this expression, $k_{1}$ is the motility of the boids, $k_{2}$ the
sensitivity to the chemical signal, $k_{3}$ the difference between the boid
birth and death rate, and $k_{4}$ and $k_{5}$ give respectively the production
and decay of the chemical signal. Hillen and Painter (2009)

When $k_{2}$ is positive, the chemical signal acts as a chemo-attractant. When
$k_{4}$ is positive and depends on $\rho$ (e.g., $k_{4}\propto\rho$), the
system is autochemotaxic. We will focus in this paper on flocks without birth
and death, which amounts to taking $k_{3}=0$. Otherwise, we will take the very
general form (II.2).

###  II.2 Flocking

Self-propelled entities in large numbers can display coherent motion over a
wide range of length scales. Inspired by the work of Viscek and the historical
progression of hydrodynamics, a hydrodynamic theory of flocking was developed
by Toner and Tu. Toner and Tu (1995); Tu  _et al._ (1998); Toner and Tu
(1998); Toner  _et al._ (2005); Toner (2012a).

The theory describes a flock composed of active, self-propelled particles
moving over a frictional substrate, so that momentum is not conserved and time
reversal symmetry is broken. In the jargon of the active matter field, this is
called a “dry” flock. The underlying dynamical rules are assumed to be
rotation and translation invariant. These symmetries, together with the
requirement that the total number of birds is conserved (that is, we ignore
birth and death “on the wing”) strongly constrain the form of the theory. In
particular, they determine which variables we need to keep in a long-
wavelength, long-time theory (which is what we mean by a “hydrodynamic”
theory).

In all hydrodynamic theories, it is sufficient to keep only the “slow”
variables; that is, the variables which evolve slowly at long length scales.
More precisely, this means the variables whose time scales of evolution
diverge as the length scale under consideration does.

There are only two generic reasons that a variable will be slow in this sense:
it must either be associated with a conserved quantity like boid number, or be
a “Goldstone mode” arising from a broken continuous symmetry. For a dry flock
as defined above, the only conserved quantity is boid number; the associated
hydrodynamic field is the coarse grained number density of boids
$\rho(\mathbf{r},t)$. Since an ordered flock spontaneously breaks the
continuous rotational symmetry of the underlying dynamics, the boid velocity
field $\mathbf{v}(\mathbf{r},t)$ (or, more precisely, its components
perpendicular to the mean velocity of the ordered state) is a Goldstone mode,
and hence, also slow.

The equations of motion for these fields keep all terms that are “relevant”,
in the renormalization group sense of affecting the scaling of the system at
long distances and times, and that are consistent with the symmetry (rotation
invariance) and conservation law (boid number) of the system. The requirement
of “relevance” is largely met by performing a gradient expansion. As shown by
Toner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998); Toner  _et al._
(2005), this forces the equations of motion to be of the form:

$\displaystyle\partial_{t}\mathbf{v}+\lambda_{1}(\mathbf{v}\cdot\nabla)\mathbf{v}+\lambda_{2}(\nabla\cdot\mathbf{v})\mathbf{v}+\lambda_{3}\nabla(|\mathbf{v}|^{2})=$
$\displaystyle\quad\quad U(|\mathbf{v}|,\rho)\mathbf{v}-\nabla
P_{1}(|\mathbf{v}|,\rho)-\mathbf{v}(\mathbf{v}\cdot\nabla
P_{2}(|\mathbf{v}|,\rho))+D_{B}\nabla(\nabla\cdot\mathbf{v})+D_{T}\nabla^{2}\mathbf{v}+D_{2}(\mathbf{v}\cdot\nabla)^{2}\mathbf{v}$
(II.3) $\displaystyle\partial_{t}\rho+\nabla\cdot(\mathbf{v}\rho)=0$ (II.4)

The significance of each of the terms in these equations is as follows:

The terms involving the parameters $\lambda_{i}$ are analogs of the convective
derivative of the coarse grained velocity field in the Navier-Stokes
equations. If our system respected Galilean invariance, we would have
$\lambda_{1}=1$ and $\lambda_{2,3}=0$. However, because our flock is on a
frictional substrate, which provides a special reference frame, we have
neither Galilean invariance nor momentum conservation.

The $U(|\mathbf{v}|,\rho)$ term, which is similar in form to a dissipative
term, clearly therefore also breaks both Galilean invariance and momentum
conservation. However, because our system is active, $U(|\mathbf{v}|,\rho)$
need not be negative for all $|\mathbf{v}|$. Indeed, if we are to model a
system in which the steady state is a moving flock, we must take it to have
the form illustrated in figure 7. In earlier literatureToner and Tu (1995); Tu
_et al._ (1998); Toner and Tu (1998); Toner  _et al._ (2005); Toner (2012b),
the special choice $U(|\mathbf{v}|,\rho)=\alpha-u|\mathbf{v}|^{2}$ is often
made. This is by no means necessary, however. Therefore, here we will make no
assumptions about the precise form of $U(|\mathbf{v}|,\rho)$, other than that
it is analytic in $|\mathbf{v}|$ and $\rho$, and that it has the qualitative
shape illustrated in (7).

Figure 7:  A plot of $U(|\mathbf{v}|,\rho)$ for fixed $\rho$. The effect of
this term is to maintain a non-zero speed, $v_{0}$, of the flock. Boids will
accelerate if their speed is below $v_{0}$ and decelerate if their speed is
larger than $v_{0}$.

The $P_{1}$ term is perfectly analogous to the isotropic pressure in the
Navier-Stokes equations. The $P_{2}$ term is an “anisotropic pressure”, and is
allowed because our system breaks rotation invariance locally, which means
that there is no reason that the response to density gradients along the local
velocity $\mathbf{v}$ should be the same as that to gradients perpendicular to
$\mathbf{v}$. Note that this term also breaks Galilean invariance.

The velocity diffusion constants $D_{B}$ and $D_{T}$ are precise analogs of
the bulk and shear viscosities, respectively, in the Navier-Stokes equations.
The diffusion constant $D_{2}$ is an anisotropic viscosity which has no analog
in the Navier-Stokes equation, because it violates Galilean invariance. Since
we lack Galilean invariance here, it is allowed in our problem. All of these
viscosities have the effect of suppressing fluctuations of the velocity away
from spatial uniformity.

The quantities $P_{1,2}(|\mathbf{v}|,\rho)$, $U(|\mathbf{v}|,\rho)$,
$D_{B,T,2}$, and $\lambda_{1,2,3}$ are in general functions of $|\mathbf{v}|$
and $\rho$. They can not depend of the direction of $\mathbf{v}$ due to
rotation invariance.

We are not interested in fluctuations here, so we have dropped a noise term
included by Toner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998);
Toner  _et al._ (2005); Toner (2012a) in the velocity equation of motion.

##  III The hydrodynamic theory of autochemotaxic flocks: Equations of motion

We now wish to develop a hydrodynamic theory of autochemotaxic flocking. This
will be a continuum model which takes as its variables the three fields
introduced above: the number density of boids $\rho(\mathbf{r},t)$, the boid
velocity field $\mathbf{v}(\mathbf{r},t)$, and the concentration
$\eta(\mathbf{r},t)$ of the chemical signal (or chemo-attractant/repellent).

Let’s begin with our modification of equation (II.2) for the dynamics of the
chemo-attractant, $\eta$, in the Keller-Segel model . First, our model is not
“Malthusian”; that is, we will not consider birth and death of the boids.
Therefore, we set $k_{3}=0$ in the equation of motion (LABEL:KSrho) for the
boid density.

In addition, because of the coherent motion of the flock, we need to
supplement the terms in the density equation of motion (LABEL:KSrho) with the
convective $\nabla\cdot(\mathbf{v}\rho)$ term in the Toner-Tu equation of
motion (II.4) for the boid density. A similar term, but with a non-unit
coefficient, can also appear in the chemo-attractant equation of motion.
Analogs of the “convective derivative” terms with coefficients
$\lambda_{1,2,3}$ that appear in the velocity equation of motion (II.3) are
allowed in the chemo-attractant equation of motion as well. We also need to
allow for anisotropic diffusive terms (i.e., anisotropic versions of the
$k_{1,2}$ terms) in the density equation, and similar terms in the chemo-
attractant equation of motion, in much the way we introduced the anisotropic
pressure $P_{2}$ in the velocity equation of motion (II.3) in the Toner-Tu
model.

In addition, all of the $k_{i}$’s can now depend on the local speed
$|\mathbf{v}|$, as well as on $\rho$ and $\eta$.

The coupling of the chemo-attractant concentration $\eta$ to the velocity
dynamics must, by symmetry, be embodied entirely by making the pressures
$P_{1,2}$, and the phenomenological coefficients $U$, $D_{B,T,2}$, and
$\lambda_{1,2,3}$, depend on $\eta$ as well as $|\mathbf{v}|$ and $\rho$. It
is the dependence of the isotropic pressure on $\eta$ that leads to the
instability we find, as we’ll show below.

This reasoning leads to the following equations of motion for the boid density
$\rho$, chemo-attractant concentration $\eta$, and the velocity
$\mathbf{v}$Tru :

velocity

$\displaystyle\partial_{t}\mathbf{v}+\lambda_{1}(\mathbf{v}\cdot\nabla)\mathbf{v}+\lambda_{2}(\nabla\cdot\mathbf{v})\mathbf{v}+\lambda_{3}\nabla(|\mathbf{v}|^{2})$
$\displaystyle=$ $\displaystyle U(|\mathbf{v}|,\rho,\eta)\mathbf{v}-\nabla
P_{1}(|\mathbf{v}|,\rho,\eta)-\mathbf{v}(\mathbf{v}\cdot\nabla
P_{2}(|\mathbf{v}|,\rho,\eta))$ (III.1) $\displaystyle+$ $\displaystyle
D_{B}\nabla(\nabla\cdot\mathbf{v})+D_{T}\nabla^{2}\mathbf{v}+D_{2}(\mathbf{v}\cdot\nabla)^{2}\mathbf{v}\,,$
$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle\nabla\cdot\bigg{(}k_{1}(|\mathbf{v}|,\rho,\eta)\nabla\rho-
k_{2}(|\mathbf{v}|,\rho,\eta)\rho\nabla\eta+k_{1a}(|\mathbf{v}|,\rho,\eta)\mathbf{v}(\mathbf{v}\cdot\nabla\rho)-k_{2a}(|\mathbf{v}|,\rho,\eta)\rho\mathbf{v}(\mathbf{v}\cdot\nabla\eta)\bigg{)}-\nabla\cdot(\mathbf{v}\rho)\,,$
$\displaystyle\partial_{t}\eta$ $\displaystyle=$ $\displaystyle
D_{\eta}\nabla^{2}\eta+k_{6}(|\mathbf{v}|,\rho,\eta)-\lambda_{\eta
1}(\mathbf{v}\cdot\nabla)\eta-\lambda_{\eta
2}(\nabla\cdot\mathbf{v})\eta{-\zeta_{\eta}\nabla\cdot(\mathbf{v}\rho)}+D_{\eta\rho}\nabla^{2}\rho+D_{\eta
v}\nabla^{2}|\mathbf{v}|$ (III.3) $\displaystyle+$ $\displaystyle D_{\eta
a}\nabla\cdot[\mathbf{v}(\mathbf{v}\cdot\nabla)\eta]+D_{\eta\rho
a}\nabla\cdot[\mathbf{v}(\mathbf{v}\cdot\nabla)\rho]+D_{\eta
va}\nabla\cdot[\mathbf{v}(\mathbf{v}\cdot\nabla)|\mathbf{v}|],$

where we’ve defined

$k_{6}(|\mathbf{v}|,\rho,\eta)\equiv
k_{4}(|\mathbf{v}|,\rho,\eta)-k_{5}(\rho,\eta)\eta\,.$ (III.4)

As just discussed, we connect the concentration of the auto-attractant to the
velocity equation by coupling them through the equations of state of the
pressures $P_{1,2}$.

We have also introduced anisotropic analogs $k_{1a}$ and $k_{2a}$ of the boid
mobility $k_{1}$ and response to the chemotactant $k_{2}$, since we expect the
response along the direction of the flock velocity $\mathbf{v}$ to be
different in general from the response perpendicular to that direction.

##  IV Derivation of the Chemo-Attractant Instability

###  IV.1 The homogeneous polar fluid state

We begin by identifying the steady state solution of our equations of motion
that corresponds to a homogeneous moving flock. In a homogeneous state, all of
our fields are constants, so we have

$\displaystyle\rho(\mathbf{r},t)$ $\displaystyle=\rho_{0}\,,$ (IV.1)
$\displaystyle\eta(\mathbf{r},t)$ $\displaystyle=\eta_{0}\,,$ (IV.2)
$\displaystyle\mathbf{v}$ $\displaystyle=v_{0}\hat{x}\,.$ (IV.3)

Note that the direction of $\mathbf{v}$ is completely arbitrary, due to the
rotation invariance of our model. We will henceforth choose our coordinate
system so that the $x$-axis is the direction of the spontaneous velocity.

Inserting these constant ansätze (IV.1), (IV.2) and (IV.3), into our equations
of motion (III.1), (III) and (III.3), it is clear that all terms involving
spatial or temporal derivatives vanish. It is easy to see that this implies
that the density equation (III) is automatically satisfied. The chemo-
attractant concentration equation in a steady, homogeneous state then reads

$k_{6}(\rho_{0},\eta_{0})=0\,,$ (IV.4)

while the velocity equation reduces to

$U(v_{0},\rho_{0},\eta_{0})=0\,.$ (IV.5)

This amounts to two scalar algebraic equations with three unknowns. We clearly
need one more condition. This can be obtained by fixing the number $N$ of
boids contained in the volume $V$ of our system (in spatial dimension $d=3$)
(or the area in $d=2$). Once we do so, we have, in a uniform system, fixed the
density $\rho_{0}$ everywhere to be $\rho_{0}=N/V$ or $\rho_{0}=N/{\rm Area}$.
Inserting that value of $\rho_{0}$ into (IV.4) then in principle determines
$\eta_{0}$. (We will assume (IV.4) has a unique solution for $\eta_{0}$ for
any value of $\rho_{0}$.)

With these values of $\rho_{0}$ and $\eta_{0}$ in hand, we can then in
principle use (IV.5) to determine the steady-state speed $v_{0}$. We will also
assume that the solution of (IV.5) is unique, which it clearly will be if
$U(v_{0},\rho_{0},\eta_{0})$ looks like figure 7.

###  IV.2 Linearization

Having found a steady state solution to our equations of motion that
corresponds to a spatially homogeneous and uniformly moving flock, we now will
determine the stability of this state. To do so, we begin by linearizing our
equations of motion about our steady state solution. That is, we will write

$\displaystyle\rho(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\rho_{0}+\delta\rho(\mathbf{r},t),$ (IV.6)
$\displaystyle\eta(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\eta_{0}+\delta\eta(\mathbf{r},t),$ (IV.7)
$\displaystyle\mathbf{v}(\mathbf{r},t)$ $\displaystyle=$ $\displaystyle
v_{0}\hat{x}+{\bf\delta v}(\mathbf{r},t)=(v_{0}+\delta
v_{x}(\mathbf{r},t))\hat{x}+\mathbf{v}_{\perp}(\mathbf{r},t),$

and then expand our equations of motion to linear order in the fluctuations
$\delta\rho(\mathbf{r},t)$, $\delta\eta(\mathbf{r},t)$, and ${\bf\delta
v}(\mathbf{r},t)$ of the fields $\rho(\mathbf{r},t)$, $\eta(\mathbf{r},t)$,
and $\mathbf{v}(\mathbf{r},t)$.

The process of linearization begins by expanding all of the $|\mathbf{v}|$,
$\rho$, and $\eta$ dependent parameters in the equations of motion to
sufficiently high order in the fluctuations $\delta\rho(\mathbf{r},t)$,
$\delta\eta(\mathbf{r},t)$, and ${\bf\delta v}(\mathbf{r},t)$ of the fields
$\rho(\mathbf{r},t)$, $\eta(\mathbf{r},t)$, and $\mathbf{v}(\mathbf{r},t)$ to
obtain all terms in the equation of motion to linear order in
$\delta\rho(\mathbf{r},t)$, $\delta\eta(\mathbf{r},t)$, and ${\bf\delta
v}(\mathbf{r},t)$. This gives us

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle D_{{}_{C}}\,,$ (IV.9)
$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle D_{{}_{C}}\mu\,,$
(IV.10) $\displaystyle k_{3}$ $\displaystyle=$ $\displaystyle 0\,,$ (IV.11)
$\displaystyle k_{6}$ $\displaystyle=$
$\displaystyle\gamma\delta\rho-{\delta\eta\over\tau}+\alpha\delta|\mathbf{v}|\,,$
(IV.12) $\displaystyle P_{1,2}(|\mathbf{v}|,\rho,\eta)$ $\displaystyle=$
$\displaystyle{\rm
constant}+\sigma_{1,2}\delta\rho-\kappa_{1,2}\delta\eta+\nu_{1,2}\delta|\mathbf{v}|\,,$
$\displaystyle U(|\mathbf{v}|,\rho,\eta)$ $\displaystyle=$
$\displaystyle-A\delta|\mathbf{v}|+{{\bar{B}}}\delta\rho+C\delta\eta\,,$
(IV.14)

where $D_{{}_{C}}$, $\mu$, $\gamma$, $\tau$, $\sigma_{1,2}$, $\kappa_{1,2}$,
$A$, ${{\bar{B}}}$, and $C$ are all phenomenological constants. We remind the
reader that boids are not being born and dying “on the wing”. This means that
the boid birth and death rate term $k_{3}$ is zero. Most of the other
parameters in the model can have either sign, but $\kappa_{1}$ and $A$ are
both positive definite. The positivity of $A$ can be seen from figure 7 ,
since $A$ is just $-1$ times the slope of $U(v)$ as it crosses the $v$ axis.
The positivity of $\kappa_{1}$ means that the chemical signal attracts the
boids ( i.e., that it is a chemo-attractant), which is the case we wish to
consider here. As we will see, it is this attraction that destabilizes the
spatially uniform flocking state.

Note that we also have ${{\bar{B}}}>0$ for systems in which the mean speed
$v_{0}$ increases with increasing density, as it does for, e.g., the Vicsek
modelVicsek  _et al._ (1995).

The parameters $\gamma$ and $\tau$ respectively set the sensitivity of the
production rate and decay rate of the chemo-attractant.

Making the substitutions (LABEL:field_expand) and (IV.9-IV.14) in our
equations of motion (III.1-III.3), and keeping only terms to linear order in
$\delta\rho$, $\delta\eta$, and $\delta\mathbf{v}$, we obtain our linearized
equations of motion:

$\displaystyle\frac{\partial\delta\rho}{\partial t}$ $\displaystyle=$
$\displaystyle D_{{}_{C}}\nabla^{2}\delta\rho-\mu
D_{{}_{C}}(\rho_{0}\nabla^{2}\delta\eta)-\rho_{0}\partial_{x}\delta
v_{x}-\rho_{0}\nabla_{\perp}\cdot\mathbf{v}_{\perp}-v_{0}\partial_{x}\delta\rho+D_{ca}v_{0}^{2}\partial_{x}^{2}\delta\rho-
D_{ca}\mu_{a}v_{0}^{2}\partial_{x}^{2}\delta\eta\,,$ (IV.15)
$\displaystyle\frac{\partial\delta\eta}{\partial t}$ $\displaystyle=$
$\displaystyle
D_{\eta}\nabla^{2}\delta\eta-\frac{\delta\eta}{\tau}+\gamma\delta\rho+\alpha\delta
v_{x}{-v_{0}\zeta_{\eta}\partial_{x}\delta\rho}-\lambda_{\eta
1}v_{0}\partial_{x}\delta\eta-{(\lambda_{\eta
2}\eta_{0}+\zeta_{\eta}\rho_{0})}[\partial_{x}\delta
v_{x}+\nabla_{\perp}\cdot\mathbf{v}_{\perp}]$ (IV.16) $\displaystyle+$
$\displaystyle D_{\eta\rho}\nabla^{2}\delta\rho+D_{\eta v}\nabla^{2}\delta
v_{x}+D_{\eta\rho a}v_{0}^{2}\partial_{x}^{2}\delta\rho+D_{\eta
a}v_{0}^{2}\partial_{x}^{2}\delta\eta+D_{\eta
va}v_{0}^{2}\partial_{x}^{2}\delta v_{x}\,,$
$\displaystyle\frac{\partial\delta v_{x}}{\partial t}$ $\displaystyle=$
$\displaystyle v_{0}(-A\delta
v_{x}+{{\bar{B}}}\delta\rho+C\delta\eta)-(\lambda_{1}+\lambda_{2}+2\lambda_{3})v_{0}\partial_{x}\delta
v_{x}-\lambda_{2}v_{0}(\nabla_{\perp}\cdot\mathbf{v}_{\perp})-(\sigma_{1}+v_{0}^{2}\sigma_{2})\partial_{x}\delta\rho$
(IV.17) $\displaystyle+$
$\displaystyle(\kappa_{1}+v_{0}^{2}\kappa_{2})\partial_{x}\delta\eta-(\nu_{1}+v_{0}^{2}\nu_{2})\partial_{x}\delta
v_{x}+D_{B}[\partial_{x}(\nabla_{\perp}\cdot\mathbf{v}_{\perp})+\partial_{x}^{2}\delta
v_{x}]+D_{T}\nabla^{2}\delta v_{x}+D_{2}v_{0}^{2}\partial_{x}^{2}\delta
v_{x}\,,$ $\displaystyle\frac{\partial\mathbf{v}_{\perp}}{\partial t}$
$\displaystyle=$
$\displaystyle-\lambda_{1}v_{0}\partial_{x}\mathbf{v}_{\perp}-2\lambda_{3}v_{0}\nabla_{\perp}\delta
v_{x}-\nabla_{\perp}(\sigma_{1}\delta\rho-\kappa_{1}\delta\eta+\nu_{1}\delta
v_{x})+D_{B}\nabla_{\perp}[\nabla_{\perp}\cdot\mathbf{v}_{\perp}+\partial_{x}\delta
v_{x}]$ (IV.18)
$\displaystyle\,\,+D_{T}\nabla^{2}_{\perp}\mathbf{v}_{\perp}+D_{x}\partial_{x}^{2}\mathbf{v}_{\perp}\,.$

where $D_{x}\equiv D_{T}+D_{2}v_{0}^{2}$.

We’ll now use the equations to investigate the stability of the uniform steady
state.

###  IV.3 Demonstration that $\delta\eta$ and $\delta v_{x}$ are Fast
Variables, and their elimination

We begin by observing that the fields $\delta v_{x}$ and $\delta\eta$ are
“fast” in the sense that their time derivatives do not vanish in the limit of
extremely slowly spatially varying fields. In contrast, the fields $\rho$ and
$\mathbf{v}_{\perp}$ are “slow”, or, to use another word, “hydrodynamic”: that
is, their time derivatives do vanish in the limit of extremely slowly
spatially varying fields.

As a result, if we are considering the most slowly evolving modes in the
system at the longest wavelengths, the “fast” fields $\delta v_{x}$ and
$\delta\eta$ will, on the time scale of those slowest modes, relax very
quickly back to values determined entirely by the slow fields. This means that
those fast variables become effectively “enslaved” to the hydrodynamic
variables in the problem - in this case the boid density $\rho$. By
“enslaved”, we mean that its value at any point in space and time is
determined entirely by the instantaneous value of $\rho$ at the same point in
space and time. It is this sort of “enslavement” that enables hydrodynamic
theories in general to eliminate “fast” variables - that is, variables that do
not relax infinitely slowly as the length scale under consideration goes to
infinity- and work only with the “slow” hydrodynamic variables.

To see how this works in our problem, note that, in the hydrodynamic limit of
long wavelengths and large time scales, derivatives of the fluctuations of the
fields are small compared to terms that contain the fluctuation with no
derivatives. We will therefore expand the equations of motion for the fast
fields $\eta$ and $\delta v_{x}$ in derivatives by writing

$\delta\eta=\delta\eta_{0}+\delta\eta_{1}\,\,\,\,\,\,,\,\,\,\,\,\,\delta
v_{x}=\delta v_{x}^{(0)}+\delta v_{x}^{(1)}\,,$ (IV.19)

where

$\delta\eta_{0}\,,\delta
v_{x}^{(0)}=O(\delta\rho)\,\,,\,\,\delta\eta_{1}\,,\delta
v_{x}^{(1)}=O(\partial_{t}\delta\rho,\nabla\delta\rho,\nabla\cdot\mathbf{v}_{\perp})\,.$
(IV.20)

Note that $\nabla\cdot v_{\perp}\ll\delta\rho$, as can be seen by the
following argument: In order to satisfy Eq (IV.15), $\nabla\cdot v_{\perp}$
must be of order $(\partial_{t}\rho)$ or $\partial_{x}\delta\rho$ (the only
other term in (IV.15), namely $\partial_{x}\delta v_{x}$, is clearly much
smaller than $\nabla\cdot v_{\perp}$, since $x$-derivatives are much smaller
than $\perp$ derivatives in the unstable regime (where
$q_{x}\ll\mathbf{q}_{\perp}$), and $\delta v_{x}$ is $\ll\mathbf{v}_{\perp}$,
since the former is massive). Therefore, $\nabla\cdot v_{\perp}$ is of order a
derivative of $\delta\rho$, and hence much less than $\delta\rho$, as we
asserted.

Therefore, to leading order in a hydrodynamic expansion, we can drop all terms
with derivatives from the equations of motion (IV.16) and (IV.17). Dropping
those terms, equations (IV.16) and (IV.17) become:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle-\frac{\delta\eta_{0}}{\tau}+\gamma\delta\rho+\alpha\delta
v_{x}^{(0)}\,,$ (IV.21) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
v_{0}(-A\delta v_{x}^{(0)}+{{\bar{B}}}\delta\rho+C\delta\eta_{0})\,.$ (IV.22)

These equations can trivially be solved to relate the “fast” variables $\delta
v_{x}$ and $\delta\eta$ to the “slow” variable $\delta\rho$:

$\displaystyle\delta v_{x}^{(0)}$ $\displaystyle=$
$\displaystyle\left(\frac{{{\bar{B}}}+\tau\gamma C}{A-\alpha\tau
C}\right)\delta\rho=K_{v}\delta\rho,\,,$ (IV.23) $\displaystyle\delta\eta_{0}$
$\displaystyle=$
$\displaystyle\left(\frac{(A\gamma+\alpha{{\bar{B}}})\tau}{A-\alpha\tau
C}\right)\delta\rho=K_{\eta}\delta\rho\,,$ (IV.24)

where we’ve defined

$K_{v}\equiv\frac{{{\bar{B}}}+\tau\gamma C}{A-\alpha\tau
C}\,\,\,\,\,\,,\,\,\,\,\,\,K_{\eta}\equiv\frac{(A\gamma+\alpha{{\bar{B}}})\tau}{A-\alpha\tau
C}\,.$ (IV.25)

The only important point to be noted about these rather complicated
expressions is that both $K_{v}$ and $K_{\eta}$ are monotonically increasing
functions of the chemo-attractant release rate $\gamma$, provided $A$ and $C$
are positive. As pointed out earlier, $A$ definitely is positive. Strictly
speaking, $C$ could, in principle, have either sign. However, we expect
intuitively that increasing the chemo-attractant concentration should increase
the speed of the boids (since most critters move faster when they smell
food!). If this is the case, then $C$ will be positive as well.

We’ll show later that the fact that both $K_{v}$ and $K_{\eta}$ are
monotonically increasing functions of the chemo-attractant release rate
$\gamma$ implies that the system will always eventually become unstable if we
increase the chemo-attractant release rate $\gamma$ without bound.

The alert reader will notice that $K_{v}$ and $K_{\eta}$ both diverge as
$A\to\alpha\tau C$ from above. The physics of this instability is simply that,
when the denominator vanishes, the decay rate of one of the “fast” modes
vanishes. This can be seen by noting that the denominator in question is
nothing but the determinant of the dynamical matrix for the fast fields if we
set the slow field $\rho=0$. That determinant is simply the product of the two
eigen-decay rates. While this instability would clearly be an interesting
topic for future research, we will not discuss it further here. Instead, we
will focus our attention on the parameter regime $A>\alpha\tau C$, where this
particular instability does not occur.

Inserting (IV.23) and (IV.24) into the linearized equation of motion (IV.17)
for $\delta v_{x}$, and gathering terms proportional to one derivative (either
spatial or temporal) of the hydrodynamic variables $\delta\rho$ and
$\mathbf{v}_{\perp}$, gives

$v_{0}(-A\delta
v_{x}^{(1)}+C\delta\eta_{1})=K_{v}\partial_{t}\delta\rho+G_{v}\partial_{x}\delta\rho+\lambda_{2}v_{0}\nabla\cdot\mathbf{v}_{\perp}\,,$
(IV.26)

where

$\displaystyle G_{v}$ $\displaystyle\equiv$
$\displaystyle\sigma_{1}+v_{0}^{2}\sigma_{2}+K_{v}(\nu_{1}+\nu_{2}v_{0}^{2}+(\lambda_{1}+\lambda_{2}+2\lambda_{3})v_{0})$
(IV.27) $\displaystyle\,-K_{\eta}(\kappa_{1}+v_{0}^{2}\kappa_{2})\,.$

Likewise, inserting (IV.23) and (IV.24) into the linearized equation of motion
(IV.15) for $\delta\eta$, and gathering terms proportional to one derivative
(either spatial or temporal) of the hydrodynamic variables $\delta\rho$ and
$\mathbf{v}_{\perp}$, gives

$-\frac{\delta\eta_{1}}{\tau}+\alpha\delta
v_{x}^{(1)}=K_{\eta}\partial_{t}\delta\rho+{(\lambda_{\eta
2}\eta_{0}+\zeta_{\eta}\rho_{0})}\nabla\cdot\mathbf{v}_{\perp}+G_{\eta}\partial_{x}\delta\rho\,,$
(IV.28)

where

$G_{\eta}\equiv\lambda_{\eta 1}v_{0}K_{\eta}+\lambda_{\eta
2}\eta_{0}K_{v}{+v_{0}\zeta_{\eta}}\,.$ (IV.29)

Solving these two simple linear equations (IV.26) and (IV.28) for $\delta
v_{x}^{(1)}$ and $\delta\eta_{1}$ gives

$\delta v_{x}^{(1)}=-K_{v1}\partial_{t}\delta\rho-
G_{v1}\partial_{x}\delta\rho-H_{v}\nabla\cdot\mathbf{v}_{\perp}$ (IV.30)

where

$K_{v1}\equiv{v_{0}C\tau K_{\eta}+K_{v}\over v_{0}(A-\alpha C\tau)}\,,$
(IV.31)

$G_{v1}\equiv{G_{v}+C\tau v_{0}G_{\eta}\over v_{0}(A-\alpha C\tau)}\,,$
(IV.32)

$H_{v}\equiv{\lambda_{2}+{(\lambda_{\eta
2}\eta_{0}+\zeta_{\eta}\rho_{0})}C\tau\over A-\alpha C\tau}\,,$ (IV.33)

and

$\delta\eta_{1}=-K_{\eta 1}\partial_{t}\delta\rho-G_{\eta
1}\partial_{x}\delta\rho-H_{\eta}\nabla\cdot\mathbf{v}_{\perp}\,,$ (IV.34)

where

$K_{\eta 1}\equiv{v_{0}A\tau K_{\eta}+\alpha\tau K_{v}\over v_{0}(A-\alpha
C\tau)}\,,$ (IV.35) $G_{\eta 1}\equiv{\tau(\alpha G_{v}+Av_{0}G_{\eta})\over
v_{0}(A-\alpha C\tau)}\,,$ (IV.36)
$H_{\eta}\equiv{\tau(\alpha\lambda_{2}+{(\lambda_{\eta
2}\eta_{0}+\zeta_{\eta}\rho_{0})}A)\over A-\alpha C\tau}\,.$ (IV.37)

Where $K_{v1}$ and $K_{\eta 1}$ give the linear response of $\delta
v_{x}^{(1)}$ and $\delta\eta_{1}$ to the time rate of change of the local
density, $\delta\rho$. Likewise, the $G_{v1}$ and $G_{\eta 1}$ give the linear
response to the spatial change in the direction of flock motion to the local
density, $\delta\rho$. Finally, $H_{\nu}$ and $H_{\eta}$ give the linear
response to the divergence of the local velocity that is perpendicular to the
direction of flock motion, $\Delta\cdot\mathbf{v}_{\perp}$.

Having eliminated the fast variables $\eta$ and $\delta v_{x}$ in favor of the
slow variable $\delta\rho$, we can now write a set of closed equations of
motion for the density fluctuation $\delta\rho$ and perpendicular velocity
field $\mathbf{v}_{\perp}$ by substituting equations (IV.23), (IV.24), (IV.30)
and (IV.34) into our equations of motion (IV.15) and (IV.18). This gives:

$\displaystyle\frac{\partial\delta\rho}{\partial t}$ $\displaystyle=$
$\displaystyle D_{\rho\perp}\nabla_{\perp}^{2}\delta\rho+D_{\rho
x}\partial_{x}^{2}\delta\rho-
v_{\rho}\partial_{x}\delta\rho-\rho_{0}\nabla_{\perp}\cdot\mathbf{v}_{\perp}+D_{\rho
v}\partial_{x}(\nabla_{\perp}\cdot\mathbf{v}_{\perp})+\phi\partial_{t}\partial_{x}\delta\rho\,,$
(IV.38) $\displaystyle\frac{\partial\mathbf{v}_{\perp}}{\partial t}$
$\displaystyle=$ $\displaystyle-
v_{v}\partial_{x}\mathbf{v}_{\perp}-\left({B\over\rho_{0}}\right)\nabla_{\perp}\delta\rho+D_{BR}\nabla_{\perp}(\nabla\cdot\mathbf{v}_{\perp})+D_{T}\nabla^{2}_{\perp}\mathbf{v}_{\perp}+D_{x}\partial_{x}^{2}\mathbf{v}_{\perp}+\nu_{x}\partial_{x}\nabla_{\perp}\delta\rho+\nu_{t}\partial_{t}\nabla_{\perp}\delta\rho\,,$

where we have defined:

$\displaystyle v_{\rho}$ $\displaystyle\equiv$
$\displaystyle\rho_{0}K_{v}+v_{0}\,,$ (IV.40) $\displaystyle D_{\rho\perp}$
$\displaystyle\equiv$ $\displaystyle D_{{}_{C}}(1-\mu K_{\eta}\rho_{0})\,,$
(IV.41) $\displaystyle D_{\rho x}$ $\displaystyle\equiv$ $\displaystyle
D_{{}_{C}}+D_{Ca}v_{0}^{2}-D_{{}_{C}}\mu\rho_{0}K_{\eta}$ (IV.42)
$\displaystyle\,\,-D_{Ca}\mu_{a}v_{0}^{2}\rho_{0}K_{\eta}+\rho_{0}G_{v1}\,,$
$\displaystyle D_{\rho v}$ $\displaystyle\equiv$
$\displaystyle\rho_{0}H_{v}\,,$ (IV.43) $\displaystyle\phi$
$\displaystyle\equiv$ $\displaystyle\rho_{0}K_{v1}\,,$ (IV.44) $\displaystyle
v_{v}$ $\displaystyle\equiv$ $\displaystyle\lambda_{1}v_{0}\,,$ (IV.45)
$\displaystyle D_{BR}$ $\displaystyle\equiv$ $\displaystyle
D_{B}+(2\lambda_{3}v_{0}+\nu_{1})H_{v}+\kappa_{1}H_{\eta}\,,$ (IV.46)
$\displaystyle\nu_{x}$ $\displaystyle\equiv$ $\displaystyle\kappa_{1}G_{\eta
1}+D_{B}K_{v}+(2\lambda_{3}v_{0}+\nu_{1})G_{v}\,,$ (IV.47)
$\displaystyle\nu_{t}$ $\displaystyle\equiv$ $\displaystyle\kappa_{1}K_{\eta
1}+(2\lambda_{3}v_{0}+\nu_{1})K_{v1}\,,$ (IV.48)

and, most importantly, the “inverse compressibility”

$\displaystyle B$ $\displaystyle\equiv$
$\displaystyle\rho_{0}[\sigma_{1}-K_{\eta}\kappa_{1}+(2\lambda_{3}v_{0}+\nu_{1})K_{v}]$
$\displaystyle=$
$\displaystyle\rho_{0}\bigg{(}\sigma_{1}-\left[\frac{\kappa_{1}(A\gamma+\alpha{{\bar{B}}})\tau{-}(2\lambda_{3}v_{0}+\nu_{1})({{\bar{B}}}+\gamma\tau
C)}{A-\alpha\tau C}\right]\bigg{)}\,.$

This last expression is the most important result of this analysis, since it
shows that autochemotaxis can make the inverse compressibility negative. One
way this can be done is by increasing the parameter $\kappa_{1}$, which
governs the sensitivity of the isotropic pressure to the chemo-attractant
concentration. As can been by inspection of (LABEL:Bfinal), increasing
$\kappa_{1}$ without bound while holding all other parameters fixed will
inevitably make $B$ become negative, provided only that

${(A\gamma+\alpha{{\bar{B}}})\tau\over A-\alpha\tau C}>0\,,$ (IV.50)

which clearly is a very large portion of the parameter space of our
hydrodynamic model.

Likewise, increasing the chemo-attractant production rate $\gamma$ with all
other parameters held fixed will also unavoidably drive the inverse
compressibility $B$ negative, provided that

${(\kappa_{1}A+(2\lambda_{3}v_{0}+\nu_{1})C)\tau\over A-\alpha\tau C}>0\,.$
(IV.51)

Finally, increasing the chemo-attractant lifetime $\tau$ with all other
parameters held fixed will also certainly drive the inverse compressibility
$B$ negative if, for example, all of the parameters $\alpha$, ${{\bar{B}}}$,
$C$, and the combination $2\lambda_{3}v_{0}+\nu_{1}$ are all positive. Indeed,
in this case, $B$ is guaranteed to turn negative as $\tau$ is increased from
zero for some $\tau<{A\over\alpha C}$.

The forms of the equations of motion (IV.38) and (LABEL:vperpEq2) are
identical to the linearized equations of motion obtained in Toner (2012a) for
stable flocks without chemotaxis. This was inevitable, since the symmetries
and conservation laws of our model are the same as those of an ordinary,
stable flock. The only long-wavelength effect of the autochemotaxis is to
provide a mechanism to drive the inverse compressibility $B$ negative, as
discussed above. Any attractive interaction, which is sufficiently strong to
drive the inverse compressibility negative, will obey the same hydrodynamic
description put forward in this paper. In the next section, we will show that
when $B$ does become negative, the homogeneous phase of the flock becomes
unstable.

###  IV.4 Mode Structure

To investigate the mode structure of our linearized model, we proceed exactly
as in Toner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998); Toner
_et al._ (2005); Toner (2012a), by seeking solutions of the linearized
equations of motion (IV.38) and (LABEL:vperpEq2) for plane-wave normal modes,
in which all fields are proportional to $\exp{[i\bm{q}\cdot\bm{r}-i\omega
t]}$. Doing so, the equations of motion then read:

$\displaystyle\bigg{[}-i(\omega-
v_{\rho}q_{x})+D_{\rho\perp}q_{\perp}^{2}+D_{\rho x}q_{x}^{2}-\phi\omega
q_{x}\bigg{]}\delta\rho$
$\displaystyle=-i\rho_{0}\mathbf{q}_{\perp}\cdot\mathbf{v}_{\perp}-D_{\rho
v}q_{x}(\mathbf{q}_{\perp}\cdot\mathbf{v}_{\perp})\,,$ (IV.52)
$\displaystyle\bigg{[}-i\omega+iv_{v}q_{x}+D_{T}q^{2}_{\perp}+D_{x}q^{2}_{x}\bigg{]}\mathbf{v}_{\perp}+D_{BR}\mathbf{q}_{\perp}(\mathbf{q}_{\perp}\cdot\mathbf{v}_{\perp})$
$\displaystyle=\bigg{[}-i\frac{B}{\rho_{0}}+\nu_{t}\omega-\nu_{x}q_{x}\bigg{]}\mathbf{q}_{\perp}\delta\rho\,.$
(IV.53)

We decouple these by projecting (IV.53) perpendicular to and along
$\mathbf{q}_{\perp}$. That is, we write

$\mathbf{v}_{\perp}=v_{L}\hat{q}_{\perp}+\mathbf{v}_{T}\,,$ (IV.54)

with the “transverse” components $\mathbf{v}_{T}$, by definition,
perpendicular to $\mathbf{q}_{\perp}$ (that is,
$\mathbf{v}_{T}\cdot\hat{q}_{\perp}=0$), while the single “longitudinal”
component
$v_{L}=\frac{\mathbf{q}_{\perp}\cdot\mathbf{v}_{\perp}}{|q_{\perp}|}$ is the
projection of $\mathbf{v}_{\perp}$ onto $\mathbf{q}_{\perp}$.

This split decouples $\delta\rho$ and $v_{L}$ from the $d-2$ transverse modes
$\mathbf{v}_{T}$, as can be seen by inserting the decomposition (IV.54) into
the $\rho$ equation of motion (IV.52), and projecting the velocity equation of
motion (IV.53) in the transverse and longitudinal directions (i.e., along and
perpendicular to $\mathbf{q}_{\perp}$). This gives:

$\displaystyle\bigg{[}-i(\omega-
v_{\rho}q_{x})+\Gamma_{\rho}(\mathbf{q})-\phi\omega q_{x}\bigg{]}\delta\rho$
$\displaystyle=-i\rho_{0}q_{\perp}v_{L}-D_{\rho v}q_{x}q_{\perp}v_{L},$
(IV.55) $\displaystyle\bigg{[}-i(\omega-
v_{v}q_{x})+\Gamma_{L}(\mathbf{q})\bigg{]}v_{L}$
$\displaystyle=\bigg{[}-i\frac{B}{\rho_{0}}+\nu_{t}\omega-\nu_{x}q_{x}\bigg{]}q_{\perp}\delta\rho,$
(IV.56) $\displaystyle\bigg{[}-i(\omega-
v_{v}q_{x})+\Gamma_{T}(\mathbf{q})\bigg{]}\mathbf{v}_{T}$ $\displaystyle={\bf
0}\,.$ (IV.57)

Here we’ve defined $\Gamma_{T}=D_{T}q^{2}_{\perp}+D_{x}q_{x}^{2}$,
$\Gamma_{L}=D_{{}_{L\perp}}q^{2}_{\perp}+D_{x}q_{x}^{2}$,
$D_{{}_{L\perp}}=D_{T}+D_{BR}$, and
$\Gamma_{\rho}=D_{\rho\perp}q_{\perp}^{2}+D_{\rho x}q_{x}^{2}$.

Altogether we have $d$ frequencies to determine. Since the transverse modes
are decoupled, we can immediately write down the transverse dispersion
relation implied by equation (IV.57):

$\displaystyle\omega=-i\Gamma_{T}+v_{v}q_{x}$ (IV.58)

All these $d-2$ modes decay, provided $D_{T}q_{\perp}^{2}+D_{x}q_{x}^{2}>0$.
Thus, there is no instability in these modes, even when the inverse
compressibility $B$ becomes negative (unsurprisingly, since $B$ does not enter
the transverse equations of motion).

The instability arises from the coupled density $\rho$ and longitudinal
velocity $v_{L}$ modes, as can be seen from equations (IV.55) and (IV.56),
which can be rewritten:

$\displaystyle-i\omega\delta\rho$ $\displaystyle=$ $\displaystyle-
D_{\rho\perp}q_{\perp}^{2}\delta\rho-D_{\rho x}q_{x}^{2}\delta\rho-
iv_{\rho}q_{x}\delta\rho-i\rho_{0}q_{\perp}v_{L}-D_{\rho
v}q_{x}q_{\perp}v_{L}+\phi\omega q_{x}\delta\rho\,,$ (IV.59)
$\displaystyle-i\omega v_{L}$ $\displaystyle=$ $\displaystyle-
iv_{v}q_{x}v_{L}-i\left({B\over\rho_{0}}\right)q_{\perp}\delta\rho-
D_{{}_{L\perp}}q^{2}_{\perp}v_{L}-D_{x}q^{2}_{x}v_{L}-\nu_{x}q_{x}q_{\perp}\delta\rho+\nu_{t}\omega
q_{\perp}\delta\rho\,\,.$ (IV.60)

As shown in Toner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998);
Toner  _et al._ (2005); Toner (2012a), the eigenfrequencies of these
equations, to linear order in $q$, are

$\displaystyle\omega_{\pm}=\frac{1}{2}\bigg{[}q_{x}v_{\rho}+q_{x}v_{v}\pm\Delta\bigg{]},$
(IV.61)
$\displaystyle\Delta=\sqrt{(v_{v}-v_{\rho})^{2}q_{x}^{2}+4Bq_{\perp}^{2}}.$
(IV.62)

In contrast to prior work, however, we now know that autochemotaxis can drive
the inverse compressibility $B$ negative. When this happens, it is clear from
(IV.62) that modes whose wavevectors $\mathbf{q}$ lie in the region in which

$\displaystyle(v_{v}-v_{\rho})^{2}q_{x}^{2}+4Bq_{\perp}^{2}<0\,,$ (IV.63)

will be unstable, thereby destabilizing the spatially uniform flocking state.
This implies the instability lies in the region

$|q_{x}|<{2\sqrt{|B|}\over|v_{v}-v_{\rho}|}|\mathbf{q}_{\perp}|=\epsilon|\mathbf{q}_{\perp}|\,,$
(IV.64)

where we’ve defined:

$\epsilon\equiv{2\sqrt{|B|}\over|v_{v}-v_{\rho}|}\,.$ (IV.65)

(Recall that in the unstable region $B<0$). If we are near the instability on
the unstable side, $|B|\ll(v_{v}-v_{\rho})^{2}$, and so (IV.64) implies
$\epsilon<<1$, and so the unstable wavevectors satisfy
$|q_{x}|\ll|\mathbf{q}_{\perp}|$, as claimed in the introduction.

The window (IV.64) gives the region in which, at sufficiently small
$\mathbf{q}$, the flock becomes unstable with a growth rate that is linear in
$q\equiv|\mathbf{q}|$. Small $\mathbf{q}$’s that lie outside this region,
however, are not necessarily stable. This is because, for
$|q_{x}|>\epsilon|\mathbf{q}_{\perp}|$, the eigenfrequencies $\omega$ to
linear order in $q$, are real, as can be seen by inspection of equation IV.61.
Since the stability or instability is determined entirely by the imaginary
part of $\omega$, we must therefore, for $\mathbf{q}$’s in the regime
$|q_{x}|>\epsilon|\mathbf{q}_{\perp}|$, work to higher (in fact, second) order
in $\mathbf{q}$ in order to determine the absolute instability limit. This
proves, as we shall see, to be have a larger opening angle at small
$\mathbf{q}$ than the limit (IV.64), but the new limit has the same
$\sqrt{|B|}$ scaling with $B$. In contrast to equation (IV.64), however, the
instability is only limited to a narrow window if $|v_{v}-v_{\rho}|$ exceeds
some finite threshold. See equation (IV.77) for more details.

It is instructive to rewrite the dispersion relation (IV.62) in terms of two
complex, direction dependent sound speeds $c_{\pm}(\theta_{\mathbf{q}})$:

$\omega=c_{\pm}(\theta_{\mathbf{q}})q\,,$ (IV.66)

where $q\equiv|\mathbf{q}|$, and $\theta_{\mathbf{q}}$ is the angle between
the $\mathbf{q}$ and the $\hat{x}$ direction (i.e., the direction of mean
flock motion $\langle\mathbf{v}\rangle$). Using $q_{x}=q\cos{\theta}$ and
$q_{\perp}=q\sin{\theta}$ in (IV.62), we obtain the expression of references
Toner and Tu (1995); Tu  _et al._ (1998); Toner and Tu (1998); Toner  _et al._
(2005); Toner (2012a) for $c_{\pm}(\theta_{\mathbf{q}})$:

$\displaystyle
c_{\pm}(\theta_{\mathbf{q}})=\frac{1}{2}\bigg{[}(v_{\rho}+v_{v})\cos(\theta_{\mathbf{q}})\pm\sqrt{(v_{v}-v_{\rho})^{2}\cos^{2}(\theta_{\mathbf{q}})+4B\sin^{2}(\theta_{\mathbf{q}})}\bigg{]}.$
(IV.67)

from this expression, we can see that $\omega_{+}$ acquires a positive
imaginary part, signaling instability, for directions of propagation
$\theta_{\mathbf{q}}$ near $\pm\pi/2$; i.e., perpendicular to the direction of
mean flock motion. Specifically, the instability occurs for

$|\theta\pm\pi/2|<{2\sqrt{|B|}\over|v_{v}-v_{\rho}|}\,.$ (IV.68)

For directions of propagation $\theta_{\mathbf{q}}$ between $0$ and ${\pi\over
2}$, this condition is equivalent to

$\theta_{\mathbf{q}}>\theta_{c}\equiv{\pi\over
2}-{2\sqrt{|B|}\over|v_{v}-v_{\rho}|}\,.$ (IV.69)

This result holds only for very small $q$. As noted earlier, even at very
small $q$, it only gives the region in which the instability growth rate is
linear in $q$; the possibility of an $O(q^{2})$ instability outside this small
angular window still exists. Furthermore, at larger $\mathbf{q}$, higher order
terms in the equation of motion become important. As a result, even
$\mathbf{q}$’s whose directions lie in the narrow angular wedge (IV.68) do not
remain unstable out to arbitrarily large $|\mathbf{q}|$.

To investigate both of these possibilities, we’ll now go beyond leading order
in $q$ for systems close to the instability threshold. Assuming that we
continue to have $q_{x}\sim\epsilon|\mathbf{q}_{\perp}|$ with $\epsilon\ll 1$,
we can greatly simplify the eigenvalue conditions (IV.38) and
(LABEL:vperpEq2). We’ll also assume, and verify a posteriori, that the typical
unstable wavevector satisfies $|\mathbf{q}_{\perp}|\propto{\epsilon}$. Then
the various terms in (IV.59) scale as follows:

$\displaystyle i\omega\delta v_{x}^{(1)}\sim\epsilon^{2}\delta
v_{x}^{(1)}\,\,\,\,,\,\,\,\,D_{\rho\perp}q_{\perp}^{2}\delta
v_{x}^{(1)}\sim\epsilon^{2}\delta v_{x}^{(1)}\,,$ $\displaystyle D_{\rho
x}q_{x}^{2}\delta v_{x}^{(1)}\sim\epsilon^{4}\delta
v_{x}^{(1)}\,\,\,\,,\,\,\,\,iv_{\rho}q_{x}\delta
v_{x}^{(1)}\sim\epsilon^{2}\delta v_{x}^{(1)}\,,$ $\displaystyle
i\rho_{0}q_{\perp}v_{L}\sim\epsilon v_{L}\,\,\,\,,\,\,\,\,D_{\rho
v}q_{x}q_{\perp}v_{L}\sim\epsilon^{3}v_{L}\,,$ $\displaystyle\phi\omega
q_{x}\delta v_{x}^{(1)}\sim\epsilon^{4}\delta v_{x}^{(1)}\,.$

Thus, as we approach the transition, so that $\epsilon\to 0$, we can neglect
the terms $D_{\rho x}q_{x}^{2}\delta\rho$ and $\phi\omega q_{x}\delta\rho$
relative to $\omega\rho$, and also drop the $D_{\rho v}q_{x}q_{\perp}v_{L}$
term relative to the $v_{\rho}q_{x}\delta\rho$ term. Doing so reduces our
$\delta\rho$ equation of motion (IV.59) to the far more manageable form:

$\displaystyle[-i(\omega-
v_{\rho}q_{x})+D_{\rho\perp}q_{\perp}^{2}]\delta\rho+i\rho_{0}q_{\perp}v_{L}$
$\displaystyle=0\,.$ (IV.70)

Performing this sort of power counting on the $v_{L}$ equation of motion
(IV.60) gives:

$\displaystyle i\omega v_{L}$
$\displaystyle\sim\epsilon^{2}v_{L}\,\,\,\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,v_{v}q_{x}v_{L}\sim\epsilon^{2}v_{L}\,,$
$\displaystyle Bq_{\perp}\delta v_{x}^{(1)}$
$\displaystyle\sim\epsilon^{3}\delta
v_{x}^{(1)}\,\,\,\,,\,\,\,\,D_{{}_{L\perp}}q_{\perp}^{2}v_{L}\sim\epsilon^{2}v_{L}\,,$
$\displaystyle D_{x}q_{x}^{2}v_{L}$
$\displaystyle\sim\epsilon^{4}v_{L}\,\,\,\,,\,\,\,\,\nu_{x}q_{x}q_{\perp}\delta
v_{x}^{(1)}\sim\epsilon^{3}\delta v_{x}^{(1)}\,,$ $\displaystyle\nu_{t}\omega
q_{\perp}\delta v_{x}^{(1)}$ $\displaystyle\sim\epsilon^{3}\delta
v_{x}^{(1)}\,,$

which shows that we can drop the $D_{x}q_{x}^{2}$ term in the $v_{L}$ equation
of motion (IV.60) relative to the $D_{{}_{L\perp}}q_{\perp}^{2}$ term. Doing
so gives:

$\displaystyle[-i(\omega-
v_{v}q_{x})+D_{{}_{L\perp}}q_{\perp}^{2}]v_{L}+\bigg{[}-i\frac{|B|}{\rho_{0}}q_{\perp}+\nu_{x}q_{x}q_{\perp}-\nu_{t}\omega
q_{\perp}\bigg{]}\delta\rho$ $\displaystyle=0\,.$ (IV.71)

These simplified equations of motion (IV.70) and (IV.71) can be combined into
a single vector equation:

$\bm{M}(\omega)\left[\begin{array}[]{cc}\delta\rho\\\\[10.0pt]
v_{L}\end{array}\right]=\left[\begin{array}[]{cc}0\\\\[10.0pt]
0\end{array}\right]\,,$ (IV.72)

where the matrix $\bm{M}$ is given by:

$\bm{M}=\left(\begin{array}[]{cc}[-i(\omega-
v_{\rho}q_{x})+D_{\rho\perp}q_{\perp}^{2}]&i\rho_{0}q_{\perp}\\\\[10.0pt]
\bigg{[}-i\frac{|B|}{\rho_{0}}q_{\perp}+\nu_{x}q_{x}q_{\perp}-\nu_{t}\omega
q_{\perp}\bigg{]}&[-i(\omega-
v_{v}q_{x})+D_{{}_{L\perp}}q_{\perp}^{2}]\end{array}\right)\;.$

The eigenvalue condition that determines the eigenfrequencies $\omega$ is
clearly simply

$\det\bm{M}(\omega)=0\,,$ (IV.73)

which leads to

Figure 8: Polar plot of the direction dependent sound speed $c_{\pm}(\theta)$
from equation (IV.67) ($B=-0.025,~{}v_{\rho}=1.50,~{}v_{v}=1$ (all units
arbitrary).). The horizontal axis, $x(\parallel)$, is the mean direction of
the flock’s motion. The solid blue line is the real part of the sound speed.
The dashed orange line is the imaginary part of the sound speed. Wavevectors
with non-zero imaginary speed values are unstable. This instability region is
highly anisotropic, with all of the unstable wavevectors very nearly
perpendicular to $y(\parallel)$. Note that the bifurcation of the real part
coincides with the opening angle of the imaginary part , $\theta_{c}$.

$\displaystyle-\omega^{2}+\omega[(v_{v}+v_{\rho})q_{x}$
$\displaystyle-i(D_{L\perp}+D_{\rho\perp}-\rho_{0}\nu_{t})q_{\perp}^{2}]$
$\displaystyle-|B|q_{\perp}^{2}+i(v_{\rho}D_{L\perp}+v_{v}D_{\rho\perp}-\rho_{0}\nu_{x})q_{x}q_{\perp}^{2}+D_{L\perp}D_{\rho\perp}q_{\perp}^{4}-v_{\rho}v_{v}q_{x}^{2}=0\,.$
(IV.74)

Figure 9: A plot of the highly anisotropic instability region in $q$-space.
There are unstable modes at all $\mathbf{q}$’s inside the two loops of this
lemniscate. The largest unstable wavevector $q_{c}$ is marked with a black
dash at the highest point on the instability boundary. The wavevectors $q_{m}$
with the fastest growing modes are identified with red crosses.

We’ll first determine the instability boundary. Instability occurs for those
$\mathbf{q}$’s for which the imaginary part ${\rm Im}\,\omega(\mathbf{q})$ of
$\omega(\mathbf{q})$ is positive ; modes with ${\rm Im}\,\omega(\mathbf{q})$
negative are stable. Therefore, the boundary between stable and unstable
regions of $\mathbf{q}$ must be the locus on which ${\rm
Im}\,\omega(\mathbf{q})=0$; i.e., the locus on which $\omega$ is purely real.
We can determine this locus by first taking the imaginary part of (IV.4) when
${\rm Im}\,\omega(\mathbf{q})=0$. This implies

$-(D_{{}_{L\perp}}+D_{\rho\perp}-\nu_{t}\rho_{0})q_{\perp}^{2}\omega+(v_{v}D_{\rho\perp}+v_{\rho}D_{{}_{L\perp}}-\rho_{0}\nu_{x})q_{\perp}^{2}q_{x}=0\,,$
(IV.75)

which can be trivially solved for the value $\omega_{inst}$ of $\omega$ at
which the instability occurs:

$\omega_{inst}=\left({v_{v}D_{\rho\perp}+v_{\rho}D_{{}_{L\perp}}-\rho_{0}\nu_{x}\over
D_{{}_{L\perp}}+D_{\rho\perp}-\rho_{0}\nu_{t}}\right)q_{x}\,.$ (IV.76)

Plugging this value of $\omega$ into the real part of (IV.4) then leads to a
linear equation for $q_{x}^{2}$, whose solution is

$\displaystyle
q_{x}^{2}=\frac{|B|q_{\perp}^{2}-D_{\rho\perp}D_{{}_{L\perp}}q_{\perp}^{4}}{V^{2}}\,,$
(IV.77)

where we’ve defined

$\displaystyle D_{3}$ $\displaystyle\equiv$ $\displaystyle
D_{{}_{L\perp}}+D_{\rho\perp}-\rho_{0}\nu_{t}\,,$ (IV.78) $\displaystyle
D_{6}$ $\displaystyle\equiv$ $\displaystyle{4\over
v_{4}}(v_{\rho}D_{{}_{L\perp}}+v_{v}D_{\rho\perp}-\nu_{x}\rho_{0})-2D_{3}\,,$
(IV.79) $\displaystyle V^{2}$ $\displaystyle\equiv$
$\displaystyle\left[{(v_{v}-v_{\rho})^{2}\over
4}-{D_{6}^{2}(v_{v}+v_{\rho})^{2}\over 16D_{3}^{2}}\right]\,.$ (IV.80)

The alert reader will recognize (IV.77) as the shape of the stability boundary
in wavevector space, as depicted in figure 1 , which we present again here,
with a new caption, as figure 9. She will also recognize this as equation (4)
of the associated short paperMiller and Toner .

Figure 10: Plot of the growth rate (IV.89). All wave vectors beyond $q_{c}$
decay. The fastest growth occurs at wavevector $q_{\perp}=q_{m}$, the position
of the maximum of this curve.

Note that the maximum value $q_{c}$ of $|\mathbf{q}_{\perp}|$ in the unstable
region is the value

$q_{c}=\sqrt{|B|\over D_{{}_{L\perp}}D_{\rho\perp}}$ (IV.81)

of $|\mathbf{q}_{\perp}|$ at which the right hand side of (IV.77) vanishes.
This is the justification for the preceding scaling argument, which assumed
that $|q_{\perp}|\lesssim\epsilon$.

Having found the instability boundary, we’d now like to determine the
wavevector of the fastest growing mode. As a first step, we’ll solve (IV.4)
for $\omega$:

$\displaystyle\omega$
$\displaystyle=\frac{1}{2}\bigg{[}v_{4}q_{x}-iD_{3}q_{\perp}^{2}\pm\sqrt{S}\bigg{]}$
(IV.82) $\displaystyle S$
$\displaystyle=(v_{\rho}-v_{v})^{2}q_{x}^{2}-D_{5}^{2}q_{\perp}^{4}-4|B|q_{\perp}^{2}+iv_{4}D_{6}q_{x}q_{\perp}^{2}$
(IV.83)

where we’ve defined

$\displaystyle v_{4}$ $\displaystyle\equiv v_{\rho}+v_{v},$ (IV.84)
$\displaystyle D_{5}^{2}$ $\displaystyle\equiv
D_{3}^{2}-4D_{{}_{L\perp}}D_{\rho\perp}.$ (IV.85)

Let’s first maximize the imaginary part of this over $q_{x}$. Note that all of
the $q_{x}$ dependence of the imaginary part comes from the square root.
Therefore, we need only maximize the imaginary part of that square root. This
implies we want to find the value of $q_{x}$ at which ${\partial{\rm
Im}\sqrt{S}\over\partial q_{x}}=0$. This condition is equivalent to saying
that ${\partial\sqrt{S}\over\partial q_{x}}={1\over 2\sqrt{S}}\left(\partial
S\over\partial{q_{x}}\right)$ is real. Squaring both sides of this condition
implies that ${1\over S}\left(\partial S\over\partial{q_{x}}\right)^{2}$ is
also real. Hence, at the wavevector of the fastest growing mode,

$\displaystyle{\left(\partial a\over\partial_{q_{x}}\right)^{2}-\left(\partial
b\over\partial_{q_{x}}\right)^{2}+2i\left(\partial
a\over\partial_{q_{x}}\right)\left(\partial b\over\partial_{q_{x}}\right)\over
S}={\rm real}\,,$ (IV.86)

where
$a\equiv(v_{\rho}-v_{v})^{2}q_{x}^{2}-D_{5}^{2}q_{\perp}^{4}-4|B|q_{\perp}^{2}$
and $b\equiv D_{6}q_{x}q_{\perp}^{2}$ are the real and imaginary parts of $S$
respectively. Calculating the required partial derivatives and performing this
algebra, we can rewrite this condition as:

$4(v_{\rho}-v_{v})^{2}\left({(v_{\rho}-v_{v})^{2}q_{x}^{2}-\left({D_{6}^{2}\over
4(v_{\rho}-v_{v})^{2}}\right)q_{\perp}^{4}+iD_{6}q_{x}q_{\perp}^{2}\over(v_{\rho}-v_{v})^{2}q_{x}^{2}-D_{5}^{2}q_{\perp}^{4}-4|B|q_{\perp}^{2}+iD_{6}q_{x}q_{\perp}^{2}}\right)={\rm
real}$ (IV.87)

Note that the imaginary parts of the numerator and denominator of the
parenthetical term in (IV.87) are equal. Therefore the condition can only be
satisfied in two ways: either

1) the real parts of the numerator and denominator are also equal, or

2) the imaginary parts of the numerator and denominator both vanish.

Note that condition 2) implies $q_{x}=0$. (We can’t have $q_{\perp}=0$ and
$q_{x}\neq 0$, since $q_{x}\ll q_{\perp}$ in the unstable region.) This proves
to be the condition at the maximum.

We will prove that condition 2) must hold by showing that condition 1) is
impossible. Assuming condition 1) is satisfied leads to the following
constraint on $q_{\perp}$:

$\displaystyle\bigg{(}-D_{5}^{2}+\frac{D_{6}^{2}}{4(v_{v}-v_{\rho})^{2}}\bigg{)}q_{\perp}^{2}$
$\displaystyle=4|B|.$ (IV.88)

This condition can never be satisfied in the unstable region, at least for
small $\nu_{x,t}$. To see this, note that when $\nu_{x}$ and $\nu_{t}$ are
both zero, the coefficient of $q_{\perp}^{2}$ on the left hand side of (IV.88)
vanishes, as can be seen from the expressions (IV.85) and (IV.80) for
$D_{5,6}$ evaluated for $\nu_{x,t}=0$. This means that for small $\nu_{x}$ and
$\nu_{t}$, the the coefficient of $q_{\perp}^{2}$ on the left hand side of
(IV.88) is guaranteed to be small compared to $D_{{}_{L\perp}}D_{\rho\perp}$.
That in turn implies that $|\mathbf{q}_{\perp}|>>\sqrt{|B|\over
D_{{}_{L\perp}}D_{\rho\perp}}=q_{c}$. However, such wavevectors lie outside
the unstable region, in which $|\mathbf{q}_{\perp}|<q_{c}$. Hence, this
solution is not acceptable, at least for small $\nu_{x}$ and $\nu_{t}$. This
is the limit on which we will henceforth focus. In this limit, we must
therefore have $q_{x}=0$ at the wavevector of fastest growth.

Setting $q_{x}=0$ in our expression (IV.82) for $\omega$ greatly simplifies
that expression. Taking ${\rm Im}~{}\omega$ gives the growth rate

$\displaystyle{{\rm
Im}~{}\omega}=-D_{3}q_{\perp}^{2}+\sqrt{4|B|q_{\perp}^{2}+D_{5}^{2}q_{\perp}^{4}}.$
(IV.89)

Maximizing this over $q_{\perp}$ gives, after a little algebra, a quadratic
equation for the square of the value $q_{m}$ of $q_{\perp}$ at which the
growth rate is maximized:

$\left({D_{5}^{2}\over
D_{{}_{L\perp}}D_{\rho\perp}}\right)q_{m}^{4}+4q_{c}^{2}q_{m}^{2}-q_{c}^{4}=0\,,$
(IV.90)

whose only real solutions are

$q_{m}=\pm
q_{c}\Bigg{(}\frac{1}{D_{5}}\sqrt{D_{{}_{L\perp}}D_{\rho\perp}\bigg{[}\sqrt{4+\left({D_{5}^{2}\over
D_{{}_{L\perp}}D_{\rho\perp}}\right)}-2\bigg{]}}\Bigg{)}\,,$ (IV.91)

whose magnitude the reader can verify for herself is always less than $q_{c}$,
as it must be.

###  IV.5 Intuitive explanation for the anisotropy of the instability

Why is the instability we just found restricted to modes with their
wavevectors $\mathbf{q}$ almost perpendicular to the mean velocity
$\langle\mathbf{v}\rangle$? The explanation is the competition between
convection and pressure forces.

By convection, we mean the $v_{\rho}$ and $v_{v}$ terms in equations (IV.38)
and (LABEL:vperpEq2) respectively. These would, in the absence of other
effects, cause any instantaneous local configuration of $\rho(\mathbf{r},t)$
or $\mathbf{v}(\mathbf{r},t)$ to simply be transported along the direction of
flock motion at speeds $v_{\rho}$ and $v_{v}$ respectively. Naively, one might
have expected $v_{\rho}$ and $v_{v}$ to both simply equal the speed of the
flock $|\langle\mathbf{v}\rangle|$ itself, since this would simply correspond
to the fluctuations moving along with the flock as it moves. However, as
explained in reference Toner and Tu (1995), due to the lack of Galilean
invariance in our system, there is no reason that either of these speeds
should be equal to $|\langle\mathbf{v}\rangle|$, since the flock itself does
not provide the only observable Galilean reference frame in the problem; there
is also the frictional substrate. Nor is there any reason that $v_{\rho}$ and
$v_{v}$ should be equal to each other. Indeed, their inequality plays a
crucial role in limiting the instability, as can be seen from, e.g., the
expression (IV.65) for $\epsilon$, which shows that $\epsilon\to\infty$ if
$v_{\rho}=v_{v}$. That is, the window of instability would have an opening
angle that covers all space, making all wavevectors unstable.

To understand why the inequality of $v_{\rho}$ and $v_{v}$ is so crucial,
consider the evolution of a small fluctuation in which both the density and
velocity depart from their mean values in a band perpendicular to the
direction of mean flock motion, as illustrated in figure 11. This corresponds
to a fluctuation with wavevectors making an angle $\theta_{\mathbf{q}}=0$ with
the mean velocity $\langle\mathbf{v}\rangle$. The fluctuation in the density
creates a pressure gradient, which, because the inverse compressibility $B$ is
negative, causes the velocity field to increase in the direction that tends to
increase the fluctuation in the density. That is, if the density is higher in
our band, the pressure forces tend to increase the velocity into the band,
making the density grow even more. The opposite will happen if the density is
lower than average in our band: pressure forces will tend to make the velocity
flow out, carrying boids out of the band and reducing the density even
further.

It is this positive feedback mechanism that is responsible for the familiar
instability of negative compressibility equilibrium systems. In our system,
however, this mechanism must compete against the aforementioned convection.
For the band direction just discussed, convection stabilizes the system by
separating the band of perturbed velocity from the band of perturbed density,
because they propagate at different speeds. This is illustrated in figure 11a.

Figure 11:  Three density and velocity fluctuations are depicted at different
angles $\theta$ relative to the mean motion of the flock. The density
fluctuation is depicted as a short rectangular region with cross hatching. The
density of the hash marks indicates the density in that region. The velocity
fluctuation is depicted as a long pink rectangular region with dashed edges.
Both regions propagate to the right, which is the in mean direction of motion
of the flock in all three panels, but at different speeds $v_{\rho}$ and
$v_{v}$ respectively, and extend infinitely in the directions parallel to
their depicted boundary. a) $\theta=0$, $c_{+}=v_{v}>v_{\rho}=c_{-}$. The
velocity fluctuation separates from the density fluctuation. The density
fluctuation begins to decay by diffusion, and the system is stable against
such fluctuations. b) $0<\theta<{\pi\over 2}$,the velocity and density bands
still separate over time, but by an a distance that’s reduced by a factor of
$\cos(\theta)$ when projected along the normal to the bands. As a result, if
$\theta>\theta_{c}$, the bands do not fully separate before the instability
sets in, and those fluctuations therefore grow, destablizing the system. c)
$\theta=\pi/2$. The density and velocity fluctuations both propagate parallel
to the band itself, so they continue to overlap for all time. This is
therefore the most unstable direction.

Once the velocity fluctuation is separated from the density fluctuation, both
simply decay away by diffusion. Therefore, this orientation of the bands is
stable.

Now consider the other extreme case, in which the band runs parallel to the
direction of mean flock motion, as illustrated in figure 11c. Now the relative
displacement of the velocity and density fluctuations makes no difference: the
two bands of each still overlap, since the bands are infinitely extended in
the direction of that relative motion. Hence, the velocity fluctuation
continues to sit right on top of the density fluctuation, allowing the
pressure induced instability to occur.

We can even derive the critical instability angle $\theta_{c}$ by this
argument. To see this, consider a fluctuation in which the initial band now
runs at an angle $\theta$ to the vertical, (which corresponds to the component
wavevectors $\mathbf{q}$ making an angle $\theta$ with the horizontal (i.e.,
with the $x$-axis)), as illustrated in 11b. Now the bands will again separate
by a distance $|v_{v}-v_{\rho}|t$ from each other in a time $t$, but that
separation is in the $x$-direction. The projection of that displacement
perpendicular to the bands (which is what determines when the bands actually
separate, of course) is therefore only $|v_{v}-v_{\rho}|t\cos\theta|$.

How long do these bands have to separate before the instability sets in? Well,
roughly speaking, the instability growth rate is $\sqrt{|B|}q$, where the
wavenumber $q\sim 1/W$, with $W$ the width of the band. Hence, the
characteristic time for the instability to grow is the inverse of this rate;
i.e., $t\sim{W\over\sqrt{|B|}}$. In this time, the bands will separate
perpendicular to the bands by a distance

$\displaystyle|v_{v}-v_{\rho}|t\cos\theta\sim{|v_{v}-v_{\rho}|\cos\theta
W\over\sqrt{|B|}}\,.$ (IV.92)

Requiring that this distance be greater than the width $W$ itself leads to the
condition:

${|v_{v}-v_{\rho}|W|\cos\theta|\over\sqrt{|B|}}>W$ (IV.93)

or, cancelling off the common factor of the band width $W$,

${|v_{v}-v_{\rho}||\cos(\theta)|\over\sqrt{|B|}}>1\,.$ (IV.94)

This is the condition for stability, because only if this condition is
satisfied will the bands separate before the fluctuation grows. This condition
can be rewritten

$|\cos(\theta)|>{\sqrt{|B|}\over|v_{v}-v_{\rho}||}\,.$ (IV.95)

For small $|B|$ (i.e., near the instability), the right hand side of this
inequality becomes very small, so that we can expand the cosine on the left
hand side for $\theta$ near $\pm{\pi\over 2}$. This leads to the stability
condition

$|\theta\pm\pi/2|>{\sqrt{|B|}\over|v_{v}-v_{\rho}|}\,,$ (IV.96)

which the alert reader will recognize as equivalent to the instability
condition (IV.68), up to a factor of $2$, which one wouldn’t expect to get
correctly from an argument as crude as this one.

###  IV.6 Summary of the stability analysis

In summary, we have shown in this section that autochemotaxis can make the
inverse compressibility negative. Next, we demonstrated that any negative
inverse compressibility, no matter how small, always destabilizes the uniform
homogeneous flock. However, in contrast to disordered systems, ordered flocks
near this instability are only unstable in a narrow region of region of
$\mathbf{q}$ space, very near the direction perpendicular to the direction of
flock motion, as illustrated in figure 1. Furthermore, the wavevector of the
most unstable mode is precisely perpendicular to the mean direction of flock
motion.

Finally, we note that the unstable wavevectors are all necessarily small near
the threshold of the instability (i.e., when the inverse compressibility $B$
is negative, but its magnitude $|B|$ is small).

We will use all of the above facts in the next section to help us determine
the final steady state the system near, but on the unstable side of, the
instability threshold.

##  V The Final state: Phase separation

###  V.1 Formulation

We’d now like to determine the final steady state of the system when the
instability found in the last section occurs. We will show in this section
that, as claimed in the introduction, this state is “phase separated”: the
system separates into two macroscopic bands running parallel to the mean
velocity, with essentially constant density within each. The interface between
these two bands has a width $\ell$ that is independent of the system width
$w$. Therefore, the ratio ${\ell\over w}$ of the interface width to the system
width vanishes as system size $w\to\infty$. The interface width $\ell$ itself
diverges as the system approaches a non-equilibrium “critical point” that we
will determine.

Since the instability is sharply focused in the directions perpendicular to
the bulk motion, we will seek solutions of our equations of motion that depend
only on one Cartesian component of position, by which we mean configurations
in which all of the fields depend on only one of the $d-1$ Cartesian
coordinates perpendicular to $\hat{x}$ ($d$ being the dimension of space).
We’ll call this direction $y$.

In addition, since the instability was in a mode in which only the
“longitudinal” component of the perpendicular velocity $\mathbf{v}_{\perp}$
(i.e., the component $v_{L}$ along the direction of the projection of the
wavevector of the instability perpendicular to the mean flock velocity vector
$\langle\mathbf{v}\rangle$) was non-zero, we will seek solutions in which only
that component of $\mathbf{v}_{\perp}$ \- which, in this geometry, means the
y-component- is non-zero.

Note that in addition there is a fluctuating and position dependent component
of the velocity in the $\hat{x}$ direction. This is due to the density
dependence of $U(|\mathbf{v}|,\rho)$, which, as the reader may recall,
maintains the non-zero speed of the flock.

Specifically, we’ll search for solutions of the form

$\mathbf{v}_{\perp}=v_{y}(y,t)\hat{y}\,\,\,\,\,,\,\,\,\,\,\rho=\rho_{ref}+\Delta\rho(y,t)\,,$
(V.1)

where $\rho_{ref}$ is some reference density. Because we wish to consider the
effect of changing the mean density $\rho_{0}$ of the system, we will chose
this reference density to in general differ from $\rho_{0}$.

The chemically sensitive reader will notice that we have not included the
chemo-attractant concentration $\eta$ in the above analysis. This is because,
as we did when we linearized our equations of motion in section IV, we can
eliminate this field. We can similarly eliminate the fluctuating component of
the velocity along $\langle\mathbf{v}\rangle$. See section IV for more
details.

This approach is justified whenever the solutions we are studying vary slowly
in space and time. Since we are looking for steady-state solutions, slow
variation in time is guaranteed (no variation is as slow as you can get!). We
will verify a posteriori that the solutions we find do indeed vary slowly in
space, for parameters near a non-equilibrium “critical point”. We will
identify this “critical point” as we proceed.

As discussed in the introduction, our system proves to be closely analogous to
equilibrium phase separation. In particular, it proves, as the analytic
solution we’ll present below shows, to have the analog of the liquid-vapor
critical pointChaikin and Lubensky (1995) of equilibrium phase separation. We
will choose our reference density $\rho_{ref}$ to be close to the critical
density of our non-equilibrium critical point. This justifies our expansion in
powers of $\Delta\rho$, since, as we’ll see, it implies that the density
$\rho(\mathbf{r})$ is never very far from $\rho_{ref}$ at any point in space
in the phase separated state, provided that the mean density $\rho_{0}$ is
close to $\rho_{ref}$.

We will not go through the process of eliminating the fast variables in detail
here. We will content ourselves with noting that the resultant equations of
motion must be of the form of the Toner-Tu equations (II.3) and (II.4). This
is an inevitable result of the fact that the symmetries and conservation laws
of our problem are the same as those considered in the formulation of the
Toner-Tu model.

Note that, in contrast to the stability analysis of the previous section, we
will not linearize our equations of motion. However, we will still truncate
them to quadratic order in the fluctuations $v_{y}(y)$, to cubic order in
$\Delta\rho(y)$ , and to bilinear order in $v_{y}\Delta\rho$. This will be
sufficient if $v_{y}(y)$ and $\Delta\rho(y)$ are small, as we will again
verify they are a posteriori.

We will further gradient expand our model. In particular, terms involving 2
derivatives of any type ($x$, $\perp$, or $t$) will be expanded only to linear
order in the fields. This is justified by our expectation that, near the
critical point, the density varies very slowly in space and time. We will
verify this assumption a posteriori.

The resultant equations readw2f :

$\displaystyle\partial_{t}\mathbf{v}_{\perp}+v_{v}\partial_{x}\mathbf{v}_{\perp}+\lambda\left(\mathbf{v}_{\perp}\cdot\nabla_{\perp}\right)\mathbf{v}_{\perp}$
$\displaystyle=-g_{1}\Delta\rho\partial_{x}\mathbf{v}_{\perp}-g_{2}\mathbf{v}_{\perp}\partial_{x}\Delta\rho-\nabla_{\perp}P_{1}+D_{B}\nabla_{\perp}(\nabla_{\perp}\cdot\mathbf{v}_{\perp})+D_{T}\nabla^{2}_{\perp}\mathbf{v}_{\perp}$
$\displaystyle~{}+D_{x}\partial_{x}^{2}\mathbf{v}_{\perp}+\nu_{x}\partial_{x}\nabla_{\perp}\Delta\rho+\nu_{t}\partial_{t}\nabla_{\perp}\Delta\rho\,,$
(V.2)
$\displaystyle\partial_{t}\Delta\rho+\rho_{ref}\nabla_{\perp}\cdot\mathbf{v}_{\perp}+\nabla_{\perp}\cdot(\mathbf{v}_{\perp}\Delta\rho)+v_{\rho}\partial_{x}\Delta\rho$
$\displaystyle=D_{\rho
x}\partial^{2}_{x}\Delta\rho+D_{\rho\perp}\nabla_{\perp}^{2}\Delta\rho+D_{\rho
v}\partial_{x}\left(\nabla_{\perp}\cdot\vec{v}_{\perp}\right)$
$\displaystyle~{}+\phi\partial_{t}\partial_{x}\Delta\rho+w_{2}\partial_{x}(\Delta\rho^{2})+{\rho_{ref}\over
2v_{0}}\partial_{x}(|\vec{v}_{\perp}|^{2})~{},$ (V.3)

where we’ve defined $\lambda\equiv\lambda_{1}(\rho_{ref},v_{0})$.

Inserting our one-dimensional ansätze (V.1) into the Toner-Tu equations (V.2)
and (V.3), and dropping terms of higher than quadratic order in $v_{y}(y)$ and
higher than cubic order in $\Delta\rho(y)$ leads to PDE’s for the one-
dimensional evolution:

$\displaystyle\partial_{t}v_{y}+\lambda
v_{y}\partial_{y}v_{y}=D_{{}_{L\perp}}\partial_{y}^{2}v_{y}-\partial_{y}P_{1}+\nu_{t}\partial_{t}\partial_{y}\Delta\rho,$
(V.4)
$\displaystyle\partial_{t}\Delta\rho=D_{{}_{C}}\partial_{y}^{2}\Delta\rho-\rho_{ref}\partial_{y}v_{y}-\partial_{y}(v_{y}\Delta\rho),$
(V.5)

where we’ve defined

$\displaystyle D_{L\perp}\equiv D_{B}+D_{T}\,\,\,\,\,,\,\,\,\,D_{{}_{C}}\equiv
D_{\rho\perp}\,.$ (V.6)

In the steady-state, these reduce to two coupled ODE’s:

$\displaystyle\lambda v_{y}{dv_{y}\over dy}=D_{{}_{L\perp}}{d^{2}v_{y}\over
dy^{2}}-{dP_{1}\over dy},$ (V.7) $\displaystyle
D_{{}_{C}}{d^{2}\Delta\rho\over dy^{2}}-\rho_{ref}{dv_{y}\over dy}-{d\over
dy}(v_{y}\Delta\rho)=0\,.$ (V.8)

We will begin by using equations (V.7) and (V.8) to analytically determine the
steady state solution. To validate our solution, we numerically solved the
dynamical equations (V.4) and (V.5), and allowed them to evolve for long
times.

In the previous sections we demonstrated the uniform state with
$\rho(\mathbf{r})=\rho_{0}$ and $v_{y}=0$ becomes unstable when the inverse
compressibility $B$ changes sign. In order to find a stable steady state, it
is therefore necessary to extend our expansion of the isotropic pressure
$P_{1}$ to third order in $\Delta\rho$:

$\displaystyle P_{1}=P_{0}+M\Delta\rho+w\Delta\rho^{2}+u\Delta\rho^{3}\,.$
(V.9)

Recall that, in the linear theory, we found that the inverse compressibility
$B$ appearing in the linearized equations of motion become negative for
sufficiently strong attractive autochemotaxis. It is therefore clear that the
expansion coefficient $M$ in the above expansion will likewise be driven from
positive to negative values by increasing autochemotaxis.

We will therefore treat $M$ as an experimentally controllable parameter, which
can be decreased by increasing the chemotaxis. The other experimentally
controllable parameter at our disposal is the mean density $\rho_{0}$.

The experimental situation is thus very similar to studying equilibrium phase
separation by tuning both the system’s mean density $\rho_{0}$ and the
temperature, with, in our problem, chemotaxis - or its proxy $M$\- playing the
role of temperature. We will therefore determine the “phase diagram” of our
system in the $\rho_{0}$-$M$ plane.

###  V.2 Stability boundary: the “spinodal” line

We’ll begin by determining the stability boundary in this $\rho_{0}$-$M$
plane. As shown in the last section, the instability occurs when the inverse
compressibility vanishes; that is, when

$B\equiv\left({dP_{1}\over d\rho}\right)\bigg{|}_{\rho_{0}}=0\,.$ (V.10)

We will define the value of $M$ that satisfies this condition to be
$M_{spinodal}$, since it is the precise analog of the spinodal line in
equilibrium phase separation: that is, the line on which the uniform density
state becomes unstable.

Inserting our expansion (V.9) into the condition (V.10) gives

$M_{\rm spinodal}+2w(\rho_{0}-\rho_{ref})+3u(\rho_{0}-\rho_{ref})^{2}=0\,,$
(V.11)

which can easily be solved for the spinodal line:

$M_{\rm spinodal}=M_{c}-3u(\rho_{0}-\rho_{c})^{2}\,,$ (V.12)

where we’ve defined the “critical” value $M_{c}$ of $M$ as

$M_{c}\equiv{w^{2}\over 3u}\,,$ (V.13)

and the “critical density”

$\rho_{c}\equiv\rho_{ref}-{w\over 3u}\,.$ (V.14)

The critical values $M_{c}$ and $\rho_{c}$ of $M$ and $\rho$ prove to be the
analogs of the critical temperature and critical density of an equilibrium
phase separating system.

Using (V.14), we can rewrite (V.7) and (V.8) in terms of the critical density
$\rho_{c}$ and a new density variable $\rho^{\prime}$ defined as the
difference between the local density and the critical density:

$\rho^{\prime}=\rho-\rho_{c}\,.$ (V.15)

Those equations read:

$\displaystyle\lambda v_{y}{dv_{y}\over dy}=D_{{}_{L\perp}}{d^{2}v_{y}\over
dy^{2}}-{dP_{1}\over dy},$ (V.16) $\displaystyle
D_{{}_{C}}{d^{2}{\rho^{\prime}}\over dy^{2}}-\rho_{c}{dv_{y}\over dy}-{d\over
dy}(v_{y}{\rho^{\prime}})=0\,.$ (V.17)

###  V.3 Two-phase coexistence boundary: the “binodal” line and the common
tangent construction

Now we turn to the question of the phase separated state of the system. In
discussing this, it is useful to define the analog $F(\rho)$ of an equilibrium
free energy for our system. We can do so by defining

$P_{1}={dF(\rho)\over d\rho}\,,$ (V.18)

which is readily seen to imply

$F(\rho)=F_{0}+P_{0}\Delta\rho+{M\over 2}\Delta\rho^{2}+{w\over
3}\Delta\rho^{3}+{u\over 4}\Delta\rho^{4}\,.$ (V.19)

We will show below that the two densities into which the system ultimately
phase separates in the steady state are determined, sufficiently close to the
critical point

$M=M_{c}\ \ \ ,\ \ \ \rho=\rho_{c}\,,$ (V.20)

by the familiar “common tangent construction” of equilibrium statistical
mechanics, applied to the “free energy” (V.19). This is illustrated in figure
12.

This analogy to equilibrium statistical mechanics only applies close to the
critical point (V.20). We will also show below that, as we move away from the
critical point, the densities are determined by an “uncommon tangent
construction”: the slope of the free energy (V.19) remains the same at both
coexisting densities, but the tangents to the free energy curve at those two
densities do not coincide, as illustrated in figure 12.

Figure 12: Illustration of the common tangent construction. The densities of
the high and low density regions are the points that share a common tangent
line. In equilibrium physics, each phase has the same thermodynamic pressure
and chemical potential.

As pointed out earlier, our calculation here is very similar to Landau theory
for equilibrium liquid-gas phase separationChaikin and Lubensky (1995). In
particular, as in equilibrium Landau theory, our truncation of the expansion
of the pressures in powers of $\Delta\rho$ at third order is only justified
very near the “critical point” (V.20) for our system, and only if the critical
density $\rho_{c}$ is close to the reference density $\rho_{ref}$ around which
we did our expansion. Of course, since we are free to choose the reference
density $\rho_{ref}$ to be anything we want, we can always choose it to be
close to the critical density $\rho_{c}$. Indeed, we could have chosen it to
be equal to $\rho_{c}$; with that choice, we would have had $w=0$. Hence, if
our reference density $\rho_{ref}$ is close to the critical density, $w$ will
be small. We will also restrict our discussion to the vicinity of the critical
point (V.20), which implies that $M$ is small as well. We will therefore make
liberal use of the assumption of small $M$ and $w$ in all of the analysis that
follows.

All of these assumptions are precisely the same as those usually made in the
Landau theory analysis of an equilibrium critical point for phase
separationChaikin and Lubensky (1995). We therefore expect that our analytic
results obtained below to be as reliable as those obtained from equilibrium
Landau theory: that is, valid (in the absence of fluctuation effects)
sufficiently close to the critical point.

For computational convenience, we will shift variable from $\Delta\rho$ to the
new density variable $\rho^{\prime}$ defined in (V.15) as the difference
between the local density and the critical density. Using our expression
(V.14) for the critical density, we see that

$\Delta\rho=\rho^{\prime}-{w\over 3u}\,.$ (V.21)

This shift eliminates the quadratic term in the pressure, or, equivalently,
the cubic term in the free energy, (V.19), enabling us to write

$\displaystyle P_{1}(\rho^{\prime})={dF(\rho^{\prime})\over d\rho^{\prime}}$
(V.22)

with

$F(\rho^{\prime})=P_{0}\rho^{\prime}-{1\over 2}m\rho^{\prime 2}+{1\over
4}u\rho^{\prime 4}+{\rm constant}\,,$ (V.23)

where we’ve defined

$m\equiv{w^{2}\over 3u}-M\,.$ (V.24)

We have not bothered evaluating the constant in (V.23), since it drops out of
the equation of motion (V.16).

Note also that since we are limiting our analysis to small $M$ and $w$, $m$ is
also small. We will make use of this fact in the perturbation theory we
develop for the non-equilibrium phase separation later.

With the expansion (V.23) in hand, we now further simplify (V.16) and (V.17)
by rescaling to dimensionless coordinates and variables, as follows:

$s\equiv{y\over\ell}\,\,\,\,\,\,,\,\,\,\,\,\tilde{v}\equiv{v_{y}\over
v_{s}}\,\,\,\,\,\,,\,\,\,\,\,\tilde{\rho}\equiv{\rho^{\prime}\over\rho_{s}}\,.$
(V.25)

with the characteristic length $\ell$, speed $v_{s}$, and density $\rho_{s}$
given by

$\displaystyle\ell$
$\displaystyle=\sqrt{\frac{D_{{}_{L\perp}}D_{{}_{C}}}{\rho_{c}m}}\,,$ (V.26)
$\displaystyle v_{s}$
$\displaystyle=\frac{m\sqrt{D_{{}_{C}}}}{\sqrt{D_{{}_{L\perp}}\rho_{c}u}}\,,$
(V.27) $\displaystyle\rho_{s}$ $\displaystyle=\sqrt{\frac{m}{u}}\,.$ (V.28)

We’ll see shortly that the length $\ell$ is roughly the width of the interface
between phase separated regions, while the difference between the densities of
the two phase separated regions is roughly $\rho_{s}$.

In terms of these new coordinates and variables, our dimensionless steady-
state equations of motion (V.16) and (V.17) read:

$\displaystyle{d^{2}\tilde{v}\over ds^{2}}$ $\displaystyle={d\over
ds}\left({d\tilde{F}(\tilde{\rho})\over
d\tilde{\rho}}\right)+\Lambda\tilde{v}{d\tilde{v}\over ds}\,,$ (V.29)
$\displaystyle{d^{2}\tilde{\rho}\over ds^{2}}$ $\displaystyle={d\tilde{v}\over
ds}+\Gamma{d\over ds}(\tilde{\rho}\tilde{v})\,,$ (V.30)

where we’ve defined the dimensionless parameters:

$\displaystyle\Lambda$ $\displaystyle=\frac{\lambda
D_{c}}{D_{{}_{L\perp}}\rho_{c}}\sqrt{\frac{m}{u}}=\frac{\lambda
D_{c}\rho_{s}}{D_{{}_{L\perp}}\rho_{c}}\,,$ (V.31) $\displaystyle\Gamma$
$\displaystyle=\frac{1}{\rho_{c}}\sqrt{\frac{m}{u}}={\rho_{s}\over\rho_{c}}\,,$
(V.32)

and the dimensionless free energy

$\tilde{F}(\tilde{\rho})\equiv{\ell\over\rho_{s}v_{s}D_{{}_{L\perp}}}F(\rho^{\prime}=\rho_{s}\tilde{\rho})=\tilde{P}_{0}\tilde{\rho}-{1\over
2}\tilde{\rho}^{2}+{1\over 4}\tilde{\rho}^{4}\,.$ (V.33)

In (V.33), we’ve defined the dimensionless pressure

$\tilde{P}_{0}\equiv{P_{0}\ell\over D_{{}_{L\perp}}v_{s}}\,.$ (V.34)

In doing this rescaling, we have implicitly assumed that $m$ is positive
($m>0$), as it obviously will be for $M<M_{c}$ (see equation (V.24).

Note that, because $m$ is small near the critical point, both of our
dimensionless parameters $\Gamma$ and $\Lambda$ are also small near the
critical point. Therefore, to leading order near the critical point, we can
set them to zero. Doing so gives, for the steady state,

$\displaystyle{d^{2}\tilde{v}\over ds^{2}}$ $\displaystyle={d\over
ds}\left({d\tilde{F}(\tilde{\rho})\over d\tilde{\rho}}\right)\,,$ (V.35)
$\displaystyle{d^{2}\tilde{\rho}\over ds^{2}}$ $\displaystyle={d\tilde{v}\over
ds}\,.$ (V.36)

Differentiating (V.36) with respect to $s$ and using equation (V.35) in the
result gives

$\displaystyle{d^{3}\tilde{\rho}\over ds^{3}}$ $\displaystyle={d\over
ds}\left({d\tilde{F}(\tilde{\rho})\over d\tilde{\rho}}\right).$ (V.37)

Integrating this once with respect to $s$ then gives

${d^{2}\tilde{\rho}\over ds^{2}}=\left({d\tilde{F}(\tilde{\rho})\over
d\tilde{\rho}}\right)+\mu\,,$ (V.38)

where $\mu$ is a constant of integration. Multiplying both sides of (V.38) by
$\frac{d\tilde{\rho}}{ds}$ and integrating with respect to $s$ again reveals
that this equation has a conserved first integral of motion. That is, the
quantity

$E\equiv-{1\over 2}\left({d\tilde{\rho}\over
ds}\right)^{2}+\tilde{F}(\tilde{\rho})+\mu\tilde{\rho}$ (V.39)

is a constant (independent of position $s$).

Assuming that the density profile in the steady state consists of a set of
“plateaus” - i.e., regions of essentially constant density - we can calculate
the possible plateau densities, which we’ll call $\tilde{\rho}_{\pm}$ (with
$\tilde{\rho}_{-}$ being the smaller of the two), from equations (V.38) and
(V.39). The result, as we’ll now show, is identical to the familiar “common
tangent construction” of equilibrium phase separation.

On a plateau, all spatial derivatives must vanish. Hence, (V.38) implies

${d\tilde{F}\over d\tilde{\rho}}\bigg{|}_{\tilde{\rho}_{\pm}}=-\mu\,.$ (V.40)

That is, at the plateau densities, the derivatives of the effective free
energy $\tilde{F}$ must be equal. Since $\tilde{F}$ and $\tilde{\rho}$ are
linearly related to to our dimensionful free energy $F(\rho)$ and $\rho$, this
condition clearly also applies to the derivatives of $F$ with respect to
$\rho$.

Now setting derivatives equal to zero in (V.39), and using the fact that the
“energy” $E$ must be the same on both plateaus, we obtain a relation between
the values of $\tilde{F}$ itself (as opposed to its derivatives) on the
plateaus:

$\tilde{F}(\tilde{\rho}_{+})+\mu\tilde{\rho}_{+}=\tilde{F}(\tilde{\rho}_{-})+\mu\tilde{\rho}_{-}\,.$
(V.41)

This expression can be reorganized to read

$\tilde{F}(\tilde{\rho}_{+})=\tilde{F}(\tilde{\rho}_{-})-\mu(\tilde{\rho}_{+}-\tilde{\rho}_{-})\,.$
(V.42)

Using the fact that (V.40) implies that ${d\tilde{F}\over
d\tilde{\rho}}|_{{}_{\tilde{\rho}_{\pm}}}=-\mu$, we see immediately that
(V.40) and (V.42) taken together imply the “common tangent construction”: a
tangent to a plot of $\tilde{F}(\tilde{\rho})$ at one of the plateau densities
$\tilde{\rho}_{\pm}$ must also be tangent to the plot at the other plateau
density. As for the common slope relation, here too, this relation must apply
to the original, unrescaled density $\rho$ and free energy $F$ as well.

This common tangent construction is illustrated in figure 12.

To solve for these two steady state densities, it is most convenient to define
a new “pseudo-Gibbs free energy” $G$ via

$G(\tilde{\rho})\equiv\tilde{F}+\mu\tilde{\rho}=(\tilde{P}_{0}+\mu)\tilde{\rho}-{1\over
2}\tilde{\rho}^{2}+{1\over 4}\tilde{\rho}^{4}\,.$ (V.43)

Keep in mind that $\mu$ is an as yet undetermined constant of integration,
which we must adjust in order to obtain plateaus in the density.

With the definition (V.43), we see that the condition (V.40) on the slope of
the effective free energy $\tilde{F}$ at the plateau densities becomes the
simpler condition

${dG\over d\tilde{\rho}}\bigg{|}_{\tilde{\rho}_{\pm}}=0\,,$ (V.44)

i.e., $G$ must be minimized (or maximized) at $\tilde{\rho}_{\pm}$.

Furthermore, the condition (V.41) implies

$G(\tilde{\rho}_{+})=G(\tilde{\rho}_{-})\,.$ (V.45)

Taking the conditions (V.44) and (V.45) together implies that the arbitrary
constant of integration $\mu$ must be chosen so that $G$ has two degenerate
minima. Inspection of (V.43) makes it fairly obvious what choice of $\mu$ will
accomplish this:

$\mu=-\tilde{P}_{0}\,.$ (V.46)

With this choice, we have

$G(\tilde{\rho})=-{1\over 2}\tilde{\rho}^{2}+{1\over 4}\tilde{\rho}^{4}\,,$
(V.47)

which clearly has two degenerate minima at

$\tilde{\rho}_{\pm}=\pm 1\,.$ (V.48)

Undoing the rescaling (V.25) and the shift (V.15), we see that this implies
the two plateau densities $\rho_{\pm}$ are given by

$\rho_{\pm}=\rho_{c}\pm\rho_{s}=\rho_{c}\pm\sqrt{\frac{m}{u}}\,.$ (V.49)

Note that, for small $m$ (i.e., close to the critical point), these densities
are close to the critical density $\rho_{c}$, which justifies our expansion in
$\Delta\rho$.

However, the two phase separated states only become accessible if the mean
density $\rho_{0}$ lies between these two phase separated densities:

$\rho_{+}>\rho_{0}>\rho_{-}\,\,\,\,\,.$ (V.50)

Otherwise, conservation of particle number forbids going to a phase separated
state (since both plateaus would have either have a lower density than the
mean density $\rho_{0}$, or a higher density). In either case, there is no way
that one can achieve the mean density $\rho_{0}$ by putting part of the system
in a plateau at $\rho_{+}$, and the remainder in a plateau of density
$\rho_{-}$. It is only in the regime (V.50) that such a separation can lead to
the mean density $\rho_{0}$.

Using our result (V.49) for $\rho_{\pm}$, we see that the condition (V.50) is
equivalent to

$\sqrt{\frac{m}{u}}>|\rho-\rho_{c}|\,.$ (V.51)

Using our expression (V.24) relating $m$ and the parameter $M$ in our
expansion (V.9) for the isotropic pressure $P_{1}$, we see that this condition
can be written as

$M<{w^{2}\over 3u}-u(\rho-\rho_{c})^{2}\equiv M_{\rm binodal}\,.$ (V.52)

Note that $M_{\rm binodal}>M_{\rm spinodal}$ for all $\rho$. In fact, the
ratio

${(M_{\rm spinodal}-M_{c})\over(M_{\rm binodal}-M_{c})}=3$ (V.53)

is universal for all $\rho$ sufficiently close to the critical point.

The result (V.52), along with our earlier result (V.12) for the “spinodal”
value $M_{\rm spinodal}$ of $M$ at which the uniform state becomes unstable,
can be summarized by the “phase diagram” in the $\rho_{0}$-$M$ plane
illustrated in figure 13. For all $M$ above the “binodal” parabola, which is
the upper parabola in 13, the system can only be in the uniform state, which
is stable. For values of $M$ between the two parabolae (the blue and orange
curves in figure 13) - that is, for $M_{\rm binodal}>M>M_{\rm spinodal}$, with
$M_{\rm spinodal}$ and $M_{\rm binodal}$ given by (V.12) and (V.52), both the
two phase state, and the homogeneous, one phase state, are stable. Finally,
for $M<M_{\rm spinodal}$, only the phase separated state is stable.

The strong similarity between our results and equilibrium liquid-vapor phase
separation is apparent from this phase diagram, which is identical to that for
an equilibrium liquid-vapor system, with $M$ playing the role of temperature.
Recall that $M$ can be increased (decreased) experimentally by decreasing
(increasing) the strength of the chemotaxis, or, more generally, by tuning the
strength of whatever attractive interactions in the flock reduce the inverse
compressiblity.

Note there is an important difference between our results and equilibrium
liquid-vapor phase separation. In equilibrium systems, a homogeneous density
state is only meta-stable in the region between the binodal and spinodal
curves. The homogeneous state exists at a local minimum in the free energy;
the global minimum is the phase separated state. In the presence of
sufficiently large fluctuations, the homogeneous phase will overcome this
energy barrier through droplet nucleationChaikin and Lubensky (1995), which
leads to phase separation. Since flocks are non-equilibrium systems, we do not
have any reason to believe that the global minimum of our pseudo-energy ($G$)
corresponds to a “preferred” stable state. All we can determine is the local
stability of each phase between the binodal and spinodal curves. Between these
two curves, both states are locally stable. We therefore refer to this region
as the bistable region.

Figure 13: Phase diagram of flock phase separation. The solid blue and dashed
orange curves are the binodal and spinodal lines respectively , analytic
expressions for which are given by equations (V.52) and (V.12). Note that
those expressions are only valid close to the critical point. In the orange
filled region under the spinodal line, which is labeled “unstable”, only the
two-phase state is stable. In the blue region between the spinodal and binodal
lines, both the two phase state, and the homogeneous, one phase state, are
stable. The analysis we perform in this paper is done close to the critical
point, where our assumption that density variations are small is valid. Figure
14: A result from evolving equation (V.29) and (V.30) for a long time, in a
periodic space. The dashed blue curve is the initial condition, and the orange
curve is the final state. This solution was iterated for 3000 units of time
and $\Lambda,\Gamma=0$. The plateaus are very weakly attracted to one another.
The long time steady state solution in a periodic space would consist one a
single high density region and a single low density region.

We turn now from the plateau densities to the form of the interface between
different plateaus. We will determine this using equation (V.39). Using our
earlier result that $\mu=-\tilde{P}_{0}$, and evaluating the right hand side
of (V.39) on a plateau, where $\tilde{\rho}=\pm 1$ and ${d\tilde{\rho}\over
ds}=0$, we see that

$E=-{1\over 4}\,.$ (V.54)

Using this in (V.39), we can solve that equation for ${d\tilde{\rho}\over
ds}$:

${d\tilde{\rho}\over ds}=\pm\sqrt{{1\over
2}(\tilde{\rho}^{4}-2\tilde{\rho}^{2}+1)}=\pm{1\over\sqrt{2}}(1-\tilde{\rho}^{2})\,,$
(V.55)

where the plus sign corresponds to the transition region from the low density
plateau at $\tilde{\rho}=-1$ on the left to the high density plateau at
$\tilde{\rho}=1$ on the right, and the minus sign to the opposite case.
Solving this first order ODE by separation of variables, we find the interface
structure:

$\tilde{\rho}(s)=\pm\tanh\left({s-s_{0}\over\sqrt{2}}\right)\,,$ (V.56)

where the constant of integration $s_{0}$ simply gives the position of the
interface.

Undoing our various shifts (V.21) and rescalings (V.25), we can use this to
write the density profile of the interface in terms of the physical variable
$\rho$ and the physical coordinate $y$:

$\displaystyle\rho(y)$ $\displaystyle=$ $\displaystyle\rho_{ref}-{w\over
3u}\pm\sqrt{\frac{{w^{2}\over
3u}-B}{u}}\tanh\left({y-y_{0}\over\ell\sqrt{2}}\right)\,,$

where as before the plus sign describes the case in which the low density
plateau is on the left (i.e., as $y\to-\infty$) and the high density on the
right, while the minus sign describes the opposite case.

Note that, as claimed earlier, the characteristic width of these interfaces is
of order $\ell$, and hence diverges like $m^{-1/2}$ as we approach the
critical point $m=0$, as can be seen from equation (V.26).

Having found the density $\rho$, we can now also determine the velocity field
$\mathbf{v}$. The $y$-component follows from equation (V.36), which can be
integrated once to give

$\tilde{v}={d\tilde{\rho}\over ds}+{\rm
constant}=\pm{1\over\sqrt{2}\cosh^{2}\left({s-s_{0}\over\sqrt{2}}\right)}+{\rm
constant}\,.$ (V.58)

The constant of integration in this expression is readily seen to be zero,
since, as can be seen from our full time-dependent 1d equations of motion
(V.7) and (V.8) the volume integral of $\tilde{v}$ is conserved. This volume
integral was zero initially, because we started with a system moving uniformly
in the $x$-direction, so $v_{y}{(y,t=0)}=0$ for all $y$. Hence, this volume
integral must be zero in the final state, as it will be if we choose the
constant of integration in (V.58) to be zero (note that a system with periodic
boundary conditions must have equal numbers of “up” and “down” interfaces, so
the small contributions to $\int dyv_{y}(y)$ coming from these regions will
cancel).

Thus we have

$\tilde{v}=\pm{1\over\sqrt{2}\cosh^{2}\left({s-s_{0}\over\sqrt{2}}\right)}\,.$
(V.59)

Undoing our rescaling (V.25), we see that the true $y$-component of the
velocity is given by

$v_{y}(y)=\pm{v_{s}\over\sqrt{2}\cosh^{2}\left({y-y_{0}\over\ell\sqrt{2}}\right)}\,,$
(V.60)

which shows that the velocity is localized to within a distance $\ell$ of the
interface, and that its characteristic velocity scale $v_{s}$ vanishes like
$m$ as the critical point is approached. The density and the velocity profiles
(V.56) and (V.59) are plotted in figure 15.

Figure 15:  Interface density and velocity profile, as described by plus
variants of (V.56) and (V.59), where $s_{0}=0$. The blue line is the density
fluctuation. The dashed orange line is the velocity field. Boids accelerate
from a low density region (on the left) and move rightward to the high density
region. They reach their maximum velocity halfway through the interface and
begin to decelerate as they move deeper into the plateau. This accumulation of
boids is balanced with diffusion from the high density region to the low
density region.

The steady state profile of the interface is stabilized by a balance between
the velocity field (V.60), which convectively carries boids from the low
density region to the high density region, and diffusion (controlled by
$D_{\rho\perp}$) which carries boids from the high density region to the low
density region.

Finally, we note that the $x$-component of the velocity must also become
position dependent; in particular, it must take on different values in the
plateaus. This is true simply because, in general, the mean local speed of the
flock should depend on the local density. Since that density is different in
the two plateaus, the speeds in each plateau should be different as well.

This can be verified by looking at the $x$-component of the Toner-Tu-Keller-
Segel equations (III.1-III.3) in each plateau. In a plateau, all spatial
derivatives must vanish. In addition, in the steady state, time derivatives
must vanish as well. Hence, the $x$-component of the Toner-Tu-Keller-Segel
equations reduces to

$U(|\mathbf{v}|,\rho,\eta)=0\,.$ (V.61)

This is an implicit equation relating the speed $|\mathbf{v}|$ to the density
$\rho$ and chemo-attractant concentration $\eta$. Since $v_{y}=0$ in the
plateaus (as we just saw), $|\mathbf{v}|=v_{x}$ there. Hence, $v_{x}$ will be
a function of $\rho$ and $\eta$. And since $\rho$ and $\eta$ take on different
values in the plus and minus plateaus, the speed $v_{x}$ will take on
different values as well.

For example, in models like the Vicsek model, in which the speed is an
increasing function of $\rho$, the birds in the high density plateaus will be
moving (parallel to the plateau boundaries) faster than those in the low
density plateaus. For systems exhibiting jamming-like behavior, in which the
local speed is a decreasing function of the local density, the opposite will
be true.

All of the above predictions for the final steady state of the system are
fully borne out by numerical solution of the time-dependent equations (V.29)
and (V.30) in a periodic space, as illustrated in figure 14, and in the
animation Movie 1, in the supplemental materials.

All of our numerical solutions started with initial conditions in which the
departure from a uniform state was small. Furthermore, we performed the analog
of a “critical quench” in an equilibrium system, by choosing our mean density
to be the critical density $\rho_{c}$.

These simulations all show that at small times, the system quickly forms
multiple well-defined plateaus, with density and velocity profiles that agree
well with our analytic solution. Adjacent plateaus merge quickly, while
distant plateaus do not merge on time scales numerically accessible to our
numerical solution (and patience).

This slowing down is a consequence of the exponential approach to the uniform
plateau state of both the density (as shown by equation (V.56)) and the
velocity (as shown by equation (V.59)). It is only the overlap of these
exponential tails that leads to any interaction at all between neighboring
interfaces. As a result, the attraction between neighboring interfaces falls
off exponentially with their separation. Hence, the time for neighboring
domains of high density to merge also grows exponentially with their
separation, which is why the merger slows down so dramatically once the near-
lying plateaus have merged.

We believe that this slowing down is an artifact of both the one-
dimensionality of our numerical solutions, and the fact that they are
noiseless. As discussed in the introduction, we expect that what will really
happen, once the plateaus form, is that the bands will begin to undulate due
to noise in the system, as illustrated in figure 5. These undulations will
grow with time until bands begin bumping into their neighbors, at which point
the bands can start to merge. The merged region will then rapidly “zipper”
along in both directions, merging the two bands. This process will then repeat
with the larger bands formed as a result of such mergers, until full phase
separation is achieved.

Again reiterating our discussion in the introduction, we expect the size of
these undulations to grow algebraically with time, as the interface
fluctuations in the KPZ equation doKardar  _et al._ (1986). Therefore, the
time scale needed for undulations of the bands to grow large enough to reach
neighboring bands will grow only algebraically with the separation between
bands, rather than exponentially. We therefore expect that the time required
for our autochemotaxic instability to completely phase separate a large system
will grow only algebraically with system size.

To make these ideas precise, a theory of the dynamics of fluctuations of an
interface between two flocks of different densities and speeds. This would be
an excellent topic for further research.

Indeed, a quantitative comparison of our analytical result, (V.56) and
(LABEL:wallreal), for the shape of the interfaces and plateau densities fits
the final steady state of these numerical solutions extremely well.

###  V.4 Two-phase coexistence boundary further from the critical point: the
Uncommon Tangent Construction

The previous section’s analysis neglected all non-linearities in the equations
of motion except for those in the isotropic pressure $P_{1}$. As we saw, this
reduces the statics of the problem to those of the equilibrium system.

While, as we also saw in the previous section, it is asymptotically exact to
neglect those other non-linearities as we approach the critical point, it must
also be the case that all non-equilibrium effects, at least on the plateau
densities and the binodal and spinodal curves, must come from those neglected
non-linearities. To see those effects, we must therefore go beyond the leading
order behavior near the critical point, and include the extra non-linearities
$\Lambda$ and $\Gamma$ in equations (V.29) and (V.30).

Recall that we are always limited in any event to analyzing the system close
to the critical point, just as one is in equilibrium, because only near the
critical point are the density variations small enough to justify our
expansion (V.9) of the isotropic pressure. Therefore, we will continue to
assume that we’re close to the critical point, where $\Lambda$ and $\Gamma$
are small. This justifies treating $\Lambda$ and $\Gamma$ as small
perturbations to the problem, and calculating their effects perturbatively.

We will now do such a perturbation theory, and show that the effect of those
terms is to change the common tangent construction of the previous section to
an uncommon tangent construction: the condition that the slopes of our
effective free energy at the two plateau densities be equal survives, but the
tangent at either of these no longer intersects the free energy curve at the
other.

A very similar “uncommon” tangent construction was foundCates and Tailleur
(2015) for “motility induced phase separation”, which occurs in some
disordered active systems.

In the steady state, our full equations (V.29) and (V.30) can be rewritten:

$\displaystyle{d^{2}\tilde{v}\over ds}$ $\displaystyle={d\over
ds}\bigg{(}\partial_{\tilde{\rho}}\tilde{F}(\tilde{\rho})\bigg{)}+{\Lambda\over
2}{d\tilde{v}^{2}\over ds}\,,$ (V.62) $\displaystyle{d^{2}\tilde{\rho}\over
ds^{2}}$ $\displaystyle={d\tilde{v}\over ds}+\Gamma{d\over
ds}(\tilde{\rho}\tilde{v})\,,$ (V.63)

where $\tilde{F}(\tilde{\rho})$ is the effective free energy we introduced in
equation (V.33).

Just as in the common tangent construction, we take the derivative of the
second equation (V.63), and substitute the first equation (V.63) into it. This
gives

$\displaystyle{d^{3}\tilde{\rho}\over ds^{3}}$ $\displaystyle={d\over
ds}\bigg{(}\partial_{\tilde{\rho}}\tilde{F}+\frac{\Lambda}{2}\tilde{v}^{2}\bigg{)}+\Gamma{d^{2}\over
ds^{2}}(\tilde{\rho}\tilde{v})\,.$ (V.64)

Next, we rearrange this equation, keeping the new nonlinearities on the right
hand side. Each term is a derivative with respect to $s$, so we’ll integrate
once, obtaining

$\displaystyle\bigg{[}{d^{2}\tilde{\rho}\over
ds^{2}}-\partial_{\tilde{\rho}}\tilde{F}(\tilde{\rho})\bigg{]}={\frac{\Lambda}{2}\tilde{v}^{2}}+{d\over
ds}\bigg{(}\Gamma(\tilde{\rho}\tilde{v})\bigg{)}+\mu\,,$ (V.65)

where $\mu$ is a constant of integration. Next, we multiply each side by
${d\tilde{\rho}\over ds}$. This enables us to rewrite (V.65) in the form

$\displaystyle{d\over ds}\bigg{[}\frac{1}{2}\bigg{(}{d\tilde{\rho}\over
ds}\bigg{)}^{2}-G(\tilde{\rho})\bigg{]}={\frac{\Lambda}{2}\bigg{(}{d\tilde{\rho}\over
ds}\bigg{)}\tilde{v}^{2}}+\Gamma\bigg{(}{d\tilde{\rho}\over ds}\bigg{)}{d\over
ds}\bigg{(}\tilde{\rho}\tilde{v}\bigg{)}\,,$ (V.66)

where $G$ is our “pseudo-Gibbs free energy”, given by equation (V.43); that is
$G(\tilde{\rho})=\tilde{F}(\tilde{\rho})+\mu\tilde{\rho}$.

Figure 16: Illustration of the uncommon tangent construction. The slope of the
pseudo-Free energy $F(\rho)$ is the same at both plateau densities
$\rho_{\pm}$, but the two tangents are distinct, rather than being the same
line, as they are in equilibrium. It is only in the limit as the system
approaches the critical point that the common tangent construction is
recovered.

Note that if we set the non-linear coefficients $\Lambda$ and $\Gamma$ to
zero, this recovers our earlier result that the term in square brackets on the
left hand side of (V.66) is a constant. That constancy then leads directly to
the common tangent construction of the last section.

Now we wish to work to next order in the non-linear coefficients $\Lambda$ and
$\Gamma$. We’ll do this by writing

$\displaystyle\tilde{\rho}$ $\displaystyle=\tilde{\rho}_{1}+\tilde{\rho}_{2},$
(V.67) $\displaystyle\tilde{v}$ $\displaystyle=\tilde{v}_{1}+\tilde{v}_{2}.$
(V.68)

where $\tilde{\rho}_{1}$ and $\tilde{v}_{1}$ are the common tangent solutions
(V.56) and (V.59) found in the previous section, and $\tilde{\rho}_{2}$ and
$\tilde{v}_{2}$ are small perturbations to those solutions; that is

$\tilde{\rho}_{1}=\pm\tanh\left({s-s_{0}\over\sqrt{2}}\right),\,$ (V.69)
$\tilde{v}_{1}=\pm{1\over\sqrt{2}\cosh^{2}\left({s-s_{0}\over\sqrt{2}}\right)}.\,$
(V.70)

We’ll assume, and verify a posteriori, that $\tilde{\rho}_{2}$ and
$\tilde{v}_{2}$ are of order $\Lambda$ and $\Gamma$; $\tilde{\rho}_{1}$ and
$\tilde{v}_{1}$ are, obviously, zeroeth order in $\Lambda$ and $\Gamma$.

In light of this perturbative hierarchy, we can, to next order in $\Lambda$
and $\Gamma$, replace $\tilde{\rho}$ and $\tilde{v}$ with $\tilde{\rho}_{1}$
and $\tilde{v}_{1}$ on the right hand side of (V.66). Doing so, using
$\tilde{v}_{1}=\partial_{s}\tilde{\rho}_{1}$, and integrating over $s$ from a
low density plateau to a high density one, and remembering once again that
within a plateau, ${d\tilde{\rho}\over ds}=0$, gives

$\displaystyle G(\tilde{\rho}_{+})-G(\tilde{\rho}_{-})$ $\displaystyle=$
$\displaystyle-\frac{\Lambda}{2}{\int_{-\infty}^{\infty}}\bigg{(}{d\tilde{\rho}_{1}\over
ds}\bigg{)}^{3}ds-\Gamma{\int_{-\infty}^{\infty}}\bigg{(}{d\tilde{\rho}_{1}\over
ds}\bigg{)}{d\over ds}\bigg{(}\tilde{\rho}_{1}\left({d\tilde{\rho}_{1}\over
ds}\right)\bigg{)}ds\,.$ (V.71)

Plugging (V.69) into this expression, and performing the resulting elementary
integrals, we find:

$\displaystyle G(\tilde{\rho}_{+})-G(\tilde{\rho}_{-})$
$\displaystyle=2\tilde{P}_{0}+2\mu=-\frac{4}{15}[\Lambda+\Gamma]$ (V.72)
$\displaystyle\mu$
$\displaystyle=-\tilde{P}_{0}-\frac{2}{15}[\Lambda+\Gamma].$ (V.73)

Since we are working in the limit in which the nonlinear coefficients are all
small, we see that (V.73) implies that $\mu+\tilde{P}_{0}$ is small. We will
make use of this fact in a moment.

Next, we evaluate equation (V.65) deep in the density plateaus, where all
spatial derivatives vanish. This implies:

$\displaystyle-\frac{\partial\tilde{F}(\tilde{\rho}_{\pm})}{\partial\tilde{\rho}}=-\tilde{P}_{0}+\rho_{\pm}-\rho_{\pm}^{3}=\mu\,.$
(V.74)

We define $\delta_{\pm}=\tilde{\rho}_{\pm}\mp 1$ as the small deviation of
these densities from those we obtained earlier by ignoring these “new” non-
linearities. From equation (V.73), we see that these obey

$\displaystyle-(\pm 1+\delta_{\pm})+(\pm 1+\delta_{\pm})^{3}$
$\displaystyle=-\mu-\tilde{P}_{0}\,,$ (V.75)

which is easily solved to linear order in the small quantities $\delta_{\pm}$
to give

$\displaystyle\delta_{\pm}=-\left(\frac{\mu+\tilde{P}_{0}}{2}\right)=\frac{1}{15}[\Lambda+\Gamma]\,.$
(V.76)

Counter-intuitively, the shift is the same for the low and high density
plateaus. We strongly suspect that this does not hold to higher orders in
perturbation theory, although we have not attempted such a calculation.

Figure 17 shows the excellent agreement between the analytically predicted
offset (V.76) and the results of simulating the time dependent equation of
motion for a long time.

An uncommon tangent construction much like that we find here is also found in
phase separation in disordered active systems (e.g., “Motility Induced Phase
Separation” (“MIPS”) Cates and Tailleur (2015)). It thus appears that this
uncommon construction is, paradoxically, quite common in active systems.

Figure 17: The result of evolving equations (V.62) and (V.63) for a long time
in a periodic space. The dashed blue curve is the initial state of the
numerical solution. The orange curve is the density profile once the final
time step is reached. We retain the plateau structure but it is shifted
upwards. The predicted offset $\delta_{\pm}$ differs from the result obtained
by direct numerical solution of the equations of motion, (V.16) and (V.17), by
less than 6%. The parameters used for this numerical solution were
$\Gamma,\Lambda=0.25$ and it was iterated for 3000 units of time.

##  VI Conclusion

We have shown that a dry active autochemotaxic flock, or, more generally, any
flock with sufficiently strong attractive interactions between the flockers,
can become unstable to the formation of density bands parallel to the mean
direction of motion of the flock. The attractive interactions must drive the
inverse compressibiliy negative. This behavior may be connected to the
formation of, e.g., unidirectional ant trails. The instability is caused by
the autochemotaxic “pressure” or other attractive mechanism overcoming the
usual repulsive isotropic pressure, causing the inverse compressibility $B$ of
the flock to be negative. We investigated that instability, and specifically
demonstrated that:

  * •

The instability is anisotropic. The direction of maximum growth rate is
perpendicular to the mean flock motion. This initially causes the formation of
a growing wavelike modulation of the local density parallel to mean flock
motion. The maximum growth rate ${\rm Im}~{}\omega_{max}\propto\sqrt{|B|}$.

  * •

This modulation quickly grows into a “phase separated” flock consisting of
high density plateaus separated by low density plateaus. The boundaries of the