Heterogeneous Nucleation and Growth of Sessile Chemically Active Droplets

Droplets are essential for spatially controlling biomolecules in cells. To work properly, cells need to control the emergence and morphology of droplets. On the one hand, driven chemical reactions can affect droplets profoundly. For instance, reactions can control how droplets nucleate and how large they grow. On the other hand, droplets coexist with various organelles and other structures inside cells, which could affect their nucleation and morphology. To understand the interplay of these two aspects, we study a continuous field theory of active phase separation. Our numerical simulations reveal that reactions suppress nucleation while attractive walls enhance it. Intriguingly, these two effects are coupled, leading to shapes that deviate substantially from the spherical caps predicted for passive systems. These distortions result from anisotropic fluxes responding to the boundary conditions dictated by the Young-Dupr\'e equation. Interestingly, an electrostatic analogy of chemical reactions confirms these effects. We thus demonstrate how driven chemical reactions affect the emergence and morphology of droplets, which could be crucial for understanding biological cells and improving technical applications, e.g., in chemical engineering.

Droplets comprised of biomolecules, also known as biomolecular condensates, are crucial for organize biological cells.For example, such droplets separate molecules, control chemical reactions, and exert forces [1][2][3][4][5] .To fulfill these functions, it is likely that cells control the nucleation, location, size, and shape of droplets.While nucleation can happen spontaneously inside the cytoplasm 6 , most droplets might be nucleated heterogeneously involving other structures as nucleation sites.Indeed, many droplets interact with other structures inside cells, like cytoskeletal elements [6][7][8] , membrane-bound organelles 1,9,10 , and the plasma membrane [11][12][13] .Such interactions of liquid-like droplets with more solid-like structures is known as wetting 14 , and directly linked to heterogeneous nucleation [15][16][17] .However, most traditional examples of heterogeneous nucleation, e.g., by dust particles in clouds 18 , concern passive systems.In contrast, biological cells use external energy input to control processes actively, but it is unclear how activity affects heterogeneous nucleation and the properties of the subsequently forming droplets attached to the solid surface (sessile droplets).
Driven chemical reactions that affect the droplet material are one crucial example for an active process that controls droplets in cells 19,20 .If such reactions take place in the entire system, droplet size can be controlled [20][21][22] , droplets can divide spontaneously 23 , and homogeneous nucleation is suppressed 24 .Moreover, if reactions are restricted to the boundary, Liese and Zhao et al. recently demonstrated modified shapes of sessile droplets 25 .However, the generic case of wetting in the presence of bulk chemical reactions has not been considered so far.
To study how boundaries affect the behavior of chemically active droplets, we consider a continuous description of phase separation with driven reactions in the bulk and a passive interaction of the droplet material with a flat wall.Our stochastic simulations 26 and analytical results from an equilibrium surrogate model reveal that reactions generally suppress nucleation, whereas attractive walls facilitate them.However, both processes exhibit a complex interplay, which is connected to substantial shape deformations: While critical nuclei maintain their stereotypic spherical cap shapes, macroscopic droplets often exhibit elongated shapes along the boundary.

II. MODEL
We consider an isothermal system of fixed volume V filled with an incompressible, binary mixture of droplet and solvent material with equal molecular volume ν.The state of the fluid is described by the concentration field c(r, t) of the droplet material, whereas the solvent concentration is ν −1 − c(r, t).We describe the interactions and entropy in the system using a free energy comprised of bulk and surface terms given by 27 where f (c) is the local free energy density, κ penalizes compositional gradients, and g(c) is the contact potential describing the interaction of the fluid with the immobile boundary ∂V of the system.For simplicity, we focus on the linear order of the expansion g(c) = g 0 + g 1 c + O(c 2 ) since g 0 merely shifts the total free energy F , but does not affect the behavior 28 .In contrast, we describe the bulk interactions by where a 2 , a 4 > 0 are phenomenological coefficients.Equilibrium states minimize the free energy F .As a necessary condition, the variation of F with respect to c must vanish, which yields the two conditions where n denotes the outward normal of the surface ∂V .The first equation describes the balance of the exchange chemical potential µ = f ′ (c) − κ∇ 2 c, where the constant is determined from material conservation 29 .This generically yields a dense phase of concentration c out ≈ (2ν) −1 − a 2 /a 4 by an interface of width w = 2κ/a 2 .The dense phase typically assumes a spherical shape to minimize the surface tension γ ds = 2 2κa 3 2 /(3a 4 ) between the two phases 29 .In contrast, Eq. (3b) describes a boundary condition, which determines the behavior of the system close to the wall.In particular, the droplet material is repelled from the wall if g 1 < 0 14 .Since the droplet still maintains its spherical shape, its geometry at the wall is fully quantified by the contact angle ϑ, which is given by the Young-Dupré equation, cos(ϑ) = (γ ws − γ wd )/γ ds 15 , where γ ws , γ wd , and γ ds denote the surface tensions between wallsolvent, wall-droplet, and droplet-solvent, respectively.Since g(c) directly quantifies surface energies, we have This equation can only be solved for ϑ if the interactions between the droplet and the wall are weak, |g 1 | < g * with g * ≈ (2a 2 2 κ)/(9a 4 ).In contrast, the droplet will be repelled from the wall if g 1 < −g * , and it will fully wet the wall if g 1 > g * .Since the interesting process of heterogeneous nucleation is related to partial wetting where ϑ is defined, we concentrate on the case |g 1 | < g * .In particular, we consider the case g 1 > 0, corresponding to an attractive wall that is the most probable site for nucleation.
To describe nucleation, we next specify the dynamics of the system.We start with the continuity equation where j denotes the diffusive exchange flux between droplet material and the solvent, and the source term s describes chemical transitions 29 .The passive diffusive flux j is driven by gradients of the exchange chemical potential, j = −Λ d ∇µ + η, where Λ d is the diffusive mobility and η denotes diffusive thermal noise, which obeys ⟨η i (r, t)⟩ = 0 and the fluctuation dissipation the- where k B T is the thermal energy [30][31][32] .The system becomes active when we drive the reactions described by s out of equilibrium 21 .We focus on driven reactions that result in size-controlled droplets, which requires production of droplet material outside the droplet, while it is degrade inside 21,29 .This behavior can be captured by the linear expression where k sets the reaction rate and c 0 denotes the stationary state of the reaction scheme.We have shown previously that this case faithfully describes homogeneous nucleation in active systems and that the thermal noise associated with the reactions can be neglected since the diffusive noise η dominates on the length scales relevant for nucleation 24 .Taken together, our system is described by the stochastic partial differential equation augmented by the no-flux boundary condition n • j = 0, and Eq.(3b) describing local equilibrium at the boundary.

