Reaction-Driven Diffusiophoresis of Liquid Condensates: Potential Mechanisms for Intracellular Organization

The cellular environment, characterized by its intricate composition and spatial organization, hosts a variety of organelles, ranging from membrane-bound ones to membraneless structures that are formed through liquid–liquid phase separation. Cells show precise control over the position of such condensates. We demonstrate that organelle movement in external concentration gradients, diffusiophoresis, is distinct from the one of colloids because fluxes can remain finite inside the liquid-phase droplets and movement of the latter arises from incompressibility. Within cellular domains diffusiophoresis naturally arises from biochemical reactions that are driven by a chemical fuel and produce waste. Simulations and analytical arguments within a minimal model of reaction-driven phase separation reveal that the directed movement stems from two contributions: Fuel and waste are refilled or extracted at the boundary, resulting in concentration gradients, which (i) induce product fluxes via incompressibility and (ii) result in an asymmetric forward reaction in the droplet’s surroundings (as well as asymmetric backward reaction inside the droplet), thereby shifting the droplet’s position. We show that the former contribution dominates and sets the direction of the movement, toward or away from fuel source and waste sink, depending on the product molecules’ affinity toward fuel and waste, respectively. The mechanism thus provides a simple means to organize condensates with different composition. Particle-based simulations and systems with more complex reaction cycles corroborate the robustness and universality of this mechanism.

