Physics of active emulsions

Within the canonical ensemble, a homogeneous binary mixture of volume $V$ can be characterised by the Helmholtz free energy F  =  E  −  TS, which combines the internal energy E and the entropy S of a system. This free energy can be expressed by the partition function (2.1) [102, 107, 108]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\be}{\beta} \newcommand{\kb}{k_{\rm B}} \displaystyle F(T, V, N_A,N_B)=-\kb T \ln Z(T, V, N_A,N_B) . \label{eqn:free_energy_Z} \nonumber \end{align} \tag{ 2.3 }$$

Derivatives of the free energy F are related to thermodynamic quantities that are relevant in our discussion of phase separation. In particular, the entropy is given as $S=-\partial F/\partial T|_{V,N_A,N_B}$, the pressure is $p=-\partial F/\partial V|_{T,N_A,N_B}$ and the chemical potentials read $\mu_A=\partial F/\partial N_A|_{T,V,N_B}$ and $\mu_B=\partial F/\partial N_B|_{T,V,N_A}$.

For simplicity, we focus on an incompressible binary system of constant volume $V=\nu M$ and constant molecular volume $\nu$ of the two components. In this case, adding an A molecule to the system corresponds to removing a B molecule. Consequently, the relevant thermodynamic quantities are the exchange chemical potential $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu$ and the osmotic pressure $\Pi$ [104, 109]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\pp}{\partial} \newcommand{\ba}{\vetr{a}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq:def_chem} \bar\mu = \left. \frac{\pp F}{\pp N_A} \right|_{T,V} =- \left. \frac{\pp F}{\pp N_B} \right|_{T,V} = \nu \left. \frac{\pp f}{\pp \phi} \right|_{T} \, , \nonumber \end{align} \tag{ 2.4a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\pp}{\partial} \newcommand{\be}{\beta} \displaystyle \label{eq:def_p} \Pi = -\left. \frac{\pp F}{\pp V} \right|_{T,N_A} = - f + \phi \, \left. \frac{\pp f}{\pp \phi} \right|_{T} \, , \nonumber \end{align} \tag{ 2.4b }$$

where the number of lattice sites M is slaved to the total volume by $M=V/\nu$. Here, we used the homogeneity of the system, $F=V f(\phi)$, where $f(\phi)$ is the free energy density as a function of the volume fraction $\phi=N_A \nu/V$ of A molecules.

The homogeneous state for a given volume is a stable thermodynamic state if it corresponds to a minimum of the free energy F. This requires that the curvature of the free energy density as a function of volume fraction is convex, i.e. $\newcommand{\g}{\gamma} f''(\phi) \geqslant 0$. The link between stability and curvature of the free energy density can be understood qualitatively: conservation of molecule numbers implies that raising the volume fraction in one spatial region requires lowering it in another. If the free energy density is convex, any such perturbation increases the overall free energy. This can be shown rigorously by considering spatially inhomogeneous perturbations that conserve molecule numbers. We will discuss this approach after introducing the free energy functional in section 2.1.5.

The stability of the homogenous state can also be shown using an ensemble where the particle number N_A is fixed and the volume $V$ can change. This ensemble is governed by the thermodynamic potential $G(N_A,\Pi)=F(N_A,V) + V \Pi$, where $\Pi$ is the osmotic pressure given by equation (2.4b) and the volume $V=\partial G /\partial \Pi$. A homogeneous state with the osmotic pressure $\Pi$ is stable if the free energy G as a function of $\Pi$ is concave, $\partial^2 G/\partial \Pi^2 <0$. The concavity of G with respect to variations of the osmotic pressure can be seen by writing $\partial^2 G/\partial \Pi^2 = \partial V /\partial \Pi = - V \kappa$, where $\kappa$ is the osmotic compressibility $\kappa = - V^{-1} \partial V /\partial \Pi$. For the homogeneous state to be stable, the osmotic pressure should increase as the volume decreases, i.e. $\kappa>0$, to push the system back to its thermodynamic state after a perturbation in volume. This condition is satisfied if the free energy density is convex, $f''(\phi)>0$, since $\kappa=(\phi^2 f''(\phi)){}^{-1}$. In the thermodynamic limit, ensembles become equivalent and thus the convexity of the free energy density determines the thermodynamic stability of the homogeneous state, not only in the ensemble $(N_A,\Pi)$ where the osmotic pressure is imposed but also in the ensemble $(N_A,V)$ where the volume is fixed.

To determine the relevant thermodynamic quantities for phase separation, we need to evaluate the free energy and the partition function given in equation (2.3). Since this is generally difficult, we discuss a mean-field approximation for the homogeneous case of the incompressible binary mixture on a lattice characterised by the Hamiltonian given in equation (2.2). This approximation neglects the spatial correlations between the molecules and is typically a good approximation far away from the critical point where correlations are very short-ranged [110, 111]. Within the mean field approximation, the probability that lattice site n is occupied by A, is given by $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \la \sigma_n \ra= N_A/M$, where $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \la \ldots \ra$ denotes the average in the canonical ensemble and $\phi = N_A/M$ is the volume fraction of A molecules in the system. Due to incompressibility the probability of the site being occupied by B is $N_B/M=1-\phi$. The partition function hence is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle \label{eq:part_function_mean_field} Z \simeq \left| \Omega \right| \exp\left(- \frac{E(\phi)}{k_{\rm B}T}\right) \nonumber \end{align} \tag{ 2.5 }$$

with the internal energy given as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \label{eq_mf_energy} E(\phi)=\frac{zM}{2}\left[e_{AA} \phi^2 +2e_{AB} \phi (1-\phi) +e_{BB} (1-\phi)^2\right] \, , \nonumber \end{align} \tag{ 2.6 }$$

where z is the number of neighbours per lattice site (e.g. z  =  6 for a cubic lattice), and zM/2 is the total number of distinct nearest neighbours. The number $\newcommand{\abs}[1]{|#1|} \abs{\Omega}$ of all possible arrangements on the lattice appearing in equation (2.5),

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \left| \Omega \right| = {M \choose N_A} = {M \choose N_B} = \frac{M !}{N_A ! N_B!} \label{eqn:possible_arrangements} \, , \nonumber \end{align} \tag{ 2.7 }$$

determines the entropy $\newcommand{\kb}{k_{\rm B}} \newcommand{\ri}{{\rm I}} S= \kb \ln \left| \Omega \right|$ for the incompressible binary mixture on a lattice, which is also referred to as mixing entropy. Using equation (2.3), we obtain the free energy density

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \fe(\phi) &\simeq \frac{z}{2\nu} \bigl[ e_{AA} \phi^2 + 2 e_{AB} \phi(1-\phi) +e_{BB} (1-\phi)^2 \bigr] \nonumber \\ &\quad + \frac{k_{\rm B}T}{\nu} \bigl[ \phi\ln \phi +(1-\phi) \ln (1-\phi)\bigr] \, , \label{eq:free_en} \nonumber \end{align} \tag{ 2.8 }$$

where we have used Stirling's approximation, $\newcommand{\s}[1]{s^{(#1)}} \ln{N!} \simeq$ $\newcommand{\s}[1]{s^{(#1)}} N\ln N -N$, to evaluate the factorials. The free energy density can also be written as $\newcommand{\fe}{f} \fe(\phi) = \phi \fe(1) \,+(1-\phi)\fe(0) +$ $\newcommand{\fe}{f} \fe_{\rm mix}(\phi)$, which separates the contribution of the pure systems from the free energy of mixing [104, 106],

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:Fmix} \fe_{\rm mix}(\phi) = \frac{k_{\rm B}T}{\nu} \bigl[ \phi\ln \phi +(1-\phi) \ln (1-\phi) + \chi \phi(1-\phi) \bigr] \nonumber \end{align} \tag{ 2.9 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \label{eq:chi} \chi=\frac{z}{2 k_{\rm B}T} \left(2e_{AB}-e_{AA}-e_{BB}\right) \nonumber \end{align} \tag{ 2.10 }$$

is the Flory–Huggins interaction parameter [2, 3]. $\newcommand{\fe}{f} \fe_{\rm mix}$ captures the competition between the mixing entropy $\newcommand{\ri}{{\rm I}} S=-k_{\rm B} (V/\nu) \left[\phi\ln \phi +(1-\phi) \ln (1-\phi) \right]$ and the molecular interactions characterised by the single parameter $\chi$. In the next section, we will see that both the free energy density (equation (2.8)) and the free energy density of mixing (equation (2.9)) lead to the same phase separation equilibrium. However, the difference will become apparent when we discuss chemical reactions in section 4.

The free energy density $f_{\rm mix}$ is a symmetric function with respect to $\phi=\frac12$. This symmetry stems from considering equal molecular volumes of components A and B and the subtraction of the free energy before mixing. Conversely, if the molecules A and B have different molecular volumes $n_A \nu$ and $n_B \nu$, the free energy of mixing is not symmetric [2, 3]:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:Fmix_asym} \tilde{\fe}_{\rm mix}(\phi) &= \frac{k_{\rm B}T}{\nu} \bigg[ \frac{\phi}{n_A} \ln \phi +\frac{(1-\phi)}{n_B} \ln (1-\phi) + \chi \phi(1-\phi) \bigg] \, , \nonumber \end{align} \tag{ 2.11 }$$

where n_A and n_B denote the non-dimensional molecular size in multiples of the volume $\nu$ of a single lattice site.

Homogeneous states governed by the free energy density $\newcommand{\fe}{f} \fe(\phi)$ are only stable when the free energy density is convex, $\newcommand{\fe}{f} \fe''(\phi)>0$; see section 2.1.1. For the expression given in equation (2.8) this is the case for all $\phi$ in the absence of interactions (e_ij  =  0 for $i,j=A,B$) and when entropic effects dominate. In the presence of interactions however, the free energy density of the homogeneous system can become concave ($\newcommand{\fe}{f} \fe''(\phi)<0$) within a range of volume fractions $\phi$; see figure 1(a). Within this range, the homogeneous state is not stable, implying that the thermodynamic equilibrium state is inhomogeneous.

Zoom In Zoom Out Reset image size

The simplest inhomogeneous state corresponds to two subsystems of different volume fractions, also referred to as phases. The associated free energy can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:freeen} F \simeq V_1 \fe(\phi_1)+V_2 \fe\left(\phi_2\right) \, , \nonumber \end{align} \tag{ 2.12 }$$

where $\newcommand{\al}{\alpha} \phi_\alpha$ and $\newcommand{\al}{\alpha} V_\alpha$ denote the volume fraction and volume of phase $\newcommand{\al}{\alpha} \alpha$, with $\newcommand{\al}{\alpha} \alpha=1,2$. The incompressibility assumption combined with conservation of particles implies $V_1+V_2=V$ and $V_1\phi_1+V_2\phi_2=V\phi$. Consequently, there are only two independent variables in the free energy above, e.g. $\phi_1$ and $V_1$. In equation (2.12) we neglected the energetic contribution of the interface region that separates the two phases. This is valid in the thermodynamic limit where the system and the volumes of the phases are infinitely large, so the energetic contribution of the interface is negligible relative to the contribution of the phases. We will have to refine equation (2.12) when discussing finite systems in section 2.2.

The inhomogeneous state is stable if it corresponds to a minimum of the free energy (2.12) consistent with the imposed constraint of particle number conservation and absence of vacancies. To find this minimum, we differentiate F with respect to $\phi_1$ and $V_1$, respectively, use the relationship $\phi_2= (\phi V-\phi_1 V_1)/(V-V_1)$, and set each expression to zero:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:equal_chem_pot} 0 = \fe'(\phi_1^{(0)})- \fe'(\phi_2^{(0)}) \, , \nonumber \end{align} \tag{ 2.13a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:equal_pressures} 0 &= \fe(\phi_1^{(0)})- \fe(\phi_2^{(0)}) +\bigl(\phi_2^{(0)}-\phi_1^{(0)}\bigr) \fe'(\phi_2^{(0)}) \, , \nonumber \end{align} \tag{ 2.13b }$$

where $\phi_1^{(0)}$ and $\phi_2^{(0)}$ denote the equilibrium volume fractions. The first equation is a balance of the exchange chemical potentials between phases, $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu_1=\bar\mu_2$, with $\newcommand{\vetr}{\boldsymbol} \newcommand{\al}{\alpha} \newcommand{\ba}{\vetr{a}} \bar \mu_\alpha = \bar{\mu} |_{\phi=\phi_\alpha^{(0)}}$ for $\newcommand{\al}{\alpha} \alpha=1,2$; see equation (2.4a)). The second equation corresponds to the balance of the osmotic pressures between the two phases, $\Pi_1=\Pi_2$, with $\newcommand{\al}{\alpha} \Pi_\alpha = \Pi |_{\phi=\phi_\alpha^{(0)}}$; see equation (2.4b).

Obviously, the equations (2.13) are satisfied for homogeneous systems with $\phi_1^{(0)} = \phi_2^{(0)}$. To see that there also exist solutions with $\phi_1^{(0)} \neq \phi_2^{(0)}$, the two conditions can be represented by a graphical tangent construction using the free energy density $\newcommand{\fe}{f} \fe(\phi)$; see figure 1(a). Here, condition (2.13a) implies that the slopes at the equilibrium volume fractions $\phi_1^{(0)}$ and $\phi_2^{(0)}$ are the same and condition (2.13b) states that they are also equal to the slope of the line connecting the points $\newcommand{\fe}{f} (\phi_1^{(0)},\fe(\phi_1^{(0)}))$ and $\newcommand{\fe}{f} (\phi_2^{(0)}, \fe(\phi_2^{(0)}))$. Taken together, these conditions can only be satisfied by a common tangent to the two points. This procedure of finding the equilibrium volume fractions is known as Maxwell's tangent construction or construction of the convex hull. The orange line in figure 1(a) shows the result of such a construction. In fact, inserting condition (2.13b) into equation (2.12) (and using conservation of A particles) shows that this line corresponds to the volume weighted average of the free energy density of the two subsystems. Consequently, the corresponding demixed system has a lower free energy than the mixed system described by the black line in figure 1(a). Clearly, a separation into two phases with volume fractions $\phi_1^{(0)}$ and $\phi_2^{(0)}$ is only possible when the average volume fraction $\phi$ obeys $\phi_1^{(0)} < \phi < \phi_2^{(0)}$. Outside this region, phase separation is not possible and only the homogeneous is stable because $f''(\phi) > 0$.

The parameter region where phase separation is possible can be determined from the solutions $\phi_1^{(0)}$ and $\phi_2^{(0)}$ as a function of the interaction parameter $\chi$; see figure 1(b). The corresponding line is called the binodal line. In the simple case of the symmetric free energy of mixing (equation (2.9)), $\newcommand{\fe}{f} \fe_{\rm mix}(\phi_1^{(0)})=\fe_{\rm mix}(\phi_2^{(0)})$ and $\newcommand{\fe}{f} \fe_{\rm mix}'(\phi_1^{(0)})= \fe_{\rm mix}'(\phi_2^{(0)})=0$. The binodal line is then given by $\chi_{\rm b}(\phi)=\ln(\phi/(1-\phi))/(2\phi-1)$ and phase separation occurs only for $\chi>\chi_{\rm b}$. In particular, the minimal interaction parameter is $\chi_{\rm b}^{\rm min} = 2$, which is obtained at the critical point $\phi=\frac12$ (figure 1(b)). Near the critical point, the equilibrium volume fractions inside each phase obey $\newcommand{\s}[1]{s^{(#1)}} \phi_{1}^{(0)}\simeq \frac12 - [\frac38(\chi-2)]^{1/2}$ and $\newcommand{\s}[1]{s^{(#1)}} \phi_{2}^{(0)}\simeq \frac12 + [\frac38(\chi-2)]^{1/2}$. Note that the same results are obtained when equation (2.8) is used instead of $\newcommand{\fe}{f} \fe_{\rm mix}$, since terms linear in $\phi$ do not alter the conditions given in equation (2.13).

The discussion of the previous section neglected the contribution of the interfacial region on the equilibrium free energy. This interfacial region is always present since the volume fraction is continuous in space and must thus interpolate between the values $\phi_1^{(0)}$ and $\phi_2^{(0)}$ in the two phases. The additional free energy contribution associated with this spatial variation can be estimated within our lattice model. For simplicity, we first consider a one-dimensional system with discrete lattice positions x_n for which the Hamiltonian given in equation (2.2) can be written as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle H(\{\sigma\}) &=\sum_n \big[ e_{AB} \bigl(\sigma_n (1-\sigma_{n+1}) + \sigma_{n+1} (1-\sigma_{n}) \bigr) \nonumber \\ & \quad + e_{AA} \sigma_{n} \sigma_{n+1} + e_{BB} (1- \sigma_{n}) (1- \sigma_{n+1}) \big] \, . \label{eq:spatvar} \nonumber \end{align} \tag{ 2.14 }$$

We proceed analogously to section 2.1 and perform a mean-field approximation after rewriting the coupling terms using $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} 2 \sigma_n (1-\sigma_{n+1}) = \left(\sigma_n-\sigma_{n+1}\right){}^2 -\sigma_n^2-\sigma_{n+1}^2+2\sigma_n$. Additionally, generalising to three dimensions and taking the continuum limit, we replace $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\dd}{{\rm d}} \sum_n \mapsto \nu^{-1}\int \dd^{3}r$, $\newcommand{\vetr}{\boldsymbol} \newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\bbr}{\vetr{r}} \mean{\sigma_n} \mapsto \phi(\bbr)$, and $\newcommand{\vetr}{\boldsymbol} \newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\s}[1]{s^{(#1)}} \mean{\sigma_{n+1} - \sigma_n} \mapsto \nu^{\frac{1}{3}} \nabla\phi(\vetr r)$. Hence,

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\ri}{{\rm I}} \newcommand{\dd}{{\rm d}} \newcommand{\bbr}{\vetr{r}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\B}{c_{B}} \newcommand{\kb}{k_{\rm B}} \newcommand{\vetr}{\boldsymbol} \displaystyle E & \simeq \frac{\kb T \chi}{\nu} \int \dd^{3}r \Bigl[ \phi(\bbr) \bigl(1-\phi(\bbr) \bigr) + \frac{\nu^{\frac23}}{2} \left| \nabla \phi(\bbr) \right|^2 \Bigr] \nonumber \\ & \quad + \frac{z}{2\nu}\int \dd^{3}r \Big[e_{AA} \, \phi(\bbr) + e_{BB} \bigl(1-\phi(\bbr) \bigr) \Big] \, , \nonumber \end{align} \tag{ 2.15 }$$

where $\chi$ is the Flory–Huggins parameter given in equation (2.10). The associated free energy F can be expressed as a functional of $\newcommand{\vetr}{\boldsymbol} \phi(\vetr r)$,

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\dd}{{\rm d}} \newcommand{\be}{\beta} \newcommand{\B}{c_{B}} \newcommand{\kb}{k_{\rm B}} \displaystyle F &=\frac{1}{\nu} \int \dd^{3}r \, \Bigl[ \frac{\kappa}{2\nu} \left|\nabla \phi \right|^2 + \frac{z}{2}\left(e_{AA} \, \phi + e_{BB} \left(1-\phi \right) \right) \nonumber \\ &\quad + \kb T \big[ \phi \ln (\phi) +(1 - \phi)\ln (1- \phi) + \chi \, \phi (1-\phi) \big] \Bigr] , \label{eq:fmix_cont_bin} \nonumber \end{align} \tag{ 2.16 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\be}{\beta} \displaystyle \label{eq:kappa_lattice} \kappa = k_{\rm B}T \chi \, \nu^{\frac53} \nonumber \end{align} \tag{ 2.17 }$$

characterises the change of free energy density due to concentration inhomogeneities. Since the interaction parameter $\chi>2$ in the phase separating regime, we have $\kappa >0$ and the gradient term thus penalises spatial inhomogeneities. We identify the integrand of the free energy given in equation (2.16) as the total free energy density $\newcommand{\fe}{f} \newcommand{\fetot}{f_{\rm tot}} \fetot$. Expressing it in terms of concentrations $c = \phi/\nu$, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\fetot}{f_{\rm tot}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:freefuncttot} \fetot(c, \nabla c) = f_{\rm mix}(c) + f_0(c) + \frac{\kappa}{2} \left|\nabla c\right|^2 \, , \nonumber \end{align} \tag{ 2.18 }$$

