# Thermodynamically consistent flocking: From discontinuous to continuous
transitions

Tal Agranov  DAMTP, Centre for Mathematical Sciences, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, United Kingdom  Robert L. Jack  DAMTP,
Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road,
Cambridge CB3 0WA, United Kingdom  Yusuf Hamied Department of Chemistry,
University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom
Michael E. Cates  DAMTP, Centre for Mathematical Sciences, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom  Étienne Fodor
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg, Luxembourg

###### Abstract

We introduce a family of lattice-gas models of flocking, whose
thermodynamically consistent dynamics admits a proper equilibrium limit at
vanishing self-propulsion. These models are amenable to an exact coarse-
graining which allows us to study their hydrodynamic behavior analytically. We
show that the equilibrium limit here belongs to the universality class of
Model C, and that it generically exhibits tricritical behavior. Self-
propulsion has a non-perturbative effect on the phase diagram, yielding novel
phase behaviors depending on the type of aligning interactions. For aligning
interaction that increase monotonically with the density, the tricritical
point diverges to infinite density reproducing the standard scenario of a
discontinuous flocking transition accompanied by traveling bands. In contrast,
for models where the aligning interaction is non-monotonic in density, the
system can exhibit either (the nonequilibrium counterpart of) an azeotropic
point, associated with a continuous flocking transition, or a state where
three phases coexist dynamically.

##  1 Introduction

Flocking is a prominent phenomenon in nonequilibrium statistical mechanics [1,
2], comprising the collective motion of a large group of aligning self-
propelled agents. It appears in systems ranging from bird flocks [3] to human
crowds [4] and synthetic self-propelled colloids [5]. Various theoretical
models lead to a qualitatively similar flocking behavior, which suggests
possible universality. A pioneering theoretical account of flocking is due to
Vicsek et al. [6], who proposed an agent-based model almost three decades ago,
followed by Toner and Tu [7] who offered a hydrodynamic description. Since
then, these models and their variants have played a leading role in the study
of active matter [8].

In the Vicsek model, locally aligning spins move with fixed speed and
individually fluctuating orientation, which drives the dynamics out of
equilibrium. Strikingly, under local alignment, this model undergoes an
ordering transition even in two dimensions, which would be precluded in
equilibrium by the Mermin-Wagner theorem [9, 10]. When increasing either
density or alignment, the system transitions from a homogeneous disordered gas
to a homogeneous collective motion with long-ranged order. Originally thought
of as a continuous transition [6], it took almost another decade to establish
flocking as a discontinuous transition [11]: in between the two homogeneous
phases, there is a region in parameter space where dense ordered bands
propagate against a dilute disordered background. Unlike equilibrium liquid-
gas transitions, the critical point is at infinite density. This feature has
been attributed to the fact that the dense and dilute phases have different
symmetries, so that they cannot be connected via a continuous transformation
[12].

A recent lattice-based model, known as the active Ising model (AIM), has led
to a series of analytical predictions for the phase diagram of the flocking
transition [13, 14]. In this simple extension of the Ising model, motile spins
with discrete symmetry can only self-propel along a single axis. Importantly,
this model is amenable to coarse-graining, yielding an exact hydrodynamic
description [12, 15, 16]. Although the phenomenology of AIM is clearly
distinct from any equilibrium analogue, it is tempting to try and rationalize
how AIM departs from an equilibrium reference model, with proper reversible
dynamics, as was done for repulsive active particles [17, 18]. Interestingly,
AIM does not reduce to equilibrium at vanishing self-propulsion, as one might
expect. This inconsistency reveals that AIM is actually not compatible with
the theory of stochastic thermodynamics [19], which assumes that the system is
in contact with equilibrium reservoirs and driven by external forces (so that
it naturally relaxes to equilibrium in the absence of drive). The same holds
for other flocking models with continuous spin symmetry [6], or with
topological interactions that involve nonreciprocal alignment [20, 21], which
also lead to irreversibility even at vanishing self-propulsion.

It is then natural to examine flocking models where the microscopic dynamics
is amended to entail an equilibrium limit in the absence of self-propulsion,
which we refer to as thermodynamic consistency. Several questions remain open:
Would novel qualitative behavior emerge from thermodynamic consistency? To
which universality class does the equilibrium limit of such flocking models
belong? Inspired by AIM, we explore these issues in the context of flocking
systems with a discrete symmetry. To make analytical progress, we choose
lattice models which, like AIM, admit an analytically exact and tractable
hydrodynamic description.

Although our models operate far from thermodynamic equilibrium, we find that
their behavior can be usefully analyzed in terms of the underlying free energy
that controls their equilibrium limit. As noted previously [22, 13, 12], the
standard scenario of a discontinuous flocking transition stems from a density
modulated coupling of the alignment interaction. In our models, we reveal a
strikingly different phenomenology between the cases of monotonic and non-
monotonic density modulations. For models where ordering increases
monotonically with density, the equilibrium limit exhibits a single
tricritical point [16], and introducing activity is a non-perturbative effect:
the tricritical point diverges to infinite density, and the flocking
transition is discontinuous all across phase space. In contrast, for models
where ordering is non-monotonic with density, the equilibrium limit exhibits a
pair of tricritical points. At finite self-propulsion, such points can either
collide into a single azeotropic point, yielding a continuous flocking
transition, or instead a new collective state emerges where three phases
coexist.

This paper is organized as follows. In Sec. 2, we set up the general
definition of our thermodynamically consistent flocking models. We propose a
microscopic dynamics with a generic type of aligning interactions, and derive
the corresponding hydrodynamic description through systematic coarse-graining.
Then, we examine in Secs. 3 and 4 the phase diagrams for monotonic and non-
monotonic dependence of the aligning strength on density, respectively.
Overall, our results illustrate how various choices of density-dependent
alignment affect the flocking transition, yielding unprecedented phenomenology
in the flocking context. Yet, some of this phenomenology echoes that of
equilibrium systems, so that our thermodynamically consistent description
allows one to anticipate and classify the topologies of nonequilibrium phase
diagrams in terms of their equilibrium counterparts.

##  2 Thermodynamic consistency: From microscopics to hydrodynamics

In this Section, we introduce a class of lattice dynamics for active spins
that obey thermodynamic consistency. We take locally aligning interactions
with mesoscopic range to capture a transition to collective motion. Such a
form of interactions is amenable to exact coarse-graining, which allows us to
derive the corresponding hydrodynamic equations for some relevant fields.
Importantly, our derivation is given in a functional form with an arbitrary
dependence of the alignment strength on density.

###  2.1 Microscopic dynamics of aligning active spins

Flocks with discrete symmetry, like in AIM [12], are defined by two species of
particles ($+$ and $-$) which hop with a bias to the right or left
(respectively) on a periodic lattice of $L\gg 1$ sites, see Fig. 1. We
consider for simplicity a one dimensional lattice, yet our fluctuating
hydrodynamics is trivially extended to higher dimensions. The dynamics does
not exclude multiple occupancy, so that each site $i$ can have any number of
particles. The total number of particles in the system is $N=\rho_{0}L$, where
$\rho_{0}$ is the mean density.

The dynamical update rules consist of diffusive hops biased by self-
propulsion, and spin flips $+\leftrightarrow-$. Differently from previous
models [12, 16], at zero self-propulsion we make both these rules derive from
the same Hamiltonian $H$ via a detailed balance constraint [23, 19]. We leave
$H$ unspecified at this stage. Self-propulsion is then added as a weak bias of
hops in the direction of the particle’s spin, yielding the following update
rules (see Fig. 1):

  1. 1.

Site hoping: Any particle jumps to a neighboring site with rate
$D_{0}e^{-\frac{\beta\Delta H}{2}+\frac{\lambda}{D_{0}L}}$ if the jump is in
the direction of its spin, and $D_{0}e^{-\frac{\beta\Delta H}{2}}$ otherwise.
Here, $\Delta H$ is the energy difference between the configurations of the
system before and after the jump.

  2. 2.

Tumbling: A $+$ ($-$) particle converts into a $-$ ($+$) particle with rate
$(\gamma/L^{2})e^{-\frac{\beta\Delta H}{2}}$.

Here, $D_{0}$ is the bare diffusion constant, $\lambda$ is the self-propulsion
strength, and $\gamma$ sets the tumble rate. The scaling of the rates with $L$
ensures that in the hydrodynamic limit ($L\to\infty$ at fixed $\rho_{0}$) all
processes occur on diffusive time scales [15, 24]. Indeed, the time it takes
for particles to traverse a macroscopic system of size $L$, either via
diffusive motion or on account of self-propulsion, scales as $L^{2}$ which is
also the time scale for tumbling events. This $L$-dependent scaling of the
rates is standard for other lattice-based diffusive models [15, 24]. The main
difference with AIM [13, 12, 22] is that our site hopping now accounts for the
energy difference $\Delta H$. This feature renders our models reversible at
vanishing self-propulsion, thus setting a proper equilibrium limit, as a
consequence of thermodynamic consistency [23, 19].

