Supporting Information: Mechanical Theory of Nonequilibrium Coexistence and Motility-Induced
Phase Separation
Ahmad K. Omar,
1, 2,
∗
Hyeongjoo Row,
3
Stewart A. Mallory,
4
and John F. Brady
3,
†
1
Department of Materials Science and Engineering, University of California, Berkeley, California 94720, USA
2
Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
3
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA
4
Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802, USA
ACTIVE BROWNIAN PARTICLE CONSERVATION EQUATIONS AND CLOSURES
We provide a derivation of the equations presented in the main text: the conservation equations for density, linear momentum,
and the other orientational moment fields for a collection of
N
interacting active Brownian particles (ABPs). We emphasize
that these equations for interacting ABPs have appeared in various forms throughout the literature [1–3]. However, in addition
to deriving the equations, we will introduce the closures and assumptions that allow us to determine the stationary states of
coexistence.
The Fokker-Planck (or Smoluchowski) equation describing the
N
-body distribution of particle positions and orientations
follows from the particle equations-of-motion and is given by:
∂f
N
(
Γ
,t
)
∂t
+
∑
α
∇
α
·
j
T
α
+
∑
α
∇
R
α
·
j
R
α
= 0
.
(S1a)
Here,
f
N
(
Γ
,t
)
is the probability density of observing a configuration
Γ
≡
(
r
1
,...,
r
N
,
q
1
,...,
q
N
)
at time
t
,
r
α
and
q
α
(
|
q
α
|
= 1)
are the position and orientation vectors of particle
α
,
j
T
α
and
j
R
α
are the translational and rotational fluxes of particle
α
, and
∇
α
=
∂/∂
r
α
and
∇
R
α
=
q
α
×
∂/∂
q
α
are translational and rotational gradient operators. The fluxes are given by:
j
T
α
=
U
0
q
α
f
N
+
1
ζ
F
α
f
N
−
D
T
∇
α
f
N
,
(S1b)
j
R
α
=
−
1
τ
R
∇
R
α
f
N
, ,
(S1c)
where
U
0
is the intrinsic active speed,
ζ
is the translational drag coefficient,
D
T
is the translational Brownian diffusivity (ne-
glected in the main text but provided here for completeness),
τ
R
is the reorientation time scale, and
F
α
is the conservative force
on particle
α
arising from a the potential energy
U
(
r
1
,...,
r
N
)
. We consider pairwise additive isotropic interparticle interac-
tions such that
U
(
r
1
,...,
r
N
) =
∑
α
U
ext
(
r
α
) +
∑
α
∑
β
6
=
α
U
2
(
r
αβ
)
/
2
. Here,
U
ext
(
x
)
is the externally imposed potential at
position
x
,
U
2
(
r
)
is the pair interaction potential, and
r
αβ
=
r
α
−
r
β
(
r
αβ
=
|
r
αβ
|
). The force on particle
α
arising from the con-
servative potentials thus has two contributions
F
α
=
F
ext
α
+
F
C
α
, where
F
ext
α
=
−
∇
α
U
ext
(
r
α
)
,
F
C
α
=
∑
β
6
=
α
F
C
αβ
, and
F
C
αβ
=
−
∇
α
U
2
(
r
αβ
)
. From Eq. (S1) the dynamical operator is defined as
L≡
∑
α
[
∇
α
·
(
−
U
0
q
α
−
ζ
−
1
F
α
+
D
T
∇
α
)+
∇
R
α
·
(
τ
−
1
R
∇
R
α
)]
such that
∂f
N
/∂t
=
L
f
N
.
The full
N
-body distribution function
f
N
allows for the determination of the ensemble average of any observable
O
:
O
=
〈
ˆ
O
(
Γ
)
〉
, where
〈·〉 ≡
∫
(
·
)
f
N
(
Γ
,t
)
d
Γ
is the ensemble average and
ˆ
O
(
Γ
)
is the microscopic definition of the observable
O
.
The time evolution of an observable is then:
∂
O
∂t
=
∫
ˆ
O
(
∂f
N
∂t
)
d
Γ
=
∫
ˆ
O
(
L
f
N
)
d
Γ
=
∫
(
L
?
ˆ
O
)
f
N
d
Γ
=
〈L
?
ˆ
O〉
,
(S2)
where, in the case of Eq. (S1), the adjoint of
L
is
L
?
≡
∑
α
[(
U
0
q
α
+
ζ
−
1
F
α
+
D
T
∇
α
)
·
∇
α
+
τ
−
1
R
∇
R
α
·
∇
R
α
]
.
We first derive an evolution equation of the (number) density field
ρ
(
x
,t
) =
〈
ˆ
ρ
(
x
)
〉
, where
ˆ
ρ
(
x
) =
∑
α
δ
(
x
−
x
α
)
is the
microscopic density of particles. The continuity equation directly follows from this procedure:
∂ρ
∂t
+
∇
·
j
ρ
= 0
,
(S3a)
where:
j
ρ
=
U
0
m
(
x
,t
) +
1
ζ
F
ext
(
x
)
ρ
+
1
ζ
∇
·
σ
C
(
x
,t
)
−
D
T
∇
ρ .
(S3b)
2
Here,
∇
≡
∂/∂
x
is the spatial gradient operator,
j
ρ
is the particle flux,
m
(
x
,t
)
≡ 〈
ˆ
m
(
x
)
〉
is the polarization density field,
ˆ
m
(
x
) =
∑
α
q
α
δ
(
x
−
r
α
)
is the microscopic density of polarization, and:
σ
C
(
x
,t
) =
〈
−
1
2
∑
α
∑
α
6
=
β
r
αβ
F
C
αβ
b
αβ
〉
,
(S4)
is the stress generated by the pairwise interparticle forces and
b
αβ
(
x
;
r
α
,
r
β
) =
∫
1
0
δ
(
x
−
r
β
−
λ
r
αβ
)
dλ
is the bond function [4,
5].
The four terms in the particle flux (S3b) correspond to the four modes of particle transport: transport driven by the active
force, external forcing, interparticle forces, and Brownian motion. Equation (S3b) is also a statement of linear momentum
conservation. To see this we rearrange and find:
0
=
∇
·
σ
(
x
,t
) +
b
(
x
,t
)
,
(S5a)
where the body forces are the terms in Eq. S3b that
generally
may not be expressed as divergences of a tensor:
b
(
x
,t
) =
−
ζ
j
ρ
(
x
,t
) +
ζU
0
m
(
x
,t
) +
F
ext
(
x
)
ρ
(
x
,t
)
,
(S5b)
and the stresses are:
σ
(
x
,t
) =
σ
C
(
x
,t
)
−
ζD
T
ρ
(
x
,t
)
I
,
(S5c)
where
−
ζD
T
ρ
I
represents the Brownian “ideal gas” stress and
I
is the second-rank identity tensor.
In order to solve Eq. (S3), we require the polarization density
m
and the stress
σ
C
. The polarization density has its own
evolution equation. By letting
ˆ
O
=
ˆ
m
in Eq. (S2) and following procedures similar to what were applied in order to obtain
Eq. (S3), we find:
∂
m
∂t
+
∇
·
j
m
+
d
−
1
τ
R
m
= 0
,
(S6a)
j
m
=
U
0
̃
Q
(
x
,t
) +
1
ζ
F
ext
m
+
1
ζ
κ
m
(
x
,t
) +
1
ζ
∇
·
Σ
m
(
x
,t
)
−
D
T
∇
m
,
(S6b)
where
j
m
is the polarization flux,
d
is the spatial dimension, and
̃
Q
(
x
,t
)
≡ 〈
ˆ
Q
(
x
)
〉
is the nematic order density field (where
ˆ
Q
(
x
) =
∑
α
q
α
q
α
δ
(
x
−
r
α
)
is the microscopic nematic density). The interparticle forces contribute to the transport of the
polarization in body-force-like (i.e., no divergence)
κ
m
and stress-like
Σ
m
manners, which are defined as:
κ
m
(
x
,t
) =
〈
1
2
∑
α
∑
α
6
=
β
F
C
αβ
(
q
α
−
q
β
)
δ
(
x
−
r
α
)
〉
,
(S7)
Σ
m
(
x
,t
) =
〈
−
1
2
∑
α
∑
α
6
=
β
r
αβ
F
C
αβ
q
α
b
αβ
〉
.
(S8)
Equation (S8) makes clear that
d
S
(
x
)
·
Σ
m
/ζ
is the average transport of polarization due to the interparticle forces acting
across
the infinitesimal area
dS
from the direction of
d
S
. From the definition [Eq. (S7)] of the body-force-like term
κ
m
, we observe that
configurations with
q
α
=
q
β
do not contribute to
κ
m
and configurations with
q
α
=
−
q
β
contribute the most to
κ
m
in magnitude.
This observation is indicative that
κ
m
is correlated with the reduction in the effective active speed
U
eff
due to interparticle
interactions — a pair of particles slow down when they collide head to head but active motion is largely unaffected when
interacting particles are oriented in the same direction. From the scaling analysis of the active pressure
p
act
∼
ρζU
0
τ
R
U
eff
[6]
and recognizing that the active pressure is proportional to the trace of the polarization flux (with
j
m
∼
ρ
〈
̇
xq
〉
) such that
p
act
∼
ρζU
0
τ
R
〈
̇
x
·
q
〉
[2, 7–10], we can indeed identify that
κ
m
is directly related to the reduction in the effective speed of
active transport of polarization. This motivates a constitutive equation
κ
m
=
−
ζ
(
U
0
−
U
m
eff
)
̃
Q
,
(S9)
3
which leads to
j
m
=
U
0
U
m
̃
Q
+
1
ζ
F
ext
m
+
1
ζ
∇
·
Σ
m
−
D
T
∇
m
,
(S10)
where
U
m
eff
=
U
0
U
m
is the effective speed of active polarization transport. The dimensionless quantity
U
m
(
∈
[0
,
1])
represents
the effective speed relative to the intrinsic speed
U
0
and is an equation-of-state depending on the system volume (area) fraction
φ
and activity
`
0
/D
.
U
m
≈
1
when particles move nearly freely at low (
φ
1
) densities and
U
m
≈
0
when particles mobility
is limited due to interparticle interactions.
To close our equations, expressions for
σ
C
,
̃
Q
,
U
m
, and
Σ
m
are required. The nematic order density
̃
Q
follows its own
evolution equation which can again be derived from Eq. (S2) with
ˆ
O
=
ˆ
Q
:
∂
̃
Q
∂t
+
∇
·
j
̃
Q
+
2
d
τ
R
(
̃
Q
−
1
d
ρ
I
)
=
0
,
(S11a)
j
̃
Q
=
U
0
̃
B
(
x
,t
) +
1
ζ
F
ext
̃
Q
+
1
ζ
κ
̃
Q
(
x
,t
) +
1
ζ
∇
·
Σ
̃
Q
(
x
,t
)
−
D
T
∇
̃
Q
.
(S11b)
Here,
j
̃
Q
is the nematic order flux and
̃
B
(
x
,t
)
≡〈
∑
α
q
α
q
α
q
α
δ
(
x
−
r
α
)
〉
. Interparticle interactions again result in body-force-
and stress-like terms in the nematic order flux:
κ
̃
Q
(
x
,t
) =
〈
1
2
∑
α
∑
α
6
=
β
F
C
αβ
(
q
α
q
α
−
q
β
q
β
)
δ
(
x
−
r
α
)
〉
,
(S12)
Σ
̃
Q
(
x
,t
) =
〈
−
1
2
∑
α
∑
α
6
=
β
r
αβ
F
C
αβ
q
α
q
α
b
αβ
〉
.
(S13)
Again, the stress-like term
Σ
̃
Q
(
x
,t
)
represents the average transport of nematic order due to interparticle forces acting across
x
and
κ
̃
Q
is related to the reduction in the effective speed of active transport of the nematic order. We again propose a constitutive
relation:
κ
̃
Q
=
−
ζ
(
U
0
−
U
̃
Q
eff
)
̃
B
.
(S14)
Consequently, the nematic order flux becomes:
j
̃
Q
=
U
0
U
̃
Q
̃
B
(
x
,t
) +
1
ζ
F
ext
̃
Q
+
1
ζ
∇
·
Σ
̃
Q
(
x
,t
)
−
D
T
∇
̃
Q
,
(S15)
where
U
̃
Q
eff
=
U
0
U
̃
Q
is the effective speed of active nematic order transport. The dimensionless quantity
U
̃
Q
is again an
equation-of-state depending on
φ
and
`
0
/D
.
It proves convenient to rewrite our equations with the traceless tensorial orientational moments to exclude the portions that are
dependent on the lower order orientational moments (i.e.,
ρ
and
m
). The traceless nematic order is defined as
Q
=
̃
Q
−
ρ
I
/d
while
B
=
̃
B
−
α
·
m
/
(
d
+ 2)
, where
α
is a fourth-rank isotropic tensor. (In indicial notation,
α
ijkl
=
δ
ij
δ
kl
+
δ
ik
δ
jl
+
δ
il
δ
jk
where
δ
ij
is the second-rank identity tensor.) The polarization flux [Eq. (S10)] and nematic order evolution equation [Eqs. (S11a)
and (S15)] become:
j
m
=
U
0
U
m
(
Q
+
1
d
ρ
I
)
+
1
ζ
F
ext
m
+
1
ζ
∇
·
Σ
m
−
D
T
∇
m
,
(S16)
∂
Q
∂t
+
∇
·
j
Q
+
2
d
τ
R
Q
=
0
,
(S17a)
where the traceless nematic flux
j
Q
=
j
̃
Q
−
1
d
j
ρ
I
is:
j
Q
=
U
0
U
̃
Q
(
B
+
1
d
+ 2
α
·
m
)
−
1
d
U
0
mI
+
1
ζ
F
ext
Q
+
1
ζ
∇
·
Σ
̃
Q
−
1
ζd
∇
·
σ
C
I
−
D
T
∇
Q
.
(S17b)
4
FIG. S1. (a) Schematic illustrating the effect of
Σ
m
to the coexistence effective pressure
P
coexist
. On the
P −
p
C
diagram, the integral
I
[Eq. (S18)] represents the area between the constant activity curve (black) and horizontal line
P
=
P
coexist
. When the stress-like contribution
Σ
m
is neglected,
I
= 0
since the equal-area construction variable is
E
(
ρ
) =
p
C
. However,
I >
0
when
Σ
m
is considered and the
corresponding coexistence effective pressure (red) should be lower than the coexistence effective pressure predicted by
I
= 0
(blue). (b)
Coexistence curves for athermal active Brownian spheres (3d) obtained by the mechanical theory with the stress-like contribution of the
interparticle forces to the polar order
Σ
m
. The parameter
ξ
represents the magnitude of
Σ
m
. The reduced
P
coexist
resulting from
Σ
m
increases the difference between the coexisting densities. However, this difference, even for extreme values of
ξ
, is not significant.
Now, we have coupled evolution equations for the density [Eq. (S3)], polar order [Eqs. (S6a) and (S16)], and nematic order
[Eq. (S17)]. Closing these equations will require expressions for the following unknowns:
B
,
Σ
̃
Q
,
Σ
m
,
σ
C
,
U
m
, and
U
̃
Q
.
Importantly, since our mechanical theory for phase coexistence requires density gradients only up to second order,
B
and
Σ
̃
Q
can be safely discarded as they will contribute third order (
O
(
k
3
)
in the Fourier space) gradient terms. We postulate that the
effective speeds of active transport for polar and nematic orders are identical, i.e.
U
(
`
0
/D,φ
)
≡
U
m
=
U
̃
Q
.
The conservative interaction stress
σ
C
is a familiar quantity. In the absence of spatial gradients, its sole contribution is the
isotropic pressure arising from interparticle interactions
σ
C
=
−
p
C
(
`
0
/D,φ
)
I
, which again will simply depend on activity and
the particle density. The stress will of course also have gradient terms (the Korteweg stresses). However, these Korteweg stresses,
which arise due to the distortion of the pair-distribution function in the presence of density gradients, are significantly smaller
than the gradient terms generated by the active stress. Indeed, the nonisotropic contribution to the Korteweg stress (which has
the same scale as the isotropic contributions, see main text) was measured in Ref. [9] in order to compute the surface tension of
phase-separated ABPs. These Korteweg stresses were found to be negligibly small, scaling as
∼
ζU
0
D
while the active stress
gradient terms scale as
∼
ζU
0
`
0
. As we are interested in scenarios where MIPS occurs (i.e.,
`
0
/D
1
), we neglect the gradient
terms arising from the conservative stress such that
σ
C
=
−
p
C
(
`
0
/D,φ
)
I
.
Finally, we show that including the stress-like contribution of the interparticle forces to the polar order (i.e.
Σ
m
in the polar
order flux [Eq. (S16)]) broadens the binodals predicted by our mechanical theory, yet only inconsequentially. We therefore close
our equations by neglecting the stress-like contribution. We first note that our first coexistence criterion for MIPS
P
(
p
liq
C
) =
P
(
p
gas
C
) =
P
coexist
is always true regardless of whether we include
Σ
m
or not, as
Σ
m
vanishes in regions of homogeneous
density. However, since density gradients generate
Σ
m
, it can in principle alter the second coexistence criterion, the equal-area
construction, by altering
E
.
From the definition Eq. (S8), it is seen that
Σ
m
represents the correlation between the interaction stress
σ
C
and orientation
m
/ρ
. We examine the effect of the correlation by considering a constitutive relation
Σ
m
=
ξ
σ
C
m
/ρ
. Here, a parameter
ξ
(
>
0)
is introduced in order to investigate effects of the magnitude of this term systematically. The effect of
Σ
m
on the binodal is most
clearly seen by considering the following integral:
I
≡
∫
p
liq
C
p
gas
C
[
P
(
p
C
)
−P
coexist
]
dp
C
=
−
`
0
d
−
1
∫
p
liq
C
p
gas
C
d
Σ
m
zzz
dz
dp
C
.
(S18)
When the correlation term is neglected, the integral
I
trivially vanishes and the corresponding equal-area construction variable
is
E
(
ρ
) =
p
C
as discussed in the main text. With the simple model
Σ
m
=
ξ
σ
C
m
/ρ
, we find that:
I
=
ξτ
R
2(
d
−
1)
ζ
∫
p
liq
C
p
gas
C
(
dp
C
dz
)
2
d
dp
C
(
p
C
ρ
)
dp
C
.
(S19)