Role of hydrodynamic flows in chemically driven droplet division

Living cells are compartmentalized in order to organize their biochemistry in space. Many cellular compartments do not possess membranes and are formed by the assembly of proteins and RNA in compact condensates [1–16]. Such condensates often have liquid like properties and resemble droplets that form by phase separation of a complex mixture [1, 11–13]. Indeed protein droplets are observed to form in vitro by phase separation in physiological buffer [13, 15, 17–19]. Such droplets can organize chemical reactions in space, and the droplet dynamics can in turn be influenced by the reactions, as has been shown both in theory [8, 20–26] and experiments [13, 15, 17–19, 27, 28]. The ubiquitous nature of RNA-protein condensates in a large variety of different cells and organisms suggests that the physical chemistry of macromolecular phase separation represents an evolutionary old mechanism for the compartmentalization of chemistry and that droplet formation could have played a key role at the origins of life and the emergence of prebiotic protocells [15, 18, 29–40].

A minimal model of a protocell consists of a droplet that turns over by a chemical reaction and is constantly supplied with droplet material by diffusion from the outside [39]. In such a scenario droplets are maintained away from thermodynamic equilibrium and can reach a non-equilibrium steady state with a radius that is set by reaction parameters [26]. An interesting possibility is that the spherical shape of active droplets becomes unstable and droplets spontaneously divide in two smaller daughters drops, providing a physical mechanism for the division of prebiotic cells [39]. Such droplet dynamics is a hydrodynamic problem because surface tension in non-spherical droplets drives hydrodynamic flows that redistribute material and deform the droplet shape [41–44].

Here we develop a hydrodynamic theory of the dynamics of chemically active droplets. We show that chemical reactions in active droplets can perform work against surface tension and flows, giving rise to a shape instability that can result in droplet division. We investigate the conditions for which droplets divide and determine hydrodynamic flow fields of dividing droplets.

We consider an incompressible liquid consisting of droplet material B and solvent component A which can phase separate. The local composition is characterized by the concentration field c(x) of component B. Volume preserving chemical reactions can transform component A into component B and back, $A\rightleftharpoons B$. For simplicity, we first discuss an effective droplet model. A single droplet characterized by high concentration c of component B coexists with the surrounding fluid that mainly consists of A and contains B at low concentration, see figure 1. Both phases are separated by a sharp interface. The concentration of B satisfies a balance equation, where the chemical reaction provide a source or sink term s_±(c),

$$\begin{eqnarray}&&{\partial }_{t}c+{\rm{\nabla }}\cdot {\boldsymbol{j}}={s}_{\pm }(c)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\boldsymbol{j}}=-{D}_{\pm }{\rm{\nabla }}c+{\boldsymbol{v}}c.\end{eqnarray} \tag{ 2 }$$

Here, the indices + and − refer to quantities outside and inside the droplet, respectively. The flux ${\boldsymbol{j}}$ consists of advection by the fluid velocity ${\boldsymbol{v}}$ and a diffusion flux, where D_± denotes the diffusion constant of the droplet material in the two phases.

Zoom In Zoom Out Reset image size

The chemical reaction is described by the concentration-dependent rate ${s}_{\pm }(c)$ which in general is a nonlinear function of c. For simplicity, we linearize the chemical reaction rates in the vicinity of reference concentrations ${c}_{\pm }^{(0)}$ in each phase:

$$\begin{eqnarray}&&{s}_{\pm }(c)\simeq -{k}_{\pm }(c-{c}_{\pm }^{(0)})\pm {\nu }_{\pm },\end{eqnarray} \tag{ 3 }$$

where ${c}_{\pm }^{(0)}$ are the equilibrium concentrations that coexist at equilibrium across a flat interface. We have defined the reaction rate ${\nu }_{\pm }=s({c}_{\pm }^{(0)})$ and the reaction constants ${k}_{\pm }={\rm{d}}s({c}_{\pm }^{(0)})/{\rm{d}}c$. We focus on the case of positive coefficients ${k}_{\pm }\gt 0$ and ${\nu }_{\pm }\gt 0$, where B is produced outside the droplet, and degraded inside, see figure 1.

The hydrodynamic flow velocity ${\boldsymbol{v}}$ obeys Stokes equation of an incompressible fluid,

$$\begin{eqnarray}&&{\eta }_{\pm }{{\rm{\nabla }}}^{2}{\boldsymbol{v}}={\rm{\nabla }}p,\end{eqnarray} \tag{ 4 }$$

which accounts for momentum conservation ∂_ασ_αβ = 0, where the stress tensor is given by σ_αβ = η_± (∂_αv_β + ∂_βv_α) − pδ_αβ. Here η_± denotes the fluid shear viscosities inside and outside of the droplet. The pressure p plays the role of a Lagrange multiplier to impose the constraint ${\rm{\nabla }}\cdot {\boldsymbol{v}}=0$.

The bulk equations (1)–(4) are connected by boundary conditions at the droplet interface which we parameterize in spherical coordinates by the radial interface position R(θ, ϕ) as a function of the polar and azimuthal angles θ and ϕ. The stress boundary conditions read

$$\begin{eqnarray}&&{\sigma }_{{nn}}^{+}(R)-{\sigma }_{{nn}}^{-}(R)=2\gamma H(R)\,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\sigma }_{{nt}}^{+}(R)-{\sigma }_{{nt}}^{-}(R)=0,\end{eqnarray} \tag{ 6 }$$

where H(R) is the local mean curvature of the interface and γ is the droplet surface tension. The stresses at the interface on the inner and outer side of the droplet are denoted by ${\sigma }_{\alpha \beta }^{\pm }(R)$. The tensor indices n and t refer to tensor components normal and tangential to the interface, respectively. The normal and tangential tensor components are defined as ${\sigma }_{{nn}}^{\pm }={n}_{\alpha }{\sigma }_{\alpha \beta }^{\pm }{n}_{\beta }$ and ${\sigma }_{{nt}}^{\pm }={n}_{\alpha }{\sigma }_{\alpha \beta }^{\pm }{t}_{\beta }$, where n_α is a unit vector normal to the surface and t_α is a unit vector tangent to the surface. Equation (6) is valid for all tangent vectors and summation over repeated indices is implied. Using no-slip boundary conditions, the velocity field is continuous at the interface,

$$\begin{eqnarray}&&{{\boldsymbol{v}}}^{+}(R)={{\boldsymbol{v}}}^{-}(R).\end{eqnarray} \tag{ 7 }$$

The concentration field c is discontinuous across the interface with values given by

$$\begin{eqnarray}&&{c}_{-}(R)={c}_{-}^{(0)}+{\beta }_{-}\gamma H(R)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{c}_{+}(R)={c}_{+}^{(0)}+{\beta }_{+}\gamma H(R)\end{eqnarray} \tag{ 9 }$$

which are set by the physics of phase coexistence and a local equilibrium assumption. The coefficients β_± describe the effects of the Laplace pressure on the equilibrium concentrations at phase coexistence. In the presence of fluxes at the interface, the interface moves in normal direction. The radial growth velocity is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{R}}{{\rm{d}}{t}}=\displaystyle \frac{{\boldsymbol{n}}}{{\boldsymbol{n}}\cdot {{\boldsymbol{e}}}_{r}}\cdot \displaystyle \frac{{{\boldsymbol{j}}}^{-}(R)-{{\boldsymbol{j}}}^{+}(R)}{{c}_{-}(R)-{c}_{+}(R)},\end{eqnarray} \tag{ 10 }$$

where ${\boldsymbol{n}}$ is a unit vector normal to the surface and ${{\boldsymbol{e}}}_{r}$ is a unit vector in radial direction. Equation (10) captures both convection of the interface by flows and droplet growth and shrinkage by addition or removal of material.

We find non-equilibrium steady state solutions to equations (1)–(10) with a spherical droplet of stationary radius $\bar{R}$ and stationary concentration field $\bar{c}(r)$, where r is the radial coordinate, see appendix A. The stationary pressure $\bar{p}$ exhibits a jump $2\gamma /\bar{R}$ across the interface and no hydrodynamic flows exist, $\bar{{\boldsymbol{v}}}=0$. An example for a stable non-equilibrium steady state with steady state concentration profile inside and outside the droplet of radius $\bar{R}$ is shown in figure 1.

We first discuss the properties of these stationary states as a function of external supersaturation $\epsilon ={\nu }_{+}/({k}_{+}{\rm{\Delta }}c)$ and the dimensionless reaction rate A = ν₋τ/Δc inside the droplet. The supersaturation is in our system generated by reactions outside the droplet and in steady state corresponds to the concentration for which s₊ = 0. Here, ${\rm{\Delta }}c={c}_{-}^{(0)}-{c}_{+}^{(0)}$ and we have introduced the time scale τ = w²/D₊, where w = 6β₊γ/Δc is a characteristic length scale. The stationary radii as a function of supersaturation are shown In figures 2(A)–(C) as solid lines for different values of A. For values of smaller than a threshold value ₀, no stationary radius exists. For values  > ₀ two steady state radii ${\bar{R}}_{c}$ and ${\bar{R}}_{s}$ exist, which become equal at ₀ where they approach a value ${\bar{R}}_{0}$. The smaller steady state radius ${\bar{R}}_{c}$ is a critical nucleation radius similar to the critical droplet radii found in passive systems. The larger radius denoted ${\bar{R}}_{s}$ stems from the interplay of phase separation and chemical reactions [26, 39]. As the supersaturation reaches a value ${\epsilon }_{\infty }=\sqrt{({D}_{-}{k}_{-})/({D}_{+}{k}_{+})}{\nu }_{-}/({k}_{-}{\rm{\Delta }}c)$, the stationary radius ${\bar{R}}_{s}$ diverges.

Zoom In Zoom Out Reset image size

We can find simple expressions for the stationary radii in the limit of small A while keeping the ratios ν₋/(k₋Δc) and k₊/k₋ of reaction parameters fixed. In this limit, the chemical reactions fluxes vanish as s_± ∝ A and the threshold value ₀ vanishes as ₀ ∝ A^1/3. The critical nucleation radius then behaves as ${\bar{R}}_{c}\simeq w/(6\epsilon )$ and the larger steady state radius ${\bar{R}}_{s}\simeq w{(3\epsilon A)}^{1/2}$ where ${\epsilon }_{0}\ll \epsilon \ll {\epsilon }_{\infty }$, see figure 2(B) and appendix A.5.