Figure 1: The active flocking model with thermodynamic consistency. Each
particle can be in either one of two states ($+$ and $-$) which determine the
direction of biased diffusion (respectively to the right and left). The
aligning Hamiltonian $H$ constrains both the change of particle states and
their diffusive hops.

We now restrict $H$ such that the interaction range spans a mesoscopic scale
$\Delta x=L^{\delta}$, with $0<\delta<1$. This scale is sub-extensive in
system size, yet contains a large number of particles, and is the scale
appearing in the coarse-graining procedure detailed below. Importantly, it
does not amount to a simple mean-field limit, for which the interaction range
would span the entire system ($\delta=1$), and it enables a phase transition
in one dimensision even in the equilibrium limit111In defiance of the usual
Landau argument, one has that the interfacial energy scale at least as
$L^{\delta}$, overcoming the entropic $\log L$ contribution..

In practice, interactions are given in terms of some local mean occupancies
for the density and ‘magnetization’ variables, respectively
$\eta_{i}^{\rho}=\eta_{i}^{+}+\eta_{i}^{-}$ and
$\eta_{i}^{m}=\eta_{i}^{+}-\eta_{i}^{-}$, where $\eta_{i}^{+}$ and
$\eta_{i}^{-}$ are the number of $+$ and $-$ particles at site $i$. The local
empirical means around site $i$ are given by

$\hat{\rho}_{i}=\frac{\sum_{|i-j|<\Delta x}\eta_{j}^{\rho}}{2\Delta
x},\quad\hat{m}_{i}=\frac{\sum_{|i-j|<\Delta x}\eta_{j}^{m}}{2\Delta x}.$ (1)

We then choose the aligning Hamiltonian to take the form

$H=-\sum_{i=1}^{L}\frac{\hat{m}_{i}^{2}}{2\hat{\rho}_{i}}f(\hat{\rho}_{i}).$
(2)

For $f(\hat{\rho}_{i})=1$, this can be viewed as a fully connected Ising
Hamiltonian within some restricted range $\Delta x$, where the coupling is
scaled with the total number of interacting spins (as encoded in the
denominator $\hat{\rho_{i}}$). The alignment strength is then set by the
density-dependent modulation $f(\hat{\rho}_{i})$, which plays a key role in
the following. For AIM, a similar modulation appears at a coarse-grained
level, as an effective renormalization due to the hydrodynamic fluctuations
[22]. In our model, since the microscopic dynamics already contains a
mesoscale interaction range, we choose to directly introduce such a modulation
at this level. Various forms of $f$ can be motivated on account of different
types of close neighbor interactions [25, 26, 27, 28].

In the following, we consider the hydrodynamic limit of this model. As a
preliminary step, note that for a single isolated particle, there are two
important length scales: $\xi_{p}=\lambda L/\gamma$ is the average
displacement of particles between tumbles due to self-propulsion, and
$\xi_{d}=L\sqrt{D_{0}/\gamma}$ is the typical diffusive displacement, over the
same time. These distances are measured in units of the lattice spacing and
they are both $\mathcal{O}(L)$, because of the $L$-dependent hopping rates of
the model.

###  2.2 Hydrodynamic equations for density and magnetization

We follow a standard coarse-graining procedure, inspired by similar lattice
models with diffusive scaling [29], which has already been deployed in some
active lattice gas models [30, 31]. To this end, we embed the $L$ sites of the
lattice into the interval $[0,\ell_{s}]$ where
$\ell_{s}=L/\xi_{d}=\sqrt{\gamma/D_{0}}$ is the system size in the
hydrodynamic representation (measured in units of $\xi_{d}$, as in [15]). Then
introduce the hydrodynamic spatial coordinate $x=i\ell_{s}/L\in[0,\ell_{s}]$.
The reason for these choices is that we will later analyse phase-separated
profiles: their interfacial widths are of order 1 in this hydrodynamic
representation and it will be convenient to take $\ell_{s}\gg 1$ so that the
systems is large compared to this interfacial width. For consistency with
these choices we make a diffusive scaling of time, setting
$t=\gamma\hat{t}/L^{2}$ with $\hat{t}$ the time variable of the microscopic
system.

The local density in the vicinity of point $x$ is defined by the mesoscopic
coarse-graining (1) as:

$\rho(x,t)=\hat{\rho}_{\frac{xL}{\ell_{s}}}\left(t\right),\quad
m(x,t)=\hat{m}_{\frac{xL}{\ell_{s}}}\left(t\right).$ (3)

Since $\rho$ is conserved and $m$ is not, the dynamics of these fields must
take the general form

$\begin{bmatrix}\partial_{t}{\rho}\\\
\partial_{t}m\end{bmatrix}=-\partial_{x}\begin{bmatrix}J_{\rho}\\\
J_{m}\end{bmatrix}-\begin{bmatrix}0\\\ 2K\end{bmatrix}$ (4)

where $J_{\rho},J_{m}$ are the conservative fluxes and $K$ is the net rate of
change of the $\rho_{-}$ density due to spin flipping. As detailed in A, to
leading order in system size we find the following expressions

$\begin{bmatrix}J_{\rho}\\\
J_{m}\end{bmatrix}=-\frac{1}{2}\mathbb{C}(\rho,m)\partial_{x}\begin{bmatrix}\frac{\delta
F}{\delta\rho}\\\ \frac{\delta F}{\delta
m}\end{bmatrix}+\text{Pe}\begin{bmatrix}m\\\ \rho\end{bmatrix},\quad
K=M(\rho,m)\sinh\left(\frac{\delta F}{\delta
m}\right),\quad\text{Pe}=\frac{\lambda}{\sqrt{D_{0}\gamma}}.$ (5)

The Péclet number $\text{Pe}=\xi_{p}/\xi_{d}$ is the standard measure of
activity in similar active systems with diffusive scaling [32]. Here,
$F=\frac{L}{\ell_{s}}\int_{0}^{\ell_{s}}dx\mathcal{F}(\rho,m)$ is the free-
energy functional of the equilibrium ($\lambda=0$) system with density

$\mathcal{F}\left(\rho,m\right)=\frac{\rho}{2}\ln\frac{\rho^{2}-m^{2}}{4}+\frac{m}{2}\ln\frac{\rho+m}{\rho-m}-\rho-\beta\frac{m^{2}}{2\rho}f\left(\rho\right),$
(6)

while $\mathbb{C}$ and $M$ are mobility coefficients given by

$\displaystyle\mathbb{C}=2\begin{bmatrix}\rho&m\\\ m&\rho\end{bmatrix},\quad
M=\sqrt{\rho^{2}-m^{2}}.$ (7)

The expressions (5-6) are exact as long as $\text{Pe}\neq 0$. Yet in order to
describe correctly the equilibrium ($\text{Pe}=0$) hydrodynamics for non
convex $\mathcal{F}$, one has to retain surface tension terms that are sub-
leading $\mathcal{O}\left(L^{-2}|\nabla\rho|^{2},L^{-2}|\nabla m|^{2}\right)$.
Still, as is usual in the thermodynamics of phase separation, the bulk term
$\mathcal{F}$ suffices for determining the binodal (coexisting) densities.

The dynamics (4-5) has a thermodynamically consistent structure. Indeed, for
$\text{Pe}=0$, all terms derive from the free energy $\mathcal{F}$, which then
serves as a Lyapunov function, as discussed in A. As shown in what follows,
this equilibrium limit belongs to the universality class of Model C [16, 33].
Note that, as written, the dynamics (4-5) neglects some sub-leading noise
terms of order $L^{-1/2}$; these can be readily deduced from the detailed-
balance condition (and are included by definition in Model C).

The general structure of (4-5) should hold rather broadly for all
thermodynamically consistent flocking models with discrete symmetry and under
diffusive scaling, including models with purely local interaction range
(albeit with a renormalized free energy depending on dimensionality). The
conservative fluxes in (5) are linear in the free-energy derivatives, as is
common when expanding close to equilibrium. That these expressions hold also
far away from equilibrium is a consequence of the $L$-dependent scaling of the
microscopic rates, and the diffusive scaling of space and time; any higher
order spatial gradients will be sub-leading in system size $L^{-1}$. This
reasoning does not apply for the flipping term $K$, since it does not involve
spatial variations of the fields. Indeed, this term goes beyond linear order
in the free energy derivatives even at the hydrodynamic level [34]. Activity
($\text{Pe}>0$) leads to the extra terms $\text{Pe}\,m$ and $\text{Pe}\,\rho$
in the conservative fluxes (5) but otherwise does not interfere with the
remaining bare equilibrium terms. Again, this is a consequence of the
diffusive scaling, where self-propulsion enters as a weak bias in the
microscopic dynamics.

