Generalized thermodynamics of phase
equilibria in scalar active matter
The MIT Faculty has made this article openly available. Please share
how this access benefits you. Your story matters.
Citation Solon, Alexandre P. et al. "Generalized thermodynamics of phase
equilibria in scalar active matter." Physical Review E 97, 2 (February
2018): 020602(R) © 2018 American Physical Society
As Published http://dx.doi.org/10.1103/PhysRevE.97.020602
Publisher American Physical Society
Version Final published version
Citable link http://hdl.handle.net/1721.1/114384
Terms of Use Article is made available in accordance with the publisher's
policy and may be subject to US copyright law. Please refer to the
publisher's site for terms of use.
PHYSICAL REVIEW E 97, 020602(R) (2018)
Rapid Communications
Generalized thermodynamics of phase equilibria in scalar active matter
Alexandre P. Solon,
1
Joakim Stenhammar,
2
Michael E. Cates,
3
Yariv Kafri,
4
and Julien Tailleur
5
1
Massachusetts Institute of Technology, Department of Physics, Cambridge, Massachusetts 02139, USA
2
Division of Physical Chemistry, Lund University, 221 00 Lund, Sweden
3
DAMTP, Centre for Mathematical Sciences, University of Cambridge, Cambridge CB3 0WA, United Kingdom
4
Department of Physics, Technion, Haifa 32000, Israel
5
Université Paris Diderot, Sorbonne Paris Cité, MSC, UMR 7057 CNRS, 75205 Paris, France
(Received 15 September 2016; revised manuscript received 21 December 2017; published 20 February 2018)
Motility-induced phase separation (MIPS) arises generically in fluids of self-propelled particles when
interactions lead to a kinetic slowdown at high densities. Starting from a continuum description of scalar active
matter akin to a generalized Cahn-Hilliard equation, we give a general prescription for the mean densities of
coexisting phases in flux-free steady states that amounts, at a hydrodynamics scale, to extremizing an effective
free energy. We illustrate our approach on two well-known models: self-propelled particles interacting either
through a density-dependent propulsion speed or via direct pairwise forces. Our theory accounts quantitatively
for their phase diagrams, providing a unified description of MIPS.
DOI: 10.1103/PhysRevE.97.020602
Active materials, composed of particles individually ca-
pable of dissipatively converting energy into motion [1–5],
display a fascinating range of large-scale properties [6–12].
Among them, motility-induced phase separation [13] (MIPS)
has recently attracted a lot of interest [5,13–30]. It arises
because self-propelled particles accumulate in regions where
they move more slowly [31]. When interactions between
particles lead to their slowing down at high density, a pos-
itive feedback leads to phase separation between a high-
density low-motility phase and a low-density high-motility
phase. Remarkably, this liquid-gas phase separation hap-
pens without the need of any attractive interactions, lead-
ing to the emergence of cohesive matter without cohesive
forces. First postulated in idealized toy models [14–20],
MIPS has since been addressed experimentally using self-
propelled colloids [5,21] and genetically engineered bacteria
[32].
The aforementioned instability mechanism leading to MIPS
is by now well understood and has been used to define
a s pinodal region where homogeneous phases are linearly
unstable [13,14]. Furthermore, this instability can be under-
stood at the ( fluctuating) hydrodynamic level [14,18,25,33,34]
where the dynamics of active particles undergoing a kinetic
slowdown at high density reduce to an equilibrium model B
[35]. On the contrary, there is no comprehensive theory
predicting the binodals: the mapping to equilibrium breaks
down at higher order in gradients [23] and the corresponding
equilibrium predictions for the coexisting binodal densities are
violated [23,36].
MIPS has been observed in two broad classes of systems.
In a first class of models [14,15,22,33], MIPS arises from
an explicit density dependence of the propulsion speed v(ρ).
This mimics the way cells adapt their motion to the local
density measured through the concentration of a chemical
signal, and we refer to such particles as “quorum-sensing active
particles” (QSAPs). There, one can define a chemical potential
μ [14] which is equal in coexisting phases, but the coexisting
pressures, whether mechanical [36] or thermodynamic [23],
are unequal. In a second class of models [16–19], particles
propelled by a constant force interact via an isotropic, repulsive
pair potential; the slowdown triggering MIPS is now due to
collisions. Contrary to QSAPs, the mechanical pressure P
of such “pairwise-force active particles” (PFAPs), defined as
the force density on a confining wall, is equal in coexisting
phases. However, an effective chemical potential defined from
the thermodynamic equilibrium relation [37]
PV = Nμ − F
with ∂F/∂N = μ takes unequal values in coexisting phases,
causing violation of the equilibrium Maxwell equal-area con-
struction [29]. For both models, we thus lack a constraint
to complement the equality of pressure (PFAPs) or chemical
potential (QSAPs) to fix the values of coexisting densities.
The difference between these two classes of models can be
shown to stem from whether or not an effective momentum
conservation holds in the steady state [38]. When such a
conservation law is present, as in PFAPs, the pressure is given
by an underlying equation of state and is equal in the two
phases [24,29,39]. On the contrary, in the absence of this
conservation law, this is generically not the case. All in all,
a comprehensive theory of the phase equilibria in MIPS, that
would in particular encompass these two different classes of
models, remains elusive.
In this Rapid Communication, we propose a unified theory
of MIPS based on phenomenological hydrodynamic equations
of motion for the scalar density field. We show how the
binodals are determined at this level from a common tangent
construction on an effective free energy density. Our formalism
encompasses equilibrium systems for which one recovers the
standard thermodynamic free energy and, in that case only, the
equality among phases of both pressure and chemical potential.
We then show how this generic formalism can be applied to
precise models of QSAPs and PFAPs, accounting for their
phase diagrams. In particular, we show that different intensive
2470-0045/2018/97(2)/020602(6) 020602-1 ©2018 American Physical Society
SOLON, STENHAMMAR, CATES, KAFRI, AND TAILLEUR PHYSICAL REVIEW E 97, 020602(R) (2018)
quantities are equal between coexisting phases in PFAPs and
QSAPs.
General framework. We consider a continuum description
of active particles with isotropic, nonaligning interactions. In
this scalar active matter, the sole hydrodynamic field is thus the
conserved density ρ(r,t ), obeying ˙ρ =−∇ · J. By symmetry,
the current J vanishes in homogeneous phases. Its expansion
in gradients of the density involves only odd terms under space
reversal. At third order, we use [40]
˙ρ = ∇ · (M∇g); g = g
0
(ρ) + λ(ρ)|∇ρ|
2
− κ(ρ)ρ . (1)
The noiseless hydrodynamic equation (1) describes the evo-
lution of the average coarse-grained density field on scales
much larger than the correlation length and time. It can thus
be used to characterize fully phase-separated profiles, away
from the critical point where noise is irrelevant [41], to predict
binodal densities. Equation (1) plays the same role as the
Cahn-Hilliard equation does for equilibrium phase-separating
systems [41] but in general does not admit an equilibrium free
energy structure. In what follows, we first start with Eq. (1) and
show how to compute analytically its phase diagram. We then
consider a microscopic model of QSAPs for which we obtain
the coefficients of Eq. (1) in terms of microscopic parameters
by coarse-graining. Finally, we show that our formalism can
also be applied to PFAPs, even though closed expressions of
the coefficients appearing in (1) are not known explicitly in
this case.
Equation (1) predicts a linear instability of a homogeneous
profile of density ρ
0
whenever g
0
(ρ
0
) < 0[42]: this is the
standard linear instability leading to MIPS [14] and defines
the spinodal region. We now proceed to establish t he cor-
responding binodals. As in equilibrium, we consider a fully
phase-separated system. A macroscopic droplet of the minority
phase has an infinite curvature radius, and hence effectively
flat interfaces, so that curvature effects are negligible. As in
equilibrium, the problem, although n-dimensional, reduces
to studying the one-dimensional profile perpendicular to the
interface [41]. We thus consider a flat interface, parallel to ˆy,
between coexisting gas and liquid phases at densities ρ
g
and
ρ
. In a steady state with vanishing current, M∇g = 0 , so that
g is constant throughout the system: g(ρ(r,t)) =
¯
g. This yields
a first equation relating ρ
g
and ρ
:
g
0
(ρ
g
) = g
0
(ρ
) =
¯
g. (2)
A second relation can now be obtained by considering
a function R(ρ) and integrating g(ρ)∂
x
R across the inter-
face. Replacing g(ρ) by its value
¯
g or its explicit expres-
sion in Eq. (1), one finds two equivalent expressions for
x
x
g
g(ρ)∂
x
Rdx:
(R
− R
g
)
¯
g = φ(R
) − φ(R
g
)
+
x
x
g
λ(∂
x
ρ)
2
− κ∂
2
x
ρ
∂
x
Rdx, (3)
where x
g
and x
lie within the bulk gas and liquid phases,
R
/g
≡ R(ρ
/g
), and φ is defined by dφ/dR = g
0
(ρ). To
simplify Eq. (3), we choose R(ρ) such that
κR
=−(2λ + κ
)R
, (4)
where (
) denotes d/dρ. Then, one has that
λ(∂
x
ρ)
2
− κ∂
2
x
ρ
∂
x
R =−∂
x
κR
2
(∂
x
ρ)
2
, (5)
the integral of which vanishes between any two bulk planes
where ∂
x
ρ = 0. Equation (3) then yields a second constraint:
h
0
(R
) = h
0
(R
g
); h
0
(R) ≡ Rφ
(R) − φ(R). (6)
Because R is nonlinear in ρ, the lever rule,
ρ
V
+ ρ
g
V
g
= ρ
0
V is nonlinear in R, but still determines
the phase volumes V
,g
. Also the densities ρ
,g
do not vary as
one moves along the “tie line” by changing the global mean
density ρ
0
. This is not true generally in nonequilibrium phase
separation [43].
Equations (2) and (6) show the coexisting densities to satisfy
a common tangent construction on an effective (bulk) free
energy φ(R) =
g
0
(ρ)dR. The mathematical similarity with
an equilibrium common tangent construction can be traced to
the fact that Eq. (1) can be written as
˙ρ = ∇ · [M[ρ]∇g]; g =
δF
δR
, (7)
with F =
dr[φ(R) +
κ
2R
(∇R)
2
]. The stationary solutions of
Eq. (1) then correspond to extrema of the “effective” free
energy F. Note that (7) holds in any dimension. This highlights
that, although the construction of the binodals (2)–(6) relies
on a single coordinate normal to the interface, our results for
the binodals are valid in any dimensions. Last, since R(ρ)
is a bijection, the spinodal region is equivalently defined by
d
2
φ/dR
2
< 0org
0
(ρ) < 0.
To see how our formalism works, let us first consider an
equilibrium case, in which g has an even simpler form
g =
δF
δρ(r)
; F[ρ] =
f (ρ) +
c(ρ)
2
(∇ρ)
2
dr. (8)
Equation (1) is then the Cahn-Hilliard equation for a system
with free energy F [ρ] and mobility M[ρ][44]. Equation (8)
is consistent with (7) since it imposes 2λ + κ
= 0 so that
R = ρ (up to an additive and a multiplicative constant which
do not affect the phase equilibria) and F = F . We recover
φ(R) = f (ρ) as the bulk free energy density, g
0
(ρ) = f
(ρ)
as the chemical potential, and h
0
(ρ) = f
(ρ)ρ − f (ρ)asthe
pressure.
Our common tangent construction on φ(R), which amounts
to extremizing F, t hus reverts to t he usual one in equilibrium,
but extends beyond this. We now show how our formalism can
be used to derive the phase diagrams of QSAPs and PFAPs.
QSAPs. We consider particles i = 1 ···N, moving at speeds
v
i
along body-fixed directions u
i
, which undergo both con-
tinuous rotational diffusion with diffusivity D
r
and complete
randomization with tumbling rate α. Each particle adapts its
speed v(˜ρ
i
) to the local density
˜ρ
i
(r) =
dr
K(r − r
)ˆρ(r
)dr
(9)
with K(r) an isotropic coarse-graining kernel, and
ˆρ(r) =
i
δ(r − r
i
) the microscopic particle density.
Deriving hydrodynamic equations from microscopics is
generally difficult, even in equilibrium [45]. For QSAPs we
020602-2
GENERALIZED THERMODYNAMICS OF PHASE … PHYSICAL REVIEW E 97, 020602(R) (2018)
0 1 2 3
2
4
6
8
10
12
v
0
/v
1
ρ/ρ
m
(a)
1d RTPs
2d RTPs
2d ABPs
0.6 0.8 1.0 1.2 1.4
R
g
R
R
φ − ¯gR (b)
FIG. 1. (a) Phase diagrams of QSAPs. The solid lines corre-
spond to common tangent constructions on φ(R) (red) or f (ρ)
(black). Dashed lines correspond to the spinodals d
2
φ/dR
2
= 0.
Data points are from simulations of run-and-tumble particles (RTPs;
α =1,D
r
=0) or active Brownian particles (ABPs; α =0,D
r
=1)
in one dimension on lattice or two dimensions in continuous
space. Black triangles correspond to Supplemental Material movies
showing nucleation or spinodal decomposition [46]. For all plots,
v(ρ) = v
0
+
v
1
−v
0
2
[1 + tanh(
ρ−ρ
m
L
f
)], K(r) = exp(−
1
1−r
2
)/Z with Z
a normalization constant, ρ
m
= 200,v
1
= 5,L
f
= 100. (b) Common
tangent construction on φ(R)forv
0
= 20.
can follow the path of [14,33,34], taking a mean-field approx-
imation of their fluctuating hydrodynamics. We first assume
a smooth density field so that the velocity can be expanded
as [46]
v(˜ρ
i
) ≈ v(ρ) +
2
v
(ρ)ρ + O(∇
3
), (10)
where ρ is evaluated at r
i
and
2
=
1
2
r
2
K(r)dr. Follow-
ing [ 33,34], the fluctuating hydrodynamics of QSAPs is then
given by ˙ρ = ∇ · (M∇g +
√
2Mρ)[46], with a unit white
noise vector and
g
0
(ρ) = log(ρv); M = ρ
τv(˜ρ)
2
d
;
κ(ρ) =−
2
v
v
; λ(ρ) = 0, (11)
where d is the number of spatial dimensions. Here, τ ≡[(d −1)
D
r
+ α]
−1
is the orientational persistence time. The mean-
field hydrodynamic equation of QSAPs is then Eq. (1) with
the coefficients in Eq. (11). This hydrodynamic description is
expected to hold whenever the correlation length is sufficiently
small for the mean-field approximation to be valid and the
interfaces are sufficiently smooth so that the gradient expansion
is justified.
To construct the phase diagram, for a given choice of v(ρ),
we first solve Eq. (4)forR(ρ) and use it to obtain both φ(R)
and h
0
(R). The binodals then follow via a common-tangent
construction on φ(R) or, equivalently, by setting equal values
of h
0
and g
0
in coexisting phases. Note that since 2λ + κ
= 0
one has R = ρ.
Figure 1 shows the phase diagrams predicted by our general-
ized thermodynamics and by QSAP simulations. As expected,
the hydrodynamic description works best fairly close to the
critical point (but outside a numerically unresolved Ginzburg
interval where fluctuations cannot be neglected). This is where
interfaces are smoothest and the gradient expansion Eq. (10)
most accurate. To determine precisely the binodals, we choose
a v(ρ) (given in the caption of Fig. 1) such that MIPS occurs
only at large densities, leading to well-separated coexisting
densities. Under these conditions, our mean-field approxima-
tion works very well: the agreement between predicted and
measured binodals is excellent. In contrast, a common tangent
construction on f (ρ) defined by f
(ρ) = g
0
(ρ) as proposed be-
fore [14,33] gives a poorer estimate since it correctly captures
the equality of g
0
in both phases but not that of h
0
. This reminds
us that gradient terms directly influence the coexisting densities
through Eq. (4)—quite unlike the equilibrium case. As an
aside, it is remarkable that for QSAPs we can quantitatively
predict the phase diagram of a microscopic model without any
fitting parameter; something rare even for equilibrium models.
Beyond the quantitative prediction of the phase diagram, our
approach provides insight into the universality of the MIPS
seen for QSAPs. For instance, the phase diagram does not
depend on the kernel K, which enters Eq. (11) only through
the constant
2
which then cancels from Eq. ( 4) defining the
nonlinear transform R(ρ). Likewise, Fig. 1 includes lattice
simulations of QSAPs in 1D where full phase separation
is replaced by alternating domains (whose densities obey
the predicted binodal values), and confirms the equivalence
of continuous (ABP) and discrete (RTP) angular relaxation
dynamics for QSAPs [33,34].
PFAPs. We now consider self-propelled particles, of di-
ameter σ , in 2D, interacting via a short-range repulsive pair
potential V (see [46] for details):
˙
r
i
=−
j
∇
i
V (|r
i
− r
j
|) +
2D
t
ξ
i
+ v
0
u
i
;
˙
θ
i
=
2D
r
η
i
.
Here a microscopic mobility multiplying the first term was set
to unity; u
i
= (cos θ
i
, sin θ
i
), and η
i
,ξ
i
are unit Gaussian white
noises. For simplicity, we only include continuous rotational
diffusion, but we expect our results to stand for tumbles as well
since this difference has been shown to have a negligible effect
on the phase equilibria [34]. MIPS occurs in this system if
the Péclet number Pe = 3v
0
/(σD
r
) exceeds a threshold value
Pe
c
∼ 60 [16–19].
We follow [29,47] to derive a fluctuating hydrodynamics for
the stochastic density ˆρ(r) =
N
i=1
δ(r − r
i
), the deterministic
limit of which gives a coarse-grained equation for the mean
density field. On time scales larger than D
−1
r
, in our phase-
separated setup with a flat interface parallel to ˆy, the dynamics
is given by ˙ρ = ∂
2
x
g [29], with
g([ρ],x) = D
t
ρ +
v
2
0
2D
r
(ρ + m
2
) +
ˆ
I
2
−
v
0
D
t
D
r
∂
x
m
1
+ P
D
;
P
D
=
x
−∞
dx
∂
x
V (r
− r) ˆρ(r
)ˆρ(r)d
2
r
;
ˆ
I
2
=−
v
0
D
r
∂
x
V (r
− r) ˆρ(r
)
ˆ
m
1
(r)d
2
r
. (12)
Here,
ˆ
m
n
=
N
i=1
δ(r − r
i
) cos(nθ
i
) and m
n
=
ˆ
m
n
, where
··· represent averages over noise realizations. The lack
of steady-state current shows g to be uniform in the phase-
separated system, equal to some constant
¯
g.
For homogeneous systems the expression for g in Eq. (12)
reduces exactly to the equation of state (EOS) found previously
for the mechanical pressure P of PFAPs [29]. Thus g is equal
020602-3
SOLON, STENHAMMAR, CATES, KAFRI, AND TAILLEUR PHYSICAL REVIEW E 97, 020602(R) (2018)
between phases, as it was for QSAPs, but now it represents
pressure, not chemical potential. Moreover, Eq. (12) general-
izes the pressure EOS of [29] to inhomogeneous situations.
It can formally be written g = g
0
(ρ(x)) + g
int
([ρ],x) where
g
0
(ρ) is the pressure in a notionally homogeneous system with
average density ρ, and the “interfacial” term g
int
represents all
nonlocal corrections to this. The form of g used in Eq. (1) can
then be viewed as a gradient expansion of Eq. (12) for PFAPs.
One way forward would be to make that expansion (or
perhaps avoid it by using a closed-form ansatz for g
int
), and
then find R( ρ) and φ(R) analytically as was done for QSAPs
above. Here, however, we proceed differently, approximating
instead the local part, g
0
(ρ), of g in Eq. (12) by a well bench-
marked, semiempirical EOS, with parameters constrained
by simulations of uniform phases at Pe = 40 < Pe
c
[46].
We thus retain the exact structure of the nonlocal terms,
g
int
(x) ≡ g([ρ],x) − g
0
(ρ(x)) in Eq. (12), but find them nu-
merically. Although less predictive than knowing such terms
algebraically, our method clearly illustrates how they select
the binodals. Furthermore, g
int
includes all orders in gradient
and hence does not rely on a gradient expansion, contrary to
Eq. (1).
We then proceed as in Eq. (3) but, instead of R,
now using the volume per particle ν ≡ ρ
−1
. The integral
x
x
g
(g − g
0
) ∂
x
νdxthen admits two equivalent expressions
ν
g
ν
[
g
0
(ν) −
¯
g
]
dν =
x
x
g
g
int
∂
x
νdx. (13)
Here g
0
(ν) is the pressure-volume EOS, so that the nonzero
value of the right-hand integral directly quantifies viola-
tion of the Maxwell construction. A fully predictive theory
would evaluate the right-hand side integral and then solve
g
0
(ν
) = g
0
(ν
g
) =
¯
g together with Eq. (13) to obtain the values
of the binodals ν
and ν
g
. In practice, we measure g(x)
numerically via Eq. (12) from which we subtract g
0
(ρ(x)) and
integrate over space to obtain the numerical value of the right-
hand side of Eq. (13). Crucially,
¯
g, ν
g
, and ν
are not inputs
here, but are found by solving Eq. (13). Concretely this is done
via a modified Maxwell construction: The binodals correspond
to the intersect between the function g
0
(ν) and a horizontal line
of (unknown) ordinate
¯
g since g
0
(ν
) = g
0
(ν
g
) =
¯
g. We then
adjust the value of
¯
g to solve Eq. (13). This construction is
illustrated in Fig. 2, and is accurately obeyed by simulations,
unlike the equilibrium Maxwell construction, which (notwith-
standing [37]) clearly fails to account for the phase equilibria
of PFAPs where interfacial terms again directly enter.
In this Rapid Commuication, we have shown how to
build a generalized theory of phase-separating scalar active
matter starting from a generalized Cahn-Hilliard description
derived on symmetry grounds. Our work accounts for the
phase equilibria of two important classes of self-propelled
particles, PFAPs and QSAPs, which each undergo MIPS.
In contrast to equilibrium systems, interfacial contributions
to pressure and/or chemical potential generically affect the
binodal densities at coexistence [23,29]. We have given in
Eqs. (2) and (6)anexplicit construction for the binodals
at leading nontrivial order in a gradient expansion. This is
quantitatively accurate for MIPS in QSAPs at high density. In
Eq. (13) we have given a more general construction that holds
1
234
ν
50
100
150
g
0
Pe = 40
Pe = 80
Pe = 120
(a)
0.2 0.4 0.6 0.8
1
1.2 1.4
ρ
40
60
80
100
120
Pe
Measured
Eq. (13)
Equal-area
(b)
FIG. 2. (a) Mechanical pressure of PFAPs. Semiempirical EOS
g
0
(ν) (line) vs numerical measurements (symbols) for various Pe.
Open symbols correspond to binodals and horizontal lines to the
pressure
¯
g predicted by Eq. (13). (b) Corresponding phase diagrams
obtained via the modified (red; see text) and the equal-area (blue)
Maxwell constructions, compared with numerically measured bin-
odals (black). Dashed lines correspond to the spinodals g
0
(ρ) = 0.
Black triangles correspond to Supplemental Material movies showing
nucleation or spinodal decomposition [46].
beyond the gradient expansion; we tested it using numerical
data on PFAPs. In practice, our results are obtained by deriving
the coexisting densities of fully phase-separated profiles in
the steady state. Extending our formalism to account for t he
dynamical convergence to this state is an exciting challenge
left for future works. Similarly, the fate of our generalized
thermodynamic formalism when more than one conserved field
is present is an open question.
Interestingly, QSAPs and PFAPs share the same mathe-
matical structure but their coexisting densities are selected by
equating intensive observables which have different physical
interpretations. In particular, the mechanical pressure is iden-
tical in coexisting phases for PFAPs, but not for QSAPs, due
to the lack of an effective momentum conservation in the latter
case [38]. This fundamental difference is well captured by our
formalism which indeed leads to different observables g for
the two models.
Beyond understanding the phase equilibria of active matter,
we hope that our approach will pave the way toward a more
general definition of intensive thermodynamic parameters
[48–50] for active systems. Building a thermodynamic theory
of active matter would further improve our understanding and
control of these intriguing systems and has become a central
question in the field [3,14,23–25,29,34,36,37,39,51–56].
020602-4