3.1 Free energy functional and mean-field theory

The simplest next step is to postulate the following free energy functional,

where $a=a(T)$ while $b$ and $\unicode[STIX]{x1D705}$ are positive and (for simplicity) independent of temperature. The bulk free energy density of a uniform state, $f(\unicode[STIX]{x1D719})=(a/2)\unicode[STIX]{x1D719}^{2}+(b/4)\unicode[STIX]{x1D719}^{4}$ , is an approximation, inspired by a Taylor expansion in small $\unicode[STIX]{x1D719}$ for weakly demixed states; for a symmetric fluid this contains only even powers. (Note that at given $\bar{\unicode[STIX]{x1D719}}$ any linear term merely adds a constant to $F$ . Were asymmetric interactions present, any cubic term could also be eliminated by an additive shift, $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{C}$ .) This expansion breaks down at large negative $a$ , where it overshoots the saturating asymptote $\unicode[STIX]{x1D719}_{b}-1\sim \exp [-u/k_{B}T]$ , with $u$ a solubilization energy, set by ideal (dilute) solution thermodynamics. A more accurate form is $f(\unicode[STIX]{x1D719})=-u\unicode[STIX]{x1D719}^{2}/2-k_{B}T[\unicode[STIX]{x1D719}\ln \unicode[STIX]{x1D719}+(1-\unicode[STIX]{x1D719})\ln (1-\unicode[STIX]{x1D719})]$ , but the quartic approximation is sufficient for our purposes, and much easier to use in calculations of both the phase diagram and the interfacial tension between phases.

These calculations can be tackled without further approximation by addressing (3.2) using field theoretic methods, including the renormalization group theory, which is essential to understanding the behaviour very close to the critical point at $T_{C}$ (Chaikin & Lubensky Reference Chaikin and Lubensky1995). Such methods take averages over the Boltzmann weight $\exp [-F/k_{B}T]$ of our fluctuating order parameter(s); we do not pursue them here. A simpler alternative is mean-field theory which considers only the most probable states found by minimizing $F$ at fixed global composition  $\bar{\unicode[STIX]{x1D719}}$ .

We first consider states of uniform $\unicode[STIX]{x1D719}(\boldsymbol{r})=\bar{\unicode[STIX]{x1D719}}$ . For such states

For $a>0$ this has a single minimum at $\bar{\unicode[STIX]{x1D719}}=0$ , with positive curvature everywhere. The latter means that, whatever $\bar{\unicode[STIX]{x1D719}}$ is chosen, one cannot lower the free energy by introducing a phase separation. On the other hand, for $a<0$ , $f$ has negative curvature between the spinodals $\pm \unicode[STIX]{x1D719}_{s}$ where $\unicode[STIX]{x1D719}_{s}=(-a/3b)^{1/2}$ (see figure 1). Also it has two symmetric minima at $\bar{\unicode[STIX]{x1D719}}=\pm \unicode[STIX]{x1D719}_{b}$ with $\unicode[STIX]{x1D719}_{b}=(-a/b)^{1/2}$ . For $|\bar{\unicode[STIX]{x1D719}}|<\unicode[STIX]{x1D719}_{b}$ , $F$ is minimized by demixing the uniform state at $\bar{\unicode[STIX]{x1D719}}$ into two coexisting states at $\unicode[STIX]{x1D719}=\pm \unicode[STIX]{x1D719}_{b}$ . A price must be paid to create an interface between these, but the interfacial area scales as $V^{1-1/d}\ll V$ , so that, for a large enough system, this price is always worth paying.

Although so far restricted to symmetric fluid pairs, this calculation is more general than it first appears. For an asymmetric fluid, one expects (to quartic order) additional linear and cubic terms in $f(\unicode[STIX]{x1D719})$ . However, any linear term in $F$ is of the form $\int \unicode[STIX]{x1D719}\,\text{d}\boldsymbol{r}=\bar{\unicode[STIX]{x1D719}}V$ , which is simply a constant set by the global composition. Such an additive constant to $F$ has no physical effects. In contrast, a cubic term $\int (c\unicode[STIX]{x1D719}^{3}/3)\,\text{d}\boldsymbol{r}$ creates an asymmetric phase diagram, which is useful in fitting the model to real fluid pairs, for which some asymmetry is always present. However, it is a simple exercise then to show this cubic contribution can be absorbed by shifts $a\rightarrow a-a_{C}$ and $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{C}$ with $a_{C}=c^{2}/3b$ and $\unicode[STIX]{x1D719}_{C}=-c/3b$ . In other words, at our chosen level of treating $f(\unicode[STIX]{x1D719})$ as a quartic polynomial, the cubic term merely shifts the mean-field critical point to a new position on the phase diagram; measuring $\unicode[STIX]{x1D719}$ and $a$ relative to this new position, nothing has changed.

3.2 Interfacial profile and tension

In equilibrium, our two bulk phases will minimize their mutual surface area; in most geometries, this requires the interface to be flat. To calculate its interfacial tension, we take the surface normal along the $x$ direction so that $\unicode[STIX]{x1D719}(\boldsymbol{r})=\unicode[STIX]{x1D719}(x)$ . The boundary conditions are that $\unicode[STIX]{x1D719}(x)$ approaches $\pm \unicode[STIX]{x1D719}_{b}$ at $x=\pm \infty$ . To find the profile, we minimize $F[\unicode[STIX]{x1D719}]-\unicode[STIX]{x1D706}\int \unicode[STIX]{x1D719}\,\text{d}\boldsymbol{r}$ with these boundary conditions. (Here $\unicode[STIX]{x1D706}$ is a Lagrange multiplier that holds the global composition fixed during the minimization.) The resulting condition,

involves the functional derivative of $F$ , which we denote by

This is called the chemical potential (or more properly the exchange chemical potential) because, up to a factor of molecular volume, it gives the free energy change on replacing an A molecule with a B molecule locally such that, if $\unicode[STIX]{x1D719}$ is incremented by $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}(\boldsymbol{r})$ , then the free energy change is $\unicode[STIX]{x1D6FF}F=\int \unicode[STIX]{x1D707}\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}\,\text{d}\boldsymbol{r}$ .

Equation (3.4) requires that $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D706}$ , which is constant in space. For our symmetric choice of $F[\unicode[STIX]{x1D719}]$ we have $\unicode[STIX]{x1D707}=\text{d}f/\text{d}\unicode[STIX]{x1D719}=0$ in the two bulk phases at density $\pm \unicode[STIX]{x1D719}_{b}$ , so it follows that $\unicode[STIX]{x1D706}=0$ . It is then a good exercise (Chaikin & Lubensky Reference Chaikin and Lubensky1995) to show that, with the boundary conditions already given, the solution for $\unicode[STIX]{x1D719}(x)$ of the ordinary differential equation $\unicode[STIX]{x1D707}(x)=0$ is

Here $\unicode[STIX]{x1D709}_{0}=(-\unicode[STIX]{x1D705}/2a)^{1/2}$ is an interfacial width parameter, $x_{0}$ marks the midpoint of the interface, and the overall sign choice depends on whether the A-rich or B-rich phase occupies the region at large positive $x$ . The interfacial profile is fixed by a trade-off between the penalty for sharp gradients (set by $\unicode[STIX]{x1D705}$ ) and the purely local free energy terms, which, on their own, would be minimized by a profile $\unicode[STIX]{x1D719}(x)$ that jumps discontinuously from one binodal value to the other. A further exercise is to show that the equilibrium interfacial tension $\unicode[STIX]{x1D6FE}_{0}$ , defined as the excess free energy per unit area of a flat interface, obeys (Chaikin & Lubensky Reference Chaikin and Lubensky1995)

3.3 Stress tensor

If the interfacial profile departs from the equilibrium one, a thermodynamic stress $\unicode[STIX]{x1D748}$ will act on the fluid. An important example is when the interface is not flat but curved; under these conditions $\unicode[STIX]{x1D707}$ cannot be zero everywhere. For use in the NSE, we require not the stress tensor directly but the thermodynamic force density $\boldsymbol{s}\equiv \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}$ . Consider now a small incompressible displacement field $\boldsymbol{u}$ : that is, $\boldsymbol{r}\rightarrow \boldsymbol{r}+\boldsymbol{u}(\boldsymbol{r})$ with $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ . Advection of the $\unicode[STIX]{x1D719}$ field by this displacement induces the change $\unicode[STIX]{x1D719}(\boldsymbol{r})\rightarrow \unicode[STIX]{x1D719}(\boldsymbol{r}-\boldsymbol{u})$ . To linear order this gives the increment $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ , from which we find the free energy increment as

where the final form follows by partial integration and incompressibility. (We consider periodic boundary conditions without loss of generality; this eliminates the boundary term.) This result can be compared with the free energy increment caused by a strain tensor $\unicode[STIX]{x1D750}=\unicode[STIX]{x1D735}\boldsymbol{u}$ ,

where the second form again follows by partial integration. Comparison of (3.8) and (3.9) shows that $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\equiv \boldsymbol{s}=-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}$ . This form, which is all we need to know about the stress tensor for the purposes of solving the NSE, could also have been derived from the following expression for the stress tensor itself:

Here $\unicode[STIX]{x1D6F1}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D707}-\mathbb{F}$ is the order parameter contribution to the (local) pressure, better known as the osmotic pressure; we have defined $\mathbb{F}=f(\unicode[STIX]{x1D719})+\unicode[STIX]{x1D705}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})^{2}/2$ , which is the local free energy density, i.e. $\int \mathbb{F}\,\text{d}\boldsymbol{r}=F$ . The non-gradient part of $\unicode[STIX]{x1D6F1}$ is in turn $\unicode[STIX]{x1D6F1}_{bulk}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D707}_{bulk}-f$ with $\unicode[STIX]{x1D707}_{bulk}\equiv \text{d}f/\text{d}\unicode[STIX]{x1D719}$ ; this is the usual thermodynamic relation between pressure and free energy density (Chaikin & Lubensky Reference Chaikin and Lubensky1995). In some cases it is possible to think of the flow resulting from the force density $\boldsymbol{s}$ as loosely arising from ‘a gradient of osmotic pressure’. This can be useful for estimating the size of terms in the NSE, but cannot be interpreted too literally, since a term $-\unicode[STIX]{x1D735}P$ already appears there to enforce incompressibility. The irrotational term in the order parameter force density $\boldsymbol{s}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}$ can be absorbed into this and so has, strictly speaking, no role in the dynamics.

5.1 Spinodal decomposition

We start by considering the spinodal instability. Ignoring advection and noise initially, we write (4.5) as

Here $f(\unicode[STIX]{x1D719})$ is the local free energy density of a uniform state as defined in (3.3), and primes denote differentiation of this function with respect to $\unicode[STIX]{x1D719}$ . Next we linearize this equation about a uniform initial composition $\bar{\unicode[STIX]{x1D719}}$ , and Fourier transform, to give

where we have defined a wavevector-dependent decay rate $r(q)$ . For $f^{\prime \prime }(\bar{\unicode[STIX]{x1D719}})>0$ , this is positive for all $q$ : all Fourier modes decay and the initial state is stable. In contrast, for $f^{\prime \prime }(\bar{\unicode[STIX]{x1D719}})<0$ , the system is unstable, with $r(q)$ negative at small and intermediate wavevectors. (Stability is restored at high enough $q$ by the $\unicode[STIX]{x1D705}$ term.) Solving $\text{d}r/\text{d}q=0$ for $q$ identifies the fastest-growing instability to be at $q^{\ast }=-f^{\prime \prime }(\bar{\unicode[STIX]{x1D719}})/2\unicode[STIX]{x1D705}$ . Because coupling to fluid flow was neglected, this is technically a model B result, but in fact such coupling (model H) does not change the result of this linear stability analysis.

Even neglecting the thermal noise in the dynamics (4.7) and (4.8), the initial condition can be assumed to have some fluctuations. Those whose wavenumbers lie near $q^{\ast }$ grow exponentially faster than the rest, so that the time-dependent composition correlator $S_{q}(t)\equiv \langle \unicode[STIX]{x1D719}_{\boldsymbol{q}}(t)\unicode[STIX]{x1D719}_{-\boldsymbol{q}}(t)\rangle$ soon develops a peak of height growing as $\exp [-r(q^{\ast })t]$ around $q^{\ast }$ . Hence, during this ‘early stage’ of spinodal decomposition, a local domain morphology is created by compositional interdiffusion with a well-defined length scale set by $\unicode[STIX]{x03C0}/q^{\ast }$ . The amplitude of these compositional fluctuations grows until local values approach the binodals $\pm \unicode[STIX]{x1D719}_{b}$ . From this emerges a domain pattern, still initially with the same length scale, consisting of phases in local coexistence, separated by relatively sharp interfaces of width $\unicode[STIX]{x1D709}_{0}$ and interfacial tension $\unicode[STIX]{x1D6FE}_{0}$ as found from the equilibrium interfacial profile (3.6).

The next stage of the phase separation depends crucially on the topology of the newly formed fluid domains. This is controlled mainly by the phase volumes $\unicode[STIX]{x1D6F7}_{A,B}$ of the A-rich and B-rich phases. Roughly speaking, if $0.3\leqslant \unicode[STIX]{x1D6F7}_{A}\leqslant 0.7$ the domains remain bicontinuous: one can trace a path through the A-rich phase from one side of the sample to the other, and likewise for the B-rich phase. Outside this window, the structure instead has droplets of A in B ( $\unicode[STIX]{x1D6F7}_{A}<0.3$ ) or B in A ( $\unicode[STIX]{x1D6F7}_{A}>0.7$ ). The values $0.3$ and $0.7$ are rule-of-thumb figures only, with details depending on many other factors (including any asymmetry in viscosity or mobility) that we do not consider here. Note also that the window of bicontinuity shrinks to zero width in two dimensions, where the slightest asymmetry in phase volume and/or material properties will generally result in a droplet geometry in which only one phase is continuous. We next assemble the tools needed to understand the late stages of phase separation before discussing their dynamics.

5.2 Laplace pressure of curved interfaces

Once the interfaces are sharp, their geometry becomes a dominant factor in the time evolution. Unless they are perfectly flat, interfaces exert forces on the fluid via the $\boldsymbol{s}=-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}$ term in (4.3), in response to which the fluid may or may not be set in motion. (No motion is implied if $\boldsymbol{s}$ remains irrotational, as discussed previously.) The physics of this term, for interfaces that are locally equilibrated but not flat, is that of Laplace pressure. Consider, for example, a spherical droplet of one fluid in another with radius $R$ and interfacial tension $\unicode[STIX]{x1D6FE}$ . Suppose that the pressure inside the droplet is greater than that outside by an amount $\unicode[STIX]{x0394}P$ . The total force on the upper half of the droplet exerted by the bottom half is then (in three dimensions)

which must vanish if the droplet is not moving. The first term comes from the vertical component of the extra pressure acting across the equatorial disc, and the second is the tension acting across its perimeter. The resulting pressure excess, known as ‘Laplace pressure’, is $\unicode[STIX]{x0394}P=2\unicode[STIX]{x1D6FE}/R$ . More generally, in three dimensions one has $\unicode[STIX]{x0394}P=\unicode[STIX]{x1D6FE}H=\unicode[STIX]{x1D6FE}({R_{1}}^{-1}+{R_{2}}^{-1})$ , where $H$ is the mean curvature and $R_{1}$ and $R_{2}$ are the principal radii of curvature. Note also that for a (hyper-)sphere in $d$ dimensions, $\unicode[STIX]{x0394}P=(d-1)\unicode[STIX]{x1D6FE}/R$ . In all these expressions we can set $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{0}$ as calculated in § 3.2, so long as curvature is weak ( $H\unicode[STIX]{x1D709}_{0}\ll 1$ ) and there is local thermodynamic equilibrium across the interface.