Our thermodynamically consistent dynamics in (4-5) is easily related to the
non-consistent one of AIM [13, 12, 22]. Our flipping term $K$ is the same as
in the hydrodynamic equations of AIM, under appropriate matching of $f$. The
conservative fluxes $J_{\rho}$ and $J_{m}$ are different from those of AIM,
since they stem from a different choice of hopping rules. Nevertheless, one
can reproduce the corresponding terms in AIM by sending the temperature that
enters these terms to infinity (namely, taking $\beta=0$ in the free energy
that enter the conservative fluxes $J_{\rho}$ and $J_{m}$), as discussed in A.
To this extent, the source of residual irreversibly in the zero self-
propulsion limit of AIM can be regarded as a mismatch between two temperatures
controlling respectively the conservative fluxes $(J_{\rho},J_{m})$ and the
flipping term $K$.

##  3 Phase diagram for monotonic alignment strength

Various forms of the density-dependent alignment strength
$f\left(\rho\right)$, entering in the microscopic alignment (2) and in the
hydrodynamic free energy (6), can be motivated to reflect different types of
close neighbor interactions [25, 26, 27, 28]. We first explore the phase
behavior for a monotonic $f(\rho)$ by considering the form appearing
previously in AIM [13, 12, 22]:

$f\left(\rho\right)=1-\frac{r}{\rho}.$ (8)

This choice is convenient mathematically, although on physical grounds one
might instead favor choices that are positive everywhere, and do not diverge
at $\rho\to 0$. Yet this does not change the phenomenology reported below.
Other monotonic forms have been recently proposed in extensions of AIM
yielding a similar phase behavior [16]. In some of these models, the alignment
term is not purely quadratic in the magnetization as assumed here (see (2)).
Yet, its expansion close to criticality reduces to this form, and the topology
of the phase diagram actually follows from the behavior close to criticality.

The resulting phase diagram, given in terms of $T=1/\beta$ and $\rho_{0}$, is
composed of three phases shown in Fig. 2: (i) a homogeneous disordered phase
(H; $m=0$), (ii) a homogeneous ordered phase (H; $m\neq 0$), which becomes a
collective motion (C.M.) for non-vanishing self-propulsion, and (iii) a phase-
separated state (P.S.) between different discrete symmetries ($m=0$ and $m\neq
0$), which becomes a traveling band (T.B.) for non-vanishing self-propulsion.
In this Section, we study the topology of this phase diagram and derive its
phase boundaries. We address first the equilibrium limit ($\text{Pe}=0$), and
then the nonequilibrium regime ($\text{Pe}>0$).

Figure 2: Phase diagrams for (a) the equilibrium case ($\text{Pe}=0$) and (b)
the nonequilibrium regime ($\text{Pe}=0.5$), with monotonic $f$ given by (8)
with $r=1$. The colored dashed and solid curves are spinodals and binodals,
respectively. The black solid curve is a critical line. The black dot marks a
tricritical point. (c) Traveling band (T.B.) profiles for $\text{Pe}=0.5$,
$\rho_{0}\simeq 3.62$, and $T\simeq 0.71$.

###  3.1 The equilibrium limit ($\text{Pe}=0$): Model C dynamics

As pointed out in [16], the equilibrium, zero-propulsion limit of
thermodynamically consistent flocks with discrete symmetry must generically
lie within the universality class of Model C [33]. In our model, this should
remain true for any non-trivial choice of $f(\rho)$, while allowing different
global phase behaviors even at equilibrium. Model C describes an equilibrium
dynamics which couples a conserved and a non-conserved scalar order parameter,
and the canonical Landau form close to criticality is given by [35, 36]

$\mathcal{F}_{\text{cri}}=\frac{a}{2}\left(\delta\rho\right)^{2}+\frac{b}{4}\left(\delta\rho\right)^{4}+\frac{c}{2}m^{2}+\frac{d}{4}m^{4}+e\,m^{2}\delta\rho,$
(9)

where higher order terms and terms purely linear in the density are omitted.
Expanding our free energy in (6), one can derive the expressions of the
different coefficients in (9), which are determined by the choice of
$f\left(\rho\right)$. Importantly, the coupling coefficient $e$ in (9) reads

$e=\beta f^{\prime}\left(\rho\right),$ (10)

which vanishes for a constant $f(\rho)$; in this special case, our model falls
outside the Model C universality class. (As discussed in Sec. 4, this also
makes the nonequilibrium $\text{Pe}>0$ regime different from the usual
flocking scenario.) As detailed in [36], the topology of the phase diagram is
set by the values of the coefficients in (9). As we have the exact form of the
free energy at our disposal (6), we follow the procedure presented in [36] to
map the phase diagram across the whole phase space, even beyond criticality.

We find the stable phases by minimizing the free energy (6) under the
constraint of conserved density $(1/\ell_{s})\int_{0}^{\ell_{s}}\rho
dx=\rho_{0}$. Since magnetization is not conserved, then its value is locked
to the local value of the density via the minimization of $F$ with respect to
the magnetization. A non-trivial minimizer $m_{0}\neq 0$ then emerges at
$\left.\partial_{mm}\mathcal{F}\right|_{m=0}=0$ which defines the curve

$T_{c}=f\left(\rho\right).$ (11)

Generally, the relation (11) only establishes linear instability of the
neutral state ($m=0$). In regions that also admit phase separation, as
discussed below, this curve will then mark a spinodal instability which we
will term the magnetization spinodal, see dashed blue curve in Fig. 2(a). Only
in the absence of phase separation, this curve marks a second-order transition
between disordered and ordered states, see solid black curve in Fig. 2(a). In
the ordered phase, the value $m_{0}\left(\rho,T\right)$ is found from
$\partial_{m}\mathcal{F}=0$, which gives the magnetization in an implicit
form:

$m_{0}=\rho x\left(\rho,T\right),\quad\frac{\text{arctanh}x}{x}=\beta
f\left(\rho\right),$ (12)

where the second equality defines the dependence of $x$ on $\rho$ and $T$. A
state with magnetization (12) minimizes the free energy only for homogeneous
profiles. To find phase separation, one has to follow the usual procedure of a
common-tangent construction over the single-variable function
$\mathcal{F}\left[\rho,m_{0}\left(\rho,T\right)\right]$ [36]. Non-convexity of
this function marks linear instability against spatial density variations, and
is equivalent to non-convexity of the two-variable free energy $\cal F$:

$\left|\text{Hess}\left(\mathcal{F}\right)\right|<0,$ (13)

where $\left|\text{Hess}\left(\mathcal{F}\right)\right|$ is the Hessian
determinant of $\mathcal{F}\left(\rho,m\right)$. An explicit expression for
(13) is presented in B.

The relation (13) defines a spinodal curve that is complementary to the
magnetization spinodal (11), see dashed red line Fig. 2(a), which we will term
the density spinodal (since instability here also involves density
modulations). The region enclosed between these two curves defines the region
of spinodal decomposition. The two spinodal curves meet at a tricritical
point, as detailed in B, whose location obeys for the density
$\rho_{\text{tri}}$

$\left[\rho_{\text{tri}}\frac{f^{\prime}\left(\rho_{\text{tri}}\right)}{f\left(\rho_{\text{tri}}\right)}\right]^{2}=\frac{2}{3}.$
(14)

Any monotonic $f(\rho)$ that vanishes at a finite density while saturating at
infinity will have at least one root of (14), yielding a phase diagram like
that in Fig. 2(a). At $T>T_{c}\left(\rho_{\text{tri}}\right)$, the curve (11)
marks a line of critical points separating the homogeneous disordered ($m=0$)
and ordered ($m\neq 0$) phases. For $T<T_{c}\left(\rho_{\text{tri}}\right)$,
the transition to the ordered phase becomes discontinuous with a region of
coexistence between an ordered liquid and a disordered gas. Within the
miscibility gap, the relation (11) marks the magnetization spinodal, and the
complementary density spinodal given by (13) lies within the symmetry-broken
phase. The location of the binodal curves follows from a common tangent
construction, as detailed in B. Note that, although coexisting bulk phases are
found within this procedure, properly resolving the interfaces would
necessitate retaining sub-leading gradient terms, which quantify surface
tension in the free energy. Note also that Model C can exhibit additional
phase diagram topologies, such as critical points and double critical end
points [36], under appropriate tuning of $f(\rho)$. These are absent for the
curves (8), but could arise for other sufficiently elaborate choices,
including strictly monotonic curves.

###  3.2 The nonequilibrium regime ($\rm{Pe}>0$): Discontinuous flocking
transition