where $f_{\rm mix}$ follows from equation (2.9) and $\newcommand{\ri}{{\rm I}} f_0(c) \,=$$\newcommand{\ri}{{\rm I}} (z/2) \left(e_{AA} \, c + e_{BB} \left(\nu^{-1}-c \right) \right)$ is the free energy of the pure system before mixing. The free energy density (2.18) captures many features of phase separation, such as an interaction term favouring phase separation that can outcompete the mixing entropy. It is sometimes useful to consider simpler forms of the free energy avoiding the logarithmic terms in the mixing entropy (equation (2.9)), which are not required to discuss the underlying basic principles. An example for such a simplification is the Ginzburg–Landau free energy.

Near the critical point, the physics of phase separation is fully captured by the shape of the free energy close to the critical concentration c_c [7, 112–114], which is obtained by expanding the free energy to fourth order in (c − c_c). Such an expansion is also very useful as a simplification away from the critical point. Expanding the total free energy at $\newcommand{\cc}{c_{\rm c}} \cc$ leads to the asymmetric Ginzburg–Landau free energy,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\dd}{{\rm d}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:GLfree} F_{\rm GL}[c]=\int \dd^{3}r \left(\fe_{\rm GL}(c)+\frac{\kappa}{2} \left|\nabla c\right|^2 \right) \; , \nonumber \end{align} \tag{ 2.19 }$$

with the corresponding free energy density

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\cc}{c_{\rm c}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:GL_free_energy_density_asym} \fe_{\rm GL}(c) &= \tilde{b} \left(c - \cc\right) -\frac{b}{2} \left(c - \cc\right)^2 \nonumber \\ \nonumber & \quad + \frac{\tilde a}{3} \left(c-\cc\right)^3 + \frac{a}{4} \left(c-\cc\right)^4 \, . \nonumber \end{align} \tag{ 2.20 }$$

This free energy density is parameterised by the coefficients a, $\tilde a$, b, $\tilde b$, and $\newcommand{\cc}{c_{\rm c}} \cc$. As discussed in section 2.1.3, phase separation equilibrium is not affected by the linear term, thus we choose $\tilde{b}=0$, and for simplicity we also consider the special case $\tilde{a}=0$. The corresponding Ginzburg–Landau free energy density then assumes a bi-quadratic form,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\cc}{c_{\rm c}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:GL_free_energy_density} \fe_{\rm GL}(c) = -\frac{b}{2} \left(c - \cc\right)^2 + \frac{a}{4} \left(c-\cc\right)^4 \, , \nonumber \end{align} \tag{ 2.21 }$$

which is symmetric around the concentration $\newcommand{\cc}{c_{\rm c}} \cc$. This free energy density has a concave region if b  >  0 and the parameter a  >  0 characterises the convex branches of the energy density. Using equation (2.13a), the equilibrium concentrations within each phase separated domain are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cc}{c_{\rm c}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle c_1^{(0)} =\cc - \sqrt{b/ a} \, , \nonumber \end{align} \tag{ 2.22a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cc}{c_{\rm c}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle c_2^{(0)} =\cc + \sqrt{b/ a} \, . \nonumber \end{align} \tag{ 2.22b }$$

In particular, expansion of $\newcommand{\fe}{f} \newcommand{\fetot}{f_{\rm tot}} \fetot - f_0(c)$ (equation (2.18)) around $\newcommand{\cc}{c_{\rm c}} \cc=1/(2\nu)$ links the parameters a and b with the molecular parameters of the lattice model,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\be}{\beta} \newcommand{\kb}{k_{\rm B}} \displaystyle \label{eq_a_b} a=\frac{16}{3} \kb T \nu^3, \quad b= 2 (\chi-2) \kb T\nu \; . \nonumber \end{align} \tag{ 2.23 }$$

Using the expressions above in equation (2.22) we consistently obtain the equilibrium concentration found in section 2.1.3.

We start by determining the stationary states $\newcommand{\vetr}{\boldsymbol} c_*(\vetr r)$ of the bi-quadratic Ginzburg–Landau free energy $F_{\rm GL}$ given in equation (2.19). Such states have an extremal free energy subjected to the constraint that the number of A and B molecules are conserved, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\dd}{{\rm d}} \newcommand{\bbr}{\vetr{r}} \newcommand{\ba}{\vetr{a}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq:particle_conservation} \bar{c} = \frac1V \int_V \dd^{3}r \, c(\bbr) , \nonumber \end{align} \tag{ 2.24 }$$

where $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar{c}=N_A/V$ denotes the mean concentration of A molecules in the mixture of the finite volume $V$. To enforce this constraint, we introduce a Lagrange multiplier $\newcommand{\la}{\langle} \lambda$ and vary the functional $\newcommand{\vetr}{\boldsymbol} \newcommand{\dd}{{\rm d}} \newcommand{\la}{\langle} F_{\rm GL} -\lambda \int_V \dd^{3}r \, c(\vetr r)$. Here, the functional derivative of the free energy is the generalization of the exchange chemical potential generalized to inhomogeneous systems,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ba}{\vetr{a}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \bar{\mu} = \frac{\delta F}{\delta c} \label{def_inhomo_chem_pot} , \nonumber \end{align} \tag{ 2.25 }$$

which reads $\newcommand{\vetr}{\boldsymbol} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \newcommand{\bbr}{\vetr{r}} \newcommand{\ri}{{\rm I}} \bar{\mu}=a \left(c(\bbr)-\cc\right){}^{3}- b \left(c(\bbr)-\cc \right)- \kappa \nabla^2 c(\bbr)$ for $F_{\rm GL}$. Consequently, the Euler–Lagrange equation for the stationary state is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ba}{\vetr{a}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq:chem_pot_GL} \bar{\mu} = \lambda \, , \nonumber \end{align} \tag{ 2.26 }$$

where we dropped a boundary term proportional to $\kappa \nabla c$, assuming appropriate boundary conditions at the system boundary. The stationary states thus correspond to a spatially uniform exchange chemical potential $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar{\mu}$, which may be realised for both homogeneous and inhomogeneous concentration profiles $\newcommand{\vetr}{\boldsymbol} c_*(\vetr r)$.

We start by considering the spatially homogeneous equilibrium states $\newcommand{\vetr}{\boldsymbol} c_*(\vetr r) = c_0$. Particle conservation implies $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} c_0 = \bar c$ (equation (2.24)) and the Euler–Langrange equations read $\newcommand{\vetr}{\boldsymbol} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \newcommand{\la}{\langle} \lambda= a(\bar{c} - \cc){}^3 - b(\bar{c} - \cc)$. The homogeneous state $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} c(\vetr r) = \bar c$ is thus an extremal state of the free energy $F_{\rm GL}$ and we will check whether it corresponds to a minimum in the next section.

We showed above that states with coexisting phases can be stable in some regions of the phase diagram shown in figure 1(b). Here, we determine the concentration profile that connects the two phases from extremising the free energy $F_{\rm GL}$. In the following, we restrict ourselves to a flat interface oriented perpendicular to the x-axis at the position x  =  0 with a concentration $\newcommand{\cc}{c_{\rm c}} c(x=0) = \cc$. For simplicity, we extend the system to infinity while keeping the position and concentration value of the interface fixed. For the symmetric free energy density (equation (2.21)), this implies $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar \mu =0$, thus $\newcommand{\la}{\langle} \lambda=0$, at the interface, and in the case of an extremal inhomogeneous state, $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar \mu =0$ at all positions. Far away from the interface, the concentration inside the phases should be governed by the equilibrium concentrations (equation (2.22)),

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \lim_{x\rightarrow -\infty} c(x) = c_1^{(0)}\, , \quad \lim_{x\rightarrow \infty} c(x) = c_2^{(0)} . \label{eq:interface_bcs} \nonumber \end{align} \tag{ 2.27 }$$

With the boundary conditions above the unique solution is [42]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\cc}{c_{\rm c}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle \label{eq:tanh} c_{\rm I}(x) = \cc + \sqrt{\frac{b}{a}}\tanh\left(\sqrt{\frac{b}{2 \kappa}}x\right) \ . \nonumber \end{align} \tag{ 2.28 }$$

Since this interfacial profile varies substantially only within the region $\newcommand{\s}[1]{s^{(#1)}} |x| \lesssim \sqrt{2\kappa/b}$, we introduce the interfacial width

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle \label{eq:Iwidth} w=\sqrt{\frac{2\kappa}{b}} \simeq \nu^{\frac13} \sqrt{\frac{\chi}{\chi-2}} \ . \nonumber \end{align} \tag{ 2.29 }$$

The right-hand side relates w to the lattice model using equations (2.17) and (2.23). In the limit of strong phase separation (large $\chi$) the interfacial width approaches the linear dimension $\nu^{1/3}$ of the molecules.

The homogeneous and inhomogeneous stationary states of the Ginzburg–Landau free energy, $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} c_*(x)=\bar c$ and $c_*(x)=c_{\rm I}(x)$, respectively, can be either stable or unstable. They are stable if they correspond to a free energy minimum, i.e. if all small concentration perturbations increase the free energy. To test this, we consider concentration profiles $\newcommand{\e}{{\rm e}} c=c_* +\epsilon$, where $\newcommand{\vetr}{\boldsymbol} \newcommand{\e}{{\rm e}} \epsilon(\vetr r)$ is a small, position dependent concentration perturbation. To quadratic order, the change in the free energy due to this perturbation is

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\cc}{c_{\rm c}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\B}{c_{B}} \newcommand{\diff}{{{\rm d}}} \displaystyle \Delta F [c_*, \epsilon] &= F_{\rm GL} (c_* + \epsilon) - F_{\rm GL} (c_*) \label{eq:triF} \nonumber \\ &\simeq \int \diff^3r \Big[\frac{\epsilon^2}{2}\bigl(3a (c_*-\cc)^2-b\bigr) \nonumber +\frac{\kappa}{2} \left(\nabla \epsilon \right)^2 \Big]. \nonumber \end{align} \tag{ 2.30 }$$

The state $\newcommand{\vetr}{\boldsymbol} c_*(\vetr r)$ is stable if all perturbations increase the free energy, i.e. if $\newcommand{\e}{{\rm e}} \Delta F[c_*, \epsilon] > 0$ for all $\newcommand{\vetr}{\boldsymbol} \newcommand{\e}{{\rm e}} \epsilon(\vetr r)$. In the case of the homogeneous state $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} c_*(\vetr r) = \bar c$ both terms in the integrand are positive if $\newcommand{\abs}[1]{|#1|} \newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \abs{\bar c - \cc} > \sqrt{b/(3a)}$, which implies $\newcommand{\e}{{\rm e}} \Delta F[c_*, \epsilon]>0$. Conversely, $\newcommand{\e}{{\rm e}} \Delta F[c_*, \epsilon]$ can be negative for $\newcommand{\abs}[1]{|#1|} \newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \abs{\bar c - \cc} < \sqrt{b/(3a)}$ in sufficiently large systems, e.g. for the perturbation $\newcommand{\e}{{\rm e}} \epsilon(x) \propto \tanh(x/w)$, which implies that the homogeneous state is unstable for these parameters. Consequently, the stationary homogeneous state can be either stable or unstable to infinitesimal perturbations. In contrast, the inhomogeneous state $c_*(x)=c_{\rm I}(x)$ given by equation (2.28) is always stable if it is a stationary state, i.e. when $\newcommand{\abs}[1]{|#1|} \newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \abs{\bar c - \cc} < \sqrt{b/a}$; see appendix A.

Taken together, we can distinguish three different parameter regimes with different stable stationary states. For mean concentrations $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar c$ far away from the symmetry point $\newcommand{\cc}{c_{\rm c}} \cc$, i.e. when $\newcommand{\abs}[1]{|#1|} \newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \abs{\bar c - \cc} > \sqrt{b/a}$, the homogeneous state is the only stable one. Conversely, when $\newcommand{\abs}[1]{|#1|} \newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \abs{\bar c - \cc} < \sqrt{b/(3a)}$, only the inhomogeneous state is stable and phase separation will thus happen spontaneously. This region is known as the spinodal decomposition region [115–121] which is enclosed by the spinodal line (dashed line in figure 1(b)). Between the spinodal region and the homogeneous region, for $\newcommand{\abs}[1]{|#1|} \newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cc}{c_{\rm c}} \newcommand{\ba}{\vetr{a}} \sqrt{b/(3a)} < \abs{\bar c - \cc} < \sqrt{b/a}$, both states are stable to infinitesimal perturbations. In this case, phases can only originate from the homogeneous state by large fluctuations, known as nucleation events. Consequently, the respective region in the phase diagram is known as the nucleation and growth regime. All three phases are shown in the phase diagram in figure 1(b).

One important difference between droplets and the condensed phases that we discussed so far is the curvature of the droplet interface, which is inevitable due to the finite size. The surface tension $\newcommand{\g}{\gamma} \gamma$ of this curved interface affects the equilibrium concentrations inside and outside the droplet, which we denote by $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqIn}{\ceq_{\rm in}} \cEqIn$ and $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqOut}{\ceq_{\rm out}} \cEqOut$, respectively. In the following, we consider the case where the surface tension $\newcommand{\g}{\gamma} \gamma$ is constant and independent of the interface curvature, which is valid for droplets large compared to the Tolman length [127, 128]. To derive how the equilibrium concentrations depend on the droplet curvature, we write equation (2.31) for a spherical droplet of radius R,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\cIn}{c_{\rm in}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\Vd}{V_{\rm d}} \newcommand{\Vsys}{V} \displaystyle F = \Vd \fe(\cIn)+(\Vsys - \Vd) \fe(c_{\rm out}) + 4\pi R^2 \gamma \, , \nonumber \end{align} \tag{ 2.42 }$$

where $\newcommand{\Vd}{V_{\rm d}} \Vd=\frac{4\pi}{3} R^3$ denotes the droplet volume and $\newcommand{\Vsys}{V} \Vsys$ is the volume of the system. Using particle conservation, $\newcommand{\vetr}{\boldsymbol} \newcommand{\Vsys}{V} \newcommand{\Vd}{V_{\rm d}} \newcommand{\cIn}{c_{\rm in}} \newcommand{\ba}{\vetr{a}} c_{\rm out}= (\Vsys \bar c -\Vd \cIn)/(\Vsys-\Vd)$, and minimising the free energy with respect to $\newcommand{\cIn}{c_{\rm in}} \cIn$ and $\newcommand{\Vd}{V_{\rm d}} \Vd$, analogously to section 2.1.3, we obtain the equilibrium conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\cEqIn}{\ceq_{\rm in}} \newcommand{\cEq}{\ceq} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \displaystyle \label{eq:shift1} 0 = \fe'(\cEqIn)- \fe'(\cEqOut) \, , \nonumber \end{align} \tag{ 2.43a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\cEqIn}{\ceq_{\rm in}} \newcommand{\cEq}{\ceq} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \displaystyle \label{eq:shift2} 0 &= \fe(\cEqIn)- \fe(\cEqOut) \nonumber \\[-3pt] & \quad +\left(\cEqOut-\cEqIn\right) \fe'(\cEqOut) +\frac{2 \gamma}{R} \; . \nonumber \end{align} \tag{ 2.43b }$$

These conditions are fulfilled at the equilibrium concentrations $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqIn}{\ceq_{\rm in}} \cEqIn$ and $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqOut}{\ceq_{\rm out}} \cEqOut$ inside and outside the droplet, respectively. Comparing these expressions to equation (2.13) in the thermodynamic limit, we find that the pressure balance (2.43b) contains an additional term $\newcommand{\g}{\gamma} 2\gamma/R$, which is known as the Laplace pressure. Graphically, the Laplace pressure corresponds to the free energy difference of the tangents in the Maxwell construction; see figure 1(a). The Laplace pressure is proportional to the interface curvature R⁻¹ and thus disappears in the thermodynamic limit ($\newcommand{\ri}{{\rm I}} R \rightarrow \infty$).

To derive approximate expressions of the equilibrium concentrations $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqIn}{\ceq_{\rm in}} \cEqIn$ and $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqOut}{\ceq_{\rm out}} \cEqOut$, we expand $c^{\rm eq}_{\rm in/out} = c^{(0)}_{\rm in/out}+ \delta c_{\rm in/out}$ in equation (2.43) to linear order in $\delta c_{\rm in/out}$. Here, $c^{(0)}_{\rm in}$ and $c^{(0)}_{\rm out}$ denote the equilibrium concentrations in the thermodynamic limit in the condensed and dilute phase, respectively, so $\delta c_{\rm in/out}$ captures the effects of Laplace pressure. We find

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \displaystyle \label{eq:phiout} \delta c_{\rm out} \simeq \frac{2 \gamma}{\bigl(\cBaseIn - \cBaseOut\bigr) \fe''\bigl(\cBaseOut\bigr) R} \, , \nonumber \end{align} \tag{ 2.44a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:phiin} \delta c_{\rm in} \simeq \frac{\fe''\bigl(c^{(0)}_{\rm out}\bigr)}{\fe''\bigl(c^{(0)}_{\rm in}\bigr)} \delta c_{\rm out} , \nonumber \end{align} \tag{ 2.44b }$$