A closely related result concerns the chemical potential $\unicode[STIX]{x1D707}$ at a representative point on a weakly curved interface. Local equilibrium fixes the interfacial profile as $\unicode[STIX]{x1D719}(\boldsymbol{r})\simeq \unicode[STIX]{x1D719}_{0}(w)\equiv \unicode[STIX]{x1D719}_{b}\tanh (w/\unicode[STIX]{x1D709}_{0})$ , where $w$ is a Cartesian coordinate normal to the interface that vanishes at the midpoint of the profile. Then $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=\unicode[STIX]{x2202}_{w}\unicode[STIX]{x1D719}(w)\hat{\boldsymbol{w}}$ with $\hat{\boldsymbol{w}}$ the unit normal (pointing from low to high $\unicode[STIX]{x1D719}$ ). It follows that $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}=\unicode[STIX]{x2202}_{w}^{2}\unicode[STIX]{x1D719}+\unicode[STIX]{x2202}_{w}\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{w}}$ . From the definition of $\unicode[STIX]{x1D707}$ as $\unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}=f^{\prime }(\unicode[STIX]{x1D719})-\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}$ we have

For the local equilibrium profile $\unicode[STIX]{x1D719}(\boldsymbol{r})=\unicode[STIX]{x1D719}_{0}(w)=\unicode[STIX]{x1D719}_{b}\tanh (w/\unicode[STIX]{x1D709}_{0})$ , the first two terms on the right cancel. Moreover, local equilibrium requires that $\unicode[STIX]{x1D707}$ is almost constant on the scale of the interfacial thickness $\unicode[STIX]{x1D709}_{0}$ . Denoting this locally constant interfacial value by $\unicode[STIX]{x1D707}_{I}$ , multiplying by $\unicode[STIX]{x2202}_{w}\unicode[STIX]{x1D719}$ and integrating through the interface gives $2\unicode[STIX]{x1D719}_{b}\unicode[STIX]{x1D707}_{I}=H\unicode[STIX]{x1D705}\int (\unicode[STIX]{x2202}_{w}\unicode[STIX]{x1D719})^{2}\,\text{d}w$ , where $H=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{w}}$ is the mean curvature, defined as above to be positive when the interface curves towards the high- $\unicode[STIX]{x1D719}$ phase. Using the previously quoted result $\unicode[STIX]{x1D6FE}_{0}=\int \unicode[STIX]{x1D705}(\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}_{0}(x))^{2}\,\text{d}x$ , we find that $\unicode[STIX]{x1D707}_{I}=\unicode[STIX]{x1D6FE}_{0}H/2\unicode[STIX]{x1D719}_{B}$ (Bray Reference Bray1994).

Thus for a large spherical droplet of the high- $\unicode[STIX]{x1D719}$ phase in the low- $\unicode[STIX]{x1D719}$ phase (so that $H=2/R$ ), the chemical potential $\unicode[STIX]{x1D707}_{I}$ at the interface has a small positive value proportional to the Laplace pressure. If the chemical potential in the surrounding phase takes its equilibrium value far away ( $\unicode[STIX]{x1D707}=0$ at coexistence), we require a gradient of $\unicode[STIX]{x1D707}$ and thus a constant outward flux of $\unicode[STIX]{x1D719}$ from the droplet surface. Stable equilibrium of a finite droplet surrounded by the opposite phase is possible, but only in a finite system that has $\unicode[STIX]{x1D707}(\boldsymbol{r})=\unicode[STIX]{x1D707}_{I}$ everywhere so that $\boldsymbol{J}=-M\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}=\mathbf{0}$ . This arises when the phase volumes $\unicode[STIX]{x1D6F7}_{A,B}$ are sufficiently asymmetric that the interfacial area of a droplet configuration is less than that of a slab geometry, which depends on the shape of the container and its boundary conditions; even if so, $R\rightarrow \infty$ and $\unicode[STIX]{x1D707}_{I}\rightarrow 0$ in the large system limit. In all other situations (non-spherical domains, spherical droplets of unequal sizes, or a droplet in an infinite bath), there is a flux onto or off the interface proportional to the local curvature.

To quantify this, consider the state of stable equilibrium for a droplet of A-rich phase ( $\unicode[STIX]{x1D719}\simeq +\unicode[STIX]{x1D719}_{b}$ ) in a finite domain of the B-rich phase ( $\unicode[STIX]{x1D719}\simeq -\unicode[STIX]{x1D719}_{b}$ ). After expanding $f(\unicode[STIX]{x1D719})=a\unicode[STIX]{x1D719}^{2}/2+b\unicode[STIX]{x1D719}^{4}/4$ for small deviations about $\pm \unicode[STIX]{x1D719}_{b}$ , equating the Laplace pressure $\unicode[STIX]{x0394}P$ to the difference in osmotic pressure $\unicode[STIX]{x1D719}\unicode[STIX]{x1D707}(\unicode[STIX]{x1D719})-f(\unicode[STIX]{x1D719})$ of the A-rich and B-rich bulk phases gives the result that $\unicode[STIX]{x1D719}_{A,B}=\pm \unicode[STIX]{x1D719}_{b}+\unicode[STIX]{x1D6FF}$ , with both phases equally shifted upwards in composition $\unicode[STIX]{x1D719}$ relative to the equivalent state of bulk coexistence with a flat interface. The upward shift obeys $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FE}_{0}/(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D719}_{b}R)$ , where $\unicode[STIX]{x1D6FC}\equiv f^{\prime \prime }(\pm \unicode[STIX]{x1D719}_{b})=-2a$ and $R$ is the droplet radius. (This $\unicode[STIX]{x1D6FF}$ would reverse in sign for a B-rich droplet in an A-rich phase.) The result can alternatively be derived by setting $\unicode[STIX]{x1D707}_{I}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FF}$ in the equation $\unicode[STIX]{x1D707}_{I}=\unicode[STIX]{x1D6FE}_{0}H/2\unicode[STIX]{x1D719}_{B}$ quoted above.

We next assume that these local equilibrium conditions apply across the curved droplet interface in an infinite B-rich reservoir whose composition far away is $-\unicode[STIX]{x1D719}_{b}+\unicode[STIX]{x1D716}$ , where $\unicode[STIX]{x1D716}<\unicode[STIX]{x1D6FF}$ . Here the far-field parameter $\unicode[STIX]{x1D716}$ is known as the ambient supersaturation. We then seek a spherically symmetric quasi-static ( $\dot{\unicode[STIX]{x1D719}}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}=0$ ) exterior solution $\unicode[STIX]{x1D719}(r)=-\unicode[STIX]{x1D719}_{b}+\tilde{\unicode[STIX]{x1D719}}(r)$ , where $\tilde{\unicode[STIX]{x1D719}}(R)=\unicode[STIX]{x1D6FF}$ and $\tilde{\unicode[STIX]{x1D719}}(\infty )=\unicode[STIX]{x1D716}$ . With $\boldsymbol{J}=-M\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}$ , this solution obeys $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D6FC}\tilde{\unicode[STIX]{x1D719}}$ , where $\unicode[STIX]{x1D6FB}^{2}\tilde{\unicode[STIX]{x1D719}}=0$ ; note that the square gradient contribution to $\unicode[STIX]{x1D707}$ vanishes in this geometry. With the given boundary conditions, the result is $\tilde{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D716}+(\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D716})R/r$ . If we now calculate the radial current $J=|\boldsymbol{J}|$ exterior to the droplet, we have $J=-\unicode[STIX]{x1D6FC}M\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D719}}/\unicode[STIX]{x2202}r=\unicode[STIX]{x1D6FC}M(\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D716})R/r^{2}$ ; just outside the droplet itself this becomes $-\unicode[STIX]{x1D6FC}M(\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D716})R$ . By mass conservation (given the compositional jump $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=2\unicode[STIX]{x1D719}_{b}$ across the interface), we then find ${\dot{R}}=-J(R)/2\unicode[STIX]{x1D719}_{b}$  or

This result will be used in § 5.4 below to discuss the process of competitive growth and shrinkage of droplets known as Ostwald ripening.

5.3 Coalescence of droplet states

If the post-spinodal structure is that of droplets, each relaxes rapidly to minimize its area at fixed volume, resulting in a spherical shape. For well-separated droplets, the thermodynamic force $\boldsymbol{s}$ in the neighbourhood of each droplet has radial symmetry and hence is the gradient of a scalar (pressure) contribution. This has no consequences for fluid flow since all such terms are subsumed into an overall pressure set by the incompressibility constraint. Accordingly, there is no net fluid motion in the absence of noise. Thermal noise allows droplets to move around by Brownian motion, resulting in collisions that cause the mean droplet radius $R$ to increase by coalescence. To estimate the time dependence of $R$ , we observe that, for A droplets in B, the mean inter-droplet distance $L$ is of order $R\unicode[STIX]{x1D6F7}_{A}^{-1/3}$ at small $\unicode[STIX]{x1D6F7}_{A}$ ; more generally, $L/R$ is a function of $\unicode[STIX]{x1D6F7}_{A}$ only. Each droplet will collide with another in a time $\unicode[STIX]{x0394}t$ of order $L^{2}/D$ , where $D\simeq k_{B}T/\unicode[STIX]{x1D702}R$ is the droplet diffusivity. Upon collision, two droplets of radius $R$ make a new one of radius approximately $2^{1/3}R$ , causing an increment $\unicode[STIX]{x0394}\ln R\simeq (\ln 2)/3$ . This gives

where the precise factor of $2^{1/3}$ has ceased to matter, and the left-hand side can be approximated as $\text{d}\ln R/\text{d}t={\dot{R}}/R$ . By integration we then obtain the scaling law

This argument assumes that coalescence is diffusion-limited, and shows that in this case Brownian motion will cause indefinite growth of the mean droplet size, culminating in total phase separation.

The assumption of independently diffusing spherical droplets should be reliable at low enough phase volumes $\unicode[STIX]{x1D6F7}_{A}$ of the dispersed phase, where it is the dominant coarsening mechanism so long as the Ostwald process (described next) is suppressed. But at larger $\unicode[STIX]{x1D6F7}_{A}$ , more complicated routes to coalescence, some involving droplet-scale or macroscopic fluid flow, are possible. One of these is so-called ‘coalescence-induced coalescence’ where the shape relaxation post-collision of a pair of droplets creates enough flow to cause another coalescence nearby (Wagner & Cates Reference Wagner and Cates2001). This gives a new scaling ( $R\sim \unicode[STIX]{x1D6FE}_{0}t/\unicode[STIX]{x1D702}$ ), which coincides with one of the regimes described later for the coarsening of bicontinuous structures ((5.16) below). Indeed, it stems from the same balance of forces as will be discussed for that case. Quite recently a new picture has emerged of hydrodynamic coarsening in moderately concentrated droplet suspensions ( $\unicode[STIX]{x1D6F7}_{A}>0.2$ ), which supersedes a longstanding view that (5.9) still holds in this regime. This picture involves mechanical (Marangoni) forces resulting from departures of the interfacial tension from its equilibrium value ( $\unicode[STIX]{x1D6FE}\neq \unicode[STIX]{x1D6FE}_{0}$ ). These departures arise because the presence of neighbours breaks rotational symmetry around any given droplet so that the thermodynamic force density is no longer a pure pressure gradient (Shimuzu & Tanaka Reference Shimuzu and Tanaka2015).

Figure 2.

Growth rate of a droplet as a function of size during the Ostwald process (a) without and (b) with trapped species.

In many droplet emulsions it is possible to inhibit the coalescence step, so that this entire route to phase separation is effectively blocked. For instance, adding charged surfactants can stabilize oil droplets in water against coalescence by creating a Coulombic barrier opposing the close approach of droplet surfaces. (Surfactants are amphiphilic molecules that have a polar or charged head-group and an oily tail; they are widely used to stabilize oil-in-water droplet emulsions.) Steric interactions between surfactant tails can likewise stabilize water-in-oil emulsions. If the rupture of a thin film of the continuous phase between droplets has a high enough free energy barrier, coalescence rates can be reduced to a manageable, if not negligible, level in both cases (Bibette et al. Reference Bibette, Leal-Calderon, Schmitt and Poulin2002).

5.4 Ostwald ripening

Sadly, however, switching off coalescence is not enough to prevent macroscopic phase separation of droplet emulsions. This is because of a process called Ostwald ripening, in which material is transported from small droplets to large ones by molecular diffusion across the intervening continuous phase. The physics of this process follows directly from (5.7). Assuming that the system has not already reached its end-point of full phase separation, the ambient supersaturation $\unicode[STIX]{x1D716}$ in that equation remains finite. The function ${\dot{R}}(R)$ is shown in figure 2(a). This exhibits an unstable fixed point at

Droplets bigger than this grow, those smaller, shrink. This is because the Laplace pressure in each droplet raises the chemical potential at its outer surface by an amount inversely proportional to its radius. This creates chemical potential gradients from small to large droplets, which shrink and grow in response to the fluxes set up by those gradients.

We can now find the scaling of the typical droplet size $\bar{R}$ by assuming this to be comparable (but not exactly equal) to $R_{\unicode[STIX]{x1D716}}$ and conversely that the ambient supersaturation $\unicode[STIX]{x1D716}$ is comparable (but not exactly equal) to the local supersaturation near a typical drop: $\unicode[STIX]{x1D716}\simeq \unicode[STIX]{x1D6FF}(\bar{R})=\unicode[STIX]{x1D6FE}_{0}/(\unicode[STIX]{x1D6FC}\bar{R}\unicode[STIX]{x1D719}_{b})$ . Substituting in (5.7) gives

which results in the scaling law

A more complete theory of Ostwald ripening, known as the Lifshitz–Slyozov–Wagner theory, is reviewed by Onuki (Reference Onuki2002). This not only confirms these scalings but gives detailed information on the droplet size distribution. Note that (5.12) has similar time dependence to (5.9) for coalescence; this stems from the fact that both mechanisms are ultimately diffusive. However, the nature of the diffusing species (droplet in one case, molecule in the other) is quite different, resulting in prefactors that involve unrelated material properties for the two mechanisms.

5.5 Preventing Ostwald ripening

It follows from (5.12) that the Ostwald process can be slowed by reducing the interfacial tension $\unicode[STIX]{x1D6FE}_{0}$ , but so long as this remains positive it cannot be stopped entirely. A more effective approach is to include within the A phase a modest concentration of some species that is effectively insoluble in B. This might be a polymer or, if A is water and B oil, a simple salt. The idea is that the trapped species in the A droplets creates an osmotic pressure, which rises as $R$ falls, hence opposing the Laplace pressure. Treating the trapped species as an ideal solution in A, equation (5.7) is replaced by (Webster & Cates Reference Webster and Cates1998)

where $\hat{\unicode[STIX]{x1D702}}$ (unrelated to the fluid viscosity $\unicode[STIX]{x1D702}$ ) is a combination of molecular parameters and physical constants, and $\unicode[STIX]{x1D708}$ is the number of trapped molecules in the droplet. The final term reflects the extra osmotic pressure of the trapped material. There is now a stable fixed point (see figure 2 b), even for zero ambient supersaturation, at finite droplet size

If we start from a uniform set of droplets of size $R_{0}$ , then $\unicode[STIX]{x1D708}=c4\unicode[STIX]{x03C0}R_{0}^{3}/3$ , with $c$ the initial concentration of the added species in the A-rich phase. Droplets that have shrunk to a size $R_{\unicode[STIX]{x1D708}}$ can coexist with a bulk A-rich phase without shrinking further: the Laplace pressure is balance by the trapped species’ osmotic pressure. Moreover, if the initial size obeys $R_{0}<R_{\unicode[STIX]{x1D708}}$ , the Ostwald process is reversed: if an emulsion of such droplets is placed in contact with a bulk A-rich phase, they will expand up to size $R_{\unicode[STIX]{x1D708}}$ . By adding trapped species one can thus make robust ‘mini-emulsions’ (Landfester Reference Landfester2003) or ‘nanoemulsions’ (Fryd & Mason Reference Fryd and Mason2012) which permanently resist coarsening by the Ostwald process. However, these are still metastable: so long as the tension $\unicode[STIX]{x1D6FE}_{0}$ is positive, the free energy can always be reduced by coalescing droplets to reduce the interfacial area.