The phase diagram for the monotonic $f(\rho)$ obeying (8) at finite
$\text{Pe}=0.5$ is shown in Fig. 2(b). The tricritical point shifts to
infinite density, turning the second-order transition line into a first order
coexistence region (miscibility gap) along its entire length. The same phase
behavior is observed for AIM and its recent variants [13, 16]. Within the
miscibility gap, the phase separated state becomes a traveling band (T.B.), as
shown in Fig. 2(c). Indeed, at non-vanishing self-propulsion, any magnetized
bulk phase attains a non-vanishing density flux given by $\text{Pe}\,m$,
according to (5). Then, simple flux balance across the interface implies that
it must propagate with velocity given by

$V=\text{Pe}\frac{\Delta m}{\Delta\rho},$ (15)

where $\Delta\rho$ and $\Delta m$ are the density and magnetization difference
between phases. Moreover, we have found numerically that the propagation
velocity is bounded below as $|V|\geq\text{Pe}$, for all models that we have
examined.

To establish more broadly the generic topology of the phase diagram in Fig.
2(b), we analyze the linear instability of the dynamics. First, the
magnetization spinodal line (11) remains unchanged even at finite
$\text{Pe}>0$, since it describes instabilities against homogeneous
variations. As the active terms in the dynamics (4-5) only couple to spatial
flux derivatives, they do not affect this instability. In contrast, the
complementary density spinodal (13) is changed. As we show in C the
instability in (13) is now modified into the following criterion

$\left|\text{Hess}\left({\cal
F}\right)\right|<\frac{\text{Pe}^{2}}{2M\left(\rho_{0},m\right)\rho_{0}}\left[\left(\frac{\partial_{\rho
m}\mathcal{F}}{\partial_{mm}\mathcal{F}}\right)^{2}-1\right],$ (16)

where the mobility $M$ is given in (7), and the magnetization $m$ is to be
evaluated at $m=m_{0}(\rho,T)$ in (12).

Based on other flocking models, it was argued that the flocking transition
should generically be discontinuous everywhere [2]. In our model, this holds
true if the conditions in (11) and (16) predict a finite spinodal gap at any
$\text{Pe}>0$. It is useful to notice that, at the magnetized state
$m=m_{0}(\rho,T)$ in (12), we get $\partial_{\rho
m}\mathcal{F}/\partial_{mm}\mathcal{F}=-\partial m_{0}/\partial\rho$, so that
the density spinodal (16) becomes

$\left|\text{Hess}\left({\cal
F}\right)\right|<\frac{\text{Pe}^{2}}{2M\left(\rho_{0},m\right)\rho_{0}}\left[\left(\frac{\partial
m_{0}}{\partial\rho}\right)^{2}-1\right].$ (17)

Moreover, close to the magnetization spinodal (11), we have
$m_{0}\sim\sqrt{T_{c}\left(\rho\right)-T}$ with $T_{c}=f(\rho)$, yielding

$\left.\frac{\partial
m_{0}}{\partial\rho}\right|_{\rho_{0}}\sim\frac{T_{c}^{\prime}\left(\rho_{0}\right)}{\sqrt{T_{c}\left(\rho_{0}\right)-T}}.$
(18)

Given that (18) diverges close to the magnetization spinodal (11), a finite
neighborhood of the magnetization spinodal remains unstable at arbitrary small
Pe no matter where one is on the phase diagram. Therefore, the non-
perturbative effect of activity on the phase diagram can be attributed to the
singularity of the order disorder transition of the equilibrium limit, as
encoded in $m_{0}$. Remarkably, however, this argument fails for any non-
monotonic $f(\rho)$ at its maximum $f^{\prime}(\rho_{0})=0$. As discussed next
in Sec. 4, the case of non-monotonic $f$ actually leads to novel phase
behaviors and a diagram very different from Fig. 2(b).

##  4 Phase diagram for non-monotonic alignment strength

Non-monotonic alignment strength functions $f(\rho)$ are likely to arise
physically in systems where alignment is restricted by visual obstruction, as
was recently proposed to model natural bird flocking [37]. Non-monotonicity
opens the door to novel phase behaviors, not encountered for the monotonic
case in Sec. 3 above. In equilibrium ($\text{Pe}=0$), when $f(\rho)$ features
a single maximum, one typically finds two roots for (14) marking a pair of
tricritical points. This scenario leads to two distinct miscibility gaps, as
shown in Fig. 3(a). The miscibility gap to the left of the maximum has a
disordered low-density phase coexisting with an ordered higher-density phase.
The miscibility gap to the right of the maximum has a disordered high-density
phase coexisting with an ordered lower-density phase.

Figure 3: Phase diagrams for (a) the equilibrium case ($\text{Pe}=0$) and (b)
the nonequilibrium regime ($\text{Pe}=1$), with non-monotonic $f$ given by
(19) with $\alpha=1/20$. The colored dashed and solid curves are spinodals and
binodals, respectively. The black solid curve is a critical line. The black
dots mark tricritical points. The black dashed line marks the maximal
magnetization curve (22). An ‘azeotropic point’ emerges in the nonequilibrium
regime. (c) Reversed band (R.B.) profiles for $\text{Pe}=1$, $\rho_{0}\simeq
15.62$, and $T\simeq 0.43$.

In the rest of this Section, we consider how self-propulsion ($\text{Pe}>0$)
affects this equilibrium picture. The regime of non-monotonic $f$ is more
elaborate compared with the monotonic case in Sec. 3. In particular, we
distinguish two qualitatively different phase behaviors, depending on the
curvature of $f(\rho)$ close to its maximum. For illustration purposes, we use
particular examples in a specific family of non monotonic-curves:

$f\left(\rho\right)=\left(1-\frac{1}{\rho}\right)e^{-\alpha\rho},$ (19)

with $\alpha>0$. In what follows, we show that the topology of the phase
diagram is controlled by local properties of $f$ close to the maximum, so that
the exponential cutoff at large densities is immaterial for our purposes.

###  4.1 Azeotropic point in the nonequilibrium case ($\rm{Pe}>0$)

For mildly curved functions $f(\rho)$, we find numerically that the
tricritical points collide at the maximum of $f$ at any finite Pe, see Fig.
3(b). As for the case of monotonic $f(\rho)$ in Sec. 3, the effect of self-
propulsion on the phase diagram is non-perturbative: the equilibrium critical
line expands almost everywhere into a finite miscibility gap, except at a
single point which remains critical. Thus, in contrast with the case of
monotonic $f(\rho)$, there does remain a critical point on the phase diagram
at finite Pe. Therefore, a non-monotonic $f(\rho)$ can allow for a continuous
transition from homogeneous disorder (H., $m=0$) to collective motion (C.M.).

For equilibrium systems showing a phase-diagram topology equivalent to Fig.
3(b), the critical point is called an azeotropic point [38], so that we also
adopt the same term here. This point marks the pairwise meeting of two
branches of first order lines, with miscibility gaps that display distinct
behavior. The left gap features the usual traveling bands (T.B.) as in Fig.
2(c). The right gap features reversed bands (R.B.), shown in Fig. 3(c). These
bands propagate in the direction opposite to the bulk magnetization. Indeed,
in the corresponding miscibility gap, the coexistence is between a magnetized
lower-density phase and a paramagnetic high-density phase. Then, according to
the flux balance condition (15), it follows that the phase boundaries
propagate in the reverse direction from usual.

Our scenario shows that the flocking transition need not necessarily be
discontinuous, at variance with a widely held view among those studying
flocking [2], although it remains an exceptional case. In order to rule out
some continuous transition scenarios, note that the disordered (H., $m=0$) and
the ordered phases (H., $m\neq 0$ at $\text{Pe}=0$ or C.M. at $\text{Pe}\neq
0$) entail different symmetries. Thus, it is not possible to pass between
these pure phases without crossing a (symmetry-breaking) phase transition of
some kind. This argument precludes the standard scenario of liquid-liquid
phase separation, where a single pair of binodals terminate at a critical
point. Indeed, one sees from the equilibrium phase diagrams in Figs. 2(a) and
3(a) that the binodals terminate at tricritical points, and any path between
the pure phases must cross a phase transition line. In the non-equilibrium
phase diagrams of Figs. 2(b) and 3(b), one observes discontinuous transitions
as well as an azeotropic point, all of which are consistent with the relevant
symmetry principles.

###  4.2 Triple phase in the nonequilibrium case ($\rm{Pe}>0$)

When $f(\rho)$ is sufficiently curved, two tricritical points in the
equilibrium limit are closer to the maximum of $f(\rho)$, see Fig. 4(a). Here,
we observe a different behavior at finite Pe compared with Sec. 4.1. Instead
of the azeotropic point reported in Fig. 3(b), there now appears a region
which we call triple phase coexistence, as shown in Fig. 4(b). This happens
when the two magnetised binodals (red solid lines) merge, at a finite
temperature $T_{\text{triple}}$. For temperatures between $T_{\text{triple}}$
and $T_{\text{max}}=\text{ max}\left[f(\rho)\right]$ , the system enters the
triple phase coexistence (see Fig. 4(c)), a strongly fluctuating dynamical
steady state whose properties can be understood in terms of the partial
restoration of a broken symmetry, as we discuss next.