with the bare fluxes With this form of coupling the bare fluxes, the incompressibility constraint is fulfilled.Now we continue with the evaluation of droplet velocities, starting with a 1D droplet in concentration gradients.In this section, we consider the dynamics in the vicinity of an individual droplet and omit the index i running over these.We assume that the droplet with radius R positioned at z cm = 0 at time t = 0 accumulates product material from its closest surroundings (i.e., 1D analog of Voronoi cell), Ω = [z Ω − , z Ω + ].Since the droplet moves with velocity, v, we exploit that the concentration profiles stay roughly constant in the (instantaneously) co-moving frame, up to small changes upon moving through the concentration gradient, and denote the stationary quantities by a prime, e. g., ϕ c (z, t) = ϕ ′ c (z ′ ) for z ′ := z − vt.We assume interfaces to be infinitely sharp.Then, the integral of ϕ P (z, t) over the droplet domain, ω = [−R+vt, R+vt] yields the total droplet material, which we assume to be instantaneously conserved in the co-moving frame.This corresponds to the assumption that the droplet instantaneously adopts the size that balances loss of product material on the inside and total production in the accumulation volume.Integrating over a domain larger than the droplet domain, in order to include the interfaces and discontinuities in the fluxes, then taking the limit we obtain At this point, we make the approximation (ICA) that product material, which is created inside the accumulation volume, instantaneously condenses onto the droplet's interface.Moreover, the droplet shape remains spherical independent of the decay of droplet material on the inside and the condensation on the outside.
In 1D, this modifies the time-evolution, Equation 2, for each droplet region, including its accumulation volume, to where δ(•) denotes Dirac's delta function and ζ(z − z cm ) is an assignment function, which accounts for the finite lifetime of product molecules, i. e., the ratio of the rate at which product molecules reach the interface and the production rate in the accumulation volume.
It is estimated for the 3D-case below.Since the driven reactions are now replaced by surface sources and sinks, we retain only the 'non-reactive' flux, i. e., the precursor-flux contribution to the product flux, −j R , is replaced by surface sources, using that outside the droplet, r / ∈ ω, the product concentration almost vanishes ϕ P ≈ 0.
With this, the above calculation becomes where the second term was integrated by parts.Taking the limit, then results in an expression that quantifies the velocity in terms of integrals over the product concentration and the fluxes.
There are two contributions: The first arises from the fluxes inside the droplet and the second stems from the reaction/condensation in the droplet's accumulation volume.The second contribution is illustrated in Figure S1.Let us generalize the above calculation to 3D.We assume a spherical droplet of radius R moving at velocity v = vê z , where êz is the unit vector, pointing in z-direction.The droplet domain is described by ω = {r| ||r − vt|| < R}, and its surface denoted by ∂ω.Again, we approximate all product that is produced within the accumulation volume to condense instantaneously onto the droplet's surface.This modifies the evolution equation to where the co-moving position is r ′ = r − vt, and r ′ , θ ′ , φ ′ denote the spherical coordinates with respect to the droplet's center.We assume the droplet morphology (shape and size) be stationary in the co-moving frame.Here and in the following, all functions of the argument r ′ are denoted by a prime.The fraction that appears in the integral is the result of the fraction of area elements from accumulation in the domain, r ′ 2 sin θ ′ dθ ′ dφ ′ , and deposition at the sphere's surface, R 2 sin θ ′ dθ ′ dφ ′ .Hence, first, we assume that molecules reach the droplet at the angle, at which they have been produced in the accumulation volume.More accurately, the molecules will perform a random motion, before reaching the droplets surface.
We will phenomenologically reintroduce this effect below and account for it in the assignment function ζ(|r − r cm |).Hence in 3D, we have where the boundary term vanishes, having formally taken a spherical integration domain with radius R → R + ϵ with a small ϵ > 0 and the subsequent limit ϵ → 0, as was explicitly shown in the 1D case.We obtain the final estimate for the velocity, where j c := j c • êz and again, θ is the angle between r − r cm and the z-axis, as illustrated in vicinity of the droplet, molecules have a high probability to reach the surface, P arrival ≈ 1, and do so at an entry angle close to the starting one, θ arrival ≈ θ 0 , i. e., cos θ arrival ≈ 1.
For starting positions further away from the droplet, the arrival probability is lower and the spread of the entry angle is larger, which also diminishes the average ⟨cos θ arrival ⟩.Comparing different reaction rates (corresponding to life times), while the arrival probability is lower for smaller life times, the entry angle, θ arrival , is spread less such that the assignment functions for all reaction rates (across a wide range of values) collapse onto a universal curve which is nicely described by an effective power law, ∼ r −2 .This motivates the assignment function Notice that while the assignment function does not explicitly depend on the reaction rate anymore, it still does implicitly because droplet sizes, R, are smaller at higher reaction rates.In this section, we compare the 3D passive droplets in external concentration gradients of Sec.Passive Droplets in External Concentration Gradients of the main text to their 1D counterpart.As we have seen, in 3D the flux can flow around the droplet in the case where the droplet is unfavorable for the F -solvent, or preferentially through the droplet, in the case where the droplet region is attractive for it.This is at the expense of the regions surrounding the droplet laterally, where the F -solvent flux is in-or decreased, respectively.In 1D such a geometric modification of the flux fields is not possible.In such a scenario the flux has to be constant in the bulk regions and jumps at the interface, if it moves with velocity v and there is a jump in concentration, for all components.Indeed this is what we observe in the 1D simulations.This is shown in Figure S3.Again, the F -solvent is depleted inside the droplet for neutral interactions, because of the asymmetry in molecular volume, while it enriched inside the droplet for favorable interactions.This is the same is in the 3D case.However, the F -flux behaves opposite to the 3D case: It is higher on the inside than on the outside for the depleted concentration and it is reversed for the enriched concentration.This relation is simply dictated by the movement of the droplet, which demands ∆j F = v∆ϕ F .This explains the higher F -flux inside the droplet for the case χ P F N P = 0 than on the outside (at the left interface: ∆ϕ F ≈ −0.05, v ≈ −0.1R e λ, ∆j F ≈ 0.005R e λ), and correspondingly the lower one inside for χ P F N P = −30 (∆ϕ F ≈ 0.14, v ≈ −0.12R e λ, ∆j F ≈ −0.018R e λ).Compared to the 3D case, here, all material of F and S necessarily must flow through the droplet.Thus, the droplet is a bottleneck in this scenario and influences the (otherwise constant) 1 flux in the whole simulation cell.That way the droplet is still a little bit faster in the attractive scenario than in the neutral one.