A. Numerical simulations reveal increased nucleation times
To investigate heterogeneous nucleation, we performed numerical simulations of Eq. ( 7) in a two-dimensional system.Here, we used finite differences to approximate derivatives and an Euler-Maruyama scheme to perform the time stepping 33 .We applied periodic boundary conditions along the x-direction and Eq.(3b) on both boundaries in the y-directions.Repeating the simulations many times, we observed that the time t nucl , quantifying when the first droplet nucleates, follows an exponential distribution; see Supplementary Fig. 9.We thus define the nucleation time τ as the ensemble average of t nucl .Fig. 1b shows that τ decreases for stronger droplet-wall attraction (larger g 1 ), as expected for nucleation of passive droplets 34 .In contrast, larger reaction rates k lead to longer nucleation times τ , indicating that active chemical reactions hinder nucleation, consistent with the results for homogeneous nucleation 24 .These numerical simulations indicate that repulsive walls (small g 1 ) and larger reactions (larger k) suppress heterogeneous nucleation, but they do not reveal the underlying principles governing the nucleation process.

B. Equilibrium surrogate model reveals trade-off between wall affinity and chemical reactions
To understand the nucleation path in detail, we next map the active system onto a surrogate equilibrium system with long-ranged interactions, which is possible for the special case of the linear reactions that we consider.Briefly, we can re-write the dynamics given by Eq. ( 7) as free energy functional where captures the energy associated with reactions [35][36][37] .Here, Ψ is the solution to the Poisson equation ∇ 2 Ψ = c 0 −c(r) and describes the long-ranged interactions, which originate from the interplay of chemical reactions and diffusion in the original model.Since the surrogate model requires mass conservation, we solve for Ψ employing Neumann boundary conditions, n.∇Ψ = 0. We used the free energy of the surrogate model, given by Eq. ( 8), to map out the minimal energy path connecting the homogeneous state with the state where a droplet wets the wall.To do this, we used constrained optimization to obtain concentration profiles at successively larger values of a nucleation coordinate X , which measures the amount of material inside the droplet; see SI, section B. Fig. 2a shows an example of such a minimal energy path, indicating that the maximal concentration inside the droplet increases alongside its size.
The minimal energy path allows us to evaluate the critical energy and the critical nucleus, which is the transition state of the nucleation process.Fig. 2b shows that chemical reactions generally increase F , potentially explaining why reactions suppress nucleation.However, reactions also increase F of the homogeneous state when the wall is attractive (g 1 > 0).This is because the boundary condition given by Eq. (3b) perturbs the homogeneous state in a layer with a thickness of roughly the interfacial width w, leading to local reactive fluxes.To see how chemical reactions affect nucleation, we thus evaluated the energy difference ∆ F between the transition state and the homogeneous initial state.Fig. 2c shows that ∆ F increases with the reaction rate k, while it decreases with more attractive walls (larger g 1 ), which is expected 24,34 .
We next quantified how chemical reactions interact with the wall affinity by decomposing the free energy barrier, where ∆F pas = ∆ F (g 1 , k = 0) quantifies the barrier for passive heterogeneous nucleation, ∆F react = ∆ F (g 1 = 0, k) is the energy barrier of chemically driven droplets at a neutral wall, and ∆F coupl (g 1 , k) denotes the energy due to the interaction of the two effects.Since we directly determined ∆ F (g 1 , k), ∆F pas (g 1 ), and ∆F react (k), we can infer ∆F coupl (g 1 , k) from Eq. (10).The data shown in Fig. 2d indicates a significant negative coupling between the wall affinity and the reactions, i.e., strong affinity to the wall can decrease the relative effect of chemical reactions.
The coupling between the wall affinity g 1 and the reaction rate k is also apparent in the contour lines of equal ∆ F , where reactions and wall affinity compensate each other; see Fig. 2c.We show in the SI, section D that concave contour lines would be expected without coupling (∆F coupl = 0), so that the slight convex shape of the observed lines is a strong indication for coupling.To analyze these contour lines in more detail, we approximate the coupling by a bilinear function, ∆F coupl (g 1 , k) = hg 1 k with pre-factor h, motivated by Fig. 2d.Moreover, we use a linear approximation for the reactive energy 24 , ∆F react = mk, and we express the energy associated with the passive case as 34 where ∆f = f (c 0 )−f (c eq in )+f ′ (c 0 )(c eq in −c 0 ) and ϑ is given by Eq. (4).Taken together, the contour line associated with energy ∆F c then satisfies We show in the SI, section D that the curvature of these contour lines is typically negative (∂ 2 k/∂g 2 1 > 0), consistent with the observed convex shape.
To explore the observed coupling further, we next ask how the energy F react , defined by Eq. ( 9), changes with the wall affinity g 1 .We thus initialized droplets of fixed volume for various affinities g 1 and numerically calculated F react inside the droplet; see SI, section D 4. Supplementary Fig. 11 shows that F react indeed decreases with larger g 1 , consistent with the negative effect we found above.However, this analysis only probes how different contact angles affect the reactions, whereas the reactions potentially also change the entire droplet shape away from a spherical cap.