The steady state solutions are independent on the fluid viscosity, however the droplet dynamics is affected by hydrodynamic effects. We now investigate the role of hydrodynamic flows on chemically driven shape instabilities that can give rise to droplet division. We perform a linear stability analysis at the stationary state given by $\bar{X}=(\bar{c},\bar{R},\bar{p},\bar{{\boldsymbol{v}}})$ for small perturbations $\delta X=(\delta c,\delta R,\delta p,\delta {\boldsymbol{v}})$. The dynamics of these perturbations can be represented using eigenmodes

$$\begin{eqnarray}&&\delta X=\displaystyle \sum _{n,l,m}{\epsilon }_{{nlm}}{X}_{{nlm}}{{\rm{e}}}^{{\mu }_{{nlm}}t},\end{eqnarray} \tag{ 11 }$$

with ${X}_{{nlm}}=({c}_{{nl}}{Y}_{{lm}},\bar{R}{Y}_{{lm}},{p}_{l}{Y}_{{lm}},{{\boldsymbol{v}}}_{{lm}})$, where ${Y}_{{lm}}(\theta ,\phi )$ are spherical harmonics with angular mode indices with l = 0, 1, ...and m = −l, ..., l. The index n = 0, 1, ... denotes radial modes. The eigenmodes exhibit an exponential time dependence with a relaxation rate given by the eigenvalue μ_nlm. The mode amplitudes are denoted _nlm. The concentration modes are characterized by the radial functions c_nl(r). The pressure modes are described by p_l(r) and the velocity modes ${{\boldsymbol{v}}}_{{lm}}(r,\theta ,\varphi )$ can be expressed as

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{lm}}={v}_{{lm}}^{r}{{\boldsymbol{Y}}}_{{lm}}+{v}_{{lm}}^{(1)}{{\boldsymbol{\Psi }}}_{{lm}}+{v}_{{lm}}^{(2)}{{\boldsymbol{\Phi }}}_{{lm}},\end{eqnarray} \tag{ 12 }$$

where ${{\boldsymbol{Y}}}_{{lm}}(\theta ,\varphi )={{\boldsymbol{e}}}_{r}{Y}_{{lm}}$, ${{\boldsymbol{\Psi }}}_{{lm}}(\theta ,\varphi )=r{\rm{\nabla }}{Y}_{{lm}}$ and ${{\boldsymbol{\Phi }}}_{{lm}}(\theta ,\varphi )={{\boldsymbol{e}}}_{r}\times {{\boldsymbol{\Psi }}}_{{lm}}$ are vector spherical harmonics [45] and the radial functions ${v}_{{lm}}^{r}(r)$, ${v}_{{lm}}^{(1)}(r)$ and ${v}_{{lm}}^{(2)}(r)$ characterize the velocity field. The radial functions can be obtained by solving the linearized dynamic equations using the corresponding boundary conditions, see appendix A. The Stokes equation can be solved for a given shape perturbation independent of the concentration field so that the velocity field and pressure field is independent of the radial mode n. The radial part of the concentration field obeys a Helmholtz equation with an inhomogeneity that stems from hydrodynamic flows. The homogeneous part is solved by modified spherical Bessel functions and the inhomogeneous solution can be found using Greens functions. Using the dynamic equation for the shape changes of the droplet equation (10), we obtain an equation for the eigenvalue μ_nlm,

$$\begin{eqnarray}&&{\mu }_{{nlm}}=\displaystyle \frac{{v}_{l}^{r}(\bar{R})}{\bar{R}}+\displaystyle \frac{{D}_{+}}{{\rm{\Delta }}c}\left(\bar{c}^{\prime\prime} ({\bar{R}}_{+})+\displaystyle \frac{{c}_{{nl}}^{{\prime} }({\bar{R}}_{+})}{\bar{R}}\right)-\displaystyle \frac{{D}_{-}}{{\rm{\Delta }}c}\left(\bar{c}^{\prime\prime} ({\bar{R}}_{-})+\displaystyle \frac{{c}_{{nl}}^{{\prime} }({\bar{R}}_{-})}{\bar{R}}\right).\end{eqnarray} \tag{ 13 }$$

Here, the primes denote radial derivatives. Note that equation (13) is an implicit equation for the eigenvalues μ_nlm because the radial concentration modes c_nl(r) depend on μ_nlm, see appendix A. Equation (13) is independent of the index m, therefore the degeneracy of an eigenvalue μ_nl is at least $2l+1$. When all μ_nl are negative, the spherical shape is stable. The modes with l = 0 correspond to changes in droplet size without flows. They are always stable for $\bar{R}={\bar{R}}_{s}$ and unstable for $\bar{R}={\bar{R}}_{c}$. Thus droplet smaller than ${\bar{R}}_{c}$ will vanish and droplets larger will grow towards the size ${\bar{R}}_{s}$. Thus we consider the stability of $\bar{R}={\bar{R}}_{s}$ in the following. The modes with l = 1 do not involve shape deformations of the droplet and are thus not associated with flows. There always exists a marginal mode with μ_l=1 = 0 corresponding to overall translations where the droplet and all concentration fields are displaced and then stay in the new position. Here we consider shape instabilities for which a mode with l > 1 becomes unstable. Because shape deformations induce flows, this instability depends on the dimensionless viscosity F = wη₋/(γτ), as well as the ratio of viscosities in the two phases, η₊/η₋.

If we increase the supersaturation while keeping the other parameters fixed, the steady state can become unstable with respect to the mode l = 2 for a critical value  = _c. In figures 2(A)–(C), the onset of instability μ = 0 for the largest eigenvalue μ of the stationary radius is shown as a red dot, and unstable radii are indicated by red lines. Different lines correspond to different supersaturations, and the panels show different values of F. In figures 2(D)–(E), the corresponding stability diagrams of stationary droplets are shown as a function of the supersaturation and the reaction amplitude for different values of F. For large A and small , no stationary radius exists (white regions), so that any droplet would shrink and disappear. For large , the stationary state diverges (gray regions). Spherical droplets are stable in the blue regions. Stationary spherical droplets are unstable inside the red region, the surrounding black line marks the shape instability with respect to the l = 2 mode. The region where spherical droplets undergo a shape instability exists for $\epsilon \geqslant {\epsilon }^{* }$, which depends on F. The value of A for which the shape instability occurs at $\epsilon ={\epsilon }^{* }$ is denoted ${A}^{* }$, see figure 2(E).

For small A, the onset of instability can be describes by simple scaling behaviors. In the limit of small A and for $\epsilon \ll {\epsilon }_{\infty }$, we find ${\epsilon }^{* }\propto {F}^{-1/2}$ and $A={A}^{* }$ with ${A}^{* }\sim {F}^{-3/2}$ (compare figures 2(E)–(F)). For $A\lt {A}^{* }$, hydrodynamic flows govern the onset of instability which occurs at a value of A which behaves as A ∝ ⁻¹ F⁻² . For A > A*, hydrodynamic flows can be neglected as compared to diffusion fluxes and the onset of instability occurs for A ∝ ³. These two scaling regimes are indicated in in figures 2(D)–(F) by dashed lines. A derivation of these results including prefactors is given in appendix A.5.

We next address the question whether the shape instability found in the linear stability analysis can indeed give rise to droplet divisions in the presence of hydrodynamic flows in the nonlinear regime of the dynamics. We use a Cahn–Hilliard model [46] for phase separation dynamics, extended to include chemical reactions and hydrodynamic flows, that can capture topological changes of the interface. We include chemical reactions via a source term linear in the concentration as well as advection by the hydrodynamic flow which is described by the incompressible Stokes equation. Using a semi-spectral method [47], we obtain numerical solutions in a cubic box with no-flux boundary conditions, see appendix B.

Starting from a weakly deformed spherical droplet, we find regimes where the droplet disappears, where it relaxes to a stable spherical shape and where it undergoes a shape instability, consistent with the linear stability analysis of the effective droplet model. The transitions between these regimes occur for parameter values close to those predicted by the linear stability analysis. In the unstable regime, droplets typically divide. This shows that the droplet division reported previously can also occurs in the presence of hydrodynamic flows. Figure 3 shows snapshots of the droplet shape together with corresponding hydrodynamic flow fields on the symmetry plane of a dividing droplet at different times. At early times when the droplet deformation is weak, the flow field is similar to the l = 2 mode obtained from the linear theory, figure 3(A). As the droplet elongates and its waistline shrinks, the flow field becomes more complex, see figures 3(B), (C). The flow field shown in figure 3(C) exhibits two additional vortex lines that form rings around the axis of rotational symmetry. Similarly, after division, two further vortex rings occur, see figure 3(D). Interestingly, for small deformations the hydrodynamic flow direction opposes the directions of interface motion at the main droplet axes, see figures 3(A), (B). For larger deformations at later times the flow switches its direction along the long droplet axis where it assists interface motion. At the waistline, the flow velocity becomes small, see figure 3(C). After division, the flow field between the daughter droplets has very small magnitude, while strong flows at the outer sides move the droplets apart figure 3(D).

Zoom In Zoom Out Reset image size

This example shows that division of active droplets can occur even if hydrodynamic flows that oppose division are taken into account. Because flows act in opposition to the initial deformation of the sphere, the linear stability analysis already provides the key information of whether droplet division can occur for a given value of dimensionless viscosity F, see figure 2. This raises the question under what experimental conditions active droplets would become unstable and division could be observed. Ignoring hydrodynamic flows, $F\to \infty$ , it was shown that oil–water droplets and soft colloidal liquids or p-granules with sizes of a few micrometers could divide in the presence of chemical reactions [39]. To address the influence of hydrodynamic flows, we have to estimate the dimensionless viscosity $F=w{\eta }_{-}/(\gamma \tau )\simeq {k}_{B}T/(6\pi \gamma {wa})$, where we have used τ = w²/D and D ≃ k_BT/(6πηa) with molecular radius a. Thus, F is an equilibrium property of the phase separating fluid. For an oil–water system, we estimate F ≈ 0.1 , see appendix C. For soft colloidal liquids or p-granules, we estimate values between F ≈ 10–10⁴ . We can discuss these values using the stability diagrams in figures 2(D)–(F). Oil–water like droplets with F ≈ 0.1 are unlikely to divide, as the unstable region in the stability diagram is very narrow. For soft colloidal systems with F ≈ 10–10⁴ , droplet division might be experimentally observable. We can estimate typical reaction rates required for division to occur based on the reaction rate A* for which the range of supersaturation is maximal. The value of A* corresponds to a reaction rate in the droplet of the order of ν₋ = 10⁻⁴ mM s^–1, see appendix C. A comparison with reported enzymatic reaction rates [48] suggests that such values can be achieved in real systems.