Gradients
For the simulations of Sec.Passive Droplets in External Concentration Gradients we can measure the droplet velocity and compare it the theoretical prediction above, which relates it to the fluxes and concentrations.
To do so, we perform the Hoshen-Kopelman cluster analysis (HKCA), 53,54 where a droplet is defined by the threshold ϕ P > 0.5.This allows to measure the droplets center of mass, r cm , and thereby its velocity using two consecutive time steps.Figure S4 depicts the measured velocities as a function of the center of mass.The velocities are approximately constant inside the bulk and approach zero once the droplet overlaps with the F -source, i. e., once its center of mass reaches z = R, with R being the droplet's radius.We compare the measured velocities to the theoretical prediction, which is worked out for the general case of reaction-driven assembly (RDA)-formed droplets in fuel and waste concentration gradients in Sec.Theoretical Prediction of Velocities, resulting in Equation S19.For the chemically passive droplets, r f = r b = 0, and assuming the concentration profiles vary slowly across the droplet volume, this simply reduces to v = j P (rcm) ϕ P (rcm) .Making use of the cluster analysis, we average the fluxes and concentrations within the droplets at every time step and use these averages to predict v.The result is plotted in Figure S4 as dashed lines and matches the measurements.
Nucleation of RDA-formed Droplets Here, we briefly demonstrate how RDA-formed droplets nucleate.We show this on the example system of Sec.Implicit Waste of the main text.The droplets move towards the sides of the simulation cell, towards higher fuel concentrations.In this setup, the collective movement creates large voids between droplets in which the forward reaction commences.
This creates more product in the solution than can diffuse into the surrounding droplets or can revert into the precursor state.Thus, the product concentration in these void begins to rise until it comes close to the spinodal concentration.When this happens, a droplet can stochastically nucleate and such a nucleation event becomes the probable the higher the super-saturation ϕ P (r) − ϕ P,binodal .Eventually, a new droplet nucleates, which quickly adsorbs the super-saturated product from its surrounding and grows to the stationary size.
Such a nucleation event is shown in Figure S5.
Note that the nucleation happens close to the (mean-field) spinodal, 51 i. e., when the free-energy barrier for the nucleation becomes very small (i.e. on the order of k B T ).The binodal was measured at ϕ P,binodal = 0.002 (with noise) and the mean-field spinodal at ϕ P,spinodal = 0.01.This is where we typically observe the nucleation process, i. e. at high supersaturation.Hence, formation of droplets occurs in our simulations at the transition between (classical) nucleation and spinodal decomposition.The simulation results for the three system are shown in Figure S6.For system I, implicit fuel with a large simulation cell, more precursor and no wall, we observe the nucleation in the center of the system and subsequent directed movement towards the closest fuel source.