IV. SESSILE ACTIVE DROPLETS SPREAD ALONG WALLS A. Droplets deviate from spherical cap shape after nucleation
The minimal energy paths that we obtained from the surrogate model not only provide energy barriers for nucleation but also the most likely droplet shape along the nucleation path.Fig. 2a shows that small droplets are shaped like a spherical cap, which corroborates the sampled critical droplet shapes shown in Fig. 1c.However, Fig. 2a also shows that larger droplets can deviate significantly from this equilibrium shape.
The droplet shape originates from a minimization of the free energy F of the surrogate equilibrium model.If reactions are absent, surface tension γ ds minimizes the interface between droplet and solvent material, resulting in a spherical shape with constant mean curvature.At the systems boundary, the shape must be compatible with the contact angle ϑ controlled by the energy balance that led to Eq. ( 4), which implies that attractive walls result in flatter spherical caps.The situation is more complicated with chemical reactions.The long-ranged interaction described by Eq. ( 9) leads to a repulsion of the droplet material, analogously to electrostatic repulsion of charged material.Consequently, different parts of the droplet repel each other 38 , which can induce spontaneous splitting for large bulk droplets 23,39 and explains why active droplets spread along the wall.

B. Interplay of reactions and wall affinity can deform macroscopic droplets
To understand how the interplay between wall affinity and driven chemical reaction influences the droplet shape, we next study the stationary shapes of sessile droplets numerically.We initialized droplets at the wall and simulated their time evolution for different reaction rates k and wall affinities g 1 until steady states were reached.To get an initial understanding, we first quantified droplet size, which is an area in our two-dimensional simulations.If droplet size is controlled by the reactiondiffusion length scale L = (D/k) 1/2 , with some diffusivity D, we expect the droplet size to scale with k −1 .Fig. 3a indeed reveals such a scaling, at least for large reaction rates.In contrast, for intermediate rates, we observe significant deviations, which increase with the wall affinity.The snapshots shown in Fig. 3b suggest that intermediate reaction rates k lead to strongly deformed droplets, which could explain the observed sizedependence.
We next quantified droplet shapes by analyzing the curve describing the interface, defined as the iso-contour c = 0.5.First, we determined the aspect ratio α, which is defined as the quotient of the lengths of the minor and major axis (see Fig. 3c).Second, we quantified the curvature K of the interface (see SI, section C) and determined its variation along the interface (see Fig. 3d) as well as the curvature at the center point (see Fig. 3e).All three quantifications shown reveal the same fundamental dependence: Droplets are essentially spherical caps when reactions are absent (k = 0), consistent with expectations 14 .Moreover, droplets are spherical for neutral walls (g 1 = 0), even when reactions are strong, because all fields are radially symmetric in this symmetric case.Interestingly, droplets also exhibit a spherical cap shape for large reaction rates, likely because strong reactions reduce droplet size so interfacial effects dominate.Between these extreme values, we observe strong deviations from a spherical shape, both in the snapshots (Fig. 3b) and in the quantifications (Fig. 3c-e).In fact, for a given wall attraction g 1 , we empirically find an intermediate rate k for which the deformations are maximally, and this rate increases with g 1 .Taken together, these quantifications reveal an interplay between wall attraction and reactions, which leads to macroscopic deformations of sessile droplets.
The macroscopic deformations induced by the reactions affect the entire interface and could thus also impact the contact angle ϑ.To test this, we evaluated the slope −∂ x c/∂ y c along the implicit curve c(x I , y I ) = 0.5 and measured the associated ϑ at the boundary.The inset in Fig. 4a shows that ϑ decreases with increasing wall affinity g 1 , consistent with the prediction from Eq. ( 4).However, the data in Fig. 4a also indicates that ϑ decreases for larger reaction rates k, indicating an effect of the reactions.For attractive walls (g 1 > 0), the contact angle first declines sharply for increasing k, although this decline becomes weaker for larger k and there is even a brief non-monotonic behavior.The value of k, where this non-monotonic behavior is observed, depends on g 1 and coincides with the maximal shape deformation of the droplets (compare Fig. 4b to Fig. 3a).Chemically active sessile droplets thus exhibit changes in apparent contact angle concomitantly with global shape deformations.