which are known as the Gibbs–Thomson relations. Since both expressions are positive, the Laplace pressure elevates the concentrations both inside and outside the droplet. This effect is stronger for smaller droplets, which becomes explicit when writing the equilibrium concentrations as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle c^{\rm eq}_{\rm out}=c^{\rm (0)}_{\rm out} \left(1+\frac{\ell_{\gamma,{\rm out}}}{R}\right) \; , \nonumber \end{align} \tag{ 2.45a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \label{eq:cineq} c^{\rm eq}_{\rm in}=c^{\rm (0)}_{\rm in} \left(1+\frac{\ell_{\gamma,{\rm in}}}{R} \right) \; , \nonumber \end{align} \tag{ 2.45b }$$

where we defined for both phases the capillary lengths

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \displaystyle \ell_{\gamma,{\rm out}} = \frac{2 \gamma}{\bigl(\cBaseIn - \cBaseOut\bigr) \fe''\bigl(\cBaseOut\bigr)\cBaseOut} \, , \nonumber \end{align} \tag{ 2.45c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \displaystyle \ell_{\gamma,{\rm in}} = \frac{\fe''(\cBaseOut)\cBaseOut }{\fe''(\cBaseIn)\cBaseIn} \, \ell_{\gamma,{\rm out}} . \nonumber \end{align} \tag{ 2.45d }$$

In the following, we are interested in the limit of strong phase separation with $\newcommand{\g}{\gamma} c^{(0)}_{\rm in} \gg c^{(0)}_{\rm out}$ and thus $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_{\gamma,{\rm in}} \ll \ell_{\gamma,{\rm out}}$. In this case, the impact of the Laplace pressure on the equilibrium concentration inside the droplet can be neglected, i.e. $\newcommand{\s}[1]{s^{(#1)}} c^{\rm eq}_{\rm in} \simeq c^{\rm (0)}_{\rm in}$. This leaves us with a single capillary length, $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma = \ell_{\gamma,{\rm out}}$. When the phase outside is dilute, i.e. when the chemical potential can be written as $\newcommand{\vetr}{\boldsymbol} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\ba}{\vetr{a}} \bar{\mu}(c) \simeq k_{\rm B}T \ln (\nu c)$ [129], we find

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\kb}{k_{\rm B}} \displaystyle \ell_{\gamma} \simeq \frac{2 \gamma}{c^{(0)}_{\rm in} \kb T} . \label{eq:lc_app} \nonumber \end{align} \tag{ 2.46 }$$

If we additionally assume that the condensed phase is highly packed such that $\newcommand{\cbase}{c^{(0)}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} \cBaseIn \simeq \nu^{-1}$, equation (2.46) provides a useful estimate of the capillary length when the surface tension $\newcommand{\g}{\gamma} \gamma$ is known [129]. Using equation (2.34) from our lattice model, the capillary length can also be expressed as $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_{\gamma} \simeq \chi^{1/2} \left(\chi -2\right){}^{3/2}\nu^{1/3}$. This expression demonstrates that for interaction parameters not too far away from the critical value, $\chi_{\rm b}^{\rm min} = 2$, it is typically on the order of the molecular length scale $\nu^{1/3}$. Consequently, we have $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma \ll R$ and the increase of the equilibrium concentrations predicted by the Gibbs Thomson relations (2.45) is actually small, supporting the validity of the linear approximation.

The dynamics of the droplet size and its shape are linked to the movement of its interface, which we assume to be thin compared to the droplet size in the following. To describe the dynamics of the interface, we consider a spherical coordinate system $(r, \varphi, \theta)$ centered on the droplet. Assuming the interface does not deviate strongly from a spherical shape, we parameterize its shape $\newcommand{\vetr}{\boldsymbol} \vetr{R}(\varphi, \theta; t)= \mathcal{R}(\varphi, \theta; t) \vetr{e}_r$ by the radial distance $\mathcal{R}(\varphi, \theta; t)$ as a function of the polar angle $\varphi$ and the azimuthal angle $\theta$. The movement of the interface is most naturally described in the local coordinate system spanned by the two tangential directions $\newcommand{\vetr}{\boldsymbol} \vetr{e}_1=\partial \vetr{R}/\partial \varphi$ and $\newcommand{\vetr}{\boldsymbol} \vetr{e}_2=\partial \vetr{R}/\partial \theta$ and the outward normal vector $\newcommand{\vetr}{\boldsymbol} \vetr{n}=\frac{\vetr{e}_1 \times \vetr{e}_2}{|\vetr{e}_1 \times \vetr{e}_2|}$. Note that the droplet shape is only affected by the normal component $v_n$ of the interfacial velocity, while the tangential components transport material along the interface. Material conservation implies that this normal component is proportional to the net material flux toward the interface,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\cEqIn}{\ceq_{\rm in}} \newcommand{\cEq}{\ceq} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq:interface_vel} v_n =\frac{\vetr{j}_{\rm in}-\vetr{j}_{\rm out}}{\cEqIn-\cEqOut} \cdot \vetr{n} \; , \nonumber \end{align} \tag{ 2.47 }$$

where $\newcommand{\vetr}{\boldsymbol} \newcommand{\e}{{\rm e}} \vetr{j}_{\rm in}=\lim_{\epsilon\to 0} \vetr{j} (\vetr{R} - \epsilon\vetr{n})$ and $\newcommand{\vetr}{\boldsymbol} \newcommand{\e}{{\rm e}} \vetr{j}_{\rm out}=\lim_{\epsilon\to 0} \vetr{j} (\vetr{R} + \epsilon\vetr{n})$ are the local material fluxes right inside and outside of the interface, respectively. Expressing the time evolution of the interface in the spherical coordinate system gives $\newcommand{\vetr}{\boldsymbol} \partial_t \vetr{R}= \partial_t \mathcal{R} \, \vetr{e}_r$, while in the local coordinate system of the interface,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \partial_t \vetr{R} = v_n \vetr{n} + v_{{t},1} \vetr{e}_1 + v_{{t},2} \vetr{e}_2 \, . \nonumber \end{align} \tag{ 2.48 }$$

We can use the connection between the local and the global coordinate system to identify the conditions $\newcommand{\vetr}{\boldsymbol} \partial_t \vetr{R} \cdot \vetr{e}_\theta = 0$ and $\newcommand{\vetr}{\boldsymbol} \partial_t \vetr{R} \cdot \vetr{e}_\varphi = 0$, which can be used to obtain the in-plane velocity components $v_{{t},1}$ and $v_{{t},2}$. The radial interface velocity then reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq_interface_evolution} \partial_t \mathcal{R} = v_n \left[ 1 + \left(\frac{\partial_\theta \vetr{R}}{R}\right)^2 + \left(\frac{\partial_\varphi \vetr{R}}{R \sin(\theta) } \right)^2 \right]^{\frac12} \; . \nonumber \end{align} \tag{ 2.49 }$$

In the case where the dynamics within the phases are described by the diffusion equation (2.40), the material flux is given by $\newcommand{\vetr}{\boldsymbol} \vetr j=-D \nabla c$ and equation (2.49) directly determines the time evolution of the interface.

Before we consider shape perturbations in the subsequent sections, we focus here on a spherical droplet, ${\mathcal R}(\varphi, \theta; t) = R(t)$, in a spherically symmetric system. To derive its growth dynamics, we employ the quasi-static approximation, which assumes that the droplet radius varies slowly such that transients in the diffusion equation (2.40) can be neglected. Within this approximation, the diffusion equation (2.40) reduces to a Laplace equation inside and outside the droplet,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\pp}{\partial} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \displaystyle \label{eq:laplace} 0 \simeq \nabla^2 c(r)=\frac{1}{r^2}\frac{\pp}{\pp r}\left(r^2 \frac{\pp c}{\pp r} \right) \, , \nonumber \end{align} \tag{ 2.50 }$$

where we have written the Laplace operator in spherical coordinates considering that there is no polar and azimuthal dependence of the concentration field. The associated boundary conditions are given by the Gibbs–Thomson relations (2.45) at the droplet interface and a differentiability condition at the droplet centre. Moreover, we consider the case where the droplet is embedded in a large system and the concentration far away is fixed to $c_\infty$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pp}{\partial} \newcommand{\be}{\beta} \displaystyle \pp_r c (r) = 0 \, \qquad {{\rm at}} \qquad r=0 \; , \nonumber \end{align} \tag{ 2.51a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \lim_{r\rightarrow \infty} c (r) = c_\infty \; . \nonumber \end{align} \tag{ 2.51b }$$

Using these boundary conditions, the solutions inside and outside the droplet read

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle \label{eq:diff} c(r) = c_\infty + \left(c_{\rm out}^{\rm eq} -c_\infty\right)\frac{R}{r} \; , \qquad r> R \; , \nonumber \end{align} \tag{ 2.52a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \displaystyle c(r) = c^{\rm eq}_{\rm in} \; ,\quad \quad r<R \, , \nonumber \end{align} \tag{ 2.52b }$$

which are illustrated in figure 2(a). These solutions imply the fluxes $\newcommand{\vetr}{\boldsymbol} \vetr j_{\rm in} = \vetr 0$ and $\newcommand{\vetr}{\boldsymbol} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqOut}{\ceq_{\rm out}} \vetr j_{\rm out} = DR^{-1}(\cEqOut - c_\infty) \vetr e_r$ inside and outside of the interface, respectively. The growth rate of the droplet then follows from equation (2.47),

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\dd}{{\rm d}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\supSat}{\varepsilon} \displaystyle \label{eq:dRdt} \frac{\dd R}{\dd t} = \frac{D \, \cBaseOut}{R \, \cBaseIn} \left(\supSat - \frac{\ell_\gamma}{R} \right) , \nonumber \end{align} \tag{ 2.53 }$$

where we consider the case of strong phase separation ($\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\g}{\gamma} \cBaseIn \gg \cBaseOut$). We have defined the supersaturation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\supSat}{\varepsilon} \displaystyle \supSat = \frac{c_\infty}{\cBaseOut} - 1 , \nonumber \end{align} \tag{ 2.54 }$$

which measures the excess concentration relative to the equilibrium concentration $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \cBaseOut$ in the dilute phase. Equation (2.53) shows that the droplet only grows in sufficiently supersaturated environments where $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \varepsilon > \ell_\gamma R^{-1}$. The droplet dynamics are thus directly linked to the concentration of droplet material in its environment. In particular, we can define the critical radius $\newcommand{\supSat}{\varepsilon} \newcommand{\R}{c_{R}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\Rcrit}{R_{\rm c}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \Rcrit=\ell_\gamma\supSat^{-1}$, below which the droplet shrinks; see figure 2(b). Although this deterministic description cannot account for the spontaneous emergence of droplets, the critical radius is key to estimating the frequency of such nucleation events [130, 131]. In essence, nucleation relies on large fluctuations that spontaneously enrich droplet material in a region of radius $\newcommand{\R}{c_{R}} \newcommand{\Rcrit}{R_{\rm c}} \Rcrit$. In this case, the resulting droplet starts growing spontaneously according to equation (2.53).

Zoom In Zoom Out Reset image size

So far, we focused on a single droplet, but most phase separated systems contain many droplets. In such emulsions, large droplets typically grow at the expense of smaller droplets, which vanish eventually. This phenomenon is referred to as Ostwald ripening [4].

In the following we consider the interactions of many droplets that are far apart from each other in a dilute system with small supersaturation. In this case, nucleation events are rare, and the surrounding of droplets can be considered to be spherically symmetric with a common concentration $c_\infty$ far away from each droplet. This implies a common supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat$, which depends on time and mediates the interactions between the droplets. As the supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat$ can be determined from the total amount of material, the state of the system is fully specified by the radii R_i of the N droplets in the system. Their dynamics follows from equation (2.53) and reads [5]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\cInfty}{c_\infty} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cIn}{c_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\difffrac}[2]{\frac{\diff #1}{\diff #2}} \newcommand{\diff}{{{\rm d}}} \displaystyle \difffrac{R_i(t)}{t} = \frac{D \, \cBaseOut}{R_i(t) \, \cBaseIn} \left(\frac{\cInfty(t)}{\cBaseOut}-1 -\frac{\ell_\gamma}{R_i(t)} \right) \, , \nonumber \end{align} \tag{ 2.55a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\ba}{\vetr{a}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\Vsys}{V} \newcommand{\cbase}{c^{(0)}} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq:multidrop_solvent} \bar c \Vsys &= \cBaseIn \sum_{i=1}^N \frac{4\pi}3 R_i(t)^3 + c_\infty(t) \left[ V- \sum_{i=1}^N \frac{4\pi}3 R_i(t)^3 \right]. \nonumber \end{align} \tag{ 2.55b }$$

Equation (2.55b) states that the material is shared between the droplets of radius R_i and the dilute phase of concentration $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInfty}{c_\infty} \cInfty(t)$. In order to neglect the spatial correlations between the droplets [132, 133], we assumed that the system volume $\newcommand{\Vsys}{V} \Vsys$ is large compared to all droplets, $\newcommand{\Vsys}{V} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\g}{\gamma} \Vsys \gg \sum_i V_i$ with $V_i = \frac{4\pi}3 R_i^3$. In this limit, equation (2.55b) can also be approximated as $\newcommand{\vetr}{\boldsymbol} \newcommand{\cbase}{c^{(0)}} \newcommand{\Vsys}{V} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\ba}{\vetr{a}} \newcommand{\cB}{{{\mathcal B}}} \bar c \Vsys \simeq \cBaseIn \sum_{i=1}^N V_i(t) + c_\infty(t) V$.

In the limit of many droplets, the system can be described by a continuous droplet size distribution. Additionally, if the supersaturation is small, Lifshitz and Slyozov [5] demonstrated that this size distribution converges to a universal form

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle P(\tilde R) &= \frac{4}{9} \, \tilde{R}^2 \biggl(1+\frac{\tilde{R}}{3} \biggr)^{-7/3} \nonumber \\&\quad \times \biggl(1-\frac{2 \tilde{R}}{3} \biggr)^{-11/3} \exp{\left(1-\frac{3}{3-2 \tilde{R}} \right)} \, , \nonumber \end{align} \tag{ 2.56 }$$

when radii are rescaled by the critical radius $\newcommand{\R}{c_{R}} \newcommand{\Rcrit}{R_{\rm c}} \Rcrit$, i.e. $\newcommand{\R}{c_{R}} \newcommand{\Rcrit}{R_{\rm c}} \tilde R=R/\Rcrit$, irrespective of the initial size distribution; see figure 2(c). In such a coarsening system where droplets grow and shrink, the critical radius scales with the average droplet radius, $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\R}{c_{R}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\Rcrit}{R_{\rm c}} \Rcrit(t) \simeq \mean{R(t)}$. Moreover, a droplet radius R_i is typically in the order of the critical radius $\newcommand{\cbase}{c^{(0)}} \newcommand{\R}{c_{R}} \newcommand{\cBase}{\cbase} \newcommand{\cIn}{c_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cInfty}{c_\infty} \newcommand{\Rcrit}{R_{\rm c}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \Rcrit(t)=\ell_\gamma/ (\frac{\cInfty(t)}{\cBaseOut}-1)$. Thus, equation (2.55) gives $\newcommand{\diff}{{{\rm d}}} \newcommand{\difffrac}[2]{\frac{\diff #1}{\diff #2}} \newcommand{\cbase}{c^{(0)}} \newcommand{\R}{c_{R}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\Rcrit}{R_{\rm c}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \difffrac{\Rcrit(t)}{t} \simeq \frac{D \ell_\gamma \cBaseOut} {\cBaseIn \Rcrit^2 }$ leading to the scaling

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\Rcrit}{R_{\rm c}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\R}{c_{R}} \newcommand{\cbase}{c^{(0)}} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\ra}{\rangle} \displaystyle \label{eq_R_t_Ostwald} \mean{R(t)} \simeq \Rcrit(t) \propto \left(\frac{D\ell_\gamma \cBaseOut}{\cBaseIn} \, t \right)^{\frac13} \, , \nonumber \end{align} \tag{ 2.57 }$$

which is the Lifshitz–Slyozov scaling law. In summary, the increasing mean droplet radius and critical radius reflect coarsening dynamics where large droplets grow at the expense of smaller ones.

Another coarsening mechanism in emulsions, besides Ostwald ripening, is the coalescence of droplets driven by their Brownian motion [116, 126]. Brownian coalescence is not included in the theory presented here, since we neglected momentum transport due to thermal fluctuations. However, the evolution of the mean droplet size due to droplet coalescence can be determined by estimating the change in radius $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \Delta \mean{R}$ for a typical fusion event and the frequency $\Delta t^{-1}$ of inter-droplet encounters. Since the droplet volume is conserved during fusion, two equally sized droplet of size $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \mean{R}$ lead to a change of the mean radius of $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\s}[1]{s^{(#1)}} \Delta \mean{R} \simeq (2^{\frac{1}{3}} - 1) \mean{R}$. The frequency of inter-droplet encounters can be estimated by the diffusion time leading to $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\la}{\langle} \newcommand{\e}{{\rm e}} \Delta t^{-1} \simeq \lambda D_R/\ell_{\rm p}^2$, where $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\kb}{k_{\rm B}} \newcommand{\etaR}{\eta_R} \newcommand{\e}{{\rm e}} D_R=\kb T/(6 \pi \etaR \mean{R})$ is the Stokes–Einstein diffusion constant of a spherical droplet with $\newcommand{\etaR}{\eta_R} \newcommand{\e}{{\rm e}} \etaR$ denoting the viscosity of the surrounding fluid experienced by the droplet of average size $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \mean{R}$. Moreover, $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\e}{{\rm e}} \ell_{\rm p}=V/(\pi N \mean{R}^2)$ is the mean free path between the droplets, where the droplet number can be estimated as $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\s}[1]{s^{(#1)}} N \simeq V_{\rm tot}/(\frac{4\pi}{3} \mean{R}^3)$ and the volume $V_{\rm tot}$ occupied by droplets is determined by particle conservation, i.e. $\newcommand{\vetr}{\boldsymbol} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqIn}{\ceq_{\rm in}} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\ba}{\vetr{a}} V_{\rm tot}/V = (\bar c - \cEqOut)/(\cEqIn-\cEqOut)$. Not every inter-droplet encounter leads to a coalescence event in particular in the presence of surfactants [134, 135]. To account for the stochastic initiation of a coalescence event we have introduced the parameter $\newcommand{\la}{\langle} \lambda \in [0,1]$ characterising the average fraction of encounters that lead to coalescence. By writing $\newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\kb}{k_{\rm B}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\etaR}{\eta_R} \newcommand{\la}{\langle} \newcommand{\e}{{\rm e}} \mean{R}^{-1} {{\rm d}}\mean{R}/{{\rm d}}t\, \simeq \,\mean{R}^{-1} \,\Delta\mean{R}/\Delta t\, \simeq \,\lambda \kb\, T V_{\rm tot}^2/(V^2 \etaR \mean{R}^3)$, skipping the numerical prefactors, we obtain a differential equation for the mean radius $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\mean}[1]{\langle #1 \rangle} \mean{R}$, where integration gives the scaling for the mean radius arising from fusion of droplets:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\ba}{\vetr{a}} \newcommand{\etaR}{\eta_R} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\cEqIn}{\ceq_{\rm in}} \newcommand{\cEq}{\ceq} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\kb}{k_{\rm B}} \newcommand{\mean}[1]{\langle #1 \rangle} \newcommand{\ra}{\rangle} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq_R_t_Coalescence} \mean{R(t)} \propto \left(\lambda \left(\frac{\bar c - \cEqOut}{\cEqIn-\cEqOut}\right)^2 \frac{\kb T}{\etaR} \, t \right)^{\frac13} \, . \nonumber \end{align} \tag{ 2.58 }$$

Remarkably, the coarsening due to droplet coalescence has the same scaling with time as the growth of droplets by Ostwald ripening described by equation (2.57).

We can use our lattice model to determine the relative contributions of the two coarsening mechanisms to the growth of droplets. In the case of Ostwald ripening, we have to estimate the molecular diffusion constant D, the capillary length $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma$ and the relative dilution of the minority phase, $c^{(0)}_{\rm out}/ c^{(0)}_{\rm in}$. We use the Stokes–Einstein relationship to express the diffusivity as $\newcommand{\kb}{k_{\rm B}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\e}{{\rm e}} D \simeq \kb T/ (6\pi \eta_{\rm m} \nu^{1/3})$, where $\newcommand{\e}{{\rm e}} \eta_{\rm m}$ denotes the fluid viscosity felt by the molecules of volume $\nu$. Moreover, from equations (2.34) and (2.46), the capillary length is $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma \simeq \frac12 \nu^{1/3} \chi^{1/2}(\chi-2){}^{3/2}$. Finally, the binodal line corresponding to our lattice model with equal-size molecules A and B (see end of section 2.1.3) can be used to estimate the relative dilution of the minority phase. It turns out that the fraction between the equilibrium concentrations is $\newcommand{\e}{{\rm e}} c^{(0)}_{\rm out}/c^{(0)}_{\rm in} \propto \exp(-\chi)$, i.e. it decreases exponentially to zero as the interaction strength $\chi$ becomes large (limit of strong phase separation), while $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} c^{(0)}_{\rm out}/c^{(0)}_{\rm in} \simeq 1-\sqrt{\left(6 (\chi-2) \right)}$ changes only weakly close to the critical point (weak phase separation). By comparing equations (2.57) to (2.58), we find that Ostwald ripening dominates coarsening if

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\ba}{\vetr{a}} \newcommand{\etaR}{\eta_R} \newcommand{\etaM}{\eta_{\rm m}} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\cEqIn}{\ceq_{\rm in}} \newcommand{\cEq}{\ceq} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\vetr}{\boldsymbol} \displaystyle \chi^{\frac12}(\chi-2)^{\frac32} \frac{c^{(0)}_{\rm out} } {c^{(0)}_{\rm in} } \left(\frac{\cEqIn-\cEqOut}{\bar c - \cEqOut} \right)^2 \frac{\etaR}{\etaM} \lambda^{-1} \gg 1 \, , \label{eqn:coarsening_ratio} \nonumber \end{align} \tag{ 2.59 }$$