Choice of Boundary Conditions
Instead of stopping in front of the source, droplet freely move to the sides of the simulation cell, where they fuse with oncoming droplets.This allows fusion to large structures that regularly reach beyond the periodic boundary conditions.Allowing the movement above the source also allows to block it which periodically leads to depletion of fuel in the system, reducing the product concentration until the source is free, again.This undesirable dynamics justifies the external boundary field that was applied in all other simulations.Additionally, the large simulation cell allows nucleation of droplets far away from the source.Here the droplets do not grow as quickly as in the case of a smaller simulation cell.Droplets stay small while they are far away from the source and only gradually grow upon approaching it.This becomes visible when looking at the droplet radii plotted against their center of mass.From the laterally averaged concentrations inside and outside of the droplets we can see that the precursor concentration is high in all of the solution.Measuring the velocity as a function of the center of mass results in a similar curve as the one from the main text, nicely matched by the theoretical prediction.Looking at the individual contribution to the prediction, the bottom panel, we realize that in this scenario, the reaction contributions are just as large as the flux contributions.
System II is closer to the system of the main text, with the difference of a doubled simulation-cell size and doubled refilling rate of the fuel, while we keep the polymer concentration low, φR + φP = 0.125.Due to the larger size of the simulation cell, we observe, again, a region of small droplet density in the center, where droplets that do nucleate stay small and gradually grow upon approaching the source.In comparison to system I, the precursor concentration becomes low in the center, while the fuel concentration does not fall off as fast.Compared to the smaller system size of the main text, the product is more concentrated at the walls, which results in larger structures that also regularly span beyond the periodic boundaries laterally.Again, the velocity prediction matches the measurements.
In this scenario the reaction contribution is comparable to the flux contribution, however smaller in the bulk.
System III features fuel and waste explicitly, with an affinity of product towards waste.
This scenario is tailored to show what happens if the simulation cell is larger than the length scale of the waste concentration profile and thus too large for droplets to approach in the center.As visible in the example snapshot and the radii plot of Figure S6, droplets nucleate close to the fuel source.From here they quickly move towards the center.At approximately z ≈ 10R e and z ≈ L z − 10R e multiple droplets are closely packed, in full correspondence to the scenario with small simulation-cell size.When droplets approach this region, their accumulation volume becomes asymmetric, with small volume closer the center and large volume closer to the source.This leads to large forward reaction contributions opposing the flux contribution, and the droplets are slowed down.In this region the droplets adopt the largest sizes, while they shrink as they move forward from here farther away from the source.
After moving through the densely packed region, the droplets distance from each other while flux contribution drives them towards the center.The droplets continuously shrink until they reach their smallest possible size, dictated by the interface width (R ≈ 0.3R e ), and at this point dissolve.In this scenario, the velocity measurements also display the 'shoulder', the curve obtaining a small gradient in the center.Here it is not only a result of the jamming and the reaction contribution opposing the flux contribution but also the flux contribution itself displays this behavior because the waste concentration becomes saturated with a vanishing gradient at the center.
These three scenarios demonstrate the robustness of the mechanism, which qualitatively does not depend on the boundary conditions of the simulation cell.Quantitatively, however, the dynamics depends on those details offering opportunities to tailor the behavior.

Reaction Contribution to Velocities
We have shown in the main text that the forward reaction contribution strongly varies, depending on the system's morphology, i. e., the arrangement of droplets and the distribution of precursor and fuel molecules.The movement that can be produced by the forward reaction is the result of an asymmetry in product that flows into the droplet.This can stem either from concentration gradients of fuel and precursor (i.e. their products) or from an asymmetry in the accumulation volume.

For the two examples of implicit waste or fuel of Sections Implicit Waste and Implicit
Fuel of the main text, let us have a look at the statistics of the accumulation volumes and the averaged concentrations.For this we perform the Voronoi tessellation of the domains for each time step and locally average the distinct domains.Assuming cylindrical symmetry of the averaged domains, the result is presented in the insets of Figure S7.Clearly, the Figure S7: Measurement results for the system of implicit waste (left) and implicit fuel (right).In the top row, we show the mean volume of the accumulation domains in front of, V (Ω − ), and behind, V (Ω + ), the center of mass.The insets show the (locally) averaged domains, assuming cylindrical symmetry (r = x 2 + y 2 ).The center of mass is marked as a cross to see the horizontal asymmetry in these.In the bottom row we show the temporally and laterally averaged product of concentrations ϕ R ϕ F , which are proportional to the product of the forward reaction (proportionality factor is, r f N P = 3.2N P λ on the left and r f N P = 0.32N P λ on the right).
accumulation volumes are asymmetric around the center of mass of the droplets: For implicit fuel and movement towards the sides of the simulation cell, the accumulation volumes are larger farther away from the source.In the opposite case, implicit fuel and movement towards the center, the accumulation volumes are larger at the side that is closer to the waste sink.
This motivates analyzing the volume of the accumulation domains in front of the center of mass, V (Ω − ) for z < z cm , and behind it, V (Ω + ) for z > z cm .The result is depicted in the top row of Figure S7.The bottom row shows the temporally and laterally averaged product of fuel and precursor concentrations outside of the droplets, ϕ R ϕ F ω.In the case of implicit waste and movement towards the side, the precursor concentration is also higher at the sides of the simulation cell, as is the fuel concentration, and there is a high gradient of the product concentration.Hence, within a infinitesimal volume around droplet, the production is higher closer to the source than further away from it.It is the opposite for the scenario where droplets move towards the center, although the fuel concentration (solvent in this scenario) still is highest at the sides.Hence the two contributions to the reaction contribution oppose each other.The result, the green lines of Figures 5 and 6, is that the reaction contribution works in the direction of the flux contribution in the first case, while it works in opposite direction of it (away from the center) in the case where droplets move away from the waste sink.
These two scenarios are representative also of the full systems, where both, fuel and waste, are treated explicitly.