C. Anisotropic fluxes cause droplet deformation
To gain further insights into shape deformations, we next analyze how reactions disturb the spherical cap shapes expected for passive droplets.We consider a spherical cap described by a radius R and contact angle ϑ in the half-space y > 0. To describe the droplet shape explicitly, we use an effective droplet model 29 , assuming a thin interface (w ≪ R), so we can approximate the dynamics of the concentration fields c i by reaction-diffusion equations inside (i = in) and outside (i = out) the droplet; see SI, section E. Here, we linearized Eq. ( 7) without noise around the concentrations c i )] quantifies surface tension effects for i = in, out 29 .For simplicity, we consider the quasi-stationary situation where the concentration fields equilibrate faster than the interfacial shape.The resulting stationary state solutions to Eq. ( 13) then depend on the reaction-diffusion length scale L = D/k, the average concentration c 0 , and the boundary conditions.The general solution of this Helmholz equation reads where I λn (z) and K λn (z) are the modified Bessel functions of first and second kind.Here, the wavelength λ n = nπ/ϑ of the polar coordinate is quantized by the boundary conditions, which can also be used to determine the series coefficients A i n and B i n as described in the SI, section E. Taken together, this approach provides approximate solutions for the concentration profiles inside and outside the droplet.
The concentration fields imply fluxes j i = −D∇c i , which can affect droplet shapes.Indeed, the interface speed v n normal to the interface reads 29 where n is the normal vector of the interface.Using Eq. ( 14), we find with This interfacial speed v n quantifies how chemical reactions would disturb the spherical cap shape.
We start by examining the shape of a droplet on a neutral wall (g 1 = 0).Fig. 5a shows that only the n = 0 mode contributes, whereas A i n = B i n = 0 for n ≥ 1, consistent with the observed spherical shapes.The zeroth mode is essentially unchanged for an attractive wall (g 1 > 0, ϑ < π/2, Fig. 5b), but the first and second modes become important for larger droplet sizes, indicating non-spherical droplet shapes.The dependence of these modes on the droplet radius R captures the behavior we observed so far: For very small radii, the zeroth mode is negative, indicating an unstable size.Consequently, droplets can only grow spontaneously when R exceed the first root of C 0 (R), which thus corresponds to the size of the critical nucleus.At the critical size, the higher modes shown in the right two panels in Fig. 5b are vanishingly small, consistent with the spherical shape of critical nuclei observed in Fig. 1c.Droplets larger than the critical size grow until they reach a stationary state, marked by the second root of C 0 (R) in Fig. 5b; see also ref. 22 .At this stationary state, the zeroth mode does not contribute to the dynamics, but higher-order modes, which are connected to shape deformations, do.Indeed, Fig. 5c shows that the interface speed of a stationary droplet with radius R * , chosen such that C 0 (R * ) = 0, is such that the droplet flattens and spreads along the walls.Taken together, the analysis of the first mode of a droplet with a stationary size indicates that fluxes caused by the driven reactions deform the droplet.
To understand droplet deformations in detail, we next focus on the first mode (n = 1).Fig. 5d shows that its magnitude |C 1 (R = R * , k)| vanishes when reaction are absent (k = 0), consistent with the expected spherical cap shape of passive droplets.Larger reaction rates k first cause |C 1 | to increase, but |C 1 | then decreases beyond a critical value of k.This non-monotonic behavior is qualitatively consistent with the shape deformations shown in Fig. 3 and might be related to the steady state size of reactive droplets: Higher reaction rates k reduce the droplet size (see Fig. 3a and ref. 22 ), leading to smaller magnitudes of the droplet deforming modes; see Fig. 5b.In summary, chemical reactions combined with symmetry breaking by an attractive wall lead to non-isotropic flows that cause droplet deformation.

D. Reactions deform sessile 3D droplets
So far, we analyzed deformed droplets only in two spatial dimensions.To check whether the observed effects persist in three dimensions, we performed a simulation for the parameter regime where droplets are deformed.Supplementary Fig. 13 shows that such droplets are circular in the xy-plane and are deformed only in the zdirection.Consequently, simulations using cylindrical symmetry are suitable to describe the problem.Fig. 6 shows that we find spherical cap shapes for small and large reaction rates k, and strongly deformed droplets between these extremes, consistent with our results for attractive walls in two-dimensional systems.Taken together, we thus expect that our results translate directly to three-dimensional systems.