5.6 Coarsening of bicontinuous states

For roughly equal volumes of the two coexisting phases, say $0.3\leqslant \unicode[STIX]{x1D6F7}_{A}\leqslant 0.7$ , the domains of A-rich and B-rich coexisting fluids remain bicontinuous. (This applies in the three-dimensional (3D) case; the two-dimensional (2D) case is special and considered later.) This allows coarsening by a process in which the interfacial stresses – or, loosely speaking, Laplace pressure gradients – pump fluid from one place to another. In all but the most viscous fluids, this process is ultimately faster than either coalescence or Ostwald ripening: we will see below that it allows $L(t)$ to grow with a larger power of the time. Here $L(t)$ is now defined as the characteristic length scale of the bicontinuous structure; for example, the inverse of its interfacial area per unit volume. The flow-mediated coarsening is captured by model H, but absent in model B. In the latter, domain growth is controlled by a version of the Ostwald process in which diffusive fluxes in both phases act to flatten out interfaces (reducing their curvature) and move material from thinner to fatter regions of the bicontinuous structure. This leads again to (5.12) for the characteristic domain size, which is now unaffected by adding insoluble species that can move diffusively throughout the relevant phase. In both cases the driving force is interfacial tension.

To describe flow-mediated coarsening, we assume that $L(t)$ is much larger than the interfacial width, and is the only relevant length in the problem, so that in estimating terms in the NSE (4.3) we may write $\unicode[STIX]{x1D735}\sim 1/L$ . The forcing term $\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}$ is then of order $\unicode[STIX]{x1D6FE}_{0}/L^{2}$ , where we have used the result found above, $\unicode[STIX]{x1D707}_{I}=\unicode[STIX]{x1D6FE}_{0}H/2\unicode[STIX]{x1D719}_{B}$ , for the chemical potential on an interface of curvature $H\sim 1/L$ . The fluid velocity is of order $\dot{L}$ , so the viscous term scales as $\unicode[STIX]{x1D702}\dot{L}/L^{2}$ . The inertial terms are $\unicode[STIX]{x1D70C}\dot{\boldsymbol{v}}\sim \unicode[STIX]{x1D70C}\ddot{L}$ and $\unicode[STIX]{x1D70C}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}\sim \unicode[STIX]{x1D70C}(\dot{L})^{2}/L$ . The $\unicode[STIX]{x1D735}P$ term, enforcing incompressibility, is slave to the other terms.

Importantly, once sharp interfaces are present so that $\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D6FE}_{0}/L^{2}$ , the NSE (4.3) involves only three parameters, $\unicode[STIX]{x1D70C},\unicode[STIX]{x1D6FE}_{0}$ and $\unicode[STIX]{x1D702}$ . From these three quantities one can make only one length, $L_{0}=\unicode[STIX]{x1D702}^{2}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FE}_{0}$ , and one time, $t_{0}=\unicode[STIX]{x1D702}^{3}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FE}_{0}^{2}$ . This means that the domain scale $L(t)$ must obey (Siggia Reference Siggia1979; Furukawa Reference Furukawa1985)

where, for given phase volumes, $f(x)$ is a function common to all fully symmetric binary fluid pairs. This scaling applies only to bicontinuous states; as previously described, in states comprising spherical droplets, the forcing term in the NSE is instead subsumed by the pressure term at leading order.

In systems showing fluid-mediated coarsening via (5.15), we expect different behaviour according to whether $L/L_{0}$ is large or small. For $L/L_{0}$ small, it is simple to confirm that the inertial terms in (4.3) are negligible. (Note also that, in any regime where $f(x)$ is a power law, the two inertial terms have the same scaling.) The primary balance in the NSE is then $\unicode[STIX]{x1D702}\dot{L}/L^{2}\sim \unicode[STIX]{x1D6FE}_{0}/L^{2}$ , resulting in the scaling law $L(t)\sim \unicode[STIX]{x1D6FE}_{0}t/\unicode[STIX]{x1D702}$ so that

This is called the viscous hydrodynamic (VH) regime and holds below some crossover value $x=x^{\ast }$ . In contrast, for $x\gg x^{\ast }$ the primary balance in the NSE is between the interfacial and inertial terms. It is a simple exercise then to show that $L(t)\sim (\unicode[STIX]{x1D6FE}_{0}/\unicode[STIX]{x1D70C})^{1/3}t^{2/3}$ so that

This is called the inertial hydrodynamic (IH) regime. In practice, this crossover from VH to IH is several decades wide, and the crossover value rather high: $x^{\ast }\simeq 10^{4}$ (Kendon et al. Reference Kendon, Cates, Pagonabarraga, Desplat and Bladon2001). The high crossover point is less surprising if one calculates a domain-scale Reynolds number,

The crossover value of $Re$ then turns out to be of order 10 (Kendon et al. Reference Kendon, Cates, Pagonabarraga, Desplat and Bladon2001), and the largeness of $x^{\ast }$ is found to stem from a modest constant of proportionality in (5.16). It means that in practice a clean observation of the inertial hydrodynamic regime has only been achieved in computer simulation: in terrestrial laboratory experiments, the domains are by then so large (millimetres to centimetres for typical fluid pairs) that the slightest density difference between A and B causes gravitational terms to dominate. These terms give a body force pulling the two fluids apart along the gravitational axis, greatly speeding phase separation.

Note finally that, in two dimensions, bicontinuity is exceptional. For fully symmetric fluid pairs, it arises only when the phase volumes $\unicode[STIX]{x1D6F7}_{A,B}$ are also symmetric, allowing percolating paths through both fluids to cross the container (but only just). More generally, it arises at only one special phase volume whose value is set by the asymmetry of the two fluids; even if one could create this state, it would typically not be sustained during coarsening (since the meticulous dynamical balance required varies with the scale parameter $x$ ). Alongside Ostwald and coalescence regimes, one can find, at least in simulations, various complex structures involving cascades of nested droplets whose details may be noise-dependent (Gonnella et al. Reference Gonnella, Orlandini and Yeomans1999).

7.1 Interfacial tension in the presence of surfactant

Consider first a solution of surfactant molecules each carrying a unit vector $\hat{\unicode[STIX]{x1D742}}_{i}$ denoting its head/tail orientation; local coarse graining creates a smooth field $\boldsymbol{p}(\boldsymbol{r})=\langle \hat{\unicode[STIX]{x1D742}}\rangle _{meso}$ . In the absence of interfaces, this fluctuating mesoscopic field is zero on average, and its small (hence Gaussian) fluctuations are governed by a free energy contribution

This is defined such that the variance of $\boldsymbol{p}$ is set by the ‘osmotic compressibility’ $\unicode[STIX]{x1D712}(\unicode[STIX]{x1D707}_{s})$ , which is an increasing function of the surfactant chemical potential $\unicode[STIX]{x1D707}_{s}(c_{s})$ , which is itself a function of surfactant concentration $c_{s}$ . We can add this to the binary fluid free energy (3.2), alongside a coupling term to represent the reduction in free energy caused when a surfactant molecule resides at the A–B interface with its orientation suitably aligned along the composition gradient $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ :

It is a simple exercise to minimize this over $\boldsymbol{p}(\boldsymbol{r})$ at fixed $\unicode[STIX]{x1D719}(\boldsymbol{r})$ , and a slightly more complicated one to explicitly integrate over the fluctuating $\boldsymbol{p}$ field by Gaussian integration to obtain $\text{e}^{-\unicode[STIX]{x1D6FD}F[\unicode[STIX]{x1D719}]}=\int \text{e}^{-\unicode[STIX]{x1D6FD}F[\unicode[STIX]{x1D719},\boldsymbol{p}]}{\mathcal{D}}\boldsymbol{p}$ . The result of either calculation is to recover the original model H free energy $F[\unicode[STIX]{x1D719}]$ as in (3.2), but with a renormalized square gradient coefficient $\unicode[STIX]{x1D705}_{r}=\unicode[STIX]{x1D705}-\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D712}$ . From this it follows that the interfacial tension is reduced by adding surfactant, according to

This calculation is based on a Gaussian approximation (7.1) that assumes $\boldsymbol{p}(\boldsymbol{r})$ to deviate only mildly from zero everywhere. It is not very realistic – for instance, the interfacial width diverges as $\unicode[STIX]{x1D6FE}_{0}$ becomes small (whereas in practice this width is set by the size of a surfactant molecule). However, it captures the main physical effects of interest here; for a fuller discussion, see Gompper & Schick (Reference Gompper, Schick, Domb and Lebowitz1994). Note that the sign of $\unicode[STIX]{x1D714}$ , which determines whether $\boldsymbol{p}$ points up or down the interfacial $\unicode[STIX]{x1D719}$ gradient, is irrelevant.

The simplest possible assumption is that the surfactant molecules form an ideal solution, with no interactions between them. In this case $\unicode[STIX]{x1D712}=\tilde{\unicode[STIX]{x1D712}}c_{s}$ , with $\tilde{\unicode[STIX]{x1D712}}$ a constant. (Linearity in $c_{s}$ of the variance parameter $\unicode[STIX]{x1D712}$ then follows from the Poisson statistics of randomly located and oriented molecules.) In this case the interfacial tension vanishes, with infinite slope, at $c_{s}=\tilde{c}=\unicode[STIX]{x1D705}/\unicode[STIX]{x1D714}^{2}\tilde{\unicode[STIX]{x1D712}}$ . But, in fact, surfactant solutions are far from ideal, due to a phenomenon called ‘micellization’ in which individual molecules self-assemble into micelles, which contain several tens of molecules. (Micelles are typically spherical; in water they have the polar heads at the exterior surface and the apolar tails in the interior of the sphere. In oil, this structure is reversed.)

The effect of micellization is to cause the surfactant chemical potential $\unicode[STIX]{x1D707}_{s}$ , and therefore $\unicode[STIX]{x1D712}$ , to rapidly saturate above a so-called ‘critical micelle concentration’ $c_{s}=c^{\ast }$ . Any surfactant added beyond this level becomes sequestered into micellar aggregates. The concentration $c_{1}$ of ‘free’ molecules remains very close to $c^{\ast }$ at all higher concentrations and, to a good approximation, $\unicode[STIX]{x1D712}$ saturates at $\tilde{\unicode[STIX]{x1D712}}c^{\ast }$ , so that $\unicode[STIX]{x1D6FE}_{0}(c_{s})$ does not fall further. For a more detailed discussion of micellization, see Safran (Reference Safran2003) or Cates (Reference Cates2012).

The outcome is to have two general classes of behaviour, depending on whether $c^{\ast }$ lies below or above $\tilde{c}$ . In case 1, $c^{\ast }<\tilde{c}$ , so that $\unicode[STIX]{x1D6FE}_{0}(c_{s})$ follows (7.3) so long as $c_{s}<c^{\ast }$ but then abruptly stops decreasing as micellization intervenes. In case 2, $c^{\ast }>\tilde{c}$ , so that $\unicode[STIX]{x1D6FE}_{0}(c_{s})$ hits zero before micelles are formed. At this point, if water and oil are both present in bulk quantities, the system can minimize its free energy by creating a macroscopic amount of interface on which the surfactant can reside. When this happens, $\unicode[STIX]{x1D707}_{s}$ again saturates: adding further surfactant simply creates more surface at fixed $\unicode[STIX]{x1D707}_{s}$ . Hence $\unicode[STIX]{x1D6FE}_{0}$ does not become finitely negative but remains stuck at an effectively zero value. This analysis is grossly simplified, ignoring (among other things) free energy contributions from interfacial curvature, which can allow interfaces to proliferate even for small positive tension. Nonetheless, the broad distinction between case 1 and case 2 is a useful one.

7.2 Finite tension: metastable emulsions and biliquid foams

Case 1 is the more common: the typical effect of surfactant is to reduce interfacial tension to one-half or one-third of its previous value. The global free energy minimum then comprises complete phase separation, just as it does without surfactant; if emulsions are formed (for instance by stirring), they are at most metastable. Since fluid-mediated coarsening is always present in bicontinuous states, metastability generally requires a droplet geometry; as already mentioned, surfactants can help prevent their coalescence by inhibiting film rupture.

Typically, such emulsions have spherical droplets (of A in B, say) but, by evaporation or drainage under gravity, for instance in a centrifuge, much of the continuous B phase can often be expelled to create a so-called biliquid foam (Bibette et al. Reference Bibette, Leal-Calderon, Schmitt and Poulin2002), in which polyhedral droplets of A (say) are separated by thin films of B. In many cases, biliquid foams can persist for hours or days, and sometimes longer. To achieve this, one must suppress not only the rupture of thin films but also the Ostwald process, which, despite the more complicated geometry, still drives diffusion of A from small (few-sided) to large (many-sided) polyhedral droplets (Weaire & Hutzler Reference Weaire and Hutzler1999). To achieve this with a trapped species requires an especially low level of solubility in the B phase so as to ensure negligible diffusion even across the thin B films present in the foam structure. So long as they remain metastable against rupture and coarsening, biliquid foams, like soap froths, are solid materials (generally amorphous, though ordered examples can be made). As such, they have an elastic modulus, and also a yield stress, both of which scale as $G\sim \unicode[STIX]{x1D6FE}_{0}/R$ , with $R$ the mean droplet size. This is an interesting example of a solid behaviour emerging solely from the spatial organization of locally fluid components – for even the surfactant on the interface is (normally) a 2D fluid film.

7.3 Near-zero tension: stable microemulsions

In case 2, added surfactant can reduce $\unicode[STIX]{x1D6FE}_{0}$ to negligible levels for $c_{s}\geqslant \tilde{c}$ . This can lead to thermodynamically stable emulsions, generally called ‘microemulsions’, in which enough A–B interface is created to accommodate all surplus surfactant, of which the concentration is $c_{s}-\tilde{c}$ . There are two broad approaches to describing the resulting structures. One avenue is to base a description on the $\unicode[STIX]{x1D719},\boldsymbol{p}$ order parameters already introduced, addressing (7.2) in the case where $\unicode[STIX]{x1D705}_{r}=\unicode[STIX]{x1D705}-\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D712}$ , so that the thermodynamic interfacial tension is negative. Unsurprisingly, this model is unstable unless further terms are added to prevent interfacial proliferation; the most natural addition to $F$ is a term in $\int (\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})^{4}\,\text{d}\boldsymbol{r}$ . The competition of this term with the effectively negative square gradient coefficient $\unicode[STIX]{x1D705}_{r}$ sets a characteristic length scale for stable domains of A-rich and B-rich fluids. This approach leads to many insights (Gompper & Schick Reference Gompper, Schick, Domb and Lebowitz1994) but is mainly appropriate for systems in which ‘weak’ surfactants (whose adsorption energy at an interface is not much larger than the thermal energy $k_{B}T$ ) are present at high concentration so that the domain size $L$ and interfacial width $\unicode[STIX]{x1D709}$ are comparable. The structure of the microemulsion can then be viewed as a smooth spatial modulation of composition $\unicode[STIX]{x1D719}(\boldsymbol{r})$ rather than as a system of well-separated, surfactant-coated interfaces.

For strong surfactants, which have small $\tilde{c}$ , one can instead treat almost all the surfactant as interfacial. The interfacial area ${\mathcal{S}}$ of the fluid film then obeys

Here $\unicode[STIX]{x1D6F4}$ is a preferred area per molecule, $\unicode[STIX]{x1D719}_{s}$ is the volume fraction of surfactant and $v_{s}$ its molecular volume. For the soluble surfactants normally used for emulsification, the area per molecule is maintained very close to $\unicode[STIX]{x1D6F4}$ by rapid adsorption and desorption at the interface; the specific interfacial area ${\mathcal{S}}/V$ is then fixed directly by $\unicode[STIX]{x1D719}_{s}$ via (7.4).