Fuel and Waste Fluxes Cancel in the Symmetric Case
When interactions and mobilities of fuel and waste are the same, their effects on the droplets cancel exactly in the statistically stationary state, such that the droplets do not move.
For this argument, consider a system where no phase separation occurs.For simplicity we only consider fuel, waste and solvent, treating precursor and product implicitly, and analyze the resulting solvent flux.We take all mobilities to be the same, λ c = λ ∀c and all interactions to vanish, χ cc ′ = 0 ∀c, c ′ .Moreover, we disregard thermal fluctuations.For the fuel component, the time evolution becomes where r f (r) ) is now a general local forward reaction rate, and in the second line we consider only small concentration gradients, allowing us to disregard the square-gradient penalty in the chemical potential.The time evolution of the waste is analogous, In the stationary state, the resulting solvent flux takes the form Since waste and fuel flux cancel exactly, the main contribution to the droplet velocity vanishes.
Indeed, such a system shows an arrest in movement and approaches a stationary state where large droplets arrange on a lattice structure close to the interface and droplet sizes decrease further away from the source.Here, the classical size-control mechanism of such RDA systems becomes visible, dictating the lattice spacing, and size of droplets then being determined by the (local) mean concentration of product, which is smaller further away from the fuel source.The simulation results and analysis are shown in Figure S8.While the flux inside the droplets vanishes entirely, the reaction contributions do not, i. e., the asymmetry in accumulation volumes does not entirely counter the asymmetric production because of concentration gradients.Here, the deviation between prediction and measurement is again prominent close to the source.

Simulation Results for Explicit Fuel and Waste
For explicit treatment of fuel and waste we vary the interactions of fuel and waste with the product, χ P F and χ P W , and the reaction rate constants, r b .Performing all measurements, as before, we obtain the simulation results in Figure S9.Panel (a) shows the measured droplet radii plotted against the droplets' center of mass.Qualitatively one can observe nucleation in the center for the attractive-fuel case.χ P F < 0, while the nucleation happens near the source in the opposite case, χ P W < 0. For large affinities one can also observe fusion events where the droplets move to, while such events are not observable for intermediate attractions (lower velocities).The droplets become larger for smaller a reaction rate constant, r b = 0.04λ.The velocity measurements and their theoretical comparisons are plotted in  panel (b).Clearly strong attractions lead to higher velocities, as well as a higher reaction rate, which is associated with higher concentration gradients.The linear fit in the bulk, 5 ≤ z ≤ 21.6R e , is indicated as dashed lines for both measurement and theory.Its slope is used as the indicator for the droplet dynamics in the main text.Panel (c) dissects the theoretical comparison into the three contributions and we obtain the same results, as in the cases of implicit treatment.The flux contribution determines the overall dynamics of the droplets which becomes stronger, the more attractive one of the component becomes.This determines the arrangement of droplets and how these share the product that is created in the surrounding solution.That way, the forward reaction contribution is enslaved to the flux contribution and acts as a small correction to the velocity.The backward reaction contribution is always negligible.