We have shown that spontaneous division of chemically active droplets involves mechanical work against surface tension as droplets deform. Active droplets thus can transduce chemical energy to mechanical work and droplet division is therefore a mechano-chemical process. The surface tension of the droplet creates pressure gradients as the droplet becomes non-spherical that lead to hydrodynamic flows. Because the flows generated act against the shape deformation, droplets divide only for sufficiently large viscosity or sufficiently small surface tension and sufficiently large reaction rates. We show that the dependence of the onset of stability on parameters is captured for small reaction fluxes by simple scaling relations. Our work shows that droplet division would be suppressed in oil–water systems due to large surface tension and low viscosity. However it could be realized in soft colloidal systems for chemical reaction parameters that could be achieved experimentally. Furthermore flux-driven droplet divisions could be observable in biological systems, as both chemical reactions and phase-separating membrane-less organelles with low surface tensions can be found within cells.

Appendix A.: Effective droplet model with hydrodynamic flows

A.1. Stationary state of a spherical active droplet

Here, we discuss stationary solutions to equations (1)–(10) in the main text with spherical symmetry and without hydrodynamic flows $\bar{{\boldsymbol{v}}}=0$, where the bar indicates a steady state value. In this case, the pressure is constant both inside and outside the droplet, with a pressure difference due to Laplace pressure between the inside and outside of the droplets,

$$\begin{eqnarray}&&{\bar{p}}_{-}={\bar{p}}_{+}+\displaystyle \frac{2\gamma }{\bar{R}}.\end{eqnarray} \tag{ A.1 }$$

The steady state concentration profiles in the presence of chemical reactions are given by [39],

$$\begin{eqnarray}&&{\bar{c}}_{+}(r)=+\displaystyle \frac{{\nu }_{+}}{{k}_{+}}+{c}_{+}^{(0)}+{A}_{+}{k}_{0}(r/{l}_{+})\,\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{\bar{c}}_{-}(r)=-\displaystyle \frac{{\nu }_{-}}{{k}_{-}}+{c}_{-}^{(0)}+{A}_{-}{i}_{0}(r/{l}_{-}),\end{eqnarray} \tag{ A.3 }$$

where ${i}_{0}(x)=2\sinh (x)/x$ and k₀(x) = e^−x/x denote modified spherical Bessel functions of order zero of the first and second kind, respectively. The characteristic length scales ${l}_{\pm }={({D}_{\pm }/{k}_{\pm })}^{1/2}$ are set by reaction rate constants and diffusion coefficients. The parameters A_± are determined by the boundary condition at the droplet interface, equations (8) and (9)in the main text,

$$\begin{eqnarray}&&{A}_{+}=\left(\displaystyle \frac{\gamma {\beta }_{+}}{\bar{R}}-\displaystyle \frac{{\nu }_{+}}{{k}_{+}}\right)\displaystyle \frac{1}{{k}_{0}(\bar{R}/{l}_{+})}\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{A}_{-}=\left(\displaystyle \frac{\gamma {\beta }_{-}}{\bar{R}}+\displaystyle \frac{{\nu }_{-}}{{k}_{-}}\right)\displaystyle \frac{1}{{i}_{0}(\bar{R}/{l}_{+})}.\end{eqnarray} \tag{ A.5 }$$

Stationarity of the droplet radius $\bar{R}$ implies

$$\begin{eqnarray}&&{D}_{+}{\bar{c}}_{+}^{{\prime} }(\bar{R})={D}_{-}{\bar{c}}_{-}^{{\prime} }(\bar{R}),\end{eqnarray} \tag{ A.6 }$$

see equation (10) in the main text. Note that this equation typically has zero, one or two solutions for a given set of parameters.

A.2. Linearized dynamics

We introduce small perturbations to the spherically symmetric stationary state, with $p=\bar{p}+\delta p$, ${\boldsymbol{v}}=\delta {\boldsymbol{v}}$, $c=\bar{c}+\delta c$ and $R=\bar{R}+\delta R$ and write the dynamics of these perturbations to linear order. The linearized dynamics reads

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\delta p={\eta }_{\pm }{\rm{\Delta }}\delta {\boldsymbol{v}}\,\end{eqnarray} \tag{ A.7 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \delta {\boldsymbol{v}}=0\,\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&{\partial }_{t}\delta c=-\delta {v}_{r}\bar{c}^{\prime} +{D}_{\pm }{{\rm{\nabla }}}^{2}\delta c-{k}_{\pm }\delta c\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}\begin{array}{rcl}\,{\partial }_{t}\delta R & = & \delta {v}_{r}(\bar{R})+\displaystyle \frac{1}{{\rm{\Delta }}c}[{D}_{+}{\bar{c}}_{+}^{{\prime\prime} }(\bar{R})-{D}_{-}{\bar{c}}_{-}^{{\prime\prime} }(\bar{R})]\delta R\\ & & +\displaystyle \frac{1}{{\rm{\Delta }}c}[{D}_{+}{\partial }_{r}\delta {c}_{+}(\bar{R})-{D}_{-}{\partial }_{r}\delta {c}_{-}(\bar{R})].\end{array}\end{eqnarray} \tag{ A.10 }$$

Here δv_r denotes the radial part of the hydrodynamic velocity. With δc₋ and δc₊ we denote perturbations of the concentration field inside and outside the droplet. The same notation holds for the other fields. In this linear analysis, boundary conditions apply at the stationary radius $\bar{R}$,

$$\begin{eqnarray}&&\delta {c}_{\pm }(\bar{R})={\beta }_{\pm }\gamma \delta H-{\bar{c}}_{\pm }^{{\prime} }(\bar{R})\delta R,\end{eqnarray} \tag{ A.11 }$$

with perturbation of the curvature $\delta H=H(R)-H(\bar{R})$.

The linearized dynamics can be decomposed in spherical harmonics, see equation (11) in the main text. The curvature perturbation then takes the form

$$\begin{eqnarray}&&\delta H=\displaystyle \sum _{{nlm}}\displaystyle \frac{{h}_{l}}{\bar{R}}{\epsilon }_{{nlm}}{Y}_{{lm}},\end{eqnarray} \tag{ A.12 }$$

with h_l = (l² + l − 2)/2.

A.3. Hydrodynamic eigenmodes of the linearized dynamics

We can expand the hydrodynamic eigenmodes using a basis of vector spherical harmonics, see equation (12) in the main text. The velocity boundary conditions equation (7) in the main text for the mode amplitudes read

$$\begin{eqnarray}&&0={v}_{{lm}}^{r+}(\bar{R})-{v}_{{lm}}^{r-}(\bar{R})\,\end{eqnarray} \tag{ A.13 }$$

$$\begin{eqnarray}&&0={v}_{{lm}}^{(1)+}(\bar{R})-{v}_{{lm}}^{(1)-}(\bar{R})\,\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}&&0={v}_{{lm}}^{(2)+}(\bar{R})-{v}_{{lm}}^{(2)-}(\bar{R}).\end{eqnarray} \tag{ A.15 }$$

The stress boundary conditions (see equations (5) and (6) in the main text) at the interface read

$$\begin{eqnarray}&&\,0=2{\eta }_{+}({v}_{{lm}}^{r+})^{\prime} (\bar{R})-{p}_{{lm}}^{+}(\bar{R})-2{\eta }_{-}({v}_{{lm}}^{r-})^{\prime} (\bar{R})+{p}_{{lm}}^{-}(\bar{R})-2\gamma {\epsilon }_{{lm}}\displaystyle \frac{{h}_{l}}{\bar{R}}\end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray}&&0={\eta }_{+}\left[({v}_{{lm}}^{(1)+})^{\prime} (\bar{R})+\displaystyle \frac{{v}_{{lm}}^{r+}(\bar{R})}{\bar{R}}-\displaystyle \frac{{v}_{{lm}}^{(1)+}(\bar{R})}{\bar{R}}\right]\,\end{eqnarray} \tag{ A.17 }$$

$$\begin{eqnarray}&&-{\eta }_{-}\left[({v}_{{lm}}^{(1)-})^{\prime} (\bar{R})+\displaystyle \frac{{v}_{{lm}}^{r-}(\bar{R})}{\bar{R}}-\displaystyle \frac{{v}_{{lm}}^{(1)-}(\bar{R})}{\bar{R}}\right]\end{eqnarray} \tag{ A.18 }$$

$$\begin{eqnarray}&&0={\eta }_{+}\left[({v}_{{lm}}^{(2)+})^{\prime} (\bar{R})-\displaystyle \frac{{v}_{{lm}}^{(2)+}(\bar{R})}{\bar{R}}\right]-{\eta }_{-}\left[({v}_{{lm}}^{(2)-})^{\prime} (\bar{R})-\displaystyle \frac{{v}_{{lm}}^{(2)-}(\bar{R})}{\bar{R}}\right].\phantom{\rule{-2.0em}{0ex}}\end{eqnarray} \tag{ A.19 }$$

We solve the radial profiles of the modes with a polynomial ansatz and exclude functions that diverge for r → 0 or $r\to \infty$ inside and outside the droplet, respectively. The pressure is then given by

$$\begin{eqnarray}&&{p}_{{lm}}^{-}(r)=\gamma {f}_{A}{\left(\displaystyle \frac{r}{R}\right)}^{l+1}\,\end{eqnarray} \tag{ A.20 }$$

$$\begin{eqnarray}&&{p}_{{lm}}^{+}(r)=-\gamma {f}_{B}{\left(\displaystyle \frac{r}{R}\right)}^{-l}.\end{eqnarray} \tag{ A.21 }$$

For the hydrodynamic flow velocity we obtain

$$\begin{eqnarray}&&{v}_{{lm}}^{r-}(r)=\displaystyle \frac{\gamma }{{\eta }_{-}}\left[{f}_{C1}{\left(\displaystyle \frac{r}{R}\right)}^{l+1}-{f}_{C3}{\left(\displaystyle \frac{r}{R}\right)}^{l-1}\right]\end{eqnarray} \tag{ A.22 }$$