Figure 4: Phase diagrams for (a) the equilibrium case ($\text{Pe}=0$) and (b)
the nonequilibrium regime ($\text{Pe}=1$), with non-monotonic $f$ given by
(19) with $\alpha=1/2$. The colored dashed and solid curves are spinodals and
binodals, respectively. The black solid curve is a critical line. The black
dots mark tricritical points. The black dashed line marks the maximal
magnetization curve (22). A triple phase coexistence emerges in the
nonequilibrium regime. (c) The triple phase profiles for $\text{Pe}=1$,
$\rho_{0}\simeq 1.8$, and $T\simeq 0.17$.

Recall that within the symmetry-broken T.B. and R.B. phases, the steady state
consists of a single magnetized band, which travels through the system. The
magnetization may be either positive or negative: this spontaneous breaking of
spin-reversal symmetry is familiar from the Ising model. In turn, this
symmetry breaking drives a particle current, which causes the band to move
either left or right. This additionally breaks the spatial reflection symmetry
of the model. Due to the periodic boundary conditions of the system, the
resulting steady state is also time-periodic, with period $\ell_{s}/V$.

Figure 5: Time series of the triple phase coexistence for $\text{Pe}=1$,
$\rho_{0}\simeq 1.8$, and $T\simeq 0.17$. (Top) The counter-propagating bands
of the $\rho_{\pm}$ fields, in cyan and magenta, barely scatter upon
encounter. (Bottom) The corresponding density $\rho=\rho_{+}+\rho_{-}$ and
magnetization $m=\rho_{+}-\rho_{-}$ display an oscillating dynamical state
with volume fraction periodically varying in time.

In the triple phase, both left- and right-moving bands are present at the same
time, as shown in Fig. 5. Since they propagate in opposite directions, the two
bands regularly meet each other. It turns out that they interact weakly (we
recall that the number of particles on each site is unbounded, in contrast
with exclusion processes). The resulting situation is most simply represented
in terms of the density of $+$ and $-$ particles

$\rho_{\pm}\left(x,t\right)=\frac{\rho\left(x,t\right)\pm
m\left(x,t\right)}{2}.$ (20)

As shown in Figs. 6(a-b), the space-time dependence of these profiles is well-
approximated by $\rho_{+}\simeq\bar{\rho}\left(x-Vt\right)$ and
$\rho_{-}\simeq\bar{\rho}\left(-x-Vt\right)$. In this case, the density
$\rho(x,t)$ is still time-periodic, with behaviour invariant under spatial
reflection, so that it is analogous to a standing-wave solution, while the
T.B. and R.B. are analogous to travelling-wave solutions. However, the band
profiles $\bar{\rho}$ are far from the sinusoidal functions that would be
found from the standard wave equations.

In fact, $\bar{\rho}$ is close to a top-hat profile, as it was in the
T.B./R.B. phases: there are extended bulk regions at densities
$(\rho_{1},\rho_{2})$, separated by interfaces whose widths are
$\mathcal{O}\left(1\right)\ll\ell_{s}$. Plotting the density profile at
generic times, as shown in Fig. 5, one typically finds up to four distinct
regions whose densities are the following. Outside of either band the density
is that of the dilute binodal $\rho_{g}=2\rho_{1}$. Inside both bands, the
density is that of the dense binodal $\rho_{l}=2\rho_{2}$. Finally, inside one
band and outside the other, the density is that of the intermediate binodal
$\rho_{i}=\rho_{1}+\rho_{2}$. Hence the intermediate binodal is always halfway
between the extreme binodals, while their difference equals twice the
magnetization $m_{0}$:

$\rho_{g}\left(T\right)=\rho_{i}\left(T\right)-m_{0}\left(T\right),\quad\rho_{l}\left(T\right)=\rho_{i}\left(T\right)+m_{0}\left(T\right).$
(21)

An important consequence of (21) is that the propagation velocity (15) attains
its limiting minimal value $|V|=\text{Pe}$. The magnetization in (21) marks
the bulk magnetization of the intermediate binodal
$m_{0}=m_{0}\left(\rho_{i},T\right)$, as given explicitly by (12). Therefore,
both binodals $\rho_{g,l}$ in (21) are fully determined by the value of the
intermediate binodal $\rho_{i}$. In practice, we find empirically that the
binodal $\rho_{i}$ coincides with the maximal magnetization curve
$\rho_{\text{max}}$, defined by $\partial_{\rho}m_{0}(\rho_{\text{max}},T)=0$,
as shown in dashed black line in Fig. 4(b). We get an implicit representation
of $\rho_{\text{max}}(T)$ by differentiating (12) with respect to $\rho$,
yielding

$\frac{x}{\left(1-x^{2}\right)\operatorname{arctanh}x}=1-\frac{\rho_{\text{max}}f^{\prime}\left(\rho_{\text{max}}\right)}{f\left(\rho_{\text{max}}\right)},\quad
T=\frac{x}{\operatorname{arctanh}x}f\left(\rho_{\text{max}}\right).$ (22)

Overall, Equations (21-22) specify all three binodals of the triple phase in
closed analytic form.

Figure 6:  (a-b) The $\pm$ traveling band profiles of the triple phase are
perfectly anti-symmetric: $\rho_{+}\simeq\bar{\rho}\left(x-Vt\right)$ and
$\rho_{-}\simeq\bar{\rho}\left(-x-Vt\right)$ coincide when reflected and
superimposed. (c) The periodic evolution of the phase volumes of the
intermediate and liquid phases $\nu_{i}$ and $\nu_{l}$ respectively, for
parameter values in Fig. 5. The complementary gas phase volume (not shown) is
simply $\nu_{g}=1-\nu_{i}-\nu_{l}$.

We now address the phase volumes of the bulk phases, defined as the fraction
of space that each phase occupies, and the corresponding highly dynamical
nature of the triple phase state. Since the profiles of $\rho_{+}$ and
$\rho_{-}$ are symmetric, each one of them carries exactly half of the total
mass $L\rho_{0}/2$. Besides, given that $\pm$ bands counter-propagate at
constant speed, their overlap region varies periodically with time. We report
in Fig. 6(c) the corresponding dynamics of the phase volumes for a mean
density below the intermediate binodal ($\rho_{0}<\rho_{i}$). The maximal
phase volume of the intermediate phase $\nu_{i}^{max}$ is set by the usual
lever rule between the intermediate and gas binodals:
$\nu_{i}^{max}=\left(\rho_{i}-\rho_{0}\right)/\left(\rho_{i}-\rho_{g}\right)$.
Then, the maximal phase volume for the liquid phase occurs at perfect overlap
of the $\pm$ profiles and is exactly $\nu_{i}^{max}/2$, as shown in Fig. 6(c).

Interestingly, one can formulate a necessary condition for the emergence of
the triple phase state. Indeed, the liquid and gas binodals must lie outside
the spinodals $\rho_{g}^{s}$ and $\rho_{l}^{s}$:

$\rho_{g}\left(T\right)<\rho_{g}^{s}\left(T\right)<\rho_{l}^{s}\left(T\right)<\rho_{l}\left(T\right).$
(23)

The condition (23) must hold across the triple phase. In particular, it should
hold in the vicinity of the maximum of $f$, where the triple phase first
emerges. Approximating $f(\rho)$ by a parabola in this region, the spinodal
densities $\rho_{g,l}^{s}$ defined in (11) are symmetrically positioned around
the maximum $\rho^{*}$:

$\rho_{g}^{s}\simeq\rho^{*}-\Delta_{s},\quad\rho_{l}^{s}\simeq\rho^{*}+\Delta_{s},\quad\Delta_{s}=\sqrt{\frac{2\epsilon
f(\rho^{*})}{\left|f^{\prime\prime}\left(\rho^{*}\right)\right|}},\quad\epsilon=\frac{f(\rho^{*})-T}{f(\rho^{*})}\ll
1.$ (24)

Performing a similar expansion for the maximal magnetization curve in (22)
yields

$\rho_{\text{max}}\left(T\right)\simeq\rho^{*},\quad
m_{0}\left(\rho_{\text{max}},T\right)\simeq\rho^{*}\sqrt{3\epsilon}.$ (25)

Substituting the expression (25) in (21), and recalling that
$\rho_{i}=\rho_{\text{max}}$ in the triple phase, we find that the liquid and
gas binodals are also symmetrically positioned around the maximum:

$\rho_{g}\simeq\rho^{*}-\Delta_{b},\quad\rho_{l}\simeq\rho^{*}+\Delta_{b},\quad\Delta_{b}=\sqrt{3\epsilon}\rho^{*}.$
(26)