where we dropped all numerical prefactors. In the simple case of a size-independent viscosity, $\newcommand{\etaM}{\eta_{\rm m}} \newcommand{\etaR}{\eta_R} \newcommand{\e}{{\rm e}} \etaM = \etaR$, our estimates from the simple binary lattice model indicate that coalescence typically dominates Oswald ripening for most interaction parameters $\chi$. In particular, the left-hand side of equation (2.59) goes to zero close to the critical point ($\chi = 2$) and in the limit of strong phase separation ($\newcommand{\ri}{{\rm I}} \chi \rightarrow \infty$). However, for intermediate $\chi$-values, Ostwald ripening could still be the dominant coarsening mechanism because coalescence events may be suppressed by surfactants ($\newcommand{\la}{\langle} \lambda \ll 1$) or because the ratio of the viscosities satisfies $\newcommand{\etaM}{\eta_{\rm m}} \newcommand{\etaR}{\eta_R} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \etaR/\etaM \gg 1$. Such different viscosities are particularly relevant for condensed phases in polymer or protein solutions, as well as droplet-like compartments in living cells. These complex, phase separated liquids can even show visco-elastic effects leading to a dramatic slowdown of the movements of large droplet-like phases [136–138]. In particular, inside cells, diffusion of very large compartments is strongly suppressed by the cytoskeleton [139], while the diffusion of small molecules may experience less hinderance. We therefore expect that Ostwald ripening relying on evaporation and condensation of small diffusing molecules is the dominant mechanism of droplet coarsening inside cells.

Here we discuss the thermodynamics of binary mixtures in the presence of external fields such as gravitation with a gravitational acceleration g, and electric or magnetic, position-dependent potentials denoted as $U(x)$. The presence of such an inhomogeneous external potential can influence the shape and the mean position x_i of the concentration profiles $\newcommand{\A}{c_{A}} \A (x)$ and $\newcommand{\B}{c_{B}} \B (x)$, which can be defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ba}{\vetr{a}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle x_{i}= \frac{1}{L} \int_0^{L} {{\rm d}}x \, x \, \frac{c_i(x)}{\bar c_i } \, , \nonumber \end{align} \tag{ 3.1 }$$

where $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar c_i = L^{-1} \int_0^{L} {{\rm d}}x \, c_i(x)$ denotes the mean concentration and L is the size of the system.

For a compressible system, the binary mixture is described by two concentration fields $\newcommand{\A}{c_{A}} \A$ and $\newcommand{\B}{c_{B}} \B$. Considering the case where the external fields vary along the x-coordinate, the total free energy density reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\fetot}{f_{\rm tot}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\B}{c_{B}} \newcommand{\A}{c_{A}} \displaystyle \nonumber \fetot&= \fe\left(\A, \B, \nabla \A, \nabla \B \right)\nonumber \\ \label{eq:free_energy_density} & \quad + \rho g x + U_A(x) \A + U_B(x) \B \, . \nonumber \end{align} \tag{ 3.2 }$$

The interactions between the components are governed by the free energy density f , $\newcommand{\A}{c_{A}} \newcommand{\B}{c_{B}} \rho=m_A \A+ m_B \B$ is the mass density, and m_A and m_B denote the molecular mass of each component. The contributions of the external potentials can be combined to $\newcommand{\A}{c_{A}} \newcommand{\B}{c_{B}} \tilde{U}_A(x) \A + \tilde{U}_B(x) \B$, where $\tilde{U}_i(x)= U_i(x)+ m_i g x$, $i=A,B$.

Thermodynamic equilibrium for systems with external fields can be defined at each position. The position-dependent equilibrium profile c_i(x) is then determined by a spatially constant generalised chemical potential, ${\mu}_{{\rm tot},i}=\mu_i(x) + \tilde{U}_i(x)$, where $\mu_{i}(x)=\delta F/\delta c_i$ with $F=\int {{\rm d}}^3 x f$. As an example of a system with such a position-dependent equilibrium profile we consider an incompressible system ($\newcommand{\A}{c_{A}} \newcommand{\B}{c_{B}} \nu_A \A+ \nu_B \B=1$) with gravitation as the only external potential and where the A-molecules are dilute, $\newcommand{\A}{c_{A}} \nu_A \A \ll 1$ [140]. The generalised exchange chemical potential then reads $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\ba}{\vetr{a}} \bar{\mu}_{\rm tot} \simeq k_{\rm B}T \ln (\nu \A) \,+$ $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\ba}{\vetr{a}} \nu_A \Delta \rho g x$ with the density difference $\Delta \rho = m_A/\nu_A-m_B/\nu_B$. At thermodynamic equilibrium, this gives the barometric height formula, $\newcommand{\A}{c_{A}} \newcommand{\e}{{\rm e}} \A(x)\propto \exp({-\Delta \rho \nu_A g x / (k_{\rm B}T) })$. Thus, gravitation always positions the molecules of highest mass density toward the lower gravitational potential. Typically, sedimentation of macromolecules on a mesosscopic length scale L is negligible because $\Delta \rho \nu_A g L \ll k_{\rm B}T$. Soft condensed phases, however, are much more prone to follow the gravitational field because the volume of the phase $V_i$ typically exceeds the molecular volumes. The corresponding condition for sedimentation to be relevant reads $\Delta \rho_i V_i g L > k_{\rm B}T$, where $\Delta \rho_i$ denotes the mass density difference between the liquid condensed phase and the surrounding solvent. Note that a micron-sized protein droplet with about twenty percent larger mass density relative to the solvent can already show sedimentation on length scales slightly above the droplet size [143].

Gravitation acts on the condensed phase and positions it into the region of lower gravitational potential. In the following sections we will discuss another positioning mechanism that relies on a chemical coupling. In this case the position of the condensed phase is affected via a gradient of macromolecules that in turn influences phase separation through molecular interactions. For a fixed gradient, this system may already change the position of the condensed phase by affecting the molecular interactions.

To explore the propensity to switch the position of condensed phases using such an additional component we propose a simple ternary model [144]. This model accounts for the demixing of two components, A and B, and a regulator component R that interacts with the other components and thereby affects phase separation between A and B. The regulator component R is influenced by an external potential $U(x)$. Interactions between the components i and j  are captured by the mean-field interaction parameters $\chi_{ij}$. The scenario of regulation of phase separation can be described by the following free energy density

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\fetot}{f_{\rm tot}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \newcommand{\R}{c_{R}} \newcommand{\B}{c_{B}} \newcommand{\A}{c_{A}} \newcommand{\kb}{k_{\rm B}} \displaystyle \nonumber \fetot &=\kb T \bigg[ \sum_{i={{\rm A,B,R}}} c_i \ln (\nu c_i) + \chi_{AB} \, \nu\,\A \B \nonumber \\ \nonumber &\quad + \R (\chi_{BR} \, \B+\chi_{AR} \A) \, \nu \bigg] + U(x) \R \nonumber \\ \label{eq:free_energy_ternary_functional_theory_hom} & \quad + \frac{\kappa_R}{2}\left| \nabla \R \right|^2+\frac{\kappa_A}{2}\left| \nabla \A \right|^2 \, , \nonumber \end{align} \tag{ 3.3 }$$

which is a Flory–Huggins free energy density [2, 3] for three components. Analogously to the case of the binary system discussed in section 2.1, the ternary free energy given above can be derived from a partition sum using a mean-field approximation (see [145, 146] and appendix of [144]). In equation (3.3) we consider an incompressible system where the molecular volumes are constant and equal to $\nu$ for all components, and the concentrations thus obey $\newcommand{\A}{c_{A}} \newcommand{\B}{c_{B}} \newcommand{\R}{c_{R}} \B=\nu^{-1}-\R-\A$. The logarithmic terms in equation (3.3) correspond to entropic contributions related to the number of possible configurations. The remaining contributions characterise the interactions between the three components with the (dimensionless) mean field interaction parameter $\chi_{ij}$, also referred to as Flory–Huggins interaction parameter. The interaction parameter between A and B, $\chi_{{AB}}$, determines the tendency of A and B to phase separate. The two terms in equation (3.3) proportional to the regulator concentration $\newcommand{\R}{c_{R}} \R$ describe the interactions between the regulator R and the demixing components A and B. To ensure that R acts as a regulator we choose these interaction parameters such that the regulator R does not demix from A or B. The terms in equation (3.3) with spatial derivatives represent contributions to the free energy associated with spatial inhomogeneities. We neglected a mixed term proportional to $\newcommand{\A}{c_{A}} \newcommand{\R}{c_{R}} \nabla \A \cdot \nabla \R$ since it has only little quantitative impact on the spatial profiles of the phase separated profiles [144].

In the following we consider a periodic system with a periodic potential $U(x)$ that varies solely along the x-coordinate. We choose a potential that affects the distribution of the regulator component of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\kb}{k_{\rm B}} \displaystyle \label{eq:potential} U(x)= - \kb T \ln \left(1 - Q \sin{\left(2 \pi x/L \right)}\right) \, , \nonumber \end{align} \tag{ 3.4 }$$

where 0  <  Q  <  1 characterises the strength of the potential and L denotes the size of the system along the x-direction.

In the case where the components A and R are dilute, $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\ba}{\vetr{a}} \newcommand{\Ahom}{\bar {c}_A} \nu \Ahom \ll 1$ and $\newcommand{\vetr}{\boldsymbol} \newcommand{\R}{c_{R}} \newcommand{\ba}{\vetr{a}} \newcommand{\Rhom}{\bar {c}_R} \nu \Rhom \ll 1$, and for weak external potentials ($\kappa_R (Q/L){}^2 \ll 1$) such that the the gradient terms in free energy can be neglected, the profile of the regulator component is solely given by the external potential $U(x)$ with the regulator profile $\newcommand{\R}{c_{R}} \R(x)$ assuming the shape of negative sine function, $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} -\sin{\left(2 \pi x/L \right)}$. Thus, the regulator profile has one minimum and one maximum in the periodic domain. The interaction of such a regulator profile with the components A and B causes a positional dependence of their concentration profiles. We would like to understand how these interactions affect the system if A and B phase separate.

For simplicity, we also consider a one-dimensional system of size L in the absence of boundaries. In this one-dimensional system the periodic boundary conditions are $c_i (0)= c_i (L)$ and $c_i' (0)= c_i' (L)$, where the primes denote spatial derivatives. For the considered case of an external potential $U(x)$ varying only along the x-coordinate, the restriction to a one dimensional, phase separating system represents a valid approximation for large system sizes, where the interface between the condensed phases becomes flat.

To find the equilibrium states in a phase separating system in the presence of a regulator gradient induced by the external potential $U(x)$, we determine the concentration profiles c_i(x) of all components $i=A,B,R$ by minimising the total free energy (equation (3.3); see [144]). Due to particle number conservation, there are two constraints for the minimisation: $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar c_i = L^{-1} \int_0^L {{\rm d}}x \, c_i(x)$ for $i=A,R$, where $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar c_i$ are the average concentrations and $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\R}{c_{R}} \newcommand{\Rhom}{\bar {c}_R} \newcommand{\Ahom}{\bar {c}_A} \newcommand{\ba}{\vetr{a}} \bar c_B =\nu^{-1}-\Ahom-\Rhom$. Variation of the total free energy with the constraints of particle number conservation implies ($i=A,R$):

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\fetot}{f_{\rm tot}} \newcommand{\fe}{f} \newcommand{\be}{\beta} \displaystyle \label{eq:minimization_2} 0= \int_{0}^{L}{{\rm d}}x \, \left(\frac{\partial \fetot}{\partial c_i} - \frac{{{\rm d}}}{{{\rm d}}x} \frac{\partial{\fetot}}{\partial c_i'} + \lambda_i \right) \delta c_i + \kappa_i \delta c_i c_i' \bigg{|}_0^L\, , \nonumber \end{align} \tag{ 3.5 }$$

where $\newcommand{\la}{\langle} \lambda_R$ and $\newcommand{\la}{\langle} \lambda_A$ are Lagrange multipliers and the $\delta c_i$ is the variation of the concentration corresponding to component i. The boundary term in equation (3.5) is zero in the case of periodic boundary conditions. Using the explicit form of the free energy density (equation (3.3)), a set of Euler–Lagrange equations can be derived from equation (3.5) (see [144]). These equations can be solved numerically using a finite difference solver but also approximately investigated analytically (see section 3.1.5). As control parameters we consider the three interaction parameters $\chi_{AR}$, $\chi_{AB}$ and $\chi_{BR}$, the strength of the external potential Q and the mean concentration of A-material, $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\ba}{\vetr{a}} \newcommand{\Ahom}{\bar {c}_A} \Ahom$. The mean concentration of the regulator material is fixed to a small volume fraction $\newcommand{\vetr}{\boldsymbol} \newcommand{\R}{c_{R}} \newcommand{\ba}{\vetr{a}} \newcommand{\Rhom}{\bar {c}_R} \nu \Rhom =0.02$ in all presented studies to avoid phase separation of the regulator. Moreover, we focus on the limit of strong phase segregation, where the interfacial width determined by $\kappa_A$ are small compared to the system size. We verified that our results depend only weakly on the specific values of $\kappa_A$ and $\kappa_R$.

Solving the Euler–Lagrange equations, we find two extremal profiles of the phase separating component A, $\newcommand{\A}{c_{A}} \newcommand{\Al}{c_{A}^{-}} \Al(x)$ and $\newcommand{\A}{c_{A}} \newcommand{\Ar}{c_{A}^{+}} \Ar(x)$, and two corresponding profiles of the regulator component R, denoted as $\newcommand{\R}{c_{R}} \newcommand{\Rl}{c_{R}^{-}} \Rl(x)$ and $\newcommand{\R}{c_{R}} \newcommand{\Rr}{c_{R}^{+}} \Rr(x)$ (the profile of B follows from number conservation). The phase separating material A can be accumulated at larger regulator concentration and correlates (+) with the concentration of the regulator material (figure 3(a)). The corresponding solutions are $\newcommand{\A}{c_{A}} \newcommand{\Ar}{c_{A}^{+}} \Ar(x)$ and $\newcommand{\R}{c_{R}} \newcommand{\Rr}{c_{R}^{+}} \Rr(x)$. Alternatively, the A-material accumulates at smaller regulator concentrations ($\newcommand{\A}{c_{A}} \newcommand{\Al}{c_{A}^{-}} \Al(x)$ and $\newcommand{\R}{c_{R}} \newcommand{\Rl}{c_{R}^{-}} \Rl(x)$) corresponding to an anti-correlation (−) with respect to the regulator profile (figure 3(b)). The free energies of the correlated and the anti-correlated states, $\newcommand{\A}{c_{A}} \newcommand{\R}{c_{R}} \newcommand{\Rr}{c_{R}^{+}} \newcommand{\Ar}{c_{A}^{+}} \newcommand{\Fr}{F^{+}} \Fr=F[\Ar, \Rr]$ and $\newcommand{\A}{c_{A}} \newcommand{\R}{c_{R}} \newcommand{\Rl}{c_{R}^{-}} \newcommand{\Al}{c_{A}^{-}} \newcommand{\Fl}{F^{-}} \Fl=F[\Al, \Rl]$, are different for most interaction parameters. The free energies only intersect at one point $\chi_{BR}=\chi_{BR}^*$ (figure 4(a)). At this point the minimal free energy exhibits a kink. This means that the system undergoes a discontinuous phase transition when switching between a spatially anti-correlated (−) and a spatially correlated (+) solution with respect to the regulator.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To study this phase transition, the appropriate set of order parameters can be defined from the changes of the free energy upon varying the interaction parameters (see figure 4):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \newcommand{\kb}{k_{\rm B}} \displaystyle \label{eq:order_parameter2} \rho_{ij} = \left(\kb T \mathcal{N}_{ij} \nu L \right)^{-1} \frac{{{\rm d}}}{{{\rm d}}\chi_{ij}} \Delta F(c_i(x),c_j(x)) \, , \nonumber \end{align} \tag{ 3.6 }$$

where $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \Delta F(c_i(x),c_j(x))=F(c_i(x),c_j(x)) - F(\bar c_i, \bar c_j)$. The normalisation $\mathcal{N}_{ij}$ is chosen such that $-1<\rho_{ij}<1$. When inserting equation (3.3), the order parameter $\rho_{ij}$ becomes the covariance,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\ba}{\vetr{a}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \label{eq:order_parameter} \rho_{ij} =\left(\mathcal{N}_{ij} L\right)^{-1} \int_0^L {{\rm d}}x \, \left(c_i(x) c_j(x) - \bar c_i \bar c_j\right) \, , \nonumber \end{align} \tag{ 3.7 }$$

which characterises the spatial correlation between the concentration profiles c_i(x) and c_j(x). If the fields are spatially correlated (+), $\rho_{ij}>0$, and if they are anti-correlated (−), $\rho_{ij}<0$, and $\rho_{ij}=\pm1$ if the concentration profiles of component i and j  follow spatially correlated or anti-correlated step functions. If the regulator is homogeneous, $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} c_R(x)=\bar c_R$, the order parameter is zero, $\rho_{iR}=0$ for $i=A,B$.

Varying the interaction parameter $\chi_{BR}$ (figure 4(b)), the order parameters $\rho_{BR}$ and $\rho_{AR}$ jump at the threshold value $\chi_{BR}^*$. The jump of both order parameters in the presence of a regulator gradient indicates that the spatial correlation of A and B with respect to R changes abruptly, which is expected in case of a first order phase transition.

By means of the order parameter $\rho_{BR}$ (equation (3.7)) we can now discuss the phase diagrams as a function of the interaction parameters. In the case of a spatial correlation (+), we have $\rho_{ij}>0$, while for an anti-correlation (−), $\rho_{ij}<0$. We thus find three regions (figure 4(c)): a mixed region, where concentration profiles are only weakly inhomogeneous, and no phase separation occurs, and two additional regions, where components A and B phase separate and A is spatially correlated or anti-correlated with the regulator R, respectively. There exists a triple point where all three states have the same free energy.

In summary, the presence of a concentration gradient in phase separating systems leads to equilibrium states of different spatial correlation with the regulator profile. The regulator gradient creates a bias in the position of the phase separated concentration profiles for almost all parameters in the phase diagram. If the external potential acting on the regulator has exactly one minimum and one maximum, there are two stationary states with different mean positions of the phase separating material. One of these stationary states corresponds to a global minimum of the free energy while the other state may only be locally stable. The parameters characterising the interactions between the regulator and the phase separating material determine which of these states corresponds to equilibrium. There is a discontinuous phase transition between both states upon changing these interaction parameters.

For simplicity we have discussed a phase separating system in the presence of an external potential restricting to inhomogeneities in one dimension. However, preliminary Monte-Carlo studies in three dimensions with phase separated droplets in a regulator gradient suggest that the position of droplets can be switched in a discontinuous manner.