$$\begin{eqnarray}&&{v}_{{lm}}^{(1)-}(r)=\displaystyle \frac{\gamma }{{\eta }_{-}}\left[\displaystyle \frac{l+3}{l(l+1)}{f}_{C1}{\left(\displaystyle \frac{r}{R}\right)}^{l+1}-\displaystyle \frac{l+1}{l(l+1)}{f}_{C3}{\left(\displaystyle \frac{r}{R}\right)}^{l-1}\right]\,\end{eqnarray} \tag{ A.23 }$$

$$\begin{eqnarray}&&{v}_{{lm}}^{(2)-}(r)=0\,\end{eqnarray} \tag{ A.24 }$$

and

$$\begin{eqnarray}&&{v}_{{lm}}^{r+}(r)=\displaystyle \frac{\gamma }{{\eta }_{-}}\left[-{f}_{C2}{\left(\displaystyle \frac{r}{R}\right)}^{-l}+{f}_{C4}{\left(\displaystyle \frac{r}{R}\right)}^{-l-2}\right]\end{eqnarray} \tag{ A.25 }$$

$$\begin{eqnarray}&&{v}_{{lm}}^{(1)+}(r)=\displaystyle \frac{\gamma }{{\eta }_{-}}\left[\displaystyle \frac{l-2}{l(l+1)}{f}_{C2}{\left(\displaystyle \frac{r}{R}\right)}^{-l}-\displaystyle \frac{1}{l+1}{f}_{C4}{\left(\displaystyle \frac{r}{R}\right)}^{-l-2}\right]\,\end{eqnarray} \tag{ A.26 }$$

$$\begin{eqnarray}&&{v}_{{lm}}^{(2)+}(r)=0.\,\end{eqnarray} \tag{ A.27 }$$

Here, we have defined

$$\begin{eqnarray}&&{f}_{A}=\displaystyle \frac{(l-1)(l+1)(l+2)(2l+3)}{{\rm{\Delta }}(2{l}^{2}+4l)+(2{l}^{2}+4l+3)}\,\end{eqnarray} \tag{ A.28 }$$

$$\begin{eqnarray}&&{f}_{B}=\displaystyle \frac{l(l-1)(l+2)(2l-1)}{(2{l}^{2}+1)+(2{l}^{2}-2)/{\rm{\Delta }}}\,\end{eqnarray} \tag{ A.29 }$$

$$\begin{eqnarray}&&{f}_{C1}=\displaystyle \frac{1}{2}\displaystyle \frac{l(l-1)(l+1)(l+2)}{{\rm{\Delta }}(2{l}^{2}+4l)+(2{l}^{2}+4l+3)}\,\end{eqnarray} \tag{ A.30 }$$

$$\begin{eqnarray}&&{f}_{C2}=\displaystyle \frac{1}{2}\displaystyle \frac{l(l-1)(l+1)(l+2)}{{\rm{\Delta }}(2{l}^{2}+1)+(2{l}^{2}-2)}\,\end{eqnarray} \tag{ A.31 }$$

$$\begin{eqnarray}&&{f}_{C3}=\displaystyle \frac{1}{2}\displaystyle \frac{l(l-1)(l+1)(l+2)({\rm{\Delta }}(2{l}^{2}+4l+3)+(2{l}^{2}+4l))}{({\rm{\Delta }}(2{l}^{2}+1)+(2{l}^{2}-2))({\rm{\Delta }}(2{l}^{2}+4l)+(2{l}^{2}+4l+3))}\,\,\end{eqnarray} \tag{ A.32 }$$

$$\begin{eqnarray}&&{f}_{C4}=\displaystyle \frac{1}{2}\displaystyle \frac{l(l-1)(l+1)(l+2)({\rm{\Delta }}(2{l}^{2}-2)+(2{l}^{2}+1))}{({\rm{\Delta }}(2{l}^{2}+1)+(2{l}^{2}-2))({\rm{\Delta }}(2{l}^{2}+4l)+(2{l}^{2}+4l+3))},\end{eqnarray} \tag{ A.33 }$$

where Δ = η₊/η₋ denotes the ratio of the viscosities inside and outside the droplet.

A.4. Concentration eigenmodes

The equation for the radial part of the concentration eigenmode is

$$\begin{eqnarray}&&\displaystyle \frac{1}{{D}_{\pm }}{v}_{l}^{r}(r)\bar{c}^{\prime} (r)=\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}}{r}^{2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}}-{\lambda }_{{nl}}^{\pm 2}-\displaystyle \frac{l(l+1)}{{r}^{2}}\right]{c}_{{nl}}(r)\end{eqnarray} \tag{ A.34 }$$

with

$$\begin{eqnarray}&&{\lambda }_{{nl}}^{\pm 2}=({k}_{\pm }+{\mu }_{{nl}})/{D}_{\pm }.\end{eqnarray} \tag{ A.35 }$$

The boundary conditions at $\bar{R}$ are

$$\begin{eqnarray}&&{c}_{{nl}}({\bar{R}}_{\pm })=\gamma {\beta }_{\pm }\displaystyle \frac{{h}_{l}}{\bar{R}}-\bar{R}\,\bar{c}^{\prime} ({\bar{R}}_{\pm }).\end{eqnarray} \tag{ A.36 }$$

The left-hand side of equation (A.34) constitutes an inhomogeneity

$$\begin{eqnarray}&&{f}_{l}^{\pm }(r)=-\displaystyle \frac{1}{{D}_{\pm }}{v}_{l}^{r}(r)\bar{c}^{\prime} (r).\end{eqnarray} \tag{ A.37 }$$

The solution ${c}_{{nl}}^{\pm }(r)$ of the inhomogeneous equation (A.34) that satisfies the boundary condition equation (A.36) can be constructed from a particular solution ${c}_{{nl},p}^{\pm }(r)$ of the inhomogeneous equation to which solutions ${c}_{{nl},h}^{\pm }(r)$ of the homogeneous equation with ${f}_{l}^{\pm }=0$ are added to satisfy the boundary conditions, equation (A.36). This can be expressed as

$$\begin{eqnarray}&&{c}_{{nl}}^{-}(r)={\alpha }_{{nl}}^{-}{c}_{{nl},h}^{-}(r)+{c}_{{nl},p}^{-}(r)\,\end{eqnarray} \tag{ A.38 }$$

$$\begin{eqnarray}&&{c}_{{nl}}^{+}(r)={\alpha }_{{nl}}^{+}{c}_{{nl},h}^{+}(r)+{c}_{{nl},p}^{+}(r),\end{eqnarray} \tag{ A.39 }$$

where the coefficients α_± read

$$\begin{eqnarray}&&{\alpha }_{{nl}}^{\pm }=\displaystyle \frac{{a}_{l}^{\pm }-{c}_{{nl},p}^{\pm }(\bar{R})}{{c}_{{nl},h}^{\pm }(\bar{R})},\end{eqnarray} \tag{ A.40 }$$

with ${a}_{l}^{\pm }={c}_{{nl}}^{\pm }(\bar{R})$.

We are especially interested in the case of unstable modes with μ_nl > 0. Therefore we focus on the solution of equation (A.34) for ${\lambda }_{{nl}}^{\pm 2}\gt 0$ and k_± > 0. In this case, the homogeneous equation with ${f}_{l}^{\pm }\,=\,0$ is a modified Helmholtz equation which is solved by modified spherical Bessel functions, ${c}_{{nl},h}^{-}(r)={i}_{l}({\lambda }_{{nl}}^{-}r)$ and ${c}_{{nl},h}^{+}(r)={k}_{l}({\lambda }_{{nl}}^{+}r)$, where i_l and k_l denote the modified spherical Bessel functions of first and second order, respectively. The particular solution of the inhomogeneous equation can be obtained by a Green's function approach,

$$\begin{eqnarray}\begin{array}{rcl}{c}_{l,p}^{-}(r) & = & {\lambda }_{{nl}}^{-}{k}_{l}({\lambda }_{{nl}}^{-}r){\displaystyle \int }_{0}^{r}[{i}_{l}({\lambda }_{{nl}}^{-}{r}_{2}){f}_{l}^{-}({r}_{2}){r}_{2}^{2}]{{\rm{d}}{r}}_{2}\\ & & +{\lambda }_{{nl}}^{-}{i}_{l}({\lambda }_{{nl}}^{-}r){\displaystyle \int }_{r}^{\bar{R}}[{k}_{l}({\lambda }_{{nl}}^{-}{r}_{2}){f}_{l}^{-}({r}_{2}){r}_{2}^{2}]{{\rm{d}}{r}}_{2}\end{array}\end{eqnarray} \tag{ A.41 }$$

$$\begin{eqnarray}\begin{array}{rcl}\,{c}_{l,p}^{+}(r) & = & {\lambda }_{{nl}}^{+}{k}_{l}({\lambda }_{{nl}}^{+}r){\displaystyle \int }_{\bar{R}}^{r}[{i}_{l}({\lambda }_{{nl}}^{+}{r}_{2}){f}_{l}^{+}({r}_{2}){r}_{2}^{2}]{{\rm{d}}{r}}_{2}\\ & & +{\lambda }_{{nl}}^{+}{i}_{l}({\lambda }_{{nl}}^{+}r){\displaystyle \int }_{r}^{\infty }[{k}_{l}({\lambda }_{{nl}}^{+}{r}_{2}){f}_{l}^{+}({r}_{2}){r}_{2}^{2}]{{\rm{d}}{r}}_{2},\end{array}\end{eqnarray} \tag{ A.42 }$$

with the radial part of the inhomogeneity ${f}_{l}^{\pm }(r)$ given by equation (A.37). The explicit calculation of these functions has to be handled with care, since the functions k_l and i_l have divergences for large and small arguments r that cancel in the final result but can still lead to numerical difficulties when evaluated directly.

The derivative of the concentration profile at $\bar{R}$ can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{c}_{{nl}}^{{\prime} }({\bar{R}}_{-}) & = & \displaystyle \frac{{a}_{l}^{-}}{\bar{R}}{g}_{l,i}({\lambda }_{{nl}}^{-}\bar{R})\\ & & +\displaystyle \frac{{c}_{l,p}^{-}(\bar{R})}{\bar{R}}\cdot [{g}_{l,k}({\lambda }_{{nl}}^{-}\bar{R})-{g}_{l,i}({\lambda }_{{nl}}^{-}\bar{R})]\end{array}\end{eqnarray} \tag{ A.43 }$$