With interfacial tension negligible, the energetics of a given structure in case 2 is determined by the cost of bending the interfacial surfactant film at fixed area. We treat this by a leading-order harmonic expansion about a state of preferred curvature set by molecular geometry. By a theorem of differential geometry (David Reference David, Nelson, Piran and Weinberg2004), at each point on the A–B interface one can uniquely define two principal radii of curvature, $R_{1}$ and $R_{2}$ , which we take positive for curvature towards A. (For a spherical droplet of A of radius $R$ , we have $R_{1}=R_{2}=R$ , whereas for a cylinder of radius $R$ , we have $R_{1}=R$ and $R_{2}=\infty$ . A saddle shape has $R_{1}$  and $R_{2}$ of opposite signs.) The harmonic bending energy for a fluid film then reads

There are three material parameters: the elastic constants $K$ and $\bar{K}$ (both with dimensions of energy) and $R_{0}$ , a length defining the preferred radius of mean curvature, whose relation to the molecular geometry of surfactants is discussed by Safran (Reference Safran2003). The integral in (7.5) is over an interface ${\mathcal{S}}$ between phases that can have disconnected parts (droplets) but must be orientable so that A is enclosed by ${\mathcal{S}}$ and B excluded. Therefore, its enclosed volume $V_{in}$ must obey

where we have partitioned the surfactant equally between A and B to allow us to define the volume $V_{in}$ as enclosed by a mathematical surface of no thickness. The phase volume of $V_{in}$ is then $\unicode[STIX]{x1D6F7}$ (with $V_{out}=1-\unicode[STIX]{x1D6F7}$ ), and a completely symmetric state has $\unicode[STIX]{x1D6F7}=1/2$ .

7.4 The physics of bending energy

The statistics of the A–B interface in case 2 is, by the above arguments, determined by the Boltzmann distribution ${\mathcal{P}}[{\mathcal{S}}]\propto \exp [-\unicode[STIX]{x1D6FD}F_{bend}]$ . Performing averages over this distribution is intractable, in general. However, some key concepts can be identified that allow the problem to be understood qualitatively.

7.5 Bicontinuous microemulsions

In so-called ‘balanced’ microemulsions, the spontaneous curvature is tuned to be small, so that $R_{0}\gg \unicode[STIX]{x1D709}_{K}$ . For simplicity, we set $R_{0}\rightarrow \infty$ as above. We also assume $\ell \ll \unicode[STIX]{x1D709}_{K}\ll 100~\unicode[STIX]{x03BC}\text{m}$ , so that the entropy of the interface (including the renormalization of $K$ ) cannot be ignored. Assuming $\unicode[STIX]{x1D6F7}$ of order 0.5 (roughly symmetric amounts of A- and B-rich fluid), we can introduce a structural length scale $\unicode[STIX]{x1D706}$ , which is then set by $\unicode[STIX]{x1D719}_{s}$ . Specifically, for a lamellar phase, comprising a one-dimensional stack of alternating layers of A and B with relatively flat interfaces between these, one has a layer spacing $\unicode[STIX]{x1D706}$ between adjacent surfactant films set by $\unicode[STIX]{x1D719}_{s}\simeq (\ell /\unicode[STIX]{x1D706}+\ell )$ . When $\unicode[STIX]{x1D706}\simeq \ell$ the system has no option but to fill space with flat parallel layers. As $\unicode[STIX]{x1D719}_{s}$ is reduced ( $\unicode[STIX]{x1D706}$ raised), the layer spacing $\unicode[STIX]{x1D706}$ becomes comparable to $\unicode[STIX]{x1D709}_{K}$ . For $\unicode[STIX]{x1D706}/\unicode[STIX]{x1D709}_{K}\leqslant 1/3$ (or so), the lamellar phase fluctuates but remains stable.

On the other hand, if $\unicode[STIX]{x1D719}_{s}$ is then decreased further so that $\unicode[STIX]{x1D706}\simeq \unicode[STIX]{x1D709}_{K}$ , these layers melt into an isotropic phase comprising (for $\unicode[STIX]{x1D6F7}\simeq 0.5$ ) bicontinuous domains of A and B fluids separated by a fluctuating surfactant film. This is the bicontinuous microemulsion and effectively represents a thermodynamic route to prevent coarsening of the transient bicontinuous structures encountered in § 5 above. If $\unicode[STIX]{x1D6F7}$ now deviates strongly from $0.5$ , then (just as found there) the structure depercolates, forming a droplet phase. This differs from the one discussed above for $R_{0}\ll \unicode[STIX]{x1D709}_{K}$ , since this one is stabilized by entropy and fluctuations, not by a preferred curvature of the droplets. Theories of the bicontinuous microemulsion were initiated by de Gennes & Taupin (Reference de Gennes and Taupin1982). Some of these theories use coarse-grained lattice models in which fluid domains are placed at random on a lattice of some scale $\unicode[STIX]{x1D709}$ ; the bending energy and area of the resulting interface can be estimated and used to calculate a phase diagram (Andelman et al. Reference Andelman, Cates, Roux and Safran1987). One specific feature is the appearance of three-phase coexistence in which a ‘middle-phase’ microemulsion coexists with excess phases of both oil and water: this roughly corresponds to a ‘double-sided’ emulsification failure in which both oil and water are expelled.

7.6 The sponge phase and vesicles

Closely related to the bicontinuous microemulsion is the ‘sponge phase’. This arises when there is a huge phase volume asymmetry between A-rich and B-rich fluids, but only in case 2 systems where the surfactant has a molecular preference to form a flat film rather than highly curved structures such as micelles. The interfacial structure that forms spontaneously at $c_{s}=\tilde{c}$ is, in the almost complete absence of B, necessarily now a bilayer with A (usually water) on both sides and a thin B layer in the middle. The lamellar state now consists of flattish bilayers alternating with A domains; for small volume fractions of bilayer, such that their separation exceeds their persistence length, this structure melts, just as described above for the microemulsion. The result is subtly different though: we now have a bilayer film that separates two randomly interpenetrating domains containing the same solvent A. This is called the ‘sponge phase’. It is not directly useful for A/B emulsification, since the minority B phase has negligible phase volume. However, the sponge phase does have remarkable phase transitions associated with the ‘in–out’ symmetry between the two A domains, which can be broken spontaneously (Huse & Leibler Reference Huse and Leibler1988; Roux et al. Reference Roux, Coulon and Cates1992). When the symmetry is strongly broken, one has discrete droplets of A separated from a continuous A phase by a bilayer. This can be viewed as a thermodynamically stabilized A-in-A emulsion (typically water-in-water), which can be useful for encapsulation. Such structures, called vesicles, can also exist in the complete absence of B, so that the bilayer contains only surfactant.

Similar vesicles can also be formed from lipids, which are biological molecules closely related to surfactants. The main difference is that for lipids $\tilde{c}$ is extremely small compared to typical values for surfactants. This renders them effectively insoluble in water, since any attempt to increase $c_{s}$ above this small value leads instead to the creation of more interface on which the lipid resides (typically organized as a lamellar phase). This insolubility means that, although in local equilibrium $\unicode[STIX]{x1D6FE}_{0}$ is effectively zero, there is no rapid equilibration of the surface area per molecule $\unicode[STIX]{x1D6F4}$ by molecular exchange between interface and bulk. As a result, if a lipid vesicle is mechanically stretched, a tension soon develops. Such tension can also arise by swelling the vesicle to an inflated spherical shape rather than a relaxed floppy shape of lesser volume. Only in the absence of such stretching can a lipid bilayer be described by the bending free energy (7.5) in which typically $R_{0}=\infty$ because there is no preferred curvature by symmetry; in this case there are large shape fluctuations so long as $K/k_{B}T$ is not large. The spontaneous curvature radius $R_{0}$ can, however, be finite if the composition of the lipid bilayer is different on its two faces. This is common in biology, where a bilayer of mixed lipids separates the interior and exterior of a cell. Note that, because of their insoluble character, the global distribution of lipids to form vesicles is almost never in equilibrium: the size of each vesicle is set by the amount of lipid present at its surface, not by thermodynamics.

8.1 Colloidal particles at fluid–fluid interfaces

Although Janus beads are often studied (Lattuada & Hatton Reference Lattuada and Hatton2011), they are rarely used for emulsion stabilization, because it is much cheaper to use ordinary colloidal spheres of homogeneous surface chemistry. Perhaps surprisingly, so long as the colloid has roughly equal affinity to the two fluids A and B, interfacial attachment energies remain vastly larger than $k_{B}T$ . The simplest case is when the two solid–fluid interfacial tensions, $\unicode[STIX]{x1D6FE}_{SA}$ and $\unicode[STIX]{x1D6FE}_{SB}$ , are the same. The energy of such a particle (of radius $a$ ) is then independent of where it resides, but the energy of the A–B interface is reduced by $\unicode[STIX]{x03C0}a^{2}\unicode[STIX]{x1D6FE}_{0}$ if the particle is placed there, because a disc of interface is now covered up by the colloid. Typically, $\unicode[STIX]{x1D6FE}_{0}=0.01~\text{N}~\text{m}^{-2}$ , giving an attachment energy of order $10^{7}k_{B}T$ for $a=1~\unicode[STIX]{x03BC}\text{m}$ and $10k_{B}T$ for $a=1$  nm. Similar remarks apply for unequal solid–fluid tensions so long as the contact angle $\unicode[STIX]{x1D703}$ , defined via $\unicode[STIX]{x1D6FE}_{0}\cos \unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FE}_{SB}-\unicode[STIX]{x1D6FE}_{SA}$ for the case of partial wetting (figure 5), is not too close to 0 or $\unicode[STIX]{x03C0}$ , whereupon the particle becomes fully wetted by one or other phase.

Figure 5.

(a) Geometry of a partially wet colloidal particle at the fluid–fluid interface; $\unicode[STIX]{x1D703}$ is the contact angle. This is usually measured through the polar phase, so the upper fluid is water and the lower is oil, as drawn here. (b) The initial locus of a fluid interface (light grey) into which particles (black circles) are inserted. These are jammed in a 2D layer but have interstitial fluid regions between them, as shown. If the volume of fluid in the droplet is now reduced, to maintain a fixed contact angle with particles that cannot move, the curvature of the interface is reversed to give the final fluid locus (dark grey arcs). This creates a negative Laplace pressure, which can switch off the Ostwald process.

Consider now placing a number of spherical colloidal particles on the surface of a spherical droplet of fluid A in fluid B, or vice versa. Each colloid can be accommodated with the required contact angle $\unicode[STIX]{x1D703}$ by cutting a small spherical cap out of the interface and slotting the particle into place there. The contact line is a perfect circle, as required for tangency at fixed angle to a sphere. To conserve the enclosed volume, the droplet radius changes slightly, but it remains perfectly spherical. There is minimal change in Laplace pressure, and the interfacial energy is independent of where the colloids are placed. Accordingly, there is no capillary force between them. These statements can only change if particles become jammed so they interact directly via particle–particle forces.

The physics is very different for spherical colloids residing on a hypothetical cylinder of fluid A in B (or vice versa). It is not possible now to insert a spherical particle onto the surface of this cylinder at fixed contact angle $\unicode[STIX]{x1D703}$ unless the fluid interface becomes deformed. This deformation costs extra interfacial energy, and can be minimized by placing two particles close together. Accordingly, there is a capillary attraction between the colloids, which will have a strong tendency to aggregate. Similar arguments apply to non-spherical particles such as ellipsoids, even on flat surfaces (Cavallaro et al. Reference Cavallaro, Botto, Lewandowski, Wang and Stebe2011).

In all these cases, particle adsorption to the interface is thermally irreversible for particles bigger than a few tens of nanometres. This non-equilibrium interfacial physics allows stabilization not only of spherical droplets (classically known as ‘Pickering’ or ‘Ramsey’ emulsions) but also of more complex structures with frozen shapes maintained by a layer of interfacial particles clamped together by interfacial tension (Binks & Horozov Reference Binks and Horozov2006; Cates & Clegg Reference Cates and Clegg2008). Closely related structures can also be made with air as one of the two fluids (Subramanian et al. Reference Subramaniam, Abkarian and Stone2005). Unless Janus particles are used (Aveyard Reference Aveyard2012), these emulsions are always metastable: their minimum free energy state comprises bulk A and B phases, separated by a flat interface, with as much of this interface as possible covered by particles, and the remainder distributed randomly in one or both solvents. However, the vast detachment energies allow the metastable state to survive almost indefinitely in many cases.

8.2 Resistance to coalescence and Ostwald ripening

One route to Pickering emulsions is to create numerous small droplets by applying a strong flow that mixes fluids A and B. The flow sweeps particles onto the interface whose initial surface area is much larger than they can cover. Coalescence initially proceeds as normal (perhaps assisted by maintaining a lower flow rate). Coalescence decreases the surface-to-volume ratio of droplets at fixed numbers of attached colloids. This proceeds until the surface particle density is high enough to prevent further coalescence. This requires a coverage comparable to, but somewhat below, that of a densely packed 2D amorphous film. Hence the resulting particle layer is usually not jammed, and droplets can relax to a spherical shape. In some cases, though, for instance if droplets remained stretched by flow as they coalesce, the final droplet is arrested in an aspherical jammed state, with particles clamped in position by interfacial tension (Cates & Clegg Reference Cates and Clegg2008).

For both spherical and aspherical structures, these clamped layers offer very strong stability against coalescence. They also resist Ostwald ripening, for the following reason. Recall that a close-packed monolayer of particles (whether ordered or amorphous) can be placed on the surface of a spherical fluid droplet of radius $R$ (say) without altering its interfacial geometry. Imagine such a droplet with the particles just in contact with one another. There is still a fluid–fluid interface at the interstices between particles, and this has the same curvature, and hence Laplace pressure $\unicode[STIX]{x1D6FE}_{0}/R$ , as the original drop.

Suppose now that this droplet is in diffusive equilibrium with one or more larger ones. According to the Ostwald mechanism, it will start to shrink. However, if the particles are already in contact, they cannot follow the droplet surface inwards as this happens. Moreover, each particle demands an unchanged contact angle with the interface, which is effectively now pinned to the particle layer. It is easy to see (figure 5) that even a small loss of volume of the droplet under these conditions will reduce the Laplace pressure to zero and then negative values. Once zero is reached, the droplet is fully resistant to Ostwald ripening. The mechanism is similar to that described in § 5.5 using a trapped species; indeed, due to their high detachment energies, interfacial particles are effectively such a species.

8.3 Morphologies of particle-stabilized emulsions

Figure 6.

Various particle-stabilized interfacial structures between fluids A and B, as imaged by confocal microscopy (which effectively views a thin slice through the material). (a) A biliquid foam with thin B films separating polyhedral A droplets; the bright regions are layers of fluorescently labelled particles and appear as lines or regions according to whether these lie oblique to the confocal plane or within it. (b) A multiple emulsion. The brightest regions (light green) are again particle-rich. Fluid A is dark grey and fluid B is dull red. (c) A bijel (see text for meaning); the fluorescent labelling is similar to the multiple emulsion. (Images courtesy of Paul Clegg.)

Alongside conventional spherical emulsion droplets, with particle stabilization a number of alternative morphologies are available (figure 6). First, drainage or centrifugation of a Pickering emulsion leads to a compressed foam structure that is like the biliquid foams discussed previously. These can be very stable thanks to the combined resistance to coalescence and Ostwald ripening provided by the particles. Second, simple manual agitation of binary immiscible solvents containing partially wettable particles often results in droplet-within-droplet structures known as multiple emulsions. Such structures require stability against both coalescence and Ostwald ripening for A droplets in B and for B droplets in A simultaneously. This is relatively difficult for surfactant formulations but seemingly quite easy with particle-stabilized ones (Clegg et al. Reference Clegg, Tavacoli and Wilde2016).