In the previous section we have considered the numerical minimisation of a set of non-linear Euler–Lagrange equations derived from the free energy density (equation (3.3)). Here we give some approximate analytic arguments to understand the minimal ingredients for the occurrence of the discontinuous phase transition. To this end, we would like to simplify the system further and consider the dilute limit of the regulator, i.e. $\newcommand{\vetr}{\boldsymbol} \newcommand{\R}{c_{R}} \newcommand{\ba}{\vetr{a}} \newcommand{\Rhom}{\bar {c}_R} \nu \Rhom \ll 1$ and thus approximate $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\Ahom}{\bar {c}_A} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\ba}{\vetr{a}} \bar c_B \simeq 1 - \nu^{-1} \Ahom$ in equation (3.3). For such dilute conditions, the equilibrium concentrations (for A component) in each phase, $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} \cBaseIn$ and $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \cBaseOut$, are then well described by the binary A-B system. In addition, for strong enough external potential $U(x)$, the regulator profile is well approximated by, $\newcommand{\R}{c_{R}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} \R(x)= \tilde A - \tilde B \sin{\left(2 \pi x/L \right)}$, apart from the peaks (see figure 3), where $\tilde A>0$ is some concentration offset and $\tilde B>0$ characterises the strength of the spatial modulations of the regulator profile. For a discussion of the relevance of the regulator peaks at the interface see [144]; here we simply neglect these peaks for simplicity. We expect that the extrema of the regulator profile $\newcommand{\R}{c_{R}} \R(x)$ at $x=L/4, 3L/4,$ determine the position of the A-rich condensed phase. As obtained from the numerical analysis presented in the last section, there are two solutions, either spatially correlated (+) or anti-correlated (−) with the regulator. These solutions can be approximated in the dilute limit of the regulator as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\A}{c_{A}} \newcommand{\cbase}{c^{(0)}} \displaystyle \nonumber \A^+(x) &\simeq \cBaseOut + \left(\cBaseIn- \cBaseOut \right) \Theta \left(x- 3L/4 + x_0/2 \right) \nonumber \\ & \quad \times \Theta \left(3L/4 + x_0/2 -x \right) \, , \nonumber \end{align} \tag{ 3.8a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\A}{c_{A}} \newcommand{\cbase}{c^{(0)}} \displaystyle \nonumber \A^-(x) &\simeq \cBaseOut + \left(\cBaseIn- \cBaseOut \right) \Theta \left(x- L/4 + x_0/2 \right) \nonumber \\ & \quad \times \Theta \left(L/4 + x_0/2 -x \right) \, , \nonumber \end{align} \tag{ 3.8b }$$

where $\Theta(\cdot)$ denotes the Heaviside step function. These solutions describe the A-rich domains either localised around the maximal (+) or minimal (−) amount of regulator. The domain size of the A-rich phase denoted as x₀ is determined by conservation of particles, i.e. $\newcommand{\vetr}{\boldsymbol} \newcommand{\A}{c_{A}} \newcommand{\ba}{\vetr{a}} \newcommand{\Ahom}{\bar {c}_A} \Ahom L= \int {{\rm d}} x \A^+(x) = \int {{\rm d}} x \A^-(x)$, leading to $\newcommand{\vetr}{\boldsymbol} \newcommand{\cbase}{c^{(0)}} \newcommand{\A}{c_{A}} \newcommand{\ba}{\vetr{a}} \newcommand{\Ahom}{\bar {c}_A} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} x_0=L (\Ahom-\cBaseOut)/(\cBaseIn-\cBaseOut)$. Using the approximate solutions above we can calculate the difference in free energy corresponding to correlated and anti-correlated states,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\Fl}{F^{-}} \newcommand{\Fr}{F^{+}} \newcommand{\kb}{k_{\rm B}} \displaystyle \label{eq_deltaF_analytic} \Fr - \Fl \simeq W \,{\kb T \, L} \left(\chi_{AR} - \chi_{BR} \right) \tilde B \, , \nonumber \end{align} \tag{ 3.9 }$$

where $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} W= 2 \pi^{-1} \sin \left(\pi \frac{x_0}{L}\right)>0$ is a positive constant. Note that all contributions apart from the A-R and B-R interaction vanish in the dilute limit and due to conservation of A and B material. Since the parameter $\tilde B$ characterises the concentration modulations of the regulator profile, the free energy difference $\newcommand{\Fr}{F^{+}} \newcommand{\Fl}{F^{-}} (\Fr - \Fl)$ consistently vanishes for zero $\tilde B$ (equation (3.9)). In the case of non-zero $\tilde B$, the free energy difference is determined by the difference in the interactions of A and B with respect to the regulator R, $\Delta \chi =\chi_{AR} - \chi_{BR}$. Most importantly, at $\Delta \chi=0$, the two solutions switch their thermodynamic stability: for $\Delta \chi>0$, the anti-correlated state is favoured, while for $\Delta \chi<0$, the correlated state is preferred; $\Delta \chi=0$ indicates the transition point. In addition, the slopes of correlated and anti-correlated free energy, $\newcommand{\Fr}{F^{+}} \Fr$ and $\newcommand{\Fl}{F^{-}} \Fl$, with respect to $\Delta \chi$ at the transition point are not equal. The difference in slopes implies that the minimal free energy exhibits a kink at the transition point $\Delta \chi=0$, which means that the system undergoes a discontinuous phase transition switching from a correlated to an anti-correlated state for increasing $\Delta \chi$. The discontinuous phase transition occurs at $\Delta \chi=0$, which agrees with the numerical predictions shown in figure 4(c) (note that the corresponding volume fraction of the regulator is rather dilute). Our approximate analytic treatment indicates that the occurrence of a discontinuous phase transition solely relies on the existence of the position-dependent profile of the regulator and the interactions of the regulator molecules with the liquid condensed phases. It seems that the peaks of regulator material at the interface of the condensed phase, which we neglected in this approximate analytic argument, are not necessary to observe the discontinuous transition. The discontinuous nature may also be preserved if the regulator gradient is driven by boundary conditions? We leave this question to future research.

Here, we discuss the ripening dynamics of droplets in a regulator gradient that varies only along the x-coordinate. To this end, we modify the concept of a common far field concentration introduced in section 2.5.3 to a concentration field $\newcommand{\cinf}{{c}_\infty} \cinf(x)$ that changes on the length scale of the system size L.

In the absence of a regulator gradient, the concentration outside approaches the 'far field' of the droplet, $\newcommand{\cinf}{{c}_\infty} \cinf$, as the distance to the droplet interface increases. The far field is created by the surrounding droplets and is well reached if the length scale $\newcommand{\e}{{\rm e}} \ell$ corresponding to the mean distance between droplets exceeds the droplet radius R, i.e. $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell \gg R$.

In the presence of a regulator gradient varying along the x-coordinate, the far field $\newcommand{\cinf}{{c}_\infty} \cinf(x)$ seen by the droplet is now also position dependent and can be approximately written as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\cOut}{c} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\cinf}{{c}_\infty} \displaystyle \cinf(x) \simeq \frac{1}{L_y L_z } \int_{0}^{L_y} {{\rm d}}y \int_{0}^{L_z} {{\rm d}}z \, \cOut_{\rm out}(x,y,z) \, , \nonumber \end{align} \tag{ 3.10 }$$

where $\newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \cOut_{\rm out}(x,y,z)$ is the concentration field outside the droplets and L_y and L_z denote the system size in the y  and z-direction, respectively. The expression above is an approximation because we have neglected weak concentration perturbations close to droplet interfaces described by the Gibbs Thomson relationship (2.45). However, for the typical case of droplet radii exceeding the capillary length ($\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} R \gg \ell_\gamma$) and the mean inter-droplet distance being larger than the droplet size ($\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell \gg R$), these concentration perturbations are very small (see section 2.5.1).

Finally, to make sure that the spatial variations of the position dependent far field are large on the system size but comparably small on the droplet scale, we consider the case where all these length scales separate, $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma \ll R \ll \ell \ll L$. This separation of length scales will allow us to investigate weak perturbations of the droplet shape parallel to the concentration gradient and also to construct the equilibrium concentration at each position x along the gradient.

The separation of length scales suggests dividing the system into independent slices of a size corresponding to the intermediate length scale $\newcommand{\e}{{\rm e}} \ell$. The phase separation dynamics can then be discussed locally for each position x corresponding to a slice element of linear length $\newcommand{\e}{{\rm e}} \ell$. For this discussion, we consider a simplified model with a free energy density given by equation (3.3) and calculate the corresponding phase diagram (figure 5(a)). If the droplet material is roughly constant, each position x maps on a single point in the phase diagram because the regulator profile is fixed. The corresponding curve in the phase diagram due to the spatial dependence of the regulator defines local values of the equilibrium concentrations inside and outside of the droplet (figure 5(b)). In other words, droplets in the slice corresponding to the position x feel the local equilibrium concentrations, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \ceqbaseIn(x)$ and $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$. For simplicity, we restrict ourselves to a special case where the equilibrium concentration inside is position independent, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} \ceqbaseIn(x)\simeq \ceqbaseIn$, so droplet growth is solely determined by the conditions outside of the droplet. Moreover, as concentration inhomogeneities of droplet material are generally weak, we do not consider weak transients of $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$ due to the space and time varying $\newcommand{\cinf}{{c}_\infty} \cinf(x)$. The actual concentration of droplet material outside, $\newcommand{\cinf}{{c}_\infty} \cinf(x)$, together with the equilibrium concentration outside, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$, determine a spatially dependent supersaturation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\cinf}{{c}_\infty} \displaystyle \label{eq:supersaturation} \varepsilon(x)=\frac{\cinf(x)}{\ceqbase(x)}-1 \, . \nonumber \end{align} \tag{ 3.11 }$$

Zoom In Zoom Out Reset image size

There is a dissolution boundary located at the position $\newcommand{\xd}{x_{\rm c}} x=\xd$ where the supersaturation $\newcommand{\xd}{x_{\rm c}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \varepsilon(\xd)=\ell_\gamma/R$. For the considered case of $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma \ll R$, this boundary approximately corresponds to a vanishing supersaturation $\newcommand{\xd}{x_{\rm c}} \newcommand{\s}[1]{s^{(#1)}} \varepsilon(\xd)\simeq 0$. In the illustration in figure 5(a), the fluid is mixed for $\newcommand{\xd}{x_{\rm c}} x<\xd$, while droplets can form ($\varepsilon>0$) for $\newcommand{\xd}{x_{\rm c}} x>\xd$. In the absence of droplets, the concentration field $\newcommand{\cinf}{{c}_\infty} \cinf(x)$ evolves in time satisfying a diffusion equation, which we will derive in the next section. When droplets are nucleated, their local dynamics of growth or shrinkage is guided by the local supersaturation $\varepsilon(x)$ as well as $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$ (see section 2.5.3). This local droplet dynamics then in turn also influences the concentration field $\newcommand{\cinf}{{c}_\infty} \cinf(x)$. As time proceeds, diffusion of droplet material occurs on length scales larger than the intermediate length scale $\newcommand{\e}{{\rm e}} \ell$. For this regime, we will derive a dynamical theory and extend the Lifschitz & Slyozov theory to concentration gradients.

A regulator concentration gradient generates a position-dependent equilibrium concentration $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$ and a position-dependent supersaturation $\varepsilon(x)$ (equation (3.11)). This supersaturation will drive the droplet dynamics and lead to a position-dependent growth, drift of droplets and even deformations of their shape. In the following, we discuss the dynamics of growth of a single droplet where the equilibrium concentration, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$, and the concentration of droplet material, $\newcommand{\cinf}{{c}_\infty} \cinf(x)$, are position dependent. The case of multiple droplets is studied in the next section.

The concentration inside the droplet can be approximated by the equilibrium concentration $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cEq}{\ceq} \newcommand{\cEqIn}{\ceq_{\rm in}} \cEqIn \simeq \ceqbaseIn$ in the limit of strong phase separation ($\newcommand{\g}{\gamma} c^{(0)}_{\rm in} \gg c^{(0)}_{\rm out}$; see section 2.5.1). This allows us to restrict the analysis to the concentration field $\newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \cOut(r,\theta,\varphi)$ outside of a droplet. Here, we use spherical coordinates centred at the droplet position x₀, with r denoting the radial distance from the centre, and $\theta$ and $\varphi$ are the azimuthal and polar angles relative to the x-axis. Within the quasi-stationary approximation (see section 2.5.2) the concentration outside but near the droplet obeys the steady state of a diffusion equation (2.50). For large r the concentration field approaches the 'far field' which in the presence of a linear gradient reads (figure 6):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\cOut}{c} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\al}{\alpha} \newcommand{\cinf}{{c}_\infty} \displaystyle \label{eq:bc1} \cinf = \lim_{r\to\infty} \cOut(r, \theta, \varphi) \simeq \al + \be \, r \cos \theta \, . \nonumber \end{align} \tag{ 3.12 }$$

Zoom In Zoom Out Reset image size

This inhomogeneous far field concentration is locally (with respect to inter-droplet distance $\newcommand{\e}{{\rm e}} \ell$) characterised by the concentration $\newcommand{\cinf}{{c}_\infty} \newcommand{\al}{\alpha} \al=\cinf(x_0)$ and the gradient $\newcommand{\cinf}{{c}_\infty} \newcommand{\be}{\beta} \be=\partial_x\cinf(x)|_{x_0}$ at the position of the droplet x₀. At the surface of the spherical droplet, r  =  R, the boundary condition is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\ri}{{\rm I}} \newcommand{\cEqOut}{\ceq_{\rm out}} \newcommand{\cEq}{\ceq} \newcommand{\cOut}{c} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \label{eq:bc2} \cOut(R,\theta, \varphi) &=\cEqOut(\theta)\nonumber \\ \nonumber &\simeq\left(\ceqbase + R \cos(\theta)\partial_x \ceqbase(x)|_{x_0} \right) \left(1+\ell_\gamma/R \right) \, . \nonumber \end{align} \tag{ 3.13 }$$

Here, $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma$ is the capillary length as introduced in section 2.5.1. Equation (3.13) corresponds to the Gibbs–Thomson relation (see section 2.5.2), which describes the increase of the local concentration at the droplet interface relative to the equilibrium concentration due to the surface tension of the droplet interface. The presence of spatial inhomogeneities on the scale of the droplet R leads to an additional term in the Gibbs–Thomson relation. To linear order, this contribution to $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqOut}{\ceq_{\rm out}} \cEqOut$ reads $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} R \cos(\theta)\partial_x \ceqbase(x)|_{x_0}$. The values of $\newcommand{\be}{\beta} \be$ and $\newcommand{\al}{\alpha} \al$ characterising the far field, $\newcommand{\cinf}{{c}_\infty} \cinf(x_0)$, together with the local equilibrium concentration at the droplet surface, $\newcommand{\ceq}{c^{\rm eq}} \newcommand{\cEq}{\ceq} \newcommand{\cEqOut}{\ceq_{\rm out}} \cEqOut(\theta)$, determine the local rates of growth or shrinkage of the drop at x  =  x₀ in a spatially inhomogeneous regulator gradient.

The solution to the Laplace equation with cylindrical symmetry is of the form $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\bi}{\boldsymbol} \cOut (r,\theta) = \sum_{i=0}^\infty \big(A_i r^i +B_i r^{-i-1}\big) P_i (\cos \theta)$, where $P_i (\cos \theta)$ are the Legendre polynomials. Using the boundary conditions (3.12) and (3.13), we find

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\ri}{{\rm I}} \newcommand{\cOut}{c} \newcommand{\be}{\beta} \newcommand{\al}{\alpha} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \label{eq:sol_conc_grad} \cOut(r, \theta)&= \al \left(1- \frac{R}{r}\right) + \be \cos \theta \left(r - \frac{R^3}{r^2} \right) \nonumber \\ & \quad + \left(\ceqbase + R \cos(\theta)\partial_x \ceqbase(x)|_{x_0} \right) \left(1+\frac{\ell_\gamma}{R}\right) \frac{R}{r} \, . \nonumber \end{align} \tag{ 3.14 }$$

The droplet could grow, drift or deform due to normal fluxes of droplet material at the interface leading to a movement of the interface. The speed normal to the interface reads $\newcommand{\vetr}{\boldsymbol} v_{\rm n}=\vetr n \cdot \vetr{v}_{\rm n}$, where $\newcommand{\vetr}{\boldsymbol} \vetr n$ denotes the normal vector to the interface. In case of a spherical droplet, $\newcommand{\vetr}{\boldsymbol} \vetr n = \vetr e_r$, where $\newcommand{\vetr}{\boldsymbol} \vetr e_r$ is the radial unit vector in spherical coordinates. In the limit of strong phase separation, i.e. $\newcommand{\g}{\gamma} c^{(0)}_{\rm in} \gg c^{(0)}_{\rm out}$, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cEq}{\ceq} \newcommand{\cEqIn}{\ceq_{\rm in}} \cEqIn \simeq \ceqbaseIn$, the velocity normal to the interface, $v_{\rm n}$ (equation (2.47)), can be expressed by $\newcommand{\vetr}{\boldsymbol} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} v_{\rm n}\simeq \vetr n \cdot (\vetr j_{\rm in}-\vetr j_{\rm out})/\ceqbaseIn$. For weak spatial variations inside and outside of the droplet, the local flux is defined as $\newcommand{\vetr}{\boldsymbol} \vetr j=-D \nabla c$. The concentration inside the droplet is approximately constant and for simplicity we consider it to be independent of the droplet position. Thus, the flux inside the droplet vanishes, $\newcommand{\vetr}{\boldsymbol} \vetr j_{\rm in}=0$, while the flux outside reads $\newcommand{\vetr}{\boldsymbol} \newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \vetr j_{\rm out}=-D \nabla \cOut|_R$.

Now we discuss how the normal speed $v_{\rm n}$ can be used to calculate the droplet growth speed $v_0$, the droplet drift velocity $v_1$, and the rate of deformations from the spherical shape, $v_2$. To this end, we parametrise the surface of the droplet in terms of Legendre polynomials as there is no dependence on the polar angle, which gives $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} \mathcal{R}(\theta, t)=\sum_i d_i(t) P_i \left(\cos\theta\right)$, where d_i(t) are the expansion coefficients characterising the shape of the interface. The corresponding interface velocity of an approximately spherical droplet is $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} v_{\rm n} \simeq\partial_t \mathcal{R}(\theta, t)=\sum_i v_i(t) P_i \left(\cos\theta\right)$, where the speeds for droplet growth, drift and deformations along the regulator gradient read $v_i(t)=\partial_t d_i(t)$. We can now identify the radius R as $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} d_0= \langle \mathcal{R}, P_0 \rangle/\langle P_0, P_0 \rangle$, the position of the droplet center x₀ as $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} d_1=\langle \mathcal{R}, P_1 \rangle/\langle P_1, P_1 \rangle$, and the deformations are characterised by $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} d_2=\langle \mathcal{R}, P_2 \rangle/\langle P_2, P_2 \rangle$ Here, the brackets indicate the scalar product $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \langle h , g \rangle = \int_0^\pi \mathrm{d}\theta \sin \theta \, h \, g$ between the functions g and h. Most importantly, we can identify $v_0= \mathrm{d}R/\mathrm{d}t$ as the rate of change of the radius and $v_1= \mathrm{d}x_0/\mathrm{d}t$ as the drift velocity of the droplet center.

Using equation (3.14), we find for the droplet growth speed

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \newcommand{\al}{\alpha} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \label{eq:LS_inhomog} \frac{\mathrm{d}R}{\mathrm{d}t }= \frac{D}{\ceqbaseIn R}\; \left[ \al -\ceqbase(x_0) \left(1 + \frac{\ell_\gamma}{R} \right)\right] \, . \nonumber \end{align} \tag{ 3.15 }$$

