Nucleation of chemically active droplets

Driven chemical reactions can control the macroscopic properties of droplets, like their size. Such active droplets are critical in structuring the interior of biological cells. Cells also need to control where and when droplets appear, so they need to control droplet nucleation. Our numerical simulations demonstrate that reactions generally suppress nucleation if they stabilize the homogeneous state. An equilibrium surrogate model reveals that reactions increase the effective energy barrier of nucleation, enabling quantitative predictions of the increased nucleation times. Moreover, the surrogate model allows us to construct a phase diagram, which summarizes how reactions affect the stability of the homogeneous phase and the droplet state. This simple picture provides accurate predictions of how driven reactions delay nucleation, which is relevant for understanding droplets in biological cells and chemical engineering.

Nucleation is a stochastic process that relies on thermal fluctuations to create a sufficiently large nucleus that can grow spontaneously [17][18][19].This is because creating the droplet interface costs energy, so tiny droplets generally dissolve.Classical nucleation theory predicts that the typical nucleation time scales exponentially with the energy barrier associated with the critical nucleus.While this theory is well-understood for passive systems, it is unclear how it can be extended to active systems, where free energies are generally unavailable.To overcome this challenge, we use an equilibrium surrogate model to reveal how driven reactions controlling droplet size suppress nucleation substantially.