Particle-based Simulations
To corroborate the continuum model results, we perform particle-based simulations using a soft, coarse-grained model.We use the Monte-Carlo simulation program SOft coarsegrained Monte carlo Acceleration (SOMA) 60 which employs the single-chain-in-mean-field (SCMF) algorithm. 59The model treats molecules as Gaussian chains, with a discretization N i , lumping several monomer repeat units into a single interaction center and employs a top-down approach to assess the interactions and dynamics.The model allows for great parallelization by splitting interactions into strong bonded interaction and weak non-bonded where m runs over all molecules, b over all bonds within the molecule, and r m,b , r m,b+1 refer to the positions of the bonded particles.The weak non-bonded interactions are expressed in terms of the normalized densities ϕ c (r), i. e., where κ 0 = 6 is the inverse isothermal compressibility which is large for the nearly incompressible system.A more detailed description of the particle-based model including the implementation of reactions is given in Refs. 52,60nveniently, the two models operate on the same set of parameters which simplifies comparison between the two models.In the particle-based simulations, we take the same parameters, with the exception that solvent, fuel and waste molecules are dimers, N S/F/W = 2, just as was done in. 52,61To comply with this, we also take a higher discretization of the precursor and product molecules, N R/P = 20.We take simulation cells of V = 10R e ×10R e × 25R e at a resolution of 8 grid cells per R e .We take a reaction rate constant of r b = 0.08τ −1 0 and vary the fuel and waste interactions with the product molecules.Since product and short fuel and waste molecules have different lengths, the invariant interaction parameter is χ P,F/W rather than χ P,F/W N P as was the case if both species had the same discretization.Therefore, we vary the interactions in the range χ P,F/W = 0, −2, −4 to match the interaction of the continuum model with half the discretization values.The results are presented in Figure S10, where we depict the temporally and laterally averaged concentrations in-and outside the droplets in (a) and the measured velocities in (b).Visibly, we obtain the same result as before.
The affinity towards one component determines the direction of the movement which can be amplified by high affinities.

Simulations of Multiple Droplet Types
The system presented in Sec.Multiple Droplet Types features two droplet molecule types, P and Q, which phase separate into two distinct droplet phases.Here, we present the In order to observe two types of droplets that do not wet each other, we take a strong repulsion of P and Q polymers, χ P Q N P = 200.We show that upon demixing from solution, the Q-rich droplets move towards the center, whereas the Prich droplets continue to move towards the fuel source.This is demonstrated in Figure 9 of the main text, where example snapshots and velocity measurements demonstrate that even the Q-rich droplets whose components are not involved in the reaction cycle, exhibit a directed movement.

Simulations of Amphiphilic Product Molecules
When introducing one component into the system that is amphiphilic, we need to extend our continuum model to include this.5][66] The amphiphilic diblock copolymer consists of two linear blocks with distinct monomer types, called A and B. For each of these blocks, we introduce a density field, as we do for a other single-block components, i. e., R, S, F, W . Introducing for molecular type index p and for the blocks of each molecule index i, the free-energy functional reads where the coefficients Ap, ij and C pij and the parameters f pi are dependent on the molecular architecture.For simple homopolymers and solvents, the theory reduces to the Flory-Huggins-de Gennes model 41 as introduced in the main text, A p = 0, C p = N P /N p and f p = 1.
For the diblock copolymer, the free-energy expansion to second order in the concentrations matches Ohta-Kawasaki free-energy. 67For these components, the coefficients become 52,61 where f B = 1 − f A = 0.7 is the block ratio and the definition The reactions are implemented as described in Ref. 52 which assumes a single reaction to convert the hydrophilic precursor into an amphiphilic molecule.This would correspond to short molecules with presumably weak amphiphilicity, where a single reaction can create such a drastic change in the molecule's interactions.Overall, we still describe the reaction cycle by a forward and backward rate constant, r f and r b , respectively.The system consists finally of a hydrophilic precursor, R, in solution, S, that reacts with a fuel molecule, F , to become amphiphilic, i. e., the tail, B, becomes hydrophobic, while the head group, A, stays hydrophilic.Waste, W , is produced within the forward reaction and the backward reaction from amphiphile to precursor happens spontaneously, modeled as a first order reaction.The chain discretizations are taken as in the main text, the amphiphilic molecules have a block ratio of f A = 1 − f B = 0.3, and we take the interaction matrix We take larger simulation-cell sizes of V = 51.2Re × 9.6R e × 9.6R e such that micelles do not grow too quickly after nucleation, but slowly grow upon moving from the bulk towards the source, since the head-groups interact favorably with the fuel.We take a backward reaction rate constant of r b = 0.04λ and refilling rate constants of r source = 0.02λR e , r sink = 10λR e .
The result is a movement of aggregates as depicted in Sec.Growing complex reaction-driven assemblies, where the slow growth allows the assembly of vesicles, that would not form if it the amphiphiles would just experience a quench into the phase separating parameter regime.
Indeed, performing simulations without fuel and waste and with either no reactions or first-order reactions between precursor and product, but taking all interactions to be the same, we observe no coarsening after the quench and initial nucleation of small micelles, see Figure S11.This also becomes visible when plotting micelle radii against time, where the sizes just stall at the initially unstable length scale and no fusion happens.1.An Extensible Markup Language (XML) input file is used specify the system architecture (including simulation cell size, grid spacing, external field) the different molecular components (including molecular volume and their initial mean concentrations), as well as the reaction dynamics.The software provides a script to create a Hierarchical Data Format version 5 (HDF5) file which is the input to the simulation.
2. The HDF5-file is given to the simulations which requires the presence of a compute unified device architecture by Nvidia® (CUDA)-graphics processing unit (GPU).The concentration fields are regularly saved according to an interval, ∆τ , as specified in the XML input file.
The HDF5 file now contains the simulation results in the form of the time-dependent (n t saved times) density fields for all n c components.In the file a n t × n c × n z × n y × n x sized array is saved.