Finally, the condition (23) for existence of the triple phase is equivalent to
the ordering $\Delta_{b}>\Delta_{s}$ of the binodal and spinodal gaps. Using
the expressions for these gaps, respectively in (24) and (26), we arrive at
the following condition for the curvature near the maximum:

$\frac{\rho^{*2}\left|f^{\prime\prime}\left(\rho^{*}\right)\right|}{f\left(\rho^{*}\right)}>\frac{2}{3}.$
(27)

Remarkably, the condition (27) is independent of the Péclet number, which
corroborates that the triple phase appears as a singular perturbation of the
underlying equilibrium phase diagram at any finite activity, no matter how
small.

##  5 Discussion

Our results show that the density-dependence of the aligning interactions can
play an unexpectedly major role in flocking models, yielding a richer flocking
phenomenology than any previously reported. These results are found by
choosing a thermodynamically consistent approach which entails a proper
equilibrium limit, pertaining to the universality class of Model C [33, 16],
at vanishing self-propulsion. Importantly, the impact of density-dependent
interactions on the equilibrium phase diagram allows one to rationalize the
emergence of exotic phase topologies in the presence of self-propulsion.
Indeed, for the nonequilibrium regime, we not only recover the standard
scenario of discontinuous flocking [2], we also reveal novel phase behaviors
where the system displays either an azeotropic point or a region of triple-
phase coexistence. Remarkably, in contrast with standard flocking models [2],
the flocking transition is continuous in the presence of an azeotropic point,
while remaining discontinuous elsewhere on the phase diagram.

Given the rich topology of these phase diagrams, it would be interesting to
examine how the entropy production rate (EPR) varies across the various phase
transitions that we have unveiled. In thermal systems [23, 19], EPR measures
both the breakdown of time-reversal symmetry and the amount of energy
dissipated by nonequilibrium currents. For flocking systems, irreversibility
measures have been provided in both microscopic models [39, 40, 41, 42] and
field theories [43, 44]. Yet, these measures should coincide with a
quantification of dissipation only for thermodynamically consistent models
[45, 46, 47]. In that respect, our approach opens the unprecedented
opportunity to properly study how the dissipation varies across the flocking
transition. For instance, it could lead to a spatial resolution of
dissipation, thus uncovering the profile of energy dissipation associated with
travelling bands.

## Acknowledgements

The authors acknowledge useful discussions with M. Esposito, G. Falasco, K.
Proesmans, and A. Tanaji Mohite. This research was funded in part by the
Luxembourg National Research Fund (FNR), grant reference 14389168. For the
purpose of open access, ÉF has applied a Creative Commons Attribution 4.0
International (CC BY 4.0) license to any Author Accepted Manuscript version
arising from this submission.

##  Appendix A Deriving hydrodynamic equations

The microscopic dynamics described in Sec. 2 consists of biased diffusive hops
along with reaction dynamics. To derive the corresponding hydrodynamics, we
follow the recipe outlined in [24], which has already been implemented in
other models involving active lattice gases [15, 30, 32, 31]. A necessary
property for the derivation is that the underlying microscopic dynamics is
dominated by symmetric diffusive hops, which have the fastest transition rate.

First, we express the Hamiltonian (2) in terms of the hydrodynamic fields (3).
Assuming that these fields are slowly varying over the macroscopic scale
(which will be verified a posteriori), then one has, to leading order at large
$L$,

$H=\frac{L}{\ell_{s}}\int_{0}^{\ell_{s}}dx\mathcal{H}\left[\rho\left(x,t\right),m\left(x,t\right)\right],$
(28)

with the Hamiltonian density

$\mathcal{H}\left(\rho,m\right)=-\frac{m^{2}}{2\rho}f\left(\rho\right).$ (29)

When a $+$ particle at position $i=Lx_{i}/\ell_{s}$ hops to the right, the
field variations follow as

$\rho\to\rho+\frac{\ell_{s}}{L}\delta\left(x-x_{i+1}\right)-\frac{\ell_{s}}{L}\delta\left(x-x_{i}\right)=\rho-\frac{\ell_{s}^{2}}{L^{2}}\partial_{x}\delta\left(x-x_{i}\right)+o\left(L^{-2}\right).$
(30)

Similarly, we have

$m\to
m-\frac{\ell_{s}^{2}}{L^{2}}\partial_{x}\delta\left(x-x_{i}\right)+o\left(L^{-2}\right),$
(31)

and analogous expressions hold for the other jumps. The expression for the
field variations under a flip of $+\to-$ is simply given by

$\rho\to\rho,\quad m\to
m-2\frac{\ell_{s}}{L}\delta\left(x-x_{i}\right)+o\left(L^{-1}\right).$ (32)

Substituting the expressions (30-31) into (28), we find that the jump of a $+$
particle to the right leads to the Hamiltonian increment

$\Delta H=\frac{\ell_{s}}{L}\partial_{x}\left(\frac{\delta
H}{\delta\rho}+\frac{\delta H}{\delta m}\right),$ (33)

with similar expressions hold for the other jumps. The Hamiltonian increment
under $+\to-$ flipping is simply

$\Delta H=-2\frac{\delta H}{\delta m}.$ (34)

Substituting (33-34) in the definition of the rates of Sec. 2, we find that
these rates are functions of the local fields to leading order:

$\displaystyle\text{rate}_{+}(i\rightarrow i+1)$
$\displaystyle=D_{0}\left[1+\frac{1}{L}\left(\frac{\lambda}{D_{0}}-\frac{\beta}{2}\ell_{s}\partial_{x}\left.\left(\frac{\delta
H}{\delta\rho}+\frac{\delta H}{\delta
m}\right)\right|_{\rho(x,t),m(x,t)}\right)\right],$ (35)
$\displaystyle\text{rate}_{+}(i\rightarrow i-1)$
$\displaystyle=D_{0}\left[1+\frac{\beta\ell_{s}}{2L}\left.\left(\frac{\delta
H}{\delta\rho}+\frac{\delta H}{\delta
m}\right)\right|_{\rho(x,t),m(x,t)}\right],$
$\displaystyle\text{rate}_{-}(i\rightarrow i+1)$
$\displaystyle=D_{0}\left[1-\frac{\beta\ell_{s}}{2L}\left.\left(\frac{\delta
H}{\delta\rho}-\frac{\delta H}{\delta
m}\right)\right|_{\rho(x,t),m(x,t)}\right],$
$\displaystyle\text{rate}_{-}(i\rightarrow i-1)$
$\displaystyle=D_{0}\left[1+\frac{1}{L}\left(\frac{\lambda}{D_{0}}+\frac{\beta}{2}\ell_{s}\partial_{x}\left.\left(\frac{\delta
H}{\delta\rho}-\frac{\delta H}{\delta
m}\right)\right|_{\rho(x,t),m(x,t)}\right)\right],$
$\displaystyle\text{rate}(+\rightarrow-)$
$\displaystyle=\frac{\gamma}{L^{2}}e^{\beta\frac{\delta H}{\delta m}}.$

Importantly, the leading order terms are symmetric hops with rate $D_{0}$.
These dominate over both the asymmetric $\mathcal{O}\left(L^{-1}\right)$ jump
rate, and over the slower $\mathcal{O}\left(L^{-2}\right)$ tumbling rate
defined in Sec. 2. This timescale separation guarantees that the measure is
controlled by the simple symmetric hop rules corresponding to noninteracting
random walkers [24]. It is then straightforward to write down the coupled
hydrodynamics of the $\rho_{+}$ and $\rho_{-}$ fields

$\begin{bmatrix}\partial_{t}{\rho_{+}}\\\
\partial_{t}{\rho_{-}}\end{bmatrix}=-\partial_{x}\begin{bmatrix}J_{\rho_{+}}\\\
J_{\rho_{-}}\end{bmatrix}+\begin{bmatrix}-K\\\ K\end{bmatrix}.$ (36)

The hydrodynamic fluxes are found by coarse-graining, using the hydrodynamic
coordinates $x=\ell_{s}i/L$ and $t=\gamma\hat{t}/L^{2}$, with $\hat{t}$ the
time variable of the microscopic dynamics:

$\displaystyle J_{\rho_{+}}$
$\displaystyle=-\partial_{x}\rho_{+}-\beta\rho_{+}\partial_{x}\left(\frac{\delta
H}{\delta\rho}+\frac{\delta H}{\delta m}\right)+\text{Pe}\rho_{+},$ (37)
$\displaystyle J_{\rho_{-}}$
$\displaystyle=-\partial_{x}\rho_{-}-\beta\rho_{-}\partial_{x}\left(\frac{\delta
H}{\delta\rho}-\frac{\delta H}{\delta m}\right)-\text{Pe}\rho_{-},$

and, following the same steps, the tumbling rate reads

$K=me^{\beta\frac{\delta H}{\delta m}}.$ (38)

