Phase Transitions and Motility-Induced Phase Separation
Equilibrium Phase Transitions: Fundamental Principles
Unlike the non-equilibrium phase separation phenomena we will discuss later, equilibrium phase transitions form the bedrock of conventional condensed matter physics. These transitions represent dramatic changes in the physical properties of a system when certain conditions, like temperature or pressure, are varied across a critical point. Phase transitions, whether in equilibrium systems or non-equilibrium active matter, are quintessential emergent phenomena where the collective interactions of enormous numbers of microscopic constituents produce entirely new states of matter with qualitatively different properties, properties that don’t exist at the microscopic level. During these transitions, we typically observe spontaneous symmetry breaking, where the system collectively “chooses” one of several equivalent states, creating order from disorder in a way that cannot be predicted from individual particle behavior. Perhaps most strikingly, systems with completely different microscopic details can exhibit identical critical behavior when they belong to the same universality class, demonstrating how macroscopic physics emerges independently of microscopic details. This non-reducibility often necessitates new theoretical frameworks that operate at higher levels of organization than the underlying microscopic laws.
Classification of Phase Transitions
Equilibrium phase transitions are traditionally classified using Ehrenfest’s approach, which focuses on the behavior of the free energy. The general expression for the Gibbs free energy is:
\[G(T,P) = H - TS = U + PV - TS\]
where \(H\) is enthalpy, \(T\) is temperature, \(S\) is entropy, \(U\) is internal energy, \(P\) is pressure, and \(V\) is volume.
First-order transitions show discontinuities in the first derivatives of the free energy. This manifests as jumps in quantities like volume (\(V = \partial G/\partial P\)) or entropy (\(S = -\partial G/\partial T\)), with familiar examples including water boiling or ice melting. These transitions involve latent heat and exhibit phase coexistence, where both phases can exist simultaneously.
In contrast, second-order (or continuous) transitions show continuous first derivatives but discontinuous second derivatives of the free energy. For example, the heat capacity \(C = -T \frac{\partial^2 G}{\partial T^2}\) may exhibit a discontinuity or divergence at the transition. These transitions are characterized by diverging correlation lengths and scale invariance near the critical point. The magnetization of a ferromagnet near its Curie temperature exemplifies this behavior, with free energy of the form
\[G(T,H) = G_0(T) - H M(T) - \frac{1}{2}\chi(T)H^2 + ...\],
where \(G(T,H)\) is the Gibbs free energy as a function of temperature \(T\) and magnetic field \(H\), \(G_0(T)\) is the field-independent part of the free energy, \(M(T)\) is the magnetization (the order parameter), \(\chi(T)\) is the magnetic susceptibility which measures how the magnetization responds to the applied field, and the susceptibility \(\chi(T)\) diverges at the critical point. Modern physics often simply distinguishes between discontinuous (first-order) and continuous transitions, with the latter revealing fascinating universal properties at the critical point.
Landau-Ginzburg Theory
Landau-Ginzburg theory offers an elegant approach to understanding continuous phase transitions through the concept of an order parameter \(\phi\), which distinguishes between different phases. For a magnet, this might be the magnetization; for a liquid-gas transition, it could be the density difference from the critical value.
The theory is built around a free energy functional:
\[F[\phi] = \int d^d\mathbf{r} \left[ \frac{1}{2}(\nabla\phi)^2 + V(\phi) \right]\]
where the first term \(\frac{1}{2}(\nabla\phi)^2\) represents the energy cost of spatial variations in the order parameter. This gradient term penalizes rapid changes in \(\phi\), reflecting the energetic cost of creating interfaces between different phases. Without this term, the system would fragment into microscopic domains; its presence ensures physically reasonable domain sizes and interface profiles.
The potential \(V(\phi)\) is typically expanded as a power series:
\[V(\phi) = r\phi^2 + u\phi^4 + \ldots\]
(Note: some texts write this as \(\frac{1}{2}r\phi^2 + \frac{1}{4}u\phi^4\); both conventions are in use, and the critical behavior is identical. The minima below \(T_c\) sit at \(\phi = \pm\sqrt{-r/(2u)}\) with the convention used here.)
These even power terms naturally appear in systems with up-down symmetry. For example, in ferromagnets, flipping all spins (\(\phi \to -\phi\)) shouldn’t change the energy. Similarly, in liquid-gas transitions, the potential depends on density fluctuations symmetrically around the critical density.
For systems without this symmetry, odd power terms can appear:
\[V(\phi) = r\phi^2 - h\phi + w\phi^3 + v\phi^4 + \ldots\]
This occurs in ferromagnets under external magnetic fields (the \(h\phi\) term) or in systems near a first-order phase transition (the \(w\phi^3\) term with \(v > 0\)), such as certain structural or metamagnetic transitions. A cubic term combined with a positive quartic coefficient generically drives the transition first-order by tilting one minimum below the other before \(r\) reaches zero.
What makes this formulation powerful is how the parameter \(r\) varies with temperature: \(r \propto (T-T_c)\). When temperature exceeds the critical temperature (\(T > T_c\)), the free energy minimum occurs at \(\phi = 0\), corresponding to the disordered phase. Below the critical temperature (\(T < T_c\)), the minima shift to non-zero values of \(\phi\), representing the ordered phase.
Code
fig, (ax1, ax_critical, ax2) = plt.subplots(1, 3, figsize=get_size(12, 5))
phi = np.linspace(-1.5, 1.5, 1000)
# Above Tc (r > 0)
r = 1.0
u = 1.0
V_above = r * phi**2 + u * phi**4
ax1.plot(phi, V_above, 'b-')
ax1.set_title('above $T_c$ ($r > 0$)')
ax1.set_xlabel('order parameter $\\phi$')
ax1.set_ylabel('potential $V(\\phi)$')
ax1.set_ylim(-1, 3)
ax1.plot(0, 0, 'ro', markersize=5) # Mark the minimum
# At Tc (r = 0)
r = 0.0
u = 1.0
V_critical = r * phi**2 + u * phi**4
ax_critical.plot(phi, V_critical, 'b-')
ax_critical.set_title('at $T_c$ ($r = 0$)')
ax_critical.set_xlabel('order parameter $\\phi$')
ax_critical.set_ylabel('potential $V(\\phi)$')
ax_critical.set_ylim(-1, 3)
ax_critical.plot(0, 0, 'ro', markersize=5) # Mark the minimum
# Below Tc (r < 0)
r = -1.0
u = 1.0
V_below = r * phi**2 + u * phi**4
ax2.plot(phi, V_below, 'b-')
ax2.set_title('below $T_c$ ($r < 0$)')
ax2.set_xlabel('order parameter $\\phi$')
ax2.set_ylabel('potential $V(\\phi)$')
ax2.set_ylim(-1, 3)
# Mark the minima
phi_min = np.sqrt(-r/(2*u))
V_min = r * phi_min**2 + u * phi_min**4
ax2.plot(phi_min, V_min, 'ro', markersize=5)
ax2.plot(-phi_min, V_min, 'ro', markersize=5)
plt.tight_layout()
plt.show()The system naturally evolves toward states that minimize this free energy, following the Euler-Lagrange equation:
\[-\nabla^2\phi + \frac{\partial V}{\partial \phi} = 0\]
Near the critical point, the system exhibits remarkable scaling behaviors that define its critical properties. Let’s explore these key behaviors in more detail:
Correlation Length (\(\xi\)): This measures how far the influence of one microscopic variable extends to others. Mathematically, it can be calculated from the spatial correlation function \(G(r) = \langle \phi(0) \phi(r) \rangle - \langle \phi \rangle^2\), which typically decays exponentially as \(G(r) \sim e^{-r/\xi}\) away from the critical point. Near the critical point, it diverges according to \(\xi \sim |T-T_c|^{-\nu}\), meaning fluctuations become correlated over macroscopic distances.
Order Parameter (\(\phi\)): This quantity distinguishes between phases - like magnetization in magnets or density difference in fluids. It scales as \(\phi \sim |T-T_c|^{\beta}\) below \(T_c\), gradually vanishing as temperature approaches the critical point.
Susceptibility (\(\chi\)): This measures how responsive the system is to external fields (like magnetic susceptibility). It diverges as \(\chi \sim |T-T_c|^{-\gamma}\), indicating the system becomes extremely sensitive to small perturbations.
Specific Heat (\(C\)): This thermodynamic quantity shows how much energy is needed to raise the temperature, scaling as \(C \sim |T-T_c|^{-\alpha}\).
These power laws are characterized by critical exponents (\(\alpha\), \(\beta\), \(\gamma\), \(\nu\)) that aren’t independent but satisfy scaling relations such as \(\alpha + 2\beta + \gamma = 2\) and the hyperscaling relation \(\nu d = 2 - \alpha\) (where \(d\) is the spatial dimension). These relations emerge from the fundamental scaling properties of the free energy. Note that the hyperscaling relation \(\nu d = 2 - \alpha\) holds only below the upper critical dimension \(d_c\) (for the Ising universality class, \(d_c = 4\)); above \(d_c\) mean-field theory applies and hyperscaling breaks down.
Universality
Perhaps the most remarkable aspect of phase transitions is universality. Seemingly different physical systems, like magnets, fluids, or alloys, can exhibit identical critical behavior if they share the same symmetries and dimensionality. This universality is reflected in matching sets of critical exponents, regardless of microscopic details.
The Ising Model
The Ising model provides perhaps the simplest theoretical framework exhibiting a continuous second order phase transition:
\[H = -J\sum_{\langle i,j \rangle} \sigma_i \sigma_j - h\sum_i \sigma_i\]
Where \(\sigma_i = \pm 1\) represents a spin at site \(i\), \(J\) is the coupling constant, and \(h\) is an external field. The model describes a ferromagnetic transition with order parameter \(M = \frac{1}{N}\sum_i \sigma_i\) (the magnetization).
The Ising model has been realized experimentally in various magnetic systems:
- Uniaxial ferromagnets such as DyAlO\(_3\) demonstrate critical behavior closely matching 3D Ising model predictions
- Thin films of ferromagnetic materials often exhibit 2D Ising behavior
- CoNb\(_2\)O\(_6\) (cobalt niobate) provides a rare example of a 1D Ising system
The 2D Ising model was solved analytically by Onsager, revealing a critical temperature \(T_c = \frac{2J}{k_B\ln(1+\sqrt{2})}\) and critical exponents: \(\beta = 1/8\), \(\gamma = 7/4\), and \(\nu = 1\). The specific heat exponent is \(\alpha = 0\), but this does not mean the specific heat is finite at \(T_c\); instead it signals a logarithmic divergence \(C \sim \ln|T-T_c|\), qualitatively different from the power-law divergence (\(\alpha > 0\)) or the jump discontinuity (\(\alpha = 0\) in mean-field theory).
Code
# Plot the susceptibility and order parameter for the 2D Ising model
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=get_size(12, 6))
# Temperature range (normalized by J/k_B)
T = np.linspace(1.8, 2.8, 1000)
T_c = 2.27 # Critical temperature for 2D Ising model
# Order parameter (magnetization) with critical exponent beta = 1/8
m = np.zeros_like(T)
beta = 1/8
for i, t in enumerate(T):
if t < T_c:
m[i] = (1 - t/T_c) ** beta
else:
m[i] = 0
# Susceptibility with critical exponent gamma = 7/4
chi = np.zeros_like(T)
gamma = 7/4
A = 10 # Amplitude factor for visualization
for i, t in enumerate(T):
if abs(t - T_c) < 0.001: # Avoid division by zero at critical point
chi[i] = A * 100 # Large but finite value at critical point
else:
chi[i] = A * (1/abs(t - T_c))**gamma
if chi[i] > A * 100: # Cap the maximum value for better visualization
chi[i] = A * 100
# Plot order parameter (magnetization)
ax1.plot(T, m, 'b-', lw=2)
ax1.set_xlabel('temperature $T$ (in units of $J/k_B$)')
ax1.set_ylabel('magnetization $|M|$')
ax1.axvline(x=T_c, color='r', linestyle='--', alpha=0.7)
ax1.text(T_c+0.05, 0.5, '$T_c$', color='r')
ax1.set_xlim(1.8, 2.8)
ax1.set_ylim(-0.05, 1.05)
# Plot susceptibility
ax2.plot(T, chi, 'g-', lw=2)
ax2.set_xlabel('temperature $T$ (in units of $J/k_B$)')
ax2.set_ylabel('susceptibility $\\chi$')
ax2.axvline(x=T_c, color='r', linestyle='--', alpha=0.7)
ax2.text(T_c+0.05, 50, '$T_c$', color='r')
ax2.set_xlim(1.8, 2.8)
ax2.set_ylim(-5, 510)
plt.tight_layout()
plt.show()
Ising Model Simulation
Below is an interactive simulation of the Ising model that allows you to explore phase transitions by adjusting key parameters. The simulation visualizes the spin configurations in real-time, demonstrating how the system behavior changes across the critical temperature.
Magnetization: 0
Energy: 0
Critical temperature (2D Ising): Tc ≈ 2.27J/kB
This interactive simulation allows you to explore the Ising model’s phase transition by adjusting three key parameters:
Temperature (T/J): Control thermal fluctuations. Notice how the system organizes into ordered domains below the critical temperature (\(T_c ≈ 2.27\)) and becomes increasingly disordered above it.
Magnetic Field (h/J): Apply an external field that biases spins. At strong fields, the system magnetizes regardless of temperature.
Coupling Strength (J): Adjust the interaction strength between neighboring spins. Higher values promote stronger ordering.
The simulation uses the Metropolis algorithm to sample configurations according to the Boltzmann distribution. The visualization shows spin states (black = up, gray = down) updating in real-time, while displaying the current magnetization and energy of the system.
Try setting the temperature near the critical point (\(T ≈ 2.27\)) and observe the large, fluctuating domains characteristic of critical phenomena. This behavior illustrates the diverging correlation length and enhanced fluctuations discussed in the Ising model section.
Motility Induced Phase Separation
When active particles are experimentally realized from colloidal systems, they can exhibit dynamical clustering phenomena such as display in the movie below and in a number of other studies.
The generic mechanism of this clustering has been termed Motility Induced Phase Separation (MIPS) and describes the spontaneous segregation of self-propelled particles into dense and dilute phases, occurring without any attractive interactions between particles. This counterintuitive phenomenon emerges purely from the interplay between self-propulsion and steric hindrance. It carries signatures of a phase transition when the control parameter, i.e. the speed of the active particles is varied.
The central mechanism of MIPS can be understood through a feedback loop: active particles move slower in crowded regions due to collisions, leading to their accumulation, which further increases local density. This positive feedback mechanism drives phase separation even with purely repulsive interactions, a stark departure from equilibrium phase transitions that typically require attractive forces.
In real experimental systems, the actual active particle interactions are not completely free of additional forces other than steric forces, though the principle phenomenon of MIPS is to some extent preserved.
Mathematical Description
At the microscopic level, MIPS is typically modeled using Active Brownian Particles (ABPs) or even discs with excluded volume interactions. These particles follow dynamics governed by:
\[\dot{\mathbf{r}}_i = v_0 \mathbf{n}_i + \mu \sum_{j \neq i} \mathbf{F}_{ij} + \sqrt{2D_t} \boldsymbol{\xi}_i(t)\]
\[\dot{\theta}_i = \sqrt{2D_r} \eta_i(t)\]
In these equations, \(\mathbf{r}_i\) represents the position of particle \(i\), and \(\mathbf{n}_i = (\cos\theta_i, \sin\theta_i)\) is its orientation unit vector. The equation includes self-propulsion with speed \(v_0\), repulsive interactions \(\mathbf{F}_{ij}\) between particles, and translational diffusion with coefficient \(D_t\) and Gaussian white noise \(\boldsymbol{\xi}_i(t)\). The particle’s orientation evolves through rotational diffusion with coefficient \(D_r = 1/\tau_p\).
The key feedback mechanism in MIPS arises from how particle speed effectively decreases in crowded regions. While \(v_0\) is constant in the equation, actual displacement reduces due to collisions. This can be captured by defining a density-dependent effective speed:
\[v_\text{eff}(\rho) = v_0 \cdot f(\rho)\]
where \(f(\rho)\) decreases with local density \(\rho\), commonly expressed as:
\[f(\rho) = (1 - \rho/\rho_{max})^\alpha\]
with \(\alpha > 0\). This establishes the crucial coupling between motility and density that drives MIPS. Later, we will explore a continuum description that more directly illuminates how this feedback mechanism leads to phase separation.
MIPS Simulation
Below is an interactive simulation of Motility-Induced Phase Separation (MIPS) that demonstrates how self-propelled particles can spontaneously form dense clusters without any attractive interactions. The key parameters you can adjust are:
Filling Fraction: Controls the density of particles in the system. Higher values promote clustering.
Orientational Noise: Determines how randomly particles change their direction. Low noise allows particles to move persistently, facilitating phase separation.
Particle Speed: Sets how quickly particles move. Higher speeds typically enhance phase separation.
Attraction Strength: An additional attraction between the particles due to phoretic fields, for example
Cluster Count: 0
Largest Cluster: 0%
Phase Separation: No
Continuum Theory
A more direct understanding of MIPS can also be obtained with the help of continuum theories such as a hydrodynamic approach that couples the density field \(\rho\) with a polarization field \(\mathbf{P}\):
\[\partial_t \rho = -\nabla \cdot \mathbf{J}\] \[\mathbf{J} = v(\rho)\mathbf{P} - D_t\nabla\rho\] \[\partial_t \mathbf{P} = -D_r \mathbf{P} - \frac{1}{d}\nabla\bigl[v(\rho)\rho\bigr] + \text{higher order terms}\]
Here \(\mathbf{P}\) is the polarization density (the coarse-grained mean orientation field), and the diffusive term \(-D_t\nabla\rho\) in the flux accounts for translational Brownian motion. The density-gradient term \(-\frac{1}{d}\nabla[v(\rho)\rho]\) in the polarization equation is essential: it couples the local swim speed to density gradients and is the term responsible for the MIPS instability.
The remaining input is the density-dependent swim speed \(v(\rho)\), which encodes the speed-density coupling at the heart of MIPS. The mechanism does not rely on a particular functional form; all it requires is that \(v(\rho)\) decreases with density. Whatever the true function looks like in a given system, at low density we can Taylor expand it around \(\rho = 0\) and keep the first-order term:
\[v(\rho) = v(0) + v'(0)\,\rho + \mathcal{O}(\rho^2) \equiv v_0\left(1-\frac{\rho}{\rho_c}\right) + \mathcal{O}(\rho^2)\]
Here \(v_0 = v(0)\) is the free swim speed, and the characteristic density \(\rho_c = -v_0/v'(0)\) is set by the initial slope; it is positive because crowding slows the particles down, \(v'(0) < 0\). The widely used linear law
\[v(\rho) = v_0(1-\rho/\rho_c)\]
is therefore nothing deeper than the first-order expansion of an unknown decreasing function, with all the microscopic physics condensed into the two constants \(v_0\) and \(\rho_c\). Two questions remain: why the speed decreases at all, and what sets the slope \(\rho_c\).
Why does crowding reduce the speed in the first place? When two active particles collide, they do not bounce off each other like billiard balls: their propulsion is persistent, so they keep pushing against each other while hardly advancing. A particle escapes such a stall only when it slides off its collision partner or when rotational diffusion turns its orientation away. During the stall its displacement is essentially zero. Collisions therefore reduce the time-averaged displacement rate, and it is this coarse-grained transport speed, not the microscopic propulsion force (which remains unchanged), that enters the continuum theory. You can observe exactly this in the interactive simulation above: colliding particles keep their orientations and simply fail to advance while they block each other.
A kinetic estimate familiar from gas theory provides the slope, and with it a microscopic interpretation of \(\rho_c\). A particle moving at speed \(v_0\) encounters neighbors at a rate \(1/\tau_\text{free} \approx \rho \sigma v_0\), where \(\sigma\) is the particle diameter (the scattering cross section in two dimensions) and \(\rho\) the number density. If every collision stalls the particle for an average time \(\tau_c\), the fraction of time lost in stalls is \(\approx \rho \sigma v_0 \tau_c\) to first order in the density, and the mean speed is reduced accordingly:
\[v(\rho) \approx v_0\,(1 - \rho \sigma v_0 \tau_c) \qquad \Rightarrow \qquad \rho_c = \frac{1}{\sigma v_0 \tau_c}\]
The stall time fixes the scale of \(\rho_c\): a collision is typically resolved by sliding around the neighbor, which takes \(\tau_c \sim \sigma/v_0\), so that \(\rho_c \sim 1/\sigma^2\), a density of the order of close packing that is independent of the propulsion speed. How far can the first-order truncation be trusted? Remarkably far: simulations of active Brownian disks find \(v(\rho)\) to be close to linear over the entire density range, extrapolating to zero near the close-packing density where motion stops entirely, see Fily & Marchetti, Phys. Rev. Lett. 108, 235702 (2012) and Stenhammar et al., Phys. Rev. Lett. 111, 145702 (2013). Model functions with curvature, such as \(f(\rho) = (1 - \rho/\rho_{max})^\alpha\) introduced in the mathematical description above, reduce to the same linear law at first order, with \(\rho_c = \rho_{max}/\alpha\).
Let’s examine the feedback loop that drives MIPS using a simple, intuitive approach. The key physics can be understood through a direct cause-and-effect chain:
Slow-down in dense regions: Active particles move more slowly in regions of higher density due to more frequent collisions. This is captured by the relation \(v(\rho) = v_0(1-\rho/\rho_c)\), where speed decreases as density increases.
Particle accumulation: When particles move from a low-density to a high-density region, they slow down. This means particles spend more time in dense regions than in dilute ones, naturally leading to further accumulation where density is already high.
Positive feedback: This accumulation further increases the local density, causing even more slowing and even more accumulation - creating a positive feedback loop.
To see how this works quantitatively, we need two ingredients: the fact that active particles accumulate where they are slow, and the density dependence of the speed. Combining the two turns the qualitative loop above into a sharp instability criterion.
Step 1: Active particles accumulate where they are slow. Consider first non-interacting active particles whose swim speed varies in space, \(v = v(\mathbf{r})\), imposed for example by a light pattern that controls the activity. In steady state, the density adjusts to be inversely proportional to the local speed:
\[\rho_\text{stat}(\mathbf{r}) = \frac{c}{v(\mathbf{r})}, \qquad c = \text{const}\]
For run-and-tumble particles in one dimension this relation is exact (the zero-flux steady state requires \(v(x)\rho(x) = \text{const}\)), and it holds for active Brownian particles up to gradient corrections. The physical reason is simple: a particle crosses fast regions quickly and lingers in slow regions, so the time spent, and with it the steady-state density, scales as \(1/v\). Since the large-scale motion of an active particle is diffusive with an active diffusivity \(D_\text{act} = v^2\tau_p/d\), the accumulation law can equivalently be written as \(\rho_\text{stat} \propto 1/\sqrt{D_\text{act}}\). Note that this behavior is impossible in equilibrium: the density of passive Brownian particles follows the Boltzmann distribution, which depends only on energies, and therefore remains uniform no matter how the mobility or diffusivity varies in space. Accumulation controlled by kinetics rather than energetics is the first genuinely non-equilibrium ingredient of MIPS.
Step 2: The speed decreases with density. In an interacting suspension, the speed landscape is not imposed from outside but generated by the particles themselves: collisions reduce the effective propulsion speed, \(v = v(\rho)\) with \(v'(\rho) < 0\), as in the linear form \(v(\rho) = v_0(1-\rho/\rho_c)\) introduced above.
Step 3: The feedback loop closes. A uniform density \(\rho_0\) must satisfy the self-consistency condition of Step 1,
\[\rho_0 = \frac{c}{v(\rho_0)}\]
Now impose a small density perturbation \(\delta\rho(\mathbf{r})\). The local speed responds as
\[v(\rho_0 + \delta\rho(\mathbf{r})) \approx v(\rho_0) + v'(\rho_0)\,\delta\rho(\mathbf{r})\]
Because \(v'(\rho) < 0\), the speed perturbation \(\delta v(\mathbf{r})\) is in antiphase with the density perturbation: the suspension is slowest where it is densest. According to Step 1, the density then relaxes toward the stationary profile dictated by this new speed landscape:
\[\rho_0 + \delta\rho' = \frac{c}{v(\rho_0)+v'(\rho_0)\,\delta\rho} \approx \frac{c}{v(\rho_0)}\left(1 - \frac{v'(\rho_0)}{v(\rho_0)}\,\delta\rho\right) = \rho_0 - \rho_0\,\frac{v'(\rho_0)}{v(\rho_0)}\,\delta\rho\]
The density response to the seed perturbation is therefore
\[\delta\rho' = -\rho_0\,\frac{v'(\rho_0)}{v(\rho_0)}\,\delta\rho\]
Since \(v' < 0\), the response \(\delta\rho'\) has the same sign as the seed \(\delta\rho\): dense regions slow down and thereby collect even more particles. If the response exceeds the seed, \(\delta\rho' > \delta\rho\), each round of the loop amplifies the perturbation and the uniform profile is unstable. This happens when
\[\frac{v'(\rho_0)}{v(\rho_0)} < -\frac{1}{\rho_0} \qquad \Longleftrightarrow \qquad \frac{\mathrm{d}\ln v}{\mathrm{d}\ln \rho} < -1\]
The uniform active fluid thus becomes unstable as soon as the speed decays faster than \(1/\rho\). The criterion can equivalently be written as \(\frac{\mathrm{d}}{\mathrm{d}\rho}[\rho\, v(\rho)] < 0\): instability sets in when the “throughput” \(\rho v\) decreases upon adding particles, which is precisely the condition under which car traffic jams.
For the linear speed law \(v(\rho) = v_0(1-\rho/\rho_c)\) we have \(v'(\rho) = -v_0/\rho_c\), and the criterion becomes
\[\rho_0 > \frac{\rho_c}{2}\]
Beyond half the jamming density the feedback wins, and the system spontaneously separates into dense (slow) and dilute (fast) regions, even though no attractive forces act between the particles. This simple argument neglects translational diffusion and the finite persistence of the orientation, which both smooth density perturbations; the collapsed note below shows that including them turns the criterion into a threshold in activity while reproducing the same instability window.
Incorporating Brownian Motion
To develop a more complete picture of MIPS, we need to incorporate rotational diffusion explicitly. Recall that active particles undergo rotational Brownian motion described by:
\[\dot{\theta}_i = \sqrt{2D_r} \eta_i(t)\]
where \(D_r\) is the rotational diffusion coefficient and \(\eta_i(t)\) is Gaussian white noise.
This rotational diffusion leads to a persistence time \(\tau_p\), the characteristic time over which a particle maintains its direction before randomizing: \(\tau_p = 1/D_r\) in two dimensions (where the orientational correlation decays as \(e^{-D_r t}\)) and \(\tau_p = 1/(2D_r)\) in three dimensions. This timescale critically affects the development of MIPS.
On timescales much longer than \(\tau_p\), a free active particle performs a persistent random walk and behaves diffusively, with an effective diffusion coefficient
\[D_\text{eff} = D_t + \frac{v_0^2 \tau_p}{d}\]
in \(d\) dimensions (in 2D, \(D_\text{eff} = D_t + v_0^2/(2D_r)\)). The second term is the active contribution \(D_\text{act}\) used in Step 1 above; it reflects the persistent motion with persistence length \(l_p = v_0\tau_p\).
Diffusion smooths density perturbations and therefore counteracts the feedback mechanism, an effect the simple argument above ignored. To include it, we return to the hydrodynamic equations of the continuum theory and eliminate the polarization adiabatically: on timescales longer than \(\tau_p\), the polarization has already relaxed, so we can set \(\partial_t \mathbf{P} \approx 0\) and obtain (in 2D, where the relaxation rate is \(D_r = 1/\tau_p\))
\[\mathbf{P} \approx -\frac{1}{d D_r}\nabla\bigl[v(\rho)\rho\bigr]\]
Substituting this into the flux \(\mathbf{J} = v(\rho)\mathbf{P} - D_t\nabla\rho\) and using \(\nabla[v(\rho)\rho] = [v(\rho) + \rho\, v'(\rho)]\nabla\rho\) gives a purely diffusive flux,
\[\mathbf{J} = -D(\rho)\,\nabla\rho, \qquad D(\rho) = D_t + \frac{v(\rho)\bigl[v(\rho) + \rho\, v'(\rho)\bigr]\tau_p}{d}\]
so that the density obeys a nonlinear diffusion equation with a density-dependent effective diffusivity:
\[\partial_t\rho = \nabla \cdot \bigl[D(\rho)\, \nabla \rho\bigr]\]
A small perturbation \(\delta\rho \propto e^{ikx + \sigma t}\) around the uniform state \(\rho_0\) grows at the rate \(\sigma = -D(\rho_0)\,k^2\), so the uniform state is linearly unstable exactly when
\[D(\rho_0) < 0\]
A negative effective diffusivity means that particles flow up density gradients. This is the mathematically precise version of the “anti-diffusion” behind MIPS. Two limits connect to what we already know. For a density-independent speed (\(v' = 0\)), \(D(\rho) = D_\text{eff} > 0\) and the uniform state is always stable. Without translational diffusion (\(D_t \to 0\)), \(D(\rho) < 0\) reduces to \(v + \rho v' < 0\), i.e. \(v'/v < -1/\rho\), which is exactly the feedback criterion derived above.
For the linear speed law \(v(\rho) = v_0(1-\rho/\rho_c)\) the effective diffusivity becomes
\[D(\rho) = D_t + \frac{v_0^2\tau_p}{d}\left(1-\frac{\rho}{\rho_c}\right)\left(1-\frac{2\rho}{\rho_c}\right)\]
The activity term is negative only in the window \(\rho_c/2 < \rho < \rho_c\), and it must in addition overcome the translational diffusion:
\[\frac{v_0^2 \tau_p}{d\, D_t} > \frac{1}{\left(\frac{2\rho}{\rho_c}-1\right)\left(1-\frac{\rho}{\rho_c}\right)}\]
The right-hand side is smallest (equal to \(8\)) at \(\rho = \tfrac{3}{4}\rho_c\), so the instability first appears at a critical activity \(v_0^2\tau_p/(d D_t) = 8\) and density \(\rho = \tfrac{3}{4}\rho_c\): the critical point of MIPS. The ratio of self-propulsion to diffusion is commonly expressed through the Péclet number Pe = \(v_0/(\sigma D_r)\), where \(\sigma\) is the particle size. In summary, the system phase separates when:
- The density lies in the unstable window: \(\rho_c/2 < \rho < \rho_c\)
- The Péclet number exceeds a critical value that depends on density
- The rotational diffusion \(D_r\) is sufficiently low to allow persistent motion
At high Péclet numbers, particles move persistently along their orientation before randomizing, so the accumulation mechanism outruns the smoothing by diffusion. When the density additionally lies in the unstable window, the motility reduction in crowded regions is pronounced enough to sustain the positive feedback loop, and the uniform suspension phase separates.
The phase separation behavior is captured in the binodal curve shown in the phase diagram below, which separates the homogeneous phase from the phase-separated region. The system exhibits a critical point where the distinction between phases vanishes, similar to equilibrium phase transitions. The universality class of the MIPS critical point is an open and actively debated question. Several numerical studies find exponents close to mean-field values (\(\beta \approx 0.5\), \(\gamma \approx 1.0\), \(\nu \approx 0.5\)), while others report Ising-like exponents in two dimensions. The difficulty is that finite-size effects and slow dynamics near the critical point make it hard to extract reliable exponents from simulations or experiments, and it remains unclear whether MIPS belongs to the Ising universality class, a new non-equilibrium class, or displays effective mean-field behavior due to long-ranged fluctuations. What is certain is that MIPS is a non-equilibrium phase transition driven by activity rather than thermodynamic interactions.
Code
v_c = 1.0 # Critical speed
rho_c = 0.2 # Critical density (mean density of the system)
rho_dense_max = 0.7 # Dense-phase density at high speed (near close packing)
rho_dilute_min = 0.02 # Residual gas density at high speed
# Binodal branches: square-root (mean-field) behavior near the critical
# point, smoothly saturating at the limiting densities far above v_c
def binodal(v):
x = np.sqrt(np.maximum(v - v_c, 0) / v_c)
high = rho_c + (rho_dense_max - rho_c) * np.tanh(1.2 * x)
low = rho_c - (rho_c - rho_dilute_min) * np.tanh(1.2 * x)
return low, high
speed_range = np.linspace(v_c, 3, 500)
low, high = binodal(speed_range)
plt.figure(figsize=get_size(8, 6))
# Homogeneous stationary density below the critical speed
plt.plot([0, v_c], [rho_c, rho_c], 'k-', lw=1)
# Coexistence region and binodal branches above the critical speed
plt.fill_between(speed_range, low, high, color='0.85', label='coexistence region')
plt.plot(speed_range, high, 'k-', lw=1, label='dense phase')
plt.plot(speed_range, low, 'k--', lw=1, label='dilute phase')
plt.plot(v_c, rho_c, 'ro', markersize=3)
plt.xlabel('particle speed $v_0$')
plt.ylabel('density $\\rho$')
plt.xlim(0, 3)
plt.ylim(0, 0.8)
plt.legend(loc='upper left',fontsize=5)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
Several key features become evident from the binodal curve shown above, which separates the homogeneous region from the phase-separated region. These features include:
Critical Point: At a critical speed \(v_c\) and density \(\rho_c\), the system undergoes a continuous phase transition.
Binodal Curve: Above \(v_c\), the system separates into coexisting phases with densities following the binodal branches.
Coexistence Region: Between the binodal branches, the system separates into regions of high and low density.
This behavior parallels liquid-gas phase transitions in equilibrium systems, but with particle speed rather than temperature serving as the control parameter. The key difference is that MIPS arises purely from the feedback between local density and particle mobility, without requiring attractive interactions.
Number density fluctuations in equilibrium
In equilibrium systems, number density fluctuations follow a characteristic scaling law that serves as a fundamental signature of thermodynamic equilibrium. This behavior emerges directly from the central limit theorem of statistics. When we observe a subsystem containing an average of \(N\) particles, the particles are either statistically independent (in non-interacting systems) or exhibit only short-range correlations. As a result, counting particles in this volume effectively sums many weakly correlated random variables.
The central limit theorem tells us that such a sum produces fluctuations with a variance that scales linearly with the system size:
\[\langle(\delta N)^2\rangle \propto \langle N\rangle\]
This means the standard deviation of these fluctuations scales as \(\sqrt{\langle N\rangle}\). When we examine relative fluctuations, the ratio of fluctuation magnitude to the mean, we find they scale as
\[\delta N/\langle N\rangle \propto 1/\sqrt{\langle N\rangle}\]
This \(1/\sqrt{N}\) scaling holds across equilibrium systems away from any critical point, appearing in ideal gases, liquids, and even interacting solids. It reflects the absence of long-range correlations that would otherwise amplify fluctuations beyond this scaling.
It is important to be precise here, because the \(\sqrt{N}\) law is not universal: it breaks down at a classical critical point. Statistical mechanics relates the number fluctuations exactly to the isothermal compressibility \(\kappa_T\),
\[\frac{\langle(\delta N)^2\rangle}{\langle N\rangle} = \rho\, k_B T\, \kappa_T\]
Away from \(T_c\) the compressibility is finite, so the variance is linear in \(\langle N\rangle\) and we recover \(\delta N/\langle N\rangle \propto 1/\sqrt{\langle N\rangle}\). But on approaching a critical point, \(\kappa_T\) diverges (this is the flat free-energy minimum at \(r = 0\) from Landau theory), the correlation length \(\xi\) grows without bound, and the number fluctuations become anomalously large and correlated over macroscopic distances. This is exactly the phenomenon of critical opalescence: a classical, equilibrium system also exhibits giant density fluctuations, but only when it is fine-tuned to its critical point.
Number density fluctuations in MIPS
Active systems exhibiting Motility-Induced Phase Separation (MIPS) demonstrate a profound departure from this equilibrium behavior. In these systems, the self-propulsion of particles creates persistent motion and induces spatial correlations that extend far beyond what’s possible in equilibrium.
The most striking manifestation of this non-equilibrium nature appears in the density fluctuations. Rather than following the equilibrium \(1/\sqrt{N}\) scaling, active systems show “giant number fluctuations” where the standard deviation scales as
\[\Delta N \propto \langle N\rangle^{a}\]
with \(a > 1/2\) (equivalently, the variance \(\langle(\delta N)^2\rangle \propto \langle N\rangle^{2a}\) grows faster than linearly). The equilibrium central-limit value is \(a = 1/2\). This super-linear scaling means that relative fluctuations actually grow with system size rather than diminish.
This anomalous behavior stems from the long-range correlations that activity creates between particle positions and velocities. As particles cluster in MIPS, they don’t distribute randomly but rather organize collectively in response to local density gradients, amplifying density variations instead of smoothing them out.
The sharper distinction from equilibrium is therefore not simply that active systems show large fluctuations, since classical systems do too at their critical point. It is that active systems exhibit these anomalous fluctuations throughout an extended region of parameter space rather than only at a single fine-tuned critical point. Activity sustains the long-range correlations generically, without the need to tune temperature or density to a critical value.
This difference in number fluctuation scaling provides researchers with a powerful experimental signature to distinguish active from passive phase separation, even without tracking individual particle trajectories. By measuring how density fluctuations scale with system size, one can immediately identify the non-equilibrium character of active matter systems, revealing one of the clearest manifestations of their departure from equilibrium physics.
MIPS exemplifies several profound concepts in non-equilibrium statistical mechanics:
Emergent Thermodynamics: Despite being driven by non-equilibrium processes, MIPS exhibits phase coexistence reminiscent of equilibrium systems. This suggests the possibility of constructing effective thermodynamic frameworks for active systems, an open question that connects to fundamental postulates of statistical mechanics.
Entropy Production: Active systems continuously dissipate energy and produce entropy. Quantifying this entropy production and its relationship to pattern formation addresses foundational questions about irreversibility and the arrow of time.
Broken Time-Reversal Symmetry: The equations governing active particles lack time-reversal symmetry, leading to novel dynamical phenomena including odd elasticity, non-reciprocal interactions, and topologically protected states, concepts that parallel developments in quantum many-body physics.
MIPS provides a conceptual bridge between statistical mechanics and several interdisciplinary fields. Cellular aggregation, bacterial biofilms, and tissue formation involve self-propelled biological entities that can exhibit MIPS-like behavior within the realm of Biological Physics. MIPS concepts help explain the behavior of synthetic active colloids, microswimmers, and other designed active materials in Soft Matter research. Pattern formation in active systems exemplifies how simple rules can generate complex, emergent behaviors, a theme that resonates across physics from condensed matter to cosmology in Complexity Science. Simulating MIPS requires numerical techniques that balance microscopic accuracy with computational efficiency, providing practical experience in scientific computing for those interested in Computational Physics.
For further reading, consider these foundational papers:
- Cates, M. E., & Tailleur, J. (2015). Motility-induced phase separation. Annual Review of Condensed Matter Physics, 6, 219-244.
- Ramaswamy, S. (2010). The mechanics and statistics of active matter. Annual Review of Condensed Matter Physics, 1, 323-345.
- Fily, Y. & Marchetti, M. C. Athermal Phase Separation of Self-Propelled Particles with No Alignment. Phys. Rev. Lett. 108, 235702 (2012).