$$\begin{eqnarray}\begin{array}{rcl}\,{c}_{{nl}}^{{\prime} }({\bar{R}}_{+}) & = & \displaystyle \frac{{a}_{l}^{+}}{\bar{R}}{g}_{l,k}({\lambda }_{{nl}}^{+}\bar{R})\\ & & +\displaystyle \frac{{c}_{l,p}^{+}(\bar{R})}{\bar{R}}\cdot [{g}_{l,i}({\lambda }_{{nl}}^{+}\bar{R})-{g}_{l,k}({\lambda }_{{nl}}^{+}\bar{R})],\end{array}\end{eqnarray} \tag{ A.44 }$$

with

$$\begin{eqnarray}&&{g}_{l,i}(x)=\displaystyle \frac{{{xi}}_{l}^{{\prime} }(x)}{{i}_{l}(x)}\end{eqnarray} \tag{ A.45 }$$

$$\begin{eqnarray}&&\,{g}_{l,k}(x)=\displaystyle \frac{{{xk}}_{l}^{{\prime} }(x)}{{k}_{l}(x)}.\end{eqnarray} \tag{ A.46 }$$

Using the equation for the shape perturbations (A.10), and using equations (A.43) and (A.44), we obtain equation (13) in the main text. This equation determines the eigenvalue μ_nlm of the hydrodynamic modes.

A.5. Scaling relations in the limit of small reaction fluxes

In the limit of small chemical reaction fluxes s_± we obtain simple scaling expressions for stationary radii and their shape instability conditions. Here we present the method and discuss the results.

A.5.1. Stationary radius

Here we discuss the stationary radius in the limit of small chemical reaction amplitude A = ν₋τ/Δc while keeping the ratios ν₋/(k₋Δc) and k₊/k₋ of reaction parameters fixed. This corresponds to the curves $\bar{R}(\epsilon )$ shown in figure 2(A) for different values of A. We can identify two regimes in the figure. The first is the region of small ,  ∼ ₀, which corresponds to the minimum of $\epsilon (\bar{R})$. The second is the region of ${\epsilon }_{\infty }$ where the stationary radius diverges. For A → 0, we see that ₀ goes to zero while ${\epsilon }_{\infty }$ stays constant, and both are connected by a straight line that indicates scaling behavior of $\bar{R}={\bar{R}}_{s}$. This increasing separation between ₀ and ${\epsilon }_{\infty }$ (and the corresponding stationary radii) in the limit of small A means that we can analyze the behavior of the stationary radius in these two regimes separately. For this we consider equations (A.2) and (A.3) for the concentration field and (A.6) for the stationary radius. We can rewrite (A.6) to obtain an expression relating the supersaturation to the stationary radius,

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{\beta }_{+}\gamma }{{\rm{\Delta }}c\bar{R}}+\left(\displaystyle \frac{{\beta }_{-}\gamma }{{\rm{\Delta }}c\bar{R}}+\displaystyle \frac{{\nu }_{-}}{{k}_{-}{\rm{\Delta }}c}\right)\displaystyle \frac{{D}_{-}}{{D}_{+}}\displaystyle \frac{\tfrac{\bar{R}}{{l}_{-}}\coth \tfrac{\bar{R}}{{l}_{-}}-1}{1+\tfrac{\bar{R}}{{l}_{+}}}.\end{eqnarray} \tag{ A.47 }$$

In this limit of small A, the characteristic length-scales of the concentration field become large with ${l}_{\pm }\propto {A}^{-1/2}$. To find scaling regimes in equation (A.47), we change variables in equation (A.47) from $(A,\bar{R})$ to $(A,\hat{R})$ with $\hat{R}=\bar{R}{A}^{a}/w$, where a is an exponent. For a = 1/3 we find the behavior of (R) close to ₀ and ${\bar{R}}_{0}$,

$$\begin{eqnarray}&&\hat{\epsilon }=\displaystyle \frac{1}{6}{\hat{R}}^{-1}+\displaystyle \frac{1}{3}{\hat{R}}^{2}+O({A}^{1/6}),\end{eqnarray} \tag{ A.48 }$$

where $\hat{\epsilon }=\epsilon {A}^{-1/3}$ becomes independent of A for small A. This function describes the supersaturation as a function of radius around the threshold value ₀. Due to the inverted presentation $\epsilon (\bar{R})$ instead of $\bar{R}(\epsilon )$ the function captures both the nucleation radius ${\bar{R}}_{c}$ and the larger radius ${\bar{R}}_{s}$. The threshold value ₀ can be obtained from equation (A.48) by minimizing $\hat{\epsilon }$ for fixed A as $\partial \hat{\epsilon }/\partial \hat{R}=0$. It behave as

$$\begin{eqnarray}&&{\epsilon }_{0}={4}^{-2/3}{A}^{1/3}+O({A}^{1/2}).\end{eqnarray} \tag{ A.49 }$$

For large and small $\hat{R}$, equation (A.48) describes the steady radii ${\bar{R}}_{s}$ and ${\bar{R}}_{c}$, respectively, for which $\epsilon \geqslant {\epsilon }_{0}$. For large , the critical radius obeys

$$\begin{eqnarray}&&{\bar{R}}_{c}\simeq \displaystyle \frac{w}{6\epsilon },\end{eqnarray} \tag{ A.50 }$$

while the larger stationary radius is

$$\begin{eqnarray}&&{\bar{R}}_{s}\simeq w{(3\epsilon A)}^{1/2}.\end{eqnarray} \tag{ A.51 }$$

In figure 2(B), the scaling behaviors given by equations (A.51) and (A.50) are indicated by dashed lines. At  = ₀ both radii meet at $\bar{R}={\bar{R}}_{0}$, where

$$\begin{eqnarray}&&{\bar{R}}_{0}=w{(4A)}^{-1/3}+O({A}^{-1/2}).\end{eqnarray} \tag{ A.52 }$$

For $a=1/2$, $\bar{R}/{l}_{\pm }$ becomes independent of A and

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{D}_{-}}{{D}_{+}}\displaystyle \frac{{\nu }_{-}}{{k}_{-}{\rm{\Delta }}c}\displaystyle \frac{\tfrac{\bar{R}}{{l}_{-}}\coth \left(\tfrac{\bar{R}}{{l}_{-}}\right)-1}{1+\tfrac{\bar{R}}{{l}_{+}}}+O({A}^{1/2}).\end{eqnarray} \tag{ A.53 }$$

For $\bar{R}/{l}_{\pm }\ll 1$, the stationary radius obeys equation (A.51) and is thus the larger stationary radius ${\bar{R}}_{s}$. For $\bar{R}/{l}_{\pm }\gg 1$, we obtain the divergence of ${\bar{R}}_{s}$ as approaches ${\epsilon }_{\infty }$ with

$$\begin{eqnarray}&&{\epsilon }_{\infty }=\sqrt{\displaystyle \frac{{D}_{-}{k}_{-}}{{D}_{+}{k}_{+}}}\displaystyle \frac{{\nu }_{-}}{{k}_{-}{\rm{\Delta }}c}.\end{eqnarray} \tag{ A.54 }$$

A.5.2. Shape instability

We now discuss scaling relations for the onset of instability in the $(A,\hat{\epsilon })$ plane in the limit of small A, which give the trends shown as dashed lines in figures 2(D)–(F). We use the scaling of the stationary radius $\bar{R}={\bar{R}}_{s}$ close to ₀ with $\hat{R}=\bar{R}{A}^{1/3}/w$, $\hat{\epsilon }=\epsilon {A}^{-1/3}$ and ${\hat{l}}_{\pm }={l}_{\pm }{A}^{1/2}$ in equation (13) to obtain

$$\begin{eqnarray}&&{\hat{\mu }}_{{nlm}}=-\displaystyle \frac{{d}_{l}}{\hat{R}}\displaystyle \frac{{A}^{-2/3}}{F}+\displaystyle \frac{2}{3}(l-1)-\displaystyle \frac{{D}_{+}}{{D}_{-}}\displaystyle \frac{(l-1){g}_{l}}{{\hat{R}}^{3}}+O({A}^{1/6}),\end{eqnarray} \tag{ A.55 }$$

where ${\hat{\mu }}_{{nlm}}={\mu }_{{nlm}}\tau /A$ and $\hat{R}$ is related to $\hat{\epsilon }$ by (A.48). Here, ${d}_{l}={f}_{C3}-{f}_{C1}$, where f_C1 and f_C3 are defined in equation (A.33) and

$$\begin{eqnarray}&&{g}_{l}=\displaystyle \frac{{h}_{l}(l+1)+\tfrac{{D}_{-}}{{D}_{+}}\tfrac{{\beta }_{-}}{{\beta }_{+}}{h}_{l}l}{l-1}\end{eqnarray} \tag{ A.56 }$$

with h_l = (l² + l − 2)/2. For large mode index l,

$$\begin{eqnarray}&&{d}_{l}=\displaystyle \frac{l}{2({\eta }_{+}/{\eta }_{-}+1)}+O(1/l).\end{eqnarray} \tag{ A.57 }$$

We now consider conditions for which ${\mu }_{{nlm}}=0$ for small A and the mode (n, l, m) becomes unstable. Using (A.48) in (A.55), we find a relation between $\hat{\epsilon }$ and $\hat{R}$ at the onset of instability μ_nlm = 0,

$$\begin{eqnarray}&&\hat{\epsilon }=\displaystyle \frac{{d}_{l}}{2(l-1)}\displaystyle \frac{1}{\hat{F}}\hat{R}+\left(\displaystyle \frac{1}{6}+\displaystyle \frac{{D}_{+}}{{D}_{-}}\displaystyle \frac{1}{2}{g}_{l}\right){\hat{R}}^{-1}+O({A}^{1/6}).\end{eqnarray} \tag{ A.58 }$$

This curve captures the scaling behavior of the onset of instability for different parameters in the $\bar{R}-\epsilon$ plane, corresponding to the red dotted line in figures 2(A)–(C).

We now focus on finding the scaling relations for the onset of stability of the stationary radius as function of A, and F, as shown in figures 2(D)–(F). At this onset, both (A.58) and (A.48) need to be satisfied. We use both equations to eliminate $\hat{R}$. We find a crossover regime with relations ${A}^{* }\sim {F}^{-3/2}$ between the region where hydrodynamic flows are relevant (A < A*) and where they can be neglected (A > A*). For A > A* we find for μ_nlm = 0 as relation between A and