Analysis
The analysis takes the time-dependent data and performs the following steps 1. Perform a HKCA on all time-steps, defining a cluster where ϕ P > 0.5.The result, a cluster map, is saved in the HDF5-file as a n t × n z × n y × n x -sized array which is zero if no cluster is detected and non-zero, indicating the cluster index, otherwise.
2. Using the HKCA, calculate the cluster radii, R, as described in Ref. 52 and the center of mass, r cm .
3. Using the center of mass, compute the Voronoi tessellation by determining for each grid cell the cluster with the closest center of mass.The result is stored to the HDF5-file.
4. If it is of interest, the Voronoi-domains can be analyzed in detail.For this the domain is separated into 51 bins (bin width of ∼ 0.5R e ).Within each bin, we can averaged each Voronoi domain within the center-of-mass frame.From this average, the averaged volume in front of (z < z cm ) or behind (z > z cm ) the droplet can be computed.
5. Using the HKCA, the concentration fields can be averaged in-and outside of the droplet phases.The cluster analysis separates the droplet from its outside exactly at the interface.Hence, simply taking averaged concentration where the cluster map is non-zero results in serious miss-evaluation of the inside concentrations of droplets, if these are small (high surface-to-volume ratio).Instead we refine the cluster, making them smaller by excluding grid cells that are (next-nearest-)neighbors of non-cluster grid cells.Averaging over the refined clusters now excludes the interface, leading to clean concentration profiles.For the analysis of the outside concentration profiles the opposite refinement is performed.
6.The velocities are measured using the center-of-mass measurement.For this, the same droplet within two time-frames are detected by finding for each cluster (droplet) a cluster with the closest center of mass (with distance smaller than |r cm,i,t −r cm,i ′ ,t+∆τ | !< R e for cluster i at time t and cluster i ′ at time t + ∆τ ).The difference in center of mass then determines the velocity.
7. Using Equation 4, the predictions are calculated using the above analysis.To obtain the flux inside the droplets, the noisy data has to be filtered (to obtain clean gradients).
To this end, we run the continuum model simulation without thermal noise for a short time, tλ = 0.1 which is long enough to properly filter the fluctuations and short enough to not alter the concentration fields significantly.The fluxes inside the droplets are measured using the cluster analysis.The reaction contribution uses the (smoothed) local concentration fields in combination with the Voronoi tessellation.

Figure S1 :
Figure S1: Sketch of the contribution of production effects in fuel and precursor concentration gradients.Left: Visualization of the second contribution of Equation S11.What is produced as excess to the average in Ω (for infinitesimal time step dt) leads to movement of the interface, vdt on the left to the left.Correspondingly, the lack of production compared to the average leads to a retraction of the interface on the right (top).In the corresponding (reactive) product flux profiles there are jumps at the interfaces corresponding to the movement, which arises from an excess production on the left (bottom).The reactive flux is set to 0 in the instantaneous condensation approximation and replaced by surface source and sink shown at the top.Right: In 3D, droplets (red) are assigned an accumulation volume, Ω, shown by dashed lines, which corresponds to the Voronoi domain.The contribution of the local production, ζ(|r − r cm |) cos(θ), to the velocity in z-direction, v z = v is indicated by the heat map (blue negative, red positive).