A third interesting morphology is that of ‘bijels’, or bicontinuous interfacially jammed emulsion gels, which are metastable analogues of the bicontinuous microemulsion: a particle layer resides at the interface between interpenetrating domains of A and B. This structure was predicted first computationally by model H simulations with added colloids (Stratford et al. Reference Stratford, Adhikari, Pagonabarraga, Desplat and Cates2005); it was confirmed in the laboratory by Herzig et al. (Reference Herzig, White, Schofield, Poon and Clegg2007). In a bijel, the interfacial film of non-detachable particles is clamped by tension into a 2D jammed layer, which imparts solidity to the whole 3D structure. This robustness can be improved further by having an interaction potential between particles with a steep barrier and then an attractive minimum at short distances. The interfacial tension pushes particles over the barrier, creating a strongly bonded interfacial film (Sanz et al. Reference Sanz, White, Clegg and Cates2009).

One recipe for making a bijel is to choose a fluid pair A + B that are miscible at high temperature (Cates & Clegg Reference Cates and Clegg2008). The colloidal particles are then dispersed within the single-phase mixture. On dropping the temperature, the fluids separate and particles are swept onto the interface. The coarsening process arrests when a jammed monolayer is formed, creating the bijel. The final structural domain size obeys $L\simeq a/\unicode[STIX]{x1D719}_{p}$ , with $a$ and $\unicode[STIX]{x1D719}_{p}$ the particle size and volume fraction, and the elastic modulus of the solid bijel scales as $G\sim \unicode[STIX]{x1D6FE}_{0}/L$ . Bijels are currently being explored for various applications in materials design (Lee et al. Reference Lee, Thijssen, Witt and Clegg2013).

9.1 Free energy of liquid-crystalline emulsions

In a polar liquid-crystalline emulsion, the hydrodynamic variables are the composition variable $\unicode[STIX]{x1D719}(\boldsymbol{r},t)$ , the polarization $\boldsymbol{p}(\boldsymbol{r},t)$ and the fluid velocity $\boldsymbol{v}(\boldsymbol{r},t)$ . These can be defined such that $\unicode[STIX]{x1D719}\simeq 1$ in the bulk liquid crystal where $|\boldsymbol{p}|>0$ , and $\unicode[STIX]{x1D719}\simeq -1$ in the bulk isotropic fluid where $|\boldsymbol{p}|=0$ . The free energy of such a system can be written as a sum of two contributions: $F[\unicode[STIX]{x1D719},\boldsymbol{p}]=F_{\unicode[STIX]{x1D719}}[\unicode[STIX]{x1D719}]+F_{\boldsymbol{p}}[\unicode[STIX]{x1D719},\boldsymbol{p}]$ . The first is the contribution from a simple binary fluid, similar to (3.2),

Here $a>0$ and $\unicode[STIX]{x1D705}>0$ are chosen such that $F$ is minimized by bulk phase separation into states with $\unicode[STIX]{x1D719}\simeq \pm 1$ . The second contribution to the free energy stems from the liquid crystallinity, which can be written as (de Gennes & Prost Reference de Gennes and Prost2002; Tjhung et al. Reference Tjhung, Marenduzzo and Cates2012)

Here $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D719})$ is a thermodynamic parameter that controls the isotropic to polar transition: $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D719})<0$ in the polar phase and $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D719})>0$ in the isotropic phase. (This notation is conventional and should be distinguishable by context from our previous use of $\unicode[STIX]{x1D6FE}$ to denote interfacial tension.) For simplicity we can assume $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D719})=-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D719}$ to obtain a bulk polar phase whenever $\unicode[STIX]{x1D719}>0$ and an isotropic phase ( $|\boldsymbol{p}|=0$ ) when $\unicode[STIX]{x1D719}<0$ . With this choice, $|\boldsymbol{p}|=1$ when $\unicode[STIX]{x1D719}=1$ , in effect setting the units for $\boldsymbol{p}$ . We require the quartic coefficient $\unicode[STIX]{x1D6FC}$ to be positive for thermodynamic stability. The gradient terms involve $K_{1}$ , $K_{2}$ and $K_{3}$ , which are the splay, twist and bend elastic constants (all positive); see figure 7(a). In the literature, one often finds the approximation $K_{1}=K_{2}=K_{3}\equiv K$ , such that all three elastic terms combine into the simpler form $(K/2)(\unicode[STIX]{x2202}_{i}p_{j})^{2}$ . The remaining terms in (9.2) describe anchoring effects. The $\unicode[STIX]{x1D6FD}_{1}$ term favours perpendicular anchoring of $\boldsymbol{p}$ at the droplet interface (see figure 7 b): if $\unicode[STIX]{x1D6FD}_{1}>0$ , $\boldsymbol{p}$ tends to point outwards (from the polar to the isotropic phase), whereas if $\unicode[STIX]{x1D6FD}_{1}<0$ , $\boldsymbol{p}$ tends to point inwards. Finally, the $\unicode[STIX]{x1D6FD}_{2}$ term promotes parallel anchoring of $\boldsymbol{p}$ at the interface (see figure 7 c).

For the case of nematic (rather than polar) emulsions, the second-rank tensor $\unicode[STIX]{x1D64C}(\boldsymbol{r},t)$ takes the place of $\boldsymbol{p}(\boldsymbol{r},t)$ as a hydrodynamic variable, alongside the composition $\unicode[STIX]{x1D719}(\boldsymbol{r},t)$ and the fluid velocity $\boldsymbol{v}(\boldsymbol{r},t)$ . The free energy can be written as $F[\unicode[STIX]{x1D719},\unicode[STIX]{x1D64C}]=F_{\unicode[STIX]{x1D719}}[\unicode[STIX]{x1D719}]+F_{\unicode[STIX]{x1D64C}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D64C}]$ , where $F_{\unicode[STIX]{x1D719}}$ is of the same form as (9.1). Here $F_{\unicode[STIX]{x1D64C}}$ comprises the so-called Landau–de Gennes free energy functional for bulk nematics (de Gennes & Prost Reference de Gennes and Prost2002), augmented with appropriate couplings to $\unicode[STIX]{x1D719}$ . In three dimensions this reads (Sulaiman et al. Reference Sulaiman, Marenduzzo and Yeomans2006)

where for simplicity we have taken a single elastic constant $K>0$ . The thermodynamic parameter $A_{0}>0$ is a scale factor for bulk free energies, while $\unicode[STIX]{x1D6FE}^{\prime }(\unicode[STIX]{x1D719})$ controls the isotropic to nematic transition; the nematic phase is stable for $\unicode[STIX]{x1D6FE}^{\prime }>2.7$ and the isotropic phase is stable for $\unicode[STIX]{x1D6FE}^{\prime }<2.7$ . This transition is discontinuous (thanks to the term in (9.3) cubic in $\unicode[STIX]{x1D64C}$ ) with a finite hysteresis window within which the phase of higher $F_{\unicode[STIX]{x1D64C}}$ remains metastable. (These statements can be checked by parametrizing $\unicode[STIX]{x1D64C}=\unicode[STIX]{x1D706}(\hat{z}\hat{z}-\unicode[STIX]{x1D644}/3)$ and then minimizing $F_{\unicode[STIX]{x1D64C}}(\unicode[STIX]{x1D706})$ with respect to $\unicode[STIX]{x1D706}$ .) Choosing $\unicode[STIX]{x1D6FE}^{\prime }(\unicode[STIX]{x1D719})=2.7+\unicode[STIX]{x1D6FE}_{1}\unicode[STIX]{x1D719}$ with some positive constant $\unicode[STIX]{x1D6FE}_{1}$ ensures a nematic phase at compositions $\unicode[STIX]{x1D719}>0$ and an isotropic phase at $\unicode[STIX]{x1D719}<0$ . For nematics, the physics of orientational anchoring at the surface of a droplet can be captured by a single term, proportional to $\unicode[STIX]{x1D6FD}_{0}$ . For $\unicode[STIX]{x1D6FD}_{0}<0$ , the director field $\hat{\boldsymbol{n}}$ will tend to align perpendicular to the interface, whereas for $\unicode[STIX]{x1D6FD}_{0}>0$ , $\hat{\boldsymbol{n}}$ will tend to align parallel to the interface. This is similar to figures 7(b) and 7(c), respectively, except that one must replace the vector with a headless vector in each case.

9.2 Topological defects

If we ignore the anchoring term $\unicode[STIX]{x1D6FD}_{1}$ in the polar free energy $F_{\boldsymbol{p}}$ , we notice that the free energy is invariant under a global inversion $\boldsymbol{p}\rightarrow -\boldsymbol{p}$ . On the other hand, for nematics, the free energy is invariant under a local inversion $\hat{\boldsymbol{n}}\rightarrow -\hat{\boldsymbol{n}}$ . As long as we are in a region without topological defects – which are places where $\boldsymbol{p}$ or $\hat{\boldsymbol{n}}$ are undefined – there should be no distinction between a global or local symmetry since any attempt to make different sign choices in different neighbourhoods will lead to $\boldsymbol{p}$ or $\hat{\boldsymbol{n}}$ being undefined wherever these neighbourhoods meet. In practice, therefore, in defect-free regions the static and dynamic properties for both $\boldsymbol{p}$ and $\hat{\boldsymbol{n}}$ are found to be very similar, so long as we ignore terms that break the $\boldsymbol{p}\leftrightarrow -\boldsymbol{p}$ symmetry in $F_{\boldsymbol{p}}$ , such as the $\unicode[STIX]{x1D6FD}_{1}$ term in (9.2).

However, differences arise when we do have topological defects. Such a defect is present when, if the state of order in some region is mapped onto a state of uniform orientation, this map cannot be made smooth everywhere in real space (Chaikin & Lubensky Reference Chaikin and Lubensky1995). In two dimensions, the lowest-order topological defects for nematics ( $\hat{\boldsymbol{n}}$ ) are point-like and have ‘topological charge’ $+1/2$ or $-1/2$ , as shown in figure 8(a). This charge denotes the number of full turns of the director along a path that makes one full turn around the defect. Because the vector $\hat{\boldsymbol{n}}$ is headless, a rotation comprising a whole number of half-turns in either direction brings it back to the same state. (Defects of charge $\pm n/2$ with $n>1$ are also possible, but in practice these rapidly dissociate for energetic reasons into $n$ half-integer defects of the appropriate sign.) On the other hand, the lowest-order defects for polar fluids have charge $\pm 1$ since a full rotation of $\boldsymbol{p}$ is needed to recover the same state, as shown in figure 8(b). (Again, defects of charge $\pm n$ with $n>1$ dissociate to reduce the elastic energy.) In both the polar and nematic cases, defects of the same sign repel one another whereas those of opposite charge attract and then annihilate. Conversely, during the isotropic to liquid crystal transition, triggered typically by a quench which alters $\unicode[STIX]{x1D6FE}$ or $\unicode[STIX]{x1D6FE}^{\prime }$ , pairs of defects with opposite signs are created.

Figure 8.

(a) Topological defects in 2D nematic liquid crystals. (b) In polar liquid crystals half-integer defects are forbidden. (c) In three dimensions, $+1/2$ and $-1/2$ point defects of the nematics become line defects. Furthermore, $+1/2$ and $-1/2$ line defects are then topologically equivalent (one can see this by flipping each molecule $180^{\circ }$ about the $y$ axis) (Chaikin & Lubensky Reference Chaikin and Lubensky1995).

In three dimensions, the $\pm 1/2$ point defects of the nematics become line defects, as shown in figure 8(c). Furthermore, in three dimensions, $+1/2$ and $-1/2$ line defects are topologically equivalent. To see this, consider the $+1/2$ line defect shown in figure 8(c), left. We can then continuously rotate the director field everywhere through $180^{\circ }$ about an axis perpendicular to $\hat{\boldsymbol{z}}$ (in this example, the $y$ axis) to get the $-1/2$ line defect shown in figure 8(c), right. Similar reasoning establishes that any two defects whose charge differs by an integer are equivalent, so that all integer defect lines are equivalent to no defect at all, and all half-integer defects are equivalent to each other. Thus in 3D nematics there is only one type of line defect; these are called ‘disclinations’. Any two disclination lines can annihilate. Note that disclinations can also form closed loops (Chaikin & Lubensky Reference Chaikin and Lubensky1995).

For polar liquid crystals there are no line defects in three dimensions. This is because cylindrical versions of the structures shown in figure 8(b) can be converted into a defect-free state by smoothly rotating all the arrows to point out of (or into) the page. The basic defects are instead the obvious 3D equivalents of the $\pm 1$ point defects shown in figure 8(b), left and right. These are called the radial hedgehog and hyperbolic hedgehog defects, respectively (Lubensky et al. Reference Lubensky, Pettey, Currier and Stark1998). These two types of point defects are also possible in nematics, but are then topologically equivalent to one another, and also equivalent, at distances much larger than its radius, to a closed disclination loop. Note that their topological equivalence does not mean that two structures are freely interconvertible; energetic considerations often favour one over the other.

We also note in passing that, at the centre of a line defect in a polar liquid crystal, $\boldsymbol{p}$ passes through zero such that its direction is undefined. For a 2D nematic, $\unicode[STIX]{x1D64C}$ likewise becomes zero at the defect core, meaning that the medium is isotropic there and $\hat{\boldsymbol{n}}$ is undefined. However, at the core of a disclination line in three dimensions, $\hat{\boldsymbol{n}}$ is undefined not because $\unicode[STIX]{x1D64C}$ is zero, but because its two largest eigenvalues are positive and degenerate, so that the orientational distribution of molecules is isotropic in the plane normal to the defect line.

Figure 9.

(a) Long chains of isotropic droplets in a bulk nematic fluid; the chains have roughly equal separation. (b) A hyperbolic hedgehog (point defect) formed by the nematic liquid crystal just outside an isotropic droplet. (c) A half-integer ring defect can also form around an isotropic droplet (saturn ring). (d) Boojum defects at the poles of a nematic droplet. (e) Radial hedgehog in the middle of a nematic droplet. (Panel (a) is adapted from Loudet et al. (Reference Loudet, Barois and Poulin2000) and (b,c) are adapted from Lubensky et al. (Reference Lubensky, Pettey, Currier and Stark1998).)

9.3 Defects in and around emulsion droplets

Emulsions comprising droplets of isotropic fluid in a bulk nematic (sometimes called ‘inverted’ nematic emulsions) can display interesting states of organization such as the one shown in figure 9 (Loudet et al. Reference Loudet, Barois and Poulin2000). Here, isotropic droplets can form parallel and very long chains, separated by roughly equal distance. This phase is obtained by quenching the mixture from an initial high-temperature phase that is uniform and isotropic, so that it phase-separates into isotropic and nematic phases along the lines discussed in § 5 for binary mixtures of simple fluids.

These long chains of isotropic droplets are stabilized in three dimensions by topological defects. In this system, surface anchoring requires $\hat{\boldsymbol{n}}$ to be radial at the droplet edge, creating in effect a radial hedgehog that cannot be annihilated because its core is effectively inside the droplet. To restore uniform ordering at large scales, there needs to be either another point defect just outside the droplet, for which the lowest-energy choice is a hyperbolic hedgehog (see figure 9 b), or a disclination loop (which is equivalent at large distances as discussed above) (Lubensky et al. Reference Lubensky, Pettey, Currier and Stark1998). This loop can be off-centred or equatorial, as shown in figure 9(c). For the asymmetric arrangements, such as that with the hyperbolic hedgehog, long-ranged elastic distortion in the director field creates an effective dipolar attraction (Poulin et al. Reference Poulin, Stark, Lubensky and Weitz1997) and this causes the droplets to come together to form chains (figure 9 a, inset). The interaction is quadrupolar for a symmetric (equatorial) disclination loop, with attraction only when the separation between particles is in a range of oblique angles to the far-field director. If perpendicular anchoring is replaced by parallel anchoring, the exterior nematic instead develops a pair of $+1$ defects at opposite poles of the included droplet. Note that half of each defect, which would lie interior to the droplet, is effectively absent. Such surface defects are called ‘boojums’. The leading-order interaction is again quadrupolar (Poulin et al. Reference Poulin, Stark, Lubensky and Weitz1997).