$$\begin{eqnarray}&&A\simeq 54\displaystyle \frac{{g}_{l}}{{\left(1+\tfrac{1}{2}{g}_{l}\right)}^{3}}{\epsilon }^{3}.\end{eqnarray} \tag{ A.59 }$$

For A < A* we find

$$\begin{eqnarray}&&A\simeq \displaystyle \frac{1}{3}{\left(\displaystyle \frac{2(l-1)}{{d}_{l}}\right)}^{2}{\epsilon }^{-1}{F}^{-2}.\end{eqnarray} \tag{ A.60 }$$

In figures 2(D)–(F), the dashed lines indicate these two scaling solutions in the limit A → 0 and $F\to \infty$ for l = 2, which we find to be the first mode to become unstable. We find that the general trends of the stability diagram is captured well, with small deviations from the full solution of equation (13) for small , and larger deviations in the regime close to ${\epsilon }_{\infty }$ where the scaling of the stationary radius ${\bar{R}}_{s}\propto {A}^{-1/3}$ breaks down.

Appendix B.: Continuum model for active droplets with flows

B.1. Continuum model for active droplets

We study an extended Cahn–Hilliard equation with chemical reactions coupled to Stokes equation for hydrodynamic flows at low Reynolds numbers. We consider an incompressible fluid containing two components A and B, with number concentration fields ${c}_{A}({\boldsymbol{r}},t)$ and $c={c}_{B}({\boldsymbol{r}},t)$ that depend on position ${\boldsymbol{r}}$ and time t, and with molecular masses m_A and m_B and molecular volumes v_A and v_B. We are interested in the case where component A forms the background fluid and B is a droplet material that forms droplets by phase separation. Additionally, chemical reactions convert the two components into each other, $A\rightleftharpoons B$. For simplicity, we consider mass and volume conserving chemical reactions with m_A/v_A = m_B/v_B, which encodes that volume is conserved in the reaction if mass is conserved. Together with incompressibility, this implies that the mass density ρ = m_Ac_A + m_Bc_B is constant. Therefore, we can describe the system by the concentration $c({\boldsymbol{r}},t)$ of the droplet material B only.

We use the following double-well free energy density [46]

$$\begin{eqnarray}&&f(c)=\displaystyle \frac{b}{2{({\rm{\Delta }}c)}^{2}}{(c-{c}_{-}^{(0)})}^{2}{(c-{c}_{+}^{(0)})}^{2}+\displaystyle \frac{\kappa }{2}{({\rm{\nabla }}c)}^{2},\end{eqnarray} \tag{ B.1 }$$

with ${\rm{\Delta }}c=| {c}_{-}^{(0)}-{c}_{+}^{(0)}|$. Here, the positive parameter b characterizes molecular interactions and entropic contributions. This free energy describes the segregation of the fluid in two coexisting phases [49]: one phase rich in droplet material with c ≈ c⁽⁰⁾₋ and a dilute phase with c ≈ c⁽⁰⁾₊. The coefficient κ is related to surface tension and the interface width [46].

The state of the system is characterized by the free energy

$$\begin{eqnarray}&&F[c]=\int {{\rm{d}}}^{3}r\,f(c),\end{eqnarray} \tag{ B.2 }$$

where the integral is over the system volume. We work with an ensemble T, ρ, c here, where T denotes temperature and the system is considered isothermal. The chemical potential $\bar{\mu }=\delta F[c]/\delta c$, governs demixing and can be split into local and nonlocal contributions, $\bar{\mu }={\bar{\mu }}_{0}-\kappa {{\rm{\nabla }}}^{2}c$ with

$$\begin{eqnarray}&&{\bar{\mu }}_{0}=\displaystyle \frac{b}{{({\rm{\Delta }}c)}^{2}}(c-{c}_{+}^{(0)})(c-{c}_{-}^{(0)})(2c-{c}_{-}^{(0)}-{c}_{+}^{(0)}).\end{eqnarray} \tag{ B.3 }$$

The dynamics of the concentration field is described by [50, 51]

$$\begin{eqnarray}&&{\partial }_{t}c=-{\rm{\nabla }}\cdot {\boldsymbol{j}}+s(c)\,\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{\boldsymbol{j}}=-M{\rm{\nabla }}\bar{\mu }+{\boldsymbol{v}}c.\end{eqnarray} \tag{ B.5 }$$

Here, M is a mobility coefficient of the droplet material and ${\boldsymbol{v}}$ is the hydrodynamic velocity. The source term s(c) describes chemical reactions, for which we choose for simplicity a linear concentration dependence,

$$\begin{eqnarray}&&s(c)=\nu -k(c-{c}_{+}^{(0)}).\end{eqnarray} \tag{ B.6 }$$

The reaction flux given in equation (B.6) does not obey detailed balance with respect to the free energy, and thus describes a situation where an external energy source maintains the system away from equilibrium [39].

The hydrodynamic velocity ${\boldsymbol{v}}$ can be calculated using momentum conservation,

$$\begin{eqnarray}&&{\partial }_{t}(\rho {v}_{\alpha })={\partial }_{\beta }{\sigma }_{\alpha \beta },\end{eqnarray} \tag{ B.7 }$$

with momentum ρv_α and stress tensor σ_αβ, where α and β number Cartesian coordinates x, y, z. We can write the stress tensor σ_αβ as

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }=-(\rho {v}_{\alpha }){v}_{\beta }+{\sigma }_{\alpha \beta }^{{\rm{eq}}}+{\sigma }_{\alpha \beta }^{d},\end{eqnarray} \tag{ B.8 }$$

where the first term describes advection of the stress tensor, ${\sigma }_{\alpha \beta }^{{\rm{eq}}}$ and ${\sigma }_{\alpha \beta }^{d}$ denote the equilibrium and dissipative stress tensors. The equilibrium stress tensor is given by

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }^{{\rm{eq}}}=-(\bar{\mu }c-f){\delta }_{\alpha \beta }-\displaystyle \frac{\partial f}{\partial ({\partial }_{\alpha }c)}{\partial }_{\beta }c-{P}_{0}{\delta }_{\alpha \beta }.\end{eqnarray} \tag{ B.9 }$$

Here, $\bar{\mu }c-f$ is the osmotic pressure of the droplet material, and δ_αβ denotes the Kronecker delta. Incompressibility is enforced by an additional partial pressure P₀. The deviatory stress tensor can be found as thermodynamic force related to momentum by writing the entropy production rate,

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }^{d}=2\eta \left({v}_{\alpha \beta }-\displaystyle \frac{1}{3}{v}_{\gamma \gamma }{\delta }_{\alpha \beta }\right)+\eta ^{\prime} {v}_{\gamma \gamma }{\delta }_{\alpha \beta },\end{eqnarray} \tag{ B.10 }$$

where η and $\eta ^{\prime}$ denote viscosities, and v_αβ = (∂_α v_β + ∂_β v_α)/2 is the symmetric strain tensor.

In the Stokes limit, the inertial terms are neglected, ${D}_{t}(\rho {v}_{\alpha })=0$, with advected derivative D_t = ∂_t + v_β∂_β, leaving $0={\partial }_{\beta }({\sigma }_{\alpha \beta }^{{\rm{eq}}}+{\sigma }_{\alpha \beta }^{d})$. This yields [52]

$$\begin{eqnarray}&&\eta {\partial }_{\beta }^{2}{v}_{\alpha }=3{\bar{\mu }}_{0}{\partial }_{\alpha }c-\kappa c{{\rm{\nabla }}}^{2}({\partial }_{\alpha }c)+{\partial }_{\alpha }{P}_{0}.\end{eqnarray} \tag{ B.11 }$$

Equations (B.3)–(B.6) and (B.11) and incompressibility ∂_αv_α = 0 define the continuum model of active droplets.

B.2. Numerical solution of the continuum model

We numerically solve the dynamic equations of the continuum model of active droplets, equations (B.3)–(B.6) and equation (B.11) with equation (B.13) and incompressibility ${\partial }_{\alpha }{v}_{\alpha }=0$.

For this we use a spectral method in a 3d rectangular box. This has the advantage that in a spectral decomposition, the spatial operators become simple multiplications with the wavenumber [47]. However, our equations contain a number of nonlinear functions, which are easier to evaluate in real space. We therefore transform forward and back in each time step.

To calculate the next timestep t_i from the fields found in timestep t_i−1, we use a semi-implicit Runge–Kutta method [53] (method (2, 3, 3)) for the concentration field. This evaluates the gradient term in $\bar{\mu }$, equation (B.3), implicitly, while evaluating the rest of $\bar{\mu }$ as well as the advection term of the fluxes, ${\boldsymbol{v}}c$, explicitly. This effectively means that the terms related to the interfacial profile are calculated implicitly, which allows for larger time steps as an explicit scheme.

For the concentration field, we choose no-flux boundary conditions (∂_nc = 0, where the derivative is in a direction normal to the simulation box), which leads to a decomposition in cosine functions in the spectral description. The Laplacian then is −k² for a mode with wave vector ${\boldsymbol{k}}$. The Stokes equation can also be solved using spectral methods. Here, no-flux conditions lead to v_n = 0. Additionally we enforce incompressibility using a reprojection method. For this, the velocity field calculated by neglecting the partial pressure, P_p = 0, can be split into two parts (Helmholtz decomposition),

$$\begin{eqnarray}&&{\boldsymbol{v}}={{\boldsymbol{v}}}_{\psi }+{{\boldsymbol{v}}}_{\phi }={\rm{\nabla }}\times {\boldsymbol{\psi }}-{\rm{\nabla }}\phi \end{eqnarray} \tag{ B.12 }$$

with vector field ${\boldsymbol{\psi }}$ and scalar field ϕ, and velocity parts ${{\boldsymbol{v}}}_{\psi }={\rm{\nabla }}\times {\boldsymbol{\psi }}$ and ${{\boldsymbol{v}}}_{\phi }=-{\rm{\nabla }}\phi$. With this, we find

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{v}}={\rm{\Delta }}\phi \end{eqnarray} \tag{ B.13 }$$

and thus, using incompressibility, ${\rm{\nabla }}\cdot {\boldsymbol{v}}=0$, we can calculate ϕ. We thus find the incompressible part of the velocity field

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{\psi }={\boldsymbol{v}}-{\rm{\nabla }}\phi .\end{eqnarray} \tag{ B.14 }$$

We can evaluate this in Fourier space using a spectral method. For a rectangular box aligned with the coordinate system, we thus find that each velocity component v_α is decomposed by sines in one direction and cosines in the other direction. Spatial derivatives convert a sine-description into cosines, and vice versa.