Figure S1 .Figure S2 :
Figure S1.The figure also presents a visualization of the second contribution to the velocity.Additionally, we account for the random, diffusive motion of molecules before reaching the droplet via the assignment ζ(|r − r cm |).Without it, molecules would be assumed to reach the droplet at the angle, at which they react to become product.In reality the random motion will blur the position of entry into the droplet.To quantify this factor, we performed simple simulations of a Wiener process, starting at position r 0 = rê z , i. e., θ 0 = 0 and measured the entry angle θ arrival at which the molecules would reach a sphere of radius R, located at the origin.This allows to calculate the average angle contribution ⟨cos θ arrival ⟩ and their probability to reach the droplet, P arrival .The assignment function is then given by ζ(r) = P arrival (r)⟨cos θ arrival ⟩(r).The result is shown in FigureS2.Starting within close

Figure S3 :
Figure S3: Simulation result of the 1D passive droplets in external concentration gradients, for two parameter combinations and without thermal noise.In both cases there is a strong repulsion χ P S N P = 30, with neutral or attractive interactions between F -solvent and droplet material, as indicated at the top.We show a snapshot at tλ = 60 after placing the droplet at the center of the simulation cell, z cm,t=0 = 12.4R e = Lz 2 − 0.4R e .Given are concentration profiles (top) and fluxes (bottom).

Figure S4 :
Figure S4: Droplet velocity in z-direction, v, plotted against center of mass, z cm , both as measured from simulation and with the theoretical prediction using the measured fluxes and concentrations inside the droplet.

Figure S5 :
Figure S5: Nucleation dynamics of RDA-formed droplets for the setup with implicit waste of Sec.Implicit Waste in the main text.The different dynamics are indicated.

Figure S6 :
Figure S6: Simulation results for various systems, numbered as described in the text.For each column, the top panel shows an example snapshot of the product concentration field.In the second row, we show the temporally and laterally averaged concentration profiles inside (solid) and outside (dashed) of the droplets.The third row displays the measured droplet radii, plotted against their center of mass.The fourth row shows measured and theoretically predicted velocities plotted against center of mass.Individual measurements are shown as dots and the local averaged is marked with errorbars, indicating the standard deviation of the local distribution.Finally, the fifth panel does the same for the three contributions of the theoretical prediction.

Figure S8 :
Figure S8: Simulation results for explicit treatment of fuel and waste with symmetric interactions, χ P F = χ P W = 0, showing (a) an example snapshot after time tλ = 1900, after the system has become stationary.(b) The averaged concentration profiles inside and outside the droplets.Due to the lattice structure, in between the droplet layers, there is also no data on the inside concentrations.(c) Measured and theoretical velocities plotted against center of mass.

Figure S9 :
Figure S9: Simulation results for explicit treatment of fuel and waste for multiple system parameters, varying fuel and waste interactions with the product, χ P F N P and χ P W N P , along the columns, and reactions rate constants, r b , along rows.Among the interactions, one interaction (not indicated in the figure) is set to 0 (nonpreferential).Each quantity is plotted against the droplet center of mass, z cm .(a) Measured droplet radius, R. (b) Measured and theoretically computed droplet velocity, v, using Equation 4. A linear fit in the bulk, 5 ≤ z ≤ 21.6R e , is shown for each curve.(c) Dissection of theoretical velocity into its three contributions.In (b) and (c) the colors and markers are the same as in Figure 5.

Figure S10 :
Figure S10: Results of the particle-based simulations for reaction rate constant r b = 0.08τ −1 0 and variation of interaction parameters.(a) Domain-averaged densities for selected components.(b) Measured individual (light blue circles) and averaged (blue squares) velocities.

Figure S11 :
Figure S11: Simple quench test of amphiphiles in equilibrium with the interactions kept the same.Left: Morphology of tail-groups, ϕ B (r) after time tλ = 500.Right: Time evolution of micelle radii.