II. MODEL
We study an isothermal fluid comprised of precursor material A that can convert into droplet material B by chemical reactions.For simplicity, we consider an incompressible fluid where both species have equal molecular volume ν, so the state of the system is characterized by the concentration c(r, t) of species B, while the concentration of A is ν −1 − c(r, t).The dynamics are governed by the continuity equation where j denotes the diffusive exchange flux and the source term s describes chemical transitions.The passive diffusive flux j is driven by the gradient of the chemical potential, j = −Λ d ∇µ + η, where Λ d is the diffusive mobility and η is the diffusive thermal noise, which obeys η i (r, t) = 0 and the fluctuation dissipation theorem η i (r, t)η j (r , t [20].The exchange chemical potential, µ = δF [c]/δc, follows from the free energy functional where f (c) is the local free energy density and κ penalizes compositional gradients.For simplicity, we consider where a 1 , a 2 , and a 4 > 0 are phenomenological coefficients.Without reactions (s = 0), Eqs. ( 1)-( 3) describe passive phase separation with a critical point at for a 2 = 0.For a 2 > 0, the spinodal line is given by c sp = 1 2ν ± a 2 /(3a 4 ) and the binodal is defined by coexisting equilibrium concentrations c out = 1 2ν − a 2 /a 4 and c in = 1 2ν + a 2 /a 4 in dilute and dense phases, respectively.
The system becomes active when phase separation is augmented by driven chemical reactions.We first consider a reaction flux s comprising passive conversion of A and B as well as an active conversion involving chemical energy ∆µ provided by a fuel [14], where Λ p r and Λ a r determine the rates of the respective reactions and η r models thermal fluctuations.Motivated by enzymes that co-localize with the droplet, we scale the rate of the active reaction with the concentration c of the droplet material.This choice allows stationary states where droplet material B turns into precursor A inside droplets, while B is replenished outside, thus controlling droplet size [14].The reactive thermal noise η r obeys η r (r, t) = 0 and η r (r, t)η r (r , t Active droplets are only stable if reactive fluxes are weak compared to diffusive fluxes [14].Consequently, the reactive noise η r is much weaker than the diffusive noise η and we neglect it in the following.Fig. 1a shows that the reaction flux s given in Eq. ( 4) is a non-monotonous function of the composition c.In particular, there are two stable homogeneous stationary state, which correspond to (meta-)stable dilute and dense systems.The main question in this paper concerns how active droplets nucleate from the dilute homogeneous state c(r) = c 0 .

A. Chemical reactions hinder nucleation
To investigate nucleation, we first perform numerical simulations of Eqs. ( 1)-( 4) in two-dimensional system with periodic boundary conditions [21]; see Fig. 1b.Repeating the simulations many times, we observe that the first droplet nucleates at random times t nucl ; see Fig. 1c.Assuming an exponential distribution of t nucl , we define the nucleation time τ as the ensemble average of t nucl .Fig. 1d shows that τ increases for stronger interactions (larger a 2 ), as expected for nucleation of passive droplets [22].More importantly, larger reaction rates Λ a r lead to longer nucleation times τ , indicating that active chemical reactions hinder nucleation.This result can be understood intuitively since the reactions stabilize the homogeneous state; see Fig. 1a.They thus help to dissolve a small accumulation of droplet material B, reducing the probability that a critical nucleus forms.

B. Surrogate equilibrium system reveals additional energy barrier
To quantitatively understand the effect of driven reactions on nucleation, we next map our system to an approximate equilibrium system.To do this, we linearize the reaction flux s around the dilute homogeneous stationary state c 0 , where k > 0 for a stable state; see Fig. 1a.Fig. 1d shows that the linearized reactions influence the nucleation time τ similarly to the full reaction flux s.The linearization allows us to map the dynamics given by Eq. ( 1) to a passive system, ∂ t c ≈ Λ d ∇ 2 δ F [c]/δc, with the augmented free energy functional where Ψ obeys the Poisson equation ∇ 2 Ψ = c 0 − c(r) and thus mediates long-ranged, Coulomb-like interactions [23][24][25].
We use the equilibrium surrogate model to investigate the energy landscape of nucleation.In particular, we use Eq. ( 6) to map the minimal energy path connecting the metastable homogeneous state with the equilibrium state containing one droplet using a proxy for droplet size as a reaction coordinate x; see Appendix I.For each value of x, we use a constrained optimization to determine the spherically symmetric composition c(r) that minimizes the energy F given by Eq. ( 6).Fig. 2a shows that the resulting profiles feature an increasing density peak, analogous to passive systems [19].However, the nucleus is also surrounded by a depletion zone originating from chemical reactions.The sequence of profiles defines the minimal energy path, from which we obtain the energy barrier ∆F as the difference between the maximal energy and the energy F (x = 0) of the homogeneous state; see Fig. 2b.Fig. 2c shows that the energy barrier ∆F depends on the reaction rate k and this dependence is approximately linear (Fig. 2d), suggesting that the long-ranged term in Eq. ( 6) could explain the suppressed nucleation caused by reactions.
We hypothesize that the increasing energy barriers explain how larger reaction rates k lead to longer nucleation times τ (see Fig. 1d).Nucleation theory predicts that τ increases exponentially with the energy barrier ∆F [22], where A is a kinetic prefactor.Fig. 1d shows that this relation explains the numerical data, particularly for larger ∆F at higher k and a 2 .The deviation at smaller k are expected since the assumptions leading to Eq. ( 7) break down for smaller ∆F [19].We conclude that the energy barriers derived from the equilibrium surrogate model explain how nucleation times increase with reaction rates.

C. Classical nucleation theory leads to phase diagram extended by chemical reactions
Motivated by the success of nucleation theory, we next approximate the minimal free energy path using the radius R of a droplet as a reaction coordinate.Assuming that the droplet with homogoneous concentration c in is embedded in an infinite system of concentration c 0 , the free energy F can be separated into contributions of bulk phases, interface, and chemical reactions, where V = πR 2 and A = 2πR in two dimensions.Classical nucleation theory implies the free energy difference g = f (c 0 )−f (c in )+µ(c 0 )(c in −c 0 ) between the phases and surface tension γ = 8κc 3 2 /(3c 4 ) [13], which is a good approximation for c 0 ≈ c out .We show in Appendix II that the free energy associated to reactions is approximately  where we neglected terms proportional to k(c out − c 0 ) 2 .Without reactions (k = 0), F (R) given by Eq. ( 8) has a single maximum at the critical radius R pas crit = γ/g with a corresponding energy barrier ∆F = πγ 2 /g; see Fig. 3a.Once nuclei exceed this critical size (by nucleation), they grow indefinitely.In contrast, Eq. ( 8) predicts that reactions (k > 0) increase the free energy of large droplets, implying a minimum at finite radius corresponding to stable droplets [10].Concomitantly, the energy barrier ∆F is elevated, consistent with increased nucleation times.Approximating the barrier by F (R pas crit ), we find which explains the linear dependence of ∆F on k observed in Fig. 2d.Taken together, this simplified picture demonstrates how large rates k gradually disfavor the droplet state until only the homogeneous state remains stable at k > k max with Finally, we use the simplified free energy of the equilibrium surrogate model to study the influence of the concentration c 0 of the homogeneous state.In particular, we determine the minimal value of c 0 beyond which the droplet state can be stable as a function of the interaction parameter a 2 .In passive systems (k = 0), the resulting line corresponds to the binodal curve.In active system (k > 0), this line is shifted to larger concentrations, thus enlarging the region where the homogeneous system is stable.The homogeneous system becomes unstable at the spinodal line, which can be determined from a linear stability analysis of Eq. ( 1); see Appendix III.Fig. 3b shows that these predictions based on Eq. ( 8) agree with numerical simulations probing the stability of the homogeneous and droplet state.Both predictions illustrate how driven chemical reactions destabilize the droplet state.

IV. DISCUSSION
In summary, we illuminated how driven chemical reactions affect the phase diagram of phase separating systems.To do this, we exploited the equilibrium surrogate of the active system to show how reactions favor the homogeneous state relative to the droplet state, which explains the suppressed nucleation qualitatively.Similar behavior was found for equilibrium system with true long-ranged interactions [25].Although the modified phase diagram was derived from the surrogate model, it is not a thermodynamic phase diagram of the phase separating system with driven reactions.For instance, the compositions of the coexisting phases at the interface are still governed by the binodal and tie lines of the passive phase diagram [13].The energy barrier associated with reactions depends linearly on their rate k, likely because reactions are weak and the system is dominated by phase separation.This implies that k decreases nucleation rates exponentially.We showed that this dependence persists for thermodynamically consistent reactions and expect that it is a general feature of phase separating systems with reactions that have a stable dilute phase.Moreover, since our derivation of the influence of the reactions is independent of the details of the free energy density, we expect that reactions suppress nucleation in a wide range of phase separating systems.
We presented results for the simple case of a binary fluid in two dimensions.While we expect that active reactions also slow nucleation in more complicated situations, it will be vital to extend our theory to three dimensions (e.g., to capture spontaneous divisions [11]) and many components (allowing for additional stable stationary states [14,15]).For better quantitative agreement, it might also be necessary to improve our treatment of nucleation theory, e.g., by describing how reactions affect the curvature of the surrogate free energy, which affects nucleation rates via the Zeldovich factor [26].However, the ultimate test of our theory will come from experiments, either from existing active droplets in biological cells [27] or in promising synthetic systems [28,29].Experiments in cells also suggest that more complex behaviors are possible, including periodic nucleation [30] and multi-step nucleation for fiber formation [31], which might involve secondary nucleation [32].In these situations, heterogeneous nucleation is likely relevant [33,34] and there are examples where nucleation is controlled by catalytically-active nucleation sites [9,35].Taken together, our approach of an equilibrium surrogate model will likely prove vital for studying nucleation in these more challenging situations.
Integrating py parts and dropping the boundary term, we obtain For the radial-symmetric system in two dimensions that we consider, we find Using Eq. (B6), we obtain where d 1 is an integration constant.Denoting the two branches by Ω in and Ω out , we find from Eq. B10 where L denotes the radius of the spherical system.Since the integrals must converge in the limit L → ∞, we conclude d 1 = 0. Hence, the second term of Eq. (B13) reduces to in the limit of large systems, L → ∞.For large droplet radii (R → ∞), this contribution scales as R and can thus be neglected compared to the contribution from inside the droplet, which scales as R 4 .Additionally, the prefactors are different, (c 0 − c in ) 2 (c 0 − c out ) 2 .Since critical radii become large close to the binodal concentration (where classical nucleation is valid), we can neglect the contribution from the outside to find Eq.9.
Appendix C: Spinodal line from linear stability analysis We perturb the dynamic equation given by Eq. 3 around the homogeneous state using the ansatz c(r, t) = c 0 + e ωt+iqr . (C1) To linear order in , we find ω = Λ d a 2 q 2 − 3Λ d a 4 q 2 (2c 0 ν − 1)

(C3)
The spinodal is given by the concentration c 0 where the ω max becomes zero.We thus solve Eq.C3 for c 0 and find )

FIG. 1 .
FIG. 1.Chemical reactions increase nucleation times.(a) Reaction flux s as a function of the concentration c of a homogeneous system for the full (solid blue line, Eq. (4)) and linearized model (dashed orange line).The spinodal concentration csp of the passive system (dotted green line), the two stable fixed points (filled disks), and the unstable fixed point (open circle) are marked.(b) Snapshots of the concentration field c of droplet material obtained from stochastic numerical simulation in two dimensions.The time between snapshots is 10/k0 and the interaction strength is a2 = 200 νkBT .(c) Distribution of measured nucleation times t nucl in the linearized model for various reaction rates k for a2 = 150 νkBT .Black lines show exponential distributions of equivalent mean τ = t nucl .(d) Nucleation time τ as a function of k for the full model (disks, k = −s (c0) ∝ Λ a r ) and the linearized reactions (triangles) for various interaction strengths a2/(νkBT ).Solid lines show predictions of Eq. (7) with A as a single fit parameter for all curves.(a-d) Additional parameters are a1ν = −1.34a2, a4 = 4 a2ν, Λ p r /Λ a r = 0.0311, ∆µν = 1.46 a2, ν = w 2 , w = 2κ/a2, and k0 = Λ d a2w −2 .
t e x i t s h a 1 _ b a s e 6 4 = " K X n W 4 c n e I 1 K y 6 J T z e 1 U 2 Z j 8 d w z r 5 J 2 v e Z d 1 i 7 u 6 p U G y e M o w g m c Q h U 8 u I I G 3 E I T W s B g C M / w C m + O c F 6 c d + d j 0 V p w 8 p l j + A P n 8 w e C R I 0 w < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 5 E t 0 P 6 C A 4 c A J z e P u 6 v l r g w 9 m x x a C 0 4 + c w x / I H z + Q O D y Y 0 x < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " Y l N k 4 M N y I z d p c a c 4 B 1 4 k u X c 3 / a w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k p S R D 0 W v H i s a G u h D W W z 2 b R L N 5 u w O x F K 6 U / w 4 k E R r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g l Q K g 6 7 7 7 R T W 1 j c 2 t 4 r b p

FIG. 2 .
FIG. 2. Reactions increase free energy barrier ∆F of surrogate equilibrium model.(a) Radial concentration profiles c(r) minimizing the free energy F given by Eq. (6) at various fixed values of reaction coordinate x (colors correspond to panel b).(b) F as a function of x with ∆F indicated.(c) F (x) for various reaction rates k (colors correspond to panel d).(d) ∆F as a function of k. (a-d) Model parameters are cν = 0.18, L = 100 w, a2 = 250 νkBT , k = 0.0025 k0 (for panels a and b) and given in Fig. 1.
t e x i t s h a 1 _ b a s e 6 4 = " r o R Q 4 FIG. 3. Extended phase diagram accounting for reactions.(a) Free energy F of the surrogate equilibrium model approximated by Eq. (8) as a function of the nucleus radius R for decreasing concentrations c0 of the homogeneous state (bottom to top) and a2ν/a4 = 0.25.(b) Extended phase diagram indicating the stability of the homogeneous and droplet state as a function of c0 and a2ν/a4 in a passive (orange lines, k = 0) and active system (blue lines, k = 10 −3 k0).The droplet state is (meta-)stable right of the solid (binodal) lines, while the homogeneous state is stable left of the dashed (spinodal) lines.Behavior of numerical simulations (symbols) for k = 10 −3 k0 corroborate the results.Inset shows zoomed in region around the spinodal.(a-b) Additional model parameters are given in Fig. 1.