We normalize concentration, length, time and energy by ${\rm{\Delta }}c={c}_{-}^{(0)}-{c}_{+}^{(0)}$, $w=2{(\kappa /b)}^{1/2}$, t₀ = w²/D and ${\hat{e}}_{0}=\kappa \hat{w}{({\rm{\Delta }}c)}^{2}/3$, respectively, where the characteristic length scale is $w=2{(\kappa /b)}^{1/2}$. The relevant dimensionless model parameters are ${c}_{+}^{(0)}/{\rm{\Delta }}c$, kt₀ , and ν₋t₀/Δc. We choose ${c}_{+}^{(0)}/{\rm{\Delta }}c=0$, kt₀ = 10⁻², $\nu {t}_{0}=2\times {10}^{-3}$ and $\eta \,{\hat{w}}^{3}/({t}_{0}{\hat{e}}_{0})=2$. Additionally, we use as box-length $L/\hat{w}=100$ in all 3 dimensions, number of grid-points in one direction N = 128 and simulation time T/t₀ = 4 × 10³. For the time step, we start with a timestep of Δt/t₀ = 10⁻⁴, and double the timestep to a final step size of Δt/t₀ = 0.01.

We start with initial conditions R = R₀ (1 +  Y_2,0), with ${R}_{0}/\hat{w}=7$ and  = 1. The concentration field at positions ${\boldsymbol{r}}$ is initialized by the function

$$\begin{eqnarray}&&c({\boldsymbol{r}})=\displaystyle \frac{{c}_{+}^{(0)}+{c}_{-}^{(0)}}{2}+\displaystyle \frac{{c}_{+}^{(0)}-{c}_{-}^{(0)}}{2}\tanh \displaystyle \frac{d({\boldsymbol{r}})}{w},\end{eqnarray} \tag{ B.15 }$$

where $d({\boldsymbol{r}})$ is the oriented distance of ${\boldsymbol{r}}$ to the nearest point on the ellipsoid. The value of $d({\boldsymbol{r}})$ is negative for points inside the droplet and positive for points outside.

B.3. Effective droplet model as a limit of the continuum model

We now discuss the relationship between the effective droplet model and the continuum model. To relate the two models, we first use the continuum model to derive jump conditions for the concentration in the effective droplet model in equilibrium. We then consider stress balance across this interface and derive stress boundary conditions in the effective droplet model. Finally we discuss the dynamical equations in the bulk and at the interface in non-equilibrium situations.

B.3.1. Derivation of jump conditions for equilibrium phase separation

First we consider the phase separation in equilibrium without chemical reactions in the continuum model.

In a one-dimensional system with a mean concentration $\bar{c}$ with ${c}_{+}^{(0)}\lt \bar{c}\lt {c}_{-}^{(0)}$, the free energy of the system in equation (B.2) is minimized by the concentration profile

$$\begin{eqnarray}&&{c}^{* }(x)=\displaystyle \frac{{c}_{-}^{(0)}+{c}_{+}^{(0)}}{2}+\displaystyle \frac{{c}_{-}^{(0)}-{c}_{+}^{(0)}}{2}\tanh \displaystyle \frac{x}{w},\end{eqnarray} \tag{ B.16 }$$

where $w=2{(\kappa /b)}^{1/2}$ denotes the interfacial width and x is the normal distance to the interface. The concentration profile describes two phases of concentration ${c}_{-}^{(0)}$ and ${c}_{+}^{(0)}$ separated by a flat interface of width w. The surface tension can be defined as

$$\begin{eqnarray}&&\gamma ={\int }_{-\infty }^{\infty }F[{c}^{* }(x)]-\displaystyle \frac{1}{2}(F[{c}_{-}^{(0)}]+F[{c}_{+}^{(0)}]){\rm{d}}{x}.\end{eqnarray} \tag{ B.17 }$$

For the free energy equation (B.2) with the concentration profile equation (B.16), this can be written as $\gamma ={\int }_{-\infty }^{\infty }\kappa {({\rm{\nabla }}{c}^{* })}^{2}{\rm{d}}{x}$ which yields $\gamma ={({\rm{\Delta }}c)}^{2}/6\,\sqrt{\kappa b}$ or [54].

This interfacial tension governs the concentration jump condition in the effective droplet model, which can be derived as follows. To describe a curved interface, we consider two homogeneous phases with concentrations c_±. For a finite volume V_s with a droplet of size V and area A the concentrations c_± can be found by minimizing the free energy F = f(c₋) V + f(c₊)(V_s − V) + γA with $\partial F/\partial {c}_{-}{| }_{V}=0$ and $\partial F/\partial V{| }_{{c}_{-}}=0$, where the concentration of both phases are related by ${V}_{s}\bar{c}={{Vc}}_{-}+({V}_{s}-V){c}_{+}$ where $\bar{c}$ denotes the average concentration in the system. Thus for two phases to be in equilibrium, their chemical potential $\bar{\mu }$ and osmotic pressure ${\rm{\Pi }}=c\bar{\mu }-f$ need to obey

$$\begin{eqnarray}&&0=\bar{\mu }({c}_{-})-\bar{\mu }({c}_{+})\,\end{eqnarray} \tag{ B.18 }$$

$$\begin{eqnarray}&&0={\rm{\Pi }}({c}_{+})-{\rm{\Pi }}({c}_{-})-2\gamma H,\end{eqnarray} \tag{ B.19 }$$

where H the mean curvature of the droplet and $2\gamma H$ is the Laplace pressure. These equations determine the concentrations in the phases c_± of coexisting phases [54].

For small Laplace pressures, we can express the equilibrium concentrations c_± of a curved interface by the concentrations of a flat interface ${c}_{\pm }^{(0)}$ plus a small perturbation,

$$\begin{eqnarray}&&{c}_{-}={c}_{-}^{(0)}+{\beta }_{-}\gamma H\,\end{eqnarray} \tag{ B.20 }$$

$$\begin{eqnarray}&&{c}_{+}={c}_{+}^{(0)}+{\beta }_{+}\gamma H,\end{eqnarray} \tag{ B.21 }$$

where ${\beta }_{\pm }=2/(f^{\prime\prime} ({c}_{\pm }^{(0)}){\rm{\Delta }}c)$. For the free energy equation (B.2), we find β_± = 2/(bΔc), which is related to the interfacial width as w = 6γβ₊/Δc.

B.3.2. Stress balance across the interface

We now consider stress balance of the continuum model across the droplet interface to derive stress jump conditions at the interface in the effective droplet model. We discuss the mechanical equilibrium in a small volume across a curved interface with a local mean curvature H corresponding to a (local) effective radius $\tilde{R}=1/H$. We focus on the case where the interface is rotationally symmetric around the considered point ${\boldsymbol{R}}$, and where the curvature does not change along the interface. We use spherical coordinates, where the radial vector ${{\boldsymbol{e}}}_{r}$ is aligned with the (outward pointing ) normal vector ${\boldsymbol{n}}$ and the tangential vectors ${\boldsymbol{t}}$ and ${\boldsymbol{s}}$ are aligned with ${{\boldsymbol{e}}}_{\theta }$ and ${{\boldsymbol{e}}}_{\phi }$, respectively (with the vector directions for ϕ = 0 in the limit θ = 0). We consider a small box enclosing ${\boldsymbol{R}}$ where the outer and inner surfaces A_out and A_in have a constant distance of δ to the interface, and the lateral surface A_lat is at a constant angle θ₀ with respect to the symmetry axis. The geometry is shown in figure B1.

Zoom In Zoom Out Reset image size

Now let us consider the balance of the stress tensor equation (B.7) across the box, taking into account the curved geometry. The stress balance ${{\rm{\partial }}}_{\beta }{\sigma }_{\alpha \beta }$ can be written as

$$\begin{eqnarray}&&0=\oint {\rm{d}}A\,{\widetilde{n}}_{\beta }{\sigma }_{\alpha \beta },\end{eqnarray} \tag{ B.22 }$$

where α and β are Cartesian coordinates and $\tilde{n}$ the (local, outward pointing) normal vector of the box-surface. We can split this in three terms,

$$\begin{eqnarray}&&0=\int {{\rm{d}}A}_{{\rm{out}}}{\sigma }_{\alpha n}-\int {{\rm{d}}A}_{{\rm{in}}}{\sigma }_{\alpha n}+\int {{\rm{d}}A}_{{\rm{lat}}}{\sigma }_{\alpha t},\end{eqnarray} \tag{ B.23 }$$

where we used that the orientation of the normal vectors of the box coincides with the normal/tangential vector of the interface.

On the inner and outer areas A_in and A_out, the stress tensor presented in equation (B.8) with equilibrium stress tensor in equation (B.9) reduces to the form of the effective droplet model given after equation (6) in the main text, as the gradient terms are negligible for $\delta \gg w$. We now consider the limit of a sharp interface w → 0 with finite surface tension γ, and consider the case of a small box of thickness δ, which remains larger than the interfacial width. The components α = x, y of equation (B.23) vanish by symmetry. For α = z we find

$$\begin{eqnarray}&&0=\pi {\widetilde{R}}^{2}{\sin }^{2}{\theta }_{0}\,{\sigma }_{{nn}}^{+}-\pi {\widetilde{R}}^{2}{\sin }^{2}{\theta }_{0}\,{\sigma }_{{nn}}^{-}-2\pi \widetilde{R}\,{\sin }^{2}{\theta }_{0}\,\gamma ,\end{eqnarray} \tag{ B.24 }$$

where ${\sigma }_{{nn}}^{\pm }$ are the stress tensor components of the effective model, equation (4), inside and outside the interface at ${\boldsymbol{R}}$. Integration over the lateral box surface A_lat yields the last term, $\int {{\rm{d}}A}_{{\rm{lat}}}{\sigma }_{\alpha t}^{}\cong 2\pi \widetilde{R}\,{\sin }^{2}{\theta }_{0}\,\gamma$. We thus find that the mechanical equilibrium of a curved interface introduces a Laplace pressure 2γH,

$$\begin{eqnarray}&&0={\sigma }_{{nn}}^{+}-{\sigma }_{{nn}}^{-}-2\gamma H.\end{eqnarray} \tag{ B.25 }$$