Interesting structures also arise in ‘direct’ nematic emulsions comprising nematic droplets immersed in an isotropic fluid. Again, the director field can be aligned either parallel to the droplet interface ( $\unicode[STIX]{x1D6FD}_{0}>0$ in the free energy equation (9.3)) or perpendicular to it ( $\unicode[STIX]{x1D6FD}_{0}<0$ ). In the case of parallel anchoring, the director field typically forms a pair of boojum defects at each pole of the droplet (now with the missing part of the defect outside rather than inside the droplet; see figure 9 d). These point defects create a cusp at each pole of the droplet (Prinsen & van der Schoot Reference Prinsen and van der Schoot2003). On the other hand, in the case of perpendicular anchoring, the director field will typically form a single radial hedgehog in the centre of the droplet (Lopez-Leon & Fernandez-Nieves Reference Lopez-Leon and Fernandez-Nieves2011), at least within the single elastic constant approximation (figure 9 e). Importantly, because nematic order exists only in the droplet interior, these effects do not lead to any long-range interaction between droplets through the isotropic continuous phase. For this reason, the self-assembly properties of direct nematic emulsions are broadly the same as for conventional emulsions of two isotropic liquids. They can be functionally useful, however, because the discrete nematic compartments can have faster switching times in response to external fields than a continuous phase in which defects can move slowly across large distances. Materials where the continuous phase is a polymer solution, from which the solvent is then evaporated to give nematic droplets in a solid matrix, are referred to as polymer-dispersed liquid crystals and are widely used in technology (Bouteiller & LeBarney Reference Bouteiller and Lebarny1996). A related structure is obtained when a binary fluid is quenched into isotropic $+$ nematic coexistence in the presence of colloidal particles. To minimize anchoring constraints, these tend to segregate into the isotropic phase. When the final phase volume of that phase is small, a biliquid foam resembling that in figure 6 is formed, with particles compressed into films that now lie between nematic cells with different directors (Anderson et al. Reference Anderson, Terentjev, Meeker, Crain and Poon2001).

10.1 Equations of motion: polar liquid crystals

Consider a patch of polar liquid-crystalline material as shown in figure 10(a). Let us initially assume that this patch rotates as a rigid body with some angular velocity $\unicode[STIX]{x1D74E}$ . In other words, the fluid velocity $\boldsymbol{v}(\boldsymbol{r},t)$ at position $\boldsymbol{r}$ at time $t$ is given by $\boldsymbol{v}=\unicode[STIX]{x1D74E}\times \boldsymbol{r}$ . At time $\unicode[STIX]{x1D6FF}t$ later, the $\boldsymbol{p}$ field can then be written as

The first term in this equation is simple advection; we displaced the material by $\boldsymbol{v}\unicode[STIX]{x1D6FF}t$ in the time interval $\unicode[STIX]{x1D6FF}t$ . However, as we can see from figure 10(a), the advective term alone is not enough; we also have to rotate $\boldsymbol{p}$ . This is given by the second term. Expanding to first order in $\unicode[STIX]{x1D6FF}t$ we obtain for the rigid rotation of $\boldsymbol{p}$ :

The second equality above follows from the fact that the angular velocity can be expressed as $\unicode[STIX]{x1D714}_{i}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D716}_{ijk}\unicode[STIX]{x1D6FA}_{jk}$ , where $\unicode[STIX]{x1D6FA}_{ij}\equiv {\textstyle \frac{1}{2}}(\unicode[STIX]{x2202}_{i}v_{j}-\unicode[STIX]{x2202}_{j}v_{i})$ is the antisymmetric part of the velocity gradient tensor.

Figure 10.

(a) Rigid-body rotation in a patch of liquid-crystalline material; (red) arrows indicate $\boldsymbol{p}$ . (b) However, most liquid crystals do not rotate like a rigid body; they tend to align with shear flow at steady state. Here $\dot{\unicode[STIX]{x1D6FE}}$ is the strain rate and $\unicode[STIX]{x1D703}_{L}$ is the Leslie angle, which is related to the phenomenological parameter $\unicode[STIX]{x1D709}$ . Horizontal (blue) arrows (left) represents fluid velocity $\boldsymbol{v}(\boldsymbol{r})$ and inclined (red) arrows (right) represent the orientational field $\boldsymbol{p}(\boldsymbol{r})$ at steady state.

In general flows, such as shear flows, liquid-crystalline materials do not rotate like a rigid body; typically $\boldsymbol{p}$ tends to align with the streamlines of the flow. This results in an additional contribution in the full equation of motion for $\boldsymbol{p}$ ; since it is advective, this term is bilinear in $\unicode[STIX]{x1D735}\boldsymbol{v}$ and $\boldsymbol{p}$ . The resulting full equation reads

Here we have introduced $\unicode[STIX]{x1D709}$ as the flow alignment parameter in the second term on the right. In this term, $\unicode[STIX]{x1D60B}_{ij}\equiv {\textstyle \frac{1}{2}}(\unicode[STIX]{x2202}_{i}v_{j}+\unicode[STIX]{x2202}_{j}v_{i})$ is the symmetric part of the velocity gradient tensor. The parameter $\unicode[STIX]{x1D709}$ is set by the molecular geometry and is independent of the terms in the free energy $F$ . If $|\unicode[STIX]{x1D709}|>1$ , the orientation of the liquid crystal tends to align in shear with the flow direction, with a steady-state orientation set by the so-called Leslie angle, $\unicode[STIX]{x1D703}_{L}=\tan ^{-1}[(\unicode[STIX]{x1D709}-1)/(\unicode[STIX]{x1D709}+1)]^{1/2}$ (see figure 10 b). On the other hand, if $|\unicode[STIX]{x1D709}|<1$ , the orientation $\boldsymbol{p}$ is not stationary in a shear flow but instead shows tumbling behaviour (Larson Reference Larson1999). Importantly, there is no corresponding molecular influence on the coefficient of the $\unicode[STIX]{x1D734}$ term: this is always unity, otherwise $\boldsymbol{p}$ would fail to evolve properly when the only motion is the slow rigid rotation of the entire sample (Beris & Edwards Reference Beris and Edwards1994).

The final ingredient in the dynamics for $\boldsymbol{p}$ is the relaxation term proportional to $1/\unicode[STIX]{x1D6E4}$ in (10.3). Here $\boldsymbol{h}(\boldsymbol{r},t)=\unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\boldsymbol{p}$ is called the molecular field, and is analogous to the chemical potential $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}$ for $\unicode[STIX]{x1D719}$ . This term, which relaxes $\boldsymbol{p}$ towards the minimum of the free energy, differs in form from the corresponding term in the equation of motion for the composition $\unicode[STIX]{x1D719}$ since $\boldsymbol{p}$ is not a locally conserved quantity. The equation for $\unicode[STIX]{x1D719}$ itself is the same as in model H for a simple binary fluid (see (4.5) above), as is the NSE (see (4.3)) for $\boldsymbol{v}$ except that a different expression, detailed below, is now needed for the thermodynamic stress  $\unicode[STIX]{x1D748}$ .

To summarize, for polar liquid crystals, at hydrodynamic (deterministic) level, the dynamics is governed by the following set of equations (Kung et al. Reference Kung, Marchetti and Saunders2006):

Here $M$ is the mobility and $\unicode[STIX]{x1D707}[\unicode[STIX]{x1D719},\boldsymbol{p}]=\unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}$ is the chemical potential. By construction, these reduce to the noise-free version of model H in the limit where $\boldsymbol{p}=\mathbf{0}$ everywhere.

It only remains to derive the total elastic stress $\unicode[STIX]{x1D748}[\unicode[STIX]{x1D719},\boldsymbol{p}]$ which enters the NSE (10.4). We know already from § 3.3 that, even in the absence of $\boldsymbol{p}$ , there is an interfacial stress contribution. This takes the form (compare (3.10))

where $\mathbb{F}_{\unicode[STIX]{x1D719}}=f(\unicode[STIX]{x1D719})+\unicode[STIX]{x1D705}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})^{2}/2$ is the local free energy density as defined in (9.1). In the final term we have generalized (3.10) to allow for the addition of anchoring terms, which also can contribute to the chemical potential and the interfacial stress; hence the full free energy density $\mathbb{F}=\mathbb{F}_{\unicode[STIX]{x1D719}}+\mathbb{F}_{\boldsymbol{p}}$ appears here. This can be confirmed by a straightforward analogue of the following argument, in which for simplicity we calculate the elastic stress contribution from $\boldsymbol{p}$ only (Beris & Edwards Reference Beris and Edwards1994).

To find this stress contribution, we consider the rate of change in $F_{\boldsymbol{p}}$ :

Substituting (10.7) for ${\dot{p}}_{i}$ , we obtain

where the underlined term has been integrated by parts. The resulting surface contribution vanishes for periodic boundary conditions, which we choose here without loss of generality (just as we did in § 3.3 when discussing model H). The underlined term in (10.10) is reminiscent of (3.8) for $\unicode[STIX]{x1D6FF}F(\unicode[STIX]{x1D719})$ in the case where only a composition variable is present. Thus $-p_{j}\unicode[STIX]{x2202}_{i}h_{j}$ is effectively a force density, playing a similar role to $-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}_{i}\unicode[STIX]{x1D707}$ there. Allowing for the fact that the small displacement $\boldsymbol{u}$ is now $\boldsymbol{v}\text{d}t$ , this term gives a contribution $\unicode[STIX]{x1D6FF}F_{\boldsymbol{p}}^{\prime }=\int (p_{j}\unicode[STIX]{x2202}_{i}h_{j})u_{i}\,\text{d}\boldsymbol{r}$ (where the prime on $F_{\boldsymbol{p}}^{\prime }$ reminds us that this is not the only contribution). Following the same procedure using (3.9) as in § 3.3, this translates into an elastic stress contribution that is the direct analogue of (10.8):

It is simple to check from this that $\unicode[STIX]{x2202}_{j}\unicode[STIX]{x1D70E}_{ij}^{\prime }=-p_{j}\unicode[STIX]{x2202}_{i}h_{j}$ .

The remaining contribution $\unicode[STIX]{x1D70E}_{ij}^{\prime \prime }$ to the elastic stress has no counterpart in model H; it stems from the rotation and alignment terms (10.10), which may be written

In the second equality, we have decomposed $h_{i}p_{j}$ into symmetric and antisymmetric parts, and used the fact that $\unicode[STIX]{x1D61F}_{ij}\unicode[STIX]{x1D620}_{ij}=0$ for any symmetric tensor $\unicode[STIX]{x1D61F}_{ij}$ and antisymmetric tensor $\unicode[STIX]{x1D620}_{ij}$ . Once again using (3.9) with $\boldsymbol{u}=\boldsymbol{v}\,\text{d}t$ we obtain

where the second form decomposes the rate-of-strain tensor $\unicode[STIX]{x2202}_{j}v_{i}$ into its symmetric and antisymmetric parts. Comparing (10.12) and (10.13), we identify ${\unicode[STIX]{x1D70E}_{ij}^{\prime \prime }}^{S}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D709}(p_{i}h_{j}+p_{j}h_{i})$ and ${\unicode[STIX]{x1D70E}_{ij}^{\prime \prime }}^{A}=-{\textstyle \frac{1}{2}}(p_{i}h_{j}-p_{j}h_{i})$ . The presence of an antisymmetric contribution is generic in liquid crystals in which rotational symmetry is spontaneously broken.

Note finally that the last term, proportional to $1/\unicode[STIX]{x1D6E4}$ , in (10.10) causes a purely dissipative loss of free energy and hence does not contribute to the elastic stress, which instead describes how the stored free energy changes with sample shape. This term is negative definite, ensuring that, in the absence of a driving force, the free energy steadily decreases towards its minimum value.

Adding the contributions $\unicode[STIX]{x1D748}^{\unicode[STIX]{x1D719}}$ and $\unicode[STIX]{x1D748}^{\boldsymbol{p}}=\unicode[STIX]{x1D748}^{\prime }+\unicode[STIX]{x1D748}^{\prime \prime S}+\unicode[STIX]{x1D748}^{\prime \prime A}$ found above, we can assemble the total elastic stress tensor $\unicode[STIX]{x1D748}[\unicode[STIX]{x1D719},\boldsymbol{p}]$ in the NSE (10.4). This reads

It is important to note that, alongside the rigid-rotation term (the final term on the right), the flow alignment parameter $\unicode[STIX]{x1D709}$ , as defined via the equation of motion (10.7) for $\boldsymbol{p}$ , also enters the elastic stress in a non-negotiable fashion (Beris & Edwards Reference Beris and Edwards1994). This is because it controls the response of $\boldsymbol{p}$ to an incremental strain and hence affects the resulting free energy increment which determines the stress.

10.2 Equations of motion: nematic liquid crystals

In principle, one can repeat the same calculation as above for nematic liquid crystals. Here we shall just quote the results (Beris & Edwards Reference Beris and Edwards1994), which are included for completeness. The dynamical equation for the order parameter $\unicode[STIX]{x1D64C}$ is

where the first term on the right describes the combined effect of rigid-body rotation and shear aligning. This takes the form

where $\unicode[STIX]{x1D709}$ is a flow aligning parameter analogous to that introduced previously for $\boldsymbol{p}$ . This equation describes similar physics to the corresponding terms for $\boldsymbol{p}$ . The quantity $\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D644}/d$ enters (rather than just the traceless $\unicode[STIX]{x1D64C}$ ) because the entire orientational distribution function, including its isotropic part, is stretched and rotated by the flow. The final term in (10.15) is a local relaxation term similar to that used in (10.7) for $\boldsymbol{p}$ , with $\unicode[STIX]{x1D6E4}$ an analogous relaxation constant. The molecular field for nematics is defined as

This is traceless by construction, so as to maintain $\text{Tr}\,\unicode[STIX]{x1D64C}=0$ under time evolution. The equations of motion for $\unicode[STIX]{x1D719}$ and $\boldsymbol{v}$ are the same as before, i.e. (10.4), (10.5) and (10.6), with the total elastic stress $\unicode[STIX]{x1D748}(\unicode[STIX]{x1D64C},\unicode[STIX]{x1D719})$ now given by

This equation set was, for example, solved numerically by Sulaiman et al. (Reference Sulaiman, Marenduzzo and Yeomans2006), who addressed droplet shapes and defect textures in equilibrium, and also switching behaviour in an applied external field $\boldsymbol{E}(t)$ (which adds a term in $\boldsymbol{E}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\boldsymbol{E}$ to the free energy density $\mathbb{F}$ , not needed here). Note, however, that for both the nematic case and the polar one of the preceding section, only recently have the relevant numerical tools (primarily involving the lattice Boltzmann method) been developed to solve the equations of motion presented above (Sulaiman et al. Reference Sulaiman, Marenduzzo and Yeomans2006; Cates et al. Reference Cates, Henrich, Marenduzzo and Stratford2009). With these methods in hand, it should be possible to understand more fully the unusual dynamical phenomena seen in direct and inverse liquid crystal emulsions, such as the kinetics of chain formation among isotropic emulsion droplets in a nematic fluid (Poulin et al. Reference Poulin, Stark, Lubensky and Weitz1997). Conversely, the recent numerical studies of rheology in such systems (Tiribocchi et al. Reference Tiribocchi, Da Re, Marenduzzo and Orlandini2016) may hopefully promote new experimental studies of their response to imposed flow.

11.1 Active models B and H

Continuum models for the description of motility-induced phase separation remain under development. Here we outline some of the interesting cases looked at so far. We start by suppressing the fluid velocity so that the relevant passive model is model B, described by (4.9) and (4.10). The simplest way to break time-reversal symmetry in this model is to retain (4.9) but add to the chemical potential in (4.10) a term that is not of the form $\unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}$ . To lowest order in gradients, this term is $|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|^{2}$ . We therefore introduce a non-equilibrium chemical potential