Note that, when sending $\beta\to 0$ in (37), one recovers the fluxes of AIM
[13]. Finally, one obtains the expressions (5-6) in terms of $\rho$ and $m$
with a few steps of algebra. Interestingly, in the equilibrium limit
($\text{Pe}=0$), the free energy $F$ serves as a Lyapunov function:

$\frac{dF}{dt}=-\frac{L}{\ell_{s}}\int_{0}^{\ell_{s}}dx\left\\{\frac{1}{2}\partial_{x}\begin{bmatrix}\frac{\delta
F}{\delta\rho}\\\ \frac{\delta F}{\delta
m}\end{bmatrix}^{\dagger}\mathbb{C}(\rho,m)\partial_{x}\begin{bmatrix}\frac{\delta
F}{\delta\rho}\\\ \frac{\delta F}{\delta m}\end{bmatrix}+2M\frac{\delta
F}{\delta m}\sinh\left(\frac{\delta F}{\delta m}\right)\right\\}\leq 0,$ (39)

where we have used that $M>0$ and $\mathbb{C}$ is positive semi-definite. The
case $d{F}/dt=0$ is achieved in steady state. In the presence of noise, this
case corresponds to the system eventually relaxing towards the minimum $\delta
F/\delta\rho=\delta F/\delta m=0$. Notice, however, that the well-posedness of
this hydrodynamic description must break in the equilibrium limit
$\text{Pe}=0$ for any non-convex free energy $\mathcal{F}$ (6). It is cured by
retaining sub-leading interfacial terms
$\mathcal{O}\left(L^{-2}|\nabla\rho|^{2},L^{-2}|\nabla m|^{2}\right)$.

##  Appendix B Spinodals and binodals in the equilibrium limit

As mentioned in Sec. 3.1, the density spinodal (13), which is complementary to
the magnetization one (11), marks the boundary of the convex region of the
single-variable function
$\mathcal{F}\left[\rho,m_{0}\left(\rho,T\right)\right]$ [36]. Indeed, this
boundary is given by

$\left.\frac{\partial^{2}\mathcal{F}}{\partial\rho^{2}}\right|_{\rho_{0},m_{0}\left(\rho,T\right)}+\left.\frac{\partial^{2}\mathcal{F}}{\partial\rho\partial
m}\right|_{\rho_{0},m_{0}\left(\rho,T\right)}\left.\frac{\partial
m_{0}\left(\rho,T\right)}{\partial\rho}\right|_{\rho_{0}}=0,$ (40)

where we used the fact that, in the ordered state
$m_{0}\left(\rho_{0},T\right)$, we have

$\left.\frac{\partial\mathcal{F}}{\partial
m}\right|_{\rho_{0},m_{0}\left(\rho,T\right)}=0.$ (41)

Differentiating (41) with respect to $\rho$ yields

$\left.\frac{\partial
m_{0}\left(\rho,T\right)}{\partial\rho}\right|_{\rho_{0}}=-\left.\left(\frac{\partial^{2}\mathcal{F}}{\partial\rho\partial
m}\bigg{/}\frac{\partial^{2}\mathcal{F}}{\partial^{2}m}\right)\right|_{\rho_{0},m_{0}\left(\rho_{0},T\right)}.$
(42)

Plugging the relation (42) into (40), it follows that finding the density
spinodal is equivalent to locating the onset of local non-convexity of the
two-variable free energy:

$\left|\text{Hess}\left(\mathcal{F}\right)\right|_{\rho_{0},m_{0}\left(\rho_{0},T\right)}=0.$
(43)

The density spinodal can be expressed in an implicit parametric form, since
all partial derivatives entering (43) are explicit functions of $\rho$ and
$m$. Indeed, we can use the parametric representation (12)

$x=\frac{m_{0}}{\rho_{0}},\quad
T\left(x,\rho_{0}\right)=\frac{x}{\text{arctanh}(x)}f\left(\rho_{0}\right),$
(44)

to explicitly write (43) as

$\left(\frac{1}{1-x^{2}}-\frac{\text{arctanh}(x)}{x}\right)\left[1-\frac{x\text{arctanh}(x)}{2}\frac{(d/d\rho)\left(\rho^{2}f^{\prime}\left(\rho\right)\right)}{f\left(\rho\right)}\right]=\text{arctanh}^{2}(x)\left(\frac{\rho
f^{\prime}\left(\rho\right)}{f\left(\rho\right)}\right)^{2},$ (45)

which gives $x=x\left(\rho_{0}\right)$ in an implicit form, from which follows
the curve $T\left(\rho_{0}\right)$ for the density spinodal. This procedure
allows us to determine the corresponding curves in Figs. 2(a), 3(a) and 4(a).

The binodals are found by the common-tangent construction over the single-
variable free energy
$\mathcal{F}\left[\rho,m_{0}\left(\rho_{0},T\right)\right]$, which can result
in several possible phase-diagram topologies [36]. The binodals correspond to
the coexistence between dense and dilute phases. For monotonic
$f\left(\rho\right)$ in Sec. 3, it is always the dense and the dilute phases
which are respectively ordered ($m=m_{0}(\rho,T)$) and disordered ($m=0$).
Yet, one can have the opposite case for non-monotonic $f(\rho)$, as explained
in Sec. 4. For illustration purposes, considering here the case of monotonic
$f(\rho)$, the condition for common tangent is given by

$\displaystyle\mathcal{F}\left[\rho_{l},m_{0}\left(\rho_{l},T\right)\right]$
$\displaystyle=\mathcal{F}\left(\rho_{g},m=0\right)+\partial_{\rho}\mathcal{F}\left(\rho_{g},m=0\right)\left(\rho_{l}-\rho_{g}\right),$
(46)
$\displaystyle\partial_{\rho}\mathcal{F}\left[\rho_{l},m_{0}\left(\rho_{l},T\right)\right]$
$\displaystyle=\partial_{\rho}\mathcal{F}\left(\rho_{g},m=0\right),$

which explicitly reads

$\displaystyle\frac{\rho_{l}}{2}\left[\ln\frac{\rho_{l}^{2}\left(1-x^{2}\right)}{4}+x\ln\frac{1+x}{1-x}\right]-\rho_{l}-\beta\frac{\rho_{l}x^{2}}{2}f\left(\rho_{l}\right)$
$\displaystyle=\rho_{g}\left(\log\frac{\rho_{g}}{2}-1\right)+\log\frac{\rho_{g}}{2}\left(\rho_{l}-\rho_{g}\right),$
(47)
$\displaystyle\log\left(\frac{\rho_{l}}{2}\sqrt{1-x^{2}}\right)-\frac{\beta
x^{2}}{2}\left[\rho_{l}f^{\prime}\left(\rho_{l}\right)-f\left(\rho_{l}\right)\right]$
$\displaystyle=\log\frac{\rho_{g}}{2},$

where we have used the parametric representation
$x=m_{0}(\rho_{l},T)/\rho_{l}$. Substituting the temperature relation (44)
into (47), we get

$1-\frac{\rho_{g}}{\rho_{l}}=\frac{x\text{arctanh}x}{2}\frac{\rho_{l}f^{\prime}\left(\rho_{l}\right)}{f\left(\rho_{l}\right)},\qquad\frac{\rho_{g}}{\rho_{l}}=\sqrt{1-x^{2}}e^{-\frac{x\text{arctanh}x}{2}\left[\frac{\rho_{l}f^{\prime}\left(\rho_{l}\right)}{f\left(\rho_{l}\right)}-1\right]}$
(48)

which finally gives the two binodals in an implicit parametric representation
for $\rho_{l}\left(x\right)$, $\rho_{g}\left(x\right)$ and
$T=T\left(x\right)$. The tricritical point is defined as the meeting point of
the two binodals (48), which is also the meeting point of the two spinodals
(45) and (11). This occurs at vanishing $x=m_{0}/\rho_{0}$. Indeed, expanding
either one of these expressions at small $x$ leads to (14).

##  Appendix C Linear instability in the nonequilibrium case

To derive the expression of the density spinodal in (16), we examine the
linear stability of the dynamics (4-5) around the homogeneous solutions
$\rho\left(x,t\right)=\rho_{0}+\delta\rho\left(x,t\right)$ and
$m\left(x,t\right)=m_{0}\left(\rho_{0},T\right)+\delta m\left(x,t\right)$,
where $m_{0}$ is given by (12). Using the Fourier convention for an arbitrary
field $X$:

$\delta X(x,t)=\sum_{n}\;\delta X_{n}(t)e^{ik_{n}x},\quad k_{n}=2\pi
n/\ell_{s},$ (49)

we arrive at the following mode dynamics

$\frac{d}{dt}\begin{bmatrix}{\delta\rho_{n}}\\\ {\delta
m_{n}}\end{bmatrix}=\mathbb{M}\begin{bmatrix}\delta\rho_{n}\\\ \delta
m_{n}\end{bmatrix},$ (50)

where the stability matrix $\mathbb{M}$ reads