V. ELECTROSTATIC ANALOGY EXPLAINS EFFECTS OF DRIVEN REACTIONS
So far, we have focused on simulating the dynamical equation (7) and analyzed the simplified reactiondiffusion equation (13).In this final section, we now investigate the equilibrium surrogate model given by Eq. ( 8) in detail, which will allow us to interpret the reactions as long-ranged electrostatic interactions.
The energy F react associated with reactions, given by Eq. ( 9), reveals that the deviation of the concentration field c(r) from the reaction equilibrium c 0 can be interpreted as a charge density.This interpretation only works for the linearized reactions given by Eq. ( 6) and when the average concentration field is equal to c 0 , e.g., when c dV = c 0 V , which implies a charge neutral system in the electrostatic interpretation.Fig. 7(a) shows a typical concentration profile of a sessile droplet, indicating that the droplet can be interpreted as a positively charged ball surrounded by a cloud of negative charges, so that the entire system is charge neutral.The reactiondiffusion length scale L = D/k controls the extent of the cloud and the reaction rate k determines the magnitude of the electrostatic interaction.Even though the sessile droplet in Fig. 7(a) looks as if it would form a dipole with the surrounding negative charges, this is not the case since image charges restore the symmetry.Distant droplets thus hardly interact with each other.
To be more quantitative, we next calculate the total energy of a single sessile droplet in the limit of a thin interface, also known as a capillary approximation; see SI, section G.The approximate expressions for small droplets in two and three dimensions read 24,37 F2D ≈ 0.029 + 0.32 log(LV where the bulk energy proportional to ∆f = f (c 0 ) − f (c in ) + f ′ (c 0 )(c in − c 0 ) scales with the droplet volume V , whereas the surface energy proportional to γ ds scales with the effective size of the interface, which depends on the contact angle ϑ: A 2D = 2V 2D (2ϑ − sin(2ϑ)) and 3 sin( ϑ 2 ) 4 .For simplicity, we neglect the influence of ϑ on the electrostatic energies associated with the charge distribution given by the respective first terms proportional to k.Instead, we use the expression for one half of a spherical droplet on a wall as an approximation; see SI, section G. Neglecting logarithmic corrections, we thus find that reactive energies scale as R 2+d , bulk energies as R d , and surface energies as R d−1 , when droplets are small and d is the space dimension.
We first use the approximate energies to investigate nucleation.The surface energy dominates for small droplets, explaining why nucleated droplets close to the critical size are spherical; see Fig. 1c.Eqs. ( 18) also show that chemical reactions raise the nucleation barrier (see Fig. 7(b,c)), consistent with suppressed nucleation.The associated energy barrier ∆ F shown in Fig. 7d exhibits a very similar dependence on the reaction rate k and the wall affinity g 1 compared to the numerically determined barriers shown in Fig. 2c.However, since we do not consider the contact angle ϑ in the simplified electrostatic energy given in Eq. ( 18), this approximation cannot capture the coupling between k and g 1 we revealed above.Moreover, the capillary approximation used here is known to overestimates the energy barrier for droplet formation 34 since it assumes a well-defined separation of the droplet and the surrounding dilute phase during nucleation.Yet, we can conclude that the electrostatics do not affect the shape of small droplets, but they oppose their formation, thus lowering nucleation rates.
To discuss shape deformations and splitting of droplets, we next approximate these two situations by considering two sessile droplets with total volume V tot that either touch each other or are well-separated; see Fig. 8a.If F (V tot ) denotes the energy of the single droplet of volume V tot , the two well-separated droplets have a total energy of 2 F ( 1 2 V tot ).In addition, the two connected droplets exhibit an electrostatic repulsion, which we approximate as the repulsion between two point charges; see SI, section G. Since the approximations given in Eqs.(18)  only hold for small droplets, we here obtain F (V ) from Eq. (G21) and Eq.(G35) in the Supplement.Taken together, we can thus compare the energies of the two droplets at different separation to the energy of a single droplet in 2D and 3D for various V tot ; see Fig. 8(b,c).Without reactions (k = 0), the energy of two droplets is always higher than that of a single droplet, consistent with a minimization of surface energies and absent droplet deformations.Moreover, F (V ) decreases monotonously beyond the critical size, implying that droplets grow until they are limited by system size.The picture changes qualitatively when reactions are enabled (and droplets become charged in the electrostatic picture): Even for arbitrarily small k, F (V ) develops a minimum at a finite size V and the electrostatic energy makes larger droplets unfavorable; see Fig. 8(b,c).Since the minimal energy is negative (implying that the droplet is favored over the homogeneous state), we observe that the states with two droplets (orange curves) have lower minima than the one-droplet states (blue curves), which implies that a large system would prefer to have many droplets.A large droplet count could emerge either from nucleating additional droplets or from splitting existing droplets.
To investigate droplet splitting, we compare the minimal energy of a single droplet with the energy of the two-droplet states at the corresponding total volume V tot .Fig. 8b shows that the two-droplet states can have lower energy for small reaction rates k, suggesting that a droplet with this volume might deform spontaneously.In contrast, the right panel shows that a single droplet has lower energy for large reaction rates k, suggesting that the single droplet is stable.The transition between these cases marks the onset of droplet deformation, which depends on the reaction rate k and the wall affinity g 1 .Fig. 8d suggests that droplets deform if reactions are not too strong and larger g 1 favor deformations; similar to Fig. 3(c-e).However, in contrast to Fig. 3(c-e), the simple scaling theory predicts that the deformed droplets (mimicked by the two adjacent droplets) are favored even when g 1 = 0.It is likely that deformations are kinetically suppressed, i.e., that small deformations are energetically unfavorable and we thus do not observe them.
The predicted transition from a single to a deformed and then split droplet can be understood as a competition between surface energies (favoring compact droplets) and effective electrostatic repulsion (favoring elongated droplets) in conjunction with the favored finite droplet size due to reactions.Droplets are spherical at large reaction rates k since large k implies small droplets where surface energies dominate.Smaller reactions imply larger droplets, where the electrostatic effects are stronger and eventually favor deformed droplets that can also split.However, if reactions are absent, we again observe spherical droplets since there are no electrostatic effects.
In summary, interpreting the reactions using electrostatics explains the delayed nucleation and droplet deformation.Nucleation is suppressed since electrostatics disfavor an accumulation of charges, but the critical droplets are still spherical since surface energies dominate.In contrast, larger droplet can minimize the total energy by deforming, which enlarges the average distance between charges and thus lowers the electrostatic energy at the expense of a larger surface energy.In this case, the deformed droplets can also grow beyond the volume predicted for spherical shapes, consistent with Fig. 3a.

VI. DISCUSSION
We showed that nucleation is accelerated by attractive walls and suppressed by active chemical reactions, consistent with our previous results 24 .However, we also found an intricate interplay between the two processes, likely stemming from shape modifications in sessile droplets.While the shape of critical nuclei is in agreement with a spherical cap, consistent with passive scenarios 14 , we found massive deformations for larger droplets at inter-mediate reaction rates.Within the surrogate model, the elongated droplets can also be interpreted as a trade-off between surface tension and effective electrostatic repulsion, which demonstrates that the active chemical reactions mediate a repulsive long-ranged interaction.
Our study of sessile active droplets provides the first step toward understanding the interplay of chemically active droplets with other structures.For simplicity, we focused on a two-component mixture, whereas cells exhibit a staggering complexity involving thousands of different components, which could be described by a multicomponent extension of our theory 40,41 .Furthermore, we considered reactions that depend linearly on composition, but realistic, thermodynamically consistent reactions are more complex 42 .While our previous work suggests that linear reactions capture nucleation quantitatively 24 , the shape of macroscopic droplets might change drastically, e.g., when surface reactions are additionally taken into account 25 .Finally, we only studied flat walls, but boundaries in cells are typically curved and deformable 10,43 .Incorporating these aspects will likely require more advanced computational methods (such as 44 ), but we expect the general behavior that we unveiled here to persist.
Active control of phase separation is a challenge in cells.Our work suggests that cells could use chemical reactions to fine-tune the rate of heterogeneous nucleation as well as the shape of sessile droplets.Such deformed droplets may offer better control of wall deformations by condensation, e.g., by using less material than a spherical cap to achieve the same deformation.Such a regulation process would be a fascinating example combining control via active chemistry and a membrane surface 45 .Moreover, our model serves as an intriguing example of boundary effects in active field theories, inviting comparisons with other active field theories 46,47 .where Γ = 10.We determine the minimal energy path by minimizing the free energy F (given in Eq. 9 of the main text) with constrained values of X .We impose a value X 0 of the reaction coordinate using a Lagrange multiplier λ and minimize by evolving the partial differential equations which corresponds to conserved and non-conserved dynamics with mobilities Λ D and Λ L , respectively.We use Λ D = 1 and Λ L = 100, which proved a good compromise between speed and stability.Using this procedure, the profiles c(r) that minimize F for each value X 0 of the constraint comprise the minimal free energy path.The profile with the largest energy F corresponds to the saddle point and thus the critical nucleus.
where ∆f = f (c 0 ) − f (c eq in ) + ∂ c f (c 0 )(c eq in − c 0 ), A ds is the droplet solvent interface, A wd is the wall droplet interface and A w is the total area of the wall.Since γ ds cos(ϑ) = γ ws − γ wd , we can rewrite this expression as The last term only gives a constant offset and will be dropped in the following.For a spherical cap geometry in two dimensions, we obtain and a critical droplet radius of R crit = γ ds /∆f .The resulting energy barrier is then given by The wetting angle ϑ can also be expressed as a function of our microscopic parameters, ϑ = arccos 3g1 a2 a4 2κ .

Augmented energy barrier by chemical reactions
The simplest assumption for how chemical reactions and wall affinities affect the energy barrier between homogeneous and droplet state is a simple superposition of the two.We consequently propose where we for simplicity assumed a linear relation for the reactive energy barrier, ∆F react.= mk, in line with our previous results 24 .We now want to look at contour lines of the energy barrier (where chemical reactions and wall affinity compensate each other) and obtain The second derivative is given by where d = 3 a2 a4 2κ .Consequently, the contour lines are concave functions of g 1 , whereas the contours in the data are slightly convex.

Additional coupling
We measured the slope of the energy barrier m for different values of g 1 and found a linear decrease (see Fig. 10).
To capture this observation we add the lowest order coupling between reactions and wall affinities, and we obtain for the contour line The second derivative in the limit of small wall affinities is then given by which for typical parameters is positive and thus shows concave behavior in agreement with the data for h < 0.

Effect of spherical cap shape on reactive energy
To test the effect of different contact angles ϑ on the reactive energy, we initialized droplets for different wall affinities g 1 .We then evaluated the reactive energy by integrating Eq. 9, neglecting the contribution from the bulk phase.The resulting reactive energy as a function of the wall affinity is depicted in Fig. 11.The diffusivity out ) is the same in both phases in our model.Note that we also neglected the fourth order derivative since we assume low concentration variations.The geometry of the problem is depicted in Fig. 12.The expected droplet shape follows a circular segment (corresponding to the spherical cap in 3D) with an opening angle of 2ϑ, where ϑ is the contact angle.The distance between the origin of the circular segment and the droplet wall interface is given by y 0 = R cos(ϑ).Since we assume fast relaxation of the concentration fields inside each phase, we seek a steady state solution of Eq. (E1).For the steady state equation (∂ t c = 0), it is instructive to define the reaction-diffusion length scale L = D/k.The solution for the concentration field in the given polar geometry is given by 50 Here, I λn (z) and K λn are modified of first and second kind, and λ n is the separation variable which has to fulfill λ n = nπ/ϑ to be consistent with the boundary conditions.We next will determine the coefficients A i n and B i n from the boundary conditions.