Slightly more generally, one can consider the case where $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D719})$ and $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D719})$ are functions of composition; time-reversal symmetry is then broken ( $\unicode[STIX]{x1D707}\neq \unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}$ for any $F$ ) whenever $\unicode[STIX]{x1D706}\neq \text{d}\unicode[STIX]{x1D705}/\text{d}\unicode[STIX]{x1D719}$ . Microscopic models (Stenhammar et al. Reference Stenhammar, Tiribocchi, Allen, Marenduzzo and Cates2013) give exactly this structure, albeit with non-polynomial $F$ and non-trivial $\unicode[STIX]{x1D719}$ dependence in $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D706}$ ; suppressing the latter dependence and restoring the simplest form for the local part of $\unicode[STIX]{x1D707}$ gives (11.1) as the prototypical model for motility-induced phase separation without fluid flow.

Active model B has some interesting properties (Wittkowski et al. Reference Wittkowski, Tiribocchi, Stenhammar, Allen, Marenduzzo and Cates2014). First, the non-TRS term (which also breaks $\unicode[STIX]{x1D719}\rightarrow -\unicode[STIX]{x1D719}$ symmetry) alters the phase boundaries at coexistence, even at mean-field (zero-noise) level. This was unexpected since in equilibrium problems the phase diagram is found by a common-tangent construction on $f(\unicode[STIX]{x1D719})$ (see (3.3)) in which no gradient terms arise. However, it turns out that this construction is valid only if the gradient terms stem from a free energy functional. If one continues to define $F[\unicode[STIX]{x1D719}]$ as the functional arising when the active term is switched off ( $\unicode[STIX]{x1D706}=0$ ), the effect of activity is to create an inequality between phases in the pressure-like quantity $P_{Th}\equiv \unicode[STIX]{x1D707}\unicode[STIX]{x1D719}-f(\unicode[STIX]{x1D719})$ , causing a shift of the binodal compositions away from their values at $\unicode[STIX]{x1D706}=0$ . (This pressure, which is defined by the usual equilibrium relation between $P,\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D719}$ , should not be confused with the mechanical one defined as the force density acting on a wall; only in passive systems are these two definitions equivalent (Solon et al. Reference Solon, Fily, Baskaran, Cates, Kafri, Kardar and Tailleur2015).) The resulting ‘pressure jump’ across the interface, which when small is linear in $\unicode[STIX]{x1D706}$ , effectively provides an anomalous active contribution to the Laplace pressure, which is finite even for a flat interface.

Because it interferes with Laplace pressure and hence with the driving force for coarsening dynamics, one might expect activity to have some important influence on the diffusive growth law (5.12) that gave a domain size $L\sim t^{1/3}$ . However, in active model B there is no conclusive numerical evidence for a change in exponent (Wittkowski et al. Reference Wittkowski, Tiribocchi, Stenhammar, Allen, Marenduzzo and Cates2014). Since coarsening is driven by Laplace pressure differences, this outcome can be rationalized by noting that the activity-induced pressure jump across interfaces is curvature-independent at leading order. For the same reason, diffusive coarsening continues indefinitely. The latter is also found to be true for the particular microscopic models (of so-called active Brownian particles, or ABPs) that inspired the form of (11.1). In several experimental systems, however, coarsening appears to saturate while clusters are at a finite size of perhaps 60–100 particles; the reasons for this are not clear, and various explanations have been suggested, for instance involving cross-phoresis mechanisms (Buttinoni et al. Reference Buttinoni, Bialke, Kummel, Lowen, Bechinger and Speck2013).

This suggests that active model B, though appealingly simple, may not capture all we need to address the phase behaviour of scalar active matter. That view is confirmed by the observation that the mathematical structure of (11.1), in combination with (4.9), enforces $\unicode[STIX]{x1D735}\times \boldsymbol{J}=\mathbf{0}$ . This rules out steady-state circulating particle currents in real space. These currents are a low-dimensional projection of a circulating probability flux in the space of configurations $\unicode[STIX]{x1D719}(\boldsymbol{r})$ . (The latter, more abstract, currents do remain present, however: one finds that in the phase-separated state there is continuous birth, in one phase, of droplets of the other, which then migrate to the interface and disappear; see Stenhammar et al. (Reference Stenhammar, Tiribocchi, Allen, Marenduzzo and Cates2013).) However, we know of situations where steadily circulating real-space currents do arise, at least in computer simulations – for example, when active Brownian particles are placed in a ratchet-like environment (Stenhammar et al. Reference Stenhammar, Wittkowski, Marenduzzo and Cates2016). Alongside the fact that active model B cannot explain cluster phases, this observation has motivated the recent introduction of an extended model, known as active model B+, in which additional time-reversal symmetry-breaking gradient terms are included in the expression for $\boldsymbol{J}$ (Nardini et al. Reference Nardini, Fodor, Tjhung, van Wijland, Tailleur and Cates2017). This model remains under investigation.

For systems in which the fluid velocity field $\boldsymbol{v}$ plays an important role, the natural starting point is model H, to which we can again add minimal TRS-breaking terms. One such term, in the chemical potential, is the same as just described for active model B. (The additional terms arising in active model B+ are yet to be addressed in this context.) Another new term enters the NSE (4.3) whose passive version contains a thermodynamic stress obeying $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}=-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}$ . This form assumes a thermodynamic relation between stress and chemical potential which only holds for equilibrium systems (in which mechanical forces and thermodynamic ones stem from the same microscopic Hamiltonian and are not independent). But in a system undergoing motility-induced phase separation, for instance, even the fact that the local ‘free energy density’ $f(\unicode[STIX]{x1D719})$ has two minima can arise purely from activity and not from attractive interactions. This means that, while the active contributions to $f(\unicode[STIX]{x1D719})$ do not break time-reversal symmetry in themselves, they have no reason to feed through via thermodynamics into the stress term in the NSE.

What matters in an incompressible fluid is the deviatoric stress, which is traceless and differs from the full stress by a pure pressure. From (3.10), this is (in $d$ dimensions)

in which $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D705}$ . In the absence of an external field that breaks rotational invariance, this form is in fact the only one possible to leading order in gradients, so that, in passing from the passive to the active case, all that is lost is the connection between $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D705}$ . Active model H thus reads (Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015)

where $\unicode[STIX]{x1D748}^{D}$ obeys (11.2), in which $\unicode[STIX]{x1D701}$ is a parameter that depends on both interaction forces and activity. Because this is no longer linked to $\unicode[STIX]{x1D705}$ (which is always postive), $\unicode[STIX]{x1D701}$ can have either sign. It includes an active contribution that is positive for ‘extensile’ swimmers, which draw fluid in along the axis of motion and expel it equatorially, and negative for ‘contractile’ ones, which do the opposite. (For a further discussion of this distinction, see § 12 and figure 12 below.)

Whenever $\unicode[STIX]{x1D701}\neq \unicode[STIX]{x1D705}$ , there are effectively two interfacial tensions in active model H, one controlling the diffusive flux (set by $\unicode[STIX]{x1D705}$ ) and one controlling the fluid flow driven by interfaces (now set by $\unicode[STIX]{x1D701}$ ). This mismatch breaks time-reversal symmetry, but if both tensions are positive, the consequences are relatively mild. (This statement is based on numerical studies in the simplest case which has $\unicode[STIX]{x1D701}\neq \unicode[STIX]{x1D705}$ but $\unicode[STIX]{x1D706}=0$ .) In particular, the physics of coarsening in the viscous hydrodynamic regime remains broadly consistent with linear scaling, $L\sim t$ . In the case when $\unicode[STIX]{x1D705}$ is finite but $\unicode[STIX]{x1D701}$ is zero, the order parameter field does not drive fluid motion and one recovers $L\sim t^{1/3}$ as for model B. On the other hand, decreasing $\unicode[STIX]{x1D701}$ below zero, as would be required to describe strongly contractile swimmers, one has in effect a negative interfacial tension in the mechanical sector (while that in the diffusive dynamics remains positive). Microscopically this arises because swimming particles tend to orient perpendicular to the interface between phases (there is a net polarization there proportional to $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ ), where their contractile swimming action pulls fluid inwards normal to the interface and pushes out sideways in the interfacial (equatorial) plane, causing the interface to stretch (figure 11 a). The spontaneous stretching motion of the interface is mechanically equivalent to a negative tension. Clearly, in this case one expects new and interesting effects to arise. One such effect is that the phase separation can arrest at a finite length scale where the diffusive shrinkage of the interfacial area is in balance with its contractile stretching (figure 11 b). This offers a hydrodynamic, rather than microscopic, mechanism for the existence of cluster phases, but only in cases where the swimming is contractile (Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015); its relation to earlier work in which the swimming particles are individually resolved (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008) is not yet established. The true character of cluster phases and their origins remains an active topic of current research (Saha et al. Reference Saha, Golestanian and Ramasawmy2014).

Figure 11.

Active model H describes a fluid of active particles with conserved fluid momentum, but no orientational interactions, so that a scalar order parameter $\unicode[STIX]{x1D719}$ is appropriate. (a) Schematic of stretching flow at an interface between low- and high-density regions ( $\unicode[STIX]{x1D719}<0$ and $\unicode[STIX]{x1D719}>0$ , respectively) of contractile swimmers. (b) Snapshot of a steady state in which this interfacial stretching balances the intrinsic coarsening dynamics of model H. (Image courtesy of A. Tiribocchi.)

12.1 Active stress in liquid crystals

As already described, the extensile swimming action of bacteria, or the contractile action of actomyosin, gives rise to an extra mechanical stress term in the NSE. To find its form, consider a rod-shaped particle to represent a single bacterium or a pair of actin filaments cross-linked by a myosin motor, as shown in figure 12(c). A coordinate $\boldsymbol{r}_{i}$ defines the centre of mass of the particle, which has extent $\ell$ along the fore–aft unit vector $\hat{\unicode[STIX]{x1D742}}_{i}$ . The activity creates a flow pattern that can be complicated in the near field, but whose far field is generically described by the lowest-order spherical harmonic compatible with global momentum conservation. The fact that there is no external force on the swimmer rules out a stokeslet contribution as would arise from a point force acting upon it. The lowest-order term is therefore a stresslet, which is equivalent to the action of a force dipole and can be represented as such. We therefore embed equal and opposite point forces at each end of the particle. The direction of the forces will determine whether the stress is contractile (forces pointing inwards, as shown in figure 12 c) or extensile (forces pointing outwards). Following Hatwalne et al. (Reference Hatwalne, Ramaswamy, Rao and Simha2004) (see also Saintillan & Shelley Reference Saintillan and Shelley2007), we write the force density acting on the fluid from $N$ such particles as

Here the summand is the force density from one active particle. Note that $\ell$ is a microscopic length scale whereas $\unicode[STIX]{x1D735}$ is the inverse of a mesoscopic or macroscopic length scale. Thus we Taylor-expanded the Dirac delta functions in the first line above in powers of $\ell \unicode[STIX]{x1D735}$ , assumed to be small. Even in a polar fluid the nematic order parameter $\unicode[STIX]{x1D64C}$ can be defined as usual: $\unicode[STIX]{x1D64C}=\langle \hat{\unicode[STIX]{x1D742}}_{i}\hat{\unicode[STIX]{x1D742}}_{i}\rangle _{meso}-\unicode[STIX]{x1D644}/d$ . More precisely, this can be expressed as

where $c(\boldsymbol{r})=\langle \sum _{i}\unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{i})\rangle _{meso}$ is the number density of the particles and the angle brackets indicate ensemble averaging. Subsituting (12.2) into (12.1), we obtain (ignoring terms $O(\ell ^{3}\unicode[STIX]{x1D6FB}^{3})$ )

Since the force density is related to the stress by $\boldsymbol{f}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}$ , we identify the active stress due to the force dipoles:

where we have introduced an activity parameter $\bar{\unicode[STIX]{x1D701}}=F\ell$ , which is positive for contractile, negative for extensile and zero in equilibrium. (Choosing $\unicode[STIX]{x1D701}=-\bar{\unicode[STIX]{x1D701}}$ matches the notation to that of § 11.1 above.) Here $c(\unicode[STIX]{x1D719})$ is the concentration of the active particles and can be taken to obey $c=c_{0}(\unicode[STIX]{x1D719}+1)/2$ so that $c\simeq c_{0}$ in the active polar phase ( $\unicode[STIX]{x1D719}=1$ ) and $c\simeq 0$ in the passive isotropic phase ( $\unicode[STIX]{x1D719}=-1$ ) for some positive constant $c_{0}$ . Also, for incompressible fluids, the isotropic part of $\unicode[STIX]{x1D748}^{active}$ can be absorbed into the isotropic pressure $P$ . Thus an equally good choice is

The form (12.5) can be used directly in the equations for an active nematic, in which the particles (or more generally their orientational statistics) are symmetric with respect to inversion $\hat{\unicode[STIX]{x1D742}}_{i}\rightarrow -\hat{\unicode[STIX]{x1D742}}_{i}$ . Self-propulsion breaks this symmetry at single-particle level, but nematic phases of self-propelled particles (known as movers) are possible in principle so long as there are equal numbers swimming up and down any chosen spatial axis. Alternatively, there are active particles (known as shakers) that set up a local circulation of fluid but do not self-advect. These can have the full nematic symmetry even at single-particle level. On the other hand, polar active liquid crystals with non-zero $\boldsymbol{p}(\boldsymbol{r},t)=\langle \hat{\unicode[STIX]{x1D742}}\rangle _{meso}$ break this head–tail symmetry by definition. In this case, one either has to carry two separate order parameter fields, $\unicode[STIX]{x1D64C}$ and $\boldsymbol{p}$ , or re-express $\unicode[STIX]{x1D64C}$ locally as a function of $\boldsymbol{p}$ so that only $\boldsymbol{p}(\boldsymbol{r},t)$ need be retained as a dynamical field variable.

There is no general relation of this kind, as exemplified by the case where molecular orientations are distributed uniformly over the unit sphere and normal to it, but with arrows pointing outwards in the upper hemisphere and inwards in the lower. This state has non-zero $\boldsymbol{p}$ but zero $\unicode[STIX]{x1D64C}$ , because, if the arrow heads are removed, the distribution of headless vectors is isotropic. Usually, though, both $\boldsymbol{p}$ and $\unicode[STIX]{x1D64C}$ are non-zero and one can assume $\boldsymbol{p}$ to point along the major axis of $\unicode[STIX]{x1D64C}$ , which is the director $\hat{\boldsymbol{n}}$ . Thus $\boldsymbol{p}=\pm p\hat{\boldsymbol{n}}$ while $\unicode[STIX]{x1D64C}=S(\hat{\boldsymbol{n}}\hat{\boldsymbol{n}}-\unicode[STIX]{x1D644}/d)$ so that the active stress can be written

where $\bar{\unicode[STIX]{x1D701}}$ has been redefined to absorb a factor of $S/p^{2}$ . This expression only differs from (12.4) or (12.5) by an isotropic term of the kind by which they already differ. Note that this form of stress cannot be derived from any free energy functional for our polar liquid crystal. In a quiescent system of non-zero uniform $\boldsymbol{p}$ , the existence of such a free energy structure, from which $\boldsymbol{p}\neq \mathbf{0}$ arises by spontaneous breaking of rotational symmetry, demands that the deviatoric stress involves gradients of $\boldsymbol{p}$ , not $\boldsymbol{p}$ itself, as in (10.14).

Absence of a free energy structure breaks time-reversal symmetry, just as it did via the chemical potential contribution $\unicode[STIX]{x1D706}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|^{2}\neq \unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}$ for active model B above. Nonetheless, just as the remaining chemical potential there was of the passive form, here we retain a free energy functional that generates (effectively passive) elastic stresses as in (10.14), which we now denote as $\unicode[STIX]{x1D70E}^{passive}$ . The chosen form is

Here $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D719})=-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D719}$ , giving the same as (9.2) for passive liquid-crystalline emulsions, except that we have taken the single elastic constant approximation, and set $\unicode[STIX]{x1D6FD}_{2}=0$ , which restricts us to cases of perpendicular anchoring. As in the passive case, this free energy will stabilize a droplet of active polar phase ( $\unicode[STIX]{x1D719}\simeq 1$ and $|\boldsymbol{p}|\simeq 1$ ) within a background fluid of the passive isotropic phase ( $\unicode[STIX]{x1D719}\simeq -1$ and $|\boldsymbol{p}|\rightarrow 0$ ), or vice versa.

