Role of hydrodynamic flows in chemically driven droplet division

We study the hydrodynamics and shape changes of chemically active droplets. In non-spherical droplets, surface tension generates hydrodynamic flows that drive liquid droplets into a spherical shape. Here we show that spherical droplets that are maintained away from thermodynamic equilibrium by chemical reactions may not remain spherical but can undergo a shape instability which can lead to spontaneous droplet division. In this case chemical activity acts against surface tension and tension-induced hydrodynamic flows. By combining low Reynolds-number hydrodynamics with phase separation dynamics and chemical reaction kinetics we determine stability diagrams of spherical droplets as a function of dimensionless viscosity and reaction parameters. We determine concentration and flow fields inside and outside the droplets during shape changes and division. Our work shows that hydrodynamic flows tends to stabilize spherical shapes but that droplet division occurs for sufficiently strong chemical driving, sufficiently large droplet viscosity or sufficiently small surface tension. Active droplets could provide simple models for prebiotic protocells that are able to proliferate. Our work captures the key hydrodynamics of droplet division that could be observable in chemically active colloidal droplets.

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 nonspherical droplets drives hydrodynamic flows that redis-
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 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 Fig. 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), Here, the indices + and − refer to quantities outside and inside the droplet, respectively. The flux j consists of advection by the fluid velocity v and a diffusion flux, where D ± denotes the diffusion constant of the droplet material in the two phases. The chemical reaction is described by the concentration-dependent rate s ± (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 (0) ± in each phase: where c (0) ± are the equilibrium concentrations that coexist at equilibrium across a flat interface. We have defined the reaction rate ν ± = s(c (0) ± ) and the reaction constants k ± = ds(c (0) ± )/dc. We focus on the case of positive coefficients k ± > 0 and ν ± > 0, where B is produced outside the droplet, and degraded inside, see Fig. 1.
The hydrodynamic flow velocity v obeys Stokes equation of an incompressible fluid, which accounts for momentum conservation ∂ α σ αβ = 0, where the stress tensor is given by 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 ∇ · v = 0. The bulk equations (1)(2)(3)(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 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 σ ± αβ (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 σ ± nn = n α σ ± αβ n β and σ ± nt = n α σ ± αβ t β , where n α is a unit vector normal to the surface and t α is a unit vector tangent to the surface. Eq (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, The concentration field c is discontinuous across the interface with values given by 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 where n is a unit vector normal to the surface and e r is a unit vector in radial direction. Eq. (10) captures both convection of the interface by flows and droplet growth and shrinkage by addition or removal of material. We find nonequilibrium steady state solutions to Equations (1-10) with a spherical droplet of stationary radius R and stationary concentration fieldc(r), where r is the radial coordinate, see Appendix A. The stationary pressurep exhibits a jump 2γ/R across the interface and no hydrodynamic flows exist,v = 0. An example for a stable non-equilibrium steady state with steady state concentration profile inside and outside the droplet of radius R is shown in Fig. 1.
We first discuss the properties of these stationary states as a function of external supersaturation = ν + /(k + ∆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, ∆c = c (0) − − c (0) + and we have introduced the time scale τ = w 2 /D + , where w = 6β + γ/∆c is a characteristic length scale. The stationary radii as a function of supersaturation are shown In Fig. 2A-C as solid lines for different values of A. For values of smaller than a threshold value 0 , no stationary radius exists. For values > 0 two steady state radiiR c and R s exist, which become equal at 0 where they approach a valueR 0 . The smaller steady state radiusR c is a critical nucleation radius similar to the critical droplet radii found in passive systems. The larger radius denotedR s stems from the interplay of phase separation and chemical reactions 26,39 . As the supersaturation reaches a value  In panel B the scaling behavior of the nucleation radiusRc and the stationary radiusRs are indicated. D-F: Stability diagram of stationary droplets of sizeRs, as a function of reaction amplitude A and supersaturation for different dimensionless viscosities F = ∞, 1000, 10 (left to right). For small supersaturation and large reaction amplitudes, no stationary radius exists (white). For large supersaturation, the stationary radius diverges (gray). In the region between these regimes, the stationary solution can be stable (blue) or unstable (red) with respect to shape perturbations of the l = 2 mode. For decreasing F , the stable regime grows, and the minimal supersaturation * at which an instability can be found increases, as well as the corresponding reaction amplitude A * . The scaling relations (dashed lines) for the regime of stable droplets and the onset of instability are indicated, with prefactors according to A 5. (Parameters: η+/η− = 1, k+/k− = 1, ν−/(k−∆c) = 1, 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 0 vanishes as 0 ∝ A 1/3 . The critical nucleation radius then behaves asR c w/(6 ) and the larger steady state radiusR s w(3 A) 1/2 where 0 ∞ , see Fig. 2B and 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 hy-drodynamic flows on chemically driven shape instabilities that can give rise to droplet division. We perform a linear stability analysis at the stationary state given byX = (c,R,p,v) for small perturbations δX = (δc, δR, δp, δv). The dynamics of these perturbations can be represented using eigenmodes δX = n,l,m nlm X nlm e µ nlm t , with X nlm = (c nl Y lm ,RY lm , p l Y lm , v lm ), where Y lm (θ, φ) 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 v lm (r, θ, ϕ) can be expressed as where Y lm (θ, ϕ) = e r Y lm , Ψ lm (θ, ϕ) = r∇Y lm and Φ lm (θ, ϕ) = e r × Ψ lm are vector spherical harmonics 45 and the radial functions v r lm (r), v lm (r) and v (2) lm (r) characterize the velocity field. The radial functions can be obtained by solving the linearized dynamic equations using the corresponding boundary conditions, see 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 Eq. (10), we obtain an equation for the eigenvalue µ nlm , Here, the primes denote radial derivatives. Note that Eq. (13) is an implicit equation for the eigenvalues µ nlm because the radial concentration modes c nl (r) depend on µ nlm , see A. Eq. (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 R =R s and unstable forR =R c . Thus droplet smaller thanR c will vanish and droplets larger will grow towards the sizeR s . Thus we consider the stability ofR =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 Fig. 2A-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 Fig. 2D-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 ≥ * , which depends on F . The value of A for which the shape instability occurs at = * is denoted A * , see Fig 2E.
For small A, the onset of instability can be describes by simple scaling behaviors. In the limit of small A and for ∞ , we find Fig. 2E-F). For A < A * , hydrodynamic flows govern the onset of instability which occurs at a value of A which behaves as A ∝ −1 F −2 . For A > A * , hydrodynamic flows can be neglected as compared to diffusion fluxes and the onset of instability occurs for A ∝ 3 . These two scaling regimes are indicated in in Fig. 2D-F by dashed lines. A derivation of these results including prefactors is given in 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 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. Fig. 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, Fig. 3 A. As the droplet elongates and its waistline shrinks, the flow field becomes more complex, see Fig. 3  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 Fig. 2. This raises the question under what experimental conditions active droplets would become unstable and division could be observed. Ignoring hydrodynamic flows, F → ∞, 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η − /(γτ ) k B T /(6πγwa), where we have used τ = w 2 /D and D k B T /(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 C. For soft colloidal liquids or p-granules, we estimate values between F ≈ 10 − 10 4 . We can discuss these values using the stability diagrams in Fig. 2D-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 4 , 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 −4 mM/s, see 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 fluxdriven 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.
Here, we discuss stationary solutions to equations (1)(2)(3)(4)(5)(6)(7)(8)(9)(10) in the main text with spherical symmetry and without hydrodynamic flowsv = 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,p The steady state concentration profiles in the presence of chemical reactions are given by 39 , where i 0 (x) = 2 sinh(x)/x and k 0 (x) = e −x /x denote modified spherical Bessel functions of order zero of the first and second kind, respectively. The characteristic length scales l ± = (D ± /k ± ) 1/2 are set by reaction rate constants and diffusion coefficients. The parameters A ± are determined by the boundary condition at the droplet interface, Eq. (8)(9) in the main text, Stationarity of the droplet radiusR implies see Eq. (10) in the main text. Note that this equation typically has zero, one or two solutions for a given set of parameters.

Linearized dynamics
We introduce small perturbations to the spherically symmetric stationary state, with p =p + δp, v = δv, c =c + δc and R =R + δR and write the dynamics of these perturbations to linear order. The linearized dynamics reads 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 radiusR, with perturbation of the curvature δH = H(R) − H(R).
The linearized dynamics can be decomposed in spherical harmonics, see Eq (11) in the main text. The curvature perturbation then takes the form with h l = (l 2 + l − 2)/2.

Hydrodynamic eigenmodes of the linearized dynamics
We can expand the hydrodynamic eigenmodes using a basis of vector spherical harmonics, see Eq. (12) in the main text. The velocity boundary conditions Eq. (7) in the main text for the mode amplitudes read The stress boundary conditions (see Eq. (5-6) in the main text) at the interface read We solve the radial profiles of the modes with a polynomial ansatz and exclude functions that diverge for r → 0 or r → ∞ inside and outside the droplet, respectively. The pressure is then given by For the hydrodynamic flow velocity we obtain Here, we have defined where ∆ = η + /η − denotes the ratio of the viscosities inside and outside the droplet.

Concentration eigenmodes
The equation for the radial part of the concentration eigenmode is The boundary conditions atR are The left-hand side of Eq. (A34) constitutes an inhomogeneity The solution c ± nl (r) of the inhomogeneous equation (A34) that satisfies the boundary condition Eq. (A36) can be constructed from a particular solution c ± nl,p (r) of the inhomogeneous equation to which solutions c ± nl,h (r) of the homogeneous equation with f ± l = 0 are added to satisfy the boundary conditions, Eq. (A36). This can be expressed as where the coefficients α ± read with a ± l = c ± nl (R). We are especially interested in the case of unstable modes with µ nl > 0. Therefore we focus on the solution of equation (A34) for λ ±2 nl > 0 and k ± > 0. In this case, the homogeneous equation with f ± l = 0 is a modified Helmholtz equation which is solved by modified spherical Bessel functions, c − nl,h (r) = i l (λ − nl r) and c + nl,h (r) = k l (λ + 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, with the radial part of the inhomogeneity f ± l (r) given by Eq. (A37). 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 atR can be expressed as with Using the equation for the shape perturbations (A10), and using Eqns (A43) and (A44), we obtain Eq. (13) in the main text. This equation determines the eigenvalue µ nlm of the hydrodynamic modes.

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. 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 curvesR( ) shown in Fig. 2A for different values of A. We can identify two regimes in the figure. The first is the region of small , ∼ 0 , which corresponds to the minimum of (R). The second is the region of ∞ where the stationary radius diverges. For A → 0, we see that 0 goes to zero while ∞ stays constant, and both are connected by a straight line that indicates scaling behavior ofR =R s . This increasing separation between 0 and ∞ (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, In this limit of small A, the characteristic length-scales of the concentration field become large with l ± ∝ A −1/2 . To find scaling regimes in equation ((A47)), we change variables in Eq. (A47) from (A,R) to (A,R) withR =RA a /w, where a is an exponent. For a = 1/3 we find the behavior of (R) close to 0 andR 0 , whereˆ = A −1/3 becomes independent of A for small A. This function describes the supersaturation as a function of radius around the threshold value 0 . Due to the inverted presentation (R) instead ofR( ) the function captures both the nucleation radiusR c and the larger radiusR s . The threshold value 0 can be obtained from Eq. (A48) by minimizingˆ for fixed A as ∂ˆ /∂R = 0. It behave as For large and smallR, Eq (A48) describes the steady radiiR s andR c , respectively, for which ≥ 0 . For large , the critical radius obeysR while the larger stationary radius isR In Fig. 2B, the scaling behaviors given by Eq. (A51) and Eq. (A50) are indicated by dashed lines. At = 0 both radii meet atR =R 0 , whereR For a = 1/2,R/l ± becomes independent of A and ForR/l ± 1, the stationary radius obeys Eq. (A51) and is thus the larger stationary radiusR s . ForR/l ± 1, we obtain the divergence ofR s as approaches ∞ with We now discuss scaling relations for the onset of instability in the (A,ˆ ) plane in the limit of small A, which give the trends shown as dashed lines in Fig. 2D-F. We use the scaling of the stationary radiusR =R s close to 0 witĥ R =RA 1/3 /w,ˆ = A −1/3 andl ± = l ± A 1/2 in Eq. (13) to obtain where f C1 and f C3 are defined in Eq. (A33) and with h l = (l 2 + l − 2)/2. For large mode index l, We now consider conditions for which µ nlm = 0 for small A and the mode (n, l, m) becomes unstable. Using (A48) in (A55), we find a relation betweenˆ andR at the onset of instability µ nlm = 0, This curve captures the scaling behavior of the onset of instability for different parameters in theR − plane, corresponding to the red dotted line in Fig. 2A-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 Fig. 2D-F. At this onset, both (A58) and (A48) need to be satisfied. We use both equations to eliminateR. We find a crossover regime with relations A * ∼ 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 In Fig. 2D-F, the dashed lines indicate these two scaling solutions in the limit A → 0 and F → ∞ 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 Eq. (13) for small , and larger deviations in the regime close to 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 (r, t) and c = c B (r, t) that depend on position 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 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 A c A + m B c B is constant. Therefore, we can describe the system by the concentration c(r, t) of the droplet material B only.
We use the following double-well free energy density 46 with ∆c = c The state of the system is characterized by the free energy 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μ = δF [c]/δc, governs demixing and can be split into local and nonlocal contributions,μ =μ 0 − κ∇ 2 c with The dynamics of the concentration field is described by 50,51 Here, M is a mobility coefficient of the droplet material and v is the hydrodynamic velocity. The source term s(c) describes chemical reactions, for which we choose for simplicity a linear concentration dependence, The reaction flux given in Eq. (B6) 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 v can be calculated using momentum conservation, with momentum ρv α and stress tensor σ αβ , where α and β number cartesian coordinates x, y, z. We can write the stress tensor σ αβ as where the first term describes advection of the stress tensor, σ eq αβ and σ d αβ denote the equilibrium and dissipative stress tensors. The equilibrium stress tensor is given by Here,μc − f is the osmotic pressure of the droplet material, and δ αβ denotes the Kronecker delta. Incompressibility is enforced by an additional partial pressure P 0 . The deviatory stress tensor can be found as thermodynamic force related to momentum by writing the entropy production rate, where η and η denote viscosities, and v αβ = (∂ α v β + ∂ β v α )/2 is the symmetric strain tensor.
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μ, Eq. (B3), implicitly, while evaluating the rest ofμ as well as the advection term of the fluxes, vc, 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 (∂ n c = 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 2 for a mode with wave vector 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), with vector field ψ and scalar field φ, and velocity parts v ψ = ∇ × ψ and v φ = −∇φ. With this, we find and thus, using incompressibility, ∇ · v = 0, we can calculate φ. We thus find the incompressible part of the velocity field 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 ∆c = c where the characteristic length scale is w = 2(κ/b) 1/2 . The relevant dimensionless model parameters are c (0) + /∆c, kt 0 , and ν − t 0 /∆c. We choose c (0) + /∆c = 0, kt 0 = 10 −2 , νt 0 = 2 · 10 −3 and ηŵ 3 /(t 0ê0 ) = 2. Additionally, we use as box-length L/ŵ = 100 in all 3 dimensions, number of grid-points in one direction N = 128 and simulation time T /t 0 = 4 · 10 3 . For the time step, we start with a timestep of ∆t/t 0 = 10 −4 , and double the timestep to a final step size of ∆t/t 0 = 0.01.
We start with initial conditions R = R 0 (1 + Y 2,0 ), with R 0 /ŵ = 7 and = 1. The concentration field at positions r is initialized by the function where d(r) is the oriented distance of r to the nearest point on the ellipsoid. The value of d(r) is negative for points inside the droplet and positive for points outside.
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.
a. 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 concentrationc with c where w = 2(κ/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 where H the mean curvature of the droplet and 2γ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 (0) ± plus a small perturbation, . For the free energy Eq. (B2), we find β ± = 2/(b∆c), which is related to the interfacial width as w = 6γβ + /∆c.

b. 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 radiusR = 1/H. We focus on the case where the interface is rotationally symmetric around the considered point R, and where the curvature does not change along the interface. We use spherical coordinates, where the radial vector e r is aligned with the (outward pointing ) normal vector n and the tangential vectors t and s are aligned with e θ and e φ , respectively (with the vector directions Geometry for the force balance. We consider a spherical cap of the droplet interface, with a box with constant distance δ to the interface inside and outside. The normal and tangential vectors n, t and s of the interface are shown, as well as the normal vectorñ of the box. The origin of the spherical coordinate system is the center of the sphere that describes the interfacial curvature, with radiusR, while θ0 gives the polar angle of the cap. for φ = 0 in the limit θ = 0). We consider a small box enclosing 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 θ 0 with respect to the symmetry axis. The geometry is shown in Fig. 4. Now let us consider the balance of the stress tensor Eq. (B7) across the box, taking into account the curved geometry. The stress balance ∂ β σ αβ can be written as where α and β are cartesian coordinates andñ the (local, outward pointing) normal vector of the box-surface. We can split this in three terms, 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 Eq. (B8) with equilibrium stress tensor in Eq. (B9) reduces to the form of the effective droplet model given after Eq. (6) in the main text, as the gradient terms are negligible for δ 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 Eq. (B23) vanish by symmetry. For α = z we find where σ ± nn are the stress tensor components of the effective model, Eq. (4), inside and outside the interface at R. Integration over the lateral box surface A lat yields the last term, dA lat σ αt ∼ = 2πR sin 2 θ 0 γ. We thus find that the mechanical equilibrium of a curved interface introduces a Laplace pressure 2γH, We therefore recover the stress jump conditions of the effective droplet model, Eq. (6). Additionally, (B25) together with (B19) implies that the partial pressure needed to satisfy incompressibility is continuous across the interface, c. 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 Eq. (B3) around the reference concentrations c (0) ± . The gradient term −κ∇ 2 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 bȳ Similarly we linearize the chemical reaction rate Eq. (B6) 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 Eq. (B26) into the equilibrium stress tensor (B9) 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 Eq. (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 v − (R) = v + (R). We thus find Eq. (7) of the effective model. The normal stress balance in Eq. (6) is derived in B 3 b.
As a last point we need to find Eq. (10) for the interface movement. We consider the concentration change in a box of width δ around the interface, see Fig. 4. We consider a box enclosing a point R on the interface at the time t aligned with the normal and tangential directions of the interface at R. The interface may move with normal movement ∂ tR (t), withR(t) = R(t) · n and normal vector n, while the box stays at a fixed position. The total change of material in the volume is given by 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) atR. Within each phase, we can express the field by the boundary values at the interface Eq. (B21) and a linear expansion, c(r, t) c − (R(t)) + ∇c − (r, t) · (r − R(t)) inside droplet c + (R(t)) + ∇c + (r, t) · (r − R(t)) outside droplet (B28) The chemical reaction is given in both phases by Eq. (B6). For small δ and θ 0 , we find for the left-hand side of Eq. (B27) that δc vanishes to lowest order and where A R is the area of the droplet interface enclosed by the box. For a spherical cap, A R = 2π(1 − cos θ 0 )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, V dV s(c) = 0 + O( ) + O(θ 0 ). The flux across the box can be expressed as − A dAñ · j = A R n · (j − (R(t)) − j + (R(t))) + O( ) + O(θ 0 ) (B30) where j ± (R(t)) denotes the flux at R inside/outside the droplet. We thus find the normal movement of the interface, ∂ tR = n · j − (R(t)) − j + (R(t)) c − (R(t)) − c + (R(t)) .
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,R(t) = R(t) · n and the radial movement R(t) = R(t) · e r . At fixed angles θ and φ, the interface movement is given by ∂ t R = ∂ t R e r . Using ∂ tR = ∂ t R(t) · n, we find a relation between the radial and normal movement, ∂ t R = ∂ tR /(n · e r ). This relation, together with Eq. (B31), yields the interfacial movement Eq. (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 1, and the chemical diffusion length, l ± /w 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 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 µ 2 . In figure 5 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.
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 B 2. We found that for F ≥ 100 droplets divide, while they are stable for F ≤ 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 = ∞, we fitted radius and spherical harmonic deformation to the concentration field using Eq. (B15). 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 Eq. (B15)) with an exponential function, Ae µ2t + B with parameters A, B and eigenvalue µ 2 to the l = 2 mode of the shape deformation.