$\mathbb{M}=\left.\begin{bmatrix}-k^{2}A&-k^{2}B-\text{Pe}ik\\\
-k^{2}C-\text{Pe}ik-2\alpha_{2}&-k^{2}D-2\alpha_{3}\end{bmatrix}\right|_{\rho_{0},m_{0}(\rho_{0},T)}$
(51)

with

$\displaystyle A=\rho\partial_{\rho\rho}\mathcal{F}+m\partial_{\rho
m}\mathcal{F},\quad B=\rho\partial_{\rho
m}\mathcal{F}+m\partial_{mm}\mathcal{F},\quad\alpha_{2}=M\partial_{\rho
m}\mathcal{F},$ (52) $\displaystyle
C=m\partial_{\rho\rho}\mathcal{F}+\rho\partial_{\rho m}\mathcal{F},\quad
D=m\partial_{\rho
m}\mathcal{F}+\rho\partial_{mm}\mathcal{F},\quad\alpha_{3}=M\partial_{mm}\mathcal{F},$

and $M\left(\rho,m\right)=\sqrt{\rho^{2}+m^{2}}$. The eigenvalues of (51) are
of the form

$\Lambda_{1,2}=\frac{1}{2}\left[\text{Tr}\pm\sqrt{\text{Tr}^{2}-4\text{Det}}\right],$
(53)

where

$\text{Tr}=-\left[2\alpha_{3}+k^{2}\left(A+D\right)\right],\quad\text{Tr}^{2}-4\text{Det}=a^{2}+ikb+k^{2}c+\mathcal{O}(k^{3}),$
(54)

with

$a=2\alpha_{3},\quad b=8\text{Pe}\,\alpha_{2},\quad
c=4\left[\alpha_{3}\left(D-A\right)+2B\alpha_{2}-\text{Pe}^{2}\right],$ (55)

yielding

$\text{Re}\left[\sqrt{\text{Tr}^{2}-4\text{Det}}\right]=a+\frac{k^{2}}{a}\left(\frac{b^{2}}{8a^{2}}+\frac{c}{2}\right)+\mathcal{O}(k^{4}),\quad\text{Im}\left[\sqrt{\text{Tr}^{2}-4\text{Det}}\right]=k\frac{b}{2a}+\mathcal{O}(k^{3}),$
(56)

where we used that $\alpha_{3}\propto\partial_{mm}\mathcal{F}>0$ at the
magnetized state $\left[\rho_{0},m_{0}\left(\rho_{0},T\right)\right]$. (Note
that the expressions $a,b,c$ in (55) are not to be confused with the Model C
coefficients (9).) Then, the eigenvalue with the larger real part is given by

$\Lambda_{+}=ik\text{Pe}\frac{\partial_{\rho
m}\mathcal{F}}{\partial_{mm}\mathcal{F}}+k^{2}\left\\{\frac{\text{Pe}^{2}}{2M\mathcal{F}_{mm}}\left[\left(\frac{\partial_{\rho
m}\mathcal{F}}{\partial_{mm}\mathcal{F}}\right)^{2}-1\right]+\frac{\rho_{0}\text{Hess}\left(\mathcal{F}\right)}{\mathcal{F}_{mm}}\right\\}.$
(57)

The condition of vanishing real part for the eigenvalue (57) yields the
ordered spinodal (16).

## References

## References

  * [1] Toner J, Tu Y and Ramaswamy S 2005 Ann. Phys. 318 170–244
  * [2] Chaté H 2020 Annu. Rev. Condens. Matter Phys. 11 189–212
  * [3] Cavagna A and Giardina I 2014 Annu. Rev. Condens. Matter Phys. 5 183–207
  * [4] Bain N and Bartolo D 2019 Science 363 46–49
  * [5] Bricard A, Caussin J B, Desreumaux N, Dauchot O and Bartolo D 2013 Nature 503 95–98
  * [6] Vicsek T, Czirók A, Ben-Jacob E, Cohen I and Shochet O 1995 Phys. Rev. Lett. 75 1226–1229
  * [7] Toner J and Tu Y 1995 Phys. Rev. Lett. 75 4326–4329
  * [8] Marchetti M C, Joanny J F, Ramaswamy S, Liverpool T B, Prost J, Rao M and Simha R A 2013 Rev. Mod. Phys. 85 1143–1189
  * [9] Mermin N D and Wagner H 1966 Phys. Rev. Lett. 17 1133
  * [10] Halperin B 2018 J. Stat. Phys. 175 521–529
  * [11] Grégoire G and Chaté H 2004 Phys. Rev. Lett. 92 025702
  * [12] Solon A P and Tailleur J 2015 Phys. Rev. E 92 042119
  * [13] Solon A P and Tailleur J 2013 Phys. Rev. Lett. 111 078101
  * [14] Solon A P, Caussin J B, Bartolo D, Chaté H and Tailleur J 2015 Phys. Rev. E 92 062111
  * [15] Kourbane-Houssene M, Erignoux C, Bodineau T and Tailleur J 2018 Phys. Rev. Lett. 120 268003
  * [16] Scandolo M, Pausch J and Cates M E 2023 Eur. Phys. J. E 46 103
  * [17] Fodor E, Nardini C, Cates M E, Tailleur J, Visco P and van Wijland F 2016 Phys. Rev. Lett. 117 038103
  * [18] Fodor E, Jack R L and Cates M E 2022 Annu. Rev. Condens. Matter Phys. 13 215–238
  * [19] Seifert U 2012 Rep. Prog. Phys. 75 126001
  * [20] Ballerini M, Cabibbo N, Candelier R, Cavagna A, Cisbani E, Giardina I, Lecomte V, Orlandi A, Parisi G, Procaccini A, Viale M and Zdravkovic V 2008 Proc. Natl. Acad. Sci. USA 105 1232–1237
  * [21] Ginelli F and Chaté H 2010 Phys. Rev. Lett. 105 168103
  * [22] Martin D, Chaté H, Nardini C, Solon A, Tailleur J and Van Wijland F 2021 Phys. Rev. Lett. 126 148001
  * [23] Lebowitz J L and Spohn H 1999 J. Stat. Phys. 95 333–365
  * [24] Bodineau T and Lagouge M 2010 J. Stat. Phys. 139 201–218
  * [25] Gorbonos D, Puckett J G, Van Der Vaart K, Sinhuber M, Ouellette N T and Gov N S 2020 J. R. Soc. Interface 17 20200367
  * [26] Bastien R and Romanczuk P 2020 Sci. Adv. 6 eaay0792
  * [27] King A E B T and Turner M S 2021 R. Soc. Open Sci. 8 201536
  * [28] Wirth T D, Dachner G C, Rio K W and Warren W H 2023 PNAS Nexus 2 pgad118
  * [29] Jona-Lasinio G, Landim C and Vares M E 1993 Prob. Theory Relat. Fields 97 339–361
  * [30] Agranov T, Ro S, Kafri Y and Lecomte V 2021 J. Stat. Mech. 2021 083208
  * [31] Agranov T, Ro S, Kafri Y and Lecomte V 2023 SciPost Phys. 14 045
  * [32] Agranov T, Cates M E and Jack R L 2022 J. Stat. Mech. 2022 123201
  * [33] Bray A J 2002 Adv. Phys. 51 481–587
  * [34] Kaiser M, Jack R L and Zimmer J 2018 J. Stat. Phys. 170 1019–1050
  * [35] Blume M, Emery V J and Griffiths R B 1971 Phys. Rev. A 4 1071–1077
  * [36] Roux D, Coulon C and Cates M 1992 J. Phys. Chem. 96 4174–4187
  * [37] Pearce D J G, Miller A M, Rowlands G and Turner M S 2014 Proc. Natl. Acad. Sci. USA 111 10422–10426
  * [38] Koynova R D, Boyanov A I and Tenchov B G 1987 Biochim. Biophys. Acta - Biomembr. 903 186–196
  * [39] Shim P S, Chun H M and Noh J D 2016 Phys. Rev. E 93 012113
  * [40] Crosato E, Prokopenko M and Spinney R E 2019 Phys. Rev. E 100 042613
  * [41] Ferretti F, Grosse-Holz S, Holmes C, Shivers J L, Giardina I, Mora T and Walczak A M 2022 Phys. Rev. E 106 034608
  * [42] Yu Q and Tu Y 2022 Phys. Rev. Lett. 129 278001
  * [43] Dadhichi L P, Maitra A and Ramaswamy S 2018 J. Stat. Mech. 2018 123201
  * [44] Borthne Ø L, Fodor É and Cates M E 2020 New J. Phys. 22 123012
  * [45] Markovich T, Fodor E, Tjhung E and Cates M E 2021 Phys. Rev. X 11 021057
  * [46] Falasco G and Esposito M 2023 arXiv e-prints arXiv:2307.12406
  * [47] K Proesmans et al., in preparation