12.2 Self-advection and equations of motion

Following the above arguments, we arrive at the dynamical equations for an active polar liquid-crystalline emulsion as follows (Kruse et al. Reference Kruse, Joanny, Julicher, Prost and Sekimoto2005):

This equation set differs from that of polar liquid-crystalline emulsions in (10.4)–(10.7), first via the active stress term in (12.8) as already discussed, and second by the self-advection terms proportional to $w$ in (12.10) and (12.11). These describe the fact that the active material is propelled through space with average local velocity $\langle w\hat{\unicode[STIX]{x1D742}}\rangle _{meso}=w\boldsymbol{p}$ , where $w$ is the swim speed of a particle (or, for actin, a suitably defined treadmilling rate). This motion is additional to the mesoscopically defined fluid velocity $\boldsymbol{v}$ . Accordingly we replace $\boldsymbol{v}\rightarrow \boldsymbol{v}+w\boldsymbol{p}$ in the advective terms of (10.4) and (10.7) to get (12.10) and (12.11) above. Equations such as (12.8)–(12.11) are often referred to as ‘active gel theory’. Here the word ‘gel’ is a slight misnomer, since liquid crystals are not strictly gels. Introducing additional polymeric degrees of freedom allows models of true active gels to be considered, but this area of study remains in its infancy (e.g. Hemingway et al. Reference Hemingway, Maitra, Banerjee, Marchetti, Ramaswamy, Fielding and Cates2015).

Turning to the case of active nematics, the dynamics for these is described by order parameter fields $(\unicode[STIX]{x1D64C},\unicode[STIX]{x1D719},\boldsymbol{v})$ instead of $(\boldsymbol{p},\unicode[STIX]{x1D719},\boldsymbol{v})$ . Since $\boldsymbol{p}$ is zero, the self-advective terms proportional to $w$ are absent, but we still have an active stress in the form (12.4). This is added to $\unicode[STIX]{x1D748}^{passive}$ obeying (10.18) in the Navier–Stokes equation for $\boldsymbol{v}$ with the equations of motion for $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D64C}$ unchanged from the passive case.

12.3 Spontaneous motion of active liquid-crystalline droplets

For simplicity, we first set both the self-advection $w$ and the anchoring term $\unicode[STIX]{x1D6FD}_{1}$ to zero. In this case the equation of motion becomes symmetric with respect to a global inversion $\boldsymbol{p}\rightarrow -\boldsymbol{p}$ . Figure 13 shows what happens when we increase the contractile stress ( $\bar{\unicode[STIX]{x1D701}}>0$ ) from some small values less than the critical activity $\bar{\unicode[STIX]{x1D701}}_{c}$ to a large value greater than $\bar{\unicode[STIX]{x1D701}}_{c}$ (Tjhung et al. Reference Tjhung, Marenduzzo and Cates2012). For $\bar{\unicode[STIX]{x1D701}}<\bar{\unicode[STIX]{x1D701}}_{c}$ , the active gel forms a uniform alignment inside the droplet. The droplet also slightly contracts along $\boldsymbol{p}$ and pumps the surrounding fluid, creating a left–right symmetric fluid flow of stresslet form (see figure 13, left). Since the fluid flow is left–right symmetric, there is no reason for this droplet to move. In other words, the droplet as a whole behaves like a large ‘shaker’ particle.

Figure 13.

For small activity $0<\bar{\unicode[STIX]{x1D701}}<\bar{\unicode[STIX]{x1D701}}_{c}$ , an active polar droplet pumps the surrounding fluid to create a quadrupolar (force-dipolar, or stresslet) fluid flow which is left–right symmetric. At higher activities $\bar{\unicode[STIX]{x1D701}}>\bar{\unicode[STIX]{x1D701}}_{c}$ , the polarization field $\boldsymbol{p}$ becomes unstable with respect to splay and this breaks the left–right symmetry, resulting in translational motion in the direction of splay (here, to the right). (See figure 7 for the definition of splay.)

However, when we increase the contractile activity $\bar{\unicode[STIX]{x1D701}}$ above some threshold $\bar{\unicode[STIX]{x1D701}}_{c}$ , the liquid-crystalline material inside the droplet becomes unstable with respect to splay deformation. This breaks the left–right symmetry in the fluid flow and causes the droplet to spontaneously swim either to the left or to the right (see figure 13, right). These phenomena are unchanged by replacing $\boldsymbol{p}\rightarrow -\boldsymbol{p}$ globally: in particular, the direction of motion corresponds to $\boldsymbol{p}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{p}$ , which is invariant under that replacement. For this reason, the sense of $\boldsymbol{p}$ is not shown in figure 13, and indeed the same figure equally depicts the corresponding phenomena that arise in contractile nematics; the direction of motion is then set by $\hat{\boldsymbol{n}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$ . Again, the director field $\hat{\boldsymbol{n}}$ inside the droplet becomes unstable with respect to splay and the droplet spontaneously swims for $\bar{\unicode[STIX]{x1D701}}>\bar{\unicode[STIX]{x1D701}}_{c}$ . This mechanism for self-propulsion has been proposed for swimming motility in some eukaryotic cells (Hawkins et al. Reference Hawkins, Poincloux, Benichou, Piel, Chavrier and Voituriez2011). Here the source of the contractile stresses comes from the actomyosin contraction in the cytoskeletal bulk of the cell. Note that extensile liquid crystals show instead an instability towards bending at large activities, resulting in banana-shaped droplets that also swim spontaneously (Tjhung et al. Reference Tjhung, Marenduzzo and Cates2012).

In striking (and somewhat related) experimental studies, Sanchez et al. (Reference Sanchez, Chen, Decamp, Heymann and Dogic2012) looked at the dynamics of suspensions of microtubules and kinesin motors in a droplet; under the conditions used, these form an extensile nematic phase. They observed spontaneous motion of the droplet driven by the activity of the kinesin motors, albeit with complex dynamics resulting in non-Brownian diffusion at the droplet scale. A further complication in this system is that the active gel phase-separates into a thin layer near the droplet surface; indeed, these experiments have fuelled theoretical and numerical studies of active nematics in quasi-2D geometries. Both experiments and simulations show complex flow patterns within the 2D film, driven primarily by defect motion. Intriguingly, defects of topological charge $-1/2$ are advected by the flow field $\boldsymbol{v}$ in a quasi-passive manner, whereas those of charge $+1/2$ have additionally an active ballistic motion, somewhat resembling the self-propulsion of contractile emulsion droplets described in the previous paragraph (Giomi et al. Reference Giomi, Bowick, Ma and Marchetti2013). This stems from the fact that the $+1/2$ defect structure breaks spatial symmetry because it looks like an arrow, unlike the $-1/2$ case which is three-fold symmetrical (see figure 8 a). The ballistic separation of defect pairs, which are created by active flow, counters the normal passive evolution in which $\pm 1/2$ defects are attracted to one another and gradually annihilate to give an increasingly ordered state; this competition allows the active system to attain a stationary chaotic flow.

Active polar and nematic emulsion models have also been used to address the biological process of cell division, known as mitosis. During mitosis, the filaments of the cytoskeleton self-organize themselves to create a structure called the ‘mitotic spindle’, which in liquid-crystalline language can be viewed as a pair of $+1$ aster defects. Experimental measurements of the fluctuations in the mitotic spindle show very good agreement with those predicted from the active gel theory outlined above (Brugues & Needleman Reference Brugues and Needleman2014). Once the two asters are created, subsequent division of the droplet can be explained mainly by equilibrium free energy minimization. Specifically, if the anchoring term $\unicode[STIX]{x1D6FD}_{1}$ is much larger than the elastic constant $K$ , division into two equally sized droplets is energetically favourable. However, if the elastic constant $K$ is much larger than $\unicode[STIX]{x1D6FD}_{1}$ , the droplet only elongates, without dividing (Leoni et al. Reference Leoni, Manyuhina, Bowick and Marchetti2017). This is one of several examples where a complex and highly regulated biophysical process can be impersonated by models involving minimal physical ingredients, prompting speculation that the biochemical machinery of the cell is sometimes used to control and exploit the autonomous function of generic active-matter building blocks, rather than to create cellular functionality from scratch. We next address another example of this type.

Figure 14.

(a) Simple model of cell crawling: a 2D active polar droplet on a substrate. When the actin treadmilling rate $w$ is larger than some critical $w_{c}$ , a thin protrusion layer is formed at the leading edge of the droplet. Black arrows indicate the orientation of actin filaments $\boldsymbol{p}$ . In this 2D model, the polarization field is confined to a thin layer near the wall and near the front of the droplet. (b) A similar model in three dimensions can also capture the shape of the crawling cell, showing a fan-shaped lamellipodium, compared to experimental observation (c). In this case there is polarization $\boldsymbol{p}$ throughout the droplet, but treadmilling occurs only in a thin layer near the wall. (See Tjhung et al. (Reference Tjhung, Tiribocchi, Marenduzzo and Cates2015); (c) is from Barnhart et al. (Reference Barnhart, Lee, Keren, Mogilner and Theriot2011).)

12.4 Active emulsion droplets and the physics of cell crawling

Some eukaryotic cell types, such as keratocytes, are able to crawl on a solid substrate such as a glass slide (Verkhovsky et al. Reference Verkhovsky, Svitkina and Borisy1999). The crawling motility is driven mainly by the actin treadmilling process. The actin filaments are found mostly at the leading edge of the crawling cell, where they have polarization $\boldsymbol{p}=\langle \hat{\unicode[STIX]{x1D742}}\rangle _{meso}$ . The filaments polymerize at their front (plus) ends and therefore self-advect in the direction of $\boldsymbol{p}$ , creating forward motion of the entire cell. An important feature is that the actin filaments communicate chemo-mechanically with the supporting wall through structures called ‘focal adhesions’. These provide anchor points for the filaments so that their treadmilling is converted into a tank-treading motion of the cellular perimeter, driving the cell forwards. This communication takes place across the plasma membrane, a lipid bilayer that encloses the cell.

To model this situation within active gel theory, we consider an active polar droplet on a solid substrate, as shown in figure 14(a). For simplicity, the confining plasma membrane is not modelled directly; instead, there is, as usual, an interfacial tension between the interior of the droplet (or now cell) and the exterior. Focal adhesions are accounted for via a no-slip or partial-slip boundary condition on the fluid flow at the wall. We introduce a free energy functional similar to (12.7), but assume either (i) that the polarization field $\boldsymbol{p}$ is confined to a thin layer of thickness $\unicode[STIX]{x1D706}$ close to the wall (and away from the rear of the droplet), with fixed treadmilling rate $w$ , or (ii) that the polarization $\boldsymbol{p}$ extends throughout the droplet but $w$ is non-zero only within a thin layer of thickness $\unicode[STIX]{x1D706}$ near the wall. Such assumptions may be implemented by having explicit spatial dependence of the model parameters, for example $w=w(0)\exp [-z/\unicode[STIX]{x1D706}]$ , where $z$ is a coordinate normal to the wall. Assumption (i) is probably closer to the real behaviour, but assumption (ii) is found to give very similar results in a 2D version of the model (Tjhung et al. Reference Tjhung, Tiribocchi, Marenduzzo and Cates2015), and is much simpler to implement in three dimensions. The confinement of actin and/or self-advection to a thin layer arises in part because the actin network not only adheres to the wall via focal adhesions, but also receives biochemical signals from these that alter its concentration and treadmilling rate.

The equations of motion of the various fields are unchanged from (12.8)–(12.11), although some models of cell crawling also neglect hydrodynamics by putting $\boldsymbol{v}=\mathbf{0}$ , in the manner of model B (Ziebert & Aranson Reference Ziebert and Aranson2013). Since we retain the NSE, we need to specify the boundary condition for $\boldsymbol{v}$ at the solid substrate. We choose a partial-slip boundary condition whereby $v_{z}(z=0)=0$ and $(1-s)\boldsymbol{v}_{\Vert }(z=0)=3\unicode[STIX]{x1D702}s(\unicode[STIX]{x2202}\boldsymbol{v}_{\Vert }/\unicode[STIX]{x2202}z)|_{z=0}$ . Here $\boldsymbol{v}_{\Vert }$ is the tangential velocity and $s$ is the slip parameter (Wolff et al. Reference Wolff, Marenduzzo and Cates2012); $s\rightarrow 0$ recovers the familiar no-slip boundary condition, whereas $s\rightarrow 1$ allows full slip. We assume a boundary condition on $\boldsymbol{p}$ at the substrate such that it lies parallel to the wall: $p_{z}(z=0)=0$ and $(\unicode[STIX]{x2202}\boldsymbol{p}_{\Vert }\unicode[STIX]{x2202}z)|_{z=0}=\mathbf{0}$ , whereas the anchoring term $\unicode[STIX]{x1D6FD}_{1}$ favours perpendicular alignment at the fluid–fluid interface.

Typical results for both 2D and 3D models are shown in figure 14. In two dimensions, for $w$ less than some critical value $w_{c}$ , the droplet crawls forwards in the direction of $\boldsymbol{p}$ . Interestingly, for larger values of $w>w_{c}$ , a thin layer of protrusion is formed at the leading edge of the droplet, which is suggestive of the lamellipodium, a fan-like protrusion that is frequently seen in real crawling cells such as keratocytes (Verkhovsky et al. Reference Verkhovsky, Svitkina and Borisy1999); see figure 14(c). The onset of this protrusion involves a dynamic transition, which becomes increasingly discontinuous when we increase the slip parameter (Tjhung et al. Reference Tjhung, Tiribocchi, Marenduzzo and Cates2015). This prediction might be tested experimentally by measuring the distribution of the projected areas of crawling cells ( $A_{\infty }$ in figure 14 a). In the case of a discontinuous transition, this should show a bimodal distribution.

Physically, the slip parameter $s$ in the model represents the inverse adhesion strength between the filaments and the substrate: if focal adhesions are few, or weak, the actin filaments are more likely to slip backwards instead of being propelled forwards. For the same treadmilling rate $w$ , the protrusion size (or projected area $A_{\infty }$ ) and the crawling speed are both found to decrease with increasing $s$ , that is, with decreasing adhesion. This agrees qualitatively with experimental observations of Barnhart et al. (Reference Barnhart, Lee, Keren, Mogilner and Theriot2011).

In 3D simulations, it is found that the active gel model supports a wide range of cell shapes. These include a ‘fried egg’ structure, with a radially symmetric protrusion and no motion; a fully 3D representation of a crawling lamellipodium as shown in figure 14(b); and a crawling finger-like protrusion, resembling a structure known to biologists as a filopodium. Transitions between these shapes are controlled in part by an interplay between the contractile activity, which promotes spontaneous splay, and the anchoring term $\unicode[STIX]{x1D6FD}_{1}$ , which favours states with $\boldsymbol{p}$ normal to the perimeter of the cell. From the active liquid crystal viewpoint, the lamellipodium is the result of a contractility-induced spontaneous splay in the interior, which then demands a fan-shaped morphology to maintain normal anchoring at the leading edge. This viewpoint is complementary to a more conventional biological one, pointing as it does to mechanistic elements that may stem from the generic properties of active fluids, rather than specific biochemical mechanisms within the cell.

Cellular motility and morphology remain active areas of research where active gel theory and its extensions can provide a physical modelling framework that captures many phenomenologies. In two recent examples, Camley et al. (Reference Camley, Zhang, Zhao, Li, Ben-Jacob, Levine and Rappel2014) modelled a pair of interacting cells that can spontaneously rotate inside a confined geometry, while Lober et al. (Reference Lober, Ziebert and Aranson2015) studied the collective dynamics of multiple motile cells and predicted that collective motion is inhibited with increasing cell–cell adhesion.