We therefore recover the stress jump conditions of the effective droplet model, equation (6). Additionally, (B.25) together with (B.19) implies that the partial pressure needed to satisfy incompressibility is continuous across the interface, ${P}_{0}^{+}={P}_{0}^{-}$.

B.3.3. Dynamics of the effective droplet model

We now consider the dynamics of a non-equilibrium system with a droplet. We show how the continuum model is related to the bulk equations and jump conditions of the effective droplet model. For this we consider a droplet with a interface that is thin compared to the dynamical length scales l_±, so that we can describe the interface by local equilibrium. In the bulk phases we focus on the case where deviations from the equilibrium concentrations are small.

In the bulk phases, we expand the chemical potential equation (B.3) around the reference concentrations ${c}_{\pm }^{(0)}$. The gradient term −κ∇²c in the chemical potential is important within the interface, but can be ignored in the bulk phases, where the length-scales on which the concentration field varies are much larger than the interfacial width. Thus we can describe the chemical potential by

$$\begin{eqnarray}&&{\bar{\mu }}_{\pm }(c)\approx {\left.\displaystyle \frac{{\rm{d}}{\bar{\mu }}_{0}}{{\rm{d}}{c}}\right|}_{{c}_{\pm }^{(0)}}(c-{c}_{\pm }^{(0)}),\end{eqnarray} \tag{ B.26 }$$

which is ${\bar{\mu }}_{\pm }(c)\approx b(c-{c}_{\pm }^{(0)})$ for our specific free energy. With this simplification, equations (B.4) and (B.5) become the reaction-diffusion-convection equations (1) and (2) with diffusion constants ${D}_{\pm }=M{({\rm{d}}{\bar{\mu }}_{0}/{\rm{d}}{c})| }_{{c}_{\pm }^{(0)}}$ or D_± = Mb. Similarly we linearize the chemical reaction rate equation (B.6) in both phases. As we already chose a linear rate for the continuum model, we only need to relate the parameters k and ν with the constants k_± and ν_± of the effective model, with k_± = k, ν₊ = ν and ν₋ = kΔc − ν . Inserting the linearized chemical potential equation (B.26) into the equilibrium stress tensor (B.9) we find that momentum conservation in the bulk phases is given by the Stokes equation (4) with viscosities η_± = η, where the pressure p is determined by the incompressibility condition ∂_αv_α = 0.

We consider the droplet interface to be in local equilibrium. We therefore obtain equation (8) for the jump of the concentration field in the effective model. The incompressibility condition ∂_αv_α = 0 implies ${v}_{n}^{-}(R)={v}_{n}^{+}(R)$ at a sharp interface, and we consider an interface without slip length, so that ${{\boldsymbol{v}}}^{-}(R)={{\boldsymbol{v}}}^{+}(R)$. We thus find equation (7) of the effective model. The normal stress balance in equation (6) is derived in B.3.2.

As a last point we need to find equation (10) for the interface movement. We consider the concentration change in a box of width δ around the interface, see figure B1. We consider a box enclosing a point ${\boldsymbol{R}}$ on the interface at the time t aligned with the normal and tangential directions of the interface at ${\boldsymbol{R}}$. The interface may move with normal movement ${\partial }_{t}\hat{R}(t)$, with $\hat{R}(t)={\boldsymbol{R}}(t)\cdot {\boldsymbol{n}}$ and normal vector ${\boldsymbol{n}}$, while the box stays at a fixed position. The total change of material in the volume is given by

$$\begin{eqnarray}&&{\partial }_{t}{\int }_{V}{\rm{d}}{V}\,c=-{\int }_{A}{\rm{d}}{A}\,\tilde{{\boldsymbol{n}}}\cdot {\boldsymbol{j}}+{\int }_{V}{\rm{d}}{V}\,s(c),\end{eqnarray} \tag{ B.27 }$$

where V denotes the volume and A the area of the box. For small w and finite δ the concentration field c makes a jump from the surface A_in to A_out given by conditions (8) and (9) at $\hat{R}$. Within each phase, we can express the field by the boundary values at the interface equation (B.21) and a linear expansion,

$$\begin{eqnarray}c({\boldsymbol{r}},t)\simeq \left\{\begin{array}{ll}{c}_{-}({\boldsymbol{R}}(t))+{\rm{\nabla }}{c}_{-}({\boldsymbol{r}},t)\cdot ({\boldsymbol{r}}-{\boldsymbol{R}}(t)) & \text{inside droplet}\\ {c}_{+}({\boldsymbol{R}}(t))+{\rm{\nabla }}{c}_{+}({\boldsymbol{r}},t)\cdot ({\boldsymbol{r}}-{\boldsymbol{R}}(t)) & \text{outside droplet.}\end{array}\right.\end{eqnarray} \tag{ B.28 }$$

The chemical reaction is given in both phases by equation (B.6). For small δ and θ₀, we find for the left-hand side of equation (B.27) that δc vanishes to lowest order and

$$\begin{eqnarray}&&{\partial }_{t}{\int }_{V}{\rm{d}}{V}\,c={A}_{R}({c}_{-}({\boldsymbol{R}}(t))-{c}_{+}({\boldsymbol{R}}(t))){\partial }_{t}\hat{R}+O(\epsilon )+O({\theta }_{0}),\end{eqnarray} \tag{ B.29 }$$

where A_R is the area of the droplet interface enclosed by the box. For a spherical cap, ${A}_{R}=2\pi (1-\cos {\theta }_{0}){\hat{R}}^{2}$. We further find that the source term due to the chemical reaction scales with the volume of the box, and thus vanishes for a small box, ${\int }_{V}{\rm{d}}{V}\,s(c)=0+O(\epsilon )+O({\theta }_{0})$. The flux across the box can be expressed as

$$\begin{eqnarray}&&-{\int }_{A}{\rm{d}}{A}\,\tilde{{\boldsymbol{n}}}\cdot {\boldsymbol{j}}={A}_{R}{\boldsymbol{n}}\cdot ({{\boldsymbol{j}}}_{-}({\boldsymbol{R}}(t))-{{\boldsymbol{j}}}_{+}({\boldsymbol{R}}(t)))+O(\epsilon )+O({\theta }_{0}),\end{eqnarray} \tag{ B.30 }$$

where ${{\boldsymbol{j}}}_{\pm }({\boldsymbol{R}}(t))$ denotes the flux at ${\boldsymbol{R}}$ inside/outside the droplet. We thus find the normal movement of the interface,

$$\begin{eqnarray}&&{\partial }_{t}\hat{R}={\boldsymbol{n}}\cdot \displaystyle \frac{{{\boldsymbol{j}}}_{-}({\boldsymbol{R}}(t))-{{\boldsymbol{j}}}_{+}({\boldsymbol{R}}(t))}{{c}_{-}({\boldsymbol{R}}(t))-{c}_{+}({\boldsymbol{R}}(t))}.\end{eqnarray} \tag{ B.31 }$$

In the main text we use spherical coordinates centered at the droplet center. For a spherical droplet, the normal and radial movement would thus be the same. For a deformed droplet, we need to consider the relation between the normal interface movement, $\hat{R}(t)={\boldsymbol{R}}(t)\cdot {\boldsymbol{n}}$ and the radial movement $R(t)={\boldsymbol{R}}(t)\cdot {{\boldsymbol{e}}}_{r}$. At fixed angles θ and ϕ, the interface movement is given by ${\partial }_{t}{\boldsymbol{R}}={\partial }_{t}R\,{{\boldsymbol{e}}}_{r}$. Using ${\partial }_{t}\hat{R}={\partial }_{t}{\boldsymbol{R}}(t)\cdot {\boldsymbol{n}}$, we find a relation between the radial and normal movement, ${\partial }_{t}R={\partial }_{t}\hat{R}/({\boldsymbol{n}}\cdot {{\boldsymbol{e}}}_{r})$. This relation, together with equation (B.31), yields the interfacial movement equation (10) presented in the main text.

We thus recover all dynamical equations of the effective droplet model from the continuum model based on irreversible thermodynamics. Note that the specific choice of the free energy leads to specific relations between parameters of the effective model such as D₊ = D₋. Our derivation shows the relation between both models in the case where the interface width w is small compared to the droplet size, $R/w\gg 1$, and the chemical diffusion length, ${l}_{\pm }/w\gg 1$. Additionally, we focused on the case where the concentrations in the phases are similar to the concentrations in equilibrium and have small concentration gradients. These conditions are not valid in all systems. Most importantly, the chemical reactions can drive concentrations far away from the equilibrium phase concentrations ${c}_{\pm }^{(0)}$. The resulting behaviors, such as the formation of new interfaces associated with instabilities of the spinodal decomposition regime, are not captured in the effective droplet model.

B.4. Comparison of the droplet dynamics in the continuum model and the effective model

Here we compare the analytical predictions of the effective model for the instability with numerical calculations of the continuous model for different values of the renormalized viscosity F. For this we numerically solved the dynamic equations of the continuous model starting with a droplet with a small initial deformation of mode l = 2. We fitted the dynamical behavior of the mode to an exponential function, with yields a numerical estimate for the eigenvalue μ₂. In figure B2 the resulting eigenvalues are shown, together with the eigenvalue of corresponding parameters of the effective model. We find that the value of F for which droplet shapes become unstable is very similar to the value predicted by the effective model. The eigenvalues are qualitatively similar to the ones of the effective model, despite working in an a parameter regime where the interfacial width and the differences of concentration within a phase cannot be considered very small, so that the models are not necessarily comparable.

Zoom In Zoom Out Reset image size

To generate the data in the figure, we initialized droplets with a small shape perturbation for different values of F. All parameters and initial conditions were chosen as described in appendix B.2. We found that for $F\geqslant 100$ droplets divide, while they are stable for $F\leqslant 1$. For F = 10, the shape deformation was very slow, so that division was not seen in the time interval T/τ = 4000. For 10 < F < 100, as well as $F=\infty$, we fitted radius and spherical harmonic deformation to the concentration field using equation (B.15). For short times, the droplet radius changes as the concentration field and droplet size go towards the stationary values. After that, the shape deformation grows until the droplet deforms so strongly that the fitting fails. By hand we chose intermediate time windows for the simulations where the size was stationary and the shape deformation small. In these windows we fitted the deformation amplitude (compare equation (B.15)) with an exponential function, ${{A}{\rm{e}}}^{{\mu }_{2}t}+B$ with parameters A, B and eigenvalue μ₂ to the l = 2 mode of the shape deformation.

Appendix C.: Estimation of parameters