Concentration field inside the droplet
The relevant boundary conditions are Translating the first boundary condition to polar coordinates gives ∂c(r(x, y), φ(x, y) ∂y We rewrite the second boundary condition as with Next, we define the projection operator for a well behaving function h m (φ), which will be defined explicitly later.Applying the projection operator to the boundary condition leads to We rewrite this expression as with coefficients Consequently, we recast the boundary condition into a system of linear equations.The number of equations is in principle infinite, but we can hope to obtain an accurate approximation by looking at finite number m of modes.We want to treat the droplet-wall boundary condition in a similar way.We can insert the solution for the concentration field into the boundary condition and obtain which we rewrite as with After applying the projection operator P , we then obtain another system of linear equations where the coefficients are given by Finally, the only task left to do is to solve the two matrix equations We can combine both equations to and invert the full matrix numerically.We choose h m (φ) = cos(πmφ/ϑ) in the projection operator and m ∈ [0, 3].The projection method we used before is not well suited for this situation since the first and the last boundary condition are defined over φ ∈ [−ϑ, ϑ], whereas the second one is defined over φ ∈ [−π/2, −ϑ] ∪ [ϑ, π/2].We can define the operators over the specified range of φ and get three equations for two unknowns, which gives an overdetermined system.One way around this is to consider a domain given by the opening angle 2ϑ.If we use the same quantization of the separation variable λ = πn/ϑ, we can solve the problem with the two boundary conditions which are defined over the opening angle.Using the same operator as defined for the inside concentration we get the two equations where with We can solve for the coefficients as before by inverting the combined matrix.Note that b out is defined as before but with a different constant value.
3. Special case of a neutral wall In the case of a neutral wall the opening angle is 2ϑ = π, and we obtain a slightly different matrix equation for the coefficient, which needs extra treatment.Since we assume constant total volume, we can relate the radii of the two smaller droplets (R 1 and R 2 ) to the radius of the single droplet by The distance between the centers of the two droplets is given by l = 2R 1 sin(ϑ) for the adjacent droplets since we assume that they touch each other.For the separated droplets we assume no interaction between them.We consider a wall of length L sys .

Energy contributions for a single droplet
We expect that the behavior of the system is governed by three energy contributions: interfacial energy, volume energy, and a contribution that captures the influence of reactions.The volume energy is given by The interfacial energy is given by E int = γ ds a + γ wd c + γ ws (L sys − c).Since γ ds cos(ϑ) = γ ws − γ wd , we generally obtain E int = γ ds (a − c cos(ϑ)) + γ ws L sys , which for the concrete systems read We can also express the interfacial energy as a function of the volume by substituting R = 2V 2D /(2ϑ − sin(2ϑ) and The reactive energy can be determined from the electrostatic analogy, but we do not have a precise expression.To derive a scaling relation, we use the electrostatic energy of a droplet on a neutral wall, which serves as an upper bound for the energy of a spherical cap.
a. Energy of single chemically active droplet in dilute phase in two-dimensional system As established in III B, we can map the chemically active system onto an electrostatic system with the energy with We want to solve this equation and consequently need an expression for the concentration field, which we can approximate as solution to reaction-diffusion equations in stationary state, i.e., we assume quasi-stationarity.The solutions, with vanishing derivative at the origin and infinity and c in out , are given by where L = D/k is the reaction-diffusion length scale as defined in the main text.These concentration profiles reproduce the main features of the solution to the full phase field equations.However, they do not conserve mass in the current form, which is easily seen by calculating the integrated reactive current inside and outside, This breaking of mass conservation is inconsistent with the initial equations Eq. ( 5) and Eq. ( 6) in the main text, which conserve the overall mass if the average density in the system is c 0 .In contrast, the thin-interface model with quasi-stationary that we here use for simplicity breaks this mass conservation.Since this is a crucial feature of our system, we correct the mass by demanding a modified equilibrium concentration outside the droplet, which we use in the following.We can then solve the energy integral by rewriting it as The function c − c 0 is not continuous, so we integrate piecewise by parts Using polar coordinates and the fact that we have Neumann conditions on Ψ at R = 0 and R → ∞, we find We next define the potential Ω = ∂ r Ψ, plug it into the Poisson equation given by Eq. (G6), and solve for Ω, where d 2 is an integration constant.The constant d 2 needs to be zero since the field needs to be divergence-free at r = 0. Similarly, we obtain We immediately see that with the modified outside concentration we have Ω in (R) = Ω out (R) and the boundary term in Eq. (G14) vanish.Importantly, the boundary terms do not vanish by themselves but only after integration over full space.We thus find We then integrate the inside contribution of the energy, where the last approximation holds for small droplet radii R and is in general a lower bound to the entire function.
For the outside contribution we get To determine the scaling behavior of this expression with R, we need to plug in the modified outside concentration Consequently the total energy scaling is given by where γ euler ≈ 0.577 is the Euler number.This expression is again valid for small droplet radii but the approximation becomes unphysical as one approaches the reaction-diffusion length scale L = D/k.The approximation then becomes negative, whereas the full expression approaches an almost linear positive increase.
b. Energy of single chemically active droplet in dilute phase in three-dimensional system We perform the same calculations as above for three dimensions, which will give slight modifications.We again start from the two reaction-diffusion equations We can write these equations in spherical coordinates and then solve for the stationary state where we expanded for small droplet radii in the last line.The outside energy can be expressed as We then get where we expanded for small droplet radii in the last line.The total energy is thus given by We now want to use the two expressions above to estimate the reactive energy F sph.cap react for sessile droplets forming a spherical cap on a wall.The volume of a spherical cap translates to its radius as R 3D (V ) = (3V 3D /(2π)) We can now insert these identities into the full expressions for the energy of a sessile droplet on an attractive wall to obtain an upper estimate for the energy of a neutral spherical cap.Note that we need to correct this energy by a factor of 1 2 since we calculated the energy of a full droplets.Hence, where F tot react is given by Eq. (G21) for 2D and by Eq. (G35) for 3D.

Energy contributions for two droplets
For the two-droplet scenarios, the energy scaling with the volume remains unchanged.For the interface contribution, we now need to consider two contributions with half of the total volume.We thus get (2/3) sin 4 (ϑ/2) . (G38b) The reactive energy can be calculated as follows b 3D kV (5/3) .(G39b)

Energy of two interacting chemically active droplets
We next consider an additional energy contribution to describe the repulsion between the two adjacent droplets.To keep it simple, we approximate the two droplets as point charges with separation l 2D = 2R 1 sin(ϑ) = √ 2R sin(ϑ), and l 3D = 2R 1 sin(ϑ) = (2/ 3 √ 2)R sin(ϑ).The total charge inside each droplet is given by q = (c 0 −c(r))dV ≈ (c 0 −c in )V /2 and the potential is where we assume that the connection between the two droplet centers is along the x-axis and one charge is located at the origin.The associated energy to bring a charge from the boundary of the system to the position l is then in 2D given by In three dimensions, we can safely assume an infinite system, such that the reference point of the potential vanishes.We consequently get (G42)

FIG. 1 .
FIG. 1. Driven reactions delay heterogeneous nucleation.(a) Schematic picture of a sessile reactive droplet.Droplet material is attracted to the wall and itself, leading to phase separation from solvent.In addition, active chemical reactions convert droplet material to solvent material inside the droplet, while the converse happens outside.This process leads to long-ranged diffusive fluxes, resulting in deformed macroscopic droplets.(b) Numerically determined nucleation times τ as a function of reaction rate k for various wall affinities g1.(c) Ensemble of critical nuclei for two reaction rates k and positive wall affinity g1 = 0.025 a2/w.(b-c) Simulations are detailed in the SI, section A; Model parameters are νc0 = 0.1, a2/(νkBT ) = 200, ν = w 2 , a4 = 4a2ν, and k0 = Λ d a2w −2 .

10 X /w 2 FIG. 2 .
FIG. 2. Reactions raise energy barrier of nucleation.(a) Cuts through concentration profile c(x, y) along the two axes (top) and shape of droplet interface (bottom) along the nucleation path (colors) for k/k0 = 0.01 and g1w/a2 = 0.05.The shape of the critical nucleus is indicated in red.(b) Energy F of the surrogate model as a function of the nucleation coordinate X for various reaction rates k for a neutral wall (g1 = 0) and an attractive wall (g1 = 0.05 a2/w).(c) Energy barrier ∆ F as a function of g1 and k.(d) Energy F coupl following from Eq. 10 as a function of g1 and k. (a-d) Model parameters are as in Fig. 1.

FIG. 3 .
FIG. 3. Active droplets spread along walls.(a) Droplet area as a function of the reaction rate k for different wall affinities g1.(b) Concentration fields of sessile droplets for two values of g1 and various k.(c) Aspect ratio α of sessile droplet as a function of g1 and k.(d) Standard deviation of the interface curvature K as a function of g1 and k.(e) K evaluated at the droplet center as a function of g1 and k. (ae) Model parameters given in Fig. 1.

FIG. 4 .
FIG. 4. Reactions affect contact angles.(a, b) Deviation of contact angle ϑ from passive case (k = 0) as a function of reaction rate k for various wall affinities g1.The inset shows ϑ as a function of g1 for k = 0. Model parameters given in Fig. 1.
out , so the diffusivity is given by D= Λ d f ′′ (c (0) in ) = Λ d f ′′ (c (0)out ) for our choice of a symmetric free energy 29 .To solve the reactiondiffusion problem, we employ polar coordinates (r, φ) centered at the sphere describing the spherical cap.We impose ∂ y c i | y=0 = 0 and ∂ r c out | r→∞ = 0 at the system boundary.The conditions at the droplet interface are governed by the local phase equilibrium and read c i

FIG. 5 .
FIG. 5. Reactions cause non-isotropic interface speeds.(a) Series coefficients Cn as a function of droplet size R for n = 0, 1, 2 at a neutral wall (g1 = 0, contact angle ϑ = 90 • ) for k = 0.0025.(b) Cn as a function of R for an attractive wall (ϑ = 60 • ).(c) Interface velocity vnn (red arrows) predicted from Eq. (16) using the first three modes for a spherical cap of R/w = 7 for ϑ = 60 • (black line).(d) |C1| as a function of the reaction rate k for ϑ = 60 • .(a-d) Additional parameters are derived from values given in Fig. 1.

FIG. 7 .
FIG. 7. Electrostatic analogy explains nucleation behavior.(a) Effective charge density c(r) − c0 of a sessile droplet for reaction rate k/k0 = 0.005.(b, c) Energy barrier F as a function of the volume V in 2D (panel b, k/k0 = 0.001) and 3D (panel c, k/k0 = 0.01) for various approximations: The full theory (solid blue lines, see Eq. (G21) and Eq.(G35)) is compared to the approximations given in Eqs.(18) (dashed lines) and the expression for k = 0 (black lines).The contact angle is ϑ = 80 • .(d) Nucleation barrier ∆ F determined from maximizing F (V ) given by Eq. (18a) for various k and g1.(a-d) Additional parameters are given in Fig. 1.

FIG. 8 .
FIG. 8. Electrostatic analogy explains droplet deformations.(a) Effective charge density c(r) − c0 of two sessile droplet that are well-separated (left) or directly adjacent (right) for for a reaction rate k/k0 = 0.005.(b, c) Energy barrier F as a function of the total volume Vtot of one droplet (I, blue lines) and two droplets that are well-separated (II, orange lines) or adjacent (II-c, green lines) in 2D (b) and 3D (c) determined using the full theory; see SI, section G. Parameters are indicated in the panels.(d) State diagram indicating which state has the lowest energy F at the total volume where a single droplet is stable as a function of reaction rate k and wall affinity g1.The single droplet has the lowest energy in the blue region and for k = 0, whereas the separated droplets have lower energy in the orange region, and both connected and separated droplets have lower energies in the green region.(a-d) Additional parameters are given in Fig. 1.

FIG. 9 .
FIG.9.Distribution of nucleation times for different wall affinities g1 and constant non-zero reaction rate k0.

g 50 FIG. 11
FIG. 10.(a) Relative increase of the energy barrier ∆ F with the reaction rate k for different wall affinities g1.(b) Fitted slope of the data in panel a as a function of the wall affinity g1.

2 .
Concentration field outside the droplet For the concentration profile outside we have in principle three boundary conditions c(R) = c

1 L 2 2 (
(c − c 0 ) = 0 .(G23)With boundary conditions ∂rc in | r=0 = 0, c(R) = c in the solution is given byc in = R(c (0) in − c 0 ) sinh r L csch R L r + c 0 .(G24)For the outside, with boundary conditions lim r→∞ ∂ r c in | r=0 = 0 and c(R) = c fields, we can next calculate the integrated reactive flux for the insideout − c 0 )R exp(R/L)L 2 .(G27)To balance the two fluxes, which ensures the solvability of the Poisson problem, we obtain a modified outside concentration,in − c 0 ) sinh(R/L) R cosh R L − L sinh R L exp(−R/L) L + c 0 , = e − R L c (0) in − c 0 L − R coth R L L + c 0 .(G28)Next, we need to calculate the potential Ω = ∂ r Ψ.We get for the inside Ω in = − r cosh(r/L) − L sinh(r/L))