In the presence of concentration gradients there also exists a net droplet drift speed

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \label{eq:drift} \frac{\mathrm{d}x_0}{\mathrm{d}t} = \frac{D}{\ceqbaseIn } \left[ 3\be - \partial_x \ceqbase (x)|_{x_0} \left(1 + \frac{\ell_\gamma}{R} \right) \right] \, , \nonumber \end{align} \tag{ 3.16 }$$

where we found in contrast to [31] an additional factor of 3 in front of the coefficient $\newcommand{\be}{\beta} \beta$. Note that both the growth rate and the drift speed are proportional to the molecular diffusion constant D of droplet material.

If the far field $\newcommand{\cinf}{{c}_\infty} \cinf(x)$ is well parametrised by a linear gradient (equation (3.12)), there are no deformations from the spherical shape, $v_2=0$. Deformations of the spherical shape can only occur if higher order polynomials $P_n (\cos \theta)$ with $\newcommand{\g}{\gamma} n \geqslant 2$ are necessary to describe the far field. If we include the quadratic order in the parametrisation of the far field (equation (3.12)), $\newcommand{\cinf}{{c}_\infty} (r \cos \theta){}^2 \partial^2_x \cinf(x)|_{x_0}$, the deformation speed reads $\newcommand{\cinf}{{c}_\infty} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} v_2=(10/3) (R D/\ceqbaseIn) \partial^2_x \cinf(x)|_{x_0}$. This quadratic contribution does not affect the droplet drift $v_1$ but the growth law $v_0$ is changed. The quadratic term gives an extra contribution inside the brackets of equation (3.15) of the form $\newcommand{\cinf}{{c}_\infty} (5/3) R^2 \partial^2_x \cinf(x)|_{x_0}$. Thus, shape deformations and their impact on the growth law in the case of a non-linear far field gradient are negligible if

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\cinf}{{c}_\infty} \displaystyle \frac{5 \partial^2_x \cinf(x)|_{x_0}}{3\cinf(x)|_{x_0}} R^2 \ll 1 \, . \nonumber \end{align} \tag{ 3.17 }$$

For the system under consideration where length scales separate, i.e. $\newcommand{\e}{{\rm e}} R \ll \ell \ll L$, deformations from the spherical shape are weak because gradients of the far field $\newcommand{\cinf}{{c}_\infty} \cinf(x)$ occurs on the length scale of the system size L. In recent numerical studies considering a continuous phase separating Flory–Huggins model in a regulator gradient maintained by sink and source terms, droplet deformations are indeed visible when droplets approach the order of the system size [91].

We can now describe the dynamics of many droplets, $i=1,\dots, N$, with positions x_i and radii R_i. If droplets are far apart from each other, the rate of growth of droplet i reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \label{eq:LS_inhomogi} \frac{\mathrm{d}R_i }{\mathrm{d}t} = \frac{D}{R_i}\; \frac{\ceqbase(x_i)}{\ceqbaseIn } \left[ \varepsilon(x_i) - \frac{\ell_\gamma}{R_i}\right] \, . \nonumber \end{align} \tag{ 3.18a }$$

The drift velocity of droplet i, ${v}_1(x_i)=\mathrm{d} x_i/\mathrm{d}t$, is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\cinf}{{c}_\infty} \displaystyle \label{eq:drifti} \frac{\mathrm{d} x_i}{\mathrm{d}t}= \frac{D}{\ceqbaseIn} \left[ 3\partial_x \cinf(x)\vert_{x_i} - \partial_x \ceqbase(x)\vert_{x_i} \left(1 + \frac{\ell_\gamma}{R_i} \right) \right] \, . \nonumber \end{align} \tag{ 3.18b }$$

If the distance between droplets is large relative to their size, droplets only interact via the far field concentration field $\newcommand{\cinf}{{c}_\infty} \cinf(x,t)$. It is governed by a diffusion equation including gain and loss terms associated with growth or shrinkage of drops:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\cinf}{{c}_\infty} \displaystyle \label{eq:beta_eq} \partial_t \cinf(x,t) &= D \frac{\partial^2}{\partial x^2} \cinf(x,t) \nonumber \\ \nonumber & \quad - H(t) \sum^{N}_{i=1}\delta(x_i-x) \frac{4\pi}{3}\frac{\mathrm{d} }{\mathrm{d} t} R_i(t)^3 \, , \nonumber \end{align} \tag{ 3.18c }$$

where the time-dependent function

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\cinf}{{c}_\infty} \displaystyle H(t)= \frac{\ceqbaseIn - \cinf(t) }{V- \sum^{N}_{i=1}\delta(x_i-x) \frac{4\pi}{3} R_i(t)^3} \nonumber \end{align} \tag{ 3.19 }$$

is approximately constant, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} H\simeq \ceqbaseIn/V$, in the limit of strong phase separation $\newcommand{\cinf}{{c}_\infty} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\g}{\gamma} \ceqbaseIn \gg \cinf(t)$ and for very large inter-droplet distances corresponding to $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\g}{\gamma} V \gg \sum^{N}_{i=1}\delta(x_i-x) \frac{4\pi}{3} R_i(t){}^3$ with $V$.

Equation (3.18c) describes the effects of large scale spatial inhomogeneities on the ripening dynamics for the case of a regulator gradient varying along the x-axis. Since large scale variations of $\newcommand{\cinf}{{c}_\infty} \cinf(x,t)$ only build up along the x-directions, derivatives of $\newcommand{\cinf}{{c}_\infty} \cinf$ along the y  and z directions do not contribute as $\newcommand{\cinf}{{c}_\infty} \cinf$ and $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase$ are constant along these directions.

In the absence of a regulator gradient, $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase$ and $\newcommand{\cinf}{{c}_\infty} \cinf$ are constant in space implying a position-independent and common supersaturation level $\varepsilon$ for all droplets (equation (3.11)). In this case equation (3.18a) gives the classical law of droplet ripening derived by Lifschitz–Slyozov [5, 6] (also referred to as Ostwald ripening) and the net drift vanishes (equation (3.18b)). In the case of Ostwald ripening, droplets larger than the critical radius $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} R_{\mathrm{c}}= \ell_{\gamma}/\varepsilon$ grow at the expense of smaller, shrinking drops, which then disappear. This competition between smaller and larger drops causes an increase of the average droplet size and a broadening of the droplet size distribution with time. Ostwald ripening is characterised by a supersaturation that decreases with time, leading to an increase of the critical droplet radius $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} R_{\rm c}=\ell_\gamma/\varepsilon(t)\propto t^{1/3}$. The droplet size distribution $P(R)$ exhibits a universal shape and is nonzero only in the interval $[0, 3R_{\rm c}/2]$ (figure 7(b), blue graph). In other words, there are no droplets larger than $3R_{\rm c}/2$ and thus also no droplets exist beyond the maximum of $\mathrm{d}R/\mathrm{d}t$ at $R=2R_{\rm c}$ where larger droplet could grow slower. Therefore, in homogeneous systems the broadening of $P(R)$ follows from larger droplets growing at a larger rate $\mathrm{d}R/\mathrm{d}t$ than smaller droplets and, because all droplets feel the same supersaturation level, droplets remain homogeneously distributed in the system. In the presence of a regulator gradient the droplet dynamics exhibits a different behaviour.

Zoom In Zoom Out Reset image size

There are two novel possibilities of how droplet material is transported: there is exchange of material between droplets at different positions along the concentration gradient due to a position dependent droplet growth, and droplets can drift along the concentration gradient.

Droplets grow or shrink with rates that vary along the gradient because the local equilibrium concentration $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$ and the far field concentration $\newcommand{\cinf}{{c}_\infty} \cinf(x)$ are position dependent (equation (3.18a)). For a supersaturation $\newcommand{\cinf}{{c}_\infty} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ri}{{\rm I}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \varepsilon(x)=\left(\cinf(x)/\ceqbase(x)-1\right) > \ell_\gamma/R$, a droplet located at position x grows, and shrinks in the opposite case. In other words, the critical droplet radius depends on position, and droplets with radii below or above $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} R_{\mathrm{c}}(x)= \ell_\gamma/\varepsilon(x)$ shrink or grow, respectively. This position dependence implies a movement of the dissolution boundary $x_{\rm c}(t)$ which is defined by $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \varepsilon(x_{\rm c}(t))=\ell_\gamma/R$ (equation (3.11)). This definition can be simplified for the case where the capillary length, which is typically in the order of the molecular size, is small relative to the droplet radii, i.e. $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma \ll R$, leading to $\newcommand{\cinf}{{c}_\infty} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} \cinf(x_{\rm c}(t)) \simeq \ceqbase(x_{\rm c}(t))$. Taking the derivative in time gives the speed of the dissolution boundary, $v_{\rm c} (t)= \frac{\mathrm{d} x_{\rm c}(t)}{\mathrm{d}t}$:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\cinf}{{c}_\infty} \displaystyle \label{eq_diss_boundary} v_{\rm c} (t) = {\frac{\mathrm{d} \cinf (x)}{\mathrm{d}t}}\bigg{|}_{x=x_{\rm c}(t)} \bigg{/}{\frac{\mathrm{d} \ceqbase(x)}{\mathrm{d}x}}\bigg{|}_{x=x_{\rm c}(t)} \, . \nonumber \end{align} \tag{ 3.20 }$$

We can now discuss the movement direction of the dissolution boundary. If droplets can grow in the system, the far field concentration should decay, $\newcommand{\cinf}{{c}_\infty} \frac{\mathrm{d} \cinf}{\mathrm{d}t}<0$. Thus, the dissolution boundary always moves toward positions corresponding to smaller values of $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase$. When the dissolution boundary moves through the system it dissolves all droplets on its way. The dissolved droplet material will diffuse and feed the growth of the remaining droplets. Thus, the moving dissolution boundary positions the phase separated material toward one boundary of the system corresponding to the lower values of the position-dependent equilibrium concentration $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$.

Another mechanism of positioning is via droplet drift (equation (3.18b)). The drift of a droplet results from an asymmetry of material flux at the interface parallel to the regulator gradient. To be more specific, we have to distinguish between the scenario of many droplets and the case of a single droplet. In the case of many droplets, the presence of these droplets keeps the position-dependent supersaturation small, thus $\newcommand{\cinf}{{c}_\infty} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\s}[1]{s^{(#1)}} \frac{\mathrm{d} \cinf}{\mathrm{d}x} \simeq \frac{\mathrm{d} \ceqbase}{\mathrm{d}x}$. Using the derived equation (3.16) for droplet drift, the drift of droplet 'd' in a system with many droplets roughly follows

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \label{eq:drift_many} \frac{\mathrm{d} x_{\rm d}}{\mathrm{d}t} \simeq 2 \frac{D}{\ceqbaseIn} \frac{\mathrm{d} \ceqbase}{\mathrm{d}x} . \nonumber \end{align} \tag{ 3.21 }$$

This equation implies that droplets drift into the opposite direction of the dissolution boundary. Thus, droplets are pulled toward dissolution making it impossible for them to escape their dissolution via drift. This picture can be different for a single droplet. When diffusion of droplet material is fast, the far field concentration is expected to be roughly homogeneous, $\newcommand{\cinf}{{c}_\infty} \newcommand{\s}[1]{s^{(#1)}} \frac{\mathrm{d} \cinf}{\mathrm{d}x} \simeq 0$. As a consequence, the droplet drifts in the same direction as the dissolution boundary in the case of single droplet:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbaseIn}{c^{(0)}_{\rm in}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \displaystyle \frac{\mathrm{d} x_0}{\mathrm{d}t}\simeq - \frac{D}{\ceqbaseIn} \frac{\mathrm{d} \ceqbase}{\mathrm{d}x} \, . \nonumber \end{align} \tag{ 3.22 }$$

The droplet may thereby escape its dissolution by drift. The drift of a single droplet in a roughly homogeneous far-field concentration should be driven by the efflux of material at the back, which diffuses to the front of the droplet (droplet front is faced into the direction of movement of the dissolution boundary). Droplet movement then arises from asymmetric droplet growth between the back and the front of the droplet. The associated time-scale of this process is roughly given by the time to diffuse along the droplet radii, i.e. R²/D. If this time-scale is smaller than the time-scale $R/v_{\rm c}$ necessary for the dissolution boundary to move the distance R, a single droplet may escape the dissolution boundary. Using equation (3.20) the condition for a single droplet to drift can thus be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\ceq}{c^{\rm eq}} \newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\cinf}{{c}_\infty} \displaystyle D \frac{\mathrm{d} \ceqbase}{\mathrm{d}x} < R \frac{\mathrm{d} \cinf}{\mathrm{d}t} \, . \nonumber \end{align} \tag{ 3.23 }$$

Otherwise, the droplet would dissolve and typically recondense via nucleation in domains corresponding to lower values of the position-dependent equilibrium concentration $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)$.

Numerically solving equation (3.18) for a large number of droplets we find that the droplet size distribution narrows during the positioning of droplets toward one boundary of the system. For details on the numerics, please refer to [31]. This narrowing of the droplet size distribution in a concentration gradient is fundamentally different from the broadening of the droplet size distributions during classical Ostwald ripening [5, 6]; see figure 7 for an illustration of the mechanism underlying the narrowing. Imagine we spatially quench the system by imposing a spatially varying equilibrium concentration $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(x)= \ceqbase(0) \, (1- m \, x)$, where $\newcommand{\ceqbase}{c^{(0)}_{\rm out}} \newcommand{\ceq}{c^{\rm eq}} \ceqbase(0)$ denotes the equilibrium concentration before the quench and m  >  0 is the slope of the 'spatial quench'. Such a spatial quench reduces the critical radius at the right boundary at x  =  L from $R_{\rm c}(0)$ (critical radius before the quench) to $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} R_{\rm c}(x= L)=\ell_\gamma/\epsilon(x= L)$ (equation (3.11)). This quench also shifts the maximum of $\mathrm{d}R/\mathrm{d}t$ for droplets at x  =  L to smaller radii since the radius corresponding to the maximum occurs at $R=2 R_{\rm c}$. After the spatial quench there are many droplets with radii $R>2 R_{\rm c}(x= L)$. According to ${{\rm d}} R/{{\rm d}}t$ (figure 7(a)) these droplets grow more slowly than those around $R=2R_{\rm c}$ which leads to a narrowing of the droplet size distribution $P(R)$ at x  =  L. The critical radius $R_{\rm c} (x = L)$ remains small because dissolution of droplets at x  <  L leads to a diffusive flux toward x  =  L and thus keeps the concentration $\newcommand{\cinf}{{c}_\infty} \cinf(L)$ at increased levels. These conditions hold for a longer time if the spatial quench has a steeper slope. As a result, the distribution narrows more for steeper quenches. For weak enough slopes of the quench, the narrowing of the droplet size distribution vanishes, but droplet positioning still occurs as long as this slope is not zero.

When the dissolution boundary reaches the boundary close to x  =  L the critical radius catches up with the mean droplet size. Concomitantly, narrowing of the droplet size distribution stops. Because all droplets are approximately of equal size, the exchange of material between droplets via Ostwald ripening is slowed down dramatically. This slowing down of inter-droplet diffusion via Ostwald ripening leads to a long phase where droplet number and size are almost constant. Close to the end of this arrest phase, the droplet distribution begins to broaden slowly, and the dynamics approaches classical Ostwald ripening.

In summary, a concentration gradient of a regulator component that affects phase separation can significantly change the dynamics of droplet coarsening. The regulator gradient causes an inhomogeneity of the equilibrium concentration and the concentration field far away from the droplet. Both induce a position-dependent ripening process where droplets can drift along the gradient and dissolve everywhere besides in a region close to one boundary of the system. During this positioning process of droplets to one boundary, the droplet size distribution can dramatically narrow for steep enough quenches, which causes a transient arrest of droplet growth. After this arrest phase, the positioned droplets are subject to a locally homogenous environment and the system recovers the dynamics of classical Ostwald ripening.

We start by considering an isolated system where the two-chemical species, the building block B and the precursor P, are converted into each other by the reaction

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \displaystyle P \rightleftharpoons B \label{eqn:reaction1} . \nonumber \end{align} \tag{ R1 }$$

At constant temperature T and volume $V$, the system is described by a free energy $F(N_P,N_B)$, where N_i are the particle numbers of type $i=P,B$. The thermodynamic equilibrium of the system corresponds to the minimum of F. The necessary condition for this minimum reads $\newcommand{\diff}{{{\rm d}}} \mu_P \diff N_P + \mu_B \diff N_B = 0$, where the chemical potentials are $\newcommand{\ri}{{\rm I}} \newcommand{\takenat}[2]{\left.{#1}\right|_{#2}} \mu_P = \takenat{\partial F/\partial N_P}{N_B}$ and $\newcommand{\ri}{{\rm I}} \newcommand{\takenat}[2]{\left.{#1}\right|_{#2}} \mu_B = \takenat{\partial F/\partial N_B}{N_P}$. Note that in contrast to phase separation without reactions, species can now be converted into each other and only the total number of particles, $M=N_P + N_B$, is conserved. This implies $\newcommand{\diff}{{{\rm d}}} \diff N_P = -\diff N_B$, such that the equilibrium condition requires $\mu_P - \mu_B=0$. Consequently, at equilibrium, the chemical reaction equalises the chemical potentials of the two species.

The difference between the chemical potentials, $\mu_P - \mu_B$, also affects the relaxation rate toward equilibrium. This rate is quantified by the total reaction flux $\newcommand{\diff}{{{\rm d}}} s=-\diff c_P/\diff t = \diff c_B/\diff t$, where $c_i = N_i/V$ denotes the concentrations in the homogeneous system for $i=P,B$. Since the reaction can proceed in both directions, the total reaction flux $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\sF}[1]{s^{(#1)}_{\rightarrow}} \newcommand{\ri}{{\rm I}} \newcommand{\sFrh}{s_{\rightarrow}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBrh}{s_{\leftarrow}} s= \sFrh - \sBrh$ is given by the difference of the forward reaction flux $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\sF}[1]{s^{(#1)}_{\rightarrow}} \newcommand{\ri}{{\rm I}} \newcommand{\sFrh}{s_{\rightarrow}} \sFrh$ associated with the conversion of P to B and the reverse flux $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBrh}{s_{\leftarrow}} \sBrh$. As a consequence of detailed balance, the ratio of the two reaction fluxes obeys (see appendix C)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ri}{{\rm I}} \newcommand{\sBrh}{s_{\leftarrow}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sFrh}{s_{\rightarrow}} \newcommand{\sF}[1]{s^{(#1)}_{\rightarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\kb}{k_{\rm B}} \displaystyle \frac{\sFrh}{\sBrh} = \exp\left(-\frac{\mu_B - \mu_P}{\kb T} \right) \label{eqn:foward_backward_ratio} , \nonumber \end{align} \tag{ 4.1 }$$

which we call detailed balance of the rates [147]. The relation shows that the net direction in which the reaction proceeds depends on the sign of the chemical potential difference $\mu_B - \mu_P$. Moreover, the net reaction flux s vanishes at chemical equilibrium ($\mu_P = \mu_B$) since $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\sF}[1]{s^{(#1)}_{\rightarrow}} \newcommand{\ri}{{\rm I}} \newcommand{\sFrh}{s_{\rightarrow}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBrh}{s_{\leftarrow}} \sFrh = \sBrh$. Close to chemical equilibrium, equation (4.1) can be linearized and the reaction flux $\newcommand{\s}[1]{s^{(#1)}} \newcommand{\sF}[1]{s^{(#1)}_{\rightarrow}} \newcommand{\ri}{{\rm I}} \newcommand{\sFrh}{s_{\rightarrow}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBrh}{s_{\leftarrow}} s= \sFrh - \sBrh$ can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\mobR}{\mobth_{\rm r}} \newcommand{\mobth}{\Lambda} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\be}{\beta} \displaystyle s = -\mobR(c_P,c_B) \, (\mu_B - \mu_P) \label{eqn:reaction_rate} , \nonumber \end{align} \tag{ 4.2 }$$

where the function $\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \mobR(c_P, c_B)$ determines the reaction rate. $\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \mobR(c_P, c_B)$ is an Onsager coefficient, which must be positive to ensure a positive entropy production rate [108, 123]; see appendix B.

To describe chemical reactions in inhomogeneous systems, we assume local thermodynamic equilibrium, i.e. there exist local volume elements where thermodynamic quantities such as concentrations and temperature can be defined. This is possible when the local volumes equilibrate quickly compared to the rates of the processes that we want to describe. In particular, the exchange with neighbouring volumes, associated with diffusive transport, and the conversion of particles into different species, associated with chemical reactions, should take place on timescales longer than the equilibration timescale of the volumes. If this is the case, a system of reacting and diffusing particles can be described by concentration fields $\newcommand{\vetr}{\boldsymbol} c_i(\vetr r)$ for all species i.

To study the interplay of the chemical reaction (R1) with phase separation, we first consider a binary, incompressible system described by the concentration of the droplet component, $\newcommand{\vetr}{\boldsymbol} c_B(\vetr r)=c(\vetr r)$, while $\newcommand{\vetr}{\boldsymbol} c_P(\vetr r) = \nu^{-1} - c(\vetr r)$ with $\nu$ denoting the molecular volume of P and B. The behavior of the system is governed by the free energy F[c], which is now a functional of the concentration field $\newcommand{\vetr}{\boldsymbol} c(\vetr r)$. For simplicity, we here consider the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ri}{{\rm I}} \newcommand{\be}{\beta} \newcommand{\diff}{{{\rm d}}} \newcommand{\abs}[1]{|#1|} \displaystyle F[c] = \int \diff^3 r \left(\,f(c) + \frac{\kappa}{2} \abs{\nabla c}^2 \right) , \nonumber \end{align} \tag{ 4.3 }$$

which combines a local contribution of the free energy density $f(c)$ with a term that accounts for the free energy costs of spatial inhomogeneities proportional to $\kappa$, analogous to section 2.1.5. The exchange chemical potential $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu = \mu_B - \mu_P$ is thus given by $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu = \delta F[c]/\delta c$. The resulting equilibrium condition of the chemical reaction is $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu(\vetr r)=0$, which includes the equilibrium condition for phase separation, $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu(\vetr r) = {{\rm const.}}$; see section 2.2.1.

The dynamical equation of the system follows from the conservation law

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \partial_t c + \nabla \cdot \vetr j = s \label{eqn:conservation_law} , \nonumber \end{align} \tag{ 4.4 }$$

where $\newcommand{\vetr}{\boldsymbol} \vetr j$ is the diffusive flux and s is the net flux of the production of species B by the reaction (R1). These two thermodynamic fluxes are driven by their respective conjugated forces $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \nabla \bar\mu$ and $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu$; see appendix B. Using linear response theory, $\newcommand{\vetr}{\boldsymbol} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\ba}{\vetr{a}} \vetr j = -\mobth \nabla \bar\mu$ and $\newcommand{\vetr}{\boldsymbol} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \newcommand{\ba}{\vetr{a}} s=-\mobR\bar\mu$, we arrive at the dynamical equation

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ba}{\vetr{a}} \newcommand{\mobR}{\mobth_{\rm r}} \newcommand{\mobth}{\Lambda} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \partial_t c = \nabla \cdot \bigl[ \mobth(c) \, \nabla \bar\mu(c) \bigr] - \mobR(c) \, \bar\mu(c) \label{eqn:active_pde_equilibrium} , \nonumber \end{align} \tag{ 4.5 }$$

which describes a binary system that exhibits phase separation and chemical reactions. Note that we recover the Cahn–Hilliard equation if chemical reactions are absent ($\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \mobR=0$); see equation (2.37). Conversely, in the limit where $\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \mobR$ is constant and diffusive fluxes vanish ($\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \mobth=0$), we obtain the Allen–Cahn model [148], which is the deterministic version of model A [124].

We study the effects of chemical reactions by first analyzing the homogeneous equilibrium states $\newcommand{\vetr}{\boldsymbol} c(\vetr r)=c_0$, which are governed by the equilibrium condition $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu(\vetr r)=0$. This condition implies vanishing chemical reaction flux, see equation (4.2), and $f'(c_0) =0$, so that c₀ is a (local) extremum of the free energy density $f(c)$. To assess the stability of these states, we consider harmonic perturbations with wave vector $\newcommand{\vetr}{\boldsymbol} \vetr q$, as described in section 2.4. In the linear regime, these perturbations grow with a rate

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\om}{\omega} \newcommand{\mobR}{\mobth_{\rm r}} \newcommand{\mobth}{\Lambda} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\be}{\beta} \newcommand{\vetr}{\boldsymbol} \displaystyle \omega(\vetr q) = -\bigl[\mobth(c_0) \vetr q^2 + \mobR(c_0)\bigr] \bigl[\,f''(c_0) + \kappa \vetr q^2\bigr]. \label{eqn:pert_growth_rate_passive} \nonumber \end{align} \tag{ 4.6 }$$

The system is stable if all perturbations decay, i.e. if $\newcommand{\vetr}{\boldsymbol} \newcommand{\om}{\omega} \omega(\vetr q) < 0$ for all wave vectors $\newcommand{\vetr}{\boldsymbol} \vetr q$. Since $\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \newcommand{\g}{\gamma} \mobth,\mobR\geqslant0$, the stability is governed by the sign of the second bracket in equation (4.6) and the homogeneous state becomes unstable when $f''(c_0) < 0$ [149]. This condition is identical to the condition for the spinodal instability in the case without chemical reactions ($\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobR}{\mobth_{\rm r}} \mobR=0$). However, in the presence of chemical reactions, only the homogeneous states with $f'(c_0)=0$ are stationary states because particles numbers of P and B are not conserved. In contrast, in the absence of chemical reactions, all homogeneous states are stationary with the homogeneous concentrations of P and B being conserved. Particle conservation also implies $\newcommand{\om}{\omega} \omega(0)=0$, while the q  =  0 mode is unstable when chemical reactions are present; see figure 8(b).

Zoom In Zoom Out Reset image size

Taken together, we showed that if the system settles in a homogeneous state it will attain minimal free energy [149, 150]; see figure 8(a). The major difference to the case without chemical reactions is that the species are not conserved individually, and the system can thus relax by altering the composition locally. In the next section, we will show that this local conversion destabilises all inhomogeneous states including the ones corresponding to coexisting phases.

We now investigate the stability of inhomogeneous states to see how chemical reactions that can relax to equilibrium ($\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu=0$) affect droplets. As an example, we first discuss the Ginzburg–Landau free energy presented in equation (2.19). When chemical reactions are present, this free energy permits two stable homogeneous solutions, corresponding to the two minima of the free energy density. Without chemical reactions, we have shown in section 2.2 that there is also a stable inhomogeneous stationary state, which consists of two bulk phases separated by an interfacial region. In the simple case of a one-dimensional system, the interfacial profile $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInt}{c_{\rm I}} \cInt(x)$ is given by equation (2.28). This interfacial profile is also a stationary state in the case with chemical reactions, since the symmetric free energy implies vanishing chemical potentials in the two bulk phases and thus satisfies the equilibrium condition $\newcommand{\vetr}{\boldsymbol} \bar\mu(\vetr r)=0$.

To scrutinise the stability of the interfacial profile in the presence of chemical reactions, we determine whether a small change $\delta c$ in the concentration profile could possibly decrease the free energy. Mathematically, this corresponds to calculating the second variation $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInt}{c_{\rm I}} \Delta F [\cInt,\delta c]$ of the free energy, which is given by equation (2.30). The stationary state is stable if $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInt}{c_{\rm I}} \Delta F [\cInt,\delta c]$ is positive for every nonzero variation $\delta c$. We show in appendix A that indeed almost all perturbations $\delta c$ increase the free energy. The only exception is $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInt}{c_{\rm I}} \delta c = \partial_x \cInt$, for which $\Delta F = 0$, implying that this perturbation does not decay in time and the state is marginally stable. This perturbation corresponds to an infinitesimal translation, indicating that interfaces can move without changing the total free energy when chemical reactions are present. Note that this perturbation does not conserve the mass of the individual species and is thus forbidden in the case without chemical reactions discussed in section 2.2.2. Taken together, we showed that the sigmoidal interface profile given by equation (2.28) is neither stable nor unstable in the special case of the Ginzburg–Landau free energy, where both minima have the same energy.

In the general case where the minima of the free energy density are at different energies (see figure 1(a)), the global free energy decreases when the interface moves such that the phase with the smaller free energy density expands [150]. Consequently, coexisting phases in asymmetric free energies are unstable due to the presence of chemical reactions. Similarly, curved interfaces are unstable (also in the case of a symmetric free energy density), since the associated Laplace pressure implies elevated concentrations on both sides of the interface; see equation (2.45). The Laplace pressure, and thus the concentrations and the free energy, decrease when the interface moves towards its concave side, implying that droplets shrink [150]. Taken together, these arguments show that inhomogeneous states are generally unstable when chemical reactions are present, and the system attains the global free energy minimum everywhere; see figure 8(a).

A hint of these dynamics is visible in the simulation results shown in the middle column in figure 10, where the chemical potential difference $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu$ is very close to zero in the bulk phases but finite at the interfaces. In fact, the chemical potential deviates from zero more strongly at interfacial regions of larger curvature. One consequence of this observations is that the local entropy production $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \Lambda_r\bar\mu^2$ by the chemical reactions is also largest at the interfaces.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Taken together, we showed that droplets are destabilised by the simple conversion reaction (R1) obeying detailed balance of the rates. In particular, the system always settles in a homogeneous state, even if droplets appear initially due to a spinodal instability or nucleation. The reason for the initial appearance of droplets is that diffusion is fast on small length scales and the formation of droplets locally reduces the free energy density. On longer time scales, the chemical reaction drives the system into a homogeneous equilibrium state. To compensate this destabilisation of droplets due to the reaction (R1) we will introduce an additional chemical reaction in the bulk phases. This chemical reaction is driven by fuel, and thereby breaks detailed balance of the rates and can prevent the chemical reaction to drive the system to a homogeneous equilibrium state. This additional chemical reaction allows us to control the behaviour of droplets. In particular, we will show that systems where droplet material is produced in one phase while it is destroyed in the other can lead to interesting phenomena such as mono-disperse emulsions or dividing droplets.

We next discuss the dynamics of active droplets in the case of strong phase separation, where the interfacial width w is small compared to the droplet radius R. In this case, the volume occupied by the interface and thus the chemical reactions inside the interfacial region are negligible. Conversely, we will show that the chemical reactions producing and destroying droplet material in the bulk phases influence the droplet dynamics significantly. We here consider a coarse-grained description, where a thin interface separates the inside of the droplet with a high concentration of droplet material B from the dilute phase outside, analogous to section 2.5.2. The dynamics in both phases is described by equation (4.11), but since the concentration variations are small within the phases, it can be approximated by a reaction-diffusion equation [74],

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sLoc}{s} \newcommand{\be}{\beta} \newcommand{\al}{\alpha} \displaystyle \partial_t c \simeq D_\alpha \nabla^2 c + \sLoc(c) \label{eqn:4_reaction_diffusion} , \nonumber \end{align} \tag{ 4.14 }$$

where $\newcommand{\cbase}{c^{(0)}} \newcommand{\al}{\alpha} \newcommand{\cBase}{\cbase} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\mobth}{\Lambda} \newcommand{\mobTot}{\mob{\mathrm{tot}}} \newcommand{\cB}{{{\mathcal B}}} D_\alpha = \mobth(\cBase_\alpha)\,f''(\cBase_\alpha) + \kappa\mobTot(\cBase_\alpha)$ denotes the diffusivity in the two phases $\newcommand{\al}{\alpha} \alpha={{\rm in}}, {{\rm out}}$. The first term in the expression for $\newcommand{\al}{\alpha} D_\alpha$ stems from the conservative fluxes and is thus equivalent to equation (2.41) while the second term captures the apparent diffusion due to chemical conversion. Note that the diffusivity $\newcommand{\DIn}{D_{\rm in}} \DIn$ inside the droplet is generally different than the diffusivity $\newcommand{\DOut}{D_{\rm out}} \DOut$ outside. The same approximation that led to equation (4.14) can also be used to linearize the reaction flux in the two phases,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\kOut}{k_{\rm out}} \newcommand{\kIn}{k_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\sLoc}{s} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \displaystyle \sLoc(c) \simeq \left\{\begin{array}{@{}ll@{}} \sBaseIn - \kIn (c - \cBaseIn) & {{\rm inside}} \nonumber \\ \sBaseOut - \kOut (c - \cBaseOut) & {{\rm outside}} \end{array}\right. \label{eqn:reaction_flux_linear} , \nonumber \end{align} \tag{ 4.15 }$$

where $\newcommand{\cbase}{c^{(0)}} \newcommand{\sLoc}{s} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} \sBaseIn = \sLoc(\cBaseIn)$ and $\newcommand{\cbase}{c^{(0)}} \newcommand{\sLoc}{s} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \sBaseOut = \sLoc(\cBaseOut)$ are the reaction fluxes when the concentrations are at their equilibrium values $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} \cBaseIn$ and $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \cBaseOut$ in the two phases, respectively. Deviations from these values are accounted for by $\newcommand{\cbase}{c^{(0)}} \newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\kIn}{k_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} \kIn = -\sLoc'(\cBaseIn)$ and $\newcommand{\cbase}{c^{(0)}} \newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\kOut}{k_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \kOut = -\sLoc'(\cBaseOut)$, which can be interpreted as elasticity coefficients of the chemical reactions [164].

The basal fluxes $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn$ and $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut$ need to have opposite sign for droplets to exist. If they had the same sign, droplet material would either be destroyed ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseIn, \sBaseOut < 0$) or produced everywhere ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseIn, \sBaseOut > 0$), which both implies unstable systems. To see under which conditions the basal fluxes have opposite sign, we express them as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ba}{\vetr{a}} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\vetr}{\boldsymbol} \displaystyle \sBaseIn = -\mob{1}(\cBaseIn) \, \bar\mu -\mob{2}(\cBaseIn)(\bar\mu - \bar\mu_2) , \nonumber \end{align} \tag{ 4.16a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ba}{\vetr{a}} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBase}{\cbase} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\vetr}{\boldsymbol} \displaystyle \sBaseOut = -\mob{1}(\cBaseOut) \, \bar\mu -\mob{2}(\cBaseOut)(\bar\mu - \bar\mu_2) , \nonumber \end{align} \tag{ 4.16b }$$

using equation (4.10). To give a concrete example, we here consider weak chemical reactions, so that the diffusive fluxes almost equilibrate the exchange chemical potential $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu$, which is therefore approximately homogeneous. If the building block B is of higher energy than the precursor P ($\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu > 0$), equation (4.16) imply that droplet material is produced outside the droplet ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut>0$) and destroyed within ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn < 0$) if the fuel supplies sufficient energy ($\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu_2 > \bar\mu > 0$) and the Onsager coefficients obey

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ba}{\vetr{a}} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\vetr}{\boldsymbol} \displaystyle \frac{\mob{1}(\cBaseOut)}{\mob{2}(\cBaseOut)} < \frac{\bar\mu_2}{\bar\mu} - 1 < \frac{\mob{1}(\cBaseIn)}{\mob{2}(\cBaseIn)} \label{eqn:reaction_onsager_condition} . \nonumber \end{align} \tag{ 4.17 }$$

For instance, if $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu$ is equal in the two phases and the external energy input is given by $\newcommand{\vetr}{\boldsymbol} \newcommand{\ba}{\vetr{a}} \bar\mu_2 = 2\bar\mu$, reaction (R1) must be faster than reaction (R2) inside the droplet, $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\cB}{{{\mathcal B}}} \mob{1}(\cBaseIn)>\mob{2}(\cBaseIn)$, while the opposite is true outside, $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \newcommand{\cB}{{{\mathcal B}}} \mob{1}(\cBaseOut) < \mob{2}(\cBaseOut)$. In this case, the conversion reaction (R1) proceeds from the building block B to the precursor P spontaneously ($\newcommand{\s}[1]{s^{(#1)}} \s{1}<0$), while reaction (R2) converts fuel to waste to produce the high-energy building block B from P ($\newcommand{\s}[1]{s^{(#1)}} \s{2}>0$); see figure 9. If the coefficients $\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \mob{1}$ and $\newcommand{\mob}[1]{\Lambda_{\rm r}^{(#1)}} \mob{2}$ obey equation (4.17) the total reaction flux $\newcommand{\sTot}{s_{\rm tot}} \newcommand{\s}[1]{s^{(#1)}} \sTot = \s{1} + \s{2}$ is then positive outside the droplet, while it is negative inside. A concrete implementation of such a system is discussed in the Supporting Information of [77].

We showed that droplet material can be produced in one phase while it is destroyed in the other when the basal fluxes $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn$ and $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut$ have opposite sign. Do similar conditions apply to the elasticity coefficients $\newcommand{\kIn}{k_{\rm in}} \kIn$ and $\newcommand{\kOut}{k_{\rm out}} \kOut$? We showed in section 4.2 that negative $\newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \sLoc'(c)$ has a stabilising effect and we would thus expect that droplets can persist if k  >  0; see also figure 8(b). Conversely, since positive $\newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \sLoc'(c)$ can destabilise homogeneous phases, the active droplets we discuss here might be unstable if k  <  0. However, there is a smallest length scale $q_{\rm max}^{-1}$ below which the instability cannot develop. This length scale can be determined from equation (4.12) and the condition $\newcommand{\om}{\omega} \omega(q_{\rm max})=0$, yielding

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \displaystyle q_{\rm max}^{-1} \simeq \left[\frac{\zeta(\cBaseIn)}{s'(\cBaseIn)} \right]^{\frac12} , \nonumber \end{align} \tag{ 4.18 }$$

for weak reactions, $\newcommand{\cbase}{c^{(0)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} s'(\cBaseIn) \ll \zeta(\cBaseIn)$. In the case where $q_{\rm max}^{-1}$ is large compared to the droplet radius, the instability is effectively suppressed, and the droplet could be stable even if $\newcommand{\kIn}{k_{\rm in}} \kIn < 0$. Conversely, the dilute phase will typically be large compared to $q_{\rm max}^{-1}$, such that an instability would develop there if $\newcommand{\kOut}{k_{\rm out}} \kOut < 0$. In the following, we thus consider all values of $\newcommand{\kIn}{k_{\rm in}} \kIn$, but restrict our discussion to $\newcommand{\kOut}{k_{\rm out}} \kOut > 0$.

We distinguish different classes of active droplets based on where droplet material is produced. If the reaction fluxes $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn$ and $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut$ have the same sign, droplet material is produced or destroyed in the entire system and such states are unstable. Consequently, stable droplets can only exist if $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn$ and $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut$ have opposite sign: we denote the case where droplet material is produced outside the droplets ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut > 0$, $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn < 0$) as externally maintained droplets, while the opposite case where droplets produce their own material ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn > 0$, $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut < 0$) is called internally maintained droplets.

We begin by studying a single, isolated active droplet surrounded by a large dilute phase. The droplet grows if there is a net flux of droplet material toward its surface, see equation (2.47). In the simple case of a spherical droplet of given radius R and volume $\newcommand{\Vd}{V_{\rm d}} \Vd= \frac{4\pi}{3} R^3$, the growth rate reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\FluxOut}{J_{\rm out}} \newcommand{\FluxIn}{J_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\Fl}{F^{-}} \newcommand{\Vd}{V_{\rm d}} \newcommand{\cbase}{c^{(0)}} \newcommand{\difffrac}[2]{\frac{\diff #1}{\diff #2}} \newcommand{\diff}{{{\rm d}}} \displaystyle \difffrac{\Vd}{t} \simeq \frac{\FluxIn - \FluxOut}{\cBaseIn - \cBaseOut} \label{eqn:droplet_growth_sym} , \nonumber \end{align} \tag{ 4.19 }$$

where we have neglected surface tension effects in the denominator and defined the integrated surface fluxes $\newcommand{\vetr}{\boldsymbol} J_{\rm in/out}=4\pi R^2 \, \vetr{n} \cdot \vetr{j}_{\rm in/out}$. These fluxes can be determined from the stationary solutions $\newcommand{\steady}[1]{#1_*} \newcommand{\s}[1]{s^{(#1)}} \steady{c}(r)$ that follow from solving equation (4.14) with the boundary conditions at the interface given in equation (2.45), see figure 11(a).

Zoom In Zoom Out Reset image size

The flux $\newcommand{\Fl}{F^{-}} \newcommand{\FluxIn}{J_{\rm in}} \FluxIn$ inside the droplet interface can be obtained in the quasi-static limit, where it equals the reaction flux $\newcommand{\vetr}{\boldsymbol} \newcommand{\diff}{{{\rm d}}} \newcommand{\steady}[1]{#1_*} \newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\ri}{{\rm I}} \newcommand{\sF}[1]{s^{(#1)}_{\rightarrow}} \newcommand{\sFluxIn}{S_{\rm in}} \sFluxIn = \int \diff^3 r \, \sLoc(\steady{c}(\vetr r))$ inside the droplet volume. In the typical case where the radius R is small compared to the length $\newcommand{\abs}[1]{|#1|} \newcommand{\DIn}{D_{\rm in}} \newcommand{\kIn}{k_{\rm in}} \newcommand{\e}{{\rm e}} \newcommand{\DiffLenIn}{\ell_{\rm in}} \DiffLenIn= (\DIn/\abs{\kIn}){}^{\frac12}$ generated by the reaction-diffusion system, the concentration inside the droplet is $\newcommand{\vetr}{\boldsymbol} \newcommand{\cbase}{c^{(0)}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cB}{{{\mathcal B}}} c_*(\vetr r) \simeq \cBaseIn$, implying

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\FluxIn}{J_{\rm in}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\Fl}{F^{-}} \newcommand{\Vd}{V_{\rm d}} \displaystyle \FluxIn \simeq \sBaseIn \Vd . \label{eqn:active_flux_inside} \nonumber \end{align} \tag{ 4.20 }$$

Droplet material is thus transported toward the interface ($\newcommand{\Fl}{F^{-}} \newcommand{\FluxIn}{J_{\rm in}} \FluxIn>0$) only in internally maintained droplets ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn > 0$).

The integrated flux $\newcommand{\Fl}{F^{-}} \newcommand{\FluxOut}{J_{\rm out}} \FluxOut$ outside the droplet interface can also be obtained in a quasi-static approximation. In contrast to the passive case discussed in section 2.5.2, the supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\cbase}{c^{(0)}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cIn}{c_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cInfty}{c_\infty} \newcommand{\cB}{{{\mathcal B}}} \supSat = (\cInfty - \cBaseOut)/\cBaseOut$ far away from droplets is now created by chemical reactions. In the typical case where the dilute phase is large compared to the length scale $\newcommand{\abs}[1]{|#1|} \newcommand{\DOut}{D_{\rm out}} \newcommand{\kOut}{k_{\rm out}} \newcommand{\e}{{\rm e}} \newcommand{\DiffLenOut}{\ell_{\rm out}} \DiffLenOut = (\DOut/\abs{\kOut}){}^{\frac12}$, the chemical reactions equilibrate far away from droplets and the composition thus reaches the value $\newcommand{\cbase}{c^{(0)}} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\cBase}{\cbase} \newcommand{\cIn}{c_{\rm in}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cInfty}{c_\infty} \newcommand{\kOut}{k_{\rm out}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \cInfty = \cBaseOut + \sBaseOut/\kOut$, so $\newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cIn}{c_{\rm in}} \newcommand{\cInfty}{c_\infty} \sLoc(\cInfty) = 0$; see figure 11(b). The associated supersaturation reads $\newcommand{\supSat}{\varepsilon} \newcommand{\cbase}{c^{(0)}} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\kOut}{k_{\rm out}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \supSat = \sBaseOut/(\kOut\cBaseOut)$ and is thus positive if $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut > 0$ since the reaction must be auto-inhibitory outside droplets ($\newcommand{\kOut}{k_{\rm out}} \kOut > 0$); see section 4.2.1. The resulting transport of droplet material can be quantified by the flux outside the droplet interface, which reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\FluxOut}{J_{\rm out}} \newcommand{\DOut}{D_{\rm out}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBase}{\cbase} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\Fl}{F^{-}} \newcommand{\cbase}{c^{(0)}} \newcommand{\supSat}{\varepsilon} \displaystyle \FluxOut \simeq 4 \pi \DOut R \cBaseOut \left(\frac{\ell_\gamma}{R} - \supSat\right) \label{eqn:active_flux_outside} \nonumber \end{align} \tag{ 4.21 }$$

in the typical case $\newcommand{\e}{{\rm e}} \newcommand{\DiffLenOut}{\ell_{\rm out}} R \ll \DiffLenOut$. The first term in the bracket captures the effect of surface tension, which is typically small. Neglecting this term, we find that droplet material is transported toward the interface ($\newcommand{\Fl}{F^{-}} \newcommand{\FluxOut}{J_{\rm out}} \FluxOut < 0$) if the dilute phase is supersaturated ($\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat > 0$), which is the case only in externally maintained droplets ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \sBaseOut>0$).

The droplet growth rate following from combining equation (4.19)–(4.21) reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\DOut}{D_{\rm out}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBaseIn}{\cbase_{\rm in}} \newcommand{\cBase}{\cbase} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\supSat}{\varepsilon} \newcommand{\difffrac}[2]{\frac{\diff #1}{\diff #2}} \newcommand{\diff}{{{\rm d}}} \displaystyle \difffrac{R}{t} \simeq \frac{\DOut\cBaseOut}{R(\cBaseIn - \cBaseOut)}\left(\supSat - \frac{\ell_\gamma}{R} + \frac{R^2\sBaseIn}{3\DOut \cBaseOut} \right) \label{eqn:active_droplet_growth} . \nonumber \end{align} \tag{ 4.22 }$$

The first term in the bracket describes droplet growth due to a supersaturated environment ($\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat > 0$) or shrinkage in undersaturated environments ($\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat < 0$). The second term, which is only relevant for small droplets, captures the reduction of growth due to surface tension $\newcommand{\g}{\gamma} \gamma$. The last term describes the growth of the droplet due to production of droplet material inside if $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn>0$ or its shrinking when $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn < 0$. Note that equation (4.22) reduces to equation (2.53) if chemical reactions are absent ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn=0$) and the supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat$ is imposed.

We begin by discussing the growth of a single droplet in an environment with constant supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat$. This corresponds for instance to an infinite system where the supersaturation reaches its equilibrium value $\newcommand{\supSat}{\varepsilon} \newcommand{\supSatEq}{\supSat_{\rm eq}} \newcommand{\cbase}{c^{(0)}} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\kOut}{k_{\rm out}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \supSatEq=\sBaseOut/(\kOut\cBaseOut)$ far away from an isolated droplet. The stationary states of this system can be determined from equation (4.22) and correspond to radii R for which the bracket vanishes. The associated cubic equation in R has at most three solutions, which we now classify. If the chemical reactions are too strong, there are no physical solutions. In particular, if $\newcommand{\abs}[1]{|#1|} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseInMax}{\sBaseIn^{\rm max}} \abs\sBaseIn > \sBaseInMax$ with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\sBaseInMax}{\sBaseIn^{\rm max}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\DOut}{D_{\rm out}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBase}{\cbase} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\cbase}{c^{(0)}} \newcommand{\supSat}{\varepsilon} \newcommand{\abs}[1]{|#1|} \displaystyle \sBaseInMax = \frac{4 \cBaseOut \DOut \abs{\supSat}^3}{9 \ell_\gamma^2} , \nonumber \end{align} \tag{ 4.23 }$$

two solutions are complex while the third one is either negative (if $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn<0$) or smaller than the interface width (if $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn>0$), which is both unphysical. Consequently, stable droplets can only exist for moderate chemical reactions, $\newcommand{\abs}[1]{|#1|} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseInMax}{\sBaseIn^{\rm max}} \abs\sBaseIn < \sBaseInMax$. In this case, the polynomial equation has three real solutions. One of the solutions is always negative and thus unphysical, while the other two read

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\g}{\gamma} \newcommand{\RsA}{R^{(1)}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\be}{\beta} \newcommand{\R}{c_{R}} \newcommand{\supSat}{\varepsilon} \displaystyle \RsA \simeq \frac{\ell_\gamma}{\supSat} \, , \nonumber \end{align} \tag{ 4.24a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\g}{\gamma} \newcommand{\cB}{{{\mathcal B}}} \newcommand{\ri}{{\rm I}} \newcommand{\RsB}{R^{(2)}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\DOut}{D_{\rm out}} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cBase}{\cbase} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\R}{c_{R}} \newcommand{\cbase}{c^{(0)}} \newcommand{\supSat}{\varepsilon} \displaystyle \RsB \simeq \left(\frac{3\DOut\supSat \cBaseOut}{-\sBaseIn}\right)^{\frac12} -\frac{\ell_\gamma}{2\supSat} \label{eqn:stationary_state_radiiB} , \nonumber \end{align} \tag{ 4.24b }$$

which is correct up to linear order in $\newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} \ell_\gamma/R$. The dynamics in the vicinity of these states follow from a linear stability analysis of equation (4.22), where the droplet radius is perturbed away from the stationary state and we determine whether the dynamics will bring the droplet back to its stationary state or not [74]. This gives us enough information to discuss the growth behavior of active droplets qualitatively.

In the case of externally maintained droplets ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn< 0$, $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat>0$), we have $\newcommand{\R}{c_{R}} \newcommand{\RsA}{R^{(1)}} \newcommand{\RsB}{R^{(2)}} 0 < \RsA < \RsB$, since $\newcommand{\abs}[1]{|#1|} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseInMax}{\sBaseIn^{\rm max}} \abs\sBaseIn < \sBaseInMax$. Here, the radius $\newcommand{\R}{c_{R}} \newcommand{\RsA}{R^{(1)}} \RsA$ corresponds to an unstable state, while $\newcommand{\R}{c_{R}} \newcommand{\RsB}{R^{(2)}} \RsB$ is stable, see figure 12. Droplets smaller than $\newcommand{\R}{c_{R}} \newcommand{\RsA}{R^{(1)}} \RsA$ dissolve and disappear, so $\newcommand{\R}{c_{R}} \newcommand{\RsA}{R^{(1)}} \RsA$ is a critical radius similar to the one discussed in section 2.5.2. Externally maintained droplets that are larger than the critical radius $\newcommand{\R}{c_{R}} \newcommand{\RsA}{R^{(1)}} \RsA$ grow until they reach the stable stationary state with radius $\newcommand{\R}{c_{R}} \newcommand{\RsB}{R^{(2)}} \RsB$. This state is not present for passive droplets in an infinite system and must thus be a consequence of the chemical reactions. This can be seen by analyzing the two fluxes $\newcommand{\Fl}{F^{-}} \newcommand{\FluxIn}{J_{\rm in}} \FluxIn$ and $\newcommand{\Fl}{F^{-}} \newcommand{\FluxOut}{J_{\rm out}} \FluxOut$, which must be equal in the stationary state. However, $\newcommand{\Fl}{F^{-}} \newcommand{\FluxIn}{J_{\rm in}} \FluxIn$ scales with the droplet volume, while $\newcommand{\Fl}{F^{-}} \newcommand{\FluxOut}{J_{\rm out}} \FluxOut$ scales with its radius, see equation (4.21).Consequently, if the droplet radius exceeds $\newcommand{\R}{c_{R}} \newcommand{\RsB}{R^{(2)}} \RsB$, the loss due to $\newcommand{\Fl}{F^{-}} \newcommand{\FluxIn}{J_{\rm in}} \FluxIn$ dominates and the droplet shrinks back to the stationary state. The chemical turnover inside the droplet thus stabilizes the stationary state. Note that the two stationary states given in equation (4.24) only exist if the chemical reactions are not too strong ($\newcommand{\abs}[1]{|#1|} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseInMax}{\sBaseIn^{\rm max}} \abs\sBaseIn < \sBaseInMax$). In the limiting case $\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sBaseInMax}{\sBaseIn^{\rm max}} \sBaseIn = -\sBaseInMax$ the situation is degenerated, and the two stationary states are identical. The corresponding radius $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\g}{\gamma} \newcommand{\e}{{\rm e}} R_{\rm min}^{\rm ext}= (3 \ell_\gamma)/(2\supSat)$ can be determined from equation (4.22) and corresponds to the smallest externally maintained droplet that can be stable.

Zoom In Zoom Out Reset image size

In the case of internally maintained droplets ($\newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \sBaseIn > 0$, $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat < 0$), the first solution given in equation (4.24) is negative and thus unphysical. The second solution is always positive, but unstable to perturbations. Consequently, $\newcommand{\R}{c_{R}} \newcommand{\RsB}{R^{(2)}} \RsB$ is the critical droplets size of internally maintained droplets, which can be significantly larger than critical sizes in externally maintained droplets. Nucleation is thus typically suppressed efficiently, but it can be promoted by catalytically active particles, which catalyze the production of droplet material at their surface [73, 74]. Internally maintained droplets larger than the critical size grow up to the system size and there is no characteristic stable size. This is because such droplets grow quicker if they become larger and this autocatalytic growth only stops when the dilute phase is depleted of material P, which can only happen in a finite system [74].

So far, we only considered systems that are large compared to the droplet size, so the supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat$ is effectively constant. In contrast, in small systems a growing droplet can deplete the surrounding dilute phase significantly. In particular, the average concentration $\newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \cOut$ of droplet material in the dilute phase evolves as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\FluxOut}{J_{\rm out}} \newcommand{\cOut}{c} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sLoc}{s} \newcommand{\be}{\beta} \newcommand{\Fl}{F^{-}} \newcommand{\Vd}{V_{\rm d}} \newcommand{\Vsys}{V} \newcommand{\difffrac}[2]{\frac{\diff #1}{\diff #2}} \newcommand{\diff}{{{\rm d}}} \displaystyle \difffrac{\cOut}{t} = \sLoc(\cOut) + \frac{\FluxOut}{\Vsys - \Vd} \label{eqn:supersaturation_dynamics_single} , \nonumber \end{align} \tag{ 4.25 }$$

where the first term on the right-hand side originates from the chemical reactions in the dilute phase and the last term accounts for the diffusive flux at the interface of the droplet of volume $\newcommand{\Vd}{V_{\rm d}} \Vd$. In large systems ($\newcommand{\Vsys}{V} \newcommand{\Vd}{V_{\rm d}} \newcommand{\g}{\gamma} \Vsys\gg \Vd$), the last term is negligible and $\newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \cOut$ relaxes to $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInfty}{c_\infty} \cInfty$, where chemical equilibrium is obeyed, $\newcommand{\sLoc}{s} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cIn}{c_{\rm in}} \newcommand{\cInfty}{c_\infty} \sLoc(\cInfty)=0$. At this point, the supersaturation $\newcommand{\supSat}{\varepsilon} \newcommand{\cbase}{c^{(0)}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\cBase}{\cbase} \newcommand{\cOut}{c} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\cO}{{{\mathcal O}}} \newcommand{\cB}{{{\mathcal B}}} \supSat= (\cOut - \cBaseOut)/\cBaseOut$ attains its equilibrium value $\newcommand{\supSat}{\varepsilon} \newcommand{\supSatEq}{\supSat_{\rm eq}} \newcommand{\cbase}{c^{(0)}} \newcommand{\sBase}{\Gamma} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\cBase}{\cbase} \newcommand{\cBaseOut}{\cbase_{\rm out}} \newcommand{\kOut}{k_{\rm out}} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\cB}{{{\mathcal B}}} \supSatEq=\sBaseOut/(\kOut\cBaseOut)$ and is thus independent of the droplet size, consistent with our assumptions in the previous subsection. In small systems, however, the last term in equation (4.25) is not negligible and $\newcommand{\supSat}{\varepsilon} \newcommand{\s}[1]{s^{(#1)}} \supSat$ depends on the droplet volume $\newcommand{\Vd}{V_{\rm d}} \Vd$. For instance, in the case of externally maintained droplets, we have $\newcommand{\Fl}{F^{-}} \newcommand{\FluxOut}{J_{\rm out}} \FluxOut < 0$ and the concentration $\newcommand{\cOut}{c} \newcommand{\cO}{{{\mathcal O}}} \cOut$ outside the droplet is thus lower than $\newcommand{\cIn}{c_{\rm in}} \newcommand{\cInfty}{c_\infty} \cInfty$. Consequently, the supersaturation is smaller for smaller systems and we would expect a reduced stationary droplet radius $\newcommand{\R}{c_{R}} \newcommand{\RsB}{R^{(2)}} \RsB$, see equation (4.24b). Indeed, when we solve for the stationary states of equation (4.22) and (4.25), we find in the simple case of large diffusion and small surface tension that there is a stationary state with volume

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\sBaseOut}{\sBase_{\rm out}} \newcommand{\sBaseIn}{\sBase_{\rm in}} \newcommand{\sB}[1]{s^{(#1)}_{\leftarrow}} \newcommand{\s}[1]{s^{(#1)}} \newcommand{\sBase}{\Gamma} \newcommand{\be}{\beta} \newcommand{\steady}[1]{#1_*} \newcommand{\Vsys}{V} \displaystyle \steady V \simeq \frac{\sBaseOut\Vsys}{\sBaseOut - \sBaseIn} \label{eqn:steady_volume_small} , \nonumber \end{align} \tag{ 4.26 }$$

which is stable for externally maintained droplets and unstable for internally maintained ones. Consequently, the size of stable stationary droplets that are externally maintained depends on system size in small systems, while it is independent of the system size in large systems, see equation (4.24).

Similar to passive droplets, we find that active droplets also have a critical radius, below which they shrink and disappear. Consequently, droplets can only be nucleated spontaneously if a concentration fluctuation creates a large enough initial droplet. In externally maintained droplets, this mechanism is similar to passive droplets, where the nucleation barrier is higher for larger surface tension and smaller supersaturation. Conversely, the critical radius is generally larger for internally maintained droplets, but it becomes smaller for larger surface tension and smaller supersaturation. This is because only a large enough droplet can produce enough droplet material to balance the efflux into the dilute phase. This efflux is smaller for larger surface tension, so surface tension actually helps to nucleate internally maintained droplets.

Once droplets exceed their critical radius they grow spontaneously by absorbing droplet material from the surrounding (passive and externally maintained droplets) or by producing more droplet material (internally maintained droplets). In the first case, the droplet growth rate is larger for smaller droplets and higher supersaturation. However, the growth of externally maintained droplets comes to a stop at a finite size, at which the loss of droplet material due to the chemical reactions inside is balanced by its influx over the surface. Conversely, internally maintained droplets grow indefinitely in an infinite system, similar to passive droplets. However, in contrast to passive droplets, this growth accelerates because of its autocatalytic nature. Consequently, in the simple case of a binary fluid in an infinite system, only externally maintained droplets reach a finite droplet size, while both passive and internally maintained droplets